pkg/gmredi/gmredi_calc_bates_k.F

0d1e4b948d Mich*0001 #include "GMREDI_OPTIONS.h"
                0002
7ea279c259 Jean*0003 C     !ROUTINE: GMREDI_CALC_BATES_K
0d1e4b948d Mich*0004 C     !INTERFACE:
7ea279c259 Jean*0005       SUBROUTINE GMREDI_CALC_BATES_K(
                0006      I                  iMin, iMax, jMin, jMax,
                0007      I                  sigmaX, sigmaY, sigmaR,
                0008      I                  bi, bj, myTime, myIter, myThid )
0d1e4b948d Mich*0009
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
7ea279c259 Jean*0012 C     | SUBROUTINE GMREDI_CALC_BATES_K
                0013 C     | o Calculates the 3D diffusivity as per Bates et al. (2014)
0d1e4b948d Mich*0014 C     *==========================================================*
                0015 C     \ev
                0016
                0017       IMPLICIT NONE
                0018
                0019 C     == Global variables ==
                0020 #include "SIZE.h"
                0021 #include "EEPARAMS.h"
                0022 #include "PARAMS.h"
a08df3a232 Jean*0023 #include "GRID.h"
                0024 #include "DYNVARS.h"
0ff5c71d8d Jean*0025 #include "FFIELDS.h"
0d1e4b948d Mich*0026 #include "GMREDI.h"
                0027
                0028 C     !INPUT/OUTPUT PARAMETERS:
                0029 C     == Routine arguments ==
                0030 C     bi, bj    :: tile indices
7ea279c259 Jean*0031 C     myTime    :: Current time in simulation
                0032 C     myIter    :: Current iteration number in simulation
0d1e4b948d Mich*0033 C     myThid    :: My Thread Id. number
                0034
                0035       INTEGER iMin,iMax,jMin,jMax
