SUBROUTINE SETCONA 1,8SETCONA1.30
C Input data SETCONA1.31
& (AK,BK,DELTA_AK,DELTA_BK,PSTAR, SETCONA1.32
& THETA,Q,RHCRIT, SETCONA1.33
& SATCON,SMVCSAT,SMVCWT, SETCONA1.34
! & C_EAG,LAND, ADM2F404.235
& C_EAG,LAND,L_LSPICE, ADM2F404.236
& RHCPT, ICE_FRAC, BL_LEVELS, L_RHCPT, ASK1F405.44
C size and control variables SETCONA1.36
& P_FIELD,P_LEVELS,P_ROWS,U_ROWS,U_FIELD,ROW_LENGTH, SETCONA1.37
& LAND_FIELD,Q_LEVELS, SETCONA1.38
C other constants SETCONA1.39
& EW_SPACE,NS_SPACE,FIRST_LAT,FIRST_LONG, SETCONA1.40
& PHI0,LAMBDA0, SETCONA1.41
C logical control switches ATJ0F403.11
& L_OLD_PWTS, ATJ0F403.12
C output data SETCONA1.42
& AKH,BKH,AK_TO_THE_KAPPA,BK_TO_THE_KAPPA,AKH_TO_THE_KAPPA, SETCONA1.43
& BKH_TO_THE_KAPPA,COS_U_LATITUDE,SEC_U_LATITUDE, SETCONA1.44
& SIN_U_LATITUDE, GJT1F304.19
& TAN_U_LATITUDE,COS_P_LATITUDE,SEC_P_LATITUDE, SETCONA1.45
& SIN_LONGITUDE,COS_LONGITUDE,TRUE_LONGITUDE, SETCONA1.46
& F1,F2,F3,F3_P,TRIGS,IFAX,P_EXNER,TWO_D_GRID_CORRECTION, SETCONA1.47
& QCL,QCF,SOILB,LAND_LIST,CLOUD_LEVELS, SETCONA1.48
& ETA_SPLIT,NUM_CLOUD_TYPES,LOW_BOT_LEVEL,LOW_TOP_LEVEL, SETCONA1.49
& MED_BOT_LEVEL,MED_TOP_LEVEL,HIGH_BOT_LEVEL,HIGH_TOP_LEVEL, SETCONA1.50
& ETA_MATRIX_INV,MATRIX_POLY_ORDER,ICODE,CMESSAGE) SETCONA1.51
SETCONA1.52
IMPLICIT NONE SETCONA1.53
SETCONA1.54
INTEGER SETCONA1.55
& P_FIELD,P_LEVELS,P_ROWS,U_ROWS,U_FIELD,ROW_LENGTH,ICODE, SETCONA1.56
& LAND_FIELD,Q_LEVELS,BL_LEVELS, ASK1F405.45
& LAND_LIST(LAND_FIELD),! List of land points SETCONA1.58
& IFAX(10), ! factors for FFT SETCONA1.59
& NUM_CLOUD_TYPES, ! No of cld types for diag cloud lyr SETCONA1.60
& CLOUD_LEVELS, ! No of cloudy levels in the model SETCONA1.61
& LOW_BOT_LEVEL, ! Bottom level of lowest cloud type SETCONA1.62
& LOW_TOP_LEVEL, ! Top " " " " " SETCONA1.63
& MED_BOT_LEVEL, ! Bottom " " med " " SETCONA1.64
& MED_TOP_LEVEL, ! Top " " " " " SETCONA1.65
& HIGH_BOT_LEVEL, ! Bottom " " top " " SETCONA1.66
& HIGH_TOP_LEVEL, ! Top " " " " " SETCONA1.67
& MATRIX_POLY_ORDER ! Order of polynomial used to SETCONA1.68
& ! calculate ETA_HALF inverse matrix SETCONA1.69
SETCONA1.70
SETCONA1.71
CHARACTER*(80) TS150793.173
& CMESSAGE ! Error message if ICODE >0 SETCONA1.73
SETCONA1.74
REAL SETCONA1.75
& AK(P_LEVELS), ! mid-layer A values for hybrid coordinates SETCONA1.76
& BK(P_LEVELS), ! mid-layer B values for hybrid coordinates SETCONA1.77
& DELTA_AK(P_LEVELS), ! layer thickness A & B values - SETCONA1.78
& DELTA_BK(P_LEVELS), ! for hybrid coordinates SETCONA1.79
& PSTAR(P_FIELD), ! SETCONA1.80
& THETA(P_FIELD,P_LEVELS), SETCONA1.81
& Q(P_FIELD,Q_LEVELS), SETCONA1.82
& RHCRIT(Q_LEVELS), SETCONA1.83
& RHCPT(P_FIELD,Q_LEVELS), ! OUT ASK1F405.46
& ICE_FRAC(P_FIELD), ! IN ASK1F405.47
& ETA_SPLIT(NUM_CLOUD_TYPES+1), SETCONA1.84
& ETA_MATRIX_INV(MATRIX_POLY_ORDER,MATRIX_POLY_ORDER,P_LEVELS) SETCONA1.85
SETCONA1.86
REAL SETCONA1.87
& EW_SPACE, ! EW (x) grid spacing in degrees SETCONA1.88
& NS_SPACE, ! NS (y) grid spacing in degrees SETCONA1.89
& FIRST_LAT, ! Latitude of first PTR row in degrees SETCONA1.90
C ! (latitude 90>-90) SETCONA1.91
& FIRST_LONG, ! Longitude of 1st PTR pt on row in degree SETCONA1.92
C ! Longitude 0>360) SETCONA1.93
& PHI0, ! Real latitude of 'pseudo' N pole in degrees SETCONA1.94
& LAMBDA0 ! Real longitude of 'pseudo' N pole in degrees SETCONA1.95
SETCONA1.96
REAL SETCONA1.97
& AKH(P_LEVELS+1), ! layer boudary A & B values SETCONA1.98
& BKH(P_LEVELS+1), ! coordinates for hybrid SETCONA1.99
& AK_TO_THE_KAPPA(P_LEVELS), ! (AK/PREF)**KAPPA SETCONA1.100
& BK_TO_THE_KAPPA(P_LEVELS), ! (BK/PREF)**KAPPA SETCONA1.101
& AKH_TO_THE_KAPPA(P_LEVELS+1), ! (AKH/PREF)**KAPPA SETCONA1.102
& BKH_TO_THE_KAPPA(P_LEVELS+1), ! (BKH/PREF)**KAPPA SETCONA1.103
& COS_U_LATITUDE(U_FIELD), ! SETCONA1.104
& SEC_U_LATITUDE(U_FIELD), ! SETCONA1.105
& SIN_U_LATITUDE(U_FIELD), ! GJT1F304.20
& TAN_U_LATITUDE(U_FIELD), ! SETCONA1.106
& COS_P_LATITUDE(P_FIELD), ! SETCONA1.107
& SEC_P_LATITUDE(P_FIELD), ! SETCONA1.108
& SIN_LONGITUDE(ROW_LENGTH), ! sin & cos of longitude SETCONA1.109
& COS_LONGITUDE(ROW_LENGTH), ! relative to coordinate pole. SETCONA1.110
C ! 1st pt long = 0 SETCONA1.111
& TRUE_LONGITUDE(P_FIELD), ! longitude relative to Greenwich SETCONA1.112
& F1(U_FIELD), ! Coriolis terms at u points as SETCONA1.113
& F2(U_FIELD), ! defined in documentation SETCONA1.114
& F3(U_FIELD), ! paper 10 equation 4 SETCONA1.115
& F3_P(P_FIELD), ! Coriolis parameter at P points SETCONA1.116
& SATCON(LAND_FIELD), ! SETCONA1.117
& SMVCSAT(LAND_FIELD), ! SETCONA1.118
& SMVCWT(LAND_FIELD), ! SETCONA1.119
& C_EAG(LAND_FIELD) ! SETCONA1.120
SETCONA1.121
REAL SETCONA1.122
& P_EXNER(P_FIELD,P_LEVELS+1), SETCONA1.123
& TRIGS(ROW_LENGTH),! coefficients for FFT SETCONA1.124
& TWO_D_GRID_CORRECTION(P_FIELD/ROW_LENGTH),! Correction SETCONA1.125
C ! used in filtering to represent SETCONA1.126
C ! 2D nature of grid. SETCONA1.127
& SOILB(LAND_FIELD), SETCONA1.128
& QCL(P_FIELD,Q_LEVELS),QCF(P_FIELD,Q_LEVELS), SETCONA1.129
& CLOUD_BOUND(NUM_CLOUD_TYPES+1) SETCONA1.130
SETCONA1.131
LOGICAL SETCONA1.132
& LAND(P_FIELD) SETCONA1.133
& ,L_LSPICE ADM2F404.237
& ,L_RHCPT ASK1F405.48
&, L_OLD_PWTS ! Switch for HADCM2 polar weights ATJ0F403.13
SETCONA1.134
C* --------------------------------------------------------------------- SETCONA1.135
SETCONA1.136
*CALL CPHYSCON 
SETCONA1.137
*CALL C_ETA_PMSL GDR3F404.13
*IF DEF,MPP GPB0F305.336
*CALL PARVARS GPB0F305.337
*IF DEF,MPP GDS1F402.398
*CALL ACPARM GDS1F402.399
*CALL MPPAC GDS1F402.400
*ENDIF GDS1F402.401
*ENDIF GPB0F305.338
*CALL TYPFLDPT ASK1F405.49
SETCONA1.138
C Local variables SETCONA1.139
REAL SETCONA1.140
& CLOUD_FRACTION(P_FIELD,Q_LEVELS), SETCONA1.141
& LS_GRID_QC(P_FIELD,Q_LEVELS), ! LOCAL Qc from LS_CLD AYY2F400.47
& LS_BS(P_FIELD,Q_LEVELS), ! LOCAL bs from LS_CLD AYY2F400.48
& LATITUDE(P_FIELD), SETCONA1.142
& LONGITUDE(P_FIELD), SETCONA1.143
& POLE_LAT, SETCONA1.144
& ETA_LEV, ! ETA at model level GDR3F404.14
& ETA(P_LEVELS), SETCONA1.145
& POLE_LONG, SETCONA1.146
& ETA_HALF(P_LEVELS+1), SETCONA1.147
& ETA_MATRIX(MATRIX_POLY_ORDER,MATRIX_POLY_ORDER,P_LEVELS) SETCONA1.148
SETCONA1.149
REAL QPOS,QNEG ADM2F404.238
INTEGER IQNEG(P_FIELD) ADM2F404.239
SETCONA1.150
*IF -DEF,GLOBAL SETCONA1.151
SETCONA1.152
& ,TRUE_LATITUDE(P_FIELD) SETCONA1.153
SETCONA1.154
*ENDIF SETCONA1.155
SETCONA1.156
INTEGER SETCONA1.157
& I,KK,II, SETCONA1.158
& LEVEL, SETCONA1.159
& ROW, SETCONA1.160
& POINT, SETCONA1.161
& FIRST_POINT, ASK1F405.50
& LAST_POINT, ASK1F405.51
& POINTS, ASK1F405.52
& JS, ASK1F405.53
& N21 SETCONA1.162
SETCONA1.163
SETCONA1.164
C External subroutines called : SETCONA1.165
EXTERNAL SETCONA1.166
& EQTOLL, SETCONA1.167
& GLUE_CLD, AYY2F400.49
& RHCRIT_CALC, ASK1F405.54
& SETFILT, SETCONA1.169
& MATINV SETCONA1.170
SETCONA1.171
REAL PU, PL SETCONA1.172
*CALL P_EXNERC SETCONA1.173
*CALL SETFLDPT ASK1F405.55
SETCONA1.174
SETCONA1.175
CL Internal structure : SETCONA1.176
CL 1. Calculate extra functions of AK,BK from dump information SETCONA1.177
C calculate layer boundary values of AK, BK SETCONA1.178
AKH(1) = 0.0 SETCONA1.179
SETCONA1.180
*IF -DEF,STRAT SETCONA1.181
BKH(1) = 1.0 SETCONA1.182
SETCONA1.183
*ELSE SETCONA1.184
BKH(1) = 0.0 SETCONA1.185
SETCONA1.186
*ENDIF SETCONA1.187
SETCONA1.188
DO LEVEL = 2, P_LEVELS + 1 SETCONA1.189
AKH(LEVEL) = AKH(LEVEL-1) + DELTA_AK(LEVEL-1) SETCONA1.190
BKH(LEVEL) = BKH(LEVEL-1) + DELTA_BK(LEVEL-1) SETCONA1.191
SETCONA1.192
IF(BK(LEVEL-1).EQ.0.0) THEN SETCONA1.193
BKH(LEVEL- 1) = 0.0 SETCONA1.194
BKH(LEVEL) = 0.0 SETCONA1.195
SETCONA1.196
END IF SETCONA1.197
END DO SETCONA1.198
SETCONA1.199
C calculate (A/PREF)**KAPPA, (B/PREF)**KAPPA for layer mid-points SETCONA1.200
DO LEVEL = 1, P_LEVELS SETCONA1.201
IF (AK(LEVEL) .GT. 0.0) THEN SETCONA1.202
AK_TO_THE_KAPPA(LEVEL) = (AK(LEVEL) / PREF)**KAPPA SETCONA1.203
SETCONA1.204
ELSE SETCONA1.205
AK_TO_THE_KAPPA(LEVEL) = 0.0 SETCONA1.206
SETCONA1.207
END IF SETCONA1.208
SETCONA1.209
IF (BK(LEVEL) .GT. 0.0) THEN SETCONA1.210
BK_TO_THE_KAPPA(LEVEL) = (BK(LEVEL) / PREF)**KAPPA SETCONA1.211
SETCONA1.212
ELSE SETCONA1.213
BK_TO_THE_KAPPA(LEVEL) = 0.0 SETCONA1.214
SETCONA1.215
END IF SETCONA1.216
END DO SETCONA1.217
SETCONA1.218
C calculate (A/PREF)**KAPPA, (B/PREF)**KAPPA, for layer boundaries SETCONA1.219
DO LEVEL = 1, P_LEVELS + 1 SETCONA1.220
IF (AKH(LEVEL) .GT. 0.0) THEN SETCONA1.221
AKH_TO_THE_KAPPA(LEVEL) = (AKH(LEVEL) / PREF)**KAPPA SETCONA1.222
SETCONA1.223
ELSE SETCONA1.224
AKH_TO_THE_KAPPA(LEVEL) = 0.0 SETCONA1.225
SETCONA1.226
END IF SETCONA1.227
SETCONA1.228
IF (BKH(LEVEL) .GT. 0.0) THEN SETCONA1.229
BKH_TO_THE_KAPPA(LEVEL) = (BKH(LEVEL)/PREF)**KAPPA SETCONA1.230
SETCONA1.231
ELSE SETCONA1.232
BKH_TO_THE_KAPPA(LEVEL) = 0.0 SETCONA1.233
SETCONA1.234
END IF SETCONA1.235
END DO SETCONA1.236
SETCONA1.237
C Calculate P_EXNER SETCONA1.238
DO LEVEL = 1, P_LEVELS + 1 SETCONA1.239
DO I = 1, P_FIELD SETCONA1.240
P_EXNER(I,LEVEL) = ((AKH(LEVEL) + BKH(LEVEL) * PSTAR(I)) SETCONA1.241
& / PREF)**KAPPA SETCONA1.242
SETCONA1.243
END DO SETCONA1.244
END DO SETCONA1.245
SETCONA1.246
CL 2. Calculate constants at u grid points SETCONA1.247
CL cos, sec, tan (u_latitude). Coriolis terms for GLOBAL model SETCONA1.248
POLE_LAT = PHI0 * PI_OVER_180 SETCONA1.249
POLE_LONG = LAMBDA0 * PI_OVER_180 SETCONA1.250
SETCONA1.251
*IF DEF,MPP GPB6F405.3
WRITE(6,*) 'Decomposition information for PE ',mype GPB6F405.4
*IF DEF,GLOBAL GPB6F405.5
WRITE(6,*) 'Running in GLOBAL configuration' GPB6F405.6
*ELSE GPB6F405.7
WRITE(6,*) 'Running in LAM configuration' GPB6F405.8
*ENDIF GPB6F405.9
WRITE(6,*) 'Overall decomposition shape is (EW x NS) ', GPB6F405.10
& nproc_x, ' x ',nproc_y GPB6F405.11
WRITE(6,*) 'Domain size (inc. halos) (Columns x Rows) ', GPB6F405.12
& ROW_LENGTH,' x ',P_ROWS GPB6F405.13
WRITE(6,*) 'Halo size (EW Halos , NS Halos) ',Offx, ' , ', Offy GPB6F405.14
WRITE(6,*) 'Position in the Logical Processor Grid is (x,y) ', GPB6F405.15
& gridpos(1), ',', gridpos(2) GPB6F405.16
WRITE(6,*) '(NB LPG co-ordinates start at zero)' GPB6F405.17
WRITE(6,*) 'Grid Co-ords of top-left corner are (x,y) ', GPB6F405.18
& datastart(1), ',', datastart(2) GPB6F405.19
WRITE(6,*) 'Location of top-left corner of T grid ', GPB6F405.20
& '(longitude,latitude) ', GPB6F405.21
& (FIRST_LONG + EW_SPACE * (datastart(1) -1.0)),',', GPB6F405.22
& (FIRST_LAT - NS_SPACE * (datastart(2) -1.0)) GPB6F405.23
*IF -DEF,GLOBAL GPB6F405.24
WRITE(6,*) '(NB Co-ordinates are for model-grid which is rotated ' GPB6F405.25
WRITE(6,*) ' for LAM models. See UMDP S1.)' GPB6F405.26
*ENDIF GPB6F405.27
*ENDIF GPB6F405.28
DO ROW = 1, U_ROWS SETCONA1.252
DO POINT = 1, ROW_LENGTH SETCONA1.253
I = POINT + (ROW - 1) * ROW_LENGTH SETCONA1.254
*IF -DEF,MPP GPB0F305.339
LATITUDE(I) = (FIRST_LAT - NS_SPACE * (ROW - 0.5)) SETCONA1.255
LONGITUDE(I) = (FIRST_LONG + EW_SPACE * (POINT - 0.5)) SETCONA1.256
*ELSE GPB0F305.340
! take account of our offset in the global data GPB0F305.341
LATITUDE(I) = (FIRST_LAT - NS_SPACE * GPB0F305.342
& ((ROW + datastart(2) - Offy - 1) - 0.5)) GPB0F305.343
LONGITUDE(I) = (FIRST_LONG + EW_SPACE * GPB0F305.344
& ((POINT + datastart(1)- Offx - 1) - 0.5)) GPB0F305.345
*ENDIF GPB0F305.346
SETCONA1.257
SETCONA1.258
LATITUDE(I) = LATITUDE(I) * PI_OVER_180 SETCONA1.259
LONGITUDE(I) = LONGITUDE(I) * PI_OVER_180 SETCONA1.260
COS_U_LATITUDE(I) = COS(LATITUDE(I)) SETCONA1.261
SEC_U_LATITUDE(I) = 1.0 / COS_U_LATITUDE(I) SETCONA1.262
SIN_U_LATITUDE(I)= SIN(LATITUDE(I)) GJT1F304.21
TAN_U_LATITUDE(I) = TAN(LATITUDE(I)) SETCONA1.263
SETCONA1.264
*IF DEF,GLOBAL SETCONA1.265
SETCONA1.266
F1(I)=0.0 SETCONA1.267
F2(I) = 2.0 * OMEGA * COS_U_LATITUDE(I) SETCONA1.268
F3(I) = 2.0 * OMEGA * SIN(LATITUDE(I)) SETCONA1.269
*ELSE SETCONA1.270
SETCONA1.271
F1(I) = -2.0 * OMEGA * SIN(LATITUDE(I)) * COS(POLE_LAT) SETCONA1.272
F2(I) = 2.0 * OMEGA * (COS(LATITUDE(I)) * SIN(POLE_LAT) - SETCONA1.273
& SIN(LATITUDE(I)) * COS(LONGITUDE(I)) * COS(POLE_LAT)) SETCONA1.274
F3(I) = 2.0 * OMEGA * (SIN(LATITUDE(I)) * SIN(POLE_LAT) + SETCONA1.275
& COS(LATITUDE(I)) * COS(LONGITUDE(I)) * COS(POLE_LAT)) SETCONA1.276
SETCONA1.277
*ENDIF SETCONA1.278
SETCONA1.279
END DO SETCONA1.280
END DO SETCONA1.281
SETCONA1.282
SETCONA1.283
CL SETCONA1.284
CL 3. Calculate constants on p grid points SETCONA1.285
CL SETCONA1.286
C calculate cos, sec (p_latitude) longitudes for global model SETCONA1.287
DO ROW = 1, P_ROWS SETCONA1.288
DO POINT = 1, ROW_LENGTH SETCONA1.289
I = POINT + (ROW - 1) * ROW_LENGTH SETCONA1.290
*IF -DEF,MPP GPB0F305.347
LATITUDE(I) = (FIRST_LAT - NS_SPACE * (ROW - 1.0)) SETCONA1.291
LONGITUDE(I) = (FIRST_LONG + EW_SPACE * (POINT - 1.0)) SETCONA1.292
*ELSE GPB0F305.348
! take account of our offset in the global data GPB0F305.349
LATITUDE(I) = (FIRST_LAT - NS_SPACE * GPB0F305.350
& ((ROW + datastart(2) - Offy - 1) - 1.0)) GPB0F305.351
LONGITUDE(I) = (FIRST_LONG + EW_SPACE * GPB0F305.352
& ((POINT + datastart(1)- Offx - 1) - 1.0)) GPB0F305.353
*ENDIF GPB0F305.354
SETCONA1.293
END DO SETCONA1.294
END DO SETCONA1.295
*IF DEF,MPP GDS1F402.402
lat_n=LATITUDE(1) GDS1F402.403
lat_s=LATITUDE(P_ROWS*ROW_LENGTH) GDS1F402.404
GDS1F402.405
long_w=LONGITUDE(1) GDS1F402.406
long_w_model=long_w GDS1F402.407
long_e=LONGITUDE(P_ROWS*ROW_LENGTH) GDS1F402.408
long_e_model=long_e GDS1F402.409
GDS1F402.410
if(long_w.gt.180.0)long_w=long_w-360.0 GDS1F402.411
if(long_e.gt.180.0)long_e=long_e-360.0 GDS1F402.412
*IF -DEF,GLOBAL GDS1F402.413
if(long_w_model.ge.360.0)long_w_model=long_w_model-360. GDS1F402.414
if(long_e_model.ge.360.0)long_e_model=long_e_model-360. GDS1F402.415
*ENDIF GDS1F402.416
*IF DEF,GLOBAL GSM5F405.2
! Without SWAPBOUNDS, west-most halo pt. is -ve equivalent of east-most GSM5F405.3
! non-halo pt. (eg -3 deg and +357 deg) causing differences in calcs GSM5F405.4
! of halo points. GSM5F405.5
CALL SWAPBOUNDS(LONGITUDE,ROW_LENGTH,P_ROWS,1,0,1) GSM5F405.6
*ENDIF GSM5F405.7
GSM5F405.8
*ENDIF GDS1F402.417
SETCONA1.296
*IF -DEF,GLOBAL SETCONA1.297
SETCONA1.298
CL calculate true longitude for limited area model using sub'tine EQTOLL SETCONA1.299
CALL EQTOLL(LATITUDE, LONGITUDE, TRUE_LATITUDE, TRUE_LONGITUDE, SETCONA1.300
& PHI0, LAMBDA0, P_FIELD) SETCONA1.301
SETCONA1.302
*ENDIF SETCONA1.303
SETCONA1.304
DO I = 1, P_FIELD SETCONA1.305
LATITUDE(I) = LATITUDE(I) * PI_OVER_180 SETCONA1.306
LONGITUDE(I) = LONGITUDE(I) * PI_OVER_180 SETCONA1.307
COS_P_LATITUDE(I) = COS(LATITUDE(I)) SETCONA1.308
SEC_P_LATITUDE(I) = 1.0 / COS_P_LATITUDE(I) SETCONA1.309
SETCONA1.310
*IF DEF,GLOBAL SETCONA1.311
TRUE_LONGITUDE(I) = LONGITUDE(I) SETCONA1.312
SETCONA1.313
F3_P(I) = 2.0 * OMEGA * SIN(LATITUDE(I)) SETCONA1.314
SETCONA1.315
*ELSE SETCONA1.316
TRUE_LONGITUDE(I) = TRUE_LONGITUDE(I) * PI_OVER_180 SETCONA1.317
SETCONA1.318
F3_P(I) = 2.0 * OMEGA * SETCONA1.319
& ( SIN(LATITUDE(I)) * SIN(POLE_LAT) SETCONA1.320
& + COS(LATITUDE(I)) * COS(LONGITUDE(I)) * COS(POLE_LAT) SETCONA1.321
& ) SETCONA1.322
SETCONA1.323
*ENDIF SETCONA1.324
SETCONA1.325
END DO SETCONA1.326
SETCONA1.327
*IF DEF,GLOBAL SETCONA1.328
C Values of COS_P_LATITUDE and SEC_P_LATITUDE at the polar points ACH1F400.3
C are used for area weighting and should reflect area of polar point. ACH1F400.4
C Geometrically this is 0.125 area of equatorward-lying neighbours. ACH1F400.5
C (But note that HADCM2 assumes 0.25 for consistency with old code). ATJ0F403.14
*IF -DEF,MPP GPB0F305.355
DO I = 1, ROW_LENGTH SETCONA1.330
ATJ0F403.15
IF (L_OLD_PWTS) THEN ATJ0F403.16
COS_P_LATITUDE(I) = 0.25 * COS_P_LATITUDE(I + ROW_LENGTH) ATJ0F403.17
SEC_P_LATITUDE(I) = 4.0 * SEC_P_LATITUDE(I + ROW_LENGTH) ATJ0F403.18
COS_P_LATITUDE(P_FIELD + I - ROW_LENGTH) = ATJ0F403.19
& 0.25 * COS_P_LATITUDE(P_FIELD + I - 2 * ROW_LENGTH) ATJ0F403.20
SEC_P_LATITUDE(P_FIELD + I - ROW_LENGTH) = ATJ0F403.21
& 4.0 * SEC_P_LATITUDE(P_FIELD + I - 2 * ROW_LENGTH) ATJ0F403.22
ELSE ATJ0F403.23
ATJ0F403.24
COS_P_LATITUDE(I) = 0.125 * COS_P_LATITUDE(I + ROW_LENGTH) ACH1F400.6
SEC_P_LATITUDE(I) = 8.0 * SEC_P_LATITUDE(I + ROW_LENGTH) ACH1F400.7
COS_P_LATITUDE(P_FIELD + I - ROW_LENGTH) = SETCONA1.333
& 0.125 * COS_P_LATITUDE(P_FIELD + I - 2 * ROW_LENGTH) ACH1F400.8
SEC_P_LATITUDE(P_FIELD + I - ROW_LENGTH) = SETCONA1.335
& 8.0 * SEC_P_LATITUDE(P_FIELD + I - 2 * ROW_LENGTH) ACH1F400.9
ATJ0F403.25
ENDIF ATJ0F403.26
SETCONA1.337
END DO SETCONA1.338
SETCONA1.339
CL set up Fourier filtering array using subroutine SETFILT SETCONA1.340
CALL SETFILT(TRIGS,IFAX,ROW_LENGTH,P_ROWS,TWO_D_GRID_CORRECTION, GPB0F401.556
& P_FIELD,COS_P_LATITUDE) GPB0F401.557
*ELSE GPB0F305.356
IF (attop) THEN GPB0F305.357
DO I = 1+Offy*ROW_LENGTH, (Offy+1)*ROW_LENGTH GPB0F305.358
ATJ0F403.27
IF (L_OLD_PWTS) THEN ATJ0F403.28
COS_P_LATITUDE(I) = 0.25 * COS_P_LATITUDE(I + ROW_LENGTH) ATJ0F403.29
SEC_P_LATITUDE(I) = 4.0 * SEC_P_LATITUDE(I + ROW_LENGTH) ATJ0F403.30
ELSE ATJ0F403.31
ATJ0F403.32
COS_P_LATITUDE(I) = 0.125 * COS_P_LATITUDE(I + ROW_LENGTH) GPB0F401.558
SEC_P_LATITUDE(I) = 8.0 * SEC_P_LATITUDE(I + ROW_LENGTH) GPB0F401.559
ATJ0F403.33
ENDIF ATJ0F403.34
ATJ0F403.35
ENDDO GPB0F305.361
ENDIF GPB0F305.362
IF (atbase) THEN GPB0F305.363
DO I = P_FIELD-((Offy+1)*ROW_LENGTH)+1, GPB0F305.364
& P_FIELD-(Offy*ROW_LENGTH)+1 GPB0F305.365
ATJ0F403.36
IF (L_OLD_PWTS) THEN ATJ0F403.37
COS_P_LATITUDE(I) = 0.25 * COS_P_LATITUDE(I - ROW_LENGTH) ATJ0F403.38
SEC_P_LATITUDE(I) = 4.0 * SEC_P_LATITUDE(I - ROW_LENGTH) ATJ0F403.39
ELSE ATJ0F403.40
ATJ0F403.41
COS_P_LATITUDE(I) = 0.125 * COS_P_LATITUDE(I - ROW_LENGTH) GPB0F401.560
SEC_P_LATITUDE(I) = 8.0 * SEC_P_LATITUDE(I - ROW_LENGTH) GPB0F401.561
ATJ0F403.42
ENDIF ATJ0F403.43
ATJ0F403.44
ENDDO GPB0F305.368
ENDIF GPB0F305.369
GPB0F305.370
! set up Fourier filtering array using subroutine SETFILT GPB0F401.562
CALL SETFILT(TRIGS,IFAX,glsize(1),P_ROWS,TWO_D_GRID_CORRECTION, GPB0F401.563
& P_FIELD,COS_P_LATITUDE) GPB0F401.564
*ENDIF GPB0F305.375
SETCONA1.343
*ENDIF SETCONA1.344
SETCONA1.345
CL SETCONA1.346
CL 4. Calculate constants depending on longitude only SETCONA1.347
CL SETCONA1.348
DO I = 1, ROW_LENGTH SETCONA1.349
SIN_LONGITUDE(I) = SIN(LONGITUDE(I)) SETCONA1.350
COS_LONGITUDE(I) = COS(LONGITUDE(I)) SETCONA1.351
SETCONA1.352
END DO SETCONA1.353
SETCONA1.354
CL Set list of land points for indexing. Initialise indices to MDIs. ARR0F403.34
C ARR0F403.35
DO I=1,LAND_FIELD ARR0F403.36
LAND_LIST(I) = IMDI ARR0F403.37
ENDDO ARR0F403.38
ARR0F403.39
I=0 GSS9F402.157
DO II=1,P_FIELD GSS9F402.158
IF(LAND(II))THEN GSS9F402.159
I = I + 1 GSS9F402.160
LAND_LIST(I) = II GSS9F402.161
END IF GSS9F402.162
END DO GSS9F402.163
C GSS9F402.164
IF (I .NE. LAND_FIELD) THEN GPB0F401.6
ICODE = 3 GPB0F401.7
CMESSAGE = 'SETCONA: Wrong number of land points' GPB0F401.8
WRITE(6,*) 'I=',I,'LAND_FIELD=',LAND_FIELD GIE0F403.619
GOTO 999 GPB0F401.10
ENDIF GPB0F401.11
SETCONA1.355
CL 5 Set additional soil and vegitation fields from those held in dummy SETCONA1.356
C SOILB parameter SETCONA1.357
DO I = 1, LAND_FIELD SETCONA1.358
IF (C_EAG(I) .LT. 0.0) THEN SETCONA1.359
ICODE = 1 SETCONA1.360
CMESSAGE = 'SETCONA: -ve numberas exponential in SOILB' SETCONA1.361
SETCONA1.362
IF (C_EAG(I) .EQ. RMDI) THEN SETCONA1.363
ICODE = 2 SETCONA1.364
CMESSAGE = 'SETCONA: missing data as exponential in SOILB' SETCONA1.365
GOTO 999 SETCONA1.366
SETCONA1.367
END IF SETCONA1.368
GOTO 999 SETCONA1.369
SETCONA1.370
END IF SETCONA1.371
IF (SATCON(I).GT.0.0) THEN AJS1F401.153
SOILB(I) = SATCON(I) / ((SMVCSAT(I) - SMVCWT(I))**C_EAG(I)) AJS1F401.154
ELSE AJS1F401.155
SOILB(I) = 0.0 AJS1F401.156
END IF AJS1F401.157
SETCONA1.373
END DO SETCONA1.374
SETCONA1.375
SETCONA1.386
CL Convert from conserved thermodynamic variables to potential SETCONA1.387
CL temperature,specific humidity, cloud water,and cloud ice. SETCONA1.388
DO LEVEL = 1, Q_LEVELS SETCONA1.389
DO I = 1, P_FIELD SETCONA1.390
PU = PSTAR(I) * BKH(LEVEL+1) + AKH(LEVEL+1) SETCONA1.391
PL = PSTAR(I) * BKH(LEVEL) + AKH(LEVEL) SETCONA1.392
THETA(I,LEVEL) = THETA(I,LEVEL) * SETCONA1.393
& P_EXNER_C(P_EXNER(I,LEVEL+1),P_EXNER(I,LEVEL),PU,PL,KAPPA) SETCONA1.394
SETCONA1.395
END DO SETCONA1.396
END DO SETCONA1.397
SETCONA1.398
! ADM2F404.240
! ASK1F405.56
IF (L_RHCPT) THEN ASK1F405.57
! Set data field with care - do not wish to use first row of P_FIELD at ASK1F405.58
! N.Pole, nor last row at S.Pole. ASK1F405.59
FIRST_POINT = FIRST_VALID_PT ASK1F405.60
LAST_POINT = LAST_P_VALID_PT ASK1F405.61
POINTS = LAST_POINT-FIRST_POINT+1 ASK1F405.62
JS = FIRST_POINT-1 ASK1F405.63
! ASK1F405.64
! Calculate critical relative humidity for all grid points ASK1F405.65
CALL RHCRIT_CALC(AK,BK,AKH,BKH,PSTAR(FIRST_POINT), ASK1F405.66
& RHCPT(FIRST_POINT,1),Q_LEVELS,POINTS,P_FIELD, ASK1F405.67
& THETA(FIRST_POINT,1),Q(FIRST_POINT,1), ASK1F405.68
& QCF(FIRST_POINT,1),ROW_LENGTH,LAND(FIRST_POINT), ASK1F405.69
& ICE_FRAC(FIRST_POINT),BL_LEVELS) ASK1F405.70
! ASK1F405.71
*IF DEF,MPP ASK1F405.72
CALL SWAPBOUNDS(RHCPT,ROW_LENGTH,P_ROWS,EW_Halo,NS_Halo, ASK1F405.73
& Q_LEVELS) ASK1F405.74
*ENDIF ASK1F405.75
! ASK1F405.76
ENDIF ! L_RHCPT ASK1F405.77
! ASK1F405.78
CALL GLUE_CLD(AK, BK, PSTAR, RHCRIT, Q_LEVELS, RHCPT, ASK1F405.79
& P_FIELD, P_FIELD, THETA, CLOUD_FRACTION, Q, QCF, QCL, AYY2F400.51
& LS_GRID_QC, LS_BS, ICODE) AYY2F400.52
SETCONA1.401
IF (ICODE .GT. 0)THEN SETCONA1.402
CMESSAGE = 'SETCONA : Error in LS_CLD' SETCONA1.403
GOTO 999 SETCONA1.404
SETCONA1.405
END IF SETCONA1.406
SETCONA1.407
DO LEVEL = 1, Q_LEVELS SETCONA1.408
DO I = 1, P_FIELD SETCONA1.409
PU = PSTAR(I) * BKH(LEVEL+1) + AKH(LEVEL+1) SETCONA1.410
PL = PSTAR(I) * BKH(LEVEL) + AKH(LEVEL) SETCONA1.411
THETA(I,LEVEL) = THETA(I,LEVEL) / SETCONA1.412
& P_EXNER_C(P_EXNER(I,LEVEL+1),P_EXNER(I,LEVEL),PU,PL,KAPPA) SETCONA1.413
SETCONA1.414
END DO SETCONA1.415
END DO SETCONA1.416
CL SETCONA1.417
CL 6 Set up cloud type boundaries for low/medium/high SETCONA1.418
CL SETCONA1.419
IF (NUM_CLOUD_TYPES .GT. 3) THEN SETCONA1.420
ICODE = 1 SETCONA1.421
CMESSAGE = 'SETCONA: Parameter NUM_CLOUD_TYPES exceeds 3' SETCONA1.422
WRITE(6,*) 'NUM_CLOUD_TYPES=',NUM_CLOUD_TYPES GIE0F403.620
GOTO 999 SETCONA1.424
SETCONA1.425
END IF SETCONA1.426
SETCONA1.427
C First calculate ETA SETCONA1.428
DO II = 1, P_LEVELS SETCONA1.429
ETA(II) = AK(II) / PREF + BK(II) SETCONA1.430
SETCONA1.431
END DO SETCONA1.432
SETCONA1.433
C The cloud amount is found by finding the max cloud cover over SETCONA1.434
C a set of levels. The boundaries of these types are defined by SETCONA1.435
C their ETA values as specified in ETA_SPLIT.The next piece of code SETCONA1.436
C finds the levels (layer mid points) which fall between the SETCONA1.437
C specfied ETA values.Thus the layers contributing to the low cloud SETCONA1.438
C have ETA values between ETA_SPLIT(1) and ETA_SPLIT(2) SETCONA1.439
SETCONA1.440
DO KK = 1, NUM_CLOUD_TYPES + 1 SETCONA1.441
LEVEL = 1 SETCONA1.442
SETCONA1.443
DO WHILE ((ETA(LEVEL) .GE. ETA_SPLIT(KK)) .AND. SETCONA1.444
& (LEVEL .LE. P_LEVELS)) SETCONA1.445
LEVEL = LEVEL + 1 SETCONA1.446
SETCONA1.447
END DO SETCONA1.448
SETCONA1.449
IF (LEVEL .GT. P_LEVELS) THEN SETCONA1.450
ICODE = 1 SETCONA1.451
CMESSAGE ='SETCONA:Error in locating levels for cloud lyrs' SETCONA1.452
GOTO 999 SETCONA1.453
SETCONA1.454
END IF SETCONA1.455
CLOUD_BOUND(KK) = LEVEL SETCONA1.456
SETCONA1.457
END DO SETCONA1.458
SETCONA1.459
LOW_BOT_LEVEL = CLOUD_BOUND(1) SETCONA1.460
LOW_TOP_LEVEL = CLOUD_BOUND(2) - 1 SETCONA1.461
MED_BOT_LEVEL = CLOUD_BOUND(2) SETCONA1.462
MED_TOP_LEVEL = CLOUD_BOUND(3) - 1 SETCONA1.463
HIGH_BOT_LEVEL = CLOUD_BOUND(3) SETCONA1.464
HIGH_TOP_LEVEL = CLOUD_BOUND(4) - 1 SETCONA1.465
SETCONA1.466
IF (LOW_TOP_LEVEL .GT. CLOUD_LEVELS) THEN SETCONA1.467
ICODE = 1 SETCONA1.468
CMESSAGE = 'SETCONA: No of cloud levels less than Top of Low' SETCONA1.469
GOTO 999 SETCONA1.470
SETCONA1.471
END IF SETCONA1.472
SETCONA1.473
IF (MED_TOP_LEVEL.GT.CLOUD_LEVELS) THEN SETCONA1.474
ICODE = 1 SETCONA1.475
CMESSAGE = 'SETCONA: No of cloud levels less than Top of Med' SETCONA1.476
GOTO 999 SETCONA1.477
SETCONA1.478
END IF SETCONA1.479
SETCONA1.480
IF (HIGH_TOP_LEVEL .GT. CLOUD_LEVELS) THEN SETCONA1.481
ICODE = 1 SETCONA1.482
CMESSAGE = 'SETCONA: No of cloud levels less than Top of High' SETCONA1.483
GOTO 999 SETCONA1.484
SETCONA1.485
END IF SETCONA1.486
SETCONA1.487
CL 7. Calculate ETA_HALF inverse matrix SETCONA1.488
DO LEVEL = 1, P_LEVELS + 1 SETCONA1.489
ETA_HALF(LEVEL) = AKH(LEVEL) / PREF + BKH(LEVEL) SETCONA1.490
SETCONA1.491
END DO SETCONA1.492
N21 = (MATRIX_POLY_ORDER + 1) / 2 SETCONA1.493
SETCONA1.494
DO LEVEL = N21, P_LEVELS - N21 + 1 SETCONA1.495
DO POINT = 1, MATRIX_POLY_ORDER SETCONA1.496
! NB ROW used as a work variable SETCONA1.497
ROW = LE VEL - N21 + POINT SETCONA1.498
DO I = 1, MATRIX_POLY_ORDER SETCONA1.499
ETA_MATRIX(I,POINT,LEVEL)=(ETA_HALF(ROW)**I-ETA_HALF(ROW+1) SETCONA1.500
& **I)/((ETA_HALF(ROW)-ETA_HALF(ROW+1))*I) SETCONA1.501
SETCONA1.502
END DO SETCONA1.503
END DO SETCONA1.504
CALL MATINV(ETA_MATRIX(1,1,LEVEL), ETA_MATRIX_INV(1,1,LEVEL), SETCONA1.505
& MATRIX_POLY_ORDER) SETCONA1.506
SETCONA1.507
END DO SETCONA1.508
SETCONA1.509
CL 8. Calculate LEVNO_PMSL_CALC for PMSL reduction calculations. GDR3F404.15
GDR3F404.16
DO LEVEL=1,P_LEVELS GDR3F404.17
ETA_LEV=AK(LEVEL)/PREF + BK(LEVEL) GDR3F404.18
IF (ETA_LEV.LT.ETA_PMSL) THEN GDR3F404.19
LEVNO_PMSL_CALC=LEVEL GDR3F404.20
WRITE (6,*) GDR3F404.21
& ' SETCONA : LEVEL FOR PMSL REDUCTION CALCULATION = ', GDR3F404.22
& LEVNO_PMSL_CALC GDR3F404.23
GO TO 555 GDR3F404.24
ENDIF GDR3F404.25
ENDDO GDR3F404.26
555 CONTINUE GDR3F404.27
GDR3F404.28
999 RETURN SETCONA1.510
END SETCONA1.511
SETCONA1.512
*ENDIF SETCONA1.513