Index: trunk/NEMO/OPA_SRC/DIA/diahth.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diahth.F90	(revision 1551)
+++ trunk/NEMO/OPA_SRC/DIA/diahth.F90	(revision 1577)
@@ -9,5 +9,5 @@
    !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
    !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
-   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning
+   !!   NEMO     3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
    !!----------------------------------------------------------------------
 
@@ -16,9 +16,10 @@
    !!   'key_diahth' :                              thermocline depth diag.
    !!----------------------------------------------------------------------
-   !!   dia_hth      : Compute diagnostics associated with the thermocline
+   !!   dia_hth      : Compute varius diagnostics associated with the mixed layer
    !!----------------------------------------------------------------------
    !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE zdf_oce         ! ocean vertical physics
    USE phycst          ! physical constants
    USE in_out_manager  ! I/O manager
@@ -32,10 +33,11 @@
 
    !! * Shared module variables
-   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      hth  ,      &  !: depth of the max vertical temperature gradient (m)
-      hd20 ,      &  !: depth of 20 C isotherm (m)
-      hd28 ,      &  !: depth of 28 C isotherm (m)
-      htc3           !: heat content of first 300 m
+   LOGICAL , PUBLIC, PARAMETER          ::   lk_diahth = .TRUE.   !: thermocline-20d depths flag
+   ! note: following variables should move to local variables once iom_put is always used 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hmld                 !: mixing layer depth (turbocline)                [m]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hth                  !: depth of the max vertical temperature gradient [m]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd20                 !: depth of 20 C isotherm                         [m]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hd28                 !: depth of 28 C isotherm                         [m]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   htc3                 !: heat content of first 300 m                    [W]
 
    !! * Substitutions
@@ -53,9 +55,17 @@
       !!                  ***  ROUTINE dia_hth  ***
       !!
-      !! ** Purpose :
-      !!      Computes the depth of strongest vertical temperature gradient
-      !!      Computes the depth of the 20 degree isotherm
-      !!      Computes the depth of the 28 degree isotherm
-      !!      Computes the heat content of first 300 m
+      !! ** Purpose : Computes
+      !!      the mixing layer depth (turbocline): avt = 5.e-4
+      !!      the depth of strongest vertical temperature gradient
+      !!      the mixed layer depth with density     criteria: rho = rho(10m or surf) + 0.03(or 0.01)
+      !!      the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2       
+      !!      the top of the thermochine: tn = tn(10m) - ztem2 
+      !!      the pycnocline depth with density criteria equivalent to a temperature variation 
+      !!                rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC) 
+      !!      the barrier layer thickness
+      !!      the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
+      !!      the depth of the 20 degree isotherm (linear interpolation)
+      !!      the depth of the 28 degree isotherm (linear interpolation)
+      !!      the heat content of first 300 m
       !!
       !! ** Method : 
@@ -66,8 +76,23 @@
       INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
       INTEGER                          ::   iid, iif, ilevel      ! temporary integers
-      INTEGER, DIMENSION(jpi,jpj)      ::   ikc                   ! levels
-      REAL(wp)                         ::   zd, zthick_0, zcoef   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zthick
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzt
+      INTEGER, DIMENSION(jpi,jpj)      ::   ik20, ik28            ! levels
+      REAL(wp)                         ::   zavt5 = 5.e-4_wp      ! Kz criterion for the turbocline depth
+      REAL(wp)                         ::   zrho3 = 0.03_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                         ::   zrho1 = 0.01_wp       ! density     criterion for mixed layer depth
+      REAL(wp)                         ::   ztem2 = 0.2_wp        ! temperature criterion for mixed layer depth
+      REAL(wp)                         ::   zthick_0, zcoef       ! temporary scalars
+      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
+      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zabs2                 ! MLD: abs( tn - tn(10m) ) = ztem2 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztm2                  ! Top of thermocline: tn = tn(10m) - ztem2     
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho10_3              ! MLD: rho = rho10m + zrho3      
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zpycn                 ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
+      REAL(wp), DIMENSION(jpi,jpj)     ::   ztinv                 ! max of temperature inversion
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdepinv               ! depth of temperature inversion
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_3               ! MLD rho = rho(surf) = 0.03
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zrho0_1               ! MLD rho = rho(surf) = 0.01
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zmaxdzT               ! max of dT/dz
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zthick                ! vertical integration thickness 
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zdelr                 ! delta rho equivalent to deltaT = 0.2
       !!----------------------------------------------------------------------
 
