SUBROUTINE IPDCOFCO 4,4IPDCOFCO.58
& ( J,JMT,IMT,IMTM1,KM,KMT,KMP,KMP1,KMP2,NT,NTMIN2, IPDCOFCO.59
& J_OFFSET, ORH7F404.55
& T,TP,TDIF, IPDCOFCO.62
& DXUR,DXU2RQ,DXT4RQ,DYUR,DYT4R,DZ2RQ,DZZ2RQ,ZDZ,DTTS, IPDCOFCO.63
& NERGY,CSR,CSTR,ITT,FM,FMP, IPDCOFCO.64
& RHOS,RHON,AH,ahi, IPDCOFCO.65
& gnu,JFT0,esav,fk1,fk2,fk3, IPDCOFCO.66
& rxp,ry,rrzp IPDCOFCO.67
& ,mld IPDCOFCO.68
& ,slope_max IPDCOFCO.69
& ) IPDCOFCO.70
C IPDCOFCO.71
C IPDCOFCO.72
IMPLICIT NONE IPDCOFCO.73
C IPDCOFCO.74
*CALL OARRYSIZ 
PXORDER.24
INTEGER IPDCOFCO.75
& I,J,K,M,N,L, IPDCOFCO.76
& JMT,IMT,IMTM1,KM,KMT,KMP,KMP1,KMP2,NT,NTMIN2, IPDCOFCO.77
& NERGY,ITT,KS,JFT0 IPDCOFCO.78
&,J_OFFSET ! IN base row of block IPDCOFCO.83
C IPDCOFCO.84
REAL IPDCOFCO.85
& T(IMT,KM,NT),TP(IMT,KM,NT),TDIF(IMT,KMP2,NTMIN2), IPDCOFCO.86
& DXUR(IMT),DXU2RQ(IMT,KM),DXT4RQ(IMT,KM), IPDCOFCO.87
& DYUR(JMT),DYT4R(JMT),DZ2RQ(IMT,KM),DZZ2RQ(IMT,KM),ZDZ(KM),DTTS, IPDCOFCO.88
& RHOS(IMT,KM),RHON(IMT,KM), IPDCOFCO.89
& FM(IMT,KM),FMP(IMT,KM),CSR(JMT),CSTR(JMT), IPDCOFCO.90
& AH, ! IN BACKGROUND horizontal diffusivity IPDCOFCO.91
C Ah in Eq.1.10 (cm2/s) IPDCOFCO.92
& ahi(KM), ! IN along-isopycnal coeff. of diffusion IPDCOFCO.93
C Ai in Eq.1.8 (cm2/s) IPDCOFCO.94
& gnu(IMT,KM), ! IN coeff. of vertical diffusion at TOP of IPDCOFCO.95
C box. Av in Eq.1.10 (cm2/s) IPDCOFCO.96
& fk1(IMT,KM,3), ! OUT } Rows 1,2,3 of diffusion matrix (Eq1.8) IPDCOFCO.97
& fk2(IMT,KM,3), ! OUT } fk1 is on EAST face, fk2 on NORTH face IPDCOFCO.98
& fk3(IMT,KM,3), ! OUT } & fk3 on TOP face of grid box (i,j,k) IPDCOFCO.99
& rxp(IMT,KM), ! INOUT delta-rho in x dirn, row J+1 (E face) IPDCOFCO.100
& ry(IMT,KM), ! INOUT delta-rho in y dirn, row J (N face) IPDCOFCO.101
& rrzp(IMT,KMP1), ! INOUT delta-rho in z dirn, row J+1 (top face IPDCOFCO.102
& esav(IMT,KM,NT) ! OUT used to save e(I,K,2) in IPDCOFCO.103
C subroutine IPDFLXCL IPDCOFCO.104
& ,mld(IMT_IPD_MIX) ! IN mixed layer depth (m) IPDCOFCO.105
& ,slope_max ! IN Maximum slope for diffusion IPDCOFCO.106
C IPDCOFCO.107
C IPDCOFCO.108
C Declare local constants and arrays IPDCOFCO.109
C IPDCOFCO.110
REAL IPDCOFCO.111
& rx(IMT,KM), ! delta-rho in x-dir, row J (E face) IPDCOFCO.112
& rym(IMT,KM), ! " " y-dir, " J-1 (N face) IPDCOFCO.113
& rrz(IMT,KMP1), ! " " z-dir, " J (top face) IPDCOFCO.114
& e(IMT,KMP1,3), ! d-rho/d-x , d-rho/d-y , d-rho/d-z IPDCOFCO.115
& slmxr(KM), ! cst. used in stability check on IPDCOFCO.116
C d-rho-dz in Eq.1.24 IPDCOFCO.117
& tempa(IMT,KMP1), ! workspace IPDCOFCO.118
& tempb(IMT,KMP1), ! workspace IPDCOFCO.119
& fxa, ! } IPDCOFCO.120
& fxb, ! } IPDCOFCO.121
& fxc, ! } local csts. IPDCOFCO.122
& fxd, ! } IPDCOFCO.123
& fxe, ! } IPDCOFCO.124
& temp ! temporary variable IPDCOFCO.125
*CALL CNTLOCN IPDCOFCO.126
C* IPDCOFCO.128
C*L---- External subroutines called ------------------------------------ IPDCOFCO.129
EXTERNAL STATEC IPDCOFCO.130
C* IPDCOFCO.131
C IPDCOFCO.132
C -------------------------------------------------------------- IPDCOFCO.133
CL Initial conditions for the isopycnal diffusion code. IPDCOFCO.134
CL These involve working out delta-rho in the x and z directions IPDCOFCO.135
CL for the second row, so that the process for calculating IPDCOFCO.136
CL these quantities for the other rows can be initiated. IPDCOFCO.137
C -------------------------------------------------------------- IPDCOFCO.138
C IPDCOFCO.139
DO 502 N=1,3 IPDCOFCO.140
DO 502 K=1,KM IPDCOFCO.141
DO 502 I=1,IMT IPDCOFCO.142
fk1(I,K,N)=0. IPDCOFCO.143
fk2(I,K,N)=0. IPDCOFCO.144
fk3(I,K,N)=0. IPDCOFCO.145
502 CONTINUE IPDCOFCO.146
fxa=1.E-25 IPDCOFCO.147
fxb=2. IPDCOFCO.148
fxc=.5 IPDCOFCO.149
fxd=1. IPDCOFCO.150
fxe=1.E10 ! Big number used to multiply all density gradients. IPDCOFCO.151
C Always cancels out in the end, so must be IPDCOFCO.152
C something to do with reducing rounding error. IPDCOFCO.153
C IPDCOFCO.154
DO 503 K=1,KM IPDCOFCO.155
C IPDCOFCO.156
C Set maximum allowable slope for the diffusion to slope_max. IPDCOFCO.157
C Also included here (commented out) are two old versions. IPDCOFCO.158
C The first value of temp implements the criterion in Cox IPDCOFCO.159
C (Ocean modelling #74, 1987). I think the 4. should be a IPDCOFCO.160
C 1. (RAW 6/5/92), but in any case this criterion gave IPDCOFCO.161
C poor results in a 2.5x3.75 GCM. The second, hardwired value IPDCOFCO.162
C of temp has worked in that GCM and in the 1x1 FOAM model. IPDCOFCO.163
C In either of these cases note that slmxr contains 2*temp, IPDCOFCO.164
C but that this is cancelled by a division by 2*dzz later on. IPDCOFCO.165
C IPDCOFCO.166
C In all cases, the maximum slope through which the diffusion IPDCOFCO.167
C tensor can be rotated is set to 2*DZ(k)/slmxr(k). IPDCOFCO.168
C IPDCOFCO.169
C temp=4.*DTTS*ahi(K) ! Cox's criterion IPDCOFCO.170
C temp=2.e13 ! Oscar's FOAM value IPDCOFCO.171
C slmxr(K)=(temp+temp)* IPDCOFCO.172
! For L_OFILTER = true -> MAX(DYUR(J),DXUR(1)*MIN(CSTR(J),CSTR(JFT0))) IPDCOFCO.173
! For L_OFILTER = false-> MAX(DYUR(J),DXUR(1)*CSTR(J)) IPDCOFCO.174
! IPDCOFCO.175
C IPDCOFCO.176
slmxr(K)=1./slope_max/DZ2RQ(2,K) ! DZ2RQ(2,.) chosen arbitraril IPDCOFCO.177
503 CONTINUE IPDCOFCO.178
C IPDCOFCO.179
IF(J+J_OFFSET.EQ.2) THEN IPDCOFCO.180
DO 505 K=1,KM IPDCOFCO.181
DO 505 I=1,IMT IPDCOFCO.182
ry(I,K)=0. IPDCOFCO.183
505 CONTINUE IPDCOFCO.184
DO 508 K=1,KMP1 IPDCOFCO.185
DO 508 I=1,IMT IPDCOFCO.186
rrzp(I,K)=0. IPDCOFCO.187
508 CONTINUE IPDCOFCO.188
DO 510 K=1,KM IPDCOFCO.189
DO 509 I=1,IMTM1 IPDCOFCO.190
rxp(I,K)=FM(I,K)*FM(I+1,K)*(RHOS(I+1,K)-RHOS(I,K))*fxe IPDCOFCO.191
509 CONTINUE IPDCOFCO.192
IF (L_OCYCLIC) THEN IPDCOFCO.193
rxp(IMT,K)=rxp(2,K) IPDCOFCO.194
ELSE IPDCOFCO.195
rxp(IMT,K)=0.0 IPDCOFCO.196
ENDIF IPDCOFCO.197
510 CONTINUE IPDCOFCO.198
DO 777 K=1,KMP1 IPDCOFCO.199
DO 777 I=1,IMT IPDCOFCO.200
tempa(I,K)=0. ! Initialise tempa IPDCOFCO.201
tempb(I,K)=0. ! Initialise tempb IPDCOFCO.202
777 CONTINUE IPDCOFCO.203
IPDCOFCO.204
IPDCOFCO.205
CALL STATEC(T,T(1,1,2),tempa,TDIF,TDIF(1,1,2),1,IMT,KM,2 IPDCOFCO.206
&,JMT ORH7F404.56
&) IPDCOFCO.209
CALL STATEC(T,T(1,1,2),tempb,TDIF,TDIF(1,1,2),2,IMT,KM,2 IPDCOFCO.210
&,JMT ORH7F404.57
&) IPDCOFCO.213
IPDCOFCO.214
IPDCOFCO.215
DO 778 I=1,IMT IPDCOFCO.216
tempa(I,KMP1)=tempa(I,KM) IPDCOFCO.217
tempb(I,KMP1)=tempb(I,KM) IPDCOFCO.218
778 CONTINUE IPDCOFCO.219
DO 520 K=2,KM,2 IPDCOFCO.220
DO 520 I=1,IMT IPDCOFCO.221
rrzp(I,K )=tempa(I,K-1)-tempa(I,K ) IPDCOFCO.222
rrzp(I,K+1)=tempb(I,K )-tempb(I,K+1) IPDCOFCO.223
520 CONTINUE IPDCOFCO.224
DO 530 K=2,KM IPDCOFCO.225
DO 530 I=1,IMT IPDCOFCO.226
rrzp(I,K)=FM(I,K)*rrzp(I,K)*fxe IPDCOFCO.227
530 CONTINUE IPDCOFCO.228
DO 540 I=1,IMT IPDCOFCO.229
rrzp(I,KMP1)=0. IPDCOFCO.230
540 CONTINUE IPDCOFCO.231
C IPDCOFCO.232
C Set e and esav to zero initially IPDCOFCO.233
C IPDCOFCO.234
DO 542 M=1,3 IPDCOFCO.235
DO 542 K=1,KMP1 IPDCOFCO.236
DO 542 I=1,IMT IPDCOFCO.237
e(I,K,M)=0. IPDCOFCO.238
542 CONTINUE IPDCOFCO.239
DO 544 M=1,NT IPDCOFCO.240
DO 544 K=1,KM IPDCOFCO.241
DO 544 I=1,IMT IPDCOFCO.242
esav(I,K,M)=0. IPDCOFCO.243
544 CONTINUE IPDCOFCO.244
ENDIF IPDCOFCO.245
C IPDCOFCO.246
C -------------------------------------------------------------- IPDCOFCO.247
CL Calculation of delta-rho in all three model directions for IPDCOFCO.248
CL all points in the slab IPDCOFCO.249
C -------------------------------------------------------------- IPDCOFCO.250
C IPDCOFCO.251
DO 546 K=1,KMP1 IPDCOFCO.252
DO 546 I=1,IMT IPDCOFCO.253
rrz(I,K)=rrzp(I,K) IPDCOFCO.254
546 CONTINUE IPDCOFCO.255
DO 548 K=1,KM IPDCOFCO.256
DO 548 I=1,IMT IPDCOFCO.257
rx(I,K)=rxp(I,K) IPDCOFCO.258
rym(I,K)=ry(I,K) IPDCOFCO.259
548 CONTINUE IPDCOFCO.260
IPDCOFCO.261
IPDCOFCO.262
CALL STATEC(TP,TP(1,1,2),tempa,TDIF,TDIF(1,1,2),1,IMT,KM,J+1 IPDCOFCO.263
&,JMT ORH7F404.58
&) IPDCOFCO.266
CALL STATEC(TP,TP(1,1,2),tempb,TDIF,TDIF(1,1,2),2,IMT,KM,J+1 IPDCOFCO.267
&,JMT ORH7F404.59
&) IPDCOFCO.270
IPDCOFCO.271
IPDCOFCO.272
DO 779 I=1,IMT IPDCOFCO.273
tempa(I,KMP1)=tempa(I,KM) IPDCOFCO.274
tempb(I,KMP1)=tempb(I,KM) IPDCOFCO.275
779 CONTINUE IPDCOFCO.276
DO 550 K=2,KM,2 IPDCOFCO.277
DO 550 I=1,IMT IPDCOFCO.278
rrzp(I,K )=tempa(I,K-1)-tempa(I,K ) IPDCOFCO.279
rrzp(I,K+1)=tempb(I,K )-tempb(I,K+1) IPDCOFCO.280
550 CONTINUE IPDCOFCO.281
DO 560 K=2,KM IPDCOFCO.282
DO 560 I=1,IMT IPDCOFCO.283
rrzp(I,K)=FMP(I,K)*rrzp(I,K)*fxe IPDCOFCO.284
560 CONTINUE IPDCOFCO.285
DO 570 I=1,IMT IPDCOFCO.286
rrzp(I,KMP1)=0. IPDCOFCO.287
570 CONTINUE IPDCOFCO.288
DO 580 K=1,KM IPDCOFCO.289
DO 575 I=1,IMTM1 IPDCOFCO.290
rxp(I,K)=FMP(I,K)*FMP(I+1,K)*(RHON(I+1,K)-RHON(I,K))*fxe IPDCOFCO.291
575 CONTINUE IPDCOFCO.292
IF (L_OCYCLIC) THEN IPDCOFCO.293
rxp(IMT,K)=rxp(2,K) IPDCOFCO.294
ELSE IPDCOFCO.295
rxp(IMT,K)=0.0 IPDCOFCO.296
ENDIF IPDCOFCO.297
DO 576 I=1,IMT IPDCOFCO.298
ry (I,K)=FMP(I,K)*FM (I ,K)*(RHON(I ,K)-RHOS(I,K))*fxe IPDCOFCO.299
576 CONTINUE IPDCOFCO.300
580 CONTINUE IPDCOFCO.301
C IPDCOFCO.302
C -------------------------------------------------------------- IPDCOFCO.303
CL In this section, (d-rho/d-x , d-rho/d-y , d-rho/d-z) (contained IPDCOFCO.304
CL in e) and the 1st row of the diffusion matrix (contained in fk1) IPDCOFCO.305
CL are calculated. IPDCOFCO.306
CL A numerical stability check on d-rho/d-z, Equation (1.23), is IPDCOFCO.307
CL also performed. IPDCOFCO.308
C -------------------------------------------------------------- IPDCOFCO.309
C IPDCOFCO.310
DO 590 K=1,KM IPDCOFCO.311
DO 585 I=1,IMTM1 IPDCOFCO.312
e(I,K,1)=(fxb*CSTR(J))*rx(I,K)*DXU2RQ(I,K) IPDCOFCO.313
e(I,K,2)=(ry(I+1,K)+ry(I,K)+rym(I+1,K)+rym(I,K))*DYT4R(J) IPDCOFCO.314
e(I,K,3)=(rrz(I,K)+rrz(I+1,K)+rrz(I,K+1)+rrz(I+1,K+1)) IPDCOFCO.315
& *fxc*DZ2RQ(I,K) IPDCOFCO.316
tempa(I,K)=-SQRT(e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2)) IPDCOFCO.317
& *(DZ2RQ(I,K)*slmxr(K)) IPDCOFCO.318
585 CONTINUE IPDCOFCO.319
DO I=1,3 IPDCOFCO.320
IF (L_OCYCLIC) THEN IPDCOFCO.321
e(IMT,K,I)=e(2,K,I) IPDCOFCO.322
ELSE IPDCOFCO.323
e(IMT,K,I)=0.0 IPDCOFCO.324
ENDIF IPDCOFCO.325
! IPDCOFCO.326
ENDDO IPDCOFCO.327
! IPDCOFCO.328
IF (L_OCYCLIC) THEN IPDCOFCO.329
tempa(IMT,K)=tempa(2,K) IPDCOFCO.330
ELSE IPDCOFCO.331
tempa(IMT,K)=0.0 IPDCOFCO.332
ENDIF IPDCOFCO.333
590 CONTINUE IPDCOFCO.334
C IPDCOFCO.335
C Performance of stability check IPDCOFCO.336
C IPDCOFCO.337
DO 600 K=1,KM IPDCOFCO.338
DO 600 I=1,IMT IPDCOFCO.339
IF(e(I,K,3).GT.tempa(I,K)) e(I,K,3)=tempa(I,K) IPDCOFCO.340
600 CONTINUE IPDCOFCO.341
C IPDCOFCO.342
630 CONTINUE IPDCOFCO.343
IF (L_OMIXLAY) THEN IPDCOFCO.344
C IPDCOFCO.345
C Set e(I,K,1) to zero in Mixed Layer IPDCOFCO.346
C IPDCOFCO.347
DO 635 K=1,KM IPDCOFCO.348
DO 635 I=1,IMT IPDCOFCO.349
IF ((mld(I)*100.) .GE. ZDZ(K)) e(I,K,1)=0. IPDCOFCO.350
635 CONTINUE IPDCOFCO.351
C IPDCOFCO.352
ENDIF IPDCOFCO.353
! IPDCOFCO.354
DO 640 K=1,KM IPDCOFCO.355
DO 640 I=1,IMT IPDCOFCO.356
tempb(I,K)=-e(I,K,1)/(e(I,K,3)*e(I,K,3)+fxa) IPDCOFCO.357
fk1(I,K,1)=fxd IPDCOFCO.358
fk1(I,K,2)=e(I,K,2)*tempb(I,K) IPDCOFCO.359
fk1(I,K,3)=e(I,K,3)*tempb(I,K) IPDCOFCO.360
640 CONTINUE IPDCOFCO.361
C IPDCOFCO.362
C -------------------------------------------------------------- IPDCOFCO.363
CL Previous section repeated for 2nd and 3rd rows IPDCOFCO.364
CL of diffusion matrix (contained in fk2 and fk3). IPDCOFCO.365
C -------------------------------------------------------------- IPDCOFCO.366
C IPDCOFCO.367
DO 650 K=1,KM IPDCOFCO.368
DO 645 I=2,IMT IPDCOFCO.369
e(I,K,1)=(rxp(I,K)+rxp(I-1,K)+rx(I,K)+rx(I-1,K)) IPDCOFCO.370
& *CSR(J)*DXT4RQ(I,K) IPDCOFCO.371
e(I,K,2)=ry(I,K)*DYUR(J) IPDCOFCO.372
e(I,K,3)=(rrzp(I,K)+rrz(I,K)+rrzp(I,K+1)+rrz(I,K+1)) IPDCOFCO.373
& *fxc*DZ2RQ(I,K) IPDCOFCO.374
tempa(I,K)=-SQRT(e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2)) IPDCOFCO.375
& *(DZ2RQ(I,K)*slmxr(K)) IPDCOFCO.376
645 CONTINUE IPDCOFCO.377
DO I=1,3 IPDCOFCO.378
IF (L_OCYCLIC) THEN IPDCOFCO.379
e(1,K,I)=e(IMTM1,K,I) IPDCOFCO.380
ELSE IPDCOFCO.381
e(1,K,I)=0.0 IPDCOFCO.382
ENDIF IPDCOFCO.383
! IPDCOFCO.384
ENDDO IPDCOFCO.385
! IPDCOFCO.386
IF (L_OCYCLIC) THEN IPDCOFCO.387
tempa(1,K)=tempa(IMTM1,K) IPDCOFCO.388
ELSE IPDCOFCO.389
tempa(1,K)=0.0 IPDCOFCO.390
ENDIF IPDCOFCO.391
650 CONTINUE IPDCOFCO.392
DO 660 K=1,KM IPDCOFCO.393
DO 660 I=1,IMT IPDCOFCO.394
IF(e(I,K,3).GT.tempa(I,K)) e(I,K,3)=tempa(I,K) IPDCOFCO.395
660 CONTINUE IPDCOFCO.396
690 CONTINUE IPDCOFCO.397
IF (L_OMIXLAY) THEN IPDCOFCO.398
C IPDCOFCO.399
C Set e(I,K,2) to zero in Mixed Layer IPDCOFCO.400
C IPDCOFCO.401
DO 695 K=1,KM IPDCOFCO.402
DO 695 I=1,IMT IPDCOFCO.403
IF ((mld(I)*100.) .GT. ZDZ(K)) e(I,K,2)=0. IPDCOFCO.404
695 CONTINUE IPDCOFCO.405
C IPDCOFCO.406
ENDIF IPDCOFCO.407
DO 700 K=1,KM IPDCOFCO.408
DO 700 I=1,IMT IPDCOFCO.409
tempb(I,K)=-e(I,K,2)/(e(I,K,3)*e(I,K,3)+fxa) IPDCOFCO.410
fk2(I,K,1)=e(I,K,1)*tempb(I,K) IPDCOFCO.411
fk2(I,K,2)=fxd IPDCOFCO.412
fk2(I,K,3)=e(I,K,3)*tempb(I,K) IPDCOFCO.413
700 CONTINUE IPDCOFCO.414
C IPDCOFCO.415
DO 708 K=2,KM IPDCOFCO.416
DO 705 I=2,IMT IPDCOFCO.417
e(I,K,1)=(rx(I,K-1)+rx(I-1,K-1)+rx(I,K)+rx(I-1,K)) IPDCOFCO.418
& *CSTR(J)*DXT4RQ(I,K) IPDCOFCO.419
e(I,K,2)=(ry(I,K-1)+rym(I,K-1)+ry(I,K)+rym(I,K))*DYT4R(J) IPDCOFCO.420
e(I,K,3)=rrz(I,K)*DZZ2RQ(I,K)*fxb IPDCOFCO.421
tempb(I,K)=e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2) IPDCOFCO.422
tempa(I,K)=-SQRT(tempb(I,K))*(DZZ2RQ(I,K)*slmxr(K)) IPDCOFCO.423
705 CONTINUE IPDCOFCO.424
IF (L_OCYCLIC) THEN IPDCOFCO.425
DO I = 1,3 IPDCOFCO.426
e(1,K,I)=e(IMTM1,K,I) IPDCOFCO.427
ENDDO IPDCOFCO.428
tempa(1,K)=tempa(IMTM1,K) IPDCOFCO.429
tempb(1,K)=tempb(IMTM1,K) IPDCOFCO.430
ELSE IPDCOFCO.431
DO I = 1,3 IPDCOFCO.432
e(1,K,I)=0.0 IPDCOFCO.433
ENDDO IPDCOFCO.434
tempa(1,K)=0.0 IPDCOFCO.435
tempb(1,K)=0.0 IPDCOFCO.436
ENDIF IPDCOFCO.437
! IPDCOFCO.438
708 CONTINUE IPDCOFCO.439
K=1 IPDCOFCO.440
DO 710 I=2,IMT IPDCOFCO.441
e(I,K,1)=(rx(I,K)+rx(I-1,K))*2.*CSTR(J)*DXT4RQ(I,K) IPDCOFCO.442
e(I,K,2)=(ry(I,K)+rym(I,K))*2.*DYT4R(J) IPDCOFCO.443
e(I,K,3)=rrz(I,K)*DZZ2RQ(I,K)*fxb IPDCOFCO.444
tempb(I,K)=e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2) IPDCOFCO.445
tempa(I,K)=-SQRT(tempb(I,K))*(DZZ2RQ(I,K)*slmxr(K)) IPDCOFCO.446
710 CONTINUE IPDCOFCO.447
IF (L_OCYCLIC) THEN IPDCOFCO.448
DO I = 1,3 IPDCOFCO.449
e(1,K,I)=e(IMTM1,K,I) IPDCOFCO.450
ENDDO IPDCOFCO.451
tempa(1,K)=tempa(IMTM1,K) IPDCOFCO.452
tempb(1,K)=tempb(IMTM1,K) IPDCOFCO.453
ELSE IPDCOFCO.454
DO I = 1,3 IPDCOFCO.455
e(1,K,I)=0.0 IPDCOFCO.456
ENDDO IPDCOFCO.457
tempa(1,K)=0.0 IPDCOFCO.458
tempb(1,K)=0.0 IPDCOFCO.459
ENDIF IPDCOFCO.460
C IPDCOFCO.461
DO 720 K=1,KM IPDCOFCO.462
DO 720 I=1,IMT IPDCOFCO.463
IF(e(I,K,3).GT.tempa(I,K)) e(I,K,3)=tempa(I,K) IPDCOFCO.464
720 CONTINUE IPDCOFCO.465
750 CONTINUE IPDCOFCO.466
IF (L_OMIXLAY) THEN IPDCOFCO.467
C IPDCOFCO.468
C Set e(I,K,3) to zero in Mixed Layer IPDCOFCO.469
C IPDCOFCO.470
DO 755 K=1,KM IPDCOFCO.471
DO 755 I=1,IMT IPDCOFCO.472
IF ((mld(I)*100.) .GE. ZDZ(K)) THEN IPDCOFCO.473
e(I,K,3)=0. IPDCOFCO.474
tempb(I,K)=0. IPDCOFCO.475
ENDIF IPDCOFCO.476
755 CONTINUE IPDCOFCO.477
C IPDCOFCO.478
ENDIF IPDCOFCO.479
DO 760 K=1,KM IPDCOFCO.480
DO 760 I=1,IMT IPDCOFCO.481
tempa(I,K)=fxd/(e(I,K,3)*e(I,K,3)+fxa) IPDCOFCO.482
C IPDCOFCO.483
C At this stage, fk3(I,K,3) only contains the delta-squared term; IPDCOFCO.484
C the epsilon term is added below. IPDCOFCO.485
C IPDCOFCO.486
fk3(I,K,3)=tempa(I,K)*tempb(I,K) IPDCOFCO.487
tempb(I,K)=-e(I,K,3)*tempa(I,K) IPDCOFCO.488
fk3(I,K,1)=e(I,K,1)*tempb(I,K) IPDCOFCO.489
fk3(I,K,2)=e(I,K,2)*tempb(I,K) IPDCOFCO.490
760 CONTINUE IPDCOFCO.491
C IPDCOFCO.492
C -------------------------------------------------------------- IPDCOFCO.493
CL Creates the diffusion matrix in its entirety IPDCOFCO.494
C -------------------------------------------------------------- IPDCOFCO.495
C IPDCOFCO.496
DO 770 L=1,3 IPDCOFCO.497
DO 770 K=1,KM IPDCOFCO.498
DO 770 I=1,IMT IPDCOFCO.499
fk1(I,K,L)=fk1(I,K,L)*ahi(K) IPDCOFCO.500
fk2(I,K,L)=fk2(I,K,L)*ahi(K) IPDCOFCO.501
fk3(I,K,L)=fk3(I,K,L)*ahi(K) IPDCOFCO.502
770 CONTINUE IPDCOFCO.503
C IPDCOFCO.504
C -------------------------------------------------------------- IPDCOFCO.505
CL To prevent horizontal roughnesses in the tracer fields IPDCOFCO.506
CL small hor. and vert. diffusion coeffs. are added to the IPDCOFCO.507
CL relevant diagonal elements of the diffusion matrix. IPDCOFCO.508
CL The epsilon term is also added in to fk3(I,K,3). IPDCOFCO.509
C -------------------------------------------------------------- IPDCOFCO.510
C IPDCOFCO.511
DO 780 K=1,KM IPDCOFCO.512
DO 780 I=1,IMT IPDCOFCO.513
fk1(I,K,1)=fk1(I,K,1)+AH IPDCOFCO.514
fk2(I,K,2)=fk2(I,K,2)+AH IPDCOFCO.515
fk3(I,K,3)=fk3(I,K,3)+gnu(I,K) IPDCOFCO.516
780 CONTINUE IPDCOFCO.517
C IPDCOFCO.518
DO K=1,KM-1 IPDCOFCO.519
DO I=1,IMT IPDCOFCO.520
FK3(I,K,3)=FK3(I,K,3)*FM(I,K) IPDCOFCO.521
END DO IPDCOFCO.522
END DO IPDCOFCO.523
C IPDCOFCO.524
790 CONTINUE IPDCOFCO.525
RETURN IPDCOFCO.526
END IPDCOFCO.527
*ENDIF IPDCOFCO.528