*IF DEF,A15_1A PVTHIN1A.2
C ******************************COPYRIGHT****************************** GTS2F400.7849
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.7850
C GTS2F400.7851
C Use, duplication or disclosure of this code is subject to the GTS2F400.7852
C restrictions as set forth in the contract. GTS2F400.7853
C GTS2F400.7854
C Meteorological Office GTS2F400.7855
C London Road GTS2F400.7856
C BRACKNELL GTS2F400.7857
C Berkshire UK GTS2F400.7858
C RG12 2SZ GTS2F400.7859
C GTS2F400.7860
C If no contract has been raised with this copy of the code, the use, GTS2F400.7861
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.7862
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.7863
C Modelling at the above address. GTS2F400.7864
C ******************************COPYRIGHT****************************** GTS2F400.7865
C GTS2F400.7866
CLL SUBROUTINE PV_THINT ----------------------------------------------- PVTHIN1A.3
CLL PVTHIN1A.4
CLL Purpose: Interpolates various fields to the desired theta level. PVTHIN1A.5
CLL Evaluates P_ON_THETA and RS_ON_THETA on the p-grid, then PVTHIN1A.6
CLL U_ON_THETA, V_ON_THETA and RS_UV_ON_THETA on the uv-grid. PVTHIN1A.7
CLL Calculates the derivative D(theta)/D(p). PVTHIN1A.8
CLL Note that routine assumes that theta is monotonic. PVTHIN1A.9
CLL It starts from the bottom of the atmosphere and PVTHIN1A.10
CLL moves up through the atmosphere until the value of theta PVTHIN1A.11
CLL is less than the value it is looking for. PVTHIN1A.12
CLL It then does the interpolation between that level PVTHIN1A.13
CLL and the level below to get the value on the desired level. PVTHIN1A.14
CLL PVTHIN1A.15
CLL Not suitable for single column use PVTHIN1A.16
CLL PVTHIN1A.17
CLL 6/10/92 Written By Simon Anderson PVTHIN1A.18
CLL PVTHIN1A.19
CLL Model Modification history from model version 3.0: PVTHIN1A.20
CLL version Date PVTHIN1A.21
CLL 3.2 28/07/93 Change subroutine name to uppercase for TS280793.10
CLL portability. Tracey Smith TS280793.11
!LL 4.3 17/02/97 Added ARGFLDPT arguments and MPP code P.Burton GSM3F403.969
CLL PVTHIN1A.22
CLL Programming standard UM DOC Paper 3, Version 4(05/02/92) A PVTHIN1A.23
CLL PVTHIN1A.24
CLL Logical Component Covered: D415 PVTHIN1A.25
CLL PVTHIN1A.26
CLL Project Task: D4 PVTHIN1A.27
CLL PVTHIN1A.28
CLL Documentation: U.M.D.P No 13. Derivation and Calculation of PVTHIN1A.29
CLL Unified Model Potential Vorticity. PVTHIN1A.30
CLL by Simon Anderson and Ian Roulstone. PVTHIN1A.31
CLL PVTHIN1A.32
CLLEND------------------------------------------------------------------ PVTHIN1A.33
PVTHIN1A.34
C*L ARGUMENTS:---------------------------------------------------------- PVTHIN1A.35
SUBROUTINE PV_THINT 1,8TS280793.12
1 (pstar,theta,rs,u,v,p_field,u_field, PVTHIN1A.37
2 p_levels,row_length, PVTHIN1A.38
*CALL ARGFLDPT 
GSM3F403.970
3 rmdi,ak,bk,des_theta, PVTHIN1A.39
4 eta_level,rs_on_theta,rs_uv_on_theta, TD141293.17
& e_levels,dth_dph,n_levels, TD141293.18
5 u_on_theta,v_on_theta,dtheta_dp) PVTHIN1A.41
PVTHIN1A.42
implicit none PVTHIN1A.43
PVTHIN1A.44
C Input variables ------------------------------------------------------ PVTHIN1A.45
PVTHIN1A.46
integer PVTHIN1A.47
& p_field !IN Points in horizontal p field. PVTHIN1A.48
&,u_field !IN Points in horizontal u field. PVTHIN1A.49
&,p_levels !IN Number of pressure levels. PVTHIN1A.50
&,row_length !IN Number of points in a row. PVTHIN1A.51
& ,n_levels !IN Number of levels of dtheta/dp TD141293.19
*CALL TYPFLDPT GSM3F403.971
PVTHIN1A.52
real PVTHIN1A.53
& pstar(p_field) !IN Primary model array for surf. press. PVTHIN1A.54
&,theta(p_field,p_levels) !IN Primary model array for theta field. PVTHIN1A.55
&,rs(p_field,p_levels) !IN Primary model array for rs. PVTHIN1A.56
&,u(u_field,p_levels) !IN Primary model array for u field. PVTHIN1A.57
&,v(u_field,p_levels) !IN Primary model array for v field. PVTHIN1A.58
PVTHIN1A.59
real PVTHIN1A.60
& rmdi !IN Real missing data indicator. PVTHIN1A.61
&,ak(p_levels) !IN A coefficient of hybrid coordinates PVTHIN1A.62
& ! at full levels. PVTHIN1A.63
&,bk(p_levels) !IN B coefficient of hybrid coordinates PVTHIN1A.64
& ! at full levels. PVTHIN1A.65
&,des_theta !IN Desired theta level we want PVTHIN1A.66
& ! variables interpolated onto. PVTHIN1A.67
& ,e_levels(n_levels) !IN half-levels over range TD141293.20
& ,dth_dph(p_field,n_levels) !IN dtheta/dp half-levels TD141293.21
PVTHIN1A.68
PVTHIN1A.69
C Output variables ----------------------------------------------------- PVTHIN1A.70
PVTHIN1A.71
real PVTHIN1A.72
& eta_level(p_field) !OUT eta value of theta level TD141293.22
&,rs_on_theta(p_field) !OUT Interpolated rs field on theta level PVTHIN1A.74
&,rs_uv_on_theta(u_field) !OUT Interpolated rs field on theta level PVTHIN1A.75
&,u_on_theta(u_field) !OUT Interpolated u field on theta level PVTHIN1A.76
&,v_on_theta(u_field) !OUT Interpolated v field on theta level PVTHIN1A.77
&,dtheta_dp(p_field) !OUT Calculated derivative D(theta)/D(p) PVTHIN1A.78
PVTHIN1A.79
C*---------------------------------------------------------------------- PVTHIN1A.80
C*L Workspace usage: PVTHIN1A.81
integer PVTHIN1A.82
& base_level_eta(p_field) ! Level pointer below desired level TD141293.23
& ! for each atmospheric column (set to PVTHIN1A.84
& ! 0 if not found, and p_levels if PVTHIN1A.85
& ! des_theta is above top level). PVTHIN1A.86
& ! Calculated at p-points. PVTHIN1A.87
&,base_level_uv(u_field) ! Base_level calculated at uv-points. PVTHIN1A.88
PVTHIN1A.89
real PVTHIN1A.90
& theta_uv(u_field,p_levels) ! Theta interpolated to uv-points. PVTHIN1A.91
PVTHIN1A.92
C*---------------------------------------------------------------------- PVTHIN1A.93
C*L External Subroutine Calls: PVTHIN1A.94
external p_to_uv ! Interpolate from p grid to u grid. PVTHIN1A.95
PVTHIN1A.96
C*---------------------------------------------------------------------- PVTHIN1A.97
C*L Local variables: PVTHIN1A.98
integer i,j,ii ! Loop counts. PVTHIN1A.99
PVTHIN1A.100
real zth1,zth2,ze1,ze2,den1,den2 TD141293.24
PVTHIN1A.105
C ---------------------------------------------------------------------- PVTHIN1A.106
CL Section 1 : Find the value of the base level. This is the largest PVTHIN1A.107
CL value of level such that theta(level)<des_theta , PVTHIN1A.108
CL calculated here on p-points. PVTHIN1A.109
C ---------------------------------------------------------------------- PVTHIN1A.110
PVTHIN1A.111
do 110 i = FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.972
base_level_eta(i)=0 TD141293.25
110 continue PVTHIN1A.114
do 120 j = 1,p_levels PVTHIN1A.115
do 130 i = FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.973
if (des_theta.gt.theta(i,j)) then PVTHIN1A.117
base_level_eta(i)=j TD141293.26
endif PVTHIN1A.119
130 continue PVTHIN1A.120
120 continue PVTHIN1A.121
C When this loop is done, base_level_p will be the value of the level PVTHIN1A.122
C with the highest value of theta LESS than the desired theta value. PVTHIN1A.123
C Base_level_p is set to 0 if no smaller value is found. PVTHIN1A.124
C Base_level_p is set to p_levels if no larger value is found. PVTHIN1A.125
PVTHIN1A.126
C ---------------------------------------------------------------------- PVTHIN1A.127
CL Section 2 : This section will interpolate variables held PVTHIN1A.128
CL on the p-grid onto the desired theta level. PVTHIN1A.129
CL Used for P_ON_THETA and RS_ON_THETA at each point. PVTHIN1A.130
C ---------------------------------------------------------------------- PVTHIN1A.131
PVTHIN1A.132
C---- Given P and RS as functions of Theta and Des_theta lying between PVTHIN1A.133
C---- Zth2 and Zth3, calculate P and RS at Des_theta linearly. PVTHIN1A.134
PVTHIN1A.135
do 210 j=FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.974
PVTHIN1A.137
if(base_level_eta(j).lt.2.or.base_level_eta(j).gt.p_levels-1) TD141293.27
& then TD141293.28
eta_level(j)=rmdi TD141293.29
rs_on_theta(j)=rmdi PVTHIN1A.140
else PVTHIN1A.141
PVTHIN1A.142
ii=base_level_eta(j) TD141293.30
PVTHIN1A.144
C Calculate theta and pressure values at the required levels. PVTHIN1A.145
PVTHIN1A.146
zth1=theta(j,ii) TD141293.31
ze1=0.00001*ak(ii)+bk(ii) TD141293.32
PVTHIN1A.149
zth2=theta(j,ii+1) TD141293.33
ze2=0.00001*ak(ii+1)+bk(ii+1) TD141293.34
PVTHIN1A.152
eta_level(j)= (des_theta-zth2)*ze1/(zth1-zth2)+ TD141293.35
& (des_theta-zth1)*ze2/(zth2-zth1) TD141293.36
rs_on_theta(j)= (des_theta-zth2)*rs(j,ii)/(zth1-zth2)+ TD141293.37
& (des_theta-zth1)*rs(j,ii+1)/(zth2-zth1) TD141293.38
if(eta_level(j).gt.e_levels(ii))then TD141293.39
base_level_eta(j)=ii-1 TD141293.40
endif TD141293.41
end if PVTHIN1A.158
PVTHIN1A.159
210 continue PVTHIN1A.160
PVTHIN1A.161
C ---------------------------------------------------------------------- PVTHIN1A.162
CL Section 3 : Interpolate values of theta and rs_on_theta from the PVTHIN1A.163
CL p_grid to the uv-grid. Theta_uv will be used below. PVTHIN1A.164
C ---------------------------------------------------------------------- PVTHIN1A.165
PVTHIN1A.166
do i=1,p_levels PVTHIN1A.167
call p_to_uv PVTHIN1A.168
1 (theta(1,i),theta_uv(1,i),p_field,u_field, PVTHIN1A.169
*IF DEF,MPP GSM3F403.975
2 row_length,P_LAST_ROW+1) GSM3F403.976
*ELSE GSM3F403.977
2 row_length,p_field/row_length) GSM3F403.978
*ENDIF GSM3F403.979
end do PVTHIN1A.171
GSM3F403.980
*IF DEF,MPP GSM3F403.981
C Initialise unused rows before p_to_uv call GSM3F403.982
do i=TOP_ROW_START-1,1,-1 GSM3F403.983
rs_on_theta(i)=rs_on_theta(i+ROW_LENGTH) GSM3F403.984
enddo GSM3F403.985
do i=LAST_P_VALID_PT+1,P_FIELD GSM3F403.986
rs_on_theta(i)=rs_on_theta(i-ROW_LENGTH) GSM3F403.987
enddo GSM3F403.988
*ENDIF GSM3F403.989
GSM3F403.990
call p_to_uv PVTHIN1A.172
1 (rs_on_theta,rs_uv_on_theta,p_field,u_field, PVTHIN1A.173
*IF DEF,MPP GSM3F403.991
2 row_length,P_LAST_ROW+1) GSM3F403.992
*ELSE GSM3F403.993
2 row_length,p_field/row_length) PVTHIN1A.174
*ENDIF GSM3F403.994
PVTHIN1A.175
PVTHIN1A.176
C ---------------------------------------------------------------------- PVTHIN1A.177
CL Section 4 : Find the value of the base level. This is the largest PVTHIN1A.178
CL value of level such that theta(level)<des_theta, PVTHIN1A.179
CL calculated here on uv-points. PVTHIN1A.180
C ---------------------------------------------------------------------- PVTHIN1A.181
do 410 i = FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.995
base_level_uv(i) = 0 PVTHIN1A.183
410 continue PVTHIN1A.184
do 420 j = 1,p_levels PVTHIN1A.185
do 430 i = FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.996
if (des_theta.gt.theta_uv(i,j)) then PVTHIN1A.187
base_level_uv(i) = j PVTHIN1A.188
endif PVTHIN1A.189
430 continue PVTHIN1A.190
420 continue PVTHIN1A.191
C When this loop is done, base_level_uv will be the value of the level PVTHIN1A.192
C with the highest value of theta LESS than the desired theta value. PVTHIN1A.193
C Base_level_uv is set to 0 if no smaller value is found. PVTHIN1A.194
C Base_level_uv is set to p_levels if no larger value is found. PVTHIN1A.195
PVTHIN1A.196
PVTHIN1A.197
C ---------------------------------------------------------------------- PVTHIN1A.198
CL Section 5 : This section will interpolate variables held PVTHIN1A.199
CL on the u-grid onto the desired theta level. PVTHIN1A.200
CL Used for U_ON_THETA and V_ON_THETA at each point. PVTHIN1A.201
C ---------------------------------------------------------------------- PVTHIN1A.202
PVTHIN1A.203
C---- Given U, V and RS_UV as functions of Theta and Des_theta lying PVTHIN1A.204
C---- between Zth2 and Zth3, calculate U, V and RS_UV at Des_theta PVTHIN1A.205
C---- by linear interpolation. PVTHIN1A.206
C---- PVTHIN1A.207
C---- If level above p_levels-1 or in bottom of boundary layer then PVTHIN1A.208
C---- we set the value to missing data. PVTHIN1A.209
C---- PVTHIN1A.210
PVTHIN1A.211
do 510 j=FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.997
PVTHIN1A.213
if(base_level_uv(j).lt.2.or. TD141293.42
& base_level_uv(j).gt.p_levels-1 ) then PVTHIN1A.215
u_on_theta(j)=rmdi PVTHIN1A.216
v_on_theta(j)=rmdi PVTHIN1A.217
eta_level(j)= rmdi TD141293.43
else PVTHIN1A.218
ii = base_level_uv(j) PVTHIN1A.219
PVTHIN1A.220
C Calculate U_ON_THETA and V_ON_THETA. PVTHIN1A.221
PVTHIN1A.222
u_on_theta(j) = ((theta_uv(j,ii+1)-des_theta)*u(j,ii) + PVTHIN1A.223
& (des_theta-theta_uv(j,ii))*u(j,ii+1)) PVTHIN1A.224
& /(theta_uv(j,ii+1)-theta_uv(j,ii)) PVTHIN1A.225
PVTHIN1A.226
v_on_theta(j) = ((theta_uv(j,ii+1)-des_theta)*v(j,ii)+ PVTHIN1A.227
& (des_theta-theta_uv(j,ii))*v(j,ii+1)) PVTHIN1A.228
& /(theta_uv(j,ii+1)-theta_uv(j,ii)) PVTHIN1A.229
PVTHIN1A.230
end if PVTHIN1A.231
PVTHIN1A.232
510 continue PVTHIN1A.233
PVTHIN1A.234
PVTHIN1A.235
C ---------------------------------------------------------------------- PVTHIN1A.236
CL Section 6 : This section will calculate derivatives held PVTHIN1A.237
CL on the p-grid on the desired theta level. PVTHIN1A.238
CL Used for DTHETA_DP at each point. PVTHIN1A.239
CL The same theta and pressure values are used as PVTHIN1A.240
CL were calculated in Section 2. PVTHIN1A.241
C ---------------------------------------------------------------------- PVTHIN1A.242
PVTHIN1A.243
C---- Given P as a function of Theta and Des_theta lying between PVTHIN1A.244
C---- Zth1 and Zth4, calculate DTHETA_DP at Des_theta by calculating PVTHIN1A.245
C---- dtheta_dp using centred finite-differences about each point PVTHIN1A.246
C---- zth1 to zth4 and then using cubic Lagrange interpolation. PVTHIN1A.247
C---- PVTHIN1A.248
C---- If level above p_levels-1 or in bottom of boundary layer then PVTHIN1A.249
C---- we set the value to missing data. PVTHIN1A.250
C---- PVTHIN1A.251
PVTHIN1A.252
do 610 j=FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.998
PVTHIN1A.254
if(base_level_eta(j).lt.2.or.base_level_eta(j).gt.p_levels-1) TD141293.44
& then TD141293.45
eta_level(j)= rmdi TD141293.46
dtheta_dp(j)=rmdi PVTHIN1A.256
elseif(base_level_eta(j).eq.p_levels-1)then TD141293.47
ii=base_level_eta(j) TD141293.48
dtheta_dp(j)=dth_dph(j,ii) TD141293.49
else TD141293.50
ii=base_level_eta(j) TD141293.51
PVTHIN1A.259
PVTHIN1A.263
C Calculate dtheta_dp and pressure values at the PVTHIN1A.264
C required four levels. PVTHIN1A.265
ze1=e_levels(ii) TD141293.52
ze2=e_levels(ii+1) TD141293.53
PVTHIN1A.281
C Calculate denominators required. PVTHIN1A.282
den1=ze1-ze2 TD141293.54
den2=ze2-ze1 TD141293.55
PVTHIN1A.287
C Calculate DTHETA_DP. PVTHIN1A.288
dtheta_dp(j) = PVTHIN1A.289
& (eta_level(j)-ze2)*dth_dph(j,ii)/den1 TD141293.56
& +(eta_level(j)-ze1)*dth_dph(j,ii+1)/den2 TD141293.57
PVTHIN1A.298
end if PVTHIN1A.299
PVTHIN1A.300
610 continue PVTHIN1A.301
PVTHIN1A.302
*IF DEF,MPP GSM3F403.999
! Swap halos on all fields GSM3F403.1000
CALL SWAPBOUNDS(eta_level,row_length,tot_P_ROWS, GSM3F403.1001
& EW_Halo,NS_Halo,1) GSM3F403.1002
CALL SWAPBOUNDS(rs_on_theta,row_length,tot_P_ROWS, GSM3F403.1003
& EW_Halo,NS_Halo,1) GSM3F403.1004
CALL SWAPBOUNDS(dtheta_dp,row_length,tot_P_ROWS, GSM3F403.1005
& EW_Halo,NS_Halo,1) GSM3F403.1006
CALL SWAPBOUNDS(rs_uv_on_theta,row_length,tot_U_ROWS, GSM3F403.1007
& EW_Halo,NS_Halo,1) GSM3F403.1008
CALL SWAPBOUNDS(u_on_theta,row_length,tot_U_ROWS, GSM3F403.1009
& EW_Halo,NS_Halo,1) GSM3F403.1010
CALL SWAPBOUNDS(v_on_theta,row_length,tot_U_ROWS, GSM3F403.1011
& EW_Halo,NS_Halo,1) GSM3F403.1012
*ENDIF GSM3F403.1013
PVTHIN1A.303
return PVTHIN1A.304
end PVTHIN1A.305
PVTHIN1A.306
*ENDIF PVTHIN1A.307