*IF DEF,A15_1A PVPIN1A.2
C ******************************COPYRIGHT****************************** GTS2F400.7831
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.7832
C GTS2F400.7833
C Use, duplication or disclosure of this code is subject to the GTS2F400.7834
C restrictions as set forth in the contract. GTS2F400.7835
C GTS2F400.7836
C Meteorological Office GTS2F400.7837
C London Road GTS2F400.7838
C BRACKNELL GTS2F400.7839
C Berkshire UK GTS2F400.7840
C RG12 2SZ GTS2F400.7841
C GTS2F400.7842
C If no contract has been raised with this copy of the code, the use, GTS2F400.7843
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.7844
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.7845
C Modelling at the above address. GTS2F400.7846
C ******************************COPYRIGHT****************************** GTS2F400.7847
C GTS2F400.7848
CLL SUBROUTINE pv_pint ------------------------------------------------ PVPIN1A.3
CLL PVPIN1A.4
CLL Purpose: Interpolates various fields to a desired pressure level. PVPIN1A.5
CLL Calculates the derivatives D(rs.u)/D(p), D(rs.v)/D(p) PVPIN1A.6
CLL and D(theta)/D(p). PVPIN1A.7
CLL Note that routine assumes that pressure is monotonic. PVPIN1A.8
CLL PVPIN1A.9
CLL Not suitable for single column use. PVPIN1A.10
CLL PVPIN1A.11
CLL Model Modification history: PVPIN1A.12
CLL Version Date PVPIN1A.13
CLL 3.1 12/11/92 Written By Simon Anderson. PVPIN1A.14
CLL 3.1 18/01/93 New deck at the release of version 3.1. PVPIN1A.15
CLL 3.2 28/07/93 Change subroutine name to uppercase for TS280793.7
CLL portability. Tracey Smith TS280793.8
!LL 4.3 17/02/97 Added ARGFLDPT arguments and MPP code P.Burton GSM3F403.930
CLL PVPIN1A.16
CLL Programming standard UM DOC Paper 3, Version 4(05/02/92) A PVPIN1A.17
CLL PVPIN1A.18
CLL Logical Component Covered: D415 PVPIN1A.19
CLL PVPIN1A.20
CLL Project Task: D4 PVPIN1A.21
CLL PVPIN1A.22
CLL Documentation: U.M.D.P. No 13. Derivation and Calculation of PVPIN1A.23
CLL Unified Model Potential Vorticity. PVPIN1A.24
CLL By Simon Anderson and Ian Roulstone. PVPIN1A.25
CLL PVPIN1A.26
CLLEND------------------------------------------------------------------ PVPIN1A.27
PVPIN1A.28
C*L ARGUMENTS:---------------------------------------------------------- PVPIN1A.29
SUBROUTINE PV_PINT 1,10TS280793.9
1 (pstar,theta,rs,u,v,p_field,u_field, PVPIN1A.31
2 p_levels,row_length, PVPIN1A.32
*CALL ARGFLDPT 
GSM3F403.931
3 rmdi,ak,bk,des_press, PVPIN1A.33
& eta_level,e_levels,dth_dph,n_levels, TD141293.120
4 theta_on_press,rs_on_press, PVPIN1A.34
5 u_on_press,v_on_press, PVPIN1A.35
6 drsu_dp,drsv_dp,dtheta_dp) PVPIN1A.36
PVPIN1A.37
implicit none PVPIN1A.38
PVPIN1A.39
C Input variables ------------------------------------------------------ PVPIN1A.40
PVPIN1A.41
integer PVPIN1A.42
& p_field !IN Points in horizontal p field. PVPIN1A.43
&,u_field !IN Points in horizontal u field. PVPIN1A.44
&,p_levels !IN Number of pressure levels. PVPIN1A.45
&,row_length !IN Number of points in a row. PVPIN1A.46
& ,n_levels !IN Number of levels of spline TD141293.121
PVPIN1A.47
*CALL TYPFLDPT GSM3F403.932
real PVPIN1A.48
& pstar(p_field) !IN Primary model array for surf. press. PVPIN1A.49
&,theta(p_field,p_levels) !IN Primary model array for theta field. PVPIN1A.50
&,rs(p_field,p_levels) !IN Primary model array for rs. PVPIN1A.51
&,u(u_field,p_levels) !IN Primary model array for u field. PVPIN1A.52
&,v(u_field,p_levels) !IN Primary model array for v field. PVPIN1A.53
& ,e_levels(n_levels) !IN half-levels over range TD141293.122
& ,dth_dph(p_field,n_levels) !IN dtheta/dp half-levels TD141293.123
PVPIN1A.54
real PVPIN1A.55
& rmdi !IN Real missing data indicator. PVPIN1A.56
&,ak(p_levels) !IN A coefficient of hybrid coordinates PVPIN1A.57
& ! at full levels. PVPIN1A.58
&,bk(p_levels) !IN B coefficient of hybrid coordinates PVPIN1A.59
& ! at full levels. PVPIN1A.60
&,des_press !IN Desired pressure level we want PVPIN1A.61
& ! variables interpolated onto. PVPIN1A.62
PVPIN1A.63
PVPIN1A.64
C Output variables ----------------------------------------------------- PVPIN1A.65
PVPIN1A.66
real PVPIN1A.67
& theta_on_press(p_field) !OUT Interpolated theta field on p level. PVPIN1A.68
&,rs_on_press(p_field) !OUT Interpolated rs field on press level PVPIN1A.69
&,u_on_press(u_field) !OUT Interpolated u field on press level. PVPIN1A.70
&,v_on_press(u_field) !OUT Interpolated v field on press level. PVPIN1A.71
&,drsu_dp(u_field) !OUT Calculated derivative D(rs.u)/D(p). PVPIN1A.72
&,drsv_dp(u_field) !OUT Calculated derivative D(rs.v)/D(p). PVPIN1A.73
&,dtheta_dp(p_field) !OUT Calculated derivative D(theta)/D(p). PVPIN1A.74
& ,eta_level(p_field) !OUT eta value of theta level TD141293.124
PVPIN1A.75
C*---------------------------------------------------------------------- PVPIN1A.76
C*L Workspace usage: 4 arrays are required. PVPIN1A.77
integer PVPIN1A.78
& base_level_eta(p_field) ! Level pointer below desired level TD141293.125
& ! level for each atmospheric column (set PVPIN1A.80
& ! to 0 if not found, and p_levels if PVPIN1A.81
& ! des_press is above top level). PVPIN1A.82
& ! Calculated at p-points. PVPIN1A.83
&,base_level_uv(u_field) ! Base_level calculated at uv-points. PVPIN1A.84
PVPIN1A.85
real PVPIN1A.86
& pstar_uv(u_field) ! Surf. Press. interpolated to uv-points. PVPIN1A.87
&,rs_uv(u_field,p_levels) ! Rs field interpolated to uv-points. PVPIN1A.88
PVPIN1A.89
C*---------------------------------------------------------------------- PVPIN1A.90
C*L External Subroutine Calls: PVPIN1A.91
external p_to_uv ! Interpolate from p-grid to uv-grid. PVPIN1A.92
PVPIN1A.93
C*---------------------------------------------------------------------- PVPIN1A.94
C*L Local variables: PVPIN1A.95
integer i,j,ii,iii ! Loop counts. PVPIN1A.96
PVPIN1A.97
integer imin,imax TD141293.126
real zth1,zth2,zp1,zp2,zp3,zp4 TD141293.127
real ze1,ze2,ze3,ze4,den1,den2 TD141293.128
PVPIN1A.101
C ---------------------------------------------------------------------- PVPIN1A.102
CL Section 1 : Find the value of the base level. This is the largest PVPIN1A.103
CL ~~~~~~~~~ value of level such that pressure(level)<des_press , PVPIN1A.104
CL calculated here on p-points. PVPIN1A.105
C ---------------------------------------------------------------------- PVPIN1A.106
do 110 i = FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.933
base_level_eta(i)=0 TD141293.129
110 continue PVPIN1A.109
do 120 j = 1,p_levels PVPIN1A.110
do 130 i = FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.934
if (des_press.lt. ak(j) + bk(j)*pstar(i)) then PVPIN1A.112
base_level_eta(i)=j TD141293.130
endif PVPIN1A.114
130 continue PVPIN1A.115
120 continue PVPIN1A.116
C When this loop is done, base_level_p will be the value of the level PVPIN1A.117
C with the highest value of p LESS than the desired pressure value. PVPIN1A.118
C Base_level_p is set to 0 if no larger value is found. PVPIN1A.119
C Base_level_p is set to p_levels if no smaller value is found. PVPIN1A.120
PVPIN1A.121
PVPIN1A.122
C ---------------------------------------------------------------------- PVPIN1A.123
CL Section 2 : This section will interpolate variables held PVPIN1A.124
CL ~~~~~~~~~ on the p-grid onto the desired pressure level. PVPIN1A.125
CL Used for THETA_ON_PRESS and RS_ON_PRESS at each point. PVPIN1A.126
C ---------------------------------------------------------------------- PVPIN1A.127
PVPIN1A.128
C---- Given Theta as a function of p and Des_press lying between PVPIN1A.129
C---- Zp1 and Zp5, calculate Theta at Des_press by fitting the PVPIN1A.130
C---- Quartic Lagrange polynomial through the data and evaluating at PVPIN1A.131
C---- Pressure=Des_press. PVPIN1A.132
C---- PVPIN1A.133
C---- If level above p_levels-1 or in bottom boundary layers then PVPIN1A.134
C---- we set the value to missing data. PVPIN1A.135
C---- PVPIN1A.136
C---- If the value calculated by the quartic is out of range PVPIN1A.137
C---- of the input values, then no quartic or single-valued quartic PVPIN1A.138
C---- is possible between these points and linear interpolation PVPIN1A.139
C---- is used between the two levels. This is usually caused by an PVPIN1A.140
C---- unstable or almost unstable model profile. PVPIN1A.141
C---- PVPIN1A.142
C---- A linear interpolation is used for rs. PVPIN1A.143
PVPIN1A.144
do 210 j= FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.935
PVPIN1A.146
if(base_level_eta(j).lt.2.or.base_level_eta(j).gt.p_levels-1) TD141293.131
& then TD141293.132
eta_level(j)=rmdi TD141293.133
theta_on_press(j)=rmdi PVPIN1A.149
rs_on_press(j)=rmdi PVPIN1A.150
else PVPIN1A.151
PVPIN1A.152
PVPIN1A.192
C Reset iii for linear interpolation. PVPIN1A.193
iii=base_level_eta(j) TD141293.134
zth1=theta(j,iii) TD141293.135
zp1=ak(iii)+bk(iii)*pstar(j) TD141293.136
ze1=0.00001*ak(iii)+bk(iii) TD141293.137
zth2=theta(j,iii+1) TD141293.138
zp2=ak(iii+1)+bk(iii+1)*pstar(j) TD141293.139
ze2=0.00001*ak(iii+1)+bk(iii+1) TD141293.140
PVPIN1A.199
C Check to see if value is in range. PVPIN1A.200
eta_level(j)= (des_press-zp2)*ze1/(zp1-zp2)+ TD141293.141
& (des_press-zp1)*ze2/(zp2-zp1) TD141293.142
theta_on_press(j)=(des_press-zp2)*zth1/(zp1-zp2)+ TD141293.143
& (des_press-zp1)*zth2/(zp2-zp1) TD141293.144
rs_on_press(j)=(des_press-zp2)*rs(j,iii)/(zp1-zp2)+ TD141293.145
& (des_press-zp1)*rs(j,iii+1)/(zp2-zp1) TD141293.146
c reset half-level pointer to level below if desired TD141293.147
c level is below full level eta TD141293.148
if(eta_level(j).gt.e_levels(iii))then TD141293.149
base_level_eta(j)=iii-1 TD141293.150
endif TD141293.151
end if PVPIN1A.210
PVPIN1A.211
210 continue PVPIN1A.212
PVPIN1A.213
C ---------------------------------------------------------------------- PVPIN1A.214
CL Section 3 : Interpolate surface pressure and Rs values from PVPIN1A.215
CL ~~~~~~~~~ the p-grid to the uv-grid. PVPIN1A.216
C ---------------------------------------------------------------------- PVPIN1A.217
PVPIN1A.218
call p_to_uv PVPIN1A.219
1 (pstar,pstar_uv,p_field,u_field, PVPIN1A.220
*IF DEF,MPP GSM3F403.936
2 row_length,P_LAST_ROW+1) GSM3F403.937
*ELSE GSM3F403.938
2 row_length,p_field/row_length) PVPIN1A.221
*ENDIF GSM3F403.939
do i=1,p_levels PVPIN1A.222
call p_to_uv PVPIN1A.223
1 (rs(1,i),rs_uv(1,i),p_field,u_field, PVPIN1A.224
*IF DEF,MPP GSM3F403.940
2 row_length,P_LAST_ROW+1) GSM3F403.941
*ELSE GSM3F403.942
2 row_length,p_field/row_length) GSM3F403.943
*ENDIF GSM3F403.944
end do PVPIN1A.226
PVPIN1A.227
C ---------------------------------------------------------------------- PVPIN1A.228
CL Section 4 : Find the value of the base level. This is the largest PVPIN1A.229
CL ~~~~~~~~~ value of level such that pressure(level)<des_press, PVPIN1A.230
CL calculated here on uv-points. PVPIN1A.231
C ---------------------------------------------------------------------- PVPIN1A.232
do 410 i = FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.945
base_level_uv(i) = 0 PVPIN1A.234
410 continue PVPIN1A.235
do 420 j = 1,p_levels PVPIN1A.236
do 430 i = FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.946
if (des_press.lt. ak(j) + bk(j)*pstar_uv(i)) then PVPIN1A.238
base_level_uv(i) = j PVPIN1A.239
endif PVPIN1A.240
430 continue PVPIN1A.241
420 continue PVPIN1A.242
C When this loop is done, base_level_uv will be the value of the level PVPIN1A.243
C with the highest value of p LESS than the desired pressure value. PVPIN1A.244
C Base_level_uv is set to 0 if no larger value is found. PVPIN1A.245
C Base_level_uv is set to p_levels if no smaller value is found. PVPIN1A.246
PVPIN1A.247
PVPIN1A.248
C ---------------------------------------------------------------------- PVPIN1A.249
CL Section 5 : This section will interpolate variables held PVPIN1A.250
CL ~~~~~~~~~ on the uv-grid onto the desired pressure level. PVPIN1A.251
CL Used for U_ON_PRESS and V_ON_PRESS at each point. PVPIN1A.252
C ---------------------------------------------------------------------- PVPIN1A.253
PVPIN1A.254
C---- Given U and V as functions of P and Des_press lying PVPIN1A.255
C---- between Zp3 and Zp4, calculate U and V at Des_press by PVPIN1A.256
C---- by linear interpolation. PVPIN1A.257
C---- PVPIN1A.258
C---- If level is above p_levels-1 or in bottom boundary layers then PVPIN1A.259
C---- we set the value to missing data. PVPIN1A.260
C---- PVPIN1A.261
PVPIN1A.262
do 510 j=FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.947
PVPIN1A.264
if(base_level_uv(j).lt.2.or. TD141293.152
& base_level_uv(j).gt.p_levels-1) then TD141293.153
eta_level(j)= rmdi TD141293.154
u_on_press(j)=rmdi PVPIN1A.267
v_on_press(j)=rmdi PVPIN1A.268
else PVPIN1A.269
PVPIN1A.270
iii=base_level_uv(j) PVPIN1A.271
PVPIN1A.272
C Calculate pressure values at the required two levels. PVPIN1A.273
zp3 =ak(iii) + bk(iii)*pstar_uv(j) PVPIN1A.274
zp4 =ak(iii+1) + bk(iii+1)*pstar_uv(j) PVPIN1A.275
PVPIN1A.276
C Calculate U_ON_PRESS and V_ON_PRESS. PVPIN1A.277
u_on_press(j) = (des_press-zp4)*u(j,iii)/(zp3-zp4) + PVPIN1A.278
& (des_press-zp3)*u(j,iii+1)/(zp4-zp3) PVPIN1A.279
v_on_press(j) = (des_press-zp4)*v(j,iii)/(zp3-zp4) + PVPIN1A.280
& (des_press-zp3)*v(j,iii+1)/(zp4-zp3) PVPIN1A.281
end if PVPIN1A.282
PVPIN1A.283
510 continue PVPIN1A.284
PVPIN1A.285
PVPIN1A.286
C ---------------------------------------------------------------------- PVPIN1A.287
CL Section 6 : This section will calculate derivatives held PVPIN1A.288
CL ~~~~~~~~~ on the p-grid on the desired pressure level. PVPIN1A.289
CL Used for DTHETA_DP at each point. PVPIN1A.290
C ---------------------------------------------------------------------- PVPIN1A.291
PVPIN1A.292
C---- Given P as a function of Theta and Des_press lying between PVPIN1A.293
C---- Zp1 and Zp4, calculate DTHETA_DP at Des_press by calculating PVPIN1A.294
C---- dtheta_dp using centred finite-differences about each point PVPIN1A.295
C---- zp1 to zp4 and then using cubic Lagrange interpolation. PVPIN1A.296
C---- PVPIN1A.297
C---- If level above p_levels-1 or in bottom of boundary layer then PVPIN1A.298
C---- we set the value to missing data. PVPIN1A.299
C---- PVPIN1A.300
PVPIN1A.301
do 610 j=FIRST_VALID_PT,LAST_P_VALID_PT GSM3F403.948
PVPIN1A.303
if(base_level_eta(j).lt.2.or.base_level_eta(j).gt.p_levels-1) TD141293.155
& then TD141293.156
eta_level(j)= rmdi TD141293.157
dtheta_dp(j)=rmdi PVPIN1A.306
elseif(base_level_eta(j).eq.p_levels-1)then TD141293.158
c in penultimate half-layer set dtheta/dp to last value TD141293.159
ii=base_level_eta(j) TD141293.160
dtheta_dp(j)=dth_dph(j,ii) TD141293.161
else TD141293.162
ii=base_level_eta(j) TD141293.163
c interpolate dtheta/dp to desired level by using TD141293.164
c half-level values of dtheta/dp TD141293.165
ze1=e_levels(ii) TD141293.166
ze2=e_levels(ii+1) TD141293.167
den1=ze1-ze2 TD141293.168
den2=ze2-ze1 TD141293.169
C Calculate DTHETA_DP. PVPIN1A.338
dtheta_dp(j) = PVPIN1A.339
& (eta_level(j)-ze2)*dth_dph(j,ii)/den1 TD141293.170
& +(eta_level(j)-ze1)*dth_dph(j,ii+1)/den2 TD141293.171
end if PVPIN1A.348
PVPIN1A.349
610 continue PVPIN1A.350
PVPIN1A.351
PVPIN1A.352
C ---------------------------------------------------------------------- PVPIN1A.353
CL Section 7 : This section will calculate derivatives held PVPIN1A.354
CL ~~~~~~~~~ on the uv-grid on the desired pressure level. PVPIN1A.355
CL Used for DRSU_DP and DRSV_DP at each point. PVPIN1A.356
C ---------------------------------------------------------------------- PVPIN1A.357
PVPIN1A.358
do 710 j=FIRST_FLD_PT,LAST_U_FLD_PT GSM3F403.949
PVPIN1A.360
if(base_level_uv(j).lt.2.or. TD141293.172
& base_level_uv(j).gt.p_levels-1) then TD141293.173
eta_level(j)= rmdi TD141293.174
drsu_dp(j)=rmdi PVPIN1A.363
drsv_dp(j)=rmdi PVPIN1A.364
else PVPIN1A.365
PVPIN1A.366
iii=base_level_uv(j) PVPIN1A.367
PVPIN1A.368
C Calculate pressure values at the required two levels. PVPIN1A.369
zp3 =ak(iii) + bk(iii)*pstar_uv(j) PVPIN1A.370
zp4 =ak(iii+1) + bk(iii+1)*pstar_uv(j) PVPIN1A.371
PVPIN1A.372
C Calculate DRSU_DP and DRSV_DP. PVPIN1A.373
drsu_dp(j) = (rs_uv(j,iii+1)*u(j,iii+1)- PVPIN1A.374
& rs_uv(j,iii)*u(j,iii))/(zp4-zp3) PVPIN1A.375
drsv_dp(j) = (rs_uv(j,iii+1)*v(j,iii+1)- PVPIN1A.376
& rs_uv(j,iii)*v(j,iii))/(zp4-zp3) PVPIN1A.377
end if PVPIN1A.378
PVPIN1A.379
710 continue PVPIN1A.380
*IF DEF,MPP GSM3F403.950
! Update halos GSM3F403.951
CALL SWAPBOUNDS(theta_on_press,row_length,tot_P_ROWS, GSM3F403.952
& EW_Halo,NS_Halo,1) GSM3F403.953
CALL SWAPBOUNDS(rs_on_press,row_length,tot_P_ROWS, GSM3F403.954
& EW_Halo,NS_Halo,1) GSM3F403.955
CALL SWAPBOUNDS(u_on_press,row_length,tot_U_ROWS, GSM3F403.956
& EW_Halo,NS_Halo,1) GSM3F403.957
CALL SWAPBOUNDS(v_on_press,row_length,tot_U_ROWS, GSM3F403.958
& EW_Halo,NS_Halo,1) GSM3F403.959
CALL SWAPBOUNDS(drsu_dp,row_length,tot_U_ROWS, GSM3F403.960
& EW_Halo,NS_Halo,1) GSM3F403.961
CALL SWAPBOUNDS(drsv_dp,row_length,tot_U_ROWS, GSM3F403.962
& EW_Halo,NS_Halo,1) GSM3F403.963
CALL SWAPBOUNDS(dtheta_dp,row_length,tot_P_ROWS, GSM3F403.964
& EW_Halo,NS_Halo,1) GSM3F403.965
CALL SWAPBOUNDS(eta_level,row_length,tot_P_ROWS, GSM3F403.966
& EW_Halo,NS_Halo,1) GSM3F403.967
*ENDIF GSM3F403.968
PVPIN1A.381
return PVPIN1A.382
end PVPIN1A.383
PVPIN1A.384
*ENDIF PVPIN1A.385