*IF DEF,A15_1A VORTI21A.2
C ******************************COPYRIGHT****************************** GTS2F400.11791
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.11792
C GTS2F400.11793
C Use, duplication or disclosure of this code is subject to the GTS2F400.11794
C restrictions as set forth in the contract. GTS2F400.11795
C GTS2F400.11796
C Meteorological Office GTS2F400.11797
C London Road GTS2F400.11798
C BRACKNELL GTS2F400.11799
C Berkshire UK GTS2F400.11800
C RG12 2SZ GTS2F400.11801
C GTS2F400.11802
C If no contract has been raised with this copy of the code, the use, GTS2F400.11803
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.11804
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.11805
C Modelling at the above address. GTS2F400.11806
C ******************************COPYRIGHT****************************** GTS2F400.11807
C GTS2F400.11808
CLL SUBROUTINE VORTIC2 ----------------------------------------------- VORTI21A.3
CLL VORTI21A.4
CLL Purpose: To compute 'vorticity2'. VORTI21A.5
CLL This is a term that is calculated here to be used later VORTI21A.6
CLL in the calculation of potential vorticity. VORTI21A.7
CLL Vorticity2 = 1/(rs*rs*cos(phi)) * VORTI21A.8
CLL (D(rs.v)/D(lambda)-D(rs.cos(phi).u)/D(phi)). VORTI21A.9
CLL VORTI21A.10
CLL Not suitable for single column use. VORTI21A.11
CLL Version for cray y-mp VORTI21A.12
CLL VORTI21A.13
CLL 7/10/92 Written by Simon Anderson VORTI21A.14
CLL VORTI21A.15
CLL Model Modification history from model version 3.0: VORTI21A.16
CLL version Date VORTI21A.17
CLL 3.1 14/01/93 Code inserted to get round FPP error in LAM mode. MM140193.1
CLL This code maybe removed when FPP is fixed. M. Mawson MM140193.2
CLL 3.2 28/07/93 Change subroutine name to uppercase for TS280793.27
CLL portability. Tracey Smith TS280793.28
!LL 4.3 18/02/97 Added ARGFLDPT arguments and MPP code P.Burton GSM3F403.1309
!LL 4.4 05/09/97 One point not initialised properly. S.D.Mullerworth GSM1F404.41
CLL VORTI21A.18
CLL Programming standard: Unified model documentation paper no. 4, VORTI21A.19
CLL standard b. version 2, dated 18/01/90 VORTI21A.20
CLL VORTI21A.21
CLL Logical components covered: D415 VORTI21A.22
CLL VORTI21A.23
CLL System task: D4 VORTI21A.24
CLL VORTI21A.25
CLL Documentation: U.M.D.P No 13. Derivation and Calculation of VORTI21A.26
CLL Unified Model Potential Vorticity. VORTI21A.27
CLL by Simon Anderson and Ian Roulstone. VORTI21A.28
CLL VORTI21A.29
CLLEND------------------------------------------------------------------ VORTI21A.30
VORTI21A.31
C*L ARGUMENTS:---------------------------------------------------------- VORTI21A.32
SUBROUTINE VORTIC2 1,1TS280793.29
1 (u_on_theta,v_on_theta,rs_on_theta, VORTI21A.34
2 rs_uv_on_theta,cos_u_latitude, VORTI21A.35
3 sec_p_latitude,vorticity2, VORTI21A.36
4 p_field,u_field,row_length, VORTI21A.37
*CALL ARGFLDPT 
GSM3F403.1310
5 latitude_step_inverse,longitude_step_inverse) VORTI21A.38
VORTI21A.39
implicit none VORTI21A.40
VORTI21A.41
C Input variables ------------------------------------------------------ VORTI21A.42
VORTI21A.43
integer VORTI21A.44
& p_field !IN Number of points in pressure field. VORTI21A.45
&,u_field !IN Number of points in velocity field. VORTI21A.46
&,row_length !IN Number of points per row. VORTI21A.47
*CALL TYPFLDPT GSM3F403.1311
VORTI21A.48
real VORTI21A.49
& u_on_theta(u_field) !IN Mass weighted u velocity. VORTI21A.50
&,v_on_theta(u_field) !IN Mass weighted v velocity* VORTI21A.51
& ! cos(latitude). VORTI21A.52
&,rs_on_theta(p_field) !IN Pseudo radius at p points. VORTI21A.53
&,rs_uv_on_theta(u_field)!IN Pseudo radius at u points. VORTI21A.54
&,cos_u_latitude(u_field)!IN Cos(lat) at u points. VORTI21A.55
&,sec_p_latitude(p_field)!IN 1/cos(lat) at p points. VORTI21A.56
&,latitude_step_inverse !IN 1/latitude increment. VORTI21A.57
&,longitude_step_inverse !IN 1/longitude increment. VORTI21A.58
VORTI21A.59
C Output variables ----------------------------------------------------- VORTI21A.60
VORTI21A.61
real VORTI21A.62
& vorticity2(p_field) !OUT Term used in potential vorticity eqn. VORTI21A.63
VORTI21A.64
C*---------------------------------------------------------------------- VORTI21A.65
C*L Workspace usage:- 3 local arrays required. VORTI21A.66
VORTI21A.67
real VORTI21A.68
& drsv_dlongitude(p_field) VORTI21A.69
&,drscosphiu_dlatitude(p_field) VORTI21A.70
&,drscosphiu_dlatitude2(u_field) VORTI21A.71
C*--------------------------------------------------------------------- VORTI21A.72
VORTI21A.73
C*L Define local variables: VORTI21A.74
integer VORTI21A.75
& i,j ! Loop counts. VORTI21A.76
*IF DEF,MPP,AND,DEF,GLOBAL GSM3F403.1312
&, info ! GCOM return code GSM3F403.1313
GSM3F403.1314
REAL GSM3F403.1315
& pole_sum(row_length) ! array for summing around pole GSM3F403.1316
*ENDIF GSM3F403.1317
VORTI21A.77
real VORTI21A.78
& scalar ! Local scalar. VORTI21A.79
VORTI21A.80
*IF DEF,GLOBAL VORTI21A.81
real sum_n,sum_s VORTI21A.82
*ENDIF VORTI21A.83
VORTI21A.84
VORTI21A.85
VORTI21A.86
CL--------------------------------------------------------------------- VORTI21A.87
CL Calculate 'vorticity2'. VORTI21A.88
CL--------------------------------------------------------------------- VORTI21A.89
VORTI21A.90
*IF -DEF,GLOBAL MM140193.3
CMIC$ PARALLEL SHARED(ROW_LENGTH,P_FIELD,VORTICITY2,SEC_P_LATITUDE, MM140193.4
CMIC$1 RS_ON_THETA,DRSCOSPHIU_DLATITUDE,U_FIELD, MM140193.5
CMIC$2 LATITUDE_STEP_INVERSE,RS_UV_ON_THETA,COS_U_LATITUDE, MM140193.6
CMIC$3 U_ON_THETA,LONGITUDE_STEP_INVERSE,V_ON_THETA,DRSV_DLONGITUDE) MM140193.7
CMIC$4 PRIVATE(J,I) MM140193.8
CMIC$ DO PARALLEL VECTOR MM140193.9
CDIR$ IVDEP MM140193.10
! Fujitsu vectorization directive GRB0F405.559
!OCL NOVREC GRB0F405.560
*ENDIF MM140193.11
C Calculate d(rsv)/d(lambda). VORTI21A.91
do 110 i=TOP_ROW_START+1,LAST_U_FLD_PT GSM3F403.1318
drsv_dlongitude(i) = longitude_step_inverse* VORTI21A.93
& (rs_uv_on_theta(i)*v_on_theta(i) - VORTI21A.94
& rs_uv_on_theta(i-1)*v_on_theta(i-1)) VORTI21A.95
110 continue VORTI21A.96
VORTI21A.97
C Calculate d(rscosphiu)/d(phi). VORTI21A.98
*IF -DEF,GLOBAL MM140193.12
CMIC$ DO PARALLEL VECTOR MM140193.13
CDIR$ IVDEP MM140193.14
! Fujitsu vectorization directive GRB0F405.561
!OCL NOVREC GRB0F405.562
*ENDIF MM140193.15
do 120 i=START_POINT_NO_HALO,LAST_U_FLD_PT GSM3F403.1319
drscosphiu_dlatitude(i) = latitude_step_inverse* VORTI21A.100
& (rs_uv_on_theta(i-row_length)* VORTI21A.101
& cos_u_latitude(i-row_length)* VORTI21A.102
& u_on_theta(i-row_length) - VORTI21A.103
& rs_uv_on_theta(i)* VORTI21A.104
& cos_u_latitude(i)* VORTI21A.105
& u_on_theta(i)) VORTI21A.106
120 continue VORTI21A.107
VORTI21A.108
*IF DEF,GLOBAL VORTI21A.109
C Calculate average of drscosphiu_dlatitude at p-points. VORTI21A.110
do 130 i=START_POINT_NO_HALO+1,END_P_POINT_NO_HALO GSM3F403.1320
drscosphiu_dlatitude2(i) = drscosphiu_dlatitude(i) + VORTI21A.112
& drscosphiu_dlatitude(i-1) VORTI21A.113
130 continue VORTI21A.114
VORTI21A.115
C Now do first point on each slice for VORTI21A.116
C drsv_dlongitude and drscosphiu_dlatitude2. VORTI21A.117
*IF -DEF,MPP GSM3F403.1321
i=FIRST_FLD_PT GSM3F403.1322
drsv_dlongitude(i) = longitude_step_inverse * VORTI21A.119
& (rs_uv_on_theta(i)*v_on_theta(i) - VORTI21A.120
& rs_uv_on_theta(i+row_length-1)* VORTI21A.121
& v_on_theta(i+row_length-1)) VORTI21A.122
do 140 i=FIRST_FLD_PT+row_length,LAST_U_FLD_PT,row_length GSM3F403.1323
drsv_dlongitude(i) = longitude_step_inverse * VORTI21A.124
& (rs_uv_on_theta(i)*v_on_theta(i) - VORTI21A.125
& rs_uv_on_theta(i+row_length-1)* VORTI21A.126
& v_on_theta(i+row_length-1)) VORTI21A.127
drscosphiu_dlatitude2(i)=drscosphiu_dlatitude(i)+ VORTI21A.128
& drscosphiu_dlatitude(i-1+row_length) VORTI21A.129
140 continue VORTI21A.130
*ELSE GSM3F403.1324
i=TOP_ROW_START GSM1F404.42
! Put a sensible number in the first element (halo) GSM3F403.1326
drsv_dlongitude(i) = drsv_dlongitude(i+1) GSM3F403.1327
i=START_POINT_NO_HALO GSM1F404.43
drscosphiu_dlatitude2(i)=drscosphiu_dlatitude2(i+1) GSM3F403.1328
*ENDIF GSM3F403.1329
VORTI21A.131
C Calculate vorticity2. VORTI21A.132
VORTI21A.133
do 150 j=START_POINT_NO_HALO,END_P_POINT_NO_HALO GSM3F403.1330
vorticity2(j)=sec_p_latitude(j)/ VORTI21A.135
& (rs_on_theta(j)*rs_on_theta(j))* VORTI21A.136
& .5*(drsv_dlongitude(j)+ VORTI21A.137
& drsv_dlongitude(j-row_length)- VORTI21A.138
& drscosphiu_dlatitude2(j)) VORTI21A.139
150 continue VORTI21A.140
VORTI21A.141
*ELSE VORTI21A.142
drsv_dlongitude(FIRST_FLD_PT) = 0. GSM3F403.1331
VORTI21A.144
C Calculate vorticity2. VORTI21A.145
VORTI21A.146
CMIC$ DO PARALLEL VECTOR MM140193.16
CDIR$ IVDEP MM140193.17
! Fujitsu vectorization directive GRB0F405.563
!OCL NOVREC GRB0F405.564
do 130 j=START_POINT_NO_HALO+1,END_P_POINT_NO_HALO-1 GSM3F403.1332
vorticity2(j)=sec_p_latitude(j)/ VORTI21A.148
& (rs_on_theta(j)*rs_on_theta(j))* VORTI21A.149
& .5*(drsv_dlongitude(j)+ VORTI21A.150
& drsv_dlongitude(j-row_length)- VORTI21A.151
& drscosphiu_dlatitude(j)- VORTI21A.152
& drscosphiu_dlatitude(j-1)) VORTI21A.153
130 continue VORTI21A.154
VORTI21A.155
C Zero vorticity2 on boundaries. VORTI21A.156
*IF DEF,MPP GSM3F403.1333
if (at_top_of_LPG) then GSM3F403.1334
*ENDIF GSM3F403.1335
! Northern boundary GSM3F403.1336
do j=TOP_ROW_START,TOP_ROW_START+row_length-1 GSM3F403.1337
vorticity2(j) = 0.0 GSM3F403.1338
enddo GSM3F403.1339
*IF DEF,MPP GSM3F403.1340
endif GSM3F403.1341
GSM3F403.1342
if (at_base_of_LPG) then GSM3F403.1343
*ENDIF GSM3F403.1344
! Southern boundary GSM3F403.1345
do j=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1 GSM3F403.1346
vorticity2(j) = 0.0 GSM3F403.1347
enddo GSM3F403.1348
*IF DEF,MPP GSM3F403.1349
endif GSM3F403.1350
GSM3F403.1351
if (at_left_of_LPG) then GSM3F403.1352
*ENDIF GSM3F403.1353
! Western boundary GSM3F403.1354
do j=TOP_ROW_START+FIRST_ROW_PT-1+row_length, GSM3F403.1355
& END_P_POINT_NO_HALO, GSM3F403.1356
& row_length GSM3F403.1357
vorticity2(j) = 0.0 GSM3F403.1358
enddo GSM3F403.1359
*IF DEF,MPP GSM3F403.1360
endif GSM3F403.1361
GSM3F403.1362
if (at_right_of_LPG) then GSM3F403.1363
*ENDIF GSM3F403.1364
! Eastern boundary GSM3F403.1365
do j=TOP_ROW_START+LAST_ROW_PT-1+row_length, GSM3F403.1366
& END_P_POINT_NO_HALO, GSM3F403.1367
& row_length GSM3F403.1368
vorticity2(j)=0.0 GSM3F403.1369
enddo GSM3F403.1370
*IF DEF,MPP GSM3F403.1371
endif GSM3F403.1372
*ENDIF GSM3F403.1373
VORTI21A.165
CMIC$ END PARALLEL MM140193.22
*ENDIF VORTI21A.166
VORTI21A.167
*IF DEF,GLOBAL VORTI21A.168
C Calculate vorticity2 at poles by summing drscosphiu/d(lat) around VORTI21A.169
C Polar circle and averaging. VORTI21A.170
scalar = latitude_step_inverse/GLOBAL_ROW_LENGTH GSM3F403.1374
*IF DEF,MPP GSM3F403.1375
if (at_top_of_LPG) then GSM3F403.1376
*ENDIF GSM3F403.1377
sum_n=0.0 GSM3F403.1378
*IF -DEF,MPP GSM3F403.1379
do i=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1 GSM3F403.1380
sum_n = sum_n - rs_uv_on_theta(i)*cos_u_latitude(i)* GSM3F403.1381
& u_on_theta(i)*scalar GSM3F403.1382
enddo GSM3F403.1383
*ELSE GSM3F403.1384
do i=1,ROW_LENGTH-2*EW_Halo GSM3F403.1385
j=TOP_ROW_START+FIRST_ROW_PT+i-2 GSM3F403.1386
pole_sum(i)=-rs_uv_on_theta(j)*cos_u_latitude(j)* GSM3F403.1387
& u_on_theta(j)*scalar GSM3F403.1388
enddo GSM3F403.1389
GSM3F403.1390
*IF DEF,REPROD GSM3F403.1391
CALL GCG_RVECSUMR( GSM3F403.1392
*ELSE GSM3F403.1393
CALL GCG_RVECSUMF( GSM3F403.1394
*ENDIF GSM3F403.1395
& ROW_LENGTH-2*EW_Halo,ROW_LENGTH-2*EW_Halo,1,1, GSM3F403.1396
& pole_sum, GSM3F403.1397
& GC_ROW_GROUP,info,sum_n) GSM3F403.1398
*ENDIF GSM3F403.1399
do i=TOP_ROW_START,TOP_ROW_START+ROW_LENGTH-1 GSM3F403.1400
vorticity2(i)=-sum_n*sec_p_latitude(i)/ GSM3F403.1401
& (rs_on_theta(i)*rs_on_theta(i)) GSM3F403.1402
enddo GSM3F403.1403
*IF DEF,MPP GSM3F403.1404
endif GSM3F403.1405
GSM3F403.1406
if (at_base_of_LPG) then GSM3F403.1407
*ENDIF GSM3F403.1408
sum_s=0.0 GSM3F403.1409
*IF -DEF,MPP GSM3F403.1410
do i=P_BOT_ROW_START-row_length,P_BOT_ROW_START-1 GSM3F403.1411
sum_s = sum_s + rs_uv_on_theta(i)*cos_u_latitude(i)* GSM3F403.1412
& u_on_theta(i)*scalar GSM3F403.1413
enddo GSM3F403.1414
*ELSE GSM3F403.1415
do i=1,ROW_LENGTH-2*EW_Halo GSM3F403.1416
j=P_BOT_ROW_START+FIRST_ROW_PT-row_length-2+i GSM3F403.1417
pole_sum(i)=rs_uv_on_theta(j)*cos_u_latitude(j)* GSM3F403.1418
& u_on_theta(j)*scalar GSM3F403.1419
enddo GSM3F403.1420
GSM3F403.1421
*IF DEF,REPROD GSM3F403.1422
CALL GCG_RVECSUMR( GSM3F403.1423
*ELSE GSM3F403.1424
CALL GCG_RVECSUMF( GSM3F403.1425
*ENDIF GSM3F403.1426
& ROW_LENGTH-2*EW_Halo,ROW_LENGTH-2*EW_Halo,1,1, GSM3F403.1427
& pole_sum, GSM3F403.1428
& GC_ROW_GROUP,info,sum_s) GSM3F403.1429
*ENDIF GSM3F403.1430
do i=P_BOT_ROW_START,P_BOT_ROW_START+row_length-1 GSM3F403.1431
vorticity2(i)=-sum_s*sec_p_latitude(i)/ GSM3F403.1432
& (rs_on_theta(i)*rs_on_theta(i)) GSM3F403.1433
enddo GSM3F403.1434
*IF DEF,MPP GSM3F403.1435
endif GSM3F403.1436
*ENDIF GSM3F403.1437
VORTI21A.191
*ENDIF VORTI21A.192
*IF DEF,MPP GSM3F403.1438
! Set rest of array to sensible values GSM3F403.1439
CALL SWAPBOUNDS(vorticity2,ROW_LENGTH,tot_P_ROWS, GSM3F403.1440
& EW_Halo,NS_Halo,1) GSM3F403.1441
*ENDIF GSM3F403.1442
VORTI21A.193
CL End of routine vortic2 VORTI21A.194
VORTI21A.195
return VORTI21A.196
end VORTI21A.197
*ENDIF VORTI21A.198