model/src/calc_oce_mxlayer.F

6aab75f3e8 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003
                0004 CBOP
                0005 C     !ROUTINE: CALC_OCE_MXLAYER
                0006 C     !INTERFACE:
                0007       SUBROUTINE CALC_OCE_MXLAYER(
                0008      I                       rhoSurf, sigmaR,
                0009      I                       bi, bj, myTime, myIter, myThid )
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
                0012 C     | S/R CALC_OCE_MXLAYER
                0013 C     | o Diagnose the Oceanic surface Mixed-Layer
0320e25227 Mart*0014 C     | Note: output "hMixLayer" is in "r" unit, i.e., in Pa
                0015 C     |       when using P-coordinate.
6aab75f3e8 Jean*0016 C     *==========================================================*
                0017 C     \ev
                0018
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021 C     == Global variables ==
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
                0026 #include "DYNVARS.h"
                0027 #ifdef ALLOW_GMREDI
                0028 # include "GMREDI.h"
                0029 #endif
                0030
                0031 C     !INPUT/OUTPUT PARAMETERS:
                0032 C     == Routine Arguments ==
                0033 C     rhoSurf   :: Surface density anomaly
                0034 C     sigmaR    :: Vertical gradient of potential density
                0035 C     bi,bj     :: tile indices
                0036 C     myTime    :: Current time in simulation
                0037 C     myIter    :: Current iteration number in simulation
                0038 C     myThid    :: my Thread Id number
                0039       _RL     rhoSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0040       _RL     sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0041       INTEGER bi, bj
                0042       _RL     myTime
                0043       INTEGER myIter
                0044       INTEGER myThid
                0045
                0046 C     === Functions ====
                0047 #ifdef ALLOW_DIAGNOSTICS
                0048       LOGICAL  DIAGNOSTICS_IS_ON
                0049       EXTERNAL DIAGNOSTICS_IS_ON
                0050 #endif /* ALLOW_DIAGNOSTICS */
                0051
                0052 C     !LOCAL VARIABLES:
                0053 C     == Local variables ==
0320e25227 Mart*0054 C     i, j, k   :: Loop counters
                0055 C     kSrf      :: surface center-level index
                0056 C     kTop      :: surface interface level index
                0057 C     kDir      :: downward k increment
                0058 C     kU, kL    :: loop ranges, Up and Low
                0059 C     kup, klw  :: interface level index above & below current level k
                0060       INTEGER i, j, k
                0061       INTEGER kSrf, kTop, kDir, deltaK
                0062       INTEGER kU, kL, klw, kup
                0063       LOGICAL calcMixLayerDepth, kIn
6aab75f3e8 Jean*0064       INTEGER method
7c01d59526 Davi*0065       _RL     rhoBigNb
6aab75f3e8 Jean*0066       _RL     rhoMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0067       _RL     rhoKm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL     rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL     tmpFac, sigmAv
                0070 CEOP
                0071
                0072       calcMixLayerDepth = .FALSE.
                0073 #ifdef ALLOW_GMREDI
050b4366e6 Jean*0074       IF ( useGMRedi .AND. .NOT.useKPP ) THEN
                0075        calcMixLayerDepth = GM_useSubMeso .OR. GM_taper_scheme.EQ.'fm07'
05118ac017 Jean*0076      &       .OR. GM_useBatesK3d
b61261c68c Davi*0077       ENDIF
6aab75f3e8 Jean*0078 #endif
                0079 #ifdef ALLOW_DIAGNOSTICS
                0080       IF ( useDiagnostics.AND. .NOT.calcMixLayerDepth ) THEN
                0081         calcMixLayerDepth = DIAGNOSTICS_IS_ON('MXLDEPTH',myThid)
                0082       ENDIF
                0083 #endif
                0084       IF ( calcMixLayerDepth ) THEN
                0085
0320e25227 Mart*0086         IF ( usingPCoords ) THEN
                0087          kTop    = Nr+1
                0088          kSrf    = Nr
                0089          kDir    = -1
                0090          deltaK  =  1
                0091         ELSE
                0092          kTop    =  1
                0093          kSrf    =  1
                0094          kDir    =  1
                0095          deltaK  =  0
                0096         ENDIF
6aab75f3e8 Jean*0097 C--   Select which "method" to use:
                0098        method = 0
                0099        IF ( hMixCriteria.LT.0. ) method = 1
                0100        IF ( hMixCriteria.GT.1. ) method = 2
                0101
                0102        IF ( method.EQ.1 ) THEN
                0103
                0104 C--   First method :
                0105 C     where the potential density (ref.lev=surface) is larger than
                0106 C       surface density plus Delta_Rho = hMixCriteria * Alpha(surf)
                0107 C     = density of water which is -hMixCriteria colder than surface water
