SUBROUTINE AI_CALC ( 1,35AI_CALC.44
*CALL ARGSIZE 
AI_CALC.45
*CALL COCAWRKA AI_CALC.46
& ,j,cstr,dyur,dxur,dxtr,dz2r,ahi,csr,csjm,dyurjm,j_1, AI_CALC.47
& slmxr,dtxsqr,dyu,dxu,dzt2r,dzw,csu,dxt4r,dyt4r,dytr,cst, AI_CALC.48
& Ai_ez,Ai_nz,Ai_bx,Ai_by,K11,K22,K33,gnu, AI_CALC.49
& delta_iso,s_plus,s_minus,athkdftu,athkdftv,athkdf) AI_CALC.50
AI_CALC.51
IMPLICIT NONE AI_CALC.52
c AI_CALC.53
c======================================================================= AI_CALC.54
c compute "Ai_ez" & "K11" at the center of eastern face of T cells. AI_CALC.55
c compute "Ai_nz" & "K22" at the center of northern face of T cells. AI_CALC.56
c compute Ai_bx, Ai_by, and "K33" at the center of the bottom face AI_CALC.57
c of T cells. AI_CALC.58
c The horizontal diffusion coefficient "ah" is added into the K11 AI_CALC.59
c and K22 components, and the vertical diffusion coefficient gnu AI_CALC.60
c is added into the K33 component. AI_CALC.61
AI_CALC.62
c Differences to Redi diffusion include: AI_CALC.63
c 1. K11, K22 and K33 are scaled with the 4 surrounding slopes AI_CALC.64
c 2. Isopycnal gradients are calculated using temperature and AI_CALC.65
c salinity expansion coefficients AI_CALC.66
c 3. The (xy) and (yx) component of the isopycnal diffusion tensor AI_CALC.67
c is assumed zero, as per the small-slope approximation AI_CALC.68
c======================================================================= AI_CALC.69
c AI_CALC.70
*CALL OARRYSIZ AI_CALC.71
*CALL TYPSIZE AI_CALC.72
*CALL CNTLOCN AI_CALC.73
*CALL COCTWRKA AI_CALC.74
*CALL UMSCALAR AI_CALC.75
*CALL OTIMER AI_CALC.76
AI_CALC.77
INTEGER i,k,ip,kr,j,jq,kp1,j_1 AI_CALC.78
AI_CALC.79
C Input variables AI_CALC.80
REAL slmxr,dtxsqr(km),dzt2r(km),dzw(0:km),dxtr(imt), AI_CALC.81
& csu(jmt),dyu(jmt),dxu(imt),dxur(imt),dyur(jmt),dz2r(km), AI_CALC.82
& cst(jmt),dytr(jmt),cstr(jmt),dxt4r(imt),dyt4r(jmt), AI_CALC.83
& gnu(imt,km),ahi(km),csr(jmt),csjm,dyurjm AI_CALC.84
AI_CALC.85
C full tensor variables AI_CALC.86
REAL delta_iso,s_minus,s_plus AI_CALC.87
AI_CALC.88
C isopycnal diffusion variables AI_CALC.89
REAL Ai_ez(imt_iso,km_iso,0:1,0:1),Ai_bx(imt_iso,km_iso,0:1,0:1), AI_CALC.90
& Ai_by(imt_iso,km_iso,0:1,0:1),Ai_nz(imt_iso,km_iso,0:1,0:1), AI_CALC.91
& K11(imt_iso,km_iso),K22(imt_iso,km_iso),K33(imt_iso,km_iso) AI_CALC.92
REAL athkdf(km),athkdftu_mom(imt,km),athkdftv_mom(imt,km) AI_CALC.93
REAL athkdftu(imt,jmt_vis),athkdftv(imt,jmt_vis) AI_CALC.94
AI_CALC.95
C Local variables AI_CALC.96
REAL sc,c1,c2,dzt4r,p5,Ai0,p25, AI_CALC.97
& sumx,sumy,sumz,c0,sxe,epsln,facty,syn,sxb, AI_CALC.98
& syb,sumxz(imt,km),sumyz(imt,km) AI_CALC.99
AI_CALC.100
! Arrays used to help speed up tanh function calculations AI_CALC.101
REAL a_sxe(imt) AI_CALC.102
& ,a_sxb1(imt),a_sxb2(imt),a_sxb3(imt),a_sxb4(imt) AI_CALC.103
& ,a_syb1(imt),a_syb2(imt),a_syb3(imt),a_syb4(imt) AI_CALC.104
AI_CALC.105
REAL csjm_loc,dyurjm_loc AI_CALC.106
AI_CALC.107
C statement functions AI_CALC.108
REAL drodxe,drodze,drodyn,drodzn,drodxb,drodyb,drodzb, AI_CALC.109
& drodye,drodxn AI_CALC.110
AI_CALC.111
c----------------------------------------------------------------------- AI_CALC.112
c compute "Ai_ez" on east face of T cell. Note re-scaling factor AI_CALC.113
c which reduces mixing coefficient "Ai" where slope "sxe" AI_CALC.114
c exceeds the critical slope "sc" for the small slope approx. For AI_CALC.115
c the full tensor, it is re-scaled if outside the stable range. AI_CALC.116
c----------------------------------------------------------------------- AI_CALC.117
c AI_CALC.118
c AI_CALC.119
c statement functions to calculate density gradients on various AI_CALC.120
c faces of tracer boxes AI_CALC.121
c AI_CALC.122
drodxe(i,k,ip) = alphai(i+ip,k,0)*ddxt(i,k,1,0) + AI_CALC.123
& betai(i+ip,k,0)*ddxt(i,k,2,0) AI_CALC.124
drodze(i,k,ip,kr) = alphai(i+ip,k,0)*ddzt(i+ip,k-1+kr,1,0) + AI_CALC.125
& betai(i+ip,k,0)*ddzt(i+ip,k-1+kr,2,0) AI_CALC.126
c AI_CALC.127
drodyn(i,k,jq) = alphai(i,k,jq)*ddyt(i,k,1,1) + AI_CALC.128
& betai(i,k,jq)*ddyt(i,k,2,1) AI_CALC.129
drodzn(i,k,jq,kr) = alphai(i,k,jq)*ddzt(i,k-1+kr,1,jq) + AI_CALC.130
& betai(i,k,jq)*ddzt(i,k-1+kr,2,jq) AI_CALC.131
c AI_CALC.132
drodxb(i,k,ip,kr) = alphai(i,k+kr,0)*ddxt(i-1+ip,k+kr,1,0) + AI_CALC.133
& betai(i,k+kr,0)*ddxt(i-1+ip,k+kr,2,0) AI_CALC.134
drodyb(i,k,jq,kr) = alphai(i,k+kr,0)*ddyt(i,k+kr,1,jq) + AI_CALC.135
& betai(i,k+kr,0)*ddyt(i,k+kr,2,jq) AI_CALC.136
drodzb(i,k,kr) = alphai(i,k+kr,0)*ddzt(i,k,1,0) + AI_CALC.137
& betai(i,k+kr,0)*ddzt(i,k,2,0) AI_CALC.138
c AI_CALC.139
c statement functions only needed with full isopycnal diffusion AI_CALC.140
c tensor (as opposed to small-slope approx) AI_CALC.141
c AI_CALC.142
drodye(i,k,ip,jq) = alphai(i+ip,k,jq)*ddyt(i+ip,k,1,jq) + AI_CALC.143
& betai(i+ip,k,jq)*ddyt(i+ip,k,2,jq) AI_CALC.144
drodxn(i,k,ip,jq) = alphai(i,k,jq)*ddxt(i-1+ip,k,1,jq) + AI_CALC.145
& betai(i,k,jq)*ddxt(i-1+ip,k,2,jq) AI_CALC.146
AI_CALC.147
c0=0. AI_CALC.148
c1=1. AI_CALC.149
c2=2. AI_CALC.150
p5=0.5 AI_CALC.151
p25=0.25 AI_CALC.152
epsln=1.0e-25 AI_CALC.153
AI_CALC.154
do k=1,km AI_CALC.155
do i=1,imt AI_CALC.156
sumxz(i,k)=c0 AI_CALC.157
sumyz(i,k)=c0 AI_CALC.158
enddo AI_CALC.159
enddo AI_CALC.160
AI_CALC.161
IF (.NOT.L_OVISBECK) THEN AI_CALC.162
do k=1,km AI_CALC.163
do i=1,imt AI_CALC.164
athkdftu_mom(i,k)=athkdf(k) AI_CALC.165
athkdftv_mom(i,k)=athkdf(k) AI_CALC.166
enddo AI_CALC.167
enddo AI_CALC.168
ELSE AI_CALC.169
do k=1,km AI_CALC.170
do i=1,imt AI_CALC.171
athkdftu_mom(i,k)=athkdftu(i,j) AI_CALC.172
athkdftv_mom(i,k)=athkdftv(i,j) AI_CALC.173
enddo AI_CALC.174
enddo AI_CALC.175
ENDIF ! L_OVISBECK AI_CALC.176
AI_CALC.177
C if we are on the initialisation row of a block, need to use AI_CALC.178
C the value of cs(j-1) and dyur(j-1) sent from the block below AI_CALC.179
if (j.lt.j_1) then AI_CALC.180
csjm_loc=csjm AI_CALC.181
dyurjm_loc=dyurjm AI_CALC.182
else AI_CALC.183
csjm_loc=csu(j-1) AI_CALC.184
dyurjm_loc=dyur(j-1) AI_CALC.185
endif AI_CALC.186
AI_CALC.187
IF (L_OISOTAPER) THEN AI_CALC.188
do k=1,km AI_CALC.189
c Ai0 = ahi(k) AI_CALC.190
Ai0 = c1 AI_CALC.191
do i=2,imtm1 AI_CALC.192
sxe = abs(drodxe(i,k,0)/(drodze(i,k,0,0)+epsln)) AI_CALC.193
a_sxe(i) =(sxe-slopec)/dslope AI_CALC.194
enddo AI_CALC.195
AI_CALC.196
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.197
AI_CALC.198
do i=2,imtm1 AI_CALC.199
Ai_ez(i,k,0,0) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.200
& *p5*(c1-a_sxe(i)) AI_CALC.201
enddo AI_CALC.202
AI_CALC.203
do i=2,imtm1 AI_CALC.204
sxe = abs(drodxe(i,k,0)/(drodze(i,k,0,1)+epsln)) AI_CALC.205
a_sxe(i) =(sxe-slopec)/dslope AI_CALC.206
enddo AI_CALC.207
AI_CALC.208
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.209
AI_CALC.210
do i=2,imtm1 AI_CALC.211
Ai_ez(i,k,0,1) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.212
& *p5*(c1-a_sxe(i)) AI_CALC.213
enddo AI_CALC.214
AI_CALC.215
do i=2,imtm1 AI_CALC.216
sxe = abs(drodxe(i,k,1)/(drodze(i,k,1,0)+epsln)) AI_CALC.217
a_sxe(i) =(sxe-slopec)/dslope AI_CALC.218
enddo AI_CALC.219
AI_CALC.220
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.221
AI_CALC.222
do i=2,imtm1 AI_CALC.223
Ai_ez(i,k,1,0) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.224
& *p5*(c1-a_sxe(i)) AI_CALC.225
enddo AI_CALC.226
AI_CALC.227
do i=2,imtm1 AI_CALC.228
sxe = abs(drodxe(i,k,1)/(drodze(i,k,1,1)+epsln)) AI_CALC.229
a_sxe(i) =(sxe-slopec)/dslope AI_CALC.230
enddo AI_CALC.231
AI_CALC.232
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.233
AI_CALC.234
do i=2,imtm1 AI_CALC.235
Ai_ez(i,k,1,1) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.236
& *p5*(c1-a_sxe(i)) AI_CALC.237
enddo AI_CALC.238
enddo AI_CALC.239
AI_CALC.240
ELSE AI_CALC.241
AI_CALC.242
do kr=0,1 AI_CALC.243
do ip=0,1 AI_CALC.244
do k=1,km AI_CALC.245
c Ai0 = ahi(k) AI_CALC.246
Ai0 = c1 AI_CALC.247
sc = c1/(slmxr*dtxsqr(k)) AI_CALC.248
do i=2,imtm1 AI_CALC.249
sumz = c0 AI_CALC.250
sxe = abs(drodxe(i,k,ip)/(drodze(i,k,ip,kr)+epsln)) AI_CALC.251
if (sxe .gt. sc) then AI_CALC.252
Ai_ez(i,k,ip,kr) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.253
& *(sc/(sxe + epsln))**2 AI_CALC.254
else AI_CALC.255
Ai_ez(i,k,ip,kr) = Ai0*fm(i,k)*fm(i+1,k) AI_CALC.256
endif AI_CALC.257
enddo AI_CALC.258
enddo AI_CALC.259
enddo AI_CALC.260
enddo AI_CALC.261
ENDIF ! TAPER AI_CALC.262
AI_CALC.263
do k=1,km AI_CALC.264
dzt4r = p5*dzt2r(k) AI_CALC.265
do i=2,imtm1 AI_CALC.266
K11(i,k) = (dzt4r*(dzw(k-1)*Ai_ez(i,k,0,0)+ AI_CALC.267
& dzw(k)*Ai_ez(i,k,0,1)+dzw(k-1)*Ai_ez(i,k,1,0)+ AI_CALC.268
& dzw(k)*Ai_ez(i,k,1,1)))*ahi(k) + ah AI_CALC.269
enddo AI_CALC.270
enddo AI_CALC.271
AI_CALC.272
call SETBCX (Ai_ez(1,1,0,0), imt, km) AI_CALC.273
call SETBCX (Ai_ez(1,1,1,0), imt, km) AI_CALC.274
call SETBCX (Ai_ez(1,1,0,1), imt, km) AI_CALC.275
call SETBCX (Ai_ez(1,1,1,1), imt, km) AI_CALC.276
call SETBCX (K11(1,1), imt, km) AI_CALC.277
AI_CALC.278
c AI_CALC.279
c AI_CALC.280
c--------------------------------------------------------------------- AI_CALC.281
c compute "Ai_nz" on north face of T cell. Note re-scaling factor AI_CALC.282
c which reduces mixing coefficient "Ai" where slope "syn" AI_CALC.283
c exceeds the critical slope "sc" for the small slope approx. For AI_CALC.284
c the full tensor, it is re-scaled if outside the stable range. AI_CALC.285
c--------------------------------------------------------------------- AI_CALC.286
c AI_CALC.287
IF (L_OISOTAPER) THEN AI_CALC.288
do k=1,km AI_CALC.289
c Ai0 = ahi(k) AI_CALC.290
Ai0 = c1 AI_CALC.291
do i=2,imtm1 AI_CALC.292
syn =abs(drodyn(i,k,0)/(drodzn(i,k,0,0)+epsln)) AI_CALC.293
a_sxe(i) =(syn-slopec)/dslope AI_CALC.294
enddo AI_CALC.295
AI_CALC.296
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.297
AI_CALC.298
do i=2,imtm1 AI_CALC.299
Ai_nz(i,k,0,0) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.300
& *p5*(c1-a_sxe(i)) AI_CALC.301
enddo AI_CALC.302
AI_CALC.303
do i=2,imtm1 AI_CALC.304
syn =abs(drodyn(i,k,0)/(drodzn(i,k,0,1)+epsln)) AI_CALC.305
a_sxe(i) =(syn-slopec)/dslope AI_CALC.306
enddo AI_CALC.307
AI_CALC.308
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.309
AI_CALC.310
do i=2,imtm1 AI_CALC.311
Ai_nz(i,k,0,1) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.312
& *p5*(c1-a_sxe(i)) AI_CALC.313
enddo AI_CALC.314
AI_CALC.315
do i=2,imtm1 AI_CALC.316
syn =abs(drodyn(i,k,1)/(drodzn(i,k,1,0)+epsln)) AI_CALC.317
a_sxe(i) =(syn-slopec)/dslope AI_CALC.318
enddo AI_CALC.319
AI_CALC.320
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.321
AI_CALC.322
do i=2,imtm1 AI_CALC.323
Ai_nz(i,k,1,0) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.324
& *p5*(c1-a_sxe(i)) AI_CALC.325
enddo AI_CALC.326
AI_CALC.327
do i=2,imtm1 AI_CALC.328
syn =abs(drodyn(i,k,1)/(drodzn(i,k,1,1)+epsln)) AI_CALC.329
a_sxe(i) =(syn-slopec)/dslope AI_CALC.330
enddo AI_CALC.331
AI_CALC.332
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.333
AI_CALC.334
do i=2,imtm1 AI_CALC.335
Ai_nz(i,k,1,1) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.336
& *p5*(c1-a_sxe(i)) AI_CALC.337
enddo AI_CALC.338
enddo AI_CALC.339
AI_CALC.340
ELSE AI_CALC.341
AI_CALC.342
do kr=0,1 AI_CALC.343
do jq=0,1 AI_CALC.344
do k=1,km AI_CALC.345
c Ai0 = ahi(k) AI_CALC.346
Ai0 = c1 AI_CALC.347
sc = c1/(slmxr*dtxsqr(k)) AI_CALC.348
dzt4r = p5*dzt2r(k) AI_CALC.349
do i=2,imtm1 AI_CALC.350
sumz = c0 AI_CALC.351
syn =abs(drodyn(i,k,jq)/(drodzn(i,k,jq,kr)+epsln)) AI_CALC.352
if (syn .gt. sc) then AI_CALC.353
Ai_nz(i,k,jq,kr) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.354
& *(sc/(syn + epsln))**2 AI_CALC.355
else AI_CALC.356
Ai_nz(i,k,jq,kr) = Ai0*fm(i,k)*fmp(i,k) AI_CALC.357
endif AI_CALC.358
enddo AI_CALC.359
enddo AI_CALC.360
enddo AI_CALC.361
enddo AI_CALC.362
ENDIF ! TAPER AI_CALC.363
AI_CALC.364
AI_CALC.365
do k=1,km AI_CALC.366
dzt4r = p5*dzt2r(k) AI_CALC.367
do i=2,imtm1 AI_CALC.368
K22(i,k) = (dzt4r*(dzw(k-1)*Ai_nz(i,k,0,0) + AI_CALC.369
& dzw(k)*Ai_nz(i,k,0,1) + dzw(k-1)*Ai_nz(i,k,1,0) + AI_CALC.370
& dzw(k)*Ai_nz(i,k,1,1)))*ahi(k) + ah AI_CALC.371
enddo AI_CALC.372
enddo AI_CALC.373
AI_CALC.374
call SETBCX (Ai_nz(1,1,0,0), imt, km) AI_CALC.375
call SETBCX (Ai_nz(1,1,1,0), imt, km) AI_CALC.376
call SETBCX (Ai_nz(1,1,0,1), imt, km) AI_CALC.377
call SETBCX (Ai_nz(1,1,1,1), imt, km) AI_CALC.378
call SETBCX (K22(1,1), imt, km) AI_CALC.379
AI_CALC.380
c--------------------------------------------------------------------- AI_CALC.381
c compute "Ai_bx", "Ai_by", & K33 on bottom face of T cell. Note AI_CALC.382
c re-scaling factor to reduce mixing coefficient "Ai" where slopes AI_CALC.383
c "sxb" and "syb" exceeds the critical slope "sc" for the small AI_CALC.384
c slope approx. For the full tensor, it is re-scaled if outside the AI_CALC.385
c stable range. AI_CALC.386
c--------------------------------------------------------------------- AI_CALC.387
c AI_CALC.388
c AI_CALC.389
c eastward slopes at the base of T cells AI_CALC.390
c AI_CALC.391
IF (L_OISOTAPER) THEN AI_CALC.392
do k=1,kmm1 AI_CALC.393
c Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.394
Ai0 = c1 AI_CALC.395
do i=2,imtm1 AI_CALC.396
a_sxb1(i) = abs(drodxb(i,k,0,0)/(drodzb(i,k,0)+epsln)) AI_CALC.397
a_sxe(i) =(a_sxb1(i)-slopec)/dslope AI_CALC.398
enddo AI_CALC.399
AI_CALC.400
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.401
AI_CALC.402
do i=2,imtm1 AI_CALC.403
Ai_bx(i,k,0,0) = Ai0*fm(i,k+1) AI_CALC.404
& *p5*(c1-a_sxe(i)) AI_CALC.405
enddo AI_CALC.406
AI_CALC.407
do i=2,imtm1 AI_CALC.408
a_sxb2(i) = abs(drodxb(i,k,0,1)/(drodzb(i,k,1)+epsln)) AI_CALC.409
a_sxe(i) =(a_sxb2(i)-slopec)/dslope AI_CALC.410
enddo AI_CALC.411
AI_CALC.412
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.413
AI_CALC.414
do i=2,imtm1 AI_CALC.415
Ai_bx(i,k,0,1) = Ai0*fm(i,k+1) AI_CALC.416
& *p5*(c1-a_sxe(i)) AI_CALC.417
enddo AI_CALC.418
AI_CALC.419
do i=2,imtm1 AI_CALC.420
a_sxb3(i) = abs(drodxb(i,k,1,0)/(drodzb(i,k,0)+epsln)) AI_CALC.421
a_sxe(i) =(a_sxb3(i)-slopec)/dslope AI_CALC.422
enddo AI_CALC.423
AI_CALC.424
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.425
AI_CALC.426
do i=2,imtm1 AI_CALC.427
Ai_bx(i,k,1,0) = Ai0*fm(i,k+1) AI_CALC.428
& *p5*(c1-a_sxe(i)) AI_CALC.429
enddo AI_CALC.430
AI_CALC.431
do i=2,imtm1 AI_CALC.432
a_sxb4(i) = abs(drodxb(i,k,1,1)/(drodzb(i,k,1)+epsln)) AI_CALC.433
a_sxe(i) =(a_sxb4(i)-slopec)/dslope AI_CALC.434
enddo AI_CALC.435
AI_CALC.436
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.437
AI_CALC.438
do i=2,imtm1 AI_CALC.439
Ai_bx(i,k,1,1) = Ai0*fm(i,k+1) AI_CALC.440
& *p5*(c1-a_sxe(i)) AI_CALC.441
enddo AI_CALC.442
AI_CALC.443
do i=2,imtm1 AI_CALC.444
sumxz(i,k) = dxu(i-1)*(Ai_bx(i,k,0,0)*a_sxb1(i)**2+ AI_CALC.445
& Ai_bx(i,k,0,1)*a_sxb2(i)**2)+ AI_CALC.446
& dxu(i)*(Ai_bx(i,k,1,0)*a_sxb3(i)**2+ AI_CALC.447
& Ai_bx(i,k,1,1)*a_sxb4(i)**2) AI_CALC.448
enddo AI_CALC.449
enddo AI_CALC.450
AI_CALC.451
ELSE AI_CALC.452
AI_CALC.453
do ip=0,1 AI_CALC.454
do kr=0,1 AI_CALC.455
do k=1,kmm1 AI_CALC.456
c Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.457
Ai0 = c1 AI_CALC.458
sc = c1/(slmxr*dtxsqr(k)) AI_CALC.459
kp1 = min(k+1,km) AI_CALC.460
do i=2,imtm1 AI_CALC.461
sxb = abs(drodxb(i,k,ip,kr)/(drodzb(i,k,kr)+epsln)) AI_CALC.462
if (sxb .gt. sc) then AI_CALC.463
Ai_bx(i,k,ip,kr) = Ai0*fm(i,k+1) AI_CALC.464
& *(sc/(sxb + epsln))**2 AI_CALC.465
else AI_CALC.466
Ai_bx(i,k,ip,kr) = Ai0*fm(i,k+1) AI_CALC.467
endif AI_CALC.468
sumxz(i,k) = sumxz(i,k) + dxu(i-1+ip)*Ai_bx(i,k,ip,kr) AI_CALC.469
& *sxb**2 AI_CALC.470
enddo AI_CALC.471
enddo AI_CALC.472
enddo AI_CALC.473
enddo AI_CALC.474
ENDIF ! taper AI_CALC.475
AI_CALC.476
c AI_CALC.477
c northward slopes at the base of T cells AI_CALC.478
c AI_CALC.479
IF (L_OISOTAPER) THEN AI_CALC.480
do k=1,kmm1 AI_CALC.481
c Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.482
Ai0 = c1 AI_CALC.483
do i=2,imtm1 AI_CALC.484
a_syb1(i) = abs(drodyb(i,k,0,0)/(drodzb(i,k,0)+epsln)) AI_CALC.485
a_sxe(i) =(a_syb1(i)-slopec)/dslope AI_CALC.486
enddo AI_CALC.487
AI_CALC.488
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.489
AI_CALC.490
do i=2,imtm1 AI_CALC.491
Ai_by(i,k,0,0) = Ai0*fm(i,k+1) AI_CALC.492
& *p5*(c1-a_sxe(i)) AI_CALC.493
enddo AI_CALC.494
AI_CALC.495
do i=2,imtm1 AI_CALC.496
a_syb2(i) = abs(drodyb(i,k,0,1)/(drodzb(i,k,1)+epsln)) AI_CALC.497
a_sxe(i) =(a_syb2(i)-slopec)/dslope AI_CALC.498
enddo AI_CALC.499
AI_CALC.500
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.501
AI_CALC.502
do i=2,imtm1 AI_CALC.503
Ai_by(i,k,0,1) = Ai0*fm(i,k+1) AI_CALC.504
& *p5*(c1-a_sxe(i)) AI_CALC.505
enddo AI_CALC.506
AI_CALC.507
do i=2,imtm1 AI_CALC.508
a_syb3(i) = abs(drodyb(i,k,1,0)/(drodzb(i,k,0)+epsln)) AI_CALC.509
a_sxe(i) =(a_syb3(i)-slopec)/dslope AI_CALC.510
enddo AI_CALC.511
AI_CALC.512
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.513
AI_CALC.514
do i=2,imtm1 AI_CALC.515
Ai_by(i,k,1,0) = Ai0*fm(i,k+1) AI_CALC.516
& *p5*(c1-a_sxe(i)) AI_CALC.517
enddo AI_CALC.518
AI_CALC.519
do i=2,imtm1 AI_CALC.520
a_syb4(i) = abs(drodyb(i,k,1,1)/(drodzb(i,k,1)+epsln)) AI_CALC.521
a_sxe(i) =(a_syb4(i)-slopec)/dslope AI_CALC.522
enddo AI_CALC.523
AI_CALC.524
call fast_tanh(imtm1-1,a_sxe(2)) AI_CALC.525
AI_CALC.526
do i=2,imtm1 AI_CALC.527
Ai_by(i,k,1,1) = Ai0*fm(i,k+1) AI_CALC.528
& *p5*(c1-a_sxe(i)) AI_CALC.529
enddo AI_CALC.530
AI_CALC.531
do i=2,imtm1 AI_CALC.532
AI_CALC.533
c sumyz(i,k) = (1./(csr(j-1)*dyur(j-1)))*(Ai_by(i,k,0,0)*syb1**2+ AI_CALC.534
sumyz(i,k) = (csjm_loc/dyurjm_loc)* AI_CALC.535
& (Ai_by(i,k,0,0)*a_syb1(i)**2+ AI_CALC.536
& Ai_by(i,k,0,1)*a_syb2(i)**2)+csu(j)*dyu(j)* AI_CALC.537
& (Ai_by(i,k,1,0)*a_syb3(i)**2+ AI_CALC.538
& Ai_by(i,k,1,1)*a_syb4(i)**2) AI_CALC.539
enddo AI_CALC.540
AI_CALC.541
enddo AI_CALC.542
AI_CALC.543
ELSE AI_CALC.544
AI_CALC.545
do jq=0,1 AI_CALC.546
facty = 1./(csr(j-1+jq)*dyur(j-1+jq)) AI_CALC.547
do kr=0,1 AI_CALC.548
do k=1,kmm1 AI_CALC.549
c Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.550
Ai0 = c1 AI_CALC.551
sc = c1/(slmxr*dtxsqr(k)) AI_CALC.552
kp1 = min(k+1,km) AI_CALC.553
do i=2,imtm1 AI_CALC.554
c Ai0 = athkdft(i)*kappa AI_CALC.555
syb = abs(drodyb(i,k,jq,kr)/(drodzb(i,k,kr)+epsln)) AI_CALC.556
if (syb .gt. sc) then AI_CALC.557
Ai_by(i,k,jq,kr) = Ai0*fm(i,k+1) AI_CALC.558
& *(sc/(syb + epsln))**2 AI_CALC.559
else AI_CALC.560
Ai_by(i,k,jq,kr) = Ai0*fm(i,k+1) AI_CALC.561
endif AI_CALC.562
sumyz(i,k) = sumyz(i,k) + facty*Ai_by(i,k,jq,kr) AI_CALC.563
& *syb**2 AI_CALC.564
enddo AI_CALC.565
enddo AI_CALC.566
enddo AI_CALC.567
enddo AI_CALC.568
ENDIF ! TAPER AI_CALC.569
AI_CALC.570
C the calculations take place offset one level in the MOM code AI_CALC.571
C from our code. So redefining K33 here to be at the TOP of tracer box AI_CALC.572
C k (HC), as opposed to the bottom of grid box k. AI_CALC.573
C Add vertical diffusion gnu into vertical part of isopycnal diff, AI_CALC.574
C for implicit solving AI_CALC.575
AI_CALC.576
do k=2,km AI_CALC.577
Ai0 = p5*(ahi(k-1)+ahi(k)) AI_CALC.578
do i=2,imtm1 AI_CALC.579
K33(i,k) = (dxt4r(i)*sumxz(i,k-1) + AI_CALC.580
& dyt4r(j)*cstr(j)*sumyz(i,k-1))*Ai0 + gnu(i,k) AI_CALC.581
AI_CALC.582
enddo AI_CALC.583
enddo AI_CALC.584
do i=2,imtm1 AI_CALC.585
K33(i,1) = gnu(i,1) AI_CALC.586
enddo AI_CALC.587
AI_CALC.588
call SETBCX (Ai_bx(1,1,1,0), imt, km) AI_CALC.589
call SETBCX (Ai_bx(1,1,0,0), imt, km) AI_CALC.590
call SETBCX (Ai_bx(1,1,1,1), imt, km) AI_CALC.591
call SETBCX (Ai_bx(1,1,0,1), imt, km) AI_CALC.592
call SETBCX (Ai_by(1,1,1,0), imt, km) AI_CALC.593
call SETBCX (Ai_by(1,1,0,0), imt, km) AI_CALC.594
call SETBCX (Ai_by(1,1,1,1), imt, km) AI_CALC.595
call SETBCX (Ai_by(1,1,0,1), imt, km) AI_CALC.596
call SETBCX (K33(1,1), imt, km) AI_CALC.597
AI_CALC.598
C the skew-diffusive tensor form for the Gent and McWilliams scheme AI_CALC.599
C is implemented by setting the K13 and K23 diffusion coefficients to AI_CALC.600
C zero, and doubling the K31 and K32 coefficients. This assumes that AI_CALC.601
C one is choosing the same value for isopycnal and thickness diffusion AI_CALC.602
AI_CALC.603
IF (L_OISOGMSKEW) THEN AI_CALC.604
do k=1,km AI_CALC.605
do i=1,imt AI_CALC.606
Ai_ez(i,k,0,0)=(ahi(k)-athkdftu_mom(i,k))*Ai_ez(i,k,0,0) AI_CALC.607
Ai_ez(i,k,0,1)=(ahi(k)-athkdftu_mom(i,k))*Ai_ez(i,k,0,1) AI_CALC.608
Ai_ez(i,k,1,0)=(ahi(k)-athkdftu_mom(i,k))*Ai_ez(i,k,1,0) AI_CALC.609
Ai_ez(i,k,1,1)=(ahi(k)-athkdftu_mom(i,k))*Ai_ez(i,k,1,1) AI_CALC.610
AI_CALC.611
Ai_nz(i,k,0,0)=(ahi(k)-athkdftv_mom(i,k))*Ai_nz(i,k,0,0) AI_CALC.612
Ai_nz(i,k,0,1)=(ahi(k)-athkdftv_mom(i,k))*Ai_nz(i,k,0,1) AI_CALC.613
Ai_nz(i,k,1,0)=(ahi(k)-athkdftv_mom(i,k))*Ai_nz(i,k,1,0) AI_CALC.614
Ai_nz(i,k,1,1)=(ahi(k)-athkdftv_mom(i,k))*Ai_nz(i,k,1,1) AI_CALC.615
enddo AI_CALC.616
enddo AI_CALC.617
do k=1,kmm1 AI_CALC.618
Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.619
do i=1,imt AI_CALC.620
Ai_bx(i,k,0,0)=(Ai0+athkdftu_mom(i,k))*Ai_bx(i,k,0,0) AI_CALC.621
Ai_bx(i,k,0,1)=(Ai0+athkdftu_mom(i,k))*Ai_bx(i,k,0,1) AI_CALC.622
Ai_bx(i,k,1,0)=(Ai0+athkdftu_mom(i,k))*Ai_bx(i,k,1,0) AI_CALC.623
Ai_bx(i,k,1,1)=(Ai0+athkdftu_mom(i,k))*Ai_bx(i,k,1,1) AI_CALC.624
AI_CALC.625
Ai_by(i,k,0,0)=(Ai0+athkdftv_mom(i,k))*Ai_by(i,k,0,0) AI_CALC.626
Ai_by(i,k,0,1)=(Ai0+athkdftv_mom(i,k))*Ai_by(i,k,0,1) AI_CALC.627
Ai_by(i,k,1,0)=(Ai0+athkdftv_mom(i,k))*Ai_by(i,k,1,0) AI_CALC.628
Ai_by(i,k,1,1)=(Ai0+athkdftv_mom(i,k))*Ai_by(i,k,1,1) AI_CALC.629
enddo AI_CALC.630
enddo AI_CALC.631
AI_CALC.632
ELSE AI_CALC.633
AI_CALC.634
do k=1,km AI_CALC.635
do i=1,imt AI_CALC.636
Ai_ez(i,k,0,0)=ahi(k)*Ai_ez(i,k,0,0) AI_CALC.637
Ai_ez(i,k,0,1)=ahi(k)*Ai_ez(i,k,0,1) AI_CALC.638
Ai_ez(i,k,1,0)=ahi(k)*Ai_ez(i,k,1,0) AI_CALC.639
Ai_ez(i,k,1,1)=ahi(k)*Ai_ez(i,k,1,1) AI_CALC.640
AI_CALC.641
Ai_nz(i,k,0,0)=ahi(k)*Ai_nz(i,k,0,0) AI_CALC.642
Ai_nz(i,k,0,1)=ahi(k)*Ai_nz(i,k,0,1) AI_CALC.643
Ai_nz(i,k,1,0)=ahi(k)*Ai_nz(i,k,1,0) AI_CALC.644
Ai_nz(i,k,1,1)=ahi(k)*Ai_nz(i,k,1,1) AI_CALC.645
enddo AI_CALC.646
enddo AI_CALC.647
do k=1,kmm1 AI_CALC.648
Ai0 = p5*(ahi(k)+ahi(k+1)) AI_CALC.649
do i=1,imt AI_CALC.650
Ai_bx(i,k,0,0)=Ai0*Ai_bx(i,k,0,0) AI_CALC.651
Ai_bx(i,k,0,1)=Ai0*Ai_bx(i,k,0,1) AI_CALC.652
Ai_bx(i,k,1,0)=Ai0*Ai_bx(i,k,1,0) AI_CALC.653
Ai_bx(i,k,1,1)=Ai0*Ai_bx(i,k,1,1) AI_CALC.654
AI_CALC.655
Ai_by(i,k,0,0)=Ai0*Ai_by(i,k,0,0) AI_CALC.656
Ai_by(i,k,0,1)=Ai0*Ai_by(i,k,0,1) AI_CALC.657
Ai_by(i,k,1,0)=Ai0*Ai_by(i,k,1,0) AI_CALC.658
Ai_by(i,k,1,1)=Ai0*Ai_by(i,k,1,1) AI_CALC.659
enddo AI_CALC.660
enddo AI_CALC.661
AI_CALC.662
ENDIF ! skew GM scheme AI_CALC.663
AI_CALC.664
c AI_CALC.665
RETURN AI_CALC.666
END AI_CALC.667
AI_CALC.668
SUBROUTINE FAST_TANH(ISIZE,ARRAY) 23AI_CALC.669
! AI_CALC.670
! Author : R.Hill AI_CALC.671
! Date : 26.10.98 AI_CALC.672
! Version : 4.5 AI_CALC.673
! AI_CALC.674
! Purpose : Calculates tanh functions (used in certain ocean AI_CALC.675
! model routines - eg isopyc_a, ai_calc) more quickly AI_CALC.676
! than using direct references to the tanh function. AI_CALC.677
! The incoming data are the size of the array for which AI_CALC.678
! tanh is to be calculated and the array itself. AI_CALC.679
! AI_CALC.680
! AI_CALC.681
! Modification History: AI_CALC.682
! Version Date Details AI_CALC.683
! ------- ------- ------------------------------------------ AI_CALC.684
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! AI_CALC.685
AI_CALC.686
IMPLICIT NONE AI_CALC.687
AI_CALC.688
INTEGER ISIZE ! IN: size of array to work on AI_CALC.689
&, I ! Local: Loop index AI_CALC.690
AI_CALC.691
REAL ARRAY(ISIZE) ! IN/OUT : The array to work with AI_CALC.692
AI_CALC.693
! First calculate the exponent part AI_CALC.694
DO I = 1, ISIZE AI_CALC.695
ARRAY(I) = -2.*ARRAY(I) AI_CALC.696
ENDDO AI_CALC.697
AI_CALC.698
*IF DEF,VECTLIB PXVECTLB.3
! Use vector exponent if available - for speed. AI_CALC.700
CALL EXP_V(ISIZE,ARRAY,ARRAY) AI_CALC.701
*ELSE AI_CALC.702
DO I = 1, ISIZE AI_CALC.703
ARRAY(I) = EXP(ARRAY(I)) AI_CALC.704
ENDDO AI_CALC.705
*ENDIF AI_CALC.706
AI_CALC.707
! Now compute the tanh values AI_CALC.708
DO I = 1, ISIZE AI_CALC.709
ARRAY(I)= (1.-ARRAY(I))/(1.+ARRAY(I)) AI_CALC.710
ENDDO AI_CALC.711
AI_CALC.712
RETURN AI_CALC.713
END AI_CALC.714
*ENDIF AI_CALC.715
AI_CALC.716
AI_CALC.717