SUBROUTINE MIXLAY (TA, WME, 2,1MIXLAY.51
+ IMT, KM, NT, MIXLAY.52
& J,JMT, ORH7F404.52
+ C2DTTS, GRAV_SI, MIXLAY.53
+ DZ, ZDZ, RZ, MIXLAY.54
+ FM, KMT, MIXLAY.55
+ LAMDA, EPSILON, MIXLAY.56
+ DELPSF, DELPSL, DECAY, MIXLAY.57
+ DRHO_SOL, DRHO_NET, MIXLAY.58
+ DRHO_ICE, MIXLAY.61
& DRHO_PSEU, ORH1F305.229
+ DIAG_MLD MIXLAY.67
+ ) MIXLAY.68
C MIXLAY.69
IMPLICIT NONE MIXLAY.70
C MIXLAY.71
C Define constant for array sizes MIXLAY.72
*CALL CNTLOCN 
ORH1F305.230
*CALL OARRYSIZ ORH1F305.231
C MIXLAY.73
C MIXLAY.74
INTEGER MIXLAY.75
+ IMT ! IN Number of points in first dimension of arrays MIXLAY.76
+, KM ! IN Number of points in vertical MIXLAY.77
+, NT ! IN Number of tracers MIXLAY.78
C MIXLAY.79
C MIXLAY.80
INTEGER MIXLAY.81
& KMT (IMT) ! IN Number of gridpoints in column ORH1F305.232
&,JMT ! IN No of rows JA121293.226
&,J ! IN ROW NUMBER JA121293.227
C MIXLAY.83
C Physical arguments MIXLAY.84
C MIXLAY.85
REAL MIXLAY.86
+ C2DTTS ! IN Timestep MIXLAY.87
+, TA (IMT,KM,NT) ! INOUT Array of tracers MIXLAY.88
+, WME (IMT) ! IN Wind mixing power (W/m2) MIXLAY.89
+, DZ (KM) ! IN Model layer thicknesses MIXLAY.90
+, ZDZ (KM) ! IN Depth of base of layers MIXLAY.91
+, RZ (KM) ! IN Layer thicknesses scaled by 1/RAT(K) MIXLAY.92
C Note that if variable timestep with MIXLAY.93
C depth not included, this is speci- MIXLAY.94
C fied as DZ in the call MIXLAY.95
+, FM (IMT, KM) ! IN Land/sea mask: 1.0=sea MIXLAY.96
+, LAMDA ! IN Reduced WME passing through surface MIXLAY.97
+, EPSILON ! IN Reduced energy after convection MIXLAY.98
+, DELPSF ! IN 'Non solar density change' (cm3/s2) MIXLAY.99
&, DELPSL (IMT,0:KM) ! IN Solar density change (cm3/s2) OJP0F404.618
+, DECAY (KM) ! IN Decay rate of WME with depth MIXLAY.101
+, DRHO_SOL (IMT) ! IN Density change at sfce: solar (g/cm3 MIXLAY.102
+, DRHO_NET (IMT) ! IN Density change at sfce: net heat (g/ MIXLAY.103
&,DRHO_ICE (IMT_ICE_MIX) ! IN Density changes at sfce: sea ice (g/ ORH1F305.237
&,DRHO_PSEU (IMT_PIC_MIX)! IN Density changes at sfce: pseudice (g ORH1F305.238
+, DIAG_MLD (IMT) ! OUT Instantaneous mixed layer depth (m) MIXLAY.112
+, GRAV_SI ! IN Acceleration due to gravity (m/s2) MIXLAY.113
C* MIXLAY.114
C*L---------- External subprograms called ------------------------------ MIXLAY.115
EXTERNAL MIXLAY.116
+ STATED ! Potential density referenced to sfc MIXLAY.117
C*---------------------------------------------------------------------- MIXLAY.118
C MIXLAY.119
C Locally defined variables and arrays MIXLAY.120
C MIXLAY.121
C MIXLAY.124
INTEGER MIXLAY.125
+ i ! Zonal index for loops MIXLAY.126
+, k ! Vertical index for loops MIXLAY.127
+, m ! Tracer number in loops MIXLAY.128
+, imixlay ! Number of points active MIXLAY.129
+, kmm1 ! KM - 1 MIXLAY.130
+, indmix(IMT) ! Compression index for mixed layer MIXLAY.131
C MIXLAY.135
REAL MIXLAY.136
+ pdz (KM) ! Weighted layer thickness MIXLAY.137
+, eqvdep (KM) ! Depth at bottom of layer (weighted) MIXLAY.138
+, wmxe (IMT) ! Depleted wind mixing energy MIXLAY.139
+, rhoref (IMT) ! Mixed layer density MIXLAY.140
+, tracer_ref (IMT, NT) ! Mixed layer tracer values MIXLAY.141
+, rho_ref (IMT) ! Mixed layer density MIXLAY.142
+, tracer_up (IMT, NT) ! Compressed mixed layer tracer MIXLAY.143
+, rho_up (IMT) ! Compressed mixed layer density MIXLAY.144
+, tracer_down (IMT,NT) ! Compressed sub-mixing tracer MIXLAY.145
+, rho_down (IMT) ! Compressed sub-mixing density MIXLAY.146
REAL MIXLAY.147
+ reng (IMT) ! Compressed energy for mixing MIXLAY.148
+, enrat_up (IMT) ! Temporary store for partial mixing MIXLAY.149
+, enrat_down (IMT) ! Temporary store for partial mixing MIXLAY.150
+, temp_mld (IMT) ! Compressed mixed layer diagnostic MIXLAY.151
+, pe (IMT) ! Potential energy difference (comp) MIXLAY.152
+, c1, c2, c3, ca, cz ! Temporary scalars MIXLAY.153
+, f1,f2 ! Temporary scalars MIXLAY.154
+, missing ! Missing data indicator MIXLAY.155
+, grav ! CGS version of GRAV_SI MIXLAY.156
C MIXLAY.157
C Workspace arrays (large arrays to permit re-use of space from MIXLAY.158
C TRACER) MIXLAY.159
C MIXLAY.160
REAL MIXLAY.161
+ worka (IMT, KM) ! Used for density calculation MIXLAY.162
+, workb (IMT, KM) ! Used for density calculation MIXLAY.163
+, rho (IMT, KM) ! Density MIXLAY.164
C MIXLAY.165
C Logical arrays MIXLAY.166
C MIXLAY.167
LOGICAL MIXLAY.168
+ mixnext (IMT) ! Mixing possible MIXLAY.169
+, fullmix (IMT) ! Full or partial layer mixing MIXLAY.170
C MIXLAY.171
C Integer arrays MIXLAY.172
C MIXLAY.173
INTEGER MIXLAY.174
+ maxmix (IMT) ! Level of maximum mixing MIXLAY.175
C MIXLAY.176
C MIXLAY.177
C Convert gravity to cgs MIXLAY.178
C MIXLAY.179
grav=GRAV_SI * 100. MIXLAY.180
C MIXLAY.181
CL 1.1 Define value of 'missing' (this is purely local) MIXLAY.182
C MIXLAY.183
missing = -1.0E30 MIXLAY.184
C MIXLAY.185
CL 1.2 Vertical constant arrays MIXLAY.186
C MIXLAY.187
DO 1200, k = 1, km MIXLAY.188
pdz (k) = RZ (k) MIXLAY.189
1200 CONTINUE MIXLAY.190
eqvdep (1) = pdz (1) MIXLAY.191
DO 1210, k = 2, km MIXLAY.192
eqvdep (k) = eqvdep (k-1) + pdz (k) MIXLAY.193
1210 CONTINUE MIXLAY.194
C MIXLAY.195
CL 1.3 Initial densities MIXLAY.196
C MIXLAY.197
JA121293.233
ORH1F305.241
CALL STATED(TA(1,1,1), TA(1,1,2), rho, worka, workb, IMT, KM,J JA121293.243
&,KM JA121293.244
&,JMT ORH7F404.53
&) JA121293.247
JA121293.248
C MIXLAY.199
CL 1.4 Default mixed layer depth MIXLAY.208
C MIXLAY.209
do 1410, i = 1,IMT @DYALLOC.4647
DIAG_MLD(i) = FM(i, 1) * ZDZ(1)*0.01 MIXLAY.227
1410 CONTINUE MIXLAY.228
C MIXLAY.229
CL 1.5 Depleteion of mixing energy as it crosses surface MIXLAY.230
C Note: 1000 converts SI WME to cgs MIXLAY.231
C MIXLAY.232
do 1510, i = 1,IMT @DYALLOC.4648
wmxe(i) = wme(i) * (C2DTTS*LAMDA*1000.0) MIXLAY.234
1510 CONTINUE MIXLAY.235
C MIXLAY.236
CL 1.6 Initialise mixed layer tracer values and pointers MIXLAY.237
C MIXLAY.238
DO 1620, m = 1, NT MIXLAY.239
do 1610, i = 1,IMT @DYALLOC.4649
tracer_ref (i, m) = TA (i, 1, m) MIXLAY.241
1610 CONTINUE MIXLAY.242
1620 CONTINUE MIXLAY.243
C MIXLAY.244
do 1630, i = 1,IMT @DYALLOC.4650
C MIXLAY.246
CL Density MIXLAY.247
C MIXLAY.248
rho_ref(i) = rho(i, 1) MIXLAY.249
C MIXLAY.250
CL Pointers for active mixing MIXLAY.251
C MIXLAY.252
mixnext(i) = FM(i, 2) .EQ. 1.0 MIXLAY.253
C MIXLAY.254
CL Pointer for maximum level mixed MIXLAY.255
maxmix(i)=1 MIXLAY.256
C MIXLAY.257
CL 1.7 Deplete initial energy due to surface layer mixing MIXLAY.258
C MIXLAY.259
IF ( DRHO_SOL(i).LT.0.0 ) THEN MIXLAY.260
c1 = DRHO_SOL(i) * DELPSL(i,KMT(i)) ORH1F305.247
ELSE MIXLAY.262
c1 = 0.0 MIXLAY.263
END IF MIXLAY.264
IF (L_SEAICE) THEN ORH1F305.249
c1 = c1 + DELPSF*DRHO_ICE(i) ORH1F305.250
ENDIF ORH1F305.251
IF (L_OPSEUDIC) THEN ORH1F305.252
c1 = c1 + DELPSF*DRHO_PSEU(i) ORH1F305.253
ENDIF ORH1F305.254
c1 = (DELPSF*DRHO_NET(i) + c1) * FM(i,1) MIXLAY.273
wmxe(i) = wmxe(i) - c1 MIXLAY.274
1630 CONTINUE MIXLAY.275
12000 CONTINUE MIXLAY.276
C MIXLAY.277
CL 2.0 Loop over layers to perform mixing MIXLAY.278
C MIXLAY.279
kmm1 = KM - 1 MIXLAY.289
DO 2999, k = 1, kmm1 MIXLAY.290
C MIXLAY.291
CL 2.1 Mixing coefficients MIXLAY.292
C MIXLAY.293
cz = 1.0 / eqvdep (K+1) MIXLAY.294
f1 = ZDZ(k+1)*ZDZ(k+1)*eqvdep(k) - MIXLAY.295
- ZDZ(k)*ZDZ(k)*eqvdep(k+1) MIXLAY.296
f2 = pdz(k+1)*ZDZ(k+1)*ZDZ(k+1) - MIXLAY.297
- (ZDZ(k)+ZDZ(k+1))*DZ(k+1)*eqvdep(k+1) MIXLAY.298
f1 = -0.5*f1*cz*grav MIXLAY.299
f2 = -0.5*f2*cz*grav MIXLAY.300
C MIXLAY.301
CL 2.2 Indicator for active points and compression index MIXLAY.302
C MIXLAY.303
do 2210, i = 1,IMT @DYALLOC.4651
mixnext (i) = (FM(i,k+1).EQ.1.0) .AND. MIXLAY.305
& (wmxe(i).GT.0.0) MIXLAY.306
2210 CONTINUE MIXLAY.307
C MIXLAY.308
C Set up index array in indmix: number of active points in imixlay MIXLAY.309
C MIXLAY.310
imixlay = 0 MIXLAY.311
do 2220, i = 1,IMT @DYALLOC.4652
IF ( mixnext(i) ) THEN MIXLAY.313
imixlay = imixlay + 1 MIXLAY.314
indmix(imixlay) = i MIXLAY.315
END IF MIXLAY.316
2220 CONTINUE MIXLAY.317
C MIXLAY.318
CL 2.3 Compress arrays to active points MIXLAY.319
C MIXLAY.320
IF ( imixlay .eq. 0) GO TO 3000 ! Stop loop if no points MIXLAY.321
C MIXLAY.322
C Tracers MIXLAY.323
C MIXLAY.324
DO 2320, m = 1, NT MIXLAY.325
DO 2310, i = 1, imixlay MIXLAY.326
tracer_up(i, m) = tracer_ref(indmix(i), m) MIXLAY.327
tracer_down(i, m) = TA(indmix(i), k+1, m) MIXLAY.328
2310 CONTINUE MIXLAY.329
2320 CONTINUE MIXLAY.330
C MIXLAY.331
C Density and wind mixing energy MIXLAY.332
C MIXLAY.333
DO 2330, i = 1, imixlay MIXLAY.334
rho_up(i) = rho_ref(indmix(i)) MIXLAY.335
rho_down(i) = rho(indmix(i), k+1) MIXLAY.336
reng(i) = wmxe(indmix(i)) MIXLAY.337
2330 CONTINUE MIXLAY.338
C MIXLAY.339
CL 2.4 Deplete wind mixing energy and calc P.E. difference MIXLAY.340
C MIXLAY.341
DO 2410, i = 1, imixlay MIXLAY.342
reng(i) = reng(i)*DECAY(k) MIXLAY.343
pe(i) = f1*rho_up(i) + f2*rho_down(i) MIXLAY.344
2410 CONTINUE MIXLAY.345
C MIXLAY.346
C Deplete energy due to convection MIXLAY.347
C MIXLAY.348
DO 2420, i = 1, imixlay MIXLAY.349
IF (pe(i) .LT. 0.0) pe(i) = pe(i)*EPSILON MIXLAY.350
2420 CONTINUE MIXLAY.351
C MIXLAY.352
CL 2.5 Mix or partially mix layers MIXLAY.353
C MIXLAY.354
c1 = eqvdep(k)*cz MIXLAY.355
c2 = (1.0 - eqvdep(k)*cz) MIXLAY.356
c3 = 2.0 / (grav*ZDZ(k)) MIXLAY.357
MIXLAY.358
DO 2510, i = 1, imixlay MIXLAY.359
C MIXLAY.360
C Indicator for full or partial mixing MIXLAY.361
C MIXLAY.362
C Added test for rho_down(i).EQ.rho_up(i), which should be equivalant ONF0F401.3
C to pe(i).EQ.0.0, but sometimes isn't due to rounding error on a 32bit ONF0F401.4
C precision machine. N.Farnon (14/2/96) ONF0F401.5
fullmix(i) = pe(i).LE.reng(i) .OR. pe(i).EQ.0.0 .OR. ONF0F401.6
& rho_down(i).EQ.rho_up(i) ONF0F401.7
C MIXLAY.364
C Coefficients for partial mixing MIXLAY.365
C MIXLAY.366
IF ( pe(i).ne.0.0 ) THEN MIXLAY.367
ca = pe(i) MIXLAY.368
ELSE MIXLAY.369
ca = 1.0 MIXLAY.370
END IF MIXLAY.371
ca = reng(i)/ca MIXLAY.372
enrat_up(i) = c2*ca MIXLAY.373
enrat_down(i) = c1*ca MIXLAY.374
2510 CONTINUE MIXLAY.375
C MIXLAY.376
DO 2530, m = 1, NT MIXLAY.377
DO 2520, i = 1, imixlay MIXLAY.378
IF ( fullmix(i) ) THEN MIXLAY.379
C MIXLAY.380
CL Complete mixing of layer tracers MIXLAY.381
C MIXLAY.382
tracer_up(i,m) = MIXLAY.383
= (tracer_up(i, m)*eqvdep(k) + MIXLAY.384
+ tracer_down(i, m)*pdz(k+1))*cz MIXLAY.385
ELSE MIXLAY.386
C MIXLAY.387
CL Partial mixing of tracers MIXLAY.388
C MIXLAY.389
ca = tracer_up(i, m) - tracer_down(i, m) MIXLAY.390
tracer_up(i, m) = tracer_up(i, m) - MIXLAY.391
+ enrat_up(i)*ca MIXLAY.392
tracer_down(i, m) = tracer_down(i, m) + MIXLAY.393
+ enrat_down(i)*ca MIXLAY.394
END IF MIXLAY.395
2520 CONTINUE MIXLAY.396
2530 CONTINUE MIXLAY.397
C MIXLAY.398
DO 2540, i = 1, imixlay MIXLAY.399
IF ( fullmix(i) ) THEN MIXLAY.400
C MIXLAY.401
CL Complete mixing of density and calculate new energy & MLD MIXLAY.402
C MIXLAY.403
rho_up(i) = MIXLAY.404
= (rho_up(i)*eqvdep(k) + MIXLAY.405
+ rho_down(i)*pdz(k+1))*cz MIXLAY.406
C MIXLAY.407
reng(i) = reng(i) - pe(i) MIXLAY.408
C MIXLAY.409
temp_mld(i) = zdz(k+1) MIXLAY.410
C MIXLAY.411
CL Set indicator to show mixing reached this layer MIXLAY.412
C MIXLAY.413
maxmix(indmix(i)) = k+1 MIXLAY.414
ELSE MIXLAY.415
C MIXLAY.416
CL Partial mixing, calculate MLD and set energy MIXLAY.417
C MIXLAY.418
temp_mld(i) = MIXLAY.419
= zdz(k) + c3*(reng(i)/(rho_down(i)- MIXLAY.420
- rho_up(i))) MIXLAY.421
C MIXLAY.422
reng(i) = missing MIXLAY.423
END IF MIXLAY.424
2540 CONTINUE MIXLAY.425
C MIXLAY.426
CL 2.6 Return all points to geographical locations MIXLAY.427
C MIXLAY.428
DO 2610, m = 1, NT MIXLAY.429
DO 2600, i = 1, imixlay MIXLAY.430
tracer_ref(indmix(i), m) = tracer_up(i, m) MIXLAY.431
2600 CONTINUE MIXLAY.432
2610 CONTINUE MIXLAY.433
DO 2620, i = 1, imixlay MIXLAY.434
rho_ref(indmix(i)) = rho_up(i) MIXLAY.435
wmxe(indmix(i)) = reng(i) MIXLAY.436
2620 CONTINUE MIXLAY.437
C MIXLAY.438
CL 2.7 Where partial mixing set tracer values below mixed layer MIXLAY.439
C MIXLAY.440
DO 2710, m = 1, NT MIXLAY.441
DO 2700, i = 1, imixlay MIXLAY.442
C MIXLAY.443
C Note: reng = missing if partial mixing MIXLAY.444
C MIXLAY.445
IF ( reng(i) .EQ. missing) THEN MIXLAY.446
TA(indmix(i), k+1, m) = MIXLAY.447
= tracer_down(i, m) MIXLAY.448
END IF MIXLAY.449
2700 CONTINUE MIXLAY.450
2710 CONTINUE MIXLAY.451
C MIXLAY.452
CL 2.8 Set diagnostic depth (instantaneous) MIXLAY.453
C MIXLAY.454
DO 2800, i = 1, imixlay MIXLAY.455
diag_mld(indmix(i)) = temp_mld(i)*0.01 MIXLAY.456
2800 CONTINUE MIXLAY.457
2999 CONTINUE MIXLAY.458
C MIXLAY.459
CL 3.0 End of mixed layer, either because bottom or no more energy MIXLAY.460
C MIXLAY.461
3000 CONTINUE MIXLAY.462
C MIXLAY.463
CL 3.1 Copy mixed layer values to main arrays MIXLAY.464
C MIXLAY.465
DO 3120, m = 1, NT MIXLAY.466
do 3110, i = 1,IMT @DYALLOC.4653
DO 3100, k = 1, maxmix(i) MIXLAY.468
TA(i, k, m) = tracer_ref(i, m) MIXLAY.469
3100 CONTINUE MIXLAY.470
3110 CONTINUE MIXLAY.471
3120 CONTINUE MIXLAY.472
ORH1F305.255
IF (L_OCYCLIC) THEN ORH1F305.256
C MIXLAY.476
CL 3.2 Wrap-around points for diagnostic MIXLAY.477
C MIXLAY.478
diag_mld(1) = diag_mld(IMT-1) ORH1F305.257
diag_mld(IMT) = diag_mld(2) ORH1F305.258
ENDIF ORH1F305.259
ORH1F305.260
C MIXLAY.483
CL Return from mixed layer calculations MIXLAY.484
C MIXLAY.485
RETURN MIXLAY.486
END MIXLAY.487
*ENDIF MIXLAY.488