*IF DEF,A15_1A THETPV1A.2
C ******************************COPYRIGHT****************************** GTS2F400.10243
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.10244
C GTS2F400.10245
C Use, duplication or disclosure of this code is subject to the GTS2F400.10246
C restrictions as set forth in the contract. GTS2F400.10247
C GTS2F400.10248
C Meteorological Office GTS2F400.10249
C London Road GTS2F400.10250
C BRACKNELL GTS2F400.10251
C Berkshire UK GTS2F400.10252
C RG12 2SZ GTS2F400.10253
C GTS2F400.10254
C If no contract has been raised with this copy of the code, the use, GTS2F400.10255
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.10256
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.10257
C Modelling at the above address. GTS2F400.10258
C ******************************COPYRIGHT****************************** GTS2F400.10259
C GTS2F400.10260
CLL Subroutine theta_pv ------------------------------------------------ THETPV1A.3
CLL THETPV1A.4
CLL Purpose: To compute Potential Temperature (Theta) on THETPV1A.5
CLL Potential Vorticity surfaces. THETPV1A.6
CLL Subroutine can cope with an array of desired pv surfaces. THETPV1A.7
CLL Includes a call to subroutine CALC_PV_P to calculate THETPV1A.8
CLL pv on pressure levels, and to output the array of THETPV1A.9
CLL values of Theta on pressure levels. THETPV1A.10
CLL This uses the Quasi-Hydrostatic equations, with complete THETPV1A.11
CLL representation of the Coriolis terms, and no metric THETPV1A.12
CLL terms omitted. THETPV1A.13
CLL The shallow atmosphere approximation is not made. THETPV1A.14
CLL Under UPDATE identifier GLOBAL, the data is THETPV1A.15
CLL assumed periodic along the rows. Note that because THETPV1A.16
CLL it is a diagnostic routine, care needs to be taken THETPV1A.17
CLL with missing data. THETPV1A.18
CLL THETPV1A.19
CLL Not suitable for single column use. THETPV1A.20
CLL THETPV1A.21
CLL Model Modification history: THETPV1A.22
CLL Version Date THETPV1A.23
CLL 3.1 21/01/93 Written by Simon Anderson. THETPV1A.24
CLL 3.1 18/01/93 New deck at the release of Version 3.1. THETPV1A.25
CLL 3.2 28/07/93 Change subroutine name to uppercase and array TS280793.13
CLL switch to switch1 for portability. Tracey Smith TS280793.14
CLL 3.4 27/05/94 Argument LLINTS added and passed to CALC_PV_P GSS1F304.216
CLL S.J.Swarbrick GSS1F304.217
CLL 4.3 21/03/97 MPP changes. S.D.Mullerworth GSM3F403.430
CLL TS280793.15
CLL THETPV1A.26
CLL Programming Standard: UM DOC Paper3, Version 4 (05/02/92) THETPV1A.27
CLL THETPV1A.28
CLL Logical Component Covered: D415 THETPV1A.29
CLL THETPV1A.30
CLL System Task: D4 THETPV1A.31
CLL THETPV1A.32
CLL Documentation: U.M.D.P No.13. Derivation and Calculation of THETPV1A.33
CLL Unified Model Potential Vorticity. THETPV1A.34
CLL By Simon Anderson and Ian Roulstone. THETPV1A.35
CLL THETPV1A.36
CLLEND THETPV1A.37
C THETPV1A.38
C*L ARGUMENTS: --------------------------------------------------------- THETPV1A.39
SUBROUTINE THETA_PV 1,5TS280793.16
1 (pstar,theta,u,v,p_field,u_field, THETPV1A.41
2 p_levels,row_length, THETPV1A.42
*CALL ARGFLDPT 
GSM3F403.431
3 rmdi,ak,bk,f3, THETPV1A.43
& e_levels,n_levels,dth_dph, TD141293.115
4 theta_pv_levs,des_pv,theta_pv_p_levs,des_p, THETPV1A.44
5 latitude_step_inverse,longitude_step_inverse, THETPV1A.45
6 cos_u_latitude,sec_p_latitude, THETPV1A.46
7 theta_on_pv,llints) GSS1F304.218
THETPV1A.48
THETPV1A.49
implicit none THETPV1A.50
logical llints GSS1F304.219
THETPV1A.51
C Input variables ------------------------------------------------------ THETPV1A.52
THETPV1A.53
integer THETPV1A.54
& p_field !IN Size of field on pressure points. THETPV1A.55
&,u_field !IN Size of field on velocity points. THETPV1A.56
&,p_levels !IN Number of pressure levels. THETPV1A.57
& ,n_levels !IN Number of half levels for dtheta/dp TD141293.116
&,row_length !IN Number of points in a row. THETPV1A.58
&,theta_pv_levs !IN Number of desired pv surfaces. THETPV1A.59
&,theta_pv_p_levs !IN Number of desired pressure levels. THETPV1A.60
THETPV1A.61
real THETPV1A.62
& pstar(p_field) !IN Surface pressure field. THETPV1A.63
&,u(u_field,p_levels) !IN Primary model array for u field. THETPV1A.64
&,v(u_field,p_levels) !IN Primary model array for v field. THETPV1A.65
&,theta(p_field,p_levels) !IN " " " theta field. THETPV1A.66
THETPV1A.67
real THETPV1A.68
& rmdi !IN Real missing data indicator. THETPV1A.69
&,ak(p_levels) !IN A coefficient of hybrid THETPV1A.70
& ! coordinates at full levels. THETPV1A.71
&,bk(p_levels) !IN B coefficient of hybrid THETPV1A.72
& ! coordinates at full levels. THETPV1A.73
&,f3(u_field) !IN Coriolis term. THETPV1A.74
& ,e_levels(n_levels) !IN half-levels over range TD141293.117
& ,dth_dph(p_field,n_levels) !IN dtheta/dp half-levels TD141293.118
&,des_pv(theta_pv_levs) !IN Value(s) of p.v. we want theta on. THETPV1A.75
&,des_p(theta_pv_p_levs) !IN Values of pressure we want pv on. THETPV1A.76
&,latitude_step_inverse !IN 1/latitude increment. THETPV1A.77
&,longitude_step_inverse !IN 1/longitude increment. THETPV1A.78
&,cos_u_latitude(u_field) !IN Cosine of latitude on uv-grid. THETPV1A.79
&,sec_p_latitude(p_field) !IN Secant of latitude on p-grid. THETPV1A.80
THETPV1A.81
THETPV1A.82
C Output variables ----------------------------------------------------- THETPV1A.83
THETPV1A.84
real THETPV1A.85
& theta_on_pv(p_field,theta_pv_levs) THETPV1A.86
& ! OUT Value of potential temperature THETPV1A.87
& ! on p.v. surface with THETPV1A.88
& ! p.v.=des_pv. THETPV1A.89
THETPV1A.90
C*---------------------------------------------------------------------- THETPV1A.91
C*L Workspace Usage: 4 arrays are required. THETPV1A.92
real THETPV1A.93
& pvort_p(p_field,theta_pv_p_levs) THETPV1A.94
& ! Calculated field of p.v. on pressure THETPV1A.95
& ! levels, from subroutine CALC_PV_P. THETPV1A.96
&,theta_on_press(p_field,theta_pv_p_levs) THETPV1A.97
& ! Calculated field of theta on pressure THETPV1A.98
& ! levels, from subroutine CALC_PV_P. THETPV1A.99
&,f3_p(p_field) ! Interpolated f3 field on p grid. THETPV1A.100
THETPV1A.101
integer THETPV1A.102
& switch1(p_field) ! Switch to determine hemispheric TS280793.17
& ! dependance. THETPV1A.104
C*---------------------------------------------------------------------- THETPV1A.105
C*L External subroutine calls: THETPV1A.106
external calc_pv_p ! Compute p.v. on pressure levels. THETPV1A.107
external th_pvint ! Interpolate theta onto pv surfaces. THETPV1A.108
external uv_to_p ! Interpolate u-grid field to p-grid fld. THETPV1A.109
THETPV1A.110
C*---------------------------------------------------------------------- THETPV1A.111
C*L Define local variable. THETPV1A.112
integer i,k ! Loop variable. THETPV1A.113
real mn ! Mean value used in computing pole values. THETPV1A.114
*CALL TYPFLDPT GSM3F403.432
*IF DEF,MPP GSM3F403.433
integer info GSM3F403.434
*ENDIF GSM3F403.435
THETPV1A.116
C ---------------------------------------------------------------------- THETPV1A.117
CL Section 1 Compute p.v. on each of the desired pressure levels. THETPV1A.118
CL ~~~~~~~~~ Set switch array to check for sign of desired pv values. THETPV1A.119
C ---------------------------------------------------------------------- THETPV1A.120
THETPV1A.121
C Loop from 1 to theta_pv_p_levs. THETPV1A.122
THETPV1A.123
do 100 k=1,theta_pv_p_levs THETPV1A.124
THETPV1A.125
C The array des_p is assumed to be ordered THETPV1A.126
C from the bottom of the atmosphere to the top. THETPV1A.127
THETPV1A.128
call calc_pv_p THETPV1A.129
1 (pstar,theta,u,v,p_field,u_field, THETPV1A.130
2 p_levels,row_length, THETPV1A.131
*CALL ARGFLDPT GSM3F403.436
3 rmdi,ak,bk,des_p(k),f3, THETPV1A.132
& e_levels,n_levels,dth_dph, TD141293.119
4 latitude_step_inverse,longitude_step_inverse, THETPV1A.133
5 cos_u_latitude,sec_p_latitude, THETPV1A.134
6 pvort_p(1,k),theta_on_press(1,k),llints) GSS1F304.220
THETPV1A.136
100 continue THETPV1A.137
THETPV1A.138
C Find f3_p on p-points. THETPV1A.139
*IF DEF,MPP GSM3F403.437
C MPP: One fewer row than normal in P_field. GSM3F403.438
call uv_to_p GSM3F403.439
& (f3,f3_p(row_length+1),u_field,p_field-row_length, GSM3F403.440
& row_length,u_field/row_length) GSM3F403.441
CALL SWAPBOUNDS(f3_p,row_length,tot_P_ROWS, GSM3F403.442
& EW_Halo,NS_Halo,1) GSM3F403.443
*ELSE GSM3F403.444
call uv_to_p THETPV1A.140
& (f3,f3_p(row_length+1),u_field,p_field, THETPV1A.141
& row_length,u_field/row_length) THETPV1A.142
*ENDIF GSM3F403.445
THETPV1A.143
C Calculate f3_p at the poles. THETPV1A.144
mn = 0. THETPV1A.145
*IF DEF,MPP GSM3F403.446
IF (at_top_of_LPG) THEN GSM3F403.447
CALL GCG_RVECSUMR(ROW_LENGTH,ROW_LENGTH-2*EW_Halo, GSM3F403.448
& EW_Halo+TOP_ROW_START,1,f3,GC_ROW_GROUP,info,mn) GSM3F403.449
mn=mn/GLOBAL_ROW_LENGTH GSM3F403.450
do i=TOP_ROW_START,START_POINT_NO_HALO-1 GSM3F403.451
f3_p(i) = mn GSM3F403.452
end do GSM3F403.453
ENDIF GSM3F403.454
*ELSE GSM3F403.455
do i=1,row_length THETPV1A.146
mn = mn + f3(i) THETPV1A.147
end do THETPV1A.148
mn = mn/row_length THETPV1A.149
GSM3F403.456
do i=1,row_length THETPV1A.150
f3_p(i) = mn THETPV1A.151
end do THETPV1A.152
*ENDIF GSM3F403.457
THETPV1A.153
mn = 0. THETPV1A.154
*IF DEF,MPP GSM3F403.458
IF (at_base_of_LPG) THEN GSM3F403.459
CALL GCG_RVECSUMR(ROW_LENGTH,ROW_LENGTH-2*EW_Halo, GSM3F403.460
& EW_Halo,1,f3(U_BOT_ROW_START),GC_ROW_GROUP,info,mn) GSM3F403.461
mn=mn/GLOBAL_ROW_LENGTH GSM3F403.462
do i=P_BOT_ROW_START,LAST_P_FLD_PT-1 GSM3F403.463
f3_p(i) = mn GSM3F403.464
end do GSM3F403.465
ENDIF GSM3F403.466
*ELSE GSM3F403.467
do i=u_field-row_length+1,u_field THETPV1A.155
mn = mn + f3(i) THETPV1A.156
end do THETPV1A.157
mn = mn/row_length THETPV1A.158
do i=p_field-row_length+1,p_field THETPV1A.159
f3_p(i) = mn THETPV1A.160
end do THETPV1A.161
*ENDIF GSM3F403.468
THETPV1A.162
C Set hemispheric switch array. THETPV1A.163
do 101 i=1,p_field THETPV1A.164
if (f3_p(i).ge.0.) then THETPV1A.165
switch1(i) = 1 TS280793.18
else THETPV1A.167
switch1(i) = -1 TS280793.19
endif THETPV1A.169
101 continue THETPV1A.170
THETPV1A.171
C Change sign of pv field in Southern hemisphere. This is required THETPV1A.172
C since user asks for values of Des_pv with Northern hemisphere sign, THETPV1A.173
C namely positive. We need to either change the sign of this Des_pv THETPV1A.174
C value depending on the hemispere we are in, or, as we do here, THETPV1A.175
C change the sign of the potential vorticity instead. THETPV1A.176
THETPV1A.177
do 102 k=1,theta_pv_p_levs THETPV1A.178
do 103 i=1,p_field THETPV1A.179
if (pvort_p(i,k).ne.rmdi) then THETPV1A.180
pvort_p(i,k) = pvort_p(i,k)*switch1(i) TS280793.20
endif THETPV1A.182
103 continue THETPV1A.183
102 continue THETPV1A.184
THETPV1A.185
C ---------------------------------------------------------------------- THETPV1A.186
CL Section 2 Compute theta on each pv surface using potential vorticity THETPV1A.187
CL ~~~~~~~~~ on pressure levels, and theta on pressure levels. THETPV1A.188
CL Note that all the missing data points at limited area THETPV1A.189
CL model boundaries, and polar rows are taken care of THETPV1A.190
CL in the CALC_PV_P subroutine. THETPV1A.191
C ---------------------------------------------------------------------- THETPV1A.192
THETPV1A.193
do 200 k=1,theta_pv_levs THETPV1A.194
THETPV1A.195
call th_pvint THETPV1A.196
1 (theta_on_press,pvort_p, THETPV1A.197
2 p_field,theta_pv_p_levs,rmdi,des_pv(k), THETPV1A.198
3 theta_on_pv(1,k)) THETPV1A.199
THETPV1A.200
200 continue THETPV1A.201
THETPV1A.202
return THETPV1A.203
end THETPV1A.204
THETPV1A.205
*ENDIF THETPV1A.206