eee0b35b17 Dimi*0108 C     (see Kara, Rochford, and Hurlburt JGR 2000 for default criterion)
6aab75f3e8 Jean*0109
                0110 c       hMixCriteria  = -0.8 _d 0
7c01d59526 Davi*0111 c       dRhoSmall = 1. _d -6
6aab75f3e8 Jean*0112         rhoBigNb  = rhoConst*1. _d 10
                0113         CALL FIND_ALPHA(
0320e25227 Mart*0114      I            bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, kSrf, kSrf,
6aab75f3e8 Jean*0115      O            rhoMxL, myThid )
                0116
0320e25227 Mart*0117         DO j=1-OLy,sNy+OLy
                0118          DO i=1-OLx,sNx+OLx
6aab75f3e8 Jean*0119            rhoKm1(i,j) = rhoSurf(i,j)
                0120            rhoMxL(i,j) = rhoSurf(i,j)
                0121      &                 + MAX( rhoMxL(i,j)*hMixCriteria, dRhoSmall )
0320e25227 Mart*0122            hMixLayer(i,j,bi,bj) = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
6aab75f3e8 Jean*0123          ENDDO
                0124         ENDDO
0320e25227 Mart*0125 C-    Z-coord: (kU,kL,kDir) = (2,Nr,1) ; P-coord: (kU,kL,kDir)= (Nr-1,1,-1)
                0126         kU = 2 + deltaK*(Nr-3)
                0127         kL = Nr - deltaK*(Nr-1)
                0128         DO k = kU,kL,kDir
6aab75f3e8 Jean*0129 C-    potential density (reference level = surface level)
0f5ba23f97 Jean*0130          CALL FIND_RHO_2D(
0320e25227 Mart*0131      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, kSrf,
0f5ba23f97 Jean*0132      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0133      O        rhoLoc,
                0134      I        k, bi, bj, myThid )
6aab75f3e8 Jean*0135
0320e25227 Mart*0136          DO j=1-OLy,sNy+OLy
                0137           DO i=1-OLx,sNx+OLx
                0138 c          kIn = k.LT.klowC(i,j,bi,bj).AND.k.GE.kSurfC(i,j,bi,bj)
                0139            kIn = k.LE.klowC(i,j,bi,bj).AND.k.GE.kSurfC(i,j,bi,bj)
                0140            IF ( kIn .AND. rhoLoc(i,j).GE.rhoMxL(i,j) ) THEN
6aab75f3e8 Jean*0141              IF ( rhoLoc(i,j).GT.rhoKm1(i,j) ) THEN
                0142               tmpFac = ( rhoMxL(i,j) - rhoKm1(i,j) )
                0143      &               / ( rhoLoc(i,j) - rhoKm1(i,j) )
                0144              ELSE
                0145               tmpFac = 0.
                0146              ENDIF
0320e25227 Mart*0147 C-    Note: sign is more related to kDir than it is to gravity orientation
                0148 c            hMixLayer(i,j,bi,bj) = ( rF(kTop)-rC(k-kDir) )*kDir
                0149              hMixLayer(i,j,bi,bj) = -gravitySign*( rF(kTop)-rC(k-kDir) )
                0150      &                            + tmpFac*drC(k+deltaK)
6aab75f3e8 Jean*0151              rhoMxL(i,j) = rhoBigNb
                0152            ELSE
                0153              rhoKm1(i,j) = rhoLoc(i,j)
                0154            ENDIF
                0155           ENDDO
                0156          ENDDO
                0157         ENDDO
                0158
                0159        ELSEIF ( method.EQ.2 ) THEN
                0160
                0161 C--   Second method :
                0162 C     where the local stratification exceed the mean stratification above
                0163 C     (from surface down to here) by factor hMixCriteria
                0164
7c01d59526 Davi*0165 c       hMixCriteria  = 1.5 _d 0
                0166 c       dRhoSmall = 1. _d -2
0320e25227 Mart*0167         DO j=1-OLy,sNy+OLy
                0168          DO i=1-OLx,sNx+OLx
                0169           IF ( klowC(i,j,bi,bj) .GT. 0 ) THEN
                0170            hMixLayer(i,j,bi,bj) = drF(kSrf)
                0171            rhoMxL(i,j) = 1.
                0172           ELSE
                0173            hMixLayer(i,j,bi,bj) = rF(kTop)
                0174            rhoMxL(i,j) = -1.
                0175           ENDIF
6aab75f3e8 Jean*0176          ENDDO
                0177         ENDDO
0320e25227 Mart*0178 C-    Z-coord: (kU,kL,kDir) = (2,Nr-1,1) ; P-coord: (kU,kL,kDir)= (Nr-1,2,-1)
                0179         kU = 2 + deltaK*(Nr-3)
                0180         kL = Nr-1 - deltaK*(Nr-3)
                0181         DO k = kU,kL,kDir
