SUBROUTINE IPDGMVEL 1IPDGMVEL.63
& ( J,JMT,IMT,IMTM1,KM,KMM1,KMP1,KMT,KMTP, IPDGMVEL.64
& DXTR,DYTR,DZ,DZ2R,DZZ2R, IPDGMVEL.65
& CSTR,CS,ITT,FM,FMP, IPDGMVEL.66
& J_OFFSET, ORH3F402.395
& ATHKDF,ATHKDFTU,ATHKDFTV,AHI, OLA2F403.362
& fk1,fk2,fk3, IPDGMVEL.68
& uisop,visopn,visops,wisop IPDGMVEL.69
& ,mld,zdz IPDGMVEL.70
& ) IPDGMVEL.71
IPDGMVEL.72
C IPDGMVEL.73
C IPDGMVEL.74
IMPLICIT NONE IPDGMVEL.75
C IPDGMVEL.76
*CALL OARRYSIZ 
PXORDER.25
INTEGER IPDGMVEL.77
& I,J,K, IPDGMVEL.78
& JMT,IMT,IMTM1,KM,KMM1,KMP1, IPDGMVEL.79
& J_OFFSET, ! Row offset to calc global J value in MPP code ORH3F402.396
& ITT IPDGMVEL.80
INTEGER IPDGMVEL.81
& KMT(IMT),KMTP(IMT) IPDGMVEL.82
C IPDGMVEL.83
REAL IPDGMVEL.84
& DXTR(IMT),DYTR(JMT),DZ2R(KM),DZZ2R(KMP1), IPDGMVEL.85
& DZ(KM), IPDGMVEL.86
& FM(IMT,KM),FMP(IMT,KM),CS(JMT),CSTR(JMT), IPDGMVEL.87
& AHI(KM), IPDGMVEL.88
& athkdf(KM), ! IN isopycnal thickness diffusivity IPDGMVEL.89
& ATHKDFTU(imt_vis,jmt_vis), ! IN isopycnal thickness diffusivity OLA2F403.363
C for Visbeck scheme (u* pts) OLA2F403.364
& ATHKDFTV(imt_vis,jmt_vis), ! IN isopycnal thickness diffusivity OLA2F403.365
C for Visbeck scheme (v* pts) OLA2F403.366
& fk1(IMT,KM,3), ! IN } Rows 1,2,3 of diffusion matrix IPDGMVEL.90
& fk2(IMT,KM,3), ! IN } fk1 is on EAST face, fk2 on NORTH face IPDGMVEL.91
& fk3(IMT,KM,3), ! IN } & fk3 on TOP face of grid box (i,j,k) IPDGMVEL.92
& mld(IMT), IPDGMVEL.93
& zdz(km) IPDGMVEL.94
C Define eddy induced velocity components IPDGMVEL.95
REAL IPDGMVEL.96
& UISOP(IMT,KM), ! OUT u* zonal velocity at E face of T gridbox IPDGMVEL.97
& VISOPN(IMT,KM), ! IN/OUT v* meridional velocity at N face IPDGMVEL.98
& VISOPS(IMT,KM), ! OUT v* meridional velocity at S face IPDGMVEL.99
& WISOP(IMT,KMP1) ! OUT w* vertical velocity at TOP of gridbox IPDGMVEL.100
C IPDGMVEL.101
C Declare local constants and arrays IPDGMVEL.102
C IPDGMVEL.103
REAL IPDGMVEL.104
& fxa,fxb IPDGMVEL.105
*CALL CNTLOCN IPDGMVEL.106
C IPDGMVEL.108
*IF DEF,OISOPYCGM IPDGMVEL.109
C* IPDGMVEL.110
C***************************************************** IPDGMVEL.111
C In our model, unlike MOM, the fk(i,k,m) values are multiplied IPDGMVEL.112
C by the along isopycnal diffusion coeff in IPDCOFCL. Hence, wherever IPDGMVEL.113
C we use the fk values in this subroutine, they need to be divided IPDGMVEL.114
C by ahi(k), so that we are just dealing with the density gradients. IPDGMVEL.115
C A better long-term option may be to change IPDCOFCL etc, but this is IPDGMVEL.116
C easier at the moment. IPDGMVEL.117
c***************************************************** IPDGMVEL.118
IPDGMVEL.119
C IPDGMVEL.120
C visopn is initialised in BLOKINIT except when j=2. IPDGMVEL.121
C Initialise arrays here for case J=2 IPDGMVEL.122
if (j+J_OFFSET .eq. 2) then ORH3F402.397
do k=1,km IPDGMVEL.124
do i=1,imt IPDGMVEL.125
uisop(i,k)=0. IPDGMVEL.126
visopn(i,k)=0. IPDGMVEL.127
visops(i,k)=0. IPDGMVEL.128
enddo IPDGMVEL.129
enddo IPDGMVEL.130
c IPDGMVEL.131
do k=1,kmp1 IPDGMVEL.132
do i=1,imt IPDGMVEL.133
wisop(i,k)=0. IPDGMVEL.134
enddo IPDGMVEL.135
enddo IPDGMVEL.136
IPDGMVEL.137
end if IPDGMVEL.138
c IPDGMVEL.139
c for all the rows INCLUDING the first row of a task, shift IPDGMVEL.140
c slabwise row of "visopn" IPDGMVEL.141
c IPDGMVEL.142
IPDGMVEL.143
fxb=2. IPDGMVEL.144
do k=1,km IPDGMVEL.145
do i=1,imt IPDGMVEL.146
visops(i,k) = visopn(i,k) IPDGMVEL.147
enddo IPDGMVEL.148
enddo IPDGMVEL.149
IPDGMVEL.150
c IPDGMVEL.151
OLA2F403.367
IF (L_OVISBECK) THEN OLA2F403.368
OLA2F403.369
c Compute the meridional component of the isopycnal mixing velocity OLA2F403.370
c at the centre of the northern face of the tracer grid box. OLA2F403.371
do k=2,kmm1 OLA2F403.372
DO I=1,IMT OOM1F404.1
visopn(i,k) = -dz2r(k)*fm(i,k)*fmp(i,k)*(fk2(i,k-1,3)* OOM1F404.2
& (ATHKDFTV(i,j)/ahi(k-1))-fk2(i,k+1,3)* OOM1F404.3
& (ATHKDFTV(i,j)/ahi(k+1))) OOM1F404.4
ENDDO OOM1F404.5
OOM1F404.6
enddo OLA2F403.378
c Consider the top and bottom levels. "fk2" is assumed to be zero OLA2F403.379
c at the ocean top and bottom. OLA2F403.380
k = 1 OLA2F403.381
DO I=1,IMT OOM1F404.7
visopn(i,k) = dz2r(k)*fm(i,k)*fmp(i,k) OOM1F404.8
& *(fk2(i,k,3)*(ATHKDFTV(i,j)/ahi(k)) + OOM1F404.9
& fk2(i,k+1,3)*(ATHKDFTV(i,j)/ahi(k+1))) OOM1F404.10
ENDDO OOM1F404.11
OOM1F404.12
do i=1,imt OLA2F403.387
visopn(i,km) = 0. OLA2F403.388
enddo OLA2F403.389
do i=1,imt OLA2F403.390
k = min(kmt(i),kmtp(i)) OLA2F403.391
if (k .ne. 0) then OLA2F403.392
visopn(i,k) = -dz2r(k)*fm(i,k)*fmp(i,k) OOM1F404.13
& *(fk2(i,k,3)*(ATHKDFTV(i,j)/ahi(k))+ OOM1F404.14
& fk2(i,k-1,3)*(ATHKDFTV(i,j)/ahi(k-1))) OOM1F404.15
OOM1F404.16
endif OLA2F403.396
enddo OLA2F403.397
c Compute the zonal component of the isopycnal mixing velocity OLA2F403.398
c at the centre of the eastern face of the tracer grid box. OLA2F403.399
do k=2,kmm1 OLA2F403.400
DO I=1,IMTM1 OOM1F404.17
uisop(i,k) = -dz2r(k)*fm(i,k)*fm(i+1,k)*(fk1(i,k-1,3)* OOM1F404.18
& (ATHKDFTU(i,j)/ahi(k-1))-fk1(i,k+1,3)* OOM1F404.19
& (ATHKDFTU(i,j)/ahi(k+1))) OOM1F404.20
ENDDO OOM1F404.21
OOM1F404.22
enddo OLA2F403.406
c Consider the top and bottom levels. "fk1" is assumed to be zero OLA2F403.407
c at the ocean top and bottom. OLA2F403.408
k = 1 OLA2F403.409
DO I=1,IMTM1 OOM1F404.23
uisop(i,k) = dz2r(k)*fm(i,k)*fm(i+1,k) OOM1F404.24
& *(fk1(i,k,3)*(ATHKDFTU(i,j)/ahi(k))+ OOM1F404.25
& fk1(i,k+1,3)*(ATHKDFTU(i,j)/ahi(k+1))) OOM1F404.26
ENDDO OOM1F404.27
OOM1F404.28
do i=1,imtm1 OLA2F403.415
uisop(i,km) = 0. OLA2F403.416
enddo OLA2F403.417
do i=1,imtm1 OLA2F403.418
k = min(kmt(i),kmt(i+1)) OLA2F403.419
if (k .ne. 0) then OLA2F403.420
uisop(i,k) = -dz2r(k)*fm(i,k)*fm(i+1,k) OOM1F404.29
& *(fk1(i,k,3)*(ATHKDFTU(i,j)/ahi(k))+ OOM1F404.30
& fk1(i,k-1,3)*(ATHKDFTU(i,j)/ahi(k-1))) OOM1F404.31
OOM1F404.32
endif OLA2F403.424
enddo OLA2F403.425
OLA2F403.426
ELSE OLA2F403.427
OLA2F403.428
C The following occurs if the Visbeck Scheme is not chosen OLA2F403.429
c compute the meridional component of the isopycnal mixing velocity IPDGMVEL.152
c at the center of the northern face of the "t" grid box. IPDGMVEL.153
c IPDGMVEL.154
c the top and bottom levels will be considered later. IPDGMVEL.155
c IPDGMVEL.156
do k=2,kmm1 IPDGMVEL.157
DO I=1,IMT OOM1F404.33
visopn(i,k) = -dz2r(k)*fm(i,k)*fmp(i,k)*(fk2(i,k-1,3)* OOM1F404.34
& (athkdf(k-1)/ahi(k-1))-fk2(i,k+1,3)* OOM1F404.35
& (athkdf(k+1)/ahi(k+1))) OOM1F404.36
ENDDO OOM1F404.37
OOM1F404.38
enddo IPDGMVEL.163
IPDGMVEL.164
c consider the top and bottom levels. "fk2" is assumed to be zero IPDGMVEL.165
c at the ocean top and bottom. IPDGMVEL.166
c IPDGMVEL.167
k = 1 IPDGMVEL.168
DO I=1,IMT OOM1F404.39
visopn(i,k) = dz2r(k)*fm(i,k)*fmp(i,k) OOM1F404.40
& *(fk2(i,k,3)*(athkdf(k)/ahi(k))+fk2(i,k+1,3)* OOM1F404.41
& (athkdf(k+1)/ahi(k+1))) OOM1F404.42
ENDDO OOM1F404.43
OOM1F404.44
IPDGMVEL.174
do i=1,imt IPDGMVEL.175
visopn(i,km) = 0. IPDGMVEL.176
enddo IPDGMVEL.177
IPDGMVEL.178
IPDGMVEL.179
do i=1,imt IPDGMVEL.180
k = min(kmt(i),kmtp(i)) IPDGMVEL.181
IF (K .NE. 0) THEN OOM1F404.45
visopn(i,k) = -dz2r(k)*fm(i,k)*fmp(i,k) OOM1F404.46
& *(fk2(i,k,3)*(athkdf(k)/ahi(k))+fk2(i,k-1,3)* OOM1F404.47
& (athkdf(k-1)/ahi(k-1))) OOM1F404.48
ENDIF OOM1F404.49
OOM1F404.50
enddo IPDGMVEL.187
c IPDGMVEL.188
c compute the zonal component of the isopycnal mixing velocity IPDGMVEL.189
c at the center of the eastern face of the "t" grid box. IPDGMVEL.190
c IPDGMVEL.191
c the top and bottom levels will be considered later. IPDGMVEL.192
c IPDGMVEL.193
do k=2,kmm1 IPDGMVEL.194
DO I=1,IMTM1 OOM1F404.51
uisop(i,k) = -dz2r(k)*fm(i,k)*fm(i+1,k)*(fk1(i,k-1,3)* OOM1F404.52
& (athkdf(k-1)/ahi(k-1)) - fk1(i,k+1,3)* OOM1F404.53
& (athkdf(k+1)/ahi(k+1))) OOM1F404.54
ENDDO OOM1F404.55
OOM1F404.56
enddo IPDGMVEL.200
IPDGMVEL.201
c consider the top and bottom levels. "fk1" is assumed to be zero IPDGMVEL.202
c at the ocean top and bottom. IPDGMVEL.203
c IPDGMVEL.204
k = 1 IPDGMVEL.205
DO I=1,IMTM1 OOM1F404.57
uisop(i,k) = dz2r(k)*fm(i,k)*fm(i+1,k) OOM1F404.58
& *(fk1(i,k,3)*(athkdf(k)/ahi(k))+fk1(i,k+1,3)* OOM1F404.59
& (athkdf(k+1)/ahi(k+1))) OOM1F404.60
ENDDO OOM1F404.61
OOM1F404.62
IPDGMVEL.211
do i=1,imtm1 IPDGMVEL.212
uisop(i,km) = 0. IPDGMVEL.213
enddo IPDGMVEL.214
IPDGMVEL.215
do i=1,imtm1 IPDGMVEL.216
k = min(kmt(i),kmt(i+1)) IPDGMVEL.217
IF (K .NE. 0) THEN OOM1F404.63
uisop(i,k) = -dz2r(k)*fm(i,k)*fm(i+1,k) OOM1F404.64
& *(fk1(i,k,3)*(athkdf(k)/ahi(k))+fk1(i,k-1,3)* OOM1F404.65
& (athkdf(k-1)/ahi(k-1))) OOM1F404.66
ENDIF OOM1F404.67
enddo IPDGMVEL.223
OLA2F403.430
ENDIF ! if L_OVISBECK OLA2F403.431
OLA2F403.432
c IPDGMVEL.224
c set the boundary conditions at "i"="imt" IPDGMVEL.225
c IPDGMVEL.226
do k=1,km IPDGMVEL.227
IF (L_OCYCLIC) THEN IPDGMVEL.228
uisop(imt,k) = uisop(2,k) IPDGMVEL.229
ELSE IPDGMVEL.230
uisop(imt,k) = 0. IPDGMVEL.231
ENDIF IPDGMVEL.232
enddo IPDGMVEL.233
IPDGMVEL.234
c compute the vertical component of the isopycnal mixing velocity IPDGMVEL.235
c at the center of the top face of the "t" grid box, using the IPDGMVEL.236
c continuity equation for the isopycnal mixing velocities IPDGMVEL.237
c IPDGMVEL.238
do i=1,imt IPDGMVEL.239
wisop(i,1) = 0. IPDGMVEL.240
enddo IPDGMVEL.241
IPDGMVEL.242
do k=1,kmm1 IPDGMVEL.243
do i=2,imt IPDGMVEL.244
wisop(i,k+1) = dz(k)*cstr(j)*((uisop(i,k)-uisop(i-1,k)) IPDGMVEL.245
& *dxtr(i)+(visopn(i,k)*cs(j) IPDGMVEL.246
& -visops(i,k)*cs(j-1))*dytr(j)) IPDGMVEL.247
enddo IPDGMVEL.248
enddo IPDGMVEL.249
IPDGMVEL.250
C integrate downward from the surface IPDGMVEL.251
do k=1,kmm1 IPDGMVEL.252
do i=2,imt IPDGMVEL.253
wisop(i,k+1) = wisop(i,k)+wisop(i,k+1) IPDGMVEL.254
enddo IPDGMVEL.255
enddo IPDGMVEL.256
IPDGMVEL.257
C set velocities at bottom to zero IPDGMVEL.258
do i=2,imt IPDGMVEL.259
wisop(i,kmt(i)+1) = 0. IPDGMVEL.260
enddo IPDGMVEL.261
IPDGMVEL.262
IPDGMVEL.263
c set the boundary conditions at "i"=1 IPDGMVEL.264
c IPDGMVEL.265
do k=1,kmp1 IPDGMVEL.266
IF (L_OCYCLIC) THEN IPDGMVEL.267
wisop(1,k) = wisop(imtm1,k) IPDGMVEL.268
ELSE IPDGMVEL.269
wisop(1,k) = 0. IPDGMVEL.270
ENDIF IPDGMVEL.271
enddo IPDGMVEL.272
*ENDIF IPDGMVEL.273
return IPDGMVEL.274
end IPDGMVEL.275
*ENDIF IPDGMVEL.276