*IF DEF,A15_1A CALCP21A.2
C ******************************COPYRIGHT****************************** GTS2F400.703
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.704
C GTS2F400.705
C Use, duplication or disclosure of this code is subject to the GTS2F400.706
C restrictions as set forth in the contract. GTS2F400.707
C GTS2F400.708
C Meteorological Office GTS2F400.709
C London Road GTS2F400.710
C BRACKNELL GTS2F400.711
C Berkshire UK GTS2F400.712
C RG12 2SZ GTS2F400.713
C GTS2F400.714
C If no contract has been raised with this copy of the code, the use, GTS2F400.715
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.716
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.717
C Modelling at the above address. GTS2F400.718
C ******************************COPYRIGHT****************************** GTS2F400.719
C GTS2F400.720
CLL Subroutine calc_pv_p ----------------------------------------------- CALCP21A.3
CLL CALCP21A.4
CLL Purpose: To compute Ertel potential vorticity CALCP21A.5
CLL on pressure levels. CALCP21A.6
CLL Uses the Quasi-Hydrostatic equations, with complete CALCP21A.7
CLL representation of the Coriolis terms, and no metric CALCP21A.8
CLL terms omitted. CALCP21A.9
CLL The shallow atmosphere approximation is not made. CALCP21A.10
CLL Under UPDATE identifier GLOBAL, the data is CALCP21A.11
CLL assumed periodic along the rows. Note that because CALCP21A.12
CLL it is a diagnostic routine care needs to be taken CALCP21A.13
CLL with missing data. CALCP21A.14
CLL CALCP21A.15
CLL Not suitable for single column use. CALCP21A.16
CLL CALCP21A.17
CLL Model Modification history: CALCP21A.18
CLL Version Date CALCP21A.19
CLL 3.1 12/11/92 Written by Simon Anderson. CALCP21A.20
CLL 3.1 18/01/93 New deck at the release of Version 3.1. CALCP21A.21
CLL 3.2 28/07/93 Change subroutine name to uppercase for TS280793.4
CLL portability. Tracey Smith TS280793.5
CLL 3.4 26/05/94 Argument llints added and passed to calc_rs GSS1F304.150
CLL S.J.Swarbrick GSS1F304.151
CLL CALCP21A.22
CLL Programming Standard: UM DOC Paper3, Version 4 (05/02/92) CALCP21A.23
CLL CALCP21A.24
CLL Logical Component Covered: D415 CALCP21A.25
CLL CALCP21A.26
CLL System Task: D4 CALCP21A.27
CLL CALCP21A.28
CLL Documentation: U.M.D.P No.13. Derivation and Calculation of CALCP21A.29
CLL Unified Model Potential Vorticity. CALCP21A.30
CLL By Simon Anderson and Ian Roulstone. CALCP21A.31
CLL CALCP21A.32
CLLEND CALCP21A.33
C CALCP21A.34
C*L ARGUMENTS: --------------------------------------------------------- CALCP21A.35
SUBROUTINE CALC_PV_P 2,9TS280793.6
1 (pstar,theta,u,v,p_field,u_field, CALCP21A.37
2 p_levels,row_length, CALCP21A.38
*CALL ARGFLDPT
GSM3F403.712
3 rmdi,ak,bk,des_press,f3, CALCP21A.39
& e_levels,n_levels,dth_dph, TD141293.175
4 latitude_step_inverse,longitude_step_inverse, CALCP21A.40
5 cos_u_latitude,sec_p_latitude, CALCP21A.41
6 pvort_p,theta_on_press,llints) GSS1F304.152
CALCP21A.43
CALCP21A.44
implicit none CALCP21A.45
logical llints GSS1F304.153
CALCP21A.46
C Input variables ------------------------------------------------------ CALCP21A.47
CALCP21A.48
integer CALCP21A.49
& p_field !IN Size of field on pressure points. CALCP21A.50
&,u_field !IN Size of field on velocity points. CALCP21A.51
&,p_levels !IN Number of pressure levels. CALCP21A.52
&,row_length !IN Number of points in a row. CALCP21A.53
& ,n_levels !IN Number of levels of spline TD141293.176
CALCP21A.54
*CALL TYPFLDPT
GSM3F403.713
GSM3F403.714
real CALCP21A.55
& pstar(p_field) !IN Surface pressure field. CALCP21A.56
&,u(u_field,p_levels) !IN Primary model array for u field. CALCP21A.57
&,v(u_field,p_levels) !IN Primary model array for v field. CALCP21A.58
&,theta(p_field,p_levels) !IN " " " theta field. CALCP21A.59
& ,dth_dph(p_field,n_levels) !IN dtheta/dp half-levels TD141293.177
CALCP21A.60
real CALCP21A.61
& rmdi !IN Real missing data indicator. CALCP21A.62
&,ak(p_levels) !IN A coefficient of hybrid CALCP21A.63
& ! coordinates at full levels. CALCP21A.64
&,bk(p_levels) !IN B coefficient of hybrid CALCP21A.65
& ! coordinates at full levels. CALCP21A.66
& ,e_levels(n_levels) !IN half-levels over range TD141293.178
&,des_press !IN Value of pressure we want pv on. CALCP21A.67
&,f3(u_field) !IN Coriolis term. CALCP21A.68
&,latitude_step_inverse !IN 1/latitude increment. CALCP21A.69
&,longitude_step_inverse !IN 1/longitude increment. CALCP21A.70
&,cos_u_latitude(u_field) !IN Cosine of latitude on u field. CALCP21A.71
&,sec_p_latitude(p_field) !IN Secant of latitude on p field. CALCP21A.72
CALCP21A.73
CALCP21A.74
C Output variables ----------------------------------------------------- CALCP21A.75
CALCP21A.76
real CALCP21A.77
& pvort_p(p_field) ! OUT Value of potential vorticity CALCP21A.78
& ! on pressure level with CALCP21A.79
& ! pressure=des_press. CALCP21A.80
&,theta_on_press(p_field) ! OUT Value of theta on pressure CALCP21A.81
& ! level with pressure=des_press. CALCP21A.82
CALCP21A.83
C*---------------------------------------------------------------------- CALCP21A.84
C*L Workspace Usage: 14 arrays are required. CALCP21A.85
logical CALCP21A.86
& mask(p_field) ! Logical mask used to make CALCP21A.87
& ! vectorising of a loop possible. CALCP21A.88
real CALCP21A.89
& rs(p_field,p_levels) ! Calculated pseudo radius of Earth. CALCP21A.90
&,u_on_press(u_field) ! Interpolated value of u on p level. CALCP21A.91
&,v_on_press(u_field) ! Interpolated value of v on p level. CALCP21A.92
&,rs_on_press(p_field) ! Interpolated value of rs on p level. CALCP21A.93
&,drsu_dp(p_field) ! D(rs.u)/D(p) on pressure level. CALCP21A.94
&,drsv_dp(p_field) ! D(rs.v)/D(p) on pressure level. CALCP21A.95
&,dtheta_dp(p_field) ! D(Theta)/D(p) on pressure level. CALCP21A.96
&,dtheta_dlatitude(p_field) ! D(Theta)/D(lat) on pressure level. CALCP21A.97
&,vorticity3(p_field) ! A calculated term in the pv equation. CALCP21A.98
&,vorticity4(p_field) ! A calculated term in the pv equation. CALCP21A.99
&,vorticity5(p_field) ! A calculated term in the pv equation. CALCP21A.100
&,f3_p(p_field) ! Interpolated f3 field on p field. CALCP21A.101
&,eta_level(p_field) ! Interpolated eta-value of theta level TD141293.179
&,work1(p_field) ! General workspace. CALCP21A.102
CALCP21A.103
C*---------------------------------------------------------------------- CALCP21A.104
C*L External subroutine calls: CALCP21A.105
external calc_rs ! Compute pseudo radius. CALCP21A.106
external pv_pint ! Interpolate variables to theta level. CALCP21A.107
external vortic3 ! Compute term 1. CALCP21A.108
external vortic4 ! Compute term 2. CALCP21A.109
external vortic5 ! Compute term 3. CALCP21A.110
external uv_to_p ! Interpolate u-grid field to p-grid fld. CALCP21A.111
CALCP21A.112
C*---------------------------------------------------------------------- CALCP21A.113
C*L Call comdecks to get required variables: CALCP21A.114
*CALL C_G
CALCP21A.115
*CALL C_OMEGA
CALCP21A.116
CALCP21A.117
C*---------------------------------------------------------------------- CALCP21A.118
C*L Define local variables. CALCP21A.119
integer i,j,k ! Loop variables. CALCP21A.120
*IF DEF,MPP GSM3F403.715
*IF -DEF,GLOBAL GSM3F403.716
integer START,END GSM3F403.717
*ENDIF GSM3F403.718
INTEGER info !GCG return code GSM3F403.719
*ENDIF GSM3F403.720
real mn ! Mean value used in computing pole values. CALCP21A.121
CALCP21A.122
CALCP21A.123
C ---------------------------------------------------------------------- CALCP21A.124
CL Section 1 Compute rs, the pseudo radius. CALCP21A.125
CL ~~~~~~~~~ CALCP21A.126
C ---------------------------------------------------------------------- CALCP21A.127
CALCP21A.128
CL Section 1.1 Call CALC_RS to get rs for level 1. CALCP21A.129
CL ~~~~~~~~~~~ Rs is returned in rs(1,k). CALCP21A.130
CL Ts is returned in work1. Rs at level k-1 is input as CALCP21A.131
CL rs(1,2) as at k-1=0 the input is not used by calc_rs. CALCP21A.132
CALCP21A.133
k=1 CALCP21A.134
call calc_rs
CALCP21A.135
1 (pstar,ak,bk,work1,rs(1,2),rs(1,k),p_field,k,p_levels, GSS1F304.154
2 llints) GSS1F304.155
CALCP21A.137
CL Section 1.2 Call CALC_RS to get rs for level k. CALCP21A.138
CL ~~~~~~~~~~~ Rs is returned in rs(1,k). CALCP21A.139
CL Ts is returned in work1. Rs at level k-1 is input as CALCP21A.140
CL rs(1,k-1). CALCP21A.141
CL Loop from 2 to p_levels. CALCP21A.142
CALCP21A.143
do 100 k=2,p_levels CALCP21A.144
CALCP21A.145
call calc_rs
CALCP21A.146
1 (pstar,ak,bk,work1,rs(1,k-1),rs(1,k),p_field,k,p_levels, GSS1F304.156
2 llints) GSS1F304.157
100 continue CALCP21A.149
CALCP21A.150
CALCP21A.151
C ---------------------------------------------------------------------- CALCP21A.152
CL Section 2 Interpolate p, u, v and rs onto desired theta level, CALCP21A.153
CL ~~~~~~~~~ and compute Dtheta/Dp, D(rs.u)/D(p) and D(rs.v)/D(p) CALCP21A.154
CL assuming cubic variation where possible. CALCP21A.155
CL Interpolate D(rs.u)/D(p) and D(rs.v)/D(p) to the p-grid. CALCP21A.156
C ---------------------------------------------------------------------- CALCP21A.157
CALCP21A.158
call pv_pint
CALCP21A.159
1 (pstar,theta,rs,u,v,p_field,u_field, CALCP21A.160
2 p_levels,row_length, CALCP21A.161
*CALL ARGFLDPT
GSM3F403.721
3 rmdi,ak,bk,des_press, CALCP21A.162
& eta_level,e_levels,dth_dph,n_levels, TD141293.180
4 theta_on_press,rs_on_press, CALCP21A.163
5 u_on_press,v_on_press, CALCP21A.164
6 drsu_dp,drsv_dp,dtheta_dp) CALCP21A.165
CALCP21A.166
CL Section 2.1 Interpolate D(rs.u)/D(p) and D(rs.v)/D(p) to the p-grid. CALCP21A.167
CL ~~~~~~~~~~~ CALCP21A.168
CALCP21A.169
do i=1,row_length GSM3F403.722
work1(i)=drsu_dp(i) GSM3F403.723
enddo GSM3F403.724
call uv_to_p
CALCP21A.170
1 (drsu_dp,work1(row_length+1),u_field,p_field, GSM3F403.725
*IF DEF,MPP GSM3F403.726
2 row_length,U_LAST_ROW+1) GSM3F403.727
*ELSE GSM3F403.728
2 row_length,u_field/row_length) CALCP21A.172
*ENDIF GSM3F403.729
do i=1,last_p_fld_pt GSM3F403.730
drsu_dp(i)=work1(i) GSM3F403.731
enddo GSM3F403.732
CALCP21A.173
do i=1,row_length GSM3F403.733
work1(i)=drsv_dp(i) GSM3F403.734
enddo GSM3F403.735
call uv_to_p
CALCP21A.174
1 (drsv_dp,work1(row_length+1),u_field,p_field, GSM3F403.736
*IF DEF,MPP GSM3F403.737
2 row_length,U_LAST_ROW+1) GSM3F403.738
*ELSE GSM3F403.739
2 row_length,u_field/row_length) CALCP21A.176
*ENDIF GSM3F403.740
do i=1,last_p_fld_pt GSM3F403.741
drsv_dp(i)=work1(i) GSM3F403.742
enddo GSM3F403.743
CALCP21A.178
C ---------------------------------------------------------------------- CALCP21A.179
CL Section 3 Compute the various terms in the potential vorticity CALCP21A.180
CL ~~~~~~~~~ equation for desired theta surface. CALCP21A.181
C ---------------------------------------------------------------------- CALCP21A.182
CALCP21A.183
CL Section 3.1 Compute 1/(rs*rs*cos(phi))*D(rs.v)/D(p)*D(theta)/D(lamda) CALCP21A.184
CL ~~~~~~~~~~~ - 1/(rs*rs)*D(rs.u)/D(p)*D(theta)/D(phi). CALCP21A.185
call vortic3
CALCP21A.186
1 (rs_on_press,theta_on_press,drsu_dp,drsv_dp, CALCP21A.187
2 sec_p_latitude, CALCP21A.188
3 vorticity3,dtheta_dlatitude, CALCP21A.189
4 p_field,row_length, CALCP21A.190
*CALL ARGFLDPT
GSM3F403.744
5 latitude_step_inverse,longitude_step_inverse) CALCP21A.191
CALCP21A.192
CALCP21A.193
CL Section 3.2 Compute 1/rs*2*omega*cos(phi) * D(theta)/D(phi). CALCP21A.194
CL ~~~~~~~~~~~ CALCP21A.195
call vortic4
CALCP21A.196
1 (rs_on_press,theta_on_press,dtheta_dlatitude, CALCP21A.197
2 sec_p_latitude, CALCP21A.198
3 vorticity4, CALCP21A.199
4 p_field,row_length, CALCP21A.200
*CALL ARGFLDPT
GSM3F403.745
5 latitude_step_inverse,longitude_step_inverse) CALCP21A.201
CALCP21A.202
CL Section 3.3 Compute 1/(rs*cos(phi))* CALCP21A.203
CL ~~~~~~~~~~~ (D(v)/D(lambda)-D(u.cos(phi))/D(phi)). CALCP21A.204
call vortic5
CALCP21A.205
1 (u_on_press,v_on_press,rs_on_press, CALCP21A.206
2 cos_u_latitude,sec_p_latitude, CALCP21A.207
3 vorticity5, CALCP21A.208
4 p_field,u_field,row_length, CALCP21A.209
*CALL ARGFLDPT
GSM3F403.746
5 latitude_step_inverse,longitude_step_inverse) CALCP21A.210
CALCP21A.211
CALCP21A.212
C ---------------------------------------------------------------------- CALCP21A.213
CL Section 4 Compute the potential vorticity using the vorticity terms, CALCP21A.214
CL ~~~~~~~~~ D(theta)/D(p) and the value of the coriolis term. CALCP21A.215
CL Care needs to be taken for rmdi. CALCP21A.216
C ---------------------------------------------------------------------- CALCP21A.217
CALCP21A.218
CL Section 4.1 Firstly, find which points we can actually calculate CALCP21A.219
CL ~~~~~~~~~~~ potential vorticity on. To do this we test for missing CALCP21A.220
CL data at all points in a 3x3 stencil around each point. CALCP21A.221
CL Special care is taken at the boundaries and poles. CALCP21A.222
CALCP21A.223
C Set logical mask default to .TRUE. CALCP21A.224
do i = 1,p_field CALCP21A.225
mask(i) = .true. CALCP21A.226
end do CALCP21A.227
CALCP21A.228
C Check points on non-polar rows: CALCP21A.229
*IF DEF,MPP GSM3F403.747
do j=FIRST_ROW,P_LAST_ROW GSM3F403.748
*ELSE GSM3F403.749
do j=2,p_field/row_length-1 CALCP21A.230
*ENDIF GSM3F403.750
C First point on each row: CALCP21A.231
i=(j-1)*row_length+1 CALCP21A.232
*IF DEF,GLOBAL CALCP21A.233
*IF -DEF,MPP GSM3F403.751
if ( eta_level(i-1 ).eq.rmdi. TD141293.181
& or.eta_level(i-row_length ).eq.rmdi. TD141293.182
& or.eta_level(i-row_length+1).eq.rmdi. TD141293.183
& or.eta_level(i+row_length-1).eq.rmdi. TD141293.184
& or.eta_level(i ).eq.rmdi. TD141293.185
& or.eta_level(i+1 ).eq.rmdi. TD141293.186
& or.eta_level(i+2*row_length-1).eq.rmdi. TD141293.187
& or.eta_level(i+row_length ).eq.rmdi. TD141293.188
& or.eta_level(i+row_length+1).eq.rmdi) mask(i)=.false. TD141293.189
*ENDIF GSM3F403.752
*ELSE GSM3F403.753
*IF DEF,MPP GSM3F403.754
IF (at_left_of_LPG) THEN GSM3F403.755
mask(i+EW_Halo)=.false. GSM3F403.756
ENDIF GSM3F403.757
*ELSE CALCP21A.243
mask(i)=.false. CALCP21A.244
*ENDIF CALCP21A.245
*ENDIF GSM3F403.758
C Inner points on each row: CALCP21A.246
*IF DEF,MPP GSM3F403.759
*IF -DEF,GLOBAL GSM3F403.760
C Borders of area done separately GSM3F403.761
START=EW_Halo+1 GSM3F403.762
END=EW_Halo GSM3F403.763
IF (at_left_of_LPG) START=START+1 GSM3F403.764
IF (at_right_of_LPG) END=END+1 GSM3F403.765
do i=(j-1)*row_length+START,(j-1)*row_length+row_length-END GSM3F403.766
*ELSE GSM3F403.767
do i=(j-1)*row_length+1+EW_Halo, GSM3F403.768
& (j-1)*row_length+row_length-EW_Halo GSM3F403.769
*ENDIF GSM3F403.770
*ELSE GSM3F403.771
do i=(j-1)*row_length+2,(j-1)*row_length+row_length-1 CALCP21A.247
*ENDIF GSM3F403.772
if ( eta_level(i-row_length-1).eq.rmdi. TD141293.190
& or.eta_level(i-row_length ).eq.rmdi. TD141293.191
& or.eta_level(i-row_length+1).eq.rmdi. TD141293.192
& or.eta_level(i-1 ).eq.rmdi. TD141293.193
& or.eta_level(i ).eq.rmdi. TD141293.194
& or.eta_level(i+1 ).eq.rmdi. TD141293.195
& or.eta_level(i+row_length-1).eq.rmdi. TD141293.196
& or.eta_level(i+row_length ).eq.rmdi. TD141293.197
& or.eta_level(i+row_length+1).eq.rmdi) mask(i)=.false. TD141293.198
end do CALCP21A.257
C Last point on each row: CALCP21A.258
i=j*row_length CALCP21A.259
*IF DEF,GLOBAL CALCP21A.260
*IF -DEF,MPP GSM3F403.773
if ( eta_level(i-row_length-1).eq.rmdi. TD141293.199
& or.eta_level(i-row_length ).eq.rmdi. TD141293.200
& or.eta_level(i-2*row_length+1).eq.rmdi. TD141293.201
& or.eta_level(i-1 ).eq.rmdi. TD141293.202
& or.eta_level(i ).eq.rmdi. TD141293.203
& or.eta_level(i-row_length+1).eq.rmdi. TD141293.204
& or.eta_level(i+row_length-1).eq.rmdi. TD141293.205
& or.eta_level(i+row_length ).eq.rmdi. TD141293.206
& or.eta_level(i+1 ).eq.rmdi) mask(i)=.false. TD141293.207
*ENDIF GSM3F403.774
*ELSE GSM3F403.775
*IF DEF,MPP GSM3F403.776
IF (at_right_of_LPG) THEN GSM3F403.777
mask(i-EW_Halo)=.false. GSM3F403.778
ENDIF GSM3F403.779
*ELSE CALCP21A.270
mask(i)=.false. CALCP21A.271
*ENDIF CALCP21A.272
*ENDIF GSM3F403.780
end do ! do j = ... GSM3F403.781
CALCP21A.274
*IF DEF,GLOBAL CALCP21A.275
C Check North pole: CALCP21A.276
*IF DEF,MPP GSM3F403.782
IF (at_top_of_LPG) THEN GSM3F403.783
*ENDIF GSM3F403.784
if (eta_level(TOP_ROW_START+LAST_ROW_PT-1) .eq. rmdi) then GSM3F403.785
j=1 GSM3F403.786
ELSE GSM3F403.787
do i=TOP_ROW_START+ROW_LENGTH, GSM3F403.788
& TOP_ROW_START+ROW_LENGTH+LAST_ROW_PT-1 GSM3F403.789
if(eta_level(i).eq.rmdi) j=1 GSM3F403.790
enddo GSM3F403.791
endif GSM3F403.792
*IF DEF,MPP GSM3F403.793
CALL GCG_IMAX(
1,GC_ROW_GROUP,info,j) GSM3F403.794
*ENDIF GSM3F403.795
if (j.eq.1) then GSM3F403.796
do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1 GSM3F403.797
mask(i)=.false. GSM3F403.798
enddo GSM3F403.799
endif GSM3F403.800
*IF DEF,MPP GSM3F403.801
endif ! if at_top_of_LPG GSM3F403.802
*ENDIF GSM3F403.803
CALCP21A.285
C Check South pole: CALCP21A.286
j=0 CALCP21A.287
*IF DEF,MPP GSM3F403.804
IF (at_base_of_LPG) THEN GSM3F403.805
*ENDIF GSM3F403.806
if (eta_level(P_BOT_ROW_START) .eq. rmdi) then GSM3F403.807
j=1 GSM3F403.808
else GSM3F403.809
do i=P_BOT_ROW_START-ROW_LENGTH, GSM3F403.810
& P_BOT_ROW_START-ROW_LENGTH+LAST_ROW_PT-1 GSM3F403.811
if (eta_level(i) .eq. rmdi) j=1 GSM3F403.812
enddo GSM3F403.813
endif GSM3F403.814
*IF DEF,MPP GSM3F403.815
CALL GCG_IMAX(
1,GC_ROW_GROUP,info,j) GSM3F403.816
*ENDIF GSM3F403.817
GSM3F403.818
if (j.eq.1) then GSM3F403.819
do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1 GSM3F403.820
mask(i)=.false. GSM3F403.821
enddo GSM3F403.822
endif GSM3F403.823
*IF DEF,MPP GSM3F403.824
endif ! if at_base_of_LPG GSM3F403.825
*ENDIF GSM3F403.826
*ELSE CALCP21A.295
*IF DEF,MPP GSM3F403.827
IF (at_top_of_LPG) THEN GSM3F403.828
*ENDIF GSM3F403.829
do i=TOP_ROW_START,TOP_ROW_START+LAST_ROW_PT-1 GSM3F403.830
mask(i)=.false. GSM3F403.831
enddo GSM3F403.832
*IF DEF,MPP GSM3F403.833
ENDIF ! if at_top_of_LPG GSM3F403.834
GSM3F403.835
IF (at_base_of_LPG) THEN GSM3F403.836
*ENDIF GSM3F403.837
do i=P_BOT_ROW_START,P_BOT_ROW_START+LAST_ROW_PT-1 GSM3F403.838
mask(i)=.false. GSM3F403.839
enddo GSM3F403.840
*IF DEF,MPP GSM3F403.841
ENDIF ! if at_base_of_LPG GSM3F403.842
*ENDIF GSM3F403.843
*ENDIF CALCP21A.303
CALCP21A.304
CALCP21A.305
CL Section 4.2 Interpolate F3 (the Coriolis term) to the p-grid. CALCP21A.306
CL ~~~~~~~~~~~ CALCP21A.307
CALCP21A.308
call uv_to_p
CALCP21A.309
1 (f3,f3_p(row_length+1),u_field,p_field, CALCP21A.310
*IF DEF,MPP GSM3F403.844
2 row_length,U_LAST_ROW+1) GSM3F403.845
*ELSE GSM3F403.846
2 row_length,u_field/row_length) CALCP21A.311
*ENDIF GSM3F403.847
CALCP21A.312
CALCP21A.313
CL Section 4.3 Calculate the potential vorticity. CALCP21A.314
CL ~~~~~~~~~~~ CALCP21A.315
CALCP21A.316
do i = START_POINT_NO_HALO,END_P_POINT_NO_HALO GSM3F403.848
if (mask(i)) then CALCP21A.318
pvort_p(i) = (vorticity3(i) + vorticity4(i) - CALCP21A.319
& dtheta_dp(i)*(f3_p(i) + vorticity5(i))) * g CALCP21A.320
else CALCP21A.321
pvort_p(i) = rmdi CALCP21A.322
end if CALCP21A.323
enddo CALCP21A.324
CALCP21A.325
CL Section 4.3.1 Compute the potential vorticity at the Northern CALCP21A.326
CL ~~~~~~~~~~~~~ and Southern boundaries. CALCP21A.327
CL Note that under UPDATE identifier GLOBAL, as pv CALCP21A.328
CL is a scaler this is defined to be the mean of CALCP21A.329
CL the near pole row. CALCP21A.330
CL If any rmdi is present in the row next to the pole, CALCP21A.331
CL then the pole has a value rmdi. CALCP21A.332
CL ELSE, the value returned is missing data. CALCP21A.333
CALCP21A.334
C Work out the Northern boundary. CALCP21A.335
*IF DEF,MPP GSM3F403.849
if (at_top_of_LPG) then GSM3F403.850
*ENDIF GSM3F403.851
if (mask(TOP_ROW_START+LAST_ROW_PT-1)) then GSM3F403.852
CALCP21A.352
mn=0 GSM3F403.853
*IF -DEF,MPP GSM3F403.854
do i=TOP_ROW_START+row_length, GSM3F403.855
& TOP_ROW_START+2*row_length-1 GSM3F403.856
mn=mn+pvort_p(i) GSM3F403.857
enddo GSM3F403.858
*ELSE GSM3F403.859
*IF DEF,REPROD GSM3F403.860
call GCG_RVECSUMR(
GSM3F403.861
*ELSE GSM3F403.862
call GCG_RVECSUMF(
GSM3F403.863
*ENDIF GSM3F403.864
& row_length-2*EW_Halo,row_length-2*EW_Halo,1,1, GSM3F403.865
& pvort_p(TOP_ROW_START+FIRST_ROW_PT-1+row_length), GSM3F403.866
& GC_ROW_GROUP,info,mn) GSM3F403.867
*ENDIF GSM3F403.868
mn=mn/GLOBAL_ROW_LENGTH GSM3F403.869
CALCP21A.370
do i=TOP_ROW_START,TOP_ROW_START+row_length-1 GSM3F403.870
pvort_p(i)=mn GSM3F403.871
enddo GSM3F403.872
else GSM3F403.873
! Missing data GSM3F403.874
do i=TOP_ROW_START,TOP_ROW_START+row_length-1 GSM3F403.875
pvort_p(i)=rmdi GSM3F403.876
enddo GSM3F403.877
endif GSM3F403.878
*IF DEF,MPP GSM3F403.879
endif ! if at_top_of_LPG GSM3F403.880
*ENDIF GSM3F403.881
GSM3F403.882
! Southern boundary GSM3F403.883
*IF DEF,MPP GSM3F403.884
GSM3F403.885
if (at_base_of_LPG) then GSM3F403.886
*ENDIF GSM3F403.887
if (mask(P_BOT_ROW_START+FIRST_ROW_PT-1)) then GSM3F403.888
GSM3F403.889
mn=0 GSM3F403.890
*IF -DEF,MPP GSM3F403.891
do i=P_BOT_ROW_START-ROW_LENGTH, GSM3F403.892
& P_BOT_ROW_START-1 GSM3F403.893
mn=mn+pvort_p(i) GSM3F403.894
enddo GSM3F403.895
*ELSE GSM3F403.896
*IF DEF,REPROD GSM3F403.897
call GCG_RVECSUMR(
GSM3F403.898
*ELSE GSM3F403.899
call GCG_RVECSUMF(
GSM3F403.900
*ENDIF GSM3F403.901
& row_length-2*EW_Halo,row_length-2*EW_Halo,1,1, GSM3F403.902
& pvort_p(P_BOT_ROW_START+FIRST_ROW_PT-1-row_length), GSM3F403.903
& GC_ROW_GROUP,info,mn) GSM3F403.904
*ENDIF GSM3F403.905
mn=mn/GLOBAL_ROW_LENGTH GSM3F403.906
GSM3F403.907
do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1 GSM3F403.908
pvort_p(i)=mn GSM3F403.909
enddo GSM3F403.910
else GSM3F403.911
! Missing data GSM3F403.912
do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1 GSM3F403.913
pvort_p(i)=rmdi GSM3F403.914
enddo GSM3F403.915
endif GSM3F403.916
*IF DEF,MPP GSM3F403.917
endif ! if at_base_of_LPG GSM3F403.918
*ENDIF GSM3F403.919
GSM3F403.920
*IF DEF,MPP GSM3F403.921
! Fill in the N+S halo areas with rmdi GSM3F403.922
do i=1,FIRST_FLD_PT-1 GSM3F403.923
pvort_p(i)=rmdi GSM3F403.924
enddo GSM3F403.925
do i=LAST_P_FLD_PT+1,p_field GSM3F403.926
pvort_p(i)=rmdi GSM3F403.927
enddo GSM3F403.928
*ENDIF GSM3F403.929
CALCP21A.371
return CALCP21A.372
end CALCP21A.373
CALCP21A.374
*ENDIF CALCP21A.375