@@ -79,63 +104,177 @@
       ENDIF
 
-      ! -------------------------- !
-      !  Depth of the thermocline  !
-      ! -------------------------- !
-      ! The depth of the thermocline is defined as the depth of the 
-      ! strongest vertical temperature gradient
-      zdzt(:,:,1) = 0.e0
-      DO jk = 2, jpk                      ! vertical gradient of temperature
-         zdzt(:,:,jk) = ( tn(:,:,jk-1) - tn(:,:,jk) ) / fse3w(:,:,jk) * tmask(:,:,jk)
-      END DO
+      ! initialization
+      ztinv  (:,:) = 0.e0_wp  
+      zdepinv(:,:) = 0.e0_wp  
+      zmaxdzT(:,:) = 0.e0_wp  
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ilevel = MAXLOC( zdzt(ji,jj,:), dim= 1 )      ! level of maximum vertical temperature gradient
-            hth(ji,jj) = fsdepw(ji,jj,ilevel)             ! depth of the thermocline
-         END DO         
-      END DO
-
-      CALL iom_put( "thermod", hth )      ! depth of the thermocline
-
-      ! ----------------------- !
-      !  Depth of 20C isotherm  !
-      ! ----------------------- !
+            zztmp = bathy(ji,jj)
+            hth     (ji,jj) = zztmp
+            hmld    (ji,jj) = zztmp
+            zabs2   (ji,jj) = zztmp
+            ztm2    (ji,jj) = zztmp
+            zrho10_3(ji,jj) = zztmp
+            zpycn   (ji,jj) = zztmp
+        END DO
+      END DO
+      IF( nla10 > 1 ) THEN 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp = bathy(ji,jj)
+               zrho0_3(ji,jj) = zztmp
+               zrho0_1(ji,jj) = zztmp
+            END DO
+         END DO
+      ENDIF
       
