Back to home page

MITgcm
 
 

    

File indexing completed on 2024-02-29 06:10:23 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!