SUBROUTINE ISOFLUX( 2,6ISOFLUX.43
*CALL ARGSIZE 
ISOFLUX.44
*CALL COCAWRKA ISOFLUX.45
& ,m,j,cstr,dyur,dxur,dzt2r,csu,csjm, ISOFLUX.46
& dytr,dxtr,cst,dxt4r,dyt4r,Ai_ez,Ai_nz,Ai_bx,Ai_by,K11,K22,K33, ISOFLUX.47
& itt,adv_vbtiso,adv_fbiso,tempx,tempa,tempb,j_1) ISOFLUX.48
ISOFLUX.49
IMPLICIT NONE ISOFLUX.50
ISOFLUX.51
*CALL OARRYSIZ ISOFLUX.52
*CALL TYPSIZE ISOFLUX.53
*CALL CNTLOCN ISOFLUX.54
*CALL COCTWRKA ISOFLUX.55
*CALL OTIMER ISOFLUX.56
c ISOFLUX.57
c======================================================================= ISOFLUX.58
c isopycnal diffusive tracer fluxes are computed. ISOFLUX.59
c======================================================================= ISOFLUX.60
c ISOFLUX.61
C Input variables ISOFLUX.62
INTEGER j,m,itt,j_1 ISOFLUX.63
REAL dzt2r(km),csu(jmt),dytr(jmt),cst(jmt),dxur(imt), ISOFLUX.64
& cstr(jmt),dxtr(imt),dyur(imt),dxt4r(imt),dyt4r(jmt), ISOFLUX.65
& csjm ISOFLUX.66
ISOFLUX.67
c compute advective tracer flux at the center of the bottom face of ISOFLUX.68
c the "T" cells ISOFLUX.69
c----------------------------------------------------------------------- ISOFLUX.70
C GM I/O variables ISOFLUX.71
REAL adv_fbiso(imt_gmm,0:km_gmm), ! vert advective tracer flux ISOFLUX.72
C ! due to GM ISOFLUX.73
& adv_vbtiso(imt_gmm,0:km_gmm) ! vertical GM velocity ISOFLUX.74
ISOFLUX.75
C Diffusion coefficient variables from AI_CALC ISOFLUX.76
REAL Ai_ez(imt_iso,km_iso,0:1,0:1),Ai_bx(imt_iso,km_iso,0:1,0:1), ISOFLUX.77
& Ai_by(imt_iso,km_iso,0:1,0:1),Ai_nz(imt_iso,km_iso,0:1,0:1), ISOFLUX.78
& K11(imt_iso,km_iso),K22(imt_iso,km_iso),K33(imt_iso,km_iso) ISOFLUX.79
ISOFLUX.80
C Heating rate variables ISOFLUX.81
REAL tempx(imt,km),tempa(imt,km),tempb(imt,km) ISOFLUX.82
ISOFLUX.83
C Local variables ISOFLUX.84
REAL c0,c1,p5,p25,fxa,fxb ISOFLUX.85
ISOFLUX.86
REAL csjm_loc ISOFLUX.87
INTEGER i,k,kr,ip,km1kr,kpkr,jq,km1kr0,kpkr0,km1kr1,kpkr1 ISOFLUX.88
REAL dzt4r,epsln,flux_x,sumz,Ai0,sumy,facty, ISOFLUX.89
& flux_y,csu_dzt4r,sumx ISOFLUX.90
ISOFLUX.91
REAL diff_tx(imt,km), ! tracer change due to zonal flux div ISOFLUX.92
& diff_ty(imt,km), ! tracer change due to merid flux div ISOFLUX.93
& diff_tz(imt,km) ! tracer change due to vert flux div ISOFLUX.94
ISOFLUX.95
REAL diff_fe(imt,km), ! zonal tracer flux ISOFLUX.96
& diff_fbiso(imt,0:km) ! vertical tracer flux ISOFLUX.97
ISOFLUX.98
C statement functions ISOFLUX.99
REAL drodxe,drodze,drodyn,drodzn,drodxb,drodyb,drodzb, ISOFLUX.100
& drodye,drodxn ISOFLUX.101
c ISOFLUX.102
c statement functions ISOFLUX.103
c ISOFLUX.104
c ISOFLUX.105
drodxe(i,k,ip) = alphai(i+ip,k,0)*ddxt(i,k,1,0) + ISOFLUX.106
& betai(i+ip,k,0)*ddxt(i,k,2,0) ISOFLUX.107
drodze(i,k,ip,kr) = alphai(i+ip,k,0)*ddzt(i+ip,k-1+kr,1,0) + ISOFLUX.108
& betai(i+ip,k,0)*ddzt(i+ip,k-1+kr,2,0) ISOFLUX.109
c ISOFLUX.110
drodyn(i,k,jq) = alphai(i,k,jq)*ddyt(i,k,1,1) + ISOFLUX.111
& betai(i,k,jq)*ddyt(i,k,2,1) ISOFLUX.112
drodzn(i,k,jq,kr) = alphai(i,k,jq)*ddzt(i,k-1+kr,1,jq) + ISOFLUX.113
& betai(i,k,jq)*ddzt(i,k-1+kr,2,jq) ISOFLUX.114
c ISOFLUX.115
drodxb(i,k,ip,kr) = alphai(i,k+kr,0)*ddxt(i-1+ip,k+kr,1,0) + ISOFLUX.116
& betai(i,k+kr,0)*ddxt(i-1+ip,k+kr,2,0) ISOFLUX.117
drodyb(i,k,jq,kr) = alphai(i,k+kr,0)*ddyt(i,k+kr,1,jq) + ISOFLUX.118
& betai(i,k+kr,0)*ddyt(i,k+kr,2,jq) ISOFLUX.119
drodzb(i,k,kr) = alphai(i,k+kr,0)*ddzt(i,k,1,0) + ISOFLUX.120
& betai(i,k+kr,0)*ddzt(i,k,2,0) ISOFLUX.121
c ISOFLUX.122
c statement functions only needed with full tensor ISOFLUX.123
c ISOFLUX.124
drodye(i,k,ip,jq) = alphai(i+ip,k,jq)*ddyt(i+ip,k,1,jq) + ISOFLUX.125
& betai(i+ip,k,jq)*ddyt(i+ip,k,2,jq) ISOFLUX.126
drodxn(i,k,ip,jq) = alphai(i,k,jq)*ddxt(i-1+ip,k,1,jq) + ISOFLUX.127
& betai(i,k,jq)*ddxt(i-1+ip,k,2,jq) ISOFLUX.128
c ISOFLUX.129
c----------------------------------------------------------------------- ISOFLUX.130
c construct total isopycnal tracer flux at east face of "T" cells ISOFLUX.131
c----------------------------------------------------------------------- ISOFLUX.132
c ISOFLUX.133
c0=0. ISOFLUX.134
c1=1. ISOFLUX.135
p5=0.5 ISOFLUX.136
p25=0.25 ISOFLUX.137
epsln=1.0e-25 ISOFLUX.138
ISOFLUX.139
C if we are on the initialisation row of a block, need to use ISOFLUX.140
C the value of cs and dyur sent from the block below ISOFLUX.141
if (j.lt.j_1) then ISOFLUX.142
csjm_loc=csjm ISOFLUX.143
else ISOFLUX.144
csjm_loc=csu(j-1) ISOFLUX.145
endif ISOFLUX.146
ISOFLUX.147
do k=1,km ISOFLUX.148
dzt4r = p5*dzt2r(k) ISOFLUX.149
km1kr0 = max(k-1,1) ISOFLUX.150
kpkr0 = min(k,km) ISOFLUX.151
km1kr1 = max(k,1) ISOFLUX.152
kpkr1 = min(k+1,km) ISOFLUX.153
do i=2,imt-1 ISOFLUX.154
sumz = - Ai_ez(i,k,0,0) ISOFLUX.155
& *(tb(i,km1kr0,m)-tb(i,kpkr0,m)) ISOFLUX.156
& *drodxe(i,k,0)/(drodze(i,k,0,0) + epsln) ISOFLUX.157
& - Ai_ez(i,k,1,0) ISOFLUX.158
& *(tb(i+1,km1kr0,m)-tb(i+1,kpkr0,m)) ISOFLUX.159
& *drodxe(i,k,1)/(drodze(i,k,1,0) + epsln) ISOFLUX.160
& - Ai_ez(i,k,0,1) ISOFLUX.161
& *(tb(i,km1kr1,m)-tb(i,kpkr1,m)) ISOFLUX.162
& *drodxe(i,k,0)/(drodze(i,k,0,1) + epsln) ISOFLUX.163
& - Ai_ez(i,k,1,1) ISOFLUX.164
& *(tb(i+1,km1kr1,m)-tb(i+1,kpkr1,m)) ISOFLUX.165
& *drodxe(i,k,1)/(drodze(i,k,1,1) + epsln) ISOFLUX.166
flux_x = dzt4r*sumz ISOFLUX.167
diff_fe(i,k) = K11(i,k)*cstr(j)*dxur(i)* ISOFLUX.168
& (tb(i+1,k,m) - tb(i,k,m)) ISOFLUX.169
& + flux_x ISOFLUX.170
enddo ISOFLUX.171
enddo ISOFLUX.172
call SETBCX (diff_fe(1,1), imt, km) ISOFLUX.173
c ISOFLUX.174
c--------------------------------------------------------------------- ISOFLUX.175
c construct total isopycnal tracer flux at north face of "T" cells ISOFLUX.176
c--------------------------------------------------------------------- ISOFLUX.177
c ISOFLUX.178
do k=1,km ISOFLUX.179
csu_dzt4r = csu(j)*p5*dzt2r(k) ISOFLUX.180
km1kr0 = max(k-1,1) ISOFLUX.181
kpkr0 = min(k,km) ISOFLUX.182
km1kr1 = max(k,1) ISOFLUX.183
kpkr1 = min(k+1,km) ISOFLUX.184
do i=2,imt-1 ISOFLUX.185
sumz = - Ai_nz(i,k,0,0) ISOFLUX.186
& *(tb(i,km1kr0,m)-tb(i,kpkr0,m)) ISOFLUX.187
& *drodyn(i,k,0)/(drodzn(i,k,0,0) + epsln) ISOFLUX.188
& - Ai_nz(i,k,1,0) ISOFLUX.189
& *(tbp(i,km1kr0,m)-tbp(i,kpkr0,m)) ISOFLUX.190
& *drodyn(i,k,1)/(drodzn(i,k,1,0) + epsln) ISOFLUX.191
& - Ai_nz(i,k,0,1) ISOFLUX.192
& *(tb(i,km1kr1,m)-tb(i,kpkr1,m)) ISOFLUX.193
& *drodyn(i,k,0)/(drodzn(i,k,0,1) + epsln) ISOFLUX.194
& - Ai_nz(i,k,1,1) ISOFLUX.195
& *(tbp(i,km1kr1,m)-tbp(i,kpkr1,m)) ISOFLUX.196
& *drodyn(i,k,1)/(drodzn(i,k,1,1) + epsln) ISOFLUX.197
flux_y = csu_dzt4r*sumz ISOFLUX.198
ISOFLUX.199
diff_fn(i,k,m,1) = K22(i,k)*csu(j)*dyur(j)* ISOFLUX.200
& (tbp(i,k,m)-tb(i,k,m)) ISOFLUX.201
& + flux_y ISOFLUX.202
enddo ISOFLUX.203
enddo ISOFLUX.204
call SETBCX (diff_fn(1,1,m,1), imt, km) ISOFLUX.205
c ISOFLUX.206
c----------------------------------------------------------------------- ISOFLUX.207
c compute the vertical tracer flux "diff_fbiso" containing the K31 ISOFLUX.208
c and K32 components which are to be solved explicitly. The K33 ISOFLUX.209
c component will be treated implicitly. Note that there are some ISOFLUX.210
c cancellations of dxu(i-1+(0:1)) and dyu(jrow-1+(0:1)) ISOFLUX.211
c----------------------------------------------------------------------- ISOFLUX.212
c ISOFLUX.213
do k=1,km-1 ISOFLUX.214
do i=2,imt-1 ISOFLUX.215
sumx = - Ai_bx(i,k,0,0)*cstr(j)* ISOFLUX.216
& (tb(i,k,m) - tb(i-1,k,m)) ISOFLUX.217
& *drodxb(i,k,0,0)/(drodzb(i,k,0) + epsln) ISOFLUX.218
& - Ai_bx(i,k,1,0)*cstr(j)* ISOFLUX.219
& (tb(i+1,k,m) - tb(i,k,m)) ISOFLUX.220
& *drodxb(i,k,1,0)/(drodzb(i,k,0) + epsln) ISOFLUX.221
& - Ai_bx(i,k,0,1)*cstr(j)* ISOFLUX.222
& (tb(i,k+1,m) - tb(i-1,k+1,m)) ISOFLUX.223
& *drodxb(i,k,0,1)/(drodzb(i,k,1) + epsln) ISOFLUX.224
& - Ai_bx(i,k,1,1)*cstr(j)* ISOFLUX.225
& (tb(i+1,k+1,m) - tb(i,k+1,m)) ISOFLUX.226
& *drodxb(i,k,1,1)/(drodzb(i,k,1) + epsln) ISOFLUX.227
ISOFLUX.228
c sumy = - Ai_by(i,k,0,0)*csu(j-1)* ISOFLUX.229
sumy = - Ai_by(i,k,0,0)*csjm_loc* ISOFLUX.230
& (tb(i,k,m)-tbm(i,k,m)) ISOFLUX.231
& *drodyb(i,k,0,0)/(drodzb(i,k,0) + epsln) ISOFLUX.232
& - Ai_by(i,k,1,0)*csu(j)* ISOFLUX.233
& (tbp(i,k,m)-tb(i,k,m)) ISOFLUX.234
& *drodyb(i,k,1,0)/(drodzb(i,k,0) + epsln) ISOFLUX.235
& - Ai_by(i,k,0,1)*csjm_loc* ISOFLUX.236
c & - Ai_by(i,k,0,1)*csu(j-1)* ISOFLUX.237
& (tb(i,k+1,m)-tbm(i,k+1,m)) ISOFLUX.238
& *drodyb(i,k,0,1)/(drodzb(i,k,1) + epsln) ISOFLUX.239
& - Ai_by(i,k,1,1)*csu(j)* ISOFLUX.240
& (tbp(i,k+1,m)-tb(i,k+1,m)) ISOFLUX.241
& *drodyb(i,k,1,1)/(drodzb(i,k,1) + epsln) ISOFLUX.242
ISOFLUX.243
diff_fbiso(i,k) = dxt4r(i)*sumx ISOFLUX.244
& + dyt4r(j)*cstr(j)*sumy ISOFLUX.245
enddo ISOFLUX.246
enddo ISOFLUX.247
ISOFLUX.248
C set top and bottom boundary conditions on vertical isopycnal fluxes ISOFLUX.249
do i=2,imt-1 ISOFLUX.250
diff_fbiso(i,0) = c0 ISOFLUX.251
diff_fbiso(i,km) = c0 ISOFLUX.252
enddo ISOFLUX.253
call SETBCX (diff_fbiso(1,0), imt, km+1) ISOFLUX.254
ISOFLUX.255
C calculate the derivatives of the fluxes to incorporate into the ISOFLUX.256
C updated tracers ISOFLUX.257
ISOFLUX.258
do k=1,km ISOFLUX.259
do i=2,imt-1 ISOFLUX.260
diff_tx(i,k)=(diff_fe(i, k)*fm(i+1,k) ISOFLUX.261
& - diff_fe(i-1,k)*fm(i-1,k))*cstr(j)*dxtr(i) ISOFLUX.262
diff_ty(i,k)=(diff_fn(i,k,m,1 )*fmp(i,k) ISOFLUX.263
& - diff_fn(i,k,m,0)*fmm(i,k))*cstr(j)*dytr(j) ISOFLUX.264
diff_tz(i,k) = (diff_fbiso(i,k-1) - diff_fbiso(i,k)) ISOFLUX.265
& *2.*dzt2r(k) ISOFLUX.266
enddo ISOFLUX.267
enddo ISOFLUX.268
ISOFLUX.269
call SETBCX (diff_tx(1,1), imt, km) ISOFLUX.270
call SETBCX (diff_ty(1,1), imt, km) ISOFLUX.271
call SETBCX (diff_tz(1,1), imt, km) ISOFLUX.272
ISOFLUX.273
C calculate the heating rate diagnostics ISOFLUX.274
do k=1,km ISOFLUX.275
do i=1,imt ISOFLUX.276
tempx(i,k)=diff_tx(i,k) ISOFLUX.277
tempa(i,k)=diff_ty(i,k) ISOFLUX.278
tempb(i,k)=diff_tz(i,k) ISOFLUX.279
enddo ISOFLUX.280
enddo ISOFLUX.281
ISOFLUX.282
C update tracers with isopycnal fluxes ISOFLUX.283
ISOFLUX.284
do k=1,km ISOFLUX.285
do i=1,imt ISOFLUX.286
ta(i,k,m)=ta(i,k,m)+(diff_tx(i,k)+diff_ty(i,k)+ ISOFLUX.287
& diff_tz(i,k))*fm(i,k) ISOFLUX.288
enddo ISOFLUX.289
enddo ISOFLUX.290
c ISOFLUX.291
IF (L_OISOGM) THEN ISOFLUX.292
c ISOFLUX.293
c----------------------------------------------------------------------- ISOFLUX.294
c compute advective tracer flux at the center of the bottom face of ISOFLUX.295
c the "T" cells ISOFLUX.296
c----------------------------------------------------------------------- ISOFLUX.297
c ISOFLUX.298
do k=1,km-1 ISOFLUX.299
do i=1,imt ISOFLUX.300
adv_fbiso(i,k) = adv_vbtiso(i,k)* ISOFLUX.301
& (tb(i,k,m) + tb(i,k+1,m)) ISOFLUX.302
enddo ISOFLUX.303
enddo ISOFLUX.304
c ISOFLUX.305
c now consider the top and bottom boundaries ISOFLUX.306
c ISOFLUX.307
do i=1,imt ISOFLUX.308
adv_fbiso(i,0) = c0 ISOFLUX.309
adv_fbiso(i,km) = c0 ISOFLUX.310
enddo ISOFLUX.311
ENDIF ! GM SCHEME ISOFLUX.312
c ISOFLUX.313
C update southern cpt of isopycnal flux ISOFLUX.314
do k=1,km ISOFLUX.315
do i=1,imt ISOFLUX.316
diff_fn(i,k,m,0)=diff_fn(i,k,m,1) ISOFLUX.317
enddo ISOFLUX.318
enddo ISOFLUX.319
ISOFLUX.320
RETURN ISOFLUX.321
END ISOFLUX.322
*ENDIF ISOFLUX.323