pkg/gmredi/gmredi_calc_geom.F

a4576c7cde Juli*0001 #include "GMREDI_OPTIONS.h"
7bb5a8a109 Jean*0002 #ifdef ALLOW_GENERIC_ADVDIFF
                0003 # include "GAD_OPTIONS.h"
                0004 #endif
a4576c7cde Juli*0005
                0006 CBOP
                0007 C !ROUTINE: GMREDI_CALC_GEOM
                0008
                0009 C !INTERFACE: ==========================================================
                0010       SUBROUTINE GMREDI_CALC_GEOM(
                0011      I             sigmaX, sigmaY, sigmaR,
                0012      I             bi, bj, myTime, myIter, myThid )
                0013
                0014 C     !DESCRIPTION:
                0015 C     *==========================================================*
                0016 C     | SUBROUTINE GMREDI_CALC_GEOM
                0017 C     | Calculate GM coefficient with the GEOMETRIC prescription
                0018 C     | GEOM_K3d is located at the grid-cell vertical interface.
                0019 C     *==========================================================*
                0020
                0021 C !USES: ===============================================================
                0022       IMPLICIT NONE
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "DYNVARS.h"
                0028 #include "GMREDI.h"
                0029
                0030 C !INPUT PARAMETERS: ===================================================
                0031 C     sigmaXYR  :: density gradient variables
                0032 C     bi, bj    :: tile indices
                0033 C     myTime    :: Current time in simulation
                0034 C     myIter    :: Current iteration number in simulation
                0035 C     myThid    :: My Thread Id. number
                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)
                0039       INTEGER bi, bj
                0040       _RL     myTime
                0041       INTEGER myIter
                0042       INTEGER myThid
                0043
                0044 #ifdef GM_GEOM_VARIABLE_K
                0045 C !INOUT PARAMETERS: ===================================================
                0046 C  GEOM_K3d     :: GEOMETRIC K-GM (GEOMETRIC in horizontal,
                0047 C                             structure function in vertical)
                0048 C  GEOM_EKE     :: GEOMETRIC parameterised energy
                0049 C              not really an output, passed around as common block
                0050
                0051 C !LOCAL VARIABLES: ====================================================
                0052 C     i,j,k           :: Loop counters
                0053 C     Slope(XY)       :: isopycnal slopes
                0054 C     dSigmaD(xyr)    :: derivative of density variables
                0055 C     dSigma(HR)      :: variables to compute isopycnal slopes
                0056 C     (MN)2loc        :: local M^2 and N^2
                0057 C     S(N)loc         :: local M^2/N^2 and M^2/N
                0058 C     (S)Nloc_zint    :: depth integrated M^2/N and N
                0059 C                       (for calculating trd_ene_gen and trd_ene_wav)
                0060 C     trd_ene_*       :: trends for energy
                0061 C                        gen, adv, lap, dis
                0062 C     ene_*           :: intermediate variables for computing trends
                0063 C     depth*          :: various depth variables for computation purposes
                0064 C     c1, c_ros[XYEN] :: long Rossby wave phase speeds for trd_ene_wav
                0065 C     UV_depth_avg    :: depth-avg flow for                trd_ene_adv
                0066 C     ab0, ab1        :: Adams-Bashforth weights
                0067
                0068       INTEGER i,j,k
                0069       INTEGER kSurf
                0070       _RL dSigmaDx
                0071       _RL dSigmaDy
                0072       _RL dSigmaDr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0073
                0074       _RL dSigmaH, dSigmaR
                0075       _RL N2loc, SNloc, Sloc
                0076       _RL SNloc3D(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0077       _RL N2loc3D(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0078
                0079       _RL recipMaxSlope
                0080       _RL c_dum, fp_im, fm_ip
                0081       _RL ab0, ab1, fCoriFac
                0082
                0083 C     metric variables
                0084       _RL recip_depthW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL recip_depthS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086 C     depth at C-points in meters
                0087       _RL depthC      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088
                0089 C     hFac of the W-control volumes
                0090       _RL hFac
                0091
                0092 C     general variables
                0093       _RL deltaTloc
                0094
                0095 C     local version of current rhs of eddy energy equation
                0096       _RL ene_rhs_now(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0097       _RL vert_struc_func(1-OLx:sNx+OLx,1-OLy:sNy+OLy, Nr)
                0098
                0099       _RL SNloc_zint (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0100       _RL  Nloc_zint (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0101       _RL trd_ene_gen(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0102       _RL trd_ene_dis(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0103       _RL trd_ene_adv(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0104       _RL trd_ene_wav(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL trd_ene_lap(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106
                0107       _RL U_depth_avg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL V_depth_avg(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109
                0110       _RL c1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111       _RL c_rosX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL c_rosY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL c_rosE
                0114
                0115       _RL ene_adv_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL ene_adv_y(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL ene_wav_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0118       _RL ene_wav_y(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0119       _RL ene_lap_x(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0120       _RL ene_lap_y(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0121
                0122 CEOP
                0123
                0124 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0125
                0126       kSurf = 1
                0127       IF ( usingPCoords ) kSurf = Nr
                0128       deltaTloc = dtTracerLev(kSurf)
                0129       recipMaxSlope = 0. _d 0
                0130       IF ( GM_maxSlope.GT.0. _d 0 ) THEN
                0131         recipMaxSlope = 1. _d 0 / GM_maxSlope
                0132       ENDIF
                0133
                0134 C--   initialise some variables for calculations
                0135
                0136 C     depths for doing depth-averaging
                0137       DO j=1-OLy,sNy+OLy
                0138        DO i=1-OLx,sNx+OLx
                0139         recip_depthW(i,j) = 0.0 _d 0
                0140         recip_depthS(i,j) = 0.0 _d 0
                0141         c_dum = rSurfW(i,j,bi,bj) - rLowW(i,j,bi,bj)
                0142         IF ( c_dum .GT. zeroRL ) recip_depthW(i,j) = 1. _d 0 / c_dum
                0143         c_dum = rSurfS(i,j,bi,bj) - rLowS(i,j,bi,bj)
                0144         IF ( c_dum .GT. zeroRL ) recip_depthS(i,j) = 1. _d 0 / c_dum
                0145 C     convert to meters in the case of p-coordinates
                0146         depthC(i,j) = ( Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) )
                0147      &              * rUnit2mass*recip_rhoConst
                0148        ENDDO
                0149       ENDDO
                0150
                0151 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0152 C   TODO: should the calculations be based on the tapered slopes?
                0153 C         here it is just the raw stuff
                0154 C         if being moved then this routine might want to be after
                0155 C         the SLOPE_LIMIT routine
                0156
                0157 C     initialisations, put this in a loop if need be
                0158       DO j=1-OLy,sNy+OLy
                0159        DO i=1-OLx,sNx+OLx
                0160         SNloc_zint(i,j) = 0.0 _d 0
                0161         Nloc_zint(i,j)  = 0.0 _d 0
                0162        ENDDO
                0163       ENDDO
                0164       DO k=1,Nr
                0165        DO j=1-OLy,sNy+OLy
                0166         DO i=1-OLx,sNx+OLx
                0167          vert_struc_func(i,j,k) = 1.0 _d 0
                0168          N2loc3D(i,j,k)         = 0.0 _d 0
                0169          SNloc3D(i,j,k)         = 0.0 _d 0
                0170         ENDDO
                0171        ENDDO
                0172       ENDDO
                0173 C-- 1st k loop : compute vertical structure to be used later
                0174       DO k=Nr,2,-1
                0175
                0176 C     For stable conditions sigmaR<0 for z-coordinates, but >0 for p-coords
                0177 C     => change sign of vertical Sigma gradient (always downward)
                0178 C     to match stratification sign (i.e., positive if stratified)
                0179 C     For p-coordinates: convert r-units to z
                0180        DO j=1-OLy,sNy+OLy
                0181         DO i=1-OLx,sNx+OLx
                0182          dSigmaDr(i,j) = MAX( gravitySign*sigmaR(i,j,k), zeroRL )
                0183      &        * mass2rUnit * rhoConst
                0184         ENDDO
                0185        ENDDO
                0186
                0187 C      NOTE: Ignores boundary cells (convenient with the 5-point filter used)
                0188        DO j=2-OLy,sNy+OLy-1
                0189         DO i=2-OLx,sNx+OLx-1
                0190          IF ( maskC(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
                0191
                0192 C      Compute -N^2 rho_0 / g via averaging, on vertical interfaces
                0193           dSigmaR = (dSigmaDr(i,j) * 4.0 _d 0
                0194      &                + maskC(i-1,j,k,bi,bj)*dSigmaDr(i-1,j)
                0195      &                + maskC(i+1,j,k,bi,bj)*dSigmaDr(i+1,j)
                0196      &                + maskC(i,j-1,k,bi,bj)*dSigmaDr(i,j-1)
                0197      &                + maskC(i,j+1,k,bi,bj)*dSigmaDr(i,j+1)
                0198      &                 ) / (4.0 _d 0
                0199      &                    + maskC(i-1,j,k,bi,bj)
                0200      &                    + maskC(i+1,j,k,bi,bj)
                0201      &                    + maskC(i,j-1,k,bi,bj)
                0202      &                    + maskC(i,j+1,k,bi,bj)
                0203      &                     )
                0204 C      Compute M^2 rho_0 / g on vertical interfaces
                0205           dSigmaDx = op25 * (sigmaX(i+1,j,k-1) + sigmaX(i,j,k-1)
                0206      &                    +  sigmaX(i+1,j,k  ) + sigmaX(i,j,k  )
                0207      &                       ) *  maskC(i,j,k,bi,bj)
                0208           dSigmaDy = op25 * (sigmaY(i,j+1,k-1) + sigmaY(i,j,k-1)
                0209      &                    +  sigmaY(i,j+1,k  ) + sigmaY(i,j,k  )
                0210      &                       ) *  maskC(i,j,k,bi,bj)
                0211           dSigmaH = SQRT(dSigmaDx * dSigmaDx
                0212      &                 + dSigmaDy * dSigmaDy)
                0213           IF ( dSigmaH .GT. 0. _d 0 ) THEN
                0214            IF ( dSigmaR .GT. dSigmaH*recipMaxSlope ) THEN
                0215             Sloc = dSigmaH / dSigmaR
                0216            ELSE
                0217             Sloc = GM_maxSlope
                0218            ENDIF
                0219            N2loc3D(i,j,k) = gravity * recip_rhoConst * dSigmaR
                0220            SNloc3D(i,j,k) = Sloc * SQRT(N2loc3D(i,j,k))
                0221           ENDIF
                0222          ENDIF
                0223         ENDDO
                0224        ENDDO
                0225       ENDDO
                0226       IF ( GEOM_vert_struc ) THEN
                0227        k=2
                0228 C      avoid division by zero
                0229        DO j=1-OLy,sNy+OLy
                0230         DO i=1-OLx,sNx+OLx
                0231          vert_struc_func(i,j,k) = MAX(N2loc3D(i,j,k),
                0232      &                                GEOM_vert_struc_min)
                0233         ENDDO
                0234        ENDDO
                0235 C      Cap the (N^2 / N^2_surf) between something (1 and 0.1 default)
                0236        DO k=Nr,2,-1
                0237         DO j=1-OLy,sNy+OLy
                0238          DO i=1-OLx,sNx+OLx
                0239           vert_struc_func(i,j,k) =
                0240      &         MAX(MIN(GEOM_vert_struc_max,
                0241      &                 N2loc3D(i,j,k)/vert_struc_func(i,j,2)),
                0242      &             GEOM_vert_struc_min)
                0243          ENDDO
                0244         ENDDO
                0245        ENDDO
                0246       ENDIF
                0247
                0248 C-- 2nd k loop : compute the rest of the GEOMETRIC stuff
                0249       DO k=Nr,2,-1
                0250
                0251 C      NOTE: Ignores boundary cells (convenient with the 5-point filter used)
                0252        DO j=2-OLy,sNy+OLy-1
                0253         DO i=2-OLx,sNx+OLx-1
                0254          IF ( maskC(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
                0255           N2loc = N2loc3D(i,j,k)
                0256           SNloc = SNloc3D(i,j,k)
                0257 C     hFac now contains the coordinate factor rUnit2mass*recip_rhoConst
                0258 C     ( = 1 for z-coordinates and 1/(gravity*rhoConst for p-coordinates)
                0259           hFac = MIN( halfRS, _hFacC(i,j,k-1,bi,bj) )
                0260      &         + MIN( halfRS, _hFacC(i,j,k  ,bi,bj) )
                0261      &                    * rUnit2mass*recip_rhoConst
                0262           SNloc_zint(i,j) = SNloc_zint(i,j)
                0263      &                    + SNloc * drC(k) * hFac
                0264      &                    * vert_struc_func(i,j,k)
                0265           Nloc_zint(i,j)  = Nloc_zint(i,j) + SQRT(N2loc)
                0266      &                    * drC(k) * hFac
                0267          ENDIF
                0268         ENDDO
                0269        ENDDO
                0270       ENDDO
                0271
                0272 C     work out the implied long Rossby phase speeds at T pts
                0273 C     average it onto UV points later
                0274
                0275       fCoriFac = 0.0 _d 0
                0276       IF ( usingCartesianGrid .AND. f0 .NE. 0.0 _d 0 )
                0277      &     fCoriFac = beta / (f0 * f0)
                0278
                0279       DO j=1-OLy,sNy+OLy
                0280        DO i=1-OLx,sNx+OLx
                0281 C       only compute the Rossby phase speed if deep enough
                0282         IF (depthC(i,j) .GT. 300.0 _d 0) THEN
                0283          c1(i,j) = MIN(10.0 _d 0, Nloc_zint(i,j) / PI)
                0284      &              * maskC(i,j,kSurf,bi,bj)
                0285
                0286 C       bound the f factors away from zero (inside Eq +/- 5degN/S),
                0287 C       compute the Rossby phase speed, and bound above by
                0288 C       1st baroclinic equatorial Rossby wave phase speed
                0289          IF ( .NOT.usingCartesianGrid ) THEN
                0290           fCoriFac = MAX( ABS(fCori(i,j,bi,bj)), 1.2676 _d -05 )
                0291           fCoriFac = fCoriCos(i,j,bi,bj) / ( rSphere * fCoriFac ** 2 )
                0292          ENDIF
                0293          c_rosE = - MIN( c1(i,j)/3.0 _d 0, c1(i,j)*c1(i,j)*fCoriFac )
                0294 C     rotating EW Rossby velocity to XY for the grid, c_rosN=0 by definition
                0295 C     (trivial rotation for Cartesian plane)
                0296          c_rosX(i,j) =  angleCosC(i,j,bi,bj)*c_rosE
                0297 c    &                + angleSinC(i,j,bi,bj)*c_rosN
                0298          c_rosY(i,j) = -angleSinC(i,j,bi,bj)*c_rosE
                0299 c    &                 +angleCosC(i,j,bi,bj)*c_rosN
                0300         ELSE
                0301          c_rosX(i,j) = 0.0 _d 0
                0302          c_rosY(i,j) = 0.0 _d 0
                0303          c1    (i,j) = 0.0 _d 0
                0304         ENDIF
                0305        ENDDO
                0306       ENDDO
                0307
                0308 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0309 C-- compute kgm as per GEOMETRIC in horizontal, then extend
                0310 C-- vertically by structure function if required
                0311 C-- do other tapering things here
                0312
                0313 C     1. compute the 2d k_GEOM first
                0314 C        bound the denominator from below (choice following NEMO)
                0315
                0316       DO j=1-OLy,sNy+OLy
                0317        DO i=1-OLx,sNx+OLx
                0318 C     for cubed sphere setups, depthC can be zero in the halo, so we
                0319 C     include it in the masking
                0320         IF ( maskC(i,j,kSurf,bi,bj)*depthC(i,j) .GT. 0. _d 0 ) THEN
                0321          GEOM_K3d(i,j,kSurf,bi,bj) = GEOM_alpha * GEOM_EKE(i,j,bi,bj)
                0322      &        / MAX( SNloc_zint(i,j), 1. _d -7 * depthC(i,j) )
                0323         ENDIF
                0324        ENDDO
                0325       ENDDO
                0326
                0327 C     2. cap k_GEOM from above! (stop it being too big)
                0328
                0329       IF ( GEOM_alpha .NE. 0.0 _d 0 ) THEN
                0330        DO j=1-OLy,sNy+OLy
                0331         DO i=1-OLx,sNx+OLx
                0332          GEOM_K3d(i,j,kSurf,bi,bj) =
                0333      &        MIN( GEOM_K3d(i,j,kSurf,bi,bj), GEOM_maxVal_K )
                0334         ENDDO
                0335        ENDDO
                0336       ENDIF
                0337
                0338 C     3. taper it according to depth to kill k_GEOM in shallow regions
                0339
                0340       DO j=1-OLy,sNy+OLy
                0341        DO i=1-OLx,sNx+OLx
                0342         GEOM_K3d(i,j,kSurf,bi,bj) = GEOM_K3d(i,j,kSurf,bi,bj)
                0343      &       * GEOM_taper(i,j,bi,bj) * maskC(i,j,kSurf,bi,bj)
                0344        ENDDO
                0345       ENDDO
                0346
                0347 C     4. extend the (tapered) k_GEOM in depth
                0348
                0349       DO k=Nr,2,-1
                0350        DO j=1-OLy,sNy+OLy
                0351         DO i=1-OLx,sNx+OLx
                0352          GEOM_K3d(i,j,k,bi,bj) = vert_struc_func(i,j,k)
                0353      &                      * GEOM_K3d(i,j,kSurf,bi,bj)
                0354      &                      * maskC(i,j,k,bi,bj)
                0355         ENDDO
                0356        ENDDO
                0357       ENDDO
                0358
                0359 C     5. cap k_GEOM from below (do not strictly need this)
                0360
                0361       IF ( GEOM_alpha .NE. 0.0 _d 0 ) THEN
                0362        DO k=1,Nr
                0363         DO j=1-OLy,sNy+OLy
                0364          DO i=1-OLx,sNx+OLx
                0365           GEOM_K3d(i,j,k,bi,bj) =
                0366      &         MAX( GEOM_K3d(i,j,k,bi,bj), GEOM_minVal_K )
                0367          ENDDO
                0368         ENDDO
                0369        ENDDO
                0370       ENDIF
                0371
                0372 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0373 C--   time step the energy
                0374 C     allocate and initialise
                0375       DO j=1-OLy,sNy+OLy
                0376        DO i=1-OLx,sNx+OLx
                0377         trd_ene_gen(i,j) = 0.0 _d 0
                0378        ENDDO
                0379       ENDDO
                0380
                0381 C     loop over k
                0382       DO k=Nr,2,-1
                0383
                0384 C      NOTE: Ignores boundary cells (convenient with the 5-point filter used)
                0385        DO j=2-OLy,sNy+OLy-1
                0386         DO i=2-OLx,sNx+OLx-1
                0387          IF ( maskC(i,j,k,bi,bj) .GT. 0. _d 0 ) THEN
                0388           hFac = MIN( halfRS, _hFacC(i,j,k-1,bi,bj) )
                0389      &         + MIN( halfRS, _hFacC(i,j,k  ,bi,bj) )
                0390           SNloc = SNloc3D(i,j,k)
                0391           trd_ene_gen(i,j) = trd_ene_gen(i,j)
                0392      &                     + GEOM_K3d(i,j,k,bi,bj)
                0393      &                     * SNloc * SNloc
                0394      &                     * drC(k) * hFac
                0395      &                     * rUnit2mass*recip_rhoConst
                0396          ENDIF
                0397         ENDDO
                0398        ENDDO
                0399       ENDDO
                0400 C ----------------------------------------------
                0401 C ---- start advection of energy ---------------
                0402 C ----------------------------------------------
                0403
                0404       DO j=1-OLy,sNy+OLy
                0405        DO i=1-OLx,sNx+OLx
                0406         trd_ene_dis(i,j) = 0.0 _d 0
                0407         trd_ene_adv(i,j) = 0.0 _d 0
                0408         trd_ene_wav(i,j) = 0.0 _d 0
                0409         trd_ene_lap(i,j) = 0.0 _d 0
                0410
                0411         U_depth_avg(i,j) = 0.0 _d 0
                0412         V_depth_avg(i,j) = 0.0 _d 0
                0413
                0414         ene_adv_x(i,j)   = 0.0 _d 0
                0415         ene_adv_y(i,j)   = 0.0 _d 0
                0416         ene_wav_x(i,j)   = 0.0 _d 0
                0417         ene_wav_y(i,j)   = 0.0 _d 0
                0418         ene_lap_x(i,j)   = 0.0 _d 0
                0419         ene_lap_y(i,j)   = 0.0 _d 0
                0420         ene_rhs_now(i,j) = 0.0 _d 0
                0421        ENDDO
                0422       ENDDO
                0423
                0424 C ----  form the energy fluxes
                0425
                0426 C     1. form the depth-AVG velocity
                0427       DO k= 1,Nr
                0428        DO j=1-OLy,sNy+OLy
                0429         DO i=1-OLx,sNx+OLx
                0430          U_depth_avg(i,j) = U_depth_avg(i,j)
                0431      &                    + uVel(i,j,k,bi,bj) * drF(k)
                0432      &                    * hFacW(i,j,k,bi,bj) * recip_depthW(i,j)
                0433          V_depth_avg(i,j) = V_depth_avg(i,j)
                0434      &                    + vVel(i,j,k,bi,bj) * drF(k)
                0435      &                    * hFacS(i,j,k,bi,bj) * recip_depthS(i,j)
                0436         ENDDO
                0437        ENDDO
                0438       ENDDO
                0439
                0440 C     2. form the advective tendency in each direction
                0441       DO j=2-OLy,sNy+OLy-1
                0442        DO i=2-OLx,sNx+OLx-1
                0443 C-    X direction:
                0444 C ---- second ordered centred difference
                0445 c       ene_adv_x(i,j) = -halfRL
                0446 c    &    *( U_depth_avg(i,j) * dyG(i,j,bi,bj)
                0447 c    &       *( GEOM_EKE( i ,j,bi,bj) - GEOM_EKE(i-1,j,bi,bj) )
                0448 c    &     + U_depth_avg(i+1,j) * dyG(i+1,j,bi,bj)
                0449 c    &       *( GEOM_EKE(i+1,j,bi,bj) - GEOM_EKE( i ,j,bi,bj) )
                0450 c    &     )
                0451 C ---- 1st order upwinding
                0452         fp_im = U_depth_avg( i ,j) + ABS( U_depth_avg( i ,j) )
                0453         fm_ip = U_depth_avg(i+1,j) - ABS( U_depth_avg(i+1,j) )
                0454         ene_adv_x(i,j) = -halfRL
                0455      &    *( fp_im * dyG(i,j,bi,bj)
                0456      &       *( GEOM_EKE( i ,j,bi,bj) - GEOM_EKE(i-1,j,bi,bj) )
                0457      &     + fm_ip * dyG(i+1,j,bi,bj)
                0458      &       *( GEOM_EKE(i+1,j,bi,bj) - GEOM_EKE( i ,j,bi,bj) )
                0459      &     )
                0460 C ---- 1st order upwinding wave stuff
                0461         c_dum = ( c_rosX(i-1,j) + c_rosX(i,j) )*halfRL
                0462         fp_im = ( c_dum + ABS( c_dum ) )*maskW( i ,j,kSurf,bi,bj)
                0463         c_dum = ( c_rosX(i,j) + c_rosX(i+1,j) )
                0464         fm_ip = ( c_dum - ABS( c_dum ) )*maskW(i+1,j,kSurf,bi,bj)
                0465         ene_wav_x(i,j) = -halfRL
                0466      &    *( fp_im * dyG(i,j,bi,bj)
                0467      &       *( GEOM_EKE( i ,j,bi,bj) - GEOM_EKE(i-1,j,bi,bj) )
                0468      &     + fm_ip * dyG(i+1,j,bi,bj)
                0469      &       *( GEOM_EKE(i+1,j,bi,bj) - GEOM_EKE( i ,j,bi,bj) )
                0470      &     )
                0471 C-    Y direction:
                0472 C ---- second ordered centred difference
                0473 c       ene_adv_y(i,j) = -halfRL
                0474 c    &    *( V_depth_avg(i,j) * dxG(i,j,bi,bj)
                0475 c    &       *( GEOM_EKE(i, j ,bi,bj) - GEOM_EKE(i,j-1,bi,bj) )
                0476 c    &     + V_depth_avg(i,j+1) * dxG(i,j+1,bi,bj)
                0477 c    &       *( GEOM_EKE(i,j+1,bi,bj) - GEOM_EKE(i, j ,bi,bj) )
                0478 c    &     )
                0479 C ---- 1st order upwinding
                0480         fp_im = V_depth_avg(i, j ) + ABS( V_depth_avg(i, j ) )
                0481         fm_ip = V_depth_avg(i,j+1) - ABS( V_depth_avg(i,j+1) )
                0482         ene_adv_y(i,j) = -halfRL
                0483      &    *( fp_im * dxG(i,j,bi,bj)
                0484      &       *( GEOM_EKE(i, j ,bi,bj) - GEOM_EKE(i,j-1,bi,bj) )
                0485      &     + fm_ip * dxG(i,j+1,bi,bj)
                0486      &       *( GEOM_EKE(i,j+1,bi,bj) - GEOM_EKE(i, j ,bi,bj) )
                0487      &     )
                0488 C ---- 1st order upwinding wave stuff
                0489         c_dum = ( c_rosY(i-1,j) + c_rosY(i,j) )*halfRL
                0490         fp_im = ( c_dum + ABS( c_dum ) )*maskS(i, j ,kSurf,bi,bj)
                0491         c_dum = ( c_rosY(i,j) + c_rosY(i,j+1) )
                0492         fm_ip = ( c_dum - ABS( c_dum ) )*maskS(i,j+1,kSurf,bi,bj)
                0493         ene_wav_y(i,j) = -halfRL
                0494      &    *( fp_im * dxG(i,j,bi,bj)
                0495      &       *( GEOM_EKE(i, j ,bi,bj) - GEOM_EKE(i,j-1,bi,bj) )
                0496      &     + fm_ip * dxG(i,j+1,bi,bj)
                0497      &       *( GEOM_EKE(i,j+1,bi,bj) - GEOM_EKE(i, j ,bi,bj) )
                0498      &     )
                0499        ENDDO
                0500       ENDDO
                0501
                0502 C     3. form the diffusive fluxes in each direction
                0503       DO j=1-OLy,sNy+OLy
                0504        DO i=2-OLx,sNx+OLx
                0505         ene_lap_x(i,j) = -GEOM_diffKh_EKE * dyG(i,j,bi,bj)
                0506      &       * ( GEOM_EKE(i,j,bi,bj) - GEOM_EKE(i-1,j,bi,bj) )
                0507      &       * recip_dxC(i,j,bi,bj)*maskW(i,j,kSurf,bi,bj)
                0508      &       * cosFacU(j,bi,bj)
                0509        ENDDO
                0510       ENDDO
                0511       DO j=2-OLy,sNy+OLy
                0512        DO i=1-OLx,sNx+OLx
                0513         ene_lap_y(i,j) = -GEOM_diffKh_EKE * dxG(i,j,bi,bj)
                0514      &       * ( GEOM_EKE(i,j,bi,bj) - GEOM_EKE(i,j-1,bi,bj) )
                0515      &       * recip_dyC(i,j,bi,bj)*maskS(i,j,kSurf,bi,bj)
                0516 #ifdef ISOTROPIC_COS_SCALING
                0517      &       *cosFacV(j,bi,bj)
                0518 #endif
                0519        ENDDO
                0520       ENDDO
                0521
                0522 C     4. sum or form tendencies on cell
                0523       DO j=2-OLy,sNy+OLy-1
                0524        DO i=2-OLx,sNx+OLx-1
                0525         trd_ene_dis(i,j) = - GEOM_lmbda * GEOM_EKE(i,j,bi,bj)
                0526      &                                          *maskInC(i,j,bi,bj)
                0527 C     Add advective contribution and divide by grid cell area
                0528         trd_ene_adv(i,j) = ( ene_adv_x(i,j) + ene_adv_y(i,j) )
                0529      &                     * recip_rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0530         trd_ene_wav(i,j) = ( ene_wav_x(i,j) + ene_wav_y(i,j) )
                0531      &                     * recip_rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0532 C     Form diffusive tendency from (minus) flux divergence
                0533         trd_ene_lap(i,j) = -( ene_lap_x(i+1,j) - ene_lap_x(i,j)
                0534      &                      + ene_lap_y(i,j+1) - ene_lap_y(i,j) )
                0535      &                      *recip_rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
                0536 C       form the RHS
                0537         ene_rhs_now(i,j) = trd_ene_gen(i,j)
                0538      &                   + trd_ene_dis(i,j)
                0539      &                   + trd_ene_adv(i,j)
                0540      &                   + trd_ene_wav(i,j)
                0541      &                   + trd_ene_lap(i,j)
                0542        ENDDO
                0543       ENDDO
                0544
                0545 C     At this point, the new ene_rhs_now is defined and correctly
                0546 C     computed for all i,j except for 1-OLx/y and sNx/y+Olx/y, and we do
                0547 C     not need any exchange. We do need to exchange the solution GEOM_EKE
                0548 C     in S/R GMREDI_DO_EXCH.
                0549
                0550 C     third: time stepping
                0551       IF ( GEOM_startAB .EQ. 0 ) THEN
                0552        ab0 =  1.0 _d 0
                0553        ab1 =  0.0 _d 0
                0554 C-    moved to gmredi_do_exch.F (outside bi,bj loop):
                0555 c       GEOM_startAB = 1
                0556       ELSE
                0557        ab0 =  1.5 _d 0 + abEps
                0558        ab1 = -0.5 _d 0 - abEps
                0559       ENDIF
                0560       DO j=2-OLy,sNy+OLy-1
                0561        DO i=2-OLx,sNx+OLx-1
                0562         GEOM_EKE(i,j,bi,bj) = GEOM_EKE(i,j,bi,bj)
                0563      &            + deltaTloc * (
                0564      &            + ab0 * ene_rhs_now(i,j)
                0565      &            + ab1 * GEOM_gEKE_Nm1(i,j,bi,bj)
                0566      &                        )
                0567        ENDDO
                0568       ENDDO
                0569 C ----------------------------------------------
                0570 C ---- end advection of energy -----------------
                0571 C ----------------------------------------------
                0572
                0573 C     update rhs fields for the next time step
                0574       DO j=1-OLy,sNy+OLy
                0575        DO i=1-OLx,sNx+OLx
                0576         GEOM_gEKE_Nm1(i,j,bi,bj) = ene_rhs_now(i,j)
                0577        ENDDO
                0578       ENDDO
                0579
                0580 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0581 C--   diagnostics and restarts
                0582
                0583 #ifdef ALLOW_DIAGNOSTICS
                0584       IF ( useDiagnostics ) THEN
                0585        CALL DIAGNOSTICS_FILL(vert_struc_func,
                0586      &                       'GEOMstru', 0,Nr,2,bi,bj,myThid)
                0587        CALL DIAGNOSTICS_FILL(trd_ene_gen,
                0588      &                       'GEOMEgen', 0,1,2,bi,bj,myThid)
                0589        CALL DIAGNOSTICS_FILL(trd_ene_dis,
                0590      &                       'GEOMEdis', 0,1,2,bi,bj,myThid)
                0591 C        advective trends are dE/dt + -(u - c) dot grad E
                0592        CALL DIAGNOSTICS_FILL(trd_ene_adv,
                0593      &                       'GEOMEadv', 0,1,2,bi,bj,myThid)
                0594        CALL DIAGNOSTICS_FILL(trd_ene_wav,
                0595      &                       'GEOMEwav', 0,1,2,bi,bj,myThid)
                0596        CALL DIAGNOSTICS_FILL(trd_ene_lap,
                0597      &                       'GEOMElap', 0,1,2,bi,bj,myThid)
                0598        CALL DIAGNOSTICS_FILL(c1,
                0599      &                       'GEOM_c1 ', 0,1,2,bi,bj,myThid)
                0600       ENDIF
                0601 #endif /* ALLOW_DIAGNOSTICS */
                0602
                0603 #endif /* GM_GEOM_VARIABLE_K */
                0604
                0605       RETURN
                0606       END