4578baf364 Jean*0036       _RL sigmaX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0037       _RL sigmaY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0038       _RL sigmaR(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0d1e4b948d Mich*0039       INTEGER bi, bj
5a6ef5c2b4 Mich*0040       _RL myTime
7ea279c259 Jean*0041       INTEGER myIter
0d1e4b948d Mich*0042       INTEGER myThid
                0043
05118ac017 Jean*0044 #ifdef GM_BATES_K3D
0d1e4b948d Mich*0045
5a6ef5c2b4 Mich*0046 C     === Functions ====
                0047       LOGICAL  DIFFERENT_MULTIPLE
                0048       EXTERNAL DIFFERENT_MULTIPLE
                0049
0d1e4b948d Mich*0050 C     !LOCAL VARIABLES:
                0051 C     == Local variables ==
74e46366ba Davi*0052       INTEGER i,j,k,kk,m,kp1
5a6ef5c2b4 Mich*0053
                0054 C     update_modes :: Whether to update the eigenmodes
                0055       LOGICAL update_modes
                0056
                0057 C     surfk  :: index of the depth of the surface layer
                0058 C     kLow_C :: Local version of the index of deepest wet grid cell on tracer grid
                0059 C     kLow_U :: Local version of the index of deepest wet grid cell on U grid
                0060 C     kLow_V :: Local version of the index of deepest wet grid cell on V grid
4578baf364 Jean*0061       INTEGER surfk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0062       INTEGER kLow_C(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       INTEGER kLow_U(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064       INTEGER kLow_V(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5a6ef5c2b4 Mich*0065
                0066 C     N2loc  :: local N**2
                0067 C     slope  :: local slope
                0068 C     Req    :: local equatorial deformation radius (m)
                0069 C     deltaH :: local thickness of Eady integration (m)
                0070 C     g_reciprho_sq :: (gravity*recip_rhoConst)**2
                0071 C     M4loc  :: local M**4
05118ac017 Jean*0072 C     maxDRhoDz :: maximum value of d(rho)/dz (derived from GM_Bates_minN2)
5a6ef5c2b4 Mich*0073 C     sigx   :: local d(rho)/dx
                0074 C     sigy   :: local d(rho)/dy
                0075 C     sigz   :: local d(rho)/dz
                0076 C     hsurf  :: local surface layer depth
                0077 C     small  :: a small number (to avoid floating point exceptions)
4578baf364 Jean*0078       _RL N2loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0d1e4b948d Mich*0079       _RL slope
4578baf364 Jean*0080       _RL slopeC(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
5a6ef5c2b4 Mich*0081       _RL Req
4578baf364 Jean*0082       _RL deltaH(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
0d1e4b948d Mich*0083       _RL g_reciprho_sq
4578baf364 Jean*0084       _RL M4loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0085       _RL M4onN2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
0d1e4b948d Mich*0086       _RL maxDRhoDz
                0087       _RL sigx, sigy, sigz
74e46366ba Davi*0088       _RL hsurf, mskp1
0d1e4b948d Mich*0089       _RL small
5a6ef5c2b4 Mich*0090
74e46366ba Davi*0091 C     dfdy    :: gradient of the Coriolis paramater, df/dy, 1/(m*s)
                0092 C     dfdx    :: gradient of the Coriolis paramater, df/dx, 1/(m*s)
0ea5811d0e Mich*0093 C     Rurms  :: a local mixing length used in calculation of urms (m)
5a6ef5c2b4 Mich*0094 C     RRhines :: The Rhines scale (m)
                0095 C     Rmix    :: Mixing length
                0096 C     N2      :: Square of the buoyancy frequency (1/s**2)
                0097 C     N2W     :: Square of the buoyancy frequency (1/s**2) averaged to west of grid cell
                0098 C     N2S     :: Square of the buoyancy frequency (1/s**2) averaged to south of grid cell
                0099 C     N       :: Buoyancy frequency, SQRT(N2)
                0100 C     BVint   :: The vertical integral of N (m/s)
                0101 C     ubar    :: Zonal velocity on a tracer point (m/s)
                0102 C     Ubaro   :: Barotropic velocity (m/s)
4578baf364 Jean*0103       _RL dfdy(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0104       _RL dfdx(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL dummy(  1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL Rurms(  1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL RRhines(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL Rmix(   1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109       _RL N2(     1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0110       _RL N2W(    1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0111       _RL N2S(    1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0112       _RL N(      1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0113       _RL BVint(  1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL Ubaro(  1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL ubar(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0116
                0117       _RL tmpU(  1-OLx:sNx+OLx,1-OLy:sNy+OLy )
                0118       _RL tmpV(  1-OLx:sNx+OLx,1-OLy:sNy+OLy )
                0119       _RL uFldX( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr )
                0120       _RL vFldY( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr )
5a6ef5c2b4 Mich*0121
                0122 C     Rmid     :: Rossby radius (m)
                0123 C     KPV      :: Diffusivity (m**2/s)
40739f94bb Mich*0124 C     Kdqdx    :: diffusivity multiplied by zonal PV gradient
                0125 C     Kdqdy    :: diffusivity multiplied by meridional PV gradient
5a6ef5c2b4 Mich*0126 C     SlopeX   :: isopycnal slope in x direction
                0127 C     SlopeY   :: isopycnal slope in y direction
                0128 C     dSigmaDx :: sigmaX averaged onto tracer grid
                0129 C     dSigmaDy :: sigmaY averaged onto tracer grid
                0130 C     tfluxX   :: thickness flux in x direction
                0131 C     tfluxY   :: thickness flux in y direction
                0132 C     fCoriU   :: Coriolis parameter averaged to U points
                0133 C     fCoriV   :: Coriolis parameter averaged to V points
5dbdc05725 Davi*0134 C     coriU    :: As fCoriU, but limited
                0135 C     coriV    :: As fCoriV, but limited
5a6ef5c2b4 Mich*0136 C     surfkz   :: Depth of surface layer (in r units)
4578baf364 Jean*0137       _RL Rmid(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0138       _RL KPV(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0139       _RL Kdqdy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0140       _RL Kdqdx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0141       _RL SlopeX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0142       _RL SlopeY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0143       _RL dSigmaDx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0144       _RL dSigmaDy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0145       _RL tfluxX(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0146       _RL tfluxY(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0147       _RL coriU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0148       _RL coriV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0149       _RL fCoriU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0150       _RL fCoriV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0151       _RL surfkz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
5a6ef5c2b4 Mich*0152
df226773af Mich*0153 C     centreX,centreY :: used for calculating averages at centre of cell
                0154 C     numerator,denominator :: of the renormalisation factor
                0155 C     uInt      :: column integral of u velocity (sum u*dz)
                0156 C     vInt      :: column integral of v velocity (sum v*dz)
                0157 C     KdqdxInt  :: column integral of K*dqdx (sum K*dqdx*dz)
                0158 C     KdqdyInt  :: column integral of K*dqdy (sum K*dqdy*dz)
                0159 C     uKdqdyInt :: column integral of u*K*dqdy (sum u*K*dqdy*dz)
                0160 C     vKdqdxInt :: column integral of v*K*dqdx (sum v*K*dqdx*dz)
                0161 C     uXiyInt   :: column integral of u*Xiy (sum u*Xiy*dz)
                0162 C     vXixInt   :: column integral of v*Xix (sum v*Xix*dz)
                0163 C     Renorm    :: renormalisation factor at the centre of a cell
                0164 C     RenormU   :: renormalisation factor at the western face of a cell
                0165 C     RenormV   :: renormalisation factor at the southern face of a cell
224bd236ac Mich*0166       _RL centreX, centreY
                0167       _RL numerator, denominator
4578baf364 Jean*0168       _RL uInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0169       _RL vInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0170       _RL KdqdxInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0171       _RL KdqdyInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0172       _RL uKdqdyInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0173       _RL vKdqdxInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0174       _RL uXiyInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0175       _RL vXixInt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0176       _RL Renorm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0177       _RL RenormU(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0178       _RL RenormV(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
224bd236ac Mich*0179
5a6ef5c2b4 Mich*0180 C     gradqx   :: Potential vorticity gradient in x direction
                0181 C     gradqy   :: Potential vorticity gradient in y direction
                0182 C     XimX     :: Vertical integral of phi_m*K*gradqx
                0183 C     XimY     :: Vertical integral of phi_m*K*gradqy
                0184 C     cDopp    :: Quasi-Doppler shifted long Rossby wave speed (m/s)
                0185 C     umc      :: ubar-c (m/s)
                0186 C     eady     :: Eady growth rate (1/s)
                0187 C     urms     :: the rms eddy velocity (m/s)
                0188 C     supp     :: The suppression factor
                0189 C     ustar    :: The eddy induced velocity in the x direction
                0190 C     vstar    :: The eddy induced velocity in the y direction
                0191 C     Xix      :: Xi in the x direction (m/s**2)
                0192 C     Xiy      :: Xi in the y direction (m/s**2)
4578baf364 Jean*0193       _RL gradqx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0194       _RL gradqy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0195       _RL XimX(GM_Bates_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0196       _RL XimY(GM_Bates_NModes,1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0197       _RL cDopp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0198       _RL umc(  1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
                0199       _RL eady( 1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0200       _RL urms( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
                0201       _RL supp( 1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
                0202       _RL ustar(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
                0203       _RL vstar(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
                0204       _RL Xix(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0205       _RL Xiy(   1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
05118ac017 Jean*0206 #ifdef GM_BATES_PASSIVE
5a6ef5c2b4 Mich*0207 C     psistar :: eddy induced streamfunction in the y direction
4578baf364 Jean*0208       _RL psistar(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1:Nr)
5a6ef5c2b4 Mich*0209 #endif
0d1e4b948d Mich*0210
                0211 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0212
                0213 C     ======================================
                0214 C     Initialise some variables
                0215 C     ======================================
                0216       small = TINY(zeroRL)
5a6ef5c2b4 Mich*0217       update_modes=.FALSE.
05118ac017 Jean*0218       IF ( DIFFERENT_MULTIPLE(GM_Bates_vecFreq,myTime,deltaTClock) )
5a6ef5c2b4 Mich*0219      &     update_modes=.TRUE.
                0220
4578baf364 Jean*0221       DO j=1-OLy,sNy+OLy
                0222        DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0223         kLow_C(i,j) = kLowC(i,j,bi,bj)
                0224        ENDDO
                0225       ENDDO
4578baf364 Jean*0226       DO j=1-OLy,sNy+OLy
                0227        DO i=1-OLx+1,sNx+OLx
0d1e4b948d Mich*0228         kLow_U(i,j) = MIN( kLow_C(i,j), kLow_C(i-1,j) )
                0229        ENDDO
                0230       ENDDO
4578baf364 Jean*0231       DO j=1-OLy+1,sNy+OLy
                0232        DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0233         kLow_V(i,j) = MIN( kLow_C(i,j), kLow_C(i,j-1) )
                0234        ENDDO
                0235       ENDDO
                0236
219a8971b6 Mich*0237 C     Dummy values for the edges. This does not affect the results
                0238 C     but avoids problems when solving for the eigenvalues.
4578baf364 Jean*0239       i=1-OLx
                0240       DO j=1-OLy,sNy+OLy
219a8971b6 Mich*0241        kLow_U(i,j) = 0
0d1e4b948d Mich*0242       ENDDO
4578baf364 Jean*0243       j=1-OLy
                0244       DO i=1-OLx,sNx+OLx
219a8971b6 Mich*0245        kLow_V(i,j) = 0
0d1e4b948d Mich*0246       ENDDO
                0247
                0248       g_reciprho_sq = (gravity*recip_rhoConst)**2
                0249 C     Gradient of Coriolis
4578baf364 Jean*0250       DO j=1-OLy+1,sNy+OLy
                0251        DO i=1-OLx+1,sNx+OLx
0d1e4b948d Mich*0252         dfdx(i,j) = ( fCori(i,j,bi,bj)-fCori(i-1,j,bi,bj) )
                0253      &              *recip_dxC(i,j,bi,bj)
                0254         dfdy(i,j) = ( fCori(i,j,bi,bj)-fCori(i,j-1,bi,bj) )
                0255      &              *recip_dyC(i,j,bi,bj)
                0256        ENDDO
                0257       ENDDO
                0258
5dbdc05725 Davi*0259 C     Coriolis at U and V points enforcing a minimum value so
0d1e4b948d Mich*0260 C     that it is defined at the equator
4578baf364 Jean*0261       DO j=1-OLy,sNy+OLy
                0262        DO i=1-OLx+1,sNx+OLx
9d5e80ac5d Mich*0263 C       Not limited
5dbdc05725 Davi*0264         fCoriU(i,j)= op5*( fCori(i,j,bi,bj)+fCori(i-1,j,bi,bj) )
                0265 C       Limited so that the inverse is defined at the equator
05118ac017 Jean*0266         coriU(i,j) = SIGN( MAX( ABS(fCoriU(i,j)),GM_Bates_minCori ),
5dbdc05725 Davi*0267      &                              fCoriU(i,j) )
9d5e80ac5d Mich*0268        ENDDO
                0269       ENDDO
4578baf364 Jean*0270       DO j=1-OLy+1,sNy+OLy
                0271        DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0272 C       Not limited
5dbdc05725 Davi*0273         fCoriV(i,j)= op5*( fCori(i,j,bi,bj)+fCori(i,j-1,bi,bj) )
                0274 C       Limited so that the inverse is defined at the equator
05118ac017 Jean*0275         coriV(i,j) = SIGN( MAX( ABS(fCoriV(i,j)),GM_Bates_minCori ),
5dbdc05725 Davi*0276      &                              fCoriV(i,j) )
0d1e4b948d Mich*0277        ENDDO
                0278       ENDDO
9d5e80ac5d Mich*0279 C     Some dummy values at the edges
4578baf364 Jean*0280       i=1-OLx
                0281       DO j=1-OLy,sNy+OLy
5dbdc05725 Davi*0282        fCoriU(i,j)= fCori(i,j,bi,bj)
05118ac017 Jean*0283        coriU(i,j) = SIGN( MAX( ABS(fCori(i,j,bi,bj)),GM_Bates_minCori ),
5dbdc05725 Davi*0284      &                             fCori(i,j,bi,bj) )
9d5e80ac5d Mich*0285       ENDDO
4578baf364 Jean*0286       j=1-OLy
                0287       DO i=1-OLx,sNx+OLx
5dbdc05725 Davi*0288        fCoriV(i,j)= fCori(i,j,bi,bj)
05118ac017 Jean*0289        coriV(i,j) = SIGN( MAX( ABS(fCori(i,j,bi,bj)),GM_Bates_minCori ),
5dbdc05725 Davi*0290      &                             fCori(i,j,bi,bj) )
9d5e80ac5d Mich*0291       ENDDO
0d1e4b948d Mich*0292
                0293 C     Zeroing some cumulative fields
4578baf364 Jean*0294       DO j=1-OLy,sNy+OLy
                0295        DO i=1-OLx,sNx+OLx
74b49ccfa8 Mich*0296         eady(i,j)   = zeroRL
                0297         BVint(i,j)  = zeroRL
                0298         Ubaro(i,j)  = zeroRL
                0299         deltaH(i,j) = zeroRL
                0300        ENDDO
                0301       ENDDO
                0302       DO k=1,Nr
4578baf364 Jean*0303        DO j=1-OLy,sNy+OLy
                0304         DO i=1-OLx,sNx+OLx
74b49ccfa8 Mich*0305          slopeC(i,j,k)=zeroRL
                0306         ENDDO
0d1e4b948d Mich*0307        ENDDO
                0308       ENDDO
                0309
9d5e80ac5d Mich*0310 C     initialise remaining 2d variables
4578baf364 Jean*0311       DO j=1-OLy,sNy+OLy
                0312        DO i=1-OLx,sNx+OLx
9d5e80ac5d Mich*0313         Rurms(i,j)=zeroRL
                0314         RRhines(i,j)=zeroRL
                0315         Rmix(i,j)=zeroRL
                0316        ENDDO
                0317       ENDDO
                0318 C     initialise remaining 3d variables
                0319       DO k=1,Nr
4578baf364 Jean*0320        DO j=1-OLy,sNy+OLy
                0321         DO i=1-OLx,sNx+OLx
05118ac017 Jean*0322          N2loc(i,j,k)=GM_Bates_minN2
                0323          N2W(i,j,k) = GM_Bates_minN2
                0324          N2S(i,j,k) = GM_Bates_minN2
9d5e80ac5d Mich*0325          M4loc(i,j,k)=zeroRL
                0326          M4onN2(i,j,k)=zeroRL
                0327          urms(i,j,k)=zeroRL
                0328          SlopeX(i,j,k)=zeroRL
                0329          SlopeY(i,j,k)=zeroRL
                0330          dSigmaDx(i,j,k)=zeroRL
                0331          dSigmaDy(i,j,k)=zeroRL
                0332          gradqx(i,j,k)=zeroRL
                0333          gradqy(i,j,k)=zeroRL
                0334         ENDDO
                0335        ENDDO
                0336       ENDDO
                0337
0d1e4b948d Mich*0338 C     Find the zonal velocity at the cell centre
0ff5c71d8d Jean*0339 #ifdef ALLOW_EDDYPSI
5853335f7f Mich*0340       IF (GM_InMomAsStress) THEN
93219a225f Mich*0341         DO k=1,Nr
4578baf364 Jean*0342          DO i = 1-OLx,sNx+OLx
                0343           DO j = 1-OLy,sNy+OLy
93219a225f Mich*0344            uFldX(i,j,k) = uEulerMean(i,j,k,bi,bj)
                0345            vFldY(i,j,k) = vEulerMean(i,j,k,bi,bj)
                0346           ENDDO
                0347          ENDDO
                0348         ENDDO
5853335f7f Mich*0349       ELSE
                0350 #endif
93219a225f Mich*0351         DO k=1,Nr
4578baf364 Jean*0352          DO i = 1-OLx,sNx+OLx
                0353           DO j = 1-OLy,sNy+OLy
93219a225f Mich*0354            uFldX(i,j,k) = uVel(i,j,k,bi,bj)
                0355            vFldY(i,j,k) = vVel(i,j,k,bi,bj)
                0356           ENDDO
                0357          ENDDO
                0358         ENDDO
0ff5c71d8d Jean*0359 #ifdef ALLOW_EDDYPSI
5853335f7f Mich*0360       ENDIF
                0361 #endif
0d1e4b948d Mich*0362
93219a225f Mich*0363 C     The following comes from rotate_uv2en_rl
                0364 C     This code does two things:
                0365 C     1) go from C grid velocity points to A grid velocity points
                0366 C     2) go from model grid directions to east/west directions
                0367       DO k = 1,Nr
4578baf364 Jean*0368        DO i = 1-OLx,sNx+OLx
                0369         j=sNy+OLy
93219a225f Mich*0370         tmpU(i,j)=zeroRL
                0371         tmpV(i,j)=zeroRL
                0372        ENDDO
4578baf364 Jean*0373        DO j = 1-OLy,sNy+OLy-1
                0374         i=sNx+OLx
93219a225f Mich*0375         tmpU(i,j)=zeroRL
                0376         tmpV(i,j)=zeroRL
4578baf364 Jean*0377         DO i = 1-OLx,sNx+OLx-1
93219a225f Mich*0378          tmpU(i,j) = 0.5 _d 0
                0379      &        *( uFldX(i+1,j,k) + uFldX(i,j,k) )
                0380          tmpV(i,j) = 0.5 _d 0
                0381      &        *( vFldY(i,j+1,k) + vFldY(i,j,k) )
74e46366ba Davi*0382
93219a225f Mich*0383          tmpU(i,j) = tmpU(i,j) * maskC(i,j,k,bi,bj)
                0384          tmpV(i,j) = tmpV(i,j) * maskC(i,j,k,bi,bj)
                0385         ENDDO
                0386        ENDDO
66842ef555 Davi*0387 C     Rotation
4578baf364 Jean*0388        DO j = 1-OLy,sNy+OLy
                0389         DO i = 1-OLx,sNx+OLx
74e46366ba Davi*0390          ubar(i,j,k) =
93219a225f Mich*0391      &        angleCosC(i,j,bi,bj)*tmpU(i,j)
                0392      &        -angleSinC(i,j,bi,bj)*tmpV(i,j)
                0393         ENDDO
                0394        ENDDO
                0395       ENDDO
74e46366ba Davi*0396
f183cca6ba Davi*0397 C     Calculate the barotropic velocity by vertically integrating
                0398 C     and the dividing by the depth of the water column
                0399 C     Note that Ubaro is at the C-point.
                0400       DO k=1,Nr
4578baf364 Jean*0401        DO j=1-OLy,sNy+OLy
                0402         DO i=1-OLx,sNx+OLx
f183cca6ba Davi*0403          Ubaro(i,j) = Ubaro(i,j) +
                0404      &        drF(k)*hfacC(i,j,k,bi,bj)*ubar(i,j,k)
                0405         ENDDO
                0406        ENDDO
                0407       ENDDO
4578baf364 Jean*0408       DO j=1-OLy,sNy+OLy
                0409        DO i=1-OLx,sNx+OLx
f183cca6ba Davi*0410         IF (kLow_C(i,j).GT.0) THEN
                0411 C         The minus sign is because r_Low<0
                0412           Ubaro(i,j) = -Ubaro(i,j)/r_Low(i,j,bi,bj)
                0413         ENDIF
                0414        ENDDO
                0415       ENDDO
                0416
0d1e4b948d Mich*0417 C     Square of the buoyancy frequency at the top of a grid cell
74e46366ba Davi*0418 C     Enforce a minimum N2
                0419 C     Mask N2, so it is zero at bottom
0d1e4b948d Mich*0420       DO k=2,Nr
4578baf364 Jean*0421        DO j=1-OLy,sNy+OLy
                0422         DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0423          N2(i,j,k) = -gravity*recip_rhoConst*sigmaR(i,j,k)
05118ac017 Jean*0424          N2(i,j,k) = MAX(N2(i,j,k),GM_Bates_minN2)*maskC(i,j,k,bi,bj)
74e46366ba Davi*0425          N(i,j,k)  = SQRT(N2(i,j,k))
0d1e4b948d Mich*0426         ENDDO
                0427        ENDDO
                0428       ENDDO
                0429 C     N2(k=1) is always zero
4578baf364 Jean*0430       DO j=1-OLy,sNy+OLy
                0431        DO i=1-OLx,sNx+OLx
74e46366ba Davi*0432         N2(i,j,1) = zeroRL
                0433         N(i,j,1)  = zeroRL
0d1e4b948d Mich*0434        ENDDO
                0435       ENDDO
                0436 C     Calculate the minimum drho/dz
05118ac017 Jean*0437       maxDRhoDz = -rhoConst*GM_Bates_minN2/gravity
0d1e4b948d Mich*0438
                0439 C     Integrate the buoyancy frequency vertically using the trapezoidal method.
74e46366ba Davi*0440 C     Assume that N(z=-H)=0
0d1e4b948d Mich*0441       DO k=1,Nr
74e46366ba Davi*0442        kp1 = min(k+1,Nr)
                0443        mskp1 = oneRL
                0444        IF ( k.EQ.Nr ) mskp1 = zeroRL
4578baf364 Jean*0445        DO j=1-OLy,sNy+OLy
                0446         DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0447            BVint(i,j) = BVint(i,j) + hFacC(i,j,k,bi,bj)*drF(k)
74e46366ba Davi*0448      &                         *op5*(N(i,j,k)+mskp1*N(i,j,kp1))
0d1e4b948d Mich*0449         ENDDO
                0450        ENDDO
                0451       ENDDO
                0452
                0453 C     Calculate the eigenvalues and eigenvectors
5a6ef5c2b4 Mich*0454       IF (update_modes) THEN
                0455         CALL GMREDI_CALC_EIGS(
                0456      I       iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
                0457      I       kLow_C, maskC(:,:,:,bi,bj),
                0458      I       hfacC(:,:,:,bi,bj), recip_hfacC(:,:,:,bi,bj),
                0459      I       R_Low(:,:,bi,bj), 1, .TRUE.,
                0460      O       Rmid, modesC(:,:,:,:,bi,bj))
                0461
                0462 C       Calculate the Rossby Radius
4578baf364 Jean*0463         DO j=1-OLy+1,sNy+OLy
                0464          DO i=1-OLx+1,sNx+OLx
f183cca6ba Davi*0465           Req  = SQRT(BVint(i,j)/(2. _d 0*pi*gradf(i,j,bi,bj)))
5a6ef5c2b4 Mich*0466           Rdef(i,j,bi,bj) = MIN(Rmid(i,j),Req)
                0467          ENDDO
                0468         ENDDO
                0469       ENDIF
0d1e4b948d Mich*0470
                0471 C     Average dsigma/dx and dsigma/dy onto the centre points
                0472       DO k=1,Nr
4578baf364 Jean*0473        DO j=1-OLy,sNy+OLy-1
                0474         DO i=1-OLx,sNx+OLx-1
0d1e4b948d Mich*0475          dSigmaDx(i,j,k) = op5*(sigmaX(i,j,k)+sigmaX(i+1,j,k))
                0476          dSigmaDy(i,j,k) = op5*(sigmaY(i,j,k)+sigmaY(i,j+1,k))
                0477         ENDDO
                0478        ENDDO
                0479       ENDDO
                0480
                0481 C     ===============================
                0482 C     Calculate the Eady growth rate
                0483 C     ===============================
                0484       DO k=1,Nr
                0485
74e46366ba Davi*0486        kp1 = min(k+1,Nr)
                0487        mskp1 = oneRL
                0488        IF ( k.EQ.Nr ) mskp1 = zeroRL
                0489
4578baf364 Jean*0490        DO j=1-OLy,sNy+OLy-1
                0491         DO i=1-OLx,sNx+OLx-1
0ff5c71d8d Jean*0492          M4loc(i,j,k) = g_reciprho_sq*( dSigmaDx(i,j,k)**2
74b49ccfa8 Mich*0493      &                                 +dSigmaDy(i,j,k)**2 )
74e46366ba Davi*0494          N2loc(i,j,k) = op5*(N2(i,j,k)+mskp1*N2(i,j,kp1))
74b49ccfa8 Mich*0495         ENDDO
                0496        ENDDO
0d1e4b948d Mich*0497 C      The bottom of the grid cell is shallower than the top
                0498 C      integration level, so, advance the depth.
05118ac017 Jean*0499        IF (-rF(k+1) .LE. GM_Bates_EadyMinDepth) CYCLE
0d1e4b948d Mich*0500
c0612052db Jean*0501 C      Do not bother going any deeper since the top of the
0d1e4b948d Mich*0502 C      cell is deeper than the bottom integration level
05118ac017 Jean*0503        IF (-rF(k).GE.GM_Bates_EadyMaxDepth) EXIT
0d1e4b948d Mich*0504
                0505 C      We are in the integration depth range
4578baf364 Jean*0506        DO j=1-OLy,sNy+OLy-1
                0507         DO i=1-OLx,sNx+OLx-1
0ff5c71d8d Jean*0508          IF ( (kLow_C(i,j).GE.k) .AND.
74b49ccfa8 Mich*0509      &        (-hMixLayer(i,j,bi,bj).LE.-rC(k)) ) THEN
0d1e4b948d Mich*0510
c2dd265de5 Mich*0511            slopeC(i,j,k) = SQRT(M4loc(i,j,k))/N2loc(i,j,k)
0d1e4b948d Mich*0512 C          Limit the slope.  Note, this is not all the Eady calculations.
40739f94bb Mich*0513            IF (slopeC(i,j,k).LE.GM_maxSlope) THEN
                0514              M4onN2(i,j,k) = M4loc(i,j,k)/N2loc(i,j,k)
0d1e4b948d Mich*0515            ELSE
40739f94bb Mich*0516              slopeC(i,j,k) = GM_maxslope
                0517              M4onN2(i,j,k) = SQRT(M4loc(i,j,k))*GM_maxslope
0d1e4b948d Mich*0518            ENDIF
0ff5c71d8d Jean*0519            eady(i,j)   = eady(i,j)
40739f94bb Mich*0520      &                   + hfacC(i,j,k,bi,bj)*drF(k)*M4onN2(i,j,k)
74b49ccfa8 Mich*0521            deltaH(i,j) = deltaH(i,j) + drF(k)
0d1e4b948d Mich*0522          ENDIF
                0523         ENDDO
                0524        ENDDO
                0525       ENDDO
                0526
4578baf364 Jean*0527       DO j=1-OLy,sNy+OLy
                0528        DO i=1-OLx,sNx+OLx
74b49ccfa8 Mich*0529 C       If the minimum depth for the integration is deeper than the ocean
0ff5c71d8d Jean*0530 C       bottom OR the mixed layer is deeper than the maximum depth of
74b49ccfa8 Mich*0531 C       integration, we set the Eady growth rate to something small
                0532 C       to avoid floating point exceptions.
                0533 C       Later, these areas will be given a small diffusivity.
                0534         IF (deltaH(i,j).EQ.zeroRL) THEN
0d1e4b948d Mich*0535           eady(i,j) = small
                0536
74b49ccfa8 Mich*0537 C       Otherwise, divide over the integration and take the square root
                0538 C       to actually find the Eady growth rate.
0d1e4b948d Mich*0539         ELSE
74b49ccfa8 Mich*0540           eady(i,j) = SQRT(eady(i,j)/deltaH(i,j))
0ff5c71d8d Jean*0541
0d1e4b948d Mich*0542         ENDIF
                0543
                0544        ENDDO
                0545       ENDDO
                0546
                0547 C     ======================================
                0548 C     Calculate the diffusivity
                0549 C     ======================================
4578baf364 Jean*0550       DO j=1-OLy+1,sNy+OLy
                0551        DO i=1-OLx+1,sNx+OLx-1
0d1e4b948d Mich*0552 C       Calculate the Visbeck velocity
05118ac017 Jean*0553         Rurms(i,j) = MIN(Rdef(i,j,bi,bj),GM_Bates_Rmax)
                0554         urms(i,j,1) = GM_Bates_Lambda*eady(i,j)*Rurms(i,j)
0d1e4b948d Mich*0555 C       Set the bottom urms to zero
                0556         k=kLow_C(i,j)
                0557         IF (k.GT.0) urms(i,j,k) = 0.0
                0558
                0559 C       Calculate the Rhines scale
66842ef555 Davi*0560         RRhines(i,j) = SQRT(urms(i,j,1)/gradf(i,j,bi,bj))
0d1e4b948d Mich*0561
                0562 C       Calculate the estimated length scale
5a6ef5c2b4 Mich*0563         Rmix(i,j) = MIN(Rdef(i,j,bi,bj), RRhines(i,j))
05118ac017 Jean*0564         Rmix(i,j) = MAX(Rmix(i,j),GM_Bates_Rmin)
0d1e4b948d Mich*0565
                0566 C       Calculate the Doppler shifted long Rossby wave speed
74e46366ba Davi*0567 C       Ubaro is at the C-point.
93219a225f Mich*0568         cDopp(i,j) = Ubaro(i,j)
66842ef555 Davi*0569      &                - gradf(i,j,bi,bj)*Rdef(i,j,bi,bj)*Rdef(i,j,bi,bj)
05118ac017 Jean*0570 C       Limit the wave speed to the namelist variable GM_Bates_maxC
                0571         IF (ABS(cDopp(i,j)).GT.GM_Bates_maxC) THEN
                0572           cDopp(i,j) = MAX(GM_Bates_maxC, cDopp(i,j))
0d1e4b948d Mich*0573         ENDIF
                0574
                0575        ENDDO
                0576       ENDDO
                0577
                0578 C     Project the surface urms to the subsurface using the first baroclinic mode
5a6ef5c2b4 Mich*0579       CALL GMREDI_CALC_URMS(
                0580      I     iMin,iMax,jMin,jMax,bi,bj,N2,myThid,
                0581      U     urms)
0d1e4b948d Mich*0582
                0583 C     Calculate the diffusivity (on the mass grid)
                0584       DO k=1,Nr
4578baf364 Jean*0585        DO j=1-OLy,sNy+OLy
                0586         DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0587          IF (k.LE.kLow_C(i,j)) THEN
74b49ccfa8 Mich*0588            IF (deltaH(i,j).EQ.zeroRL) THEN
05118ac017 Jean*0589              GM_BatesK3d(i,j,k,bi,bj) = GM_Bates_smallK
0d1e4b948d Mich*0590            ELSE
                0591              IF (urms(i,j,k).EQ.0.0) THEN
05118ac017 Jean*0592                GM_BatesK3d(i,j,k,bi,bj) = GM_Bates_smallK
0d1e4b948d Mich*0593              ELSE
74e46366ba Davi*0594               umc(i,j,k) =ubar(i,j,k) - cDopp(i,j)
05118ac017 Jean*0595               supp(i,j,k) = 1./
                0596      &               ( 1. + GM_Bates_b1*umc(i,j,k)**2/urms(i,j,1)**2 )
40739f94bb Mich*0597 C              2*Rmix gives the diameter
05118ac017 Jean*0598               GM_BatesK3d(i,j,k,bi,bj) = GM_Bates_gamma*urms(i,j,k)
74e46366ba Davi*0599      &                           *2.*Rmix(i,j)*supp(i,j,k)
0d1e4b948d Mich*0600              ENDIF
                0601
                0602 C            Enforce lower and upper bounds on the diffusivity
05118ac017 Jean*0603              GM_BatesK3d(i,j,k,bi,bj) =
                0604      &               MIN( GM_BatesK3d(i,j,k,bi,bj), GM_Bates_maxK )
                0605              GM_BatesK3d(i,j,k,bi,bj) =
                0606      &               MAX( GM_BatesK3d(i,j,k,bi,bj), GM_Bates_smallK )
0d1e4b948d Mich*0607            ENDIF
                0608          ENDIF
                0609         ENDDO
                0610        ENDDO
                0611       ENDDO
                0612
                0613 C     ======================================
                0614 C     Find the PV gradient
                0615 C     ======================================
4e65b6d0de Mich*0616 C     Calculate the surface layer thickness.
0ff5c71d8d Jean*0617 C     Use hMixLayer (calculated in model/src/calc_oce_mxlayer)
4e65b6d0de Mich*0618 C     for the mixed layer depth.
                0619
                0620 C     Enforce a minimum surface layer depth
4578baf364 Jean*0621       DO j=1-OLy,sNy+OLy
                0622        DO i=1-OLx,sNx+OLx
05118ac017 Jean*0623         surfkz(i,j) = MIN(-GM_Bates_surfMinDepth,-hMixLayer(i,j,bi,bj))
4e65b6d0de Mich*0624         surfkz(i,j) = MAX(surfkz(i,j),R_low(i,j,bi,bj))
                0625         IF(maskC(i,j,1,bi,bj).EQ.0.0) surfkz(i,j)=0.0
                0626         surfk(i,j) = 0
0d1e4b948d Mich*0627        ENDDO
                0628       ENDDO
5a6ef5c2b4 Mich*0629       DO k=1,Nr
4578baf364 Jean*0630        DO j=1-OLy,sNy+OLy
                0631         DO i=1-OLx,sNx+OLx
0ff5c71d8d Jean*0632          IF (rF(k).GT.surfkz(i,j) .AND. surfkz(i,j).GE.rF(k+1))
4e65b6d0de Mich*0633      &        surfk(i,j) = k
0d1e4b948d Mich*0634         ENDDO
                0635        ENDDO
                0636       ENDDO
4e65b6d0de Mich*0637
0d1e4b948d Mich*0638 C     Calculate the isopycnal slope
4578baf364 Jean*0639       DO j=1-OLy,sNy+OLy-1
                0640        DO i=1-OLx,sNx+OLx-1
0d1e4b948d Mich*0641         SlopeX(i,j,1) = zeroRL
                0642         SlopeY(i,j,1) = zeroRL
                0643        ENDDO
                0644       ENDDO
                0645       DO k=2,Nr
4578baf364 Jean*0646        DO j=1-OLy+1,sNy+OLy
                0647         DO i=1-OLx+1,sNx+OLx
0d1e4b948d Mich*0648          IF(surfk(i,j).GE.kLowC(i,j,bi,bj)) THEN
                0649 C          If the surface layer is thinner than the water column
                0650 C          the set the slope to zero to avoid problems.
                0651            SlopeX(i,j,k) = zeroRL
                0652            SlopeY(i,j,k) = zeroRL
                0653
                0654          ELSE
                0655 C          Calculate the zonal slope at the western cell face (U grid)
5a6ef5c2b4 Mich*0656            sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i-1,j,k)), maxDRhoDz )
0d1e4b948d Mich*0657            sigx  = op5*( sigmaX(i,j,k)+sigmaX(i,j,k-1) )
                0658            slope = sigx/sigz
                0659            IF(ABS(slope).GT.GM_maxSlope)
                0660      &          slope = SIGN(GM_maxSlope,slope)
                0661            SlopeX(i,j,k)=-maskW(i,j,k-1,bi,bj)*maskW(i,j,k,bi,bj)*slope
0ff5c71d8d Jean*0662
0d1e4b948d Mich*0663 C          Calculate the meridional slope at the southern cell face (V grid)
5a6ef5c2b4 Mich*0664            sigz  = MIN( op5*(sigmaR(i,j,k)+sigmaR(i,j-1,k)), maxDRhoDz )
0d1e4b948d Mich*0665            sigy  = op5*( sigmaY(i,j,k) + sigmaY(i,j,k-1) )
                0666            slope = sigy/sigz
                0667            IF(ABS(slope).GT.GM_maxSlope)
                0668      &          slope = SIGN(GM_maxSlope,slope)
                0669            SlopeY(i,j,k)=-maskS(i,j,k-1,bi,bj)*maskS(i,j,k,bi,bj)*slope
                0670          ENDIF
                0671         ENDDO
                0672        ENDDO
                0673       ENDDO
                0674
c8602656d9 Davi*0675 C     Calculate gradients of PV stretching term and PV diffusivity.
                0676 C     These may be altered later depending on namelist options.
0d1e4b948d Mich*0677 C     Enforce a zero slope bottom boundary condition for the bottom most cells (k=Nr)
                0678       k=Nr
4578baf364 Jean*0679       DO j=1-OLy,sNy+OLy
                0680        DO i=1-OLx,sNx+OLx
c8602656d9 Davi*0681 C          Zonal gradient of PV stretching term: at the western cell face
0d1e4b948d Mich*0682            tfluxX(i,j,k) = -fCoriU(i,j)*SlopeX(i,j,k)
                0683      &                      *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
c8602656d9 Davi*0684 C     Meridional gradient of PV stretching term: at the southern cell face
0d1e4b948d Mich*0685            tfluxY(i,j,k) = -fCoriV(i,j)*SlopeY(i,j,k)
                0686      &                     *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
40312237e2 Mich*0687
0ff5c71d8d Jean*0688 C          Use the interior diffusivity. Note: if we are using a
40312237e2 Mich*0689 C          constant diffusivity KPV is overwritten later
05118ac017 Jean*0690            KPV(i,j,k) = GM_BatesK3d(i,j,k,bi,bj)
40312237e2 Mich*0691
0d1e4b948d Mich*0692        ENDDO
                0693       ENDDO
                0694
c8602656d9 Davi*0695 C     Calculate gradients of PV stretching term and PV diffusivity: for other cells (k<Nr)
0d1e4b948d Mich*0696       DO k=Nr-1,1,-1
4578baf364 Jean*0697        DO j=1-OLy,sNy+OLy
                0698         DO i=1-OLx,sNx+OLx
c8602656d9 Davi*0699 C          Zonal gradient of PV stretching term: at the western cell face
40312237e2 Mich*0700          tfluxX(i,j,k)=-fCoriU(i,j)*(SlopeX(i,j,k)-SlopeX(i,j,k+1))
                0701      &        *recip_drF(k)*recip_hFacW(i,j,k,bi,bj)
                0702      &        *maskW(i,j,k,bi,bj)
c8602656d9 Davi*0703 C     Meridional gradient of PV stretching term: at the southern cell face
40312237e2 Mich*0704          tfluxY(i,j,k)=-fCoriV(i,j)*(SlopeY(i,j,k)-SlopeY(i,j,k+1))
                0705      &        *recip_drF(k)*recip_hFacS(i,j,k,bi,bj)
                0706      &        *maskS(i,j,k,bi,bj)
0ff5c71d8d Jean*0707
                0708 C        Use the interior diffusivity. Note: if we are using a
40312237e2 Mich*0709 C        constant diffusivity KPV is overwritten later
05118ac017 Jean*0710          KPV(i,j,k) = GM_BatesK3d(i,j,k,bi,bj)
40312237e2 Mich*0711         ENDDO
                0712        ENDDO
                0713       ENDDO
                0714
0ff5c71d8d Jean*0715 C     Special treatment for the surface layer (if necessary), which overwrites
40312237e2 Mich*0716 C     values in the previous loops.
05118ac017 Jean*0717       IF (GM_Bates_ThickSheet .OR. GM_Bates_surfK) THEN
40312237e2 Mich*0718         DO k=Nr-1,1,-1
4578baf364 Jean*0719          DO j=1-OLy,sNy+OLy
                0720           DO i=1-OLx,sNx+OLx
0ff5c71d8d Jean*0721            IF(k.LE.surfk(i,j)) THEN
                0722 C            We are in the surface layer.  Change the thickness flux
40312237e2 Mich*0723 C            and diffusivity as necessary.
                0724
05118ac017 Jean*0725              IF (GM_Bates_ThickSheet) THEN
40312237e2 Mich*0726 C              We are in the surface layer, so set the thickness flux
                0727 C              based on the average slope over the surface layer
                0728 C              If we are on the edge of a "cliff" the surface layer at the
0ff5c71d8d Jean*0729 C              centre of the grid point could be deeper than the U or V point.
40312237e2 Mich*0730 C              So, we ensure that we always take a sensible slope.
                0731                IF(kLow_U(i,j).LT.surfk(i,j)) THEN
                0732                  kk=kLow_U(i,j)
                0733                  hsurf = -rLowW(i,j,bi,bj)
                0734                ELSE
                0735                  kk=surfk(i,j)
                0736                  hsurf = -surfkz(i,j)
                0737                ENDIF
                0738                IF(kk.GT.0) THEN
                0739                  IF(kk.EQ.Nr) THEN
                0740                    tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
                0741      &                  *SlopeX(i,j,kk)/hsurf
                0742                  ELSE
                0743                    tfluxX(i,j,k) = -fCoriU(i,j)*maskW(i,j,k,bi,bj)
                0744      &                  *( SlopeX(i,j,kk)-SlopeX(i,j,kk+1) )/hsurf
                0745                  ENDIF
                0746                ELSE
                0747                  tfluxX(i,j,k) = zeroRL
                0748                ENDIF
0ff5c71d8d Jean*0749
40312237e2 Mich*0750                IF(kLow_V(i,j).LT.surfk(i,j)) THEN
                0751                  kk=kLow_V(i,j)
                0752                  hsurf = -rLowS(i,j,bi,bj)
                0753                ELSE
                0754                  kk=surfk(i,j)
                0755                  hsurf = -surfkz(i,j)
                0756                ENDIF
                0757                IF(kk.GT.0) THEN
                0758                  IF(kk.EQ.Nr) THEN
                0759                    tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
                0760      &                  *SlopeY(i,j,kk)/hsurf
                0761                  ELSE
                0762                    tfluxY(i,j,k) = -fCoriV(i,j)*maskS(i,j,k,bi,bj)
                0763      &                  *( SlopeY(i,j,kk)-SlopeY(i,j,kk+1) )/hsurf
                0764                  ENDIF
                0765                ELSE
                0766                  tfluxY(i,j,k) = zeroRL
                0767                ENDIF
5a6ef5c2b4 Mich*0768              ENDIF
0d1e4b948d Mich*0769
05118ac017 Jean*0770              IF (GM_Bates_surfK) THEN
40312237e2 Mich*0771 C              Use a constant K in the surface layer.
05118ac017 Jean*0772                KPV(i,j,k) = GM_Bates_constK
5a6ef5c2b4 Mich*0773              ENDIF
0d1e4b948d Mich*0774            ENDIF
40312237e2 Mich*0775           ENDDO
                0776          ENDDO
0d1e4b948d Mich*0777         ENDDO
40312237e2 Mich*0778       ENDIF
0d1e4b948d Mich*0779
                0780 C     Calculate gradq
05118ac017 Jean*0781       IF (GM_Bates_beta_eq_0) THEN
ecb6b844fc Mich*0782 C       Ignore beta in the calculation of grad(q)
ce74227608 Mich*0783         DO k=1,Nr
4578baf364 Jean*0784          DO j=1-OLy+1,sNy+OLy
                0785           DO i=1-OLx+1,sNx+OLx
ce74227608 Mich*0786            gradqx(i,j,k) = maskW(i,j,k,bi,bj)*tfluxX(i,j,k)
                0787            gradqy(i,j,k) = maskS(i,j,k,bi,bj)*tfluxY(i,j,k)
                0788           ENDDO
                0789          ENDDO
0d1e4b948d Mich*0790         ENDDO
0ff5c71d8d Jean*0791
ce74227608 Mich*0792       ELSE
                0793 C       Do not ignore beta
                0794         DO k=1,Nr
4578baf364 Jean*0795          DO j=1-OLy+1,sNy+OLy
                0796           DO i=1-OLx+1,sNx+OLx
ce74227608 Mich*0797            gradqx(i,j,k) = maskW(i,j,k,bi,bj)*(dfdx(i,j)+tfluxX(i,j,k))
                0798            gradqy(i,j,k) = maskS(i,j,k,bi,bj)*(dfdy(i,j)+tfluxY(i,j,k))
                0799           ENDDO
                0800          ENDDO
                0801         ENDDO
                0802       ENDIF
0d1e4b948d Mich*0803
                0804 C     ======================================
                0805 C     Find Xi and the eddy induced velocities
                0806 C     ======================================
                0807 C     Find the buoyancy frequency at the west and south faces of a cell
                0808 C     This is necessary to find the eigenvectors at those points
                0809       DO k=1,Nr
4578baf364 Jean*0810        DO j=1-OLy+1,sNy+OLy
                0811         DO i=1-OLx+1,sNx+OLx
0d1e4b948d Mich*0812          N2W(i,j,k) = maskW(i,j,k,bi,bj)
                0813      &                *( N2(i,j,k)+N2(i-1,j,k) )
                0814          N2S(i,j,k) = maskS(i,j,k,bi,bj)
                0815      &                *( N2(i,j,k)+N2(i,j-1,k) )
                0816         ENDDO
                0817        ENDDO
                0818       ENDDO
                0819
05118ac017 Jean*0820 C     If GM_Bates_use_constK=T, the diffusivity for the eddy transport is
                0821 C     set to a constant equal to GM_Bates_constK.
952f8e8fd1 Mich*0822 C     If the diffusivity is constant the method here is the same as GM.
05118ac017 Jean*0823 C     If GM_Bates_constRedi=T, BatesK3d will be set to GM_Bates_constK later.
                0824       IF(GM_Bates_use_constK) THEN
ce74227608 Mich*0825         DO k=1,Nr
4578baf364 Jean*0826          DO j=1-OLy,sNy+OLy
                0827           DO i=1-OLx,sNx+OLx
05118ac017 Jean*0828            KPV(i,j,k) = GM_Bates_constK
ce74227608 Mich*0829           ENDDO
                0830          ENDDO
                0831         ENDDO
4e65b6d0de Mich*0832       ENDIF
ce74227608 Mich*0833
05118ac017 Jean*0834       IF (.NOT. GM_Bates_smooth) THEN
4e65b6d0de Mich*0835 C       Do not expand K grad(q) => no smoothing
0ff5c71d8d Jean*0836 C       May only be done with a constant K, otherwise the
4e65b6d0de Mich*0837 C       integral constraint is violated.
ce74227608 Mich*0838         DO k=1,Nr
4578baf364 Jean*0839          DO j=1-OLy,sNy+OLy
                0840           DO i=1-OLx,sNx+OLx
4e65b6d0de Mich*0841            Xix(i,j,k) = -maskW(i,j,k,bi,bj)*KPV(i,j,k)*gradqx(i,j,k)
                0842            Xiy(i,j,k) = -maskS(i,j,k,bi,bj)*KPV(i,j,k)*gradqy(i,j,k)
ce74227608 Mich*0843           ENDDO
                0844          ENDDO
                0845         ENDDO
                0846
                0847       ELSE
4e65b6d0de Mich*0848 C       Expand K grad(q) in terms of baroclinic modes to smooth
                0849 C       and satisfy the integral constraint
ce74227608 Mich*0850
0d1e4b948d Mich*0851 C       Start with the X direction
                0852 C       ------------------------------
                0853 C       Calculate the eigenvectors at the West face of a cell
5a6ef5c2b4 Mich*0854         IF (update_modes) THEN
                0855           CALL GMREDI_CALC_EIGS(
05118ac017 Jean*0856      I         iMin,iMax,jMin,jMax,bi,bj, N2W, myThid,
                0857      I         kLow_U, maskW(:,:,:,bi,bj),
                0858      I         hfacW(:,:,:,bi,bj), recip_hfacW(:,:,:,bi,bj),
                0859      I         rLowW(:,:,bi,bj), GM_Bates_NModes, .FALSE.,
                0860      O         dummy, modesW(:,:,:,:,bi,bj) )
5a6ef5c2b4 Mich*0861         ENDIF
0ff5c71d8d Jean*0862
0d1e4b948d Mich*0863 C       Calculate Xi_m at the west face of a cell
4578baf364 Jean*0864         DO j=1-OLy,sNy+OLy
                0865          DO i=1-OLx,sNx+OLx
05118ac017 Jean*0866           DO m=1,GM_Bates_NModes
0d1e4b948d Mich*0867            XimX(m,i,j) = zeroRL
                0868           ENDDO
                0869          ENDDO
                0870         ENDDO
                0871         DO k=1,Nr
4578baf364 Jean*0872          DO j=1-OLy,sNy+OLy
                0873           DO i=1-OLx,sNx+OLx
05118ac017 Jean*0874            DO m=1,GM_Bates_NModes
40739f94bb Mich*0875             Kdqdx(i,j,k) = KPV(i,j,k)*gradqx(i,j,k)
0ff5c71d8d Jean*0876             XimX(m,i,j) = XimX(m,i,j)
40739f94bb Mich*0877      &           - maskW(i,j,k,bi,bj)*drF(k)*hfacW(i,j,k,bi,bj)
                0878      &           *Kdqdx(i,j,k)*modesW(m,i,j,k,bi,bj)
0d1e4b948d Mich*0879            ENDDO
                0880           ENDDO
                0881          ENDDO
                0882         ENDDO
0ff5c71d8d Jean*0883
0d1e4b948d Mich*0884 C     Calculate Xi in the X direction at the west face
                0885         DO k=1,Nr
4578baf364 Jean*0886          DO j=1-OLy,sNy+OLy
                0887           DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0888            Xix(i,j,k) = zeroRL
                0889           ENDDO
                0890          ENDDO
                0891         ENDDO
                0892         DO k=1,Nr
4578baf364 Jean*0893          DO j=1-OLy,sNy+OLy
                0894           DO i=1-OLx,sNx+OLx
05118ac017 Jean*0895            DO m=1,GM_Bates_NModes
0ff5c71d8d Jean*0896             Xix(i,j,k) = Xix(i,j,k)
5a6ef5c2b4 Mich*0897      &           + maskW(i,j,k,bi,bj)*XimX(m,i,j)*modesW(m,i,j,k,bi,bj)
0d1e4b948d Mich*0898            ENDDO
                0899           ENDDO
                0900          ENDDO
                0901         ENDDO
                0902
                0903 C     Now the Y direction
                0904 C     ------------------------------
                0905 C     Calculate the eigenvectors at the West face of a cell
5a6ef5c2b4 Mich*0906         IF (update_modes) THEN
                0907           CALL GMREDI_CALC_EIGS(
05118ac017 Jean*0908      I         iMin,iMax,jMin,jMax,bi,bj, N2S, myThid,
                0909      I         kLow_V, maskS(:,:,:,bi,bj),
                0910      I         hfacS(:,:,:,bi,bj), recip_hfacS(:,:,:,bi,bj),
                0911      I         rLowS(:,:,bi,bj), GM_Bates_NModes, .FALSE.,
                0912      O         dummy, modesS(:,:,:,:,bi,bj) )
5a6ef5c2b4 Mich*0913         ENDIF
                0914
4578baf364 Jean*0915         DO j=1-OLy,sNy+OLy
                0916          DO i=1-OLx,sNx+OLx
05118ac017 Jean*0917           DO m=1,GM_Bates_NModes
0d1e4b948d Mich*0918            XimY(m,i,j) = zeroRL
                0919           ENDDO
                0920          ENDDO
                0921         ENDDO
                0922         DO k=1,Nr
4578baf364 Jean*0923          DO j=1-OLy,sNy+OLy
                0924           DO i=1-OLx,sNx+OLx
05118ac017 Jean*0925            DO m=1,GM_Bates_NModes
40739f94bb Mich*0926             Kdqdy(i,j,k) = KPV(i,j,k)*gradqy(i,j,k)
4e65b6d0de Mich*0927             XimY(m,i,j) = XimY(m,i,j)
                0928      &           - drF(k)*hfacS(i,j,k,bi,bj)
40739f94bb Mich*0929      &           *Kdqdy(i,j,k)*modesS(m,i,j,k,bi,bj)
0d1e4b948d Mich*0930            ENDDO
                0931           ENDDO
                0932          ENDDO
                0933         ENDDO
0ff5c71d8d Jean*0934
0d1e4b948d Mich*0935 C     Calculate Xi for Y direction at the south face
                0936         DO k=1,Nr
4578baf364 Jean*0937          DO j=1-OLy,sNy+OLy
                0938           DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*0939            Xiy(i,j,k) = zeroRL
                0940           ENDDO
                0941          ENDDO
                0942         ENDDO
                0943         DO k=1,Nr
4578baf364 Jean*0944          DO j=1-OLy,sNy+OLy
                0945           DO i=1-OLx,sNx+OLx
05118ac017 Jean*0946            DO m=1,GM_Bates_NModes
0ff5c71d8d Jean*0947             Xiy(i,j,k) = Xiy(i,j,k)
5a6ef5c2b4 Mich*0948      &           + maskS(i,j,k,bi,bj)*XimY(m,i,j)*modesS(m,i,j,k,bi,bj)
0d1e4b948d Mich*0949            ENDDO
                0950           ENDDO
                0951          ENDDO
                0952         ENDDO
                0953
05118ac017 Jean*0954 C     ENDIF (.NOT. GM_Bates_smooth)
0d1e4b948d Mich*0955       ENDIF
                0956
224bd236ac Mich*0957 C     Calculate the renormalisation factor
4578baf364 Jean*0958       DO j=1-OLy,sNy+OLy
                0959        DO i=1-OLx,sNx+OLx
224bd236ac Mich*0960         uInt(i,j)=zeroRL
                0961         vInt(i,j)=zeroRL
                0962         KdqdyInt(i,j)=zeroRL
                0963         KdqdxInt(i,j)=zeroRL
                0964         uKdqdyInt(i,j)=zeroRL
                0965         vKdqdxInt(i,j)=zeroRL
                0966         uXiyInt(i,j)=zeroRL
                0967         vXixInt(i,j)=zeroRL
df226773af Mich*0968         Renorm(i,j)=oneRL
                0969         RenormU(i,j)=oneRL
                0970         RenormV(i,j)=oneRL
224bd236ac Mich*0971        ENDDO
                0972       ENDDO
                0973       DO k=1,Nr
4578baf364 Jean*0974        DO j=1-OLy,sNy+OLy-1
                0975         DO i=1-OLx,sNx+OLx-1
224bd236ac Mich*0976          centreX = op5*(uVel(i,j,k,bi,bj)+uVel(i+1,j,k,bi,bj))
                0977          centreY = op5*(Kdqdy(i,j,k)     +Kdqdy(i,j+1,k)     )
                0978 C        For the numerator
0ff5c71d8d Jean*0979          uInt(i,j) = uInt(i,j)
224bd236ac Mich*0980      &        + centreX*hfacC(i,j,k,bi,bj)*drF(k)
                0981          KdqdyInt(i,j) = KdqdyInt(i,j)
                0982      &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
                0983          uKdqdyInt(i,j) = uKdqdyInt(i,j)
                0984      &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
                0985 C        For the denominator
                0986          centreY = op5*(Xiy(i,j,k) + Xiy(i,j+1,k))
                0987          uXiyInt(i,j) = uXiyInt(i,j)
                0988      &        + centreX*centreY*hfacC(i,j,k,bi,bj)*drF(k)
                0989
                0990          centreX = op5*(Kdqdx(i,j,k)     +Kdqdx(i+1,j,k))
                0991          centreY = op5*(vVel(i,j,k,bi,bj)+vVel(i,j+1,k,bi,bj) )
                0992 C        For the numerator
                0993          vInt(i,j) = vInt(i,j)
                0994      &        + centreY*hfacC(i,j,k,bi,bj)*drF(k)
                0995          KdqdxInt(i,j) = KdqdxInt(i,j)
                0996      &        + CentreX*hfacC(i,j,k,bi,bj)*drF(k)
                0997          vKdqdxInt(i,j) = vKdqdxInt(i,j)
                0998      &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
                0999 C        For the denominator
                1000          centreX = op5*(Xix(i,j,k) + Xix(i+1,j,k))
                1001          vXixInt(i,j) = vXixInt(i,j)
                1002      &        + centreY*centreX*hfacC(i,j,k,bi,bj)*drF(k)
                1003
                1004         ENDDO
                1005        ENDDO
                1006       ENDDO
                1007
4578baf364 Jean*1008       DO j=1-OLy,sNy+OLy-1
                1009        DO i=1-OLx,sNx+OLx-1
224bd236ac Mich*1010         IF (kLowC(i,j,bi,bj).GT.0) THEN
0ff5c71d8d Jean*1011           numerator =
224bd236ac Mich*1012      &         (uKdqdyInt(i,j)-uInt(i,j)*KdqdyInt(i,j)/R_low(i,j,bi,bj))
                1013      &        -(vKdqdxInt(i,j)-vInt(i,j)*KdqdxInt(i,j)/R_low(i,j,bi,bj))
                1014           denominator = uXiyInt(i,j) - vXixInt(i,j)
0ff5c71d8d Jean*1015 C         We can have troubles with floating point exceptions if the denominator
                1016 C         of the renormalisation if the ocean is resting (e.g. intial conditions).
df226773af Mich*1017 C         So we make the renormalisation factor one if the denominator is very small
                1018 C         The renormalisation factor is supposed to correct the error in the extraction of
                1019 C         potential energy associated with the truncation of the expansion. Thus, we
                1020 C         enforce a minimum value for the renormalisation factor.
                1021 C         We also enforce a maximum renormalisation factor.
                1022           IF (denominator.GT.small) THEN
224bd236ac Mich*1023             Renorm(i,j) = ABS(numerator/denominator)
05118ac017 Jean*1024             Renorm(i,j) = MAX(Renorm(i,j),GM_Bates_minRenorm)
                1025             Renorm(i,j) = MIN(Renorm(i,j),GM_Bates_maxRenorm)
224bd236ac Mich*1026           ENDIF
                1027         ENDIF
                1028        ENDDO
                1029       ENDDO
                1030 C     Now put it back on to the velocity grids
4578baf364 Jean*1031       DO j=1-OLy+1,sNy+OLy-1
                1032        DO i=1-OLx+1,sNx+OLx-1
224bd236ac Mich*1033         RenormU(i,j) = op5*(Renorm(i-1,j)+Renorm(i,j))
                1034         RenormV(i,j) = op5*(Renorm(i,j-1)+Renorm(i,j))
                1035        ENDDO
                1036       ENDDO
                1037
0d1e4b948d Mich*1038 C     Calculate the eddy induced velocity in the X direction at the west face
                1039       DO k=1,Nr
4578baf364 Jean*1040        DO j=1-OLy+1,sNy+OLy
                1041         DO i=1-OLx+1,sNx+OLx
df226773af Mich*1042          ustar(i,j,k) = -RenormU(i,j)*Xix(i,j,k)/coriU(i,j)
0d1e4b948d Mich*1043         ENDDO
                1044        ENDDO
0ff5c71d8d Jean*1045       ENDDO
0d1e4b948d Mich*1046
                1047 C     Calculate the eddy induced velocity in the Y direction at the south face
                1048       DO k=1,Nr
4578baf364 Jean*1049        DO j=1-OLy+1,sNy+OLy
                1050         DO i=1-OLx+1,sNx+OLx
df226773af Mich*1051          vstar(i,j,k) = -RenormV(i,j)*Xiy(i,j,k)/coriV(i,j)
0d1e4b948d Mich*1052         ENDDO
                1053        ENDDO
0ff5c71d8d Jean*1054       ENDDO
0d1e4b948d Mich*1055
                1056 C     ======================================
                1057 C     Calculate the eddy induced overturning streamfunction
                1058 C     ======================================
05118ac017 Jean*1059 #ifdef GM_BATES_PASSIVE
0d1e4b948d Mich*1060       k=Nr
4578baf364 Jean*1061       DO j=1-OLy,sNy+OLy
                1062        DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*1063         psistar(i,j,Nr) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
                1064        ENDDO
                1065       ENDDO
                1066       DO k=Nr-1,1,-1
4578baf364 Jean*1067        DO j=1-OLy,sNy+OLy
                1068         DO i=1-OLx,sNx+OLx
0d1e4b948d Mich*1069            psistar(i,j,k) = psistar(i,j,k+1)
                1070      &          - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
                1071         ENDDO
                1072        ENDDO
                1073       ENDDO
0ff5c71d8d Jean*1074
0d1e4b948d Mich*1075 #else
                1076
                1077       IF (GM_AdvForm) THEN
                1078         k=Nr
4578baf364 Jean*1079         DO j=1-OLy+1,sNy+1
                1080          DO i=1-OLx+1,sNx+1
0d1e4b948d Mich*1081           GM_PsiX(i,j,k,bi,bj) = -hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
                1082           GM_PsiY(i,j,k,bi,bj) = -hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
                1083          ENDDO
                1084         ENDDO
                1085         DO k=Nr-1,1,-1
4578baf364 Jean*1086          DO j=1-OLy+1,sNy+1
                1087           DO i=1-OLx+1,sNx+1
0d1e4b948d Mich*1088            GM_PsiX(i,j,k,bi,bj) = GM_PsiX(i,j,k+1,bi,bj)
                1089      &          - hfacW(i,j,k,bi,bj)*drF(k)*ustar(i,j,k)
                1090            GM_PsiY(i,j,k,bi,bj) = GM_PsiY(i,j,k+1,bi,bj)
                1091      &          - hfacS(i,j,k,bi,bj)*drF(k)*vstar(i,j,k)
                1092           ENDDO
                1093          ENDDO
                1094         ENDDO
                1095
                1096       ENDIF
                1097 #endif
                1098
                1099 #ifdef ALLOW_DIAGNOSTICS
                1100 C     Diagnostics
                1101       IF ( useDiagnostics ) THEN
05118ac017 Jean*1102         CALL DIAGNOSTICS_FILL( GM_BatesK3d, 'GM_BaK  ',
                1103      I                         0, Nr, 1, bi, bj, myThid )
0ff5c71d8d Jean*1104         CALL DIAGNOSTICS_FILL(KPV,    'GM_KPV  ',0,Nr,2,bi,bj,myThid)
                1105         CALL DIAGNOSTICS_FILL(urms,   'GM_URMS ',0,Nr,2,bi,bj,myThid)
                1106         CALL DIAGNOSTICS_FILL(Rdef,   'GM_RDEF ',0, 1,1,bi,bj,myThid)
                1107         CALL DIAGNOSTICS_FILL(Rurms,  'GM_RURMS',0, 1,2,bi,bj,myThid)
                1108         CALL DIAGNOSTICS_FILL(RRhines,'GM_RRHNS',0, 1,2,bi,bj,myThid)
                1109         CALL DIAGNOSTICS_FILL(Rmix,   'GM_RMIX ',0, 1,2,bi,bj,myThid)
                1110         CALL DIAGNOSTICS_FILL(supp,   'GM_SUPP ',0,Nr,2,bi,bj,myThid)
                1111         CALL DIAGNOSTICS_FILL(Xix,    'GM_Xix  ',0,Nr,2,bi,bj,myThid)
                1112         CALL DIAGNOSTICS_FILL(Xiy,    'GM_Xiy  ',0,Nr,2,bi,bj,myThid)
                1113         CALL DIAGNOSTICS_FILL(cDopp,  'GM_C    ',0, 1,2,bi,bj,myThid)
                1114         CALL DIAGNOSTICS_FILL(Ubaro,  'GM_UBARO',0, 1,2,bi,bj,myThid)
                1115         CALL DIAGNOSTICS_FILL(eady,   'GM_EADY ',0, 1,2,bi,bj,myThid)
                1116         CALL DIAGNOSTICS_FILL(SlopeX, 'GM_Sx   ',0,Nr,2,bi,bj,myThid)
                1117         CALL DIAGNOSTICS_FILL(SlopeY, 'GM_Sy   ',0,Nr,2,bi,bj,myThid)
                1118         CALL DIAGNOSTICS_FILL(tfluxX, 'GM_TFLXX',0,Nr,2,bi,bj,myThid)
                1119         CALL DIAGNOSTICS_FILL(tfluxY, 'GM_TFLXY',0,Nr,2,bi,bj,myThid)
                1120         CALL DIAGNOSTICS_FILL(gradqx, 'GM_dqdx ',0,Nr,2,bi,bj,myThid)
                1121         CALL DIAGNOSTICS_FILL(gradqy, 'GM_dqdy ',0,Nr,2,bi,bj,myThid)
                1122         CALL DIAGNOSTICS_FILL(Kdqdy,  'GM_Kdqdy',0,Nr,2,bi,bj,myThid)
                1123         CALL DIAGNOSTICS_FILL(Kdqdx,  'GM_Kdqdx',0,Nr,2,bi,bj,myThid)
                1124         CALL DIAGNOSTICS_FILL(surfkz, 'GM_SFLYR',0, 1,2,bi,bj,myThid)
                1125         CALL DIAGNOSTICS_FILL(ustar,  'GM_USTAR',0,Nr,2,bi,bj,myThid)
                1126         CALL DIAGNOSTICS_FILL(vstar,  'GM_VSTAR',0,Nr,2,bi,bj,myThid)
                1127         CALL DIAGNOSTICS_FILL(umc,    'GM_UMC  ',0,Nr,2,bi,bj,myThid)
                1128         CALL DIAGNOSTICS_FILL(ubar,   'GM_UBAR ',0,Nr,2,bi,bj,myThid)
                1129         CALL DIAGNOSTICS_FILL(modesC, 'GM_MODEC',0,Nr,1,bi,bj,myThid)
                1130         CALL DIAGNOSTICS_FILL(M4loc,  'GM_M4   ',0,Nr,2,bi,bj,myThid)
                1131         CALL DIAGNOSTICS_FILL(N2loc,  'GM_N2   ',0,Nr,2,bi,bj,myThid)
                1132         CALL DIAGNOSTICS_FILL(M4onN2, 'GM_M4_N2',0,Nr,2,bi,bj,myThid)
                1133         CALL DIAGNOSTICS_FILL(slopeC, 'GM_SLOPE',0,Nr,2,bi,bj,myThid)
                1134         CALL DIAGNOSTICS_FILL(Renorm, 'GM_RENRM',0, 1,2,bi,bj,myThid)
05118ac017 Jean*1135 #ifdef GM_BATES_PASSIVE
74e46366ba Davi*1136         CALL DIAGNOSTICS_FILL(psistar,'GM_PSTAR',0,Nr,2,bi,bj,myThid)
                1137 #endif
0d1e4b948d Mich*1138       ENDIF
                1139 #endif
                1140
05118ac017 Jean*1141 C     For the Redi diffusivity, we set GM_BatesK3d to a constant if
                1142 C     GM_Bates_constRedi=.TRUE.
                1143       IF (GM_Bates_constRedi) THEN
952f8e8fd1 Mich*1144         DO k=1,Nr
4578baf364 Jean*1145          DO j=1-OLy,sNy+OLy
                1146           DO i=1-OLx,sNx+OLx
05118ac017 Jean*1147            GM_BatesK3d(i,j,k,bi,bj) = GM_Bates_constK
952f8e8fd1 Mich*1148           ENDDO
                1149          ENDDO
                1150         ENDDO
                1151       ENDIF
                1152
40312237e2 Mich*1153 #ifdef ALLOW_DIAGNOSTICS
0ff5c71d8d Jean*1154       IF ( useDiagnostics )
05118ac017 Jean*1155      &  CALL DIAGNOSTICS_FILL( GM_BatesK3d, 'GM_BaK_T',
                1156      I                         0, Nr, 1, bi, bj, myThid )
40312237e2 Mich*1157 #endif
                1158
05118ac017 Jean*1159 #endif /* GM_BATES_K3D */
0d1e4b948d Mich*1160       RETURN
                1161       END