SUBROUTINE ISOPYC_A ( 1,14ISOPYC_A.42
*CALL ARGSIZE 
ISOPYC_A.43
*CALL COCAWRKA ISOPYC_A.44
& ,j,dtxsqr,slmxr,dz2r,csu,dzt,athkdftu,athkdftv,athkdf, ISOPYC_A.45
& cstr,dxtr,dytr,adv_vetiso,adv_vbtiso,slopec,dslope,dxur,dyur, ISOPYC_A.46
& cst,csr,J_OFFSET,athkdf_bi,srho) ISOPYC_A.47
c ISOPYC_A.48
c======================================================================= ISOPYC_A.49
c compute isopycnal transport velocities. ISOPYC_A.50
c======================================================================= ISOPYC_A.51
c ISOPYC_A.52
IMPLICIT NONE ISOPYC_A.53
ISOPYC_A.54
*CALL OARRYSIZ ISOPYC_A.55
*CALL TYPSIZE ISOPYC_A.56
*CALL CNTLOCN ISOPYC_A.57
*CALL COCTWRKA ISOPYC_A.58
ISOPYC_A.59
C Input variables ISOPYC_A.60
INTEGER j,J_OFFSET ISOPYC_A.61
REAL dtxsqr(km),slmxr,dz2r(km),csu(jmt),dzt(km), ISOPYC_A.62
& cstr(jmt),dxtr(imt),dytr(jmt),slopec,dslope,dxur(imt), ISOPYC_A.63
& dyur(jmt),cst(jmt),csr(jmt),kappa ISOPYC_A.64
REAL athkdftu(imt_vis,jmt_vis),athkdftv(imt_vis,jmt_vis) ISOPYC_A.65
REAL athkdf(km),athkdftu_bi,athkdftv_bi,athkdf_bi ISOPYC_A.66
REAL athkdftu_mom(imt,km),athkdftv_mom(imt,km) ISOPYC_A.67
c REAL fmpp(imt,km) ISOPYC_A.68
ISOPYC_A.69
C output variables ISOPYC_A.70
REAL adv_vetiso(imt_gmm,km_gmm), ! GM eastern velocity ISOPYC_A.71
& adv_vbtiso(imt_gmm,0:km_gmm) ! GM vertical velocity ISOPYC_A.72
REAL srho(imt,km) ISOPYC_A.73
ISOPYC_A.74
C Local variables ISOPYC_A.75
INTEGER k,km1,kp1,i ISOPYC_A.76
REAL sc,Ath0,at,bt,epsln,ab,bb,absstn,abssbn, ISOPYC_A.77
& c0,c1,p5,absste,abssbe, ISOPYC_A.78
& ath_t(imt,km), ! scaling factor for eastern top slopes ISOPYC_A.79
& ath_b(imt,km) ! scaling factor for eastern bottom slopes ISOPYC_A.80
REAL rhoz(imt,km),rhotmp,rhotmp1 ISOPYC_A.81
ISOPYC_A.82
C top and bottom boundary conditions ISOPYC_A.83
REAL top_bc(km), bot_bc(km) ISOPYC_A.84
ISOPYC_A.85
C ratio of zonal to vertical density gradients ISOPYC_A.86
REAL ste(imt,km),sbe(imt,km) ISOPYC_A.87
ISOPYC_A.88
C grad^2 (ste,stn) etc for biharmonic GM scheme ISOPYC_A.89
REAL ste_d2(imt,km),sbe_d2(imt,km), ISOPYC_A.90
& stn_d2(imt,km),sbn_d2(imt,km) ISOPYC_A.91
ISOPYC_A.92
! Temporary values used in optimised tanh calculations ISOPYC_A.93
REAL tanh_temp(imt*2) ISOPYC_A.94
ISOPYC_A.95
REAL part1, part2, Ath0_j, Ath0_jp1 ISOPYC_A.96
ISOPYC_A.97
c0=0. ISOPYC_A.98
c1=1. ISOPYC_A.99
p5=.5 ISOPYC_A.100
epsln=1.0e-25 ISOPYC_A.101
c do k=1,km ISOPYC_A.102
c do i=1,imt ISOPYC_A.103
c fmpp(i,k)=0. ISOPYC_A.104
c enddo ISOPYC_A.105
c enddo ISOPYC_A.106
c ISOPYC_A.107
C coefficient for biharmonic GM scheme ISOPYC_A.108
athkdftu_bi=athkdf_bi ISOPYC_A.109
athkdftv_bi=athkdf_bi ISOPYC_A.110
ISOPYC_A.111
IF (.NOT.L_OVISBECK) THEN ISOPYC_A.112
do k=1,km ISOPYC_A.113
do i=1,imt ISOPYC_A.114
athkdftu_mom(i,k)=athkdf(k) ISOPYC_A.115
athkdftv_mom(i,k)=athkdf(k) ISOPYC_A.116
enddo ISOPYC_A.117
enddo ISOPYC_A.118
ELSE ISOPYC_A.119
do k=1,km ISOPYC_A.120
do i=1,imt ISOPYC_A.121
athkdftu_mom(i,k)=athkdftu(i,j) ISOPYC_A.122
athkdftv_mom(i,k)=athkdftv(i,j) ISOPYC_A.123
enddo ISOPYC_A.124
enddo ISOPYC_A.125
ENDIF ! L_OVISBECK ISOPYC_A.126
ISOPYC_A.127
IF (L_OISOGM) THEN ISOPYC_A.128
ISOPYC_A.129
do k=1,km ISOPYC_A.130
top_bc(k) = c1 ISOPYC_A.131
bot_bc(k) = c1 ISOPYC_A.132
enddo ISOPYC_A.133
top_bc(1) = c0 ISOPYC_A.134
bot_bc(km) = c0 ISOPYC_A.135
ISOPYC_A.136
IF (J+J_OFFSET.eq.2) THEN ISOPYC_A.137
ISOPYC_A.138
do k=1,km ISOPYC_A.139
km1 = max(k-1,1) ISOPYC_A.140
kp1 = min(k+1,km) ISOPYC_A.141
do i=1,imt ISOPYC_A.142
at = ((alphai(i,k,0) + alphai(i,k,1)) + alphai(i,km1,0)) ISOPYC_A.143
& + alphai(i,km1,1) ISOPYC_A.144
bt = ((betai(i,k,0) + betai(i,k,1)) + betai(i,km1,0)) ISOPYC_A.145
& + betai(i,km1,1) ISOPYC_A.146
stn(i,k,1) = -(at*(ddyt(i,k,1,1) + ddyt(i,km1,1,1)) ISOPYC_A.147
& + bt*(ddyt(i,k,2,1) + ddyt(i,km1,2,1))) ISOPYC_A.148
& / (at*(ddzt(i,km1,1,0) + ddzt(i,km1,1,1)) ISOPYC_A.149
& + bt*(ddzt(i,km1,2,0) + ddzt(i,km1,2,1))+epsln) ISOPYC_A.150
c ISOPYC_A.151
ab = ((alphai(i,k,0) + alphai(i,k,1)) + alphai(i,kp1,0)) ISOPYC_A.152
& + alphai(i,kp1,1) ISOPYC_A.153
bb = ((betai(i,k,0) + betai(i,k,1)) + betai(i,kp1,0)) ISOPYC_A.154
& + betai(i,kp1,1) ISOPYC_A.155
sbn(i,k,1) = -(ab*(ddyt(i,k,1,1) + ddyt(i,kp1,1,1)) ISOPYC_A.156
& + bb*(ddyt(i,k,2,1) + ddyt(i,kp1,2,1))) ISOPYC_A.157
& / (ab*(ddzt(i,k,1,0) + ddzt(i,k,1,1)) ISOPYC_A.158
& + bb*(ddzt(i,k,2,0) + ddzt(i,k,2,1))+epsln) ISOPYC_A.159
ISOPYC_A.160
stn(i,k,0)=c0 ISOPYC_A.161
sbn(i,k,0)=c0 ISOPYC_A.162
enddo ISOPYC_A.163
enddo ISOPYC_A.164
ISOPYC_A.165
IF (L_OISOTAPER) THEN ISOPYC_A.166
do k=1,km ISOPYC_A.167
kp1 = min(k+1,km) ISOPYC_A.168
do i=1,imt ISOPYC_A.169
Ath0 = c1 ISOPYC_A.170
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.171
ath_tn(i,k,0) = c0 ISOPYC_A.172
ath_bn(i,k,0) = c0 ISOPYC_A.173
absstn = abs(stn(i,k,1)) ISOPYC_A.174
abssbn = abs(sbn(i,k,1)) ISOPYC_A.175
tanh_temp(i) = (absstn-slopec)/dslope ISOPYC_A.176
tanh_temp(i+imt) = (abssbn-slopec)/dslope ISOPYC_A.177
enddo ISOPYC_A.178
ISOPYC_A.179
call fast_tanh(imt*2,tanh_temp) ISOPYC_A.180
ISOPYC_A.181
do i=1,imt ISOPYC_A.182
ath_tn(i,k,1) = (((Ath0*fm(i,k))*fmp(i,k)) ISOPYC_A.183
& *p5)*(c1-tanh_temp(i)) ISOPYC_A.184
ath_bn(i,k,1) = (((Ath0*fm(i,kp1))*fmp(i,kp1)) ISOPYC_A.185
& *p5)*(c1-tanh_temp(i+imt)) ISOPYC_A.186
enddo ISOPYC_A.187
enddo ISOPYC_A.188
ISOPYC_A.189
ELSE ISOPYC_A.190
ISOPYC_A.191
do k=1,km ISOPYC_A.192
sc = c1/(slmxr*dtxsqr(k)) ISOPYC_A.193
kp1 = min(k+1,km) ISOPYC_A.194
do i=1,imt ISOPYC_A.195
Ath0 = c1 ISOPYC_A.196
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.197
absstn = abs(stn(i,k,1)) ISOPYC_A.198
abssbn = abs(sbn(i,k,1)) ISOPYC_A.199
if (absstn .gt. sc) then ISOPYC_A.200
ath_tn(i,k,1) = ((Ath0*fm(i,k))*fmp(i,k)) ISOPYC_A.201
& *(sc/(absstn + epsln))**2 ISOPYC_A.202
else ISOPYC_A.203
ath_tn(i,k,1) = (Ath0*fm(i,k))*fmp(i,k) ISOPYC_A.204
endif ISOPYC_A.205
if (abssbn .gt. sc) then ISOPYC_A.206
ath_bn(i,k,1) = ((Ath0*fm(i,kp1))*fmp(i,kp1)) ISOPYC_A.207
& *(sc/(abssbn + epsln))**2 ISOPYC_A.208
else ISOPYC_A.209
ath_bn(i,k,1) = (Ath0*fm(i,kp1))*fmp(i,kp1) ISOPYC_A.210
endif ISOPYC_A.211
enddo ISOPYC_A.212
enddo ISOPYC_A.213
ENDIF ! taper ISOPYC_A.214
ISOPYC_A.215
ELSE ! j.ne.2 ISOPYC_A.216
ISOPYC_A.217
do k=1,km ISOPYC_A.218
do i=1,imt ISOPYC_A.219
stn(i,k,0)=stn(i,k,1) ISOPYC_A.220
ath_tn(i,k,0)=ath_tn(i,k,1) ISOPYC_A.221
sbn(i,k,0)=sbn(i,k,1) ISOPYC_A.222
ath_bn(i,k,0)=ath_bn(i,k,1) ISOPYC_A.223
enddo ISOPYC_A.224
enddo ISOPYC_A.225
ISOPYC_A.226
IF (L_OBIHARMGM) THEN ISOPYC_A.227
do k=1,km ISOPYC_A.228
do i=1,imt ISOPYC_A.229
stn(i,k,1)=stn(i,k,2) ISOPYC_A.230
stn(i,k,2)=c0 ISOPYC_A.231
ath_tn(i,k,1)=ath_tn(i,k,2) ISOPYC_A.232
ath_tn(i,k,2)=c0 ISOPYC_A.233
sbn(i,k,1)=sbn(i,k,2) ISOPYC_A.234
sbn(i,k,2)=c0 ISOPYC_A.235
ath_bn(i,k,1)=ath_bn(i,k,2) ISOPYC_A.236
ath_bn(i,k,2)=c0 ISOPYC_A.237
enddo ISOPYC_A.238
enddo ISOPYC_A.239
ENDIF ISOPYC_A.240
ISOPYC_A.241
ENDIF ! j=2 ISOPYC_A.242
ISOPYC_A.243
C update old meridional GM velocity ISOPYC_A.244
do k=1,km ISOPYC_A.245
do i=1,imt ISOPYC_A.246
adv_vntiso(i,k,0)=adv_vntiso(i,k,1) ISOPYC_A.247
enddo ISOPYC_A.248
enddo ISOPYC_A.249
c ISOPYC_A.250
c----------------------------------------------------------------------- ISOPYC_A.251
c compute the meridional component of the isopycnal mixing velocity ISOPYC_A.252
c at the center of the northern face of the "T" cells. ISOPYC_A.253
c----------------------------------------------------------------------- ISOPYC_A.254
c ISOPYC_A.255
IF (L_OBIHARMGM) THEN ISOPYC_A.256
do k=1,km ISOPYC_A.257
km1 = max(k-1,1) ISOPYC_A.258
kp1 = min(k+1,km) ISOPYC_A.259
do i=1,imt ISOPYC_A.260
at = ((alphai(i,k,1) + alphai(i,k,2)) + alphai(i,km1,1)) ISOPYC_A.261
& + alphai(i,km1,2) ISOPYC_A.262
bt = ((betai(i,k,1) + betai(i,k,2)) + betai(i,km1,1)) ISOPYC_A.263
& + betai(i,km1,2) ISOPYC_A.264
stn(i,k,2) = -(at*(ddyt(i,k,1,2) + ddyt(i,km1,1,2)) ISOPYC_A.265
& + bt*(ddyt(i,k,2,2) + ddyt(i,km1,2,2))) ISOPYC_A.266
& / (at*(ddzt(i,km1,1,1) + ddzt(i,km1,1,2)) ISOPYC_A.267
& + bt*(ddzt(i,km1,2,1) + ddzt(i,km1,2,2))+epsln) ISOPYC_A.268
ISOPYC_A.269
ab = ((alphai(i,k,1) + alphai(i,k,2)) + alphai(i,kp1,1)) ISOPYC_A.270
& + alphai(i,kp1,2) ISOPYC_A.271
bb = ((betai(i,k,1) + betai(i,k,2)) + betai(i,kp1,1)) ISOPYC_A.272
& + betai(i,kp1,2) ISOPYC_A.273
sbn(i,k,2) = -(ab*(ddyt(i,k,1,2) + ddyt(i,kp1,1,2)) ISOPYC_A.274
& + bb*(ddyt(i,k,2,2) + ddyt(i,kp1,2,2))) ISOPYC_A.275
& / (ab*(ddzt(i,k,1,1) + ddzt(i,k,1,2)) ISOPYC_A.276
& + bb*(ddzt(i,k,2,1) + ddzt(i,k,2,2))+epsln) ISOPYC_A.277
enddo ISOPYC_A.278
enddo ISOPYC_A.279
ISOPYC_A.280
ELSE ISOPYC_A.281
ISOPYC_A.282
do k=1,km ISOPYC_A.283
km1 = max(k-1,1) ISOPYC_A.284
kp1 = min(k+1,km) ISOPYC_A.285
do i=1,imt ISOPYC_A.286
at = ((alphai(i,k,0) + alphai(i,k,1))+alphai(i,km1,0)) ISOPYC_A.287
& + alphai(i,km1,1) ISOPYC_A.288
bt = ((betai(i,k,0) + betai(i,k,1)) + betai(i,km1,0)) ISOPYC_A.289
& + betai(i,km1,1) ISOPYC_A.290
stn(i,k,1) = -(at*(ddyt(i,k,1,1) + ddyt(i,km1,1,1)) ISOPYC_A.291
& + bt*(ddyt(i,k,2,1) + ddyt(i,km1,2,1))) ISOPYC_A.292
& / (at*(ddzt(i,km1,1,0) + ddzt(i,km1,1,1)) ISOPYC_A.293
& + bt*(ddzt(i,km1,2,0) + ddzt(i,km1,2,1))+epsln) ISOPYC_A.294
ISOPYC_A.295
ab = ((alphai(i,k,0) + alphai(i,k,1))+alphai(i,kp1,0)) ISOPYC_A.296
& + alphai(i,kp1,1) ISOPYC_A.297
bb = ((betai(i,k,0) + betai(i,k,1)) + betai(i,kp1,0)) ISOPYC_A.298
& + betai(i,kp1,1) ISOPYC_A.299
sbn(i,k,1) = -(ab*(ddyt(i,k,1,1) + ddyt(i,kp1,1,1)) ISOPYC_A.300
& + bb*(ddyt(i,k,2,1) + ddyt(i,kp1,2,1))) ISOPYC_A.301
& / (ab*(ddzt(i,k,1,0) + ddzt(i,k,1,1)) ISOPYC_A.302
& + bb*(ddzt(i,k,2,0) + ddzt(i,k,2,1))+epsln) ISOPYC_A.303
enddo ISOPYC_A.304
enddo ISOPYC_A.305
ENDIF ! L_OBIHARMGM ISOPYC_A.306
ISOPYC_A.307
do k=1,km ISOPYC_A.308
km1 = max(k-1,1) ISOPYC_A.309
kp1 = min(k+1,km) ISOPYC_A.310
do i=1,imt ISOPYC_A.311
at = (alphai(i,k,0) + alphai(i,km1,0)) ISOPYC_A.312
bt = (betai(i,k,0) + betai(i,km1,0)) ISOPYC_A.313
rhoz(i,k)=-(at*(ddzt(i,km1,1,0)) ISOPYC_A.314
& + bt*(ddzt(i,km1,2,0))+epsln)/4. ISOPYC_A.315
enddo ISOPYC_A.316
enddo ISOPYC_A.317
ISOPYC_A.318
IF (L_OISOTAPER) THEN ISOPYC_A.319
ISOPYC_A.320
IF (L_OBIHARMGM) THEN ISOPYC_A.321
do k=1,km ISOPYC_A.322
kp1 = min(k+1,km) ISOPYC_A.323
do i=1,imt ISOPYC_A.324
Ath0 = c1 ISOPYC_A.325
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.326
absstn = abs(stn(i,k,2)) ISOPYC_A.327
abssbn = abs(sbn(i,k,2)) ISOPYC_A.328
tanh_temp(i) = (absstn-slopec)/dslope ISOPYC_A.329
tanh_temp(i+imt) = (abssbn-slopec)/dslope ISOPYC_A.330
enddo ISOPYC_A.331
ISOPYC_A.332
call fast_tanh(imt*2,tanh_temp) ISOPYC_A.333
ISOPYC_A.334
do i=1,imt ISOPYC_A.335
ath_tn(i,k,2) = (((Ath0*fmp(i,k))*fmpp(i,k)) ISOPYC_A.336
& *p5)*(c1-tanh_temp(i)) ISOPYC_A.337
ath_bn(i,k,2) = (((Ath0*fmp(i,kp1))*fmpp(i,kp1)) ISOPYC_A.338
& *p5)*(c1-tanh_temp(i+imt)) ISOPYC_A.339
enddo ISOPYC_A.340
enddo ISOPYC_A.341
ELSE ISOPYC_A.342
do k=1,km ISOPYC_A.343
kp1 = min(k+1,km) ISOPYC_A.344
do i=1,imt ISOPYC_A.345
Ath0 = c1 ISOPYC_A.346
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.347
absstn = abs(stn(i,k,1)) ISOPYC_A.348
abssbn = abs(sbn(i,k,1)) ISOPYC_A.349
ISOPYC_A.350
tanh_temp(i) = (absstn-slopec)/dslope ISOPYC_A.351
tanh_temp(i+imt) = (abssbn-slopec)/dslope ISOPYC_A.352
enddo ISOPYC_A.353
ISOPYC_A.354
call fast_tanh(imt*2,tanh_temp) ISOPYC_A.355
ISOPYC_A.356
do i=1,imt ISOPYC_A.357
ath_tn(i,k,1) = (((Ath0*fm(i,k))*fmp(i,k)) ISOPYC_A.358
& *p5)*(c1-tanh_temp(i)) ISOPYC_A.359
ath_bn(i,k,1) = (((Ath0*fm(i,kp1))*fmp(i,kp1)) ISOPYC_A.360
& *p5)*(c1-tanh_temp(i+imt)) ISOPYC_A.361
enddo ISOPYC_A.362
enddo ISOPYC_A.363
ENDIF !L_OBIHARMGM ISOPYC_A.364
ISOPYC_A.365
ELSE ISOPYC_A.366
ISOPYC_A.367
do k=1,km ISOPYC_A.368
sc = c1/(slmxr*dtxsqr(k)) ISOPYC_A.369
kp1 = min(k+1,km) ISOPYC_A.370
do i=1,imt ISOPYC_A.371
IF (L_OBIHARMGM) THEN ISOPYC_A.372
Ath0 = c1 ISOPYC_A.373
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.374
absstn = abs(stn(i,k,2)) ISOPYC_A.375
abssbn = abs(sbn(i,k,2)) ISOPYC_A.376
if (absstn .gt. sc) then ISOPYC_A.377
ath_tn(i,k,2) = ((Ath0*fmp(i,k))*fmpp(i,k)) ISOPYC_A.378
& *(sc/(absstn + epsln))**2 ISOPYC_A.379
else ISOPYC_A.380
ath_tn(i,k,2) = (Ath0*fmp(i,k))*fmpp(i,k) ISOPYC_A.381
endif ISOPYC_A.382
if (abssbn .gt. sc) then ISOPYC_A.383
ath_bn(i,k,2) = ((Ath0*fmp(i,kp1))*fmpp(i,kp1)) ISOPYC_A.384
& *(sc/(abssbn + epsln))**2 ISOPYC_A.385
else ISOPYC_A.386
ath_bn(i,k,2) = Ath0*fmp(i,kp1)*fmpp(i,kp1) ISOPYC_A.387
endif ISOPYC_A.388
ELSE ISOPYC_A.389
Ath0 = c1 ISOPYC_A.390
c Ath0 = athkdftv_bi(i,k) ISOPYC_A.391
absstn = abs(stn(i,k,1)) ISOPYC_A.392
abssbn = abs(sbn(i,k,1)) ISOPYC_A.393
if (absstn .gt. sc) then ISOPYC_A.394
ath_tn(i,k,1) = ((Ath0*fm(i,k))*fmp(i,k)) ISOPYC_A.395
& *(sc/(absstn + epsln))**2 ISOPYC_A.396
else ISOPYC_A.397
ath_tn(i,k,1) = (Ath0*fm(i,k))*fmp(i,k) ISOPYC_A.398
endif ISOPYC_A.399
if (abssbn .gt. sc) then ISOPYC_A.400
ath_bn(i,k,1) = ((Ath0*fm(i,kp1))*fmp(i,kp1)) ISOPYC_A.401
& *(sc/(abssbn + epsln))**2 ISOPYC_A.402
else ISOPYC_A.403
ath_bn(i,k,1) = (Ath0*fm(i,kp1))*fmp(i,kp1) ISOPYC_A.404
endif ISOPYC_A.405
ENDIF ! L_OBIHARMGM ISOPYC_A.406
enddo ISOPYC_A.407
enddo ISOPYC_A.408
ENDIF ! taper ISOPYC_A.409
ISOPYC_A.410
do k=1,km ISOPYC_A.411
do i=1,imt ISOPYC_A.412
Ath0 = (athkdftv_mom(i,k)*(2.*dz2r(k)))*csu(j) ISOPYC_A.413
adv_vntiso(i,k,1) = -((ath_tn(i,k,1)*stn(i,k,1))*top_bc(k)- ISOPYC_A.414
& (ath_bn(i,k,1)*sbn(i,k,1))*bot_bc(k))* ISOPYC_A.415
& Ath0 ISOPYC_A.416
enddo ISOPYC_A.417
enddo ISOPYC_A.418
ISOPYC_A.419
IF (L_OBIHARMGM) THEN ISOPYC_A.420
ISOPYC_A.421
do k=1,km ISOPYC_A.422
do i=1,imt ISOPYC_A.423
Ath0_j = (athkdftv_bi*dytr(j))*cst(j) ISOPYC_A.424
Ath0_jp1 = (athkdftv_bi*dytr(j+1))*cst(j+1) ISOPYC_A.425
ISOPYC_A.426
part1 = (ath_tn(i,k,2)*stn(i,k,2)) - ISOPYC_A.427
& (ath_tn(i,k,1)*stn(i,k,1))*Ath0_jp1 ISOPYC_A.428
part2 = (ath_tn(i,k,1)*stn(i,k,1)) - ISOPYC_A.429
& (ath_tn(i,k,0)*stn(i,k,0))*Ath0_j ISOPYC_A.430
stn_d2(i,k) = ((((part2-part1) ISOPYC_A.431
& *dyur(j))*csr(j))*fm(i,k))*fmp(i,k) ISOPYC_A.432
ISOPYC_A.433
ISOPYC_A.434
ISOPYC_A.435
part1 = (ath_bn(i,k,2)*sbn(i,k,2) ISOPYC_A.436
& - ath_bn(i,k,1)*sbn(i,k,1))*Ath0_jp1 ISOPYC_A.437
part2 = (ath_bn(i,k,1)*sbn(i,k,1) ISOPYC_A.438
& - ath_bn(i,k,0)*sbn(i,k,0))*Ath0_j ISOPYC_A.439
sbn_d2(i,k) = ((((part2-part1) ISOPYC_A.440
& *dyur(j))*csr(j))*fm(i,k))*fmp(i,k) ISOPYC_A.441
enddo ISOPYC_A.442
enddo ISOPYC_A.443
ISOPYC_A.444
do k=1,km ISOPYC_A.445
do i=1,imt ISOPYC_A.446
adv_vntiso(i,k,1) = adv_vntiso(i,k,1) - ISOPYC_A.447
& (((stn_d2(i,k)*top_bc(k))-(sbn_d2(i,k)*bot_bc(k)))* ISOPYC_A.448
& (2.*dz2r(k)))*csu(j) ISOPYC_A.449
enddo ISOPYC_A.450
enddo ISOPYC_A.451
ISOPYC_A.452
ENDIF ! L_OBIHARMGM ISOPYC_A.453
c ISOPYC_A.454
c----------------------------------------------------------------------- ISOPYC_A.455
c compute the zonal component of the isopycnal mixing velocity ISOPYC_A.456
c at the center of the eastern face of the "T" grid cell. ISOPYC_A.457
c----------------------------------------------------------------------- ISOPYC_A.458
c ISOPYC_A.459
do k=1,km ISOPYC_A.460
km1 = max(k-1,1) ISOPYC_A.461
kp1 = min(k+1,km) ISOPYC_A.462
do i=1,imt-1 ISOPYC_A.463
c ISOPYC_A.464
at = (alphai(i,k,0) + alphai(i+1,k,0) + alphai(i,km1,0) ISOPYC_A.465
& + alphai(i+1,km1,0)) ISOPYC_A.466
bt = (betai(i,k,0) + betai(i+1,k,0) + betai(i,km1,0) ISOPYC_A.467
& + betai(i+1,km1,0)) ISOPYC_A.468
ste(i,k) =-(at*(ddxt(i,k,1,0) + ddxt(i,km1,1,0)) ISOPYC_A.469
& + bt*(ddxt(i,k,2,0) + ddxt(i,km1,2,0))) ISOPYC_A.470
& / (at*(ddzt(i,km1,1,0) + ddzt(i+1,km1,1,0)) ISOPYC_A.471
& + bt*(ddzt(i,km1,2,0) + ddzt(i+1,km1,2,0))+epsln) ISOPYC_A.472
c ISOPYC_A.473
ab = (alphai(i,k,0) + alphai(i+1,k,0) + alphai(i,kp1,0) ISOPYC_A.474
& + alphai(i+1,kp1,0)) ISOPYC_A.475
bb = (betai(i,k,0) + betai(i+1,k,0) + betai(i,kp1,0) ISOPYC_A.476
& + betai(i+1,kp1,0)) ISOPYC_A.477
sbe(i,k) =-(ab*(ddxt(i,k,1,0) + ddxt(i,kp1,1,0)) ISOPYC_A.478
& + bb*(ddxt(i,k,2,0) + ddxt(i,kp1,2,0))) ISOPYC_A.479
& / (ab*(ddzt(i,k,1,0) + ddzt(i+1,k,1,0)) ISOPYC_A.480
& + bb*(ddzt(i,k,2,0) + ddzt(i+1,k,2,0))+epsln) ISOPYC_A.481
ISOPYC_A.482
enddo ISOPYC_A.483
enddo ISOPYC_A.484
ISOPYC_A.485
call SETBCX(ste(1,1),imt,km) ISOPYC_A.486
call SETBCX(sbe(1,1),imt,km) ISOPYC_A.487
ISOPYC_A.488
IF (L_OISOTAPER) THEN ISOPYC_A.489
do k=1,km ISOPYC_A.490
kp1 = min(k+1,km) ISOPYC_A.491
do i=2,imt-1 ISOPYC_A.492
Ath0 = c1 ISOPYC_A.493
c Ath0 = athkdftu_bi(i,k) ISOPYC_A.494
absste = abs(ste(i,k)) ISOPYC_A.495
abssbe = abs(sbe(i,k)) ISOPYC_A.496
tanh_temp(i) = (absste-slopec)/dslope ISOPYC_A.497
tanh_temp(i+imt) = (abssbe-slopec)/dslope ISOPYC_A.498
enddo ISOPYC_A.499
ISOPYC_A.500
call fast_tanh(imt-2,tanh_temp(2)) ISOPYC_A.501
call fast_tanh(imt-2,tanh_temp(2+imt)) ISOPYC_A.502
ISOPYC_A.503
do i= 2, imt-1 ISOPYC_A.504
ath_t(i,k) = Ath0*fm(i,k)*fm(i+1,k) ISOPYC_A.505
& *p5*(c1-tanh_temp(i)) ISOPYC_A.506
ath_b(i,k) = Ath0*fm(i,kp1)*fm(i+1,kp1) ISOPYC_A.507
& *p5*(c1-tanh_temp(i+imt)) ISOPYC_A.508
enddo ISOPYC_A.509
enddo ISOPYC_A.510
ISOPYC_A.511
ELSE ISOPYC_A.512
ISOPYC_A.513
do k=1,km ISOPYC_A.514
sc = c1/(slmxr*dtxsqr(k)) ISOPYC_A.515
kp1 = min(k+1,km) ISOPYC_A.516
do i=2,imt-1 ISOPYC_A.517
Ath0 = c1 ISOPYC_A.518
c Ath0 = athkdftu_bi(i,k) ISOPYC_A.519
absste = abs(ste(i,k)) ISOPYC_A.520
abssbe = abs(sbe(i,k)) ISOPYC_A.521
if (absste .gt. sc) then ISOPYC_A.522
ath_t(i,k) = Ath0*fm(i,k)*fm(i+1,k) ISOPYC_A.523
& *(sc/(absste + epsln))**2 ISOPYC_A.524
else ISOPYC_A.525
ath_t(i,k) = Ath0*fm(i,k)*fm(i+1,k) ISOPYC_A.526
endif ISOPYC_A.527
if (abssbe .gt. sc) then ISOPYC_A.528
ath_b(i,k) = Ath0*fm(i,kp1)*fm(i+1,kp1) ISOPYC_A.529
& *(sc/(abssbe + epsln))**2 ISOPYC_A.530
else ISOPYC_A.531
ath_b(i,k) = Ath0*fm(i,kp1)*fm(i+1,kp1) ISOPYC_A.532
endif ISOPYC_A.533
enddo ISOPYC_A.534
enddo ISOPYC_A.535
ENDIF ! taper ISOPYC_A.536
ISOPYC_A.537
call SETBCX(ath_t(1,1),imt,km) ISOPYC_A.538
call SETBCX(ath_b(1,1),imt,km) ISOPYC_A.539
ISOPYC_A.540
do k=1,km ISOPYC_A.541
do i=2,imt-1 ISOPYC_A.542
rhotmp=(((ste(i,k)+ste(i-1,k))/2.)**2.+ ISOPYC_A.543
& ((stn(i,k,1)+stn(i,k,0))/2.)**2.) ISOPYC_A.544
rhotmp1=sqrt(rhotmp) ISOPYC_A.545
if (rhoz(i,k).lt.0.) rhoz(i,k)=0. ISOPYC_A.546
srho(i,k)=fm(i,k)*rhotmp1*p5*sqrt(rhoz(i,k))* ISOPYC_A.547
& (c1-tanh((rhotmp1-slopec)/dslope)) ISOPYC_A.548
enddo ISOPYC_A.549
enddo ISOPYC_A.550
call SETBCX(srho(1,1),imt,km) ISOPYC_A.551
ISOPYC_A.552
do k=1,km ISOPYC_A.553
do i=2,imt-1 ISOPYC_A.554
Ath0 = athkdftu_mom(i,k) ISOPYC_A.555
adv_vetiso(i,k) = -((ath_t(i,k)*ste(i,k))*Ath0*top_bc(k) - ISOPYC_A.556
& (ath_b(i,k)*sbe(i,k))*Ath0*bot_bc(k))*(2.*dz2r(k)) ISOPYC_A.557
enddo ISOPYC_A.558
enddo ISOPYC_A.559
ISOPYC_A.560
IF (L_OBIHARMGM) THEN ISOPYC_A.561
ISOPYC_A.562
do k=1,km ISOPYC_A.563
do i=2,imt-1 ISOPYC_A.564
Ath0=athkdftu_bi ISOPYC_A.565
ste_d2(i,k)=(-((ath_t(i+1,k)*ste(i+1,k) ISOPYC_A.566
& - ath_t(i,k)*ste(i,k))*Ath0*dxtr(i+1) ISOPYC_A.567
& - (ath_t(i,k)*ste(i,k) - ath_t(i-1,k)*ste(i-1,k))*Ath0* ISOPYC_A.568
& dxtr(i))*dxur(i))*fm(i,k)*fm(i+1,k) ISOPYC_A.569
ISOPYC_A.570
sbe_d2(i,k)=(-((ath_b(i+1,k)*sbe(i+1,k) ISOPYC_A.571
& - ath_b(i,k)*sbe(i,k))*Ath0*dxtr(i+1) ISOPYC_A.572
& - (ath_b(i,k)*sbe(i,k) - ath_b(i-1,k)*sbe(i-1,k))*Ath0* ISOPYC_A.573
& dxtr(i))*dxur(i))*fm(i,k)*fm(i+1,k) ISOPYC_A.574
enddo ISOPYC_A.575
enddo ISOPYC_A.576
ISOPYC_A.577
call SETBCX(ste_d2(1,1),imt,km) ISOPYC_A.578
call SETBCX(sbe_d2(1,1),imt,km) ISOPYC_A.579
ISOPYC_A.580
do k=1,km ISOPYC_A.581
do i=2,imt-1 ISOPYC_A.582
c adv_vetiso(i,k) = -(ste_d2(i,k)*top_bc(k) - ISOPYC_A.583
adv_vetiso(i,k) = adv_vetiso(i,k) -(ste_d2(i,k)*top_bc(k) - ISOPYC_A.584
& sbe_d2(i,k)*bot_bc(k))*(2.*dz2r(k)) ISOPYC_A.585
enddo ISOPYC_A.586
enddo ISOPYC_A.587
ISOPYC_A.588
ENDIF !L_OBIHARMGM) ISOPYC_A.589
c ISOPYC_A.590
c set the boundary conditions ISOPYC_A.591
c ISOPYC_A.592
call SETBCX (adv_vetiso(1,1), imt, km) ISOPYC_A.593
ISOPYC_A.594
c ISOPYC_A.595
c----------------------------------------------------------------------- ISOPYC_A.596
c compute the vertical component of the isopycnal mixing velocity ISOPYC_A.597
c at the center of the bottom face of the "T" cells, using the ISOPYC_A.598
c continuity equation for the isopycnal mixing velocities ISOPYC_A.599
c----------------------------------------------------------------------- ISOPYC_A.600
c ISOPYC_A.601
do i=1,imt ISOPYC_A.602
adv_vbtiso(i,0) = c0 ISOPYC_A.603
enddo ISOPYC_A.604
c ISOPYC_A.605
do k=1,km-1 ISOPYC_A.606
do i=2,imt ISOPYC_A.607
adv_vbtiso(i,k) = dzt(k)*cstr(j)*( ISOPYC_A.608
& (adv_vetiso(i,k) - adv_vetiso(i-1,k))*dxtr(i) + ISOPYC_A.609
& (adv_vntiso(i,k,1) - adv_vntiso(i,k,0))*dytr(j)) ISOPYC_A.610
enddo ISOPYC_A.611
enddo ISOPYC_A.612
c ISOPYC_A.613
do k=1,km-1 ISOPYC_A.614
do i=2,imt ISOPYC_A.615
adv_vbtiso(i,k) = adv_vbtiso(i,k) + adv_vbtiso(i,k-1) ISOPYC_A.616
enddo ISOPYC_A.617
enddo ISOPYC_A.618
c ISOPYC_A.619
do i=2,imt ISOPYC_A.620
adv_vbtiso(i,kmt(i)) = c0 ISOPYC_A.621
enddo ISOPYC_A.622
c ISOPYC_A.623
c set the boundary conditions ISOPYC_A.624
c ISOPYC_A.625
call SETBCX (adv_vbtiso(1,0), imt, km+1) ISOPYC_A.626
ISOPYC_A.627
ENDIF ! GM ISOPYC_A.628
ISOPYC_A.629
RETURN ISOPYC_A.630
END ISOPYC_A.631
*ENDIF ISOPYC_A.632