6aab75f3e8 Jean*0182 C-    potential density (reference level = surface level)
0f5ba23f97 Jean*0183          CALL FIND_RHO_2D(
0320e25227 Mart*0184      I        1-OLx, sNx+OLx, 1-OLy, sNy+OLy, kSrf,
0f5ba23f97 Jean*0185      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0186      O        rhoLoc,
                0187      I        k, bi, bj, myThid )
6aab75f3e8 Jean*0188
0320e25227 Mart*0189          kup = k+deltaK
                0190          klw = kup+kDir
                0191          DO j=1-OLy,sNy+OLy
                0192           DO i=1-OLx,sNx+OLx
                0193            kIn = k.LT.klowC(i,j,bi,bj).AND.k.GT.kSurfC(i,j,bi,bj)
                0194            IF ( kIn .AND. rhoMxL(i,j).GE.0. ) THEN
39de785bb5 Jean*0195              sigmAv = ( rhoLoc(i,j)-rhoSurf(i,j)+dRhoSmall )
0320e25227 Mart*0196      &              / ( rC(k)-rC(kSrf) ) * gravitySign
                0197              IF ( gravitySign * sigmaR(i,j,klw)
                0198      &            .GT.sigmAv*hMixCriteria ) THEN
6aab75f3e8 Jean*0199                tmpFac = 0. _d 0
                0200                IF ( sigmAv.GT.0. _d 0 ) THEN
0320e25227 Mart*0201                  tmpFac = hMixCriteria*sigmaR(i,j,kup)/sigmaR(i,j,klw)
39de785bb5 Jean*0202                  IF ( tmpFac .GT. 1. _d 0 ) THEN
                0203                    tmpFac = 1. _d 0
                0204      &             + ( tmpFac - 1. _d 0 )/( hMixCriteria - 1. _d 0 )
                0205                  ENDIF
                0206                  tmpFac = MAX( 0. _d 0, MIN( tmpFac, 2. _d 0 ) )
6aab75f3e8 Jean*0207                ENDIF
0320e25227 Mart*0208 C-    Note: sign is more related to kDir than it is to gravity orientation
                0209                hMixLayer(i,j,bi,bj) = -gravitySign *
                0210      &                               ( rF(kTop)-rF(klw) )
                0211      &                               - drF(k)*tmpFac*0.5 _d 0
6aab75f3e8 Jean*0212                rhoMxL(i,j) = -1.
                0213              ENDIF
                0214            ENDIF
                0215           ENDDO
                0216          ENDDO
                0217         ENDDO
                0218
                0219        ELSE
                0220         STOP 'S/R CALC_OCE_MXLAYER: invalid method'
                0221        ENDIF
                0222
d1b0368d70 Davi*0223        IF ( hMixSmooth .GT. 0. _d 0 ) THEN
                0224         tmpFac = (1. _d 0 - hMixSmooth ) / 4. _d 0
0320e25227 Mart*0225         DO j=1-OLy+1,sNy+OLy-1
                0226          DO i=1-OLx+1,sNx+OLx-1
d1b0368d70 Davi*0227             rhoLoc(i,j)=(hMixSmooth *   hMixLayer(i,j,bi,bj)   +
                0228      &                       tmpFac * ( hMixLayer(i-1,j,bi,bj) +
                0229      &                                  hMixLayer(i+1,j,bi,bj) +
                0230      &                                  hMixLayer(i,j-1,bi,bj) +
                0231      &                                  hMixLayer(i,j+1,bi,bj) )
                0232      &                  )
                0233      &                 /(hMixSmooth +
0320e25227 Mart*0234      &                       tmpFac * ( maskC(i-1,j,kSrf,bi,bj) +
                0235      &                                  maskC(i+1,j,kSrf,bi,bj) +
                0236      &                                  maskC(i,j-1,kSrf,bi,bj) +
                0237      &                                  maskC(i,j+1,kSrf,bi,bj) )
                0238      &                  ) * maskC(i,j,kSrf,bi,bj)
d1b0368d70 Davi*0239          ENDDO
                0240         ENDDO
0320e25227 Mart*0241         DO j=1-OLy+1,sNy+OLy-1
                0242          DO i=1-OLx+1,sNx+OLx-1
d1b0368d70 Davi*0243             hMixLayer(i,j,bi,bj) = rhoLoc(i,j)
                0244          ENDDO
                0245         ENDDO
                0246        ENDIF
                0247
6aab75f3e8 Jean*0248 #ifdef ALLOW_DIAGNOSTICS
                0249        IF ( useDiagnostics ) THEN
                0250         CALL DIAGNOSTICS_FILL( hMixLayer, 'MXLDEPTH',
                0251      &                         0, 1, 1, bi, bj, myThid )
                0252        ENDIF
                0253 #endif /* ALLOW_DIAGNOSTICS */
                0254
                0255 C--   end if calcMixLayerDepth
                0256       ENDIF
                0257
                0258       RETURN
                0259       END