MODULE zdfmxl_crs
   !!======================================================================
   !!                       ***  MODULE  zdfmxl  ***
   !! Ocean physics: mixed layer depth 
   !!======================================================================
   !! History :  1.0  ! 2003-08  (G. Madec)  original code
   !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop
   !!            3.7  ! 2012-03  (G. Madec)  make public the density criteria for trdmxl 
   !!             -   ! 2014-02  (F. Roquet)  mixed layer depth calculated using N2 instead of rhop 
   !!----------------------------------------------------------------------
   !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
   !!----------------------------------------------------------------------
   !USE oce             ! ocean dynamics and tracers variables
   !USE dom_oce         ! ocean space and time domain variables
   !USE oce_trc
   USE zdf_oce         ! ocean vertical physics
   USE in_out_manager  ! I/O manager
   USE prtctl          ! Print control
   USE phycst          ! physical constants
   USE iom             ! I/O library
   USE lib_mpp         ! MPP library
   USE wrk_nemo        ! work arrays
   USE timing          ! Timing
   USE trc_oce, ONLY : lk_offline ! offline flag
   USE crs

   IMPLICIT NONE
   PRIVATE

   PUBLIC   zdf_mxl_crs       ! called by step.F90

   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth

   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
   !! $Id: zdfmxl.F90 4990 2014-12-15 16:42:49Z timgraham $ 
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE zdf_mxl_crs( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zdfmxl  ***
      !!                   
      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
      !!              with density criteria.
      !!
      !! ** Method  :   The mixed layer depth is the shallowest W depth with 
      !!      the density of the corresponding T point (just bellow) bellow a
      !!      given value defined locally as rho(10m) + rho_c
      !!               The turbocline depth is the depth at which the vertical
      !!      eddy diffusivity coefficient (resulting from the vertical physics
      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
      !!      value defined locally (avt_c here taken equal to 5 cm/s2 by default)
      !!
      !! ** Action  :   nmln, hmld, hmlp, hmlpt
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
      !
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      INTEGER  ::   iiki,iikn   ! local integer
      REAL(wp) ::   zN2_c        ! local scalar
      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
      REAL(wp), POINTER, DIMENSION(:,:) ::  z2d   ! 2D workspace
      !!----------------------------------------------------------------------
      !
      IF( nn_timing == 1 )  CALL timing_start('zdf_mxl_crs')
      !
      CALL wrk_alloc( jpi_crs,jpj_crs, imld )
      CALL wrk_alloc( jpi_crs,jpj_crs, z2d )

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'zdf_mxl_crs : mixed layer depth'
         IF(lwp) WRITE(numout,*) '~~~~~~~ '
      ENDIF

      ! w-level of the turbocline
      imld(:,:)=0
      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10 
         DO jj = 1, jpj_crs
            DO ji = 1, jpi_crs
               IF( avt_crs (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( jk, 1 )      ! Turbocline 
            END DO
         END DO
      END DO

      ! depth of the mixing and mixed layers
      hmld_crs(:,:) = 0._wp
      hmlpt_crs(:,:) = 0._wp
      DO jj = 1, jpj_crs
         DO ji = 1, jpi_crs
            iiki = imld(ji,jj)
            iikn = nmln_crs(ji,jj)
            IF( iiki .NE. 0 ) hmld_crs (ji,jj) = ( gdepw_crs(ji,jj,iiki  ) - gdepw_crs(ji,jj,1 )            )   * tmask_crs(ji,jj,1)  ! Turbocline depth 
            IF( iiki .NE. 0 ) hmlpt_crs(ji,jj) = ( gdept_crs(ji,jj,iikn-1) - gdepw_crs(ji,jj,1 )            )   * tmask_crs(ji,jj,1)  ! depth of the last T-point inside the mixed layer
         END DO
      END DO
      !
      z2d=REAL(nmln_crs,wp)
      CALL iom_put("nmln_crs",z2d)
      CALL iom_put("hmlpt_crs",hmlpt_crs)
      !
      CALL wrk_dealloc( jpi_crs,jpj_crs, imld )
      CALL wrk_dealloc( jpi_crs,jpj_crs, z2d )
      !
      IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl_crs')
      !
   END SUBROUTINE zdf_mxl_crs

   !!======================================================================
END MODULE zdfmxl_crs