SUBROUTINE IPDCOFCL 4,5IPDCOFCL.41
& ( J,JMT,IMT,IMTM1,KM,KMP,KMP1,KMP2,NT,NTMIN2, OLA0F401.367
& J_OFFSET, ORH3F402.436
& T,TP,TDIF, IPDCOFCL.43
& DZZ,KMT,KMTP,KMTPP, OLA0F401.368
& DXUR,DXU2RQ,DXT4RQ,DYUR,DYT4R,DZ2RQ,DZZ2RQ,ZDZ,DTTS, IPDCOFCL.44
& NERGY,CSR,CSTR,ITT,FM,FMP,FMM, OLA0F401.369
& RHOS,RHON,AH,ahi, IPDCOFCL.46
& gnu,JFT0,esav,fk1,fk2,fk3, IPDCOFCL.47
& rxp,ry,rrzp IPDCOFCL.48
& ,drhob1p,drhob2p OLA0F401.370
& ,mld IPDCOFCL.51
& ,slope_max IPDCOFCL.54
&,dslope,slopec OLA0F401.371
& ) IPDCOFCL.55
C IPDCOFCL.56
C IPDCOFCL.57
IMPLICIT NONE IPDCOFCL.58
C IPDCOFCL.59
*CALL OARRYSIZ 
ORH6F401.18
INTEGER IPDCOFCL.60
& I,J,K,M,N,L, IPDCOFCL.61
& JMT,IMT,IMTM1,KM,KMP,KMP1,KMP2,NT,NTMIN2, OLA0F401.372
& NERGY,ITT,KS,JFT0 IPDCOFCL.63
& ,KMT(IMT),KMTP(IMT),KMTPP(IMT) OLA0F401.373
&,J_OFFSET ! IN Offset used in working out global equivalent ORH3F402.437
! values of row number. (0 for non mpp version) ORH3F402.438
INTEGER kminj(imt) OLA0F401.374
C IPDCOFCL.64
REAL IPDCOFCL.65
& T(IMT,KM,NT),TP(IMT,KM,NT),TDIF(IMT,KMP2,NTMIN2), IPDCOFCL.66
& DZZ(KMP1), OLA0F401.375
& DXUR(IMT),DXU2RQ(IMT,KM),DXT4RQ(IMT,KM), IPDCOFCL.67
& DYUR(JMT),DYT4R(JMT),DZ2RQ(IMT,KM),DZZ2RQ(IMT,KM),ZDZ(KM),DTTS, IPDCOFCL.68
& RHOS(IMT,KM),RHON(IMT,KM), IPDCOFCL.69
& FM(IMT,KM),FMP(IMT,KM),FMM(IMT,KM),CSR(JMT),CSTR(JMT), OLA0F401.376
& AH, ! IN BACKGROUND horizontal diffusivity IPDCOFCL.71
C Ah in Eq.1.10 (cm2/s) IPDCOFCL.72
& ahi(KM), ! IN along-isopycnal coeff. of diffusion IPDCOFCL.73
C Ai in Eq.1.8 (cm2/s) IPDCOFCL.74
& gnu(IMT,KM), ! IN coeff. of vertical diffusion at TOP of IPDCOFCL.75
C box. Av in Eq.1.10 (cm2/s) IPDCOFCL.76
& fk1(IMT,KM,3), ! OUT } Rows 1,2,3 of diffusion matrix (Eq1.8) IPDCOFCL.77
& fk2(IMT,KM,3), ! OUT } fk1 is on EAST face, fk2 on NORTH face IPDCOFCL.78
& fk3(IMT,KM,3), ! OUT } & fk3 on TOP face of grid box (i,j,k) IPDCOFCL.79
& rxp(IMT,KM), ! INOUT delta-rho in x dirn, row J+1 (E face) IPDCOFCL.80
& ry(IMT,KM), ! INOUT delta-rho in y dirn, row J (N face) IPDCOFCL.81
& rrzp(IMT,KMP1), ! INOUT delta-rho in z dirn, row J+1 (top face IPDCOFCL.82
C If L_EXTRAP is true, then rrz for (*,1) at level of tracer points OLA0F401.384
& esav(IMT,KM,NT) ! OUT used to save e(I,K,2) in IPDCOFCL.83
C subroutine IPDFLXCL IPDCOFCL.84
& ,mld(IMT_IPD_MIX) ! IN mixed layer depth (m) ORH1F305.4767
& ,slope_max ! IN Maximum slope for diffusion IPDCOFCL.90
REAL OLA0F401.377
& drhob1p(IMT) ! extrapolated density gradient at the bottom of OLA0F401.378
! column i+1/2, row j+1 (relative to T grid) OLA0F401.379
REAL OLA0F401.380
& drhob2p(IMT,2) ! normalised density for row j+1, OLA0F401.381
C (*,1) at level one less than min(kmtp(i),kmtpp(i)) OLA0F401.382
C (*,2) at level min(kmtp(i),kmtpp(i)) OLA0F401.383
C IPDCOFCL.91
C IPDCOFCL.92
C Declare local constants and arrays IPDCOFCL.93
C IPDCOFCL.94
REAL IPDCOFCL.95
& rx(IMT,KM), ! delta-rho in x-dir, row J (E face) IPDCOFCL.96
& rym(IMT,KM), ! " " y-dir, " J-1 (N face) IPDCOFCL.97
& rrz(IMT,KMP1), ! " " z-dir, " J (top face) IPDCOFCL.98
C If L_EXTRAP is true, then rrz for (*,1) at level of tracer points OLA0F401.385
& e(IMT,KMP1,3), ! d-rho/d-x , d-rho/d-y , d-rho/d-z IPDCOFCL.99
& slmxr(KM), ! cst. used in stability check on IPDCOFCL.100
C d-rho-dz in Eq.1.24 IPDCOFCL.101
& tempa(IMT,KMP1), ! workspace IPDCOFCL.102
& tempb(IMT,KMP1), ! workspace IPDCOFCL.103
& fxa, ! } IPDCOFCL.104
& fxb, ! } IPDCOFCL.105
& fxc, ! } local csts. IPDCOFCL.106
& fxd, ! } IPDCOFCL.107
& fxe, ! } IPDCOFCL.108
& temp ! temporary variable IPDCOFCL.109
& ,chkslp(imt,km) ! local variable ORW1F400.4
*CALL CNTLOCN ORH1F305.4761
REAL slope,slopec,dslope, OLA0F401.386
& rloc1,rloc2,rloc3,rloc4 OLA0F401.387
REAL dciso1(imt,km),dciso2(imt,km) OLA0F401.388
REAL OLA0F401.389
& drhob1(imt) ! extrapolated density gradient at the bottom of OLA0F401.390
! column i+1/2, row j (relative to T grid) OLA0F401.391
REAL OLA0F401.392
& drhob2(imt,3) ! normalised density for row j, OLA0F401.393
C (*,1) at level one less than min(kmtp(i),kmtpp(i)) OLA0F401.394
C (*,2) at level min(kmtp(i),kmtpp(i)) OLA0F401.395
C (*,3) is extrapolated density gradient at bottom OLA0F401.396
C of column i,row j+1/2 OLA0F401.397
c OLA0F401.398
c REMARK: throughout the computations, "fk3(.,1,.)" computed at the OLA0F401.399
c center of the top face of the first vertical level (i.e. OLA0F401.400
c at the ocean surface) will stay as zero, representing no OLA0F401.401
c diffusion of tracers through the ocean surface. OLA0F401.402
c OLA0F401.403
c OLA0F401.404
C* IPDCOFCL.110
C*L---- External subroutines called ------------------------------------ IPDCOFCL.111
C STATEC - calculates density IPDCOFCL.112
C MATRIX - prints out data IPDCOFCL.113
C* IPDCOFCL.114
C IPDCOFCL.115
C -------------------------------------------------------------- IPDCOFCL.116
CL Initial conditions for the isopycnal diffusion code. IPDCOFCL.117
CL These involve working out delta-rho in the x and z directions IPDCOFCL.118
CL for the second row, so that the process for calculating IPDCOFCL.119
CL these quantities for the other rows can be initiated. IPDCOFCL.120
C -------------------------------------------------------------- IPDCOFCL.121
*IF DEF,OISOPYC ORH1F305.457
C IPDCOFCL.122
DO 502 N=1,3 IPDCOFCL.123
DO 502 K=1,KM IPDCOFCL.124
DO 502 I=1,IMT IPDCOFCL.125
fk1(I,K,N)=0. IPDCOFCL.126
fk2(I,K,N)=0. IPDCOFCL.127
fk3(I,K,N)=0. IPDCOFCL.128
502 CONTINUE IPDCOFCL.129
IF (L_OISOTAPER) THEN OLA0F401.405
do k=1,km OLA0F401.406
do i=1,imt OLA0F401.407
dciso1(i,k) = 0. OLA0F401.408
dciso2(i,k) = 0. OLA0F401.409
enddo OLA0F401.410
enddo OLA0F401.411
ENDIF OLA0F401.412
OLA0F401.413
C need local var kminj that contains the value of min(kmt(i),kmtp(i)) OLA0F401.414
IF (L_OEXTRAP) THEN OLA0F401.415
do i=1,imt OLA0F401.416
kminj(i)=min(kmt(i),kmtp(i)) OLA0F401.417
drhob1(i)=0. OLA0F401.418
drhob2(i,1)=0. OLA0F401.419
drhob2(i,2)=0. OLA0F401.420
drhob2(i,3)=0. OLA0F401.421
enddo OLA0F401.422
ENDIF OLA0F401.423
OLA0F401.424
fxa=1.E-25 IPDCOFCL.130
fxb=2. IPDCOFCL.131
fxc=.5 IPDCOFCL.132
fxd=1. IPDCOFCL.133
fxe=1.E10 ! Big number used to multiply all density gradients. IPDCOFCL.134
C Always cancels out in the end, so must be IPDCOFCL.135
C something to do with reducing rounding error. IPDCOFCL.136
C IPDCOFCL.137
DO 503 K=1,KM IPDCOFCL.138
C IPDCOFCL.139
C Set maximum allowable slope for the diffusion to slope_max. IPDCOFCL.140
C Also included here (commented out) are two old versions. IPDCOFCL.141
C The first value of temp implements the criterion in Cox IPDCOFCL.142
C (Ocean modelling #74, 1987). I think the 4. should be a IPDCOFCL.143
C 1. (RAW 6/5/92), but in any case this criterion gave IPDCOFCL.144
C poor results in a 2.5x3.75 GCM. The second, hardwired value IPDCOFCL.145
C of temp has worked in that GCM and in the 1x1 FOAM model. IPDCOFCL.146
C In either of these cases note that slmxr contains 2*temp, IPDCOFCL.147
C but that this is cancelled by a division by 2*dzz later on. IPDCOFCL.148
C IPDCOFCL.149
C In all cases, the maximum slope through which the diffusion IPDCOFCL.150
C tensor can be rotated is set to 2*DZ(k)/slmxr(k). IPDCOFCL.151
C ORW1F400.5
C In old (< vn 4.0) versions of the code, if the slope ORW1F400.6
C exceeded slope_max then e(i,k,3) was reset (e.g. in the ORW1F400.7
C DO 600 loop etc) to give a diffusivity with a slope of ORW1F400.8
C slope_max. This implied a diapycnal component to the ORW1F400.9
C diffusive flux. For vn>=4.0, the sheme of Gerdes, Koeberle ORW1F400.10
C and Willebrand (J. Climate 5:211-226, 1991) is used. This ORW1F400.11
C scheme scales down the isopycnal diffusivity if slope > ORW1F400.12
C slope_max, while preserving the isopycnal orientation of ORW1F400.13
C diffusive flux. ORW1F400.14
C IPDCOFCL.152
C temp=4.*DTTS*ahi(K) ! Cox's criterion IPDCOFCL.153
C temp=2.e13 ! Oscar's FOAM value IPDCOFCL.154
C slmxr(K)=(temp+temp)* IPDCOFCL.155
! For L_OFILTER = true -> MAX(DYUR(J),DXUR(1)*MIN(CSTR(J),CSTR(JFT0))) ORH1F305.4768
! For L_OFILTER = false-> MAX(DYUR(J),DXUR(1)*CSTR(J)) ORH1F305.4769
! ORH1F305.4770
C IPDCOFCL.164
slmxr(K)=1./slope_max/DZ2RQ(2,K) ! DZ2RQ(2,.) chosen arbitraril IPDCOFCL.165
503 CONTINUE IPDCOFCL.166
C IPDCOFCL.167
IF(J+J_OFFSET.EQ.2) THEN ORH3F402.439
DO 505 K=1,KM IPDCOFCL.169
DO 505 I=1,IMT IPDCOFCL.170
ry(I,K)=0. IPDCOFCL.171
505 CONTINUE IPDCOFCL.172
DO 508 K=1,KMP1 IPDCOFCL.173
DO 508 I=1,IMT IPDCOFCL.174
rrzp(I,K)=0. IPDCOFCL.175
508 CONTINUE IPDCOFCL.176
DO 510 K=1,KM IPDCOFCL.177
DO 509 I=1,IMTM1 IPDCOFCL.178
rxp(I,K)=FM(I,K)*FM(I+1,K)*(RHOS(I+1,K)-RHOS(I,K))*fxe IPDCOFCL.179
509 CONTINUE IPDCOFCL.180
IF (L_OCYCLIC) THEN ORH1F305.4771
rxp(IMT,K)=rxp(2,K) ORH1F305.4772
ELSE ORH1F305.4773
rxp(IMT,K)=0.0 ORH1F305.4774
ENDIF ORH1F305.4775
510 CONTINUE IPDCOFCL.189
DO 777 K=1,KMP1 IPDCOFCL.190
DO 777 I=1,IMT IPDCOFCL.191
tempa(I,K)=0. ! Initialise tempa IPDCOFCL.192
tempb(I,K)=0. ! Initialise tempb IPDCOFCL.193
777 CONTINUE IPDCOFCL.194
JA121293.301
JA121293.303
CALL STATEC(T,T(1,1,2),tempa,TDIF,TDIF(1,1,2),1,IMT,KM,2 JA121293.304
&,JMT ORH7F404.61
&) JA121293.307
CALL STATEC(T,T(1,1,2),tempb,TDIF,TDIF(1,1,2),2,IMT,KM,2 JA121293.308
&,JMT ORH7F404.62
&) JA121293.311
JA121293.312
JA121293.317
DO 778 I=1,IMT IPDCOFCL.197
tempa(I,KMP1)=tempa(I,KM) IPDCOFCL.198
tempb(I,KMP1)=tempb(I,KM) IPDCOFCL.199
778 CONTINUE IPDCOFCL.200
DO 520 K=2,KM,2 IPDCOFCL.201
DO 520 I=1,IMT IPDCOFCL.202
rrzp(I,K )=tempa(I,K-1)-tempa(I,K ) IPDCOFCL.203
rrzp(I,K+1)=tempb(I,K )-tempb(I,K+1) IPDCOFCL.204
520 CONTINUE IPDCOFCL.205
IF (L_OEXTRAP) THEN OLA0F401.425
CALL EXTRAP OLA0F401.426
& (imt,imtm1,kmp1,km,fxe, OLA0F401.427
& kmt,kmtp,dzz,dzz2rq,dz2rq, OLA0F401.428
& tempa,tempb,rrzp,drhob1p,drhob2p ORH3F403.41
& ) OLA0F401.430
ENDIF OLA0F401.431
OLA0F401.432
IF (.NOT. L_OEXTRAP) THEN OLA0F401.433
C k loop from 2 to KM when setting rrz(k=1) to zero OLA0F401.434
DO 530 K=2,KM IPDCOFCL.206
DO 530 I=1,IMT IPDCOFCL.207
rrzp(I,K)=FM(I,K)*rrzp(I,K)*fxe IPDCOFCL.208
530 CONTINUE IPDCOFCL.209
ENDIF OLA0F401.435
IF (L_OEXTRAP) THEN OLA0F401.436
C If extrapolation required, value of rrzp at k=1 is redefined to OLA0F401.437
C be the extrapolated density gradient at the level of the tracer OLA0F401.438
C points (rather than at the top of the tracer gridbox as the OLA0F401.439
C other values (k=2,km)) OLA0F401.440
DO K=1,KM OLA0F401.441
DO I=1,IMT OLA0F401.442
rrzp(I,K)=FM(I,K)*rrzp(I,K)*fxe OLA0F401.443
ENDDO OLA0F401.444
ENDDO OLA0F401.445
ENDIF OLA0F401.446
DO 540 I=1,IMT IPDCOFCL.210
rrzp(I,KMP1)=0. IPDCOFCL.211
540 CONTINUE IPDCOFCL.212
C IPDCOFCL.213
C Set e and esav to zero initially IPDCOFCL.214
C IPDCOFCL.215
DO 542 M=1,3 IPDCOFCL.216
DO 542 K=1,KMP1 IPDCOFCL.217
DO 542 I=1,IMT IPDCOFCL.218
e(I,K,M)=0. IPDCOFCL.219
542 CONTINUE IPDCOFCL.220
DO 544 M=1,NT IPDCOFCL.221
DO 544 K=1,KM IPDCOFCL.222
DO 544 I=1,IMT IPDCOFCL.223
esav(I,K,M)=0. IPDCOFCL.224
544 CONTINUE IPDCOFCL.225
ENDIF IPDCOFCL.226
C IPDCOFCL.227
C -------------------------------------------------------------- IPDCOFCL.228
CL Calculation of delta-rho in all three model directions for IPDCOFCL.229
CL all points in the slab IPDCOFCL.230
C -------------------------------------------------------------- IPDCOFCL.231
C IPDCOFCL.232
DO 546 K=1,KMP1 IPDCOFCL.233
DO 546 I=1,IMT IPDCOFCL.234
rrz(I,K)=rrzp(I,K) IPDCOFCL.235
546 CONTINUE IPDCOFCL.236
IF (L_OEXTRAP) THEN OLA0F401.447
C Move extrapolated data down a row OLA0F401.448
do i=1,imt OLA0F401.449
drhob1(i)=drhob1p(i) OLA0F401.450
drhob2(i,1)=drhob2p(i,1) OLA0F401.451
drhob2(i,2)=drhob2p(i,2) OLA0F401.452
enddo OLA0F401.453
ENDIF OLA0F401.454
DO 548 K=1,KM IPDCOFCL.237
DO 548 I=1,IMT IPDCOFCL.238
rx(I,K)=rxp(I,K) IPDCOFCL.239
rym(I,K)=ry(I,K) IPDCOFCL.240
548 CONTINUE IPDCOFCL.241
JA121293.318
JA121293.320
CALL STATEC(TP,TP(1,1,2),tempa,TDIF,TDIF(1,1,2),1,IMT,KM,J+1 JA121293.321
&,JMT ORH7F404.63
&) JA121293.324
CALL STATEC(TP,TP(1,1,2),tempb,TDIF,TDIF(1,1,2),2,IMT,KM,J+1 JA121293.325
&,JMT ORH7F404.64
&) JA121293.328
JA121293.329
JA121293.334
DO 779 I=1,IMT IPDCOFCL.244
tempa(I,KMP1)=tempa(I,KM) IPDCOFCL.245
tempb(I,KMP1)=tempb(I,KM) IPDCOFCL.246
779 CONTINUE IPDCOFCL.247
DO 550 K=2,KM,2 IPDCOFCL.248
DO 550 I=1,IMT IPDCOFCL.249
rrzp(I,K )=tempa(I,K-1)-tempa(I,K ) IPDCOFCL.250
rrzp(I,K+1)=tempb(I,K )-tempb(I,K+1) IPDCOFCL.251
550 CONTINUE IPDCOFCL.252
IF (L_OEXTRAP) THEN OLA0F401.455
C extrapolate to find density gradient at k=1 on tracer levels OLA0F401.456
k = 1 OLA0F401.457
rloc1 = dzz(k)+dzz(k+1) OLA0F401.458
do i=1,imt OLA0F401.459
rloc2 = 2.*dzz2rq(i,k+1)*(tempa(i,k)*rloc1-tempa(i,k+1)*dzz(k)) OLA0F401.460
rrzp(i,k) = (rloc2-0.5*(tempa(i,k+1)+tempa(i,k))) OLA0F401.461
enddo OLA0F401.462
do i=1,imtm1 OLA0F401.463
k=min(kmtp(i),kmtp(i+1)) OLA0F401.464
if (k .ne. 0) then OLA0F401.465
rloc1=dzz(k)+dzz(k+1) OLA0F401.466
if (mod(k,2).eq.0) then OLA0F401.467
rloc2=0.5*(tempa(i,k-1)+tempa(i+1,k-1)) OLA0F401.468
rloc3=0.5*(tempa(i,k) +tempa(i+1,k)) OLA0F401.469
else OLA0F401.470
rloc2=0.5*(tempb(i,k-1)+tempb(i+1,k-1)) OLA0F401.471
rloc3=0.5*(tempb(i,k) +tempb(i+1,k)) OLA0F401.472
endif OLA0F401.473
rloc4=2.*dzz2rq(i,k)*(rloc3*rloc1-rloc2*dzz(k+1)) OLA0F401.474
drhob1p(i)=2.*dz2rq(i,k)*(0.5*(rloc2+rloc3)-rloc4)*fxe OLA0F401.475
endif OLA0F401.476
enddo OLA0F401.477
do i=1,imt OLA0F401.478
k=kminj(i) OLA0F401.479
if (k.ne.0) then OLA0F401.480
rloc1=dzz(k)+dzz(k+1) OLA0F401.481
if (mod(k,2).eq.0) then OLA0F401.482
rloc2=0.5*(drhob2(i,1)+tempa(i,k-1)) OLA0F401.483
rloc3=0.5*(drhob2(i,2)+tempa(i,k)) OLA0F401.484
else OLA0F401.485
rloc2=0.5*(drhob2(i,1)+tempb(i,k-1)) OLA0F401.486
rloc3=0.5*(drhob2(i,2)+tempb(i,k)) OLA0F401.487
endif OLA0F401.488
rloc4=2.*dzz2rq(i,k)*(rloc3*rloc1-rloc2*dzz(k+1)) OLA0F401.489
drhob2(i,3)=2.*dz2rq(i,k)*(0.5*(rloc2+rloc3)-rloc4)*fxe OLA0F401.490
endif OLA0F401.491
enddo OLA0F401.492
do i=1,imt OLA0F401.493
k=min(kmtp(i),kmtpp(i)) OLA0F401.494
if (k.ne.0) then OLA0F401.495
if (mod(k,2).eq.0) then OLA0F401.496
drhob2p(i,1)=tempa(i,k-1) OLA0F401.497
drhob2p(i,2)=tempa(i,k) OLA0F401.498
else OLA0F401.499
drhob2p(i,1)=tempb(i,k-1) OLA0F401.500
drhob2p(i,2)=tempb(i,k) OLA0F401.501
endif OLA0F401.502
endif OLA0F401.503
enddo OLA0F401.504
OLA0F401.505
ENDIF OLA0F401.506
OLA0F401.507
IF (.NOT. L_OEXTRAP) THEN OLA0F401.508
C k loop from 2 to KM when setting rrz(k=1) to zero OLA0F401.509
DO 560 K=2,KM IPDCOFCL.253
DO 560 I=1,IMT IPDCOFCL.254
rrzp(I,K)=FMP(I,K)*rrzp(I,K)*fxe IPDCOFCL.255
560 CONTINUE IPDCOFCL.256
ENDIF OLA0F401.510
IF (L_OEXTRAP) THEN OLA0F401.511
C If extrapolation required, value of rrzp at k=1 is redefined to OLA0F401.512
C be the extrapolated density gradient at the level of the tracer OLA0F401.513
C points (rather than at the top of the tracer gridbox as the OLA0F401.514
C other values (k=2,km)) OLA0F401.515
DO K=1,KM OLA0F401.516
DO I=1,IMT OLA0F401.517
rrzp(I,K)=FMP(I,K)*rrzp(I,K)*fxe OLA0F401.518
ENDDO OLA0F401.519
ENDDO OLA0F401.520
ENDIF OLA0F401.521
DO 570 I=1,IMT IPDCOFCL.257
rrzp(I,KMP1)=0. IPDCOFCL.258
570 CONTINUE IPDCOFCL.259
DO 580 K=1,KM IPDCOFCL.260
DO 575 I=1,IMTM1 IPDCOFCL.261
rxp(I,K)=FMP(I,K)*FMP(I+1,K)*(RHON(I+1,K)-RHON(I,K))*fxe IPDCOFCL.262
575 CONTINUE IPDCOFCL.263
IF (L_OCYCLIC) THEN ORH1F305.4776
rxp(IMT,K)=rxp(2,K) ORH1F305.4777
ELSE ORH1F305.4778
rxp(IMT,K)=0.0 ORH1F305.4779
ENDIF ORH1F305.4780
DO 576 I=1,IMT IPDCOFCL.272
ry (I,K)=FMP(I,K)*FM (I ,K)*(RHON(I ,K)-RHOS(I,K))*fxe IPDCOFCL.273
576 CONTINUE IPDCOFCL.274
580 CONTINUE IPDCOFCL.275
C IPDCOFCL.276
C -------------------------------------------------------------- IPDCOFCL.277
CL In this section, (d-rho/d-x , d-rho/d-y , d-rho/d-z) (contained IPDCOFCL.278
CL in e) and the 1st row of the diffusion matrix (contained in fk1) IPDCOFCL.279
CL are calculated. IPDCOFCL.280
CL A numerical stability check on d-rho/d-z, Equation (1.23), is IPDCOFCL.281
CL also performed. IPDCOFCL.282
C -------------------------------------------------------------- IPDCOFCL.283
C IPDCOFCL.284
c OLA0F401.522
c in order to compute "fk1(.,.,3)", (xz) component of the tensor, OLA0F401.523
c evaluate the density gradients. both "fk1" and the density OLA0F401.524
c gradients are obtained at the center of the eastern face of the OLA0F401.525
c "t" grid box as indicated by "X" below. OLA0F401.526
c OLA0F401.527
c lat/lon view depth/lon view OLA0F401.528
c u,v j -> +-------+ +-------+ <- zw(k-1) OLA0F401.529
c | | | | OLA0F401.530
c t,s j -> | i X | i X <- zt(k) OLA0F401.531
c | | | | OLA0F401.532
c u,v j-1 -> +-------+ +-------+ <- zw(k) OLA0F401.533
c OLA0F401.534
c first compute "e(.,.,3)". the top and bottom levels will be OLA0F401.535
c considered later. note that the gradients are multiplied by OLA0F401.536
c fxe=1.e10 to keep the exponents in range. note that no masks OLA0F401.537
c are used in the computation of "e(.,.,3)", because it will be OLA0F401.538
c later multiplied by "e(.,.,1)" which involves masks. OLA0F401.539
c OLA0F401.540
IF (L_OEXTRAP) THEN OLA0F401.541
DO K=2,KM-1 OLA0F401.542
DO I=1,IMTM1 OLA0F401.543
e(I,K,3)=(rrz(I,K)+rrz(I+1,K)+rrz(I,K+1)+rrz(I+1,K+1)) OLA0F401.544
& *fxc*DZ2RQ(I,K) OLA0F401.545
ENDDO OLA0F401.546
ENDDO OLA0F401.547
C if using GM scheme, need the interpolated density at k=1 OLA0F401.548
C hence use the previously calculated rrz(k=1) OLA0F401.549
k=1 OLA0F401.550
do i=1,imtm1 OLA0F401.551
e(i,k,3)=dz2rq(i,k)*(rrz(i,k)+rrz(i+1,k)) OLA0F401.552
enddo OLA0F401.553
OLA0F401.554
C consider bottom level OLA0F401.555
do i=1,imtm1 OLA0F401.556
e(i,km,3)=0. OLA0F401.557
enddo OLA0F401.558
OLA0F401.559
do i=1,imtm1 OLA0F401.560
k=min(kmt(i),kmt(i+1)) OLA0F401.561
if (k.ne.0) then OLA0F401.562
e(i,k,3)=drhob1(i) OLA0F401.563
endif OLA0F401.564
enddo OLA0F401.565
OLA0F401.566
ELSE OLA0F401.567
OLA0F401.568
DO K=1,KM OLA0F401.569
DO I=1,IMTM1 OLA0F401.570
e(I,K,3)=(rrz(I,K)+rrz(I+1,K)+rrz(I,K+1)+rrz(I+1,K+1)) OLA0F401.571
& *fxc*DZ2RQ(I,K) OLA0F401.572
ENDDO OLA0F401.573
ENDDO OLA0F401.574
ENDIF OLA0F401.575
OLA0F401.576
do k=1,km ORW1F400.15
do i=1,imt ORW1F400.16
chkslp(i,k)=0. ORW1F400.17
enddo ORW1F400.18
enddo ORW1F400.19
DO 590 K=1,KM IPDCOFCL.285
DO 585 I=1,IMTM1 IPDCOFCL.286
e(I,K,1)=(fxb*CSTR(J))*rx(I,K)*DXU2RQ(I,K) IPDCOFCL.287
e(I,K,2)=(ry(I+1,K)+ry(I,K)+rym(I+1,K)+rym(I,K))*DYT4R(J) IPDCOFCL.288
tempa(I,K)=-SQRT(e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2)) IPDCOFCL.291
& *(DZ2RQ(I,K)*slmxr(K)) IPDCOFCL.292
585 CONTINUE IPDCOFCL.293
DO I=1,3 IPDCOFCL.294
IF (L_OCYCLIC) THEN ORH1F305.4781
e(IMT,K,I)=e(2,K,I) ORH1F305.4782
ELSE ORH1F305.4783
e(IMT,K,I)=0.0 ORH1F305.4784
ENDIF ORH1F305.4785
! ORH1F305.4786
ENDDO ORH1F305.4787
! ORH1F305.4788
IF (L_OCYCLIC) THEN ORH1F305.4789
tempa(IMT,K)=tempa(2,K) ORH1F305.4790
ELSE ORH1F305.4791
tempa(IMT,K)=0.0 ORH1F305.4792
ENDIF ORH1F305.4793
590 CONTINUE IPDCOFCL.312
C IPDCOFCL.313
C Performance of stability check IPDCOFCL.314
C IPDCOFCL.315
DO 600 K=1,KM IPDCOFCL.316
DO 600 I=1,IMT IPDCOFCL.317
chkslp(i,k) = amax1(abs(e(i,k,3)),abs(tempa(i,k))) OLE1F403.1
600 CONTINUE IPDCOFCL.319
C IPDCOFCL.320
630 CONTINUE IPDCOFCL.321
IF (L_OMIXLAY) THEN ORH1F305.4794
C IPDCOFCL.324
C Set e(I,K,1) to zero in Mixed Layer IPDCOFCL.325
C IPDCOFCL.326
DO 635 K=1,KM IPDCOFCL.327
DO 635 I=1,IMT IPDCOFCL.328
IF ((mld(I)*100.) .GE. ZDZ(K)) e(I,K,1)=0. IPDCOFCL.329
635 CONTINUE IPDCOFCL.330
C IPDCOFCL.331
ENDIF ORH1F305.4795
! ORH1F305.4796
c OLA0F401.577
c compute "fk1" OLA0F401.578
c OLA0F401.579
do k=1,km OLA0F401.580
do i=1,imt OLA0F401.581
IF (L_OISOTAPER) THEN OLA0F401.582
rloc1 = 0. OLA0F401.583
rloc2 = sign(1.,e(i,k,3))/(abs(e(i,k,3))+fxa) OLA0F401.584
slope = rloc2*sqrt(e(i,k,1)**2+e(i,k,2)**2) OLA0F401.585
if ((slope.le.0.0) .and. OLA0F401.586
& (slope.ge.(-1./(slmxr(k)*DZ2RQ(I,K))))) then OLA0F401.587
rloc1 = 0.5*(1.+tanh((slope+slopec)/dslope)) OLA0F401.588
endif OLA0F401.589
dciso1(i,k) = rloc1 OLA0F401.590
fk1(i,k,3) = -rloc2*e(i,k,1)*dciso1(i,k) OLA0F401.591
C dciso scales the isopycnal diffusion coefficient dependent on OLA0F401.592
C the slope of the isopycnals (see GM references from subroutine OLA0F401.593
C IPDGMVEL for details) OLA0F401.594
fk1(i,k,1)= 1.*dciso1(i,k) OLA0F401.595
fk1(i,k,2)= -e(i,k,1)*e(i,k,2)*rloc2*rloc2*dciso1(i,k) OLA0F401.596
ELSE OLA0F401.597
tempb(i,k)=fxd / (chkslp(i,k)*chkslp(i,k) +fxa) OLA0F401.598
fk1(i,k,1)= e(i,k,3)*e(i,k,3)*tempb(i,k) OLA0F401.599
fk1(i,k,2)= -e(i,k,1)*e(i,k,2)*tempb(I,K) OLA0F401.600
fk1(i,k,3)= -e(i,k,1)*e(i,k,3)*tempb(I,K) OLA0F401.601
ENDIF OLA0F401.602
enddo OLA0F401.603
enddo OLA0F401.604
OLA0F401.605
C IPDCOFCL.341
C -------------------------------------------------------------- IPDCOFCL.342
CL Previous section repeated for 2nd and 3rd rows IPDCOFCL.343
CL of diffusion matrix (contained in fk2 and fk3). IPDCOFCL.344
C -------------------------------------------------------------- IPDCOFCL.345
C IPDCOFCL.346
c OLA0F401.606
c in order to compute "fk2(.,.,3)", (yz) component of the tensor, OLA0F401.607
c evaluate the density gradients. both "fk2" and the density OLA0F401.608
c gradients are obtained at the center of the northern face of the OLA0F401.609
c "t" grid box as indicated by "X" below. OLA0F401.610
c OLA0F401.611
c lat/lon view depth/lon view OLA0F401.612
c u,v j -> +---X---+ +-------+ <- zw(k-1) OLA0F401.613
c | | | | OLA0F401.614
c t,s j -> | i | | X | <- zt(k) OLA0F401.615
c | | | | OLA0F401.616
c u,v j-1 -> +-------+ +-------+ <- zw(k) OLA0F401.617
c OLA0F401.618
c first compute "e(.,.,3)". the top and bottom levels will be OLA0F401.619
c considered later. note that the gradients are multiplied by OLA0F401.620
c fxe=1.e10 to keep the exponents in range. note that no masks OLA0F401.621
c are used in the computation of "e(.,.,3)", because it will be OLA0F401.622
c later multiplied by "e(.,.,2)" which involves masks. OLA0F401.623
c OLA0F401.624
IF (L_OEXTRAP) THEN OLA0F401.625
DO K=2,KM OLA0F401.626
DO I=2,IMT OLA0F401.627
e(I,K,3)=(rrzp(I,K)+rrz(I,K)+rrzp(I,K+1)+rrz(I,K+1)) OLA0F401.628
& *fxc*DZ2RQ(I,K) OLA0F401.629
ENDDO OLA0F401.630
ENDDO OLA0F401.631
OLA0F401.632
C consider top level OLA0F401.633
k=1 OLA0F401.634
do i=2,imt OLA0F401.635
e(i,k,3)=dz2rq(i,k)*(rrzp(i,k)+rrz(i,k)) OLA0F401.636
enddo OLA0F401.637
OLA0F401.638
C consider bottom level OLA0F401.639
do i=2,imt OLA0F401.640
e(i,km,3)=0. OLA0F401.641
enddo OLA0F401.642
OLA0F401.643
do i=2,imt OLA0F401.644
k=kminj(i) OLA0F401.645
if (k.ne.0) then OLA0F401.646
e(i,k,3)=drhob2(i,3) OLA0F401.647
endif OLA0F401.648
enddo OLA0F401.649
OLA0F401.650
ELSE OLA0F401.651
DO K=1,KM OLA0F401.652
DO I=2,IMT OLA0F401.653
e(I,K,3)=(rrzp(I,K)+rrz(I,K)+rrzp(I,K+1)+rrz(I,K+1)) OLA0F401.654
& *fxc*DZ2RQ(I,K) OLA0F401.655
ENDDO OLA0F401.656
ENDDO OLA0F401.657
ENDIF OLA0F401.658
OLA0F401.659
DO 650 K=1,KM IPDCOFCL.347
DO 645 I=2,IMT IPDCOFCL.348
e(I,K,1)=(rxp(I,K)+rxp(I-1,K)+rx(I,K)+rx(I-1,K)) IPDCOFCL.349
& *CSR(J)*DXT4RQ(I,K) IPDCOFCL.350
e(I,K,2)=ry(I,K)*DYUR(J) IPDCOFCL.351
tempa(I,K)=-SQRT(e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2)) IPDCOFCL.354
& *(DZ2RQ(I,K)*slmxr(K)) IPDCOFCL.355
645 CONTINUE IPDCOFCL.356
DO I=1,3 IPDCOFCL.357
IF (L_OCYCLIC) THEN ORH1F305.4797
e(1,K,I)=e(IMTM1,K,I) ORH1F305.4798
ELSE ORH1F305.4799
e(1,K,I)=0.0 ORH1F305.4800
ENDIF ORH1F305.4801
! ORH1F305.4802
ENDDO ORH1F305.4803
! ORH1F305.4804
IF (L_OCYCLIC) THEN ORH1F305.4805
tempa(1,K)=tempa(IMTM1,K) ORH1F305.4806
ELSE ORH1F305.4807
tempa(1,K)=0.0 ORH1F305.4808
ENDIF ORH1F305.4809
650 CONTINUE IPDCOFCL.375
DO 660 K=1,KM IPDCOFCL.376
DO 660 I=1,IMT IPDCOFCL.377
chkslp(i,k) = amax1(abs(e(i,k,3)),abs(tempa(i,k))) OLE1F403.2
660 CONTINUE IPDCOFCL.379
690 CONTINUE IPDCOFCL.380
IF (L_OMIXLAY) THEN ORH1F305.4810
C IPDCOFCL.383
C Set e(I,K,2) to zero in Mixed Layer IPDCOFCL.384
C IPDCOFCL.385
DO 695 K=1,KM IPDCOFCL.386
DO 695 I=1,IMT IPDCOFCL.387
IF ((mld(I)*100.) .GT. ZDZ(K)) e(I,K,2)=0. IPDCOFCL.388
695 CONTINUE IPDCOFCL.389
C IPDCOFCL.390
ENDIF ORH1F305.4811
c OLA0F401.660
c compute "fk2" OLA0F401.661
c OLA0F401.662
DO K=1,KM OLA0F401.663
DO I=1,IMT OLA0F401.664
IF (L_OISOTAPER) THEN OLA0F401.665
rloc1 = 0. OLA0F401.666
rloc2 = sign(1.,e(i,k,3))/(abs(e(i,k,3))+fxa) OLA0F401.667
slope = rloc2*sqrt(e(i,k,1)**2+e(i,k,2)**2) OLA0F401.668
if ((slope.le.0.0) .and. OLA0F401.669
& (slope.ge.(-1./(slmxr(k)*dz2rq(i,k))))) then OLA0F401.670
rloc1 = 0.5*(1.+tanh((slope+slopec)/dslope)) OLA0F401.671
endif OLA0F401.672
dciso2(i,k) = rloc1 OLA0F401.673
fk2(i,k,3) = -rloc2*e(i,k,2)*dciso2(i,k) OLA0F401.674
C dciso scales the isopycnal diffusion coefficient dependent on OLA0F401.675
C the slope of the isopycnals (see GM references from subroutine OLA0F401.676
C IPDGMVEL for details) OLA0F401.677
fk2(i,k,2)=1.*dciso2(i,k) OLA0F401.678
fk2(i,k,1)=-e(i,k,1)*e(i,k,2)*rloc2*rloc2*dciso2(i,k) OLA0F401.679
OLA0F401.680
ELSE OLA0F401.681
tempb(i,k)=fxd / (chkslp(i,k)*chkslp(i,k) + fxa) OLA0F401.682
fk2(i,k,1)=-e(I,K,1)*e(i,k,2)*tempb(I,K) OLA0F401.683
fk2(i,k,2)=e(i,k,3)*e(i,k,3)*tempb(i,k) OLA0F401.684
fk2(i,k,3)=-e(I,K,2)*e(I,K,3)*tempb(I,K) OLA0F401.685
ENDIF OLA0F401.686
ENDDO OLA0F401.687
ENDDO OLA0F401.688
OLA0F401.689
C IPDCOFCL.400
c OLA0F401.690
c in order to compute "fk3(.,.,1:3)", (zx), (zy), and (zz) components OLA0F401.691
c of the tensor, evaluate the density gradients. both "fk3" and the OLA0F401.692
c density gradients are obtained at the center of the top face of the OLA0F401.693
c "t" grid box as indicated by "X" below. OLA0F401.694
c OLA0F401.695
c lat/lon view depth/lon view OLA0F401.696
c u,v j -> +-------+ +---X---+ <- zw(k-1) OLA0F401.697
c | | | | OLA0F401.698
c t,s j -> | X | | i | <- zt(k) OLA0F401.699
c | | | | OLA0F401.700
c u,v j-1 -> +-------+ +-------+ <- zw(k) OLA0F401.701
c OLA0F401.702
c recall that the top level has been already initialized. OLA0F401.703
c note that the gradients are multiplied by fxe=1.e10 to keep the OLA0F401.704
c exponents in range. OLA0F401.705
c OLA0F401.706
DO 708 K=2,KM IPDCOFCL.401
DO 705 I=2,IMT IPDCOFCL.402
e(I,K,1)=(rx(I,K-1)+rx(I-1,K-1)+rx(I,K)+rx(I-1,K)) IPDCOFCL.403
& *CSTR(J)*DXT4RQ(I,K) IPDCOFCL.404
e(I,K,2)=(ry(I,K-1)+rym(I,K-1)+ry(I,K)+rym(I,K))*DYT4R(J) IPDCOFCL.405
e(I,K,3)=rrz(I,K)*DZZ2RQ(I,K)*fxb IPDCOFCL.406
tempb(I,K)=e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2) IPDCOFCL.407
tempa(I,K)=-SQRT(tempb(I,K))*(DZZ2RQ(I,K)*slmxr(K)) IPDCOFCL.408
705 CONTINUE IPDCOFCL.409
IF (L_OCYCLIC) THEN ORH1F305.4812
DO I = 1,3 ORH1F305.4813
e(1,K,I)=e(IMTM1,K,I) ORH1F305.4814
ENDDO ORH1F305.4815
tempa(1,K)=tempa(IMTM1,K) ORH1F305.4816
tempb(1,K)=tempb(IMTM1,K) ORH1F305.4817
ELSE ORH1F305.4818
DO I = 1,3 ORH1F305.4819
e(1,K,I)=0.0 ORH1F305.4820
ENDDO ORH1F305.4821
tempa(1,K)=0.0 ORH1F305.4822
tempb(1,K)=0.0 ORH1F305.4823
ENDIF ORH1F305.4824
! ORH1F305.4825
708 CONTINUE IPDCOFCL.434
K=1 IPDCOFCL.435
DO 710 I=2,IMT IPDCOFCL.436
e(I,K,1)=(rx(I,K)+rx(I-1,K))*2.*CSTR(J)*DXT4RQ(I,K) IPDCOFCL.437
e(I,K,2)=(ry(I,K)+rym(I,K))*2.*DYT4R(J) IPDCOFCL.438
IF (L_OEXTRAP) THEN OLA0F401.707
e(i,k,3)=0. OLA0F401.708
ELSE OLA0F401.709
e(I,K,3)=rrz(I,K)*DZZ2RQ(I,K)*fxb OLA0F401.710
ENDIF OLA0F401.711
tempb(I,K)=e(I,K,1)*e(I,K,1)+e(I,K,2)*e(I,K,2) IPDCOFCL.440
tempa(I,K)=-SQRT(tempb(I,K))*(DZZ2RQ(I,K)*slmxr(K)) IPDCOFCL.441
710 CONTINUE IPDCOFCL.442
IF (L_OCYCLIC) THEN ORH1F305.4826
DO I = 1,3 ORH1F305.4827
e(1,K,I)=e(IMTM1,K,I) ORH1F305.4828
ENDDO ORH1F305.4829
tempa(1,K)=tempa(IMTM1,K) ORH1F305.4830
tempb(1,K)=tempb(IMTM1,K) ORH1F305.4831
ELSE ORH1F305.4832
DO I = 1,3 ORH1F305.4833
e(1,K,I)=0.0 ORH1F305.4834
ENDDO ORH1F305.4835
tempa(1,K)=0.0 ORH1F305.4836
tempb(1,K)=0.0 ORH1F305.4837
ENDIF ORH1F305.4838
C IPDCOFCL.467
DO 720 K=1,KM IPDCOFCL.468
DO 720 I=1,IMT IPDCOFCL.469
chkslp(i,k) = amax1(abs(e(i,k,3)),abs(tempa(i,k))) OLE1F403.3
tempa(i,k)=fxd / (chkslp(i,k)*chkslp(i,k) + fxa) ORW1F400.31
720 CONTINUE IPDCOFCL.471
750 CONTINUE IPDCOFCL.472
IF (L_OMIXLAY) THEN ORH1F305.4839
C IPDCOFCL.475
C Set e(I,K,3) to zero in Mixed Layer IPDCOFCL.476
C IPDCOFCL.477
DO 755 K=1,KM IPDCOFCL.478
DO 755 I=1,IMT IPDCOFCL.479
IF ((mld(I)*100.) .GE. ZDZ(K)) THEN IPDCOFCL.480
e(I,K,3)=0. IPDCOFCL.481
tempb(I,K)=0. IPDCOFCL.482
ENDIF IPDCOFCL.483
755 CONTINUE IPDCOFCL.484
C IPDCOFCL.485
ENDIF ORH1F305.4840
c OLA0F401.712
c compute "fk3" OLA0F401.713
c OLA0F401.714
IF (L_OISOTAPER) THEN OLA0F401.715
DO K=2,KM OLA0F401.716
DO I=1,IMT OLA0F401.717
rloc1 = 0. OLA0F401.718
rloc2 = sign(1.,e(i,k,3))/(abs(e(i,k,3))+fxa) OLA0F401.719
slope = rloc2*sqrt(e(i,k,1)**2+e(i,k,2)**2) OLA0F401.720
if ((slope.le.0.0) .and. OLA0F401.721
& (slope.ge.(-1./(slmxr(k)*dz2rq(i,k))))) OLA0F401.722
& rloc1 = 0.5*(1.+tanh((slope+slopec)/dslope)) OLA0F401.723
rloc3 = rloc2*rloc1 OLA0F401.724
fk3(i,k,1) = -rloc3*e(i,k,1) OLA0F401.725
fk3(i,k,2) = -rloc3*e(i,k,2) OLA0F401.726
fk3(i,k,3) = rloc2*rloc2*(e(i,k,1)**2+e(i,k,2)**2)*rloc1 OLA0F401.727
ENDDO OLA0F401.728
ENDDO OLA0F401.729
ELSE OLA0F401.730
DO K=1,KM OLA0F401.731
DO I=1,IMT OLA0F401.732
fk3(i,k,3) = tempa(i,k)*tempb(i,k) OLA0F401.733
tempb(i,k) = -e(i,k,3)*tempa(i,k) OLA0F401.734
fk3(i,k,1) = e(i,k,1)*tempb(i,k) OLA0F401.735
fk3(i,k,2) = e(i,k,2)*tempb(i,k) OLA0F401.736
ENDDO OLA0F401.737
ENDDO OLA0F401.738
ENDIF OLA0F401.739
do n=1,3 ORW1F400.36
do k=1,km ORW1F400.37
do i=1,imt ORW1F400.38
if (fk1(i,k,n).gt.1.) then ORW1F400.39
fk1(i,k,n)=1. ORW1F400.40
else if (fk1(i,k,n).lt.-1.) then OLE1F403.4
fk1(i,k,n)=-1. OLE1F403.5
endif ORW1F400.41
enddo ORW1F400.42
do i=1,imt ORW1F400.43
if (fk2(i,k,n).gt.1.) then ORW1F400.44
fk2(i,k,n)=1. ORW1F400.45
else if (fk2(i,k,n).lt.-1.) then OLE1F403.6
fk2(i,k,n)=-1. OLE1F403.7
endif ORW1F400.46
enddo ORW1F400.47
do i=1,imt ORW1F400.48
if (fk3(i,k,n).gt.1.) then ORW1F400.49
fk3(i,k,n)=1. ORW1F400.50
else if (fk3(i,k,n).lt.-1.) then OLE1F403.8
fk3(i,k,n)=-1. OLE1F403.9
endif ORW1F400.51
enddo ORW1F400.52
enddo ORW1F400.53
enddo ORW1F400.54
C IPDCOFCL.500
C -------------------------------------------------------------- IPDCOFCL.501
CL Creates the diffusion matrix in its entirety IPDCOFCL.502
C -------------------------------------------------------------- IPDCOFCL.503
C IPDCOFCL.504
DO 770 L=1,3 IPDCOFCL.505
DO 770 K=1,KM IPDCOFCL.506
DO 770 I=1,IMT IPDCOFCL.507
fk1(I,K,L)=fk1(I,K,L)*ahi(K) IPDCOFCL.508
fk2(I,K,L)=fk2(I,K,L)*ahi(K) IPDCOFCL.509
fk3(I,K,L)=fk3(I,K,L)*ahi(K) IPDCOFCL.510
770 CONTINUE IPDCOFCL.511
C IPDCOFCL.512
C -------------------------------------------------------------- IPDCOFCL.513
CL To prevent horizontal roughnesses in the tracer fields IPDCOFCL.514
CL small hor. and vert. diffusion coeffs. are added to the IPDCOFCL.515
CL relevant diagonal elements of the diffusion matrix. IPDCOFCL.516
CL The epsilon term is also added in to fk3(I,K,3). IPDCOFCL.517
C -------------------------------------------------------------- IPDCOFCL.518
C IPDCOFCL.519
DO 780 K=1,KM IPDCOFCL.520
DO 780 I=1,IMT IPDCOFCL.521
fk1(I,K,1)=fk1(I,K,1)+AH IPDCOFCL.522
fk2(I,K,2)=fk2(I,K,2)+AH IPDCOFCL.523
fk3(I,K,3)=fk3(I,K,3)+gnu(I,K) IPDCOFCL.524
780 CONTINUE IPDCOFCL.525
C IPDCOFCL.526
DO K=1,KM-1 IPDCOFCL.527
DO I=1,IMT IPDCOFCL.528
FK3(I,K,3)=FK3(I,K,3)*FM(I,K) IPDCOFCL.529
END DO IPDCOFCL.530
END DO IPDCOFCL.531
C IPDCOFCL.532
790 CONTINUE IPDCOFCL.533
*ENDIF ORH1F305.458
RETURN IPDCOFCL.534
END IPDCOFCL.535
*ENDIF IPDCOFCL.536