C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/calc_oce_mxlayer.F,v 1.1 2007/06/03 22:50:47 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: CALC_OCE_MXLAYER
C     !INTERFACE:
      SUBROUTINE CALC_OCE_MXLAYER(
     I                       rhoSurf, sigmaR,
     I                       bi, bj, myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R CALC_OCE_MXLAYER
C     | o Diagnose the Oceanic surface Mixed-Layer
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#ifdef ALLOW_GMREDI
# include "GMREDI.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine Arguments ==
C     rhoSurf   :: Surface density anomaly
C     sigmaR    :: Vertical gradient of potential density
C     bi,bj     :: tile indices
C     myTime    :: Current time in simulation
C     myIter    :: Current iteration number in simulation
C     myThid    :: my Thread Id number
      _RL     rhoSurf(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     sigmaR (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      INTEGER bi, bj
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C     === Functions ====
#ifdef ALLOW_DIAGNOSTICS
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
#endif /* ALLOW_DIAGNOSTICS */

C     !LOCAL VARIABLES:
C     == Local variables ==
C     i,j :: Loop counters
      INTEGER i,j,k
      LOGICAL calcMixLayerDepth
      INTEGER method
      _RL     dRhoSmall, rhoBigNb
      _RL     rhoMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     rhoKm1(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     tmpFac, sigmAv
CEOP

      calcMixLayerDepth = .FALSE.
#ifdef ALLOW_GMREDI
      IF ( useGMRedi ) calcMixLayerDepth = GM_taper_scheme.EQ.'fm07'
#endif
#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics.AND. .NOT.calcMixLayerDepth ) THEN
        calcMixLayerDepth = DIAGNOSTICS_IS_ON('MXLDEPTH',myThid)
      ENDIF
#endif
      IF ( calcMixLayerDepth ) THEN

C--   Select which "method" to use:
       method = 0
       IF ( hMixCriteria.LT.0. ) method = 1
       IF ( hMixCriteria.GT.1. ) method = 2

       IF ( method.EQ.1 ) THEN

C--   First method :
C     where the potential density (ref.lev=surface) is larger than
C       surface density plus Delta_Rho = hMixCriteria * Alpha(surf)
C     = density of water which is -hMixCriteria colder than surface water

c       hMixCriteria  = -0.8 _d 0
        dRhoSmall = 1. _d -6
        rhoBigNb  = rhoConst*1. _d 10
        CALL FIND_ALPHA(
     I            bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, 1, 1,
     O            rhoMxL, myThid )

        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
           rhoKm1(i,j) = rhoSurf(i,j)
           rhoMxL(i,j) = rhoSurf(i,j)
     &                 + MAX( rhoMxL(i,j)*hMixCriteria, dRhoSmall )
           hMixLayer(i,j,bi,bj) = rF(1)-R_low(I,J,bi,bj)
         ENDDO
        ENDDO
        DO k = 2,Nr
C-    potential density (reference level = surface level)
         CALL FIND_RHO(
     I             bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, 1,
     I             theta, salt,
     O             rhoLoc, myThid )

         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           IF ( k.LE.klowC(i,j,bi,bj) .AND.
     &          rhoLoc(i,j).GE.rhoMxL(i,j) ) THEN
             IF ( rhoLoc(i,j).GT.rhoKm1(i,j) ) THEN
              tmpFac = ( rhoMxL(i,j) - rhoKm1(i,j) )
     &               / ( rhoLoc(i,j) - rhoKm1(i,j) )
             ELSE
              tmpFac = 0.
             ENDIF
             hMixLayer(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
             rhoMxL(i,j) = rhoBigNb
           ELSE
             rhoKm1(i,j) = rhoLoc(i,j)
           ENDIF
          ENDDO
         ENDDO
        ENDDO

       ELSEIF ( method.EQ.2 ) THEN

C--   Second method :
C     where the local stratification exceed the mean stratification above
C     (from surface down to here) by factor hMixCriteria

c       hMixCriteria  = 2. _d 0
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
           IF ( klowC(i,j,bi,bj) .GT. 0 ) THEN
            hMixLayer(i,j,bi,bj) = drF(1)
            rhoMxL(i,j) = 1.
           ELSE
            hMixLayer(i,j,bi,bj) = rF(1)
            rhoMxL(i,j) = -1.
           ENDIF
         ENDDO
        ENDDO
        DO k = 3,Nr
C-    potential density (reference level = surface level)
         CALL FIND_RHO(
     I             bi, bj, 1-OLx, sNx+OLx, 1-OLy, sNy+OLy, K, 1,
     I             theta, salt,
     O             rhoLoc, myThid )

         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           IF ( k.LE.klowC(i,j,bi,bj)
     &          .AND. rhoMxL(i,j).GE.0. ) THEN
             sigmAv = ( rhoLoc(i,j)-rhoSurf(i,j) ) / ( rC(1)-rC(k) )
             IF ( -sigmaR(i,j,k).GE.sigmAv*hMixCriteria ) THEN
               tmpFac = 0. _d 0
               IF ( sigmAv.GT.0. _d 0 ) THEN
                 tmpFac = sigmaR(i,j,k-1)/sigmaR(i,j,k)
                 tmpFac = ( tmpFac*hMixCriteria - 1. _d 0 )
     &                         / ( hMixCriteria - 1. _d 0 )
                 tmpFac = MAX( 0. _d 0, MIN( tmpFac, 1. _d 0 ) )
               ENDIF
               hMixLayer(i,j,bi,bj) = rF(1)-rC(k-1)+tmpFac*drC(k)
               rhoMxL(i,j) = -1.
             ENDIF
           ENDIF
          ENDDO
         ENDDO
        ENDDO

       ELSE
        STOP 'S/R CALC_OCE_MXLAYER: invalid method'
       ENDIF

#ifdef ALLOW_DIAGNOSTICS
       IF ( useDiagnostics ) THEN
        CALL DIAGNOSTICS_FILL( hMixLayer, 'MXLDEPTH',
     &                         0, 1, 1, bi, bj, myThid )
       ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C--   end if calcMixLayerDepth
      ENDIF

      RETURN
      END