-      ! search last level above 20C (beware temperature is not always decreasing with depth)
-      ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 20., dim = 3 )   !   ikc = 0 if mask always equal to false
-      ! Depth of 20C isotherm, linear interpolation
+      ! Preliminary computation
+      ! computation of zdelr = (dr/dT)(T,S,10m)*(-0.2 degC)
+      DO jj=1, jpj
+         DO ji=1, jpi
+            IF( tmask(ji,jj,nla10) == 1. ) THEN
+               zu  =  1779.50 + 11.250*tn(ji,jj,nla10) - 3.80*sn(ji,jj,nla10) - 0.0745*tn(ji,jj,nla10)*tn(ji,jj,nla10)   &
+                  &                                                           - 0.0100*tn(ji,jj,nla10)*sn(ji,jj,nla10)
+               zv  =  5891.00 + 38.000*tn(ji,jj,nla10) + 3.00*sn(ji,jj,nla10) - 0.3750*tn(ji,jj,nla10)*tn(ji,jj,nla10)
+               zut =    11.25 -  0.149*tn(ji,jj,nla10) - 0.01*sn(ji,jj,nla10)
+               zvt =    38.00 -  0.750*tn(ji,jj,nla10)
+               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
+               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+            ELSE
+               zdelr(ji,jj) = 0.e0
+            ENDIF
+         END DO
+      END DO
+
+      ! ------------------------------------------------------------- !
+      ! thermocline depth: strongest vertical gradient of temperature !
+      ! turbocline depth (mixing layer depth): avt = zavt5            !
+      ! MLD: rho = rho(1) + zrho3                                     !
+      ! MLD: rho = rho(1) + zrho1                                     !
+      ! ------------------------------------------------------------- !
+      DO jk = jpkm1, 2, -1   ! loop from bottom to 2
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zzdep = fsdepw(ji,jj,jk)
+               zztmp = ( tn(ji,jj,jk-1) - tn(ji,jj,jk) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
+               zzdep = zzdep * tmask(ji,jj,1)
+
+               IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
+                  zmaxdzT(ji,jj) = zztmp   ;   hth    (ji,jj) = zzdep                ! max and depth of dT/dz
+               ENDIF
+               
+               IF( avt (ji,jj,jk) < zavt5 )    hmld   (ji,jj) = zzdep                ! avt < zavt5
+
+               IF( nla10 > 1 ) THEN 
+                  zztmp = rhop(ji,jj,jk) - rhop(ji,jj,1)                             ! delta rho(1)
+                  IF( zztmp > zrho3 )          zrho0_3(ji,jj) = zzdep                ! > 0.03
+                  IF( zztmp > zrho1 )          zrho0_1(ji,jj) = zzdep                ! > 0.01
+               ENDIF
+
+            END DO
+         END DO
+      END DO
+      
+      CALL iom_put( "mlddzt", hth )            ! depth of the thermocline
+      CALL iom_put( "mldkz5", hmld )           ! turbocline depth
+      IF( nla10 > 1 ) THEN 
+         CALL iom_put( "mldr0_3", zrho0_3 )   ! MLD delta rho(surf) = 0.03
+         CALL iom_put( "mldr0_1", zrho0_1 )   ! MLD delta rho(surf) = 0.01
+      ENDIF
+
+      ! ------------------------------------------------------------- !
+      ! MLD: abs( tn - tn(10m) ) = ztem2                              !
+      ! Top of thermocline: tn = tn(10m) - ztem2                      !
+      ! MLD: rho = rho10m + zrho3                                     !
+      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)       !
+      ! temperature inversion: max( 0, max of tn - tn(10m) )          !
+      ! depth of temperature inversion                                !
+      ! ------------------------------------------------------------- !
+      DO jk = jpkm1, nlb10, -1   ! loop from bottom to nlb10
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+
+               zzdep = fsdepw(ji,jj,jk) * tmask(ji,jj,1)
+
+               zztmp = tn(ji,jj,nla10) - tn(ji,jj,jk)                  ! - delta T(10m)
+               IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
+               IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
+               zztmp = -zztmp                                          ! delta T(10m)
+               IF( zztmp >  ztinv(ji,jj) ) THEN                        ! temperature inversion
+                  ztinv(ji,jj) = zztmp   ;   zdepinv (ji,jj) = zzdep   ! max value and depth
+               ENDIF
+
+               zztmp = rhop(ji,jj,jk) - rhop(ji,jj,nla10)              ! delta rho(10m)
+               IF( zztmp > zrho3        )    zrho10_3(ji,jj) = zzdep   ! > 0.03
+               IF( zztmp > zdelr(ji,jj) )    zpycn   (ji,jj) = zzdep   ! > equi. delta T(10m) - 0.2
+
+            END DO
+         END DO
+      END DO
+
+      CALL iom_put( "mld|dt|" , zabs2        )   ! MLD abs(delta t) - 0.2
+      CALL iom_put( "topthdep", ztm2         )   ! T(10) - 0.2
+      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03
+      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2
+      CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness
+      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref) 
+      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref) 
+
+
+      ! ----------------------------------- !
+      ! search deepest level above 20C/28C  !
+      ! ----------------------------------- !
+      ik20(:,:) = 1
+      ik28(:,:) = 1
+      DO jk = 1, jpkm1   ! beware temperature is not always decreasing with depth => loop from top to bottom
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp = tn(ji,jj,jk)
+               IF( zztmp >= 20. )   ik20(ji,jj) = jk
+               IF( zztmp >= 28. )   ik28(ji,jj) = jk
+            END DO
+         END DO
+      END DO
+
+      ! --------------------------- !
+      !  Depth of 20C/28C isotherm  !
+      ! --------------------------- !
       DO jj = 1, jpj
          DO ji = 1, jpi
