*IF DEF,A15_1A THPVIN1A.2
C ******************************COPYRIGHT****************************** GTS2F400.10279
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.10280
C GTS2F400.10281
C Use, duplication or disclosure of this code is subject to the GTS2F400.10282
C restrictions as set forth in the contract. GTS2F400.10283
C GTS2F400.10284
C Meteorological Office GTS2F400.10285
C London Road GTS2F400.10286
C BRACKNELL GTS2F400.10287
C Berkshire UK GTS2F400.10288
C RG12 2SZ GTS2F400.10289
C GTS2F400.10290
C If no contract has been raised with this copy of the code, the use, GTS2F400.10291
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.10292
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.10293
C Modelling at the above address. GTS2F400.10294
C ******************************COPYRIGHT****************************** GTS2F400.10295
C GTS2F400.10296
CLL SUBROUTINE th_pvint ----------------------------------------------- THPVIN1A.3
CLL THPVIN1A.4
CLL Purpose: Interpolates fields onto a desired pv surface. THPVIN1A.5
CLL Note that routine assumes that pv is monotonic. THPVIN1A.6
CLL THPVIN1A.7
CLL Suitable for single column use. THPVIN1A.8
CLL THPVIN1A.9
CLL Model Modification history: THPVIN1A.10
CLL Version Date THPVIN1A.11
CLL 3.1 21/01/93 Written By Simon Anderson. THPVIN1A.12
CLL 3.1 18/01/93 New deck at the release of Version 3.1. THPVIN1A.13
CLL 3.2 28/07/93 Change subroutine name to uppercase for TS280793.21
CLL portability. Tracey Smith TS280793.22
CLL 3.3 14/12/93 Change to use linear interpolation only TD141293.1
CLL in vertical. Terry Davies TD141293.2
CLL THPVIN1A.14
CLL Programming standard UM DOC Paper 3, Version 4(05/02/92) A THPVIN1A.15
CLL THPVIN1A.16
CLL Logical Component Covered: D415 THPVIN1A.17
CLL THPVIN1A.18
CLL Project Task: D4 THPVIN1A.19
CLL THPVIN1A.20
CLL Documentation: U.M.D.P No.13. Derivation and Calculation of THPVIN1A.21
CLL Unified Model Potential Vorticity. THPVIN1A.22
CLL By Simon Anderson and Ian Roulstone. THPVIN1A.23
CLL THPVIN1A.24
CLLEND------------------------------------------------------------------ THPVIN1A.25
THPVIN1A.26
C*L ARGUMENTS:---------------------------------------------------------- THPVIN1A.27
SUBROUTINE TH_PVINT 1TS280793.23
1 (theta_on_press,pvort_p, THPVIN1A.29
2 p_field,theta_pv_p_levs,rmdi,des_pv, THPVIN1A.30
3 theta_on_pv) THPVIN1A.31
THPVIN1A.32
implicit none THPVIN1A.33
THPVIN1A.34
C Input variables ------------------------------------------------------ THPVIN1A.35
THPVIN1A.36
integer THPVIN1A.37
& p_field !IN Points in horizontal p field. THPVIN1A.38
&,theta_pv_p_levs !IN Number of desired pressure levels. THPVIN1A.39
THPVIN1A.40
real THPVIN1A.41
& theta_on_press(p_field,theta_pv_p_levs) THPVIN1A.42
& !IN Calculated field of THPVIN1A.43
& ! theta on pressure levels. THPVIN1A.44
&,pvort_p(p_field,theta_pv_p_levs) THPVIN1A.45
& !IN Calculated field of potential THPVIN1A.46
& ! vorticity on pressure levels. THPVIN1A.47
THPVIN1A.48
real THPVIN1A.49
& rmdi !IN Real missing data indicator. THPVIN1A.50
&,des_pv !IN Desired pv surface in pv units we THPVIN1A.51
& ! want variables interpolated onto. THPVIN1A.52
THPVIN1A.53
THPVIN1A.54
C Output variables ----------------------------------------------------- THPVIN1A.55
THPVIN1A.56
real THPVIN1A.57
& theta_on_pv(p_field) !OUT Interpolated theta field on THPVIN1A.58
& ! pv=des_pv surface. THPVIN1A.59
THPVIN1A.60
C*---------------------------------------------------------------------- THPVIN1A.61
C*L Workspace usage: THPVIN1A.62
integer THPVIN1A.63
& base_level_p(p_field) ! The pressure level below the desired pv THPVIN1A.64
& ! surface for each atmospheric column THPVIN1A.65
& ! (set to theta_pv_p_levs if Des_pv THPVIN1A.66
& ! is above the top level, or if THPVIN1A.67
& ! level is not found). THPVIN1A.68
& ! Calculated at p-points. THPVIN1A.69
real THPVIN1A.70
& des_pv_mks ! Desired pv surface in mks units we THPVIN1A.71
& ! want variables interpolated onto. THPVIN1A.72
THPVIN1A.76
C*---------------------------------------------------------------------- THPVIN1A.77
C*L External subroutine calls: None. THPVIN1A.78
THPVIN1A.79
C*---------------------------------------------------------------------- THPVIN1A.80
C*L Local variables: THPVIN1A.81
integer i,j,iii,level ! Loop counts. THPVIN1A.82
THPVIN1A.83
real zthp1,zthp2,zpv1,zpv2 TD141293.3
THPVIN1A.87
logical l_down ! Logical used to control Do While. THPVIN1A.88
THPVIN1A.89
C ---------------------------------------------------------------------- THPVIN1A.90
CL Section 1 : Find the value of the base level. This is the largest THPVIN1A.91
CL ~~~~~~~~~ value of the pressure level such that pv(level) < Des_pv, THPVIN1A.92
CL calculated on p-points. THPVIN1A.93
C ---------------------------------------------------------------------- THPVIN1A.94
THPVIN1A.95
C Firstly, we must divide the desired pv field by 1,000,000 in THPVIN1A.96
C order to get actual values of potential vorticity. THPVIN1A.97
THPVIN1A.98
des_pv_mks = des_pv * 1.0E-06 THPVIN1A.99
THPVIN1A.100
do 110 i = 1,p_field THPVIN1A.101
base_level_p(i) = theta_pv_p_levs THPVIN1A.102
110 continue THPVIN1A.104
THPVIN1A.105
C Find first level, searching down from above, above which Des_pv THPVIN1A.106
C value can be found. Set to top level if none found. THPVIN1A.107
THPVIN1A.108
do 120 i=1,p_field THPVIN1A.109
l_down = .true. THPVIN1A.110
level = theta_pv_p_levs THPVIN1A.111
do while (l_down) THPVIN1A.112
level = level - 1 THPVIN1A.113
if(pvort_p(i,level).ne.rmdi .and. THPVIN1A.114
& pvort_p(i,level+1).ne.rmdi ) then THPVIN1A.115
if(pvort_p(i,level).lt.des_pv_mks.and. THPVIN1A.116
& pvort_p(i,level+1).gt.des_pv_mks ) then THPVIN1A.117
base_level_p(i) = level THPVIN1A.118
l_down = .false. THPVIN1A.119
end if THPVIN1A.120
end if THPVIN1A.121
if (level.le.1) l_down = .false. THPVIN1A.122
end do THPVIN1A.123
120 continue THPVIN1A.124
THPVIN1A.125
C When this loop is done, base_level_p will be the pressure level THPVIN1A.126
C with the highest value of pv LESS than the desired pv value. THPVIN1A.127
C If for any reason, the desired pv value does not lie within any two THPVIN1A.128
C desired pressure levels then base_level_p is set to theta_pv_p_levs. THPVIN1A.129
THPVIN1A.130
THPVIN1A.131
C ---------------------------------------------------------------------- THPVIN1A.132
CL Section 2 : This section will interpolate variables held on THPVIN1A.133
CL ~~~~~~~~~ the p-grid onto the desired potential vorticity surface. THPVIN1A.134
CL Used for THETA_ON_PV at each point. THPVIN1A.135
C ---------------------------------------------------------------------- THPVIN1A.136
THPVIN1A.137
C---- Given Theta as a function of pv and Des_pv lying between THPVIN1A.138
C---- Zpv1 and Zpv2, calculate theta at DEs_pv by fitting TD141293.4
C---- linear Lagrange polynomial and evaluating TD141293.5
C---- Potential Vorticity = Des_pv. THPVIN1A.141
C---- THPVIN1A.142
C---- If surface is set to theta_pv_p_levs, then we set the value THPVIN1A.143
C---- to missing data. THPVIN1A.144
C---- THPVIN1A.145
do 210 i=1,p_field THPVIN1A.151
THPVIN1A.152
if(base_level_p(i).eq.theta_pv_p_levs) then THPVIN1A.153
theta_on_pv(i)=rmdi THPVIN1A.154
else THPVIN1A.155
THPVIN1A.156
iii=base_level_p(i) THPVIN1A.209
zpv1=pvort_p(i,iii) TD141293.6
zthp1=theta_on_press(i,iii) TD141293.7
zpv2=pvort_p(i,iii+1) TD141293.8
zthp2=theta_on_press(i,iii+1) TD141293.9
THPVIN1A.214
theta_on_pv(i)=(des_pv_mks-zpv2)*zthp1/(zpv1-zpv2)+ TD141293.10
& (des_pv_mks-zpv1)*zthp2/(zpv2-zpv1) TD141293.11
THPVIN1A.223
end if THPVIN1A.224
THPVIN1A.225
210 continue THPVIN1A.226
THPVIN1A.227
return THPVIN1A.228
end THPVIN1A.229
THPVIN1A.230
*ENDIF THPVIN1A.231