-            iid = MAX( 1, ikc(ji,jj) )    ! make sure that iid /= 0 (these points are masked later)
-            zd = fsdept(ji,jj,iid) + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                           )   &
-               &                   * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                           )   &
-               &                   / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) )
-            ! bound by the ocean depth, minimum value, first T-point depth
+
             iif = mbathy(ji,jj)
-            hd20(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) )
-         END DO
-      END DO
-      WHERE(ikc == 0 )   hd20 = 0.e0      ! set to 0 points where tn is always < 20 (not very good for temporal mean...)
-      CALL iom_put( "20d", hd20 )         ! depth of the 20 isotherm
-
-      ! ----------------------- !
-      !  Depth of 28C isotherm  ! 
-      ! ----------------------- !
-      
-      ! search last level above 28C (beware temperature is not always decreasing with depth)
-      ikc(:,:) = MAXLOC( fsdept(:,:,1:jpkm1), mask = tn(:,:,1:jpkm1) >= 28., dim = 3 )   !   ikc = 0 if mask always equal to false
-      ! Depth of 28C isotherm, linear interpolation
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            iid = MAX( 1, ikc(ji,jj) )    ! make sure that iid /= 0 (these points are masked later)
-            zd = fsdept(ji,jj,iid) + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                           )   &
-               &                   * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                           )   &
-               &                   / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + ( 1. - tmask(ji,jj,1) ) )
-            ! bound by the ocean depth, minimum value, first T-point depth
-            iif = mbathy(ji,jj)
-            hd28(ji,jj) = MIN( zd*tmask(ji,jj,1), fsdepw(ji,jj,iif) )
-         END DO
-      END DO
-      WHERE(ikc == 0 )   hd28 = 0.e0      ! set to 0 points where tn is always < 28 (not very good for temporal mean...)
-      CALL iom_put( "28d", hd28 )         ! depth of the 28 isotherm
+            zzdep = fsdepw(ji,jj,iif)
+
+            iid = ik20(ji,jj)
+            IF( iid /= 1 ) THEN 
+               ! linear interpolation
+               zztmp =      fsdept(ji,jj,iid  )   &
+                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
+                  &  * ( 20.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
+                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
+               ! bound by the ocean depth, minimum value, first T-point depth
+               hd20(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep)
+            ELSE 
+               hd20(ji,jj)=0.
+            ENDIF
+
+            iid = ik28(ji,jj)
+            IF( iid /= 1 ) THEN 
+               ! linear interpolation
+               zztmp =      fsdept(ji,jj,iid  )   &
+                  &  + (    fsdept(ji,jj,iid+1) - fsdept(ji,jj,iid)                       )   &
+                  &  * ( 28.*tmask(ji,jj,iid+1) -     tn(ji,jj,iid)                       )   &
+                  &  / (        tn(ji,jj,iid+1) -     tn(ji,jj,iid) + (1.-tmask(ji,jj,1)) )
+               ! bound by the ocean depth, minimum value, first T-point depth
+               hd28(ji,jj) = MIN( zztmp*tmask(ji,jj,1), zzdep )
+            ELSE 
+               hd28(ji,jj) = 0.
+            ENDIF
+
+         END DO
+      END DO
+      CALL iom_put( "20d", hd20 )   ! depth of the 20 isotherm
+      CALL iom_put( "28d", hd28 )   ! depth of the 28 isotherm
 
       ! ----------------------------- !
@@ -145,5 +284,5 @@
       ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to do this search...)
       ilevel = 0
-      zthick_0 = 0.e0
+      zthick_0 = 0.e0_wp
       DO jk = 1, jpkm1                      
          zthick_0 = zthick_0 + e3t_0(jk)
@@ -151,5 +290,5 @@
       END DO
       ! surface boundary condition
-      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0        ;   htc3(:,:) = 0.e0                                     
+      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0.e0_wp     ;   htc3(:,:) = 0.e0_wp                                   
       ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tn(:,:,jk) * sshn(:,:) * tmask(:,:,jk)   
       ENDIF
@@ -169,5 +308,5 @@
       zcoef = rau0 * rcp
       htc3(:,:) = zcoef * htc3(:,:)
-      CALL iom_put( "hc300", htc3 )      ! first 300m heaat content
+      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
 
 
