Index: /NEMO/releases/release-4.0/src/OCE/DIA/diahth.F90
===================================================================
--- /NEMO/releases/release-4.0/src/OCE/DIA/diahth.F90	(revision 11090)
+++ /NEMO/releases/release-4.0/src/OCE/DIA/diahth.F90	(revision 11091)
@@ -5,11 +5,10 @@
    !!======================================================================
    !! History :  OPA  !  1994-09  (J.-P. Boulanger)  Original code
-   !!                 !  1996-11  (E. Guilyardi)  OPA8 
+   !!                 !  1996-11  (E. Guilyardi)  OPA8
    !!                 !  1997-08  (G. Madec)  optimization
-   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content 
+   !!                 !  1999-07  (E. Guilyardi)  hd28 + heat content
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!            3.2  !  2009-07  (S. Masson) hc300 bugfix + cleaning + add new diag
    !!----------------------------------------------------------------------
-#if   defined key_diahth
    !!----------------------------------------------------------------------
    !!   'key_diahth' :                              thermocline depth diag.
@@ -24,5 +23,4 @@
    USE lib_mpp         ! MPP library
    USE iom             ! I/O library
-   USE timing          ! preformance summary
 
    IMPLICIT NONE
@@ -30,7 +28,8 @@
 
    PUBLIC   dia_hth       ! routine called by step.F90
-
-   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .TRUE.    !: thermocline-20d depths flag
-   
+   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90
+
+   LOGICAL, PUBLIC ::   ln_diahth   !: Compute further diagnostics of ML and thermocline depth
+
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -49,8 +48,8 @@
      !!      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 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) )
@@ -59,5 +58,5 @@
      !!      the heat content of first 300 m
      !!
-     !! ** Method : 
+     !! ** Method :
      !!-------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
@@ -65,5 +64,5 @@
      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
      INTEGER                          ::   iid, ilevel           ! temporary integers
-     INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ik20, ik28  ! levels
+     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
@@ -73,79 +72,20 @@
      REAL(wp)                         ::   zztmp, zzdep          ! temporary scalars inside do loop
      REAL(wp)                         ::   zu, zv, zw, zut, zvt  ! temporary workspace
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zabs2      ! MLD: abs( tn - tn(10m) ) = ztem2 
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztm2       ! Top of thermocline: tn = tn(10m) - ztem2     
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho10_3   ! MLD: rho = rho10m + zrho3      
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zpycn      ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ztinv      ! max of temperature inversion
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdepinv    ! depth of temperature inversion
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_3    ! MLD rho = rho(surf) = 0.03
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zrho0_1    ! MLD rho = rho(surf) = 0.01
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zmaxdzT    ! max of dT/dz
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zthick     ! vertical integration thickness 
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zdelr      ! delta rho equivalent to deltaT = 0.2
-   ! note: following variables should move to local variables once iom_put is always used 
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhth    !: depth of the max vertical temperature gradient [m]
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhd20   !: depth of 20 C isotherm                         [m]
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhd28   !: depth of 28 C isotherm                         [m]
-     REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zhtc3   !: heat content of first 300 m                    [W]
-
-     !!----------------------------------------------------------------------
-     IF( ln_timing )   CALL timing_start('dia_hth')
-
-     IF( kt == nit000 ) THEN
-        !                                      ! allocate dia_hth array
-        IF(.NOT. ALLOCATED(ik20) ) THEN
-           ALLOCATE(ik20(jpi,jpj), ik28(jpi,jpj), &
-                &      zhth(jpi,jpj),   &
-                &      zhd20(jpi,jpj),   &
-                &      zhd28(jpi,jpj),   &
-                &      zhtc3(jpi,jpj),   &
-                &      zabs2(jpi,jpj),   &
-                &      ztm2(jpi,jpj),    &
-                &      zrho10_3(jpi,jpj),&
-                &      zpycn(jpi,jpj),   &
-                &      ztinv(jpi,jpj),   &
-                &      zdepinv(jpi,jpj), &
-                &      zrho0_3(jpi,jpj), &
-                &      zrho0_1(jpi,jpj), &
-                &      zmaxdzT(jpi,jpj), &
-                &      zthick(jpi,jpj),  &
-                &      zdelr(jpi,jpj), STAT=ji)
-           IF( lk_mpp  )   CALL mpp_sum(ji)
-           IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard ocean arrays' )
-        END IF
-
-        IF(lwp) WRITE(numout,*)
-        IF(lwp) WRITE(numout,*) 'dia_hth : diagnostics of the thermocline depth'
-        IF(lwp) WRITE(numout,*) '~~~~~~~ '
-        IF(lwp) WRITE(numout,*)
-     ENDIF
-
-     IF(iom_use("mlddzt").OR.iom_use("mldr0_3").OR.iom_use("mldr0_1").OR.iom_use("mld_dt02") &
-          & .OR.iom_use("topthdep").OR.iom_use("mldr10_3").OR.iom_use("pycndep").OR.iom_use("tinv").OR.iom_use("depti")) THEN
-        ! initialization
-        ztinv  (:,:) = 0._wp  
-        zdepinv(:,:) = 0._wp  
-        zmaxdzT(:,:) = 0._wp  
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-              zhth     (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 = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
-                 zrho0_3(ji,jj) = zztmp
-                 zrho0_1(ji,jj) = zztmp
-              END DO
-           END DO
-        ENDIF
-     ENDIF
+     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
+   ! note: following variables should move to local variables once iom_put is always used
+     REAL(wp), DIMENSION(jpi,jpj) ::   zhth    !: depth of the max vertical temperature gradient [m]
+     REAL(wp), DIMENSION(jpi,jpj) ::   zhd20   !: depth of 20 C isotherm                         [m]
+     REAL(wp), DIMENSION(jpi,jpj) ::   zhd28   !: depth of 28 C isotherm                         [m]
+     REAL(wp), DIMENSION(jpi,jpj) ::   zhtc3   !: heat content of first 300 m                    [W]
 
      IF (iom_use("mlddzt").OR.iom_use("mldr0_3").OR.iom_use("mldr0_1")) THEN
@@ -156,4 +96,25 @@
         ! MLD: rho = rho(1) + zrho1                                     !
         ! ------------------------------------------------------------- !
+        zmaxdzT(:,:) = 0._wp
+        IF( nla10 > 1 ) THEN
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
+                 zrho0_3(ji,jj) = zztmp
+                 zrho0_1(ji,jj) = zztmp
+                 zhth(ji,jj) = zztmp
+              END DO
+           END DO
+        ELSE IF (iom_use("mlddzt")) THEN
+           DO jj = 1, jpj
+              DO ji = 1, jpi
+                 zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
+                 zhth(ji,jj) = zztmp
+              END DO
+           END DO
+        ELSE
+           zhth(:,:) = 0._wp
+        ENDIF
+
         DO jk = jpkm1, 2, -1   ! loop from bottom to 2
            DO jj = 1, jpj
@@ -164,9 +125,9 @@
                  zzdep = zzdep * tmask(ji,jj,1)
 
-                 IF( zztmp > zmaxdzT(ji,jj) ) THEN                        
+                 IF( zztmp > zmaxdzT(ji,jj) ) THEN
                     zmaxdzT(ji,jj) = zztmp   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz
                  ENDIF
 
-                 IF( nla10 > 1 ) THEN 
+                 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
@@ -179,5 +140,5 @@
 
         IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline
-        IF( nla10 > 1 ) THEN 
+        IF( nla10 > 1 ) THEN
            IF (iom_use("mldr0_3")) CALL iom_put( "mldr0_3", zrho0_3*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.03
            IF (iom_use("mldr0_1")) CALL iom_put( "mldr0_1", zrho0_1*tmask(:,:,1) )   ! MLD delta rho(surf) = 0.01
@@ -187,4 +148,16 @@
      IF (iom_use("mld_dt02").OR.iom_use("topthdep").OR.iom_use("mldr10_3").OR. &
           & iom_use("pycndep").OR.iom_use("tinv").OR.iom_use("depti")) THEN
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1)
+              zabs2   (ji,jj) = zztmp
+              ztm2    (ji,jj) = zztmp
+              zrho10_3(ji,jj) = zztmp
+              zpycn   (ji,jj) = zztmp
+           END DO
+        END DO
+        ztinv  (:,:) = 0._wp
+        zdepinv(:,:) = 0._wp
+
         IF (iom_use("pycndep")) THEN
            ! Preliminary computation
@@ -207,4 +180,6 @@
               END DO
            END DO
+        ELSE
+           zdelr(:,:) = 0._wp
         ENDIF
 
@@ -243,6 +218,6 @@
         IF (iom_use("mldr10_3")) CALL iom_put( "mldr10_3", zrho10_3*tmask(:,:,1)     )   ! MLD delta rho(10m) = 0.03
         IF (iom_use("pycndep")) CALL iom_put( "pycndep" , zpycn*tmask(:,:,1)        )   ! MLD delta rho equi. delta T(10m) = 0.2
-        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref) 
-        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref) 
+        IF (iom_use("tinv")) CALL iom_put( "tinv"    , ztinv*tmask(:,:,1)        )   ! max. temp. inv. (t10 ref)
+        IF (iom_use("depti")) CALL iom_put( "depti"   , zdepinv*tmask(:,:,1)      )   ! depth of max. temp. inv. (t10 ref)
      ENDIF
 
@@ -272,5 +247,5 @@
               !
               iid = ik20(ji,jj)
-              IF( iid /= 1 ) THEN 
+              IF( iid /= 1 ) THEN
                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
@@ -278,10 +253,10 @@
                       &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
                  zhd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
-              ELSE 
+              ELSE
                  zhd20(ji,jj) = 0._wp
               ENDIF
               !
               iid = ik28(ji,jj)
-              IF( iid /= 1 ) THEN 
+              IF( iid /= 1 ) THEN
                  zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
                       &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
@@ -289,5 +264,5 @@
                       &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
                  zhd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
-              ELSE 
+              ELSE
                  zhd28(ji,jj) = 0._wp
               ENDIF
@@ -306,11 +281,11 @@
         ilevel   = 0
         zthick_0 = 0._wp
-        DO jk = 1, jpkm1                      
+        DO jk = 1, jpkm1
            zthick_0 = zthick_0 + e3t_1d(jk)
            IF( zthick_0 < 300. )   ilevel = jk
         END DO
         ! surface boundary condition
-        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
-        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp                                   
+        IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   zhtc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)
+        ELSE                   ;   zthick(:,:) = 0._wp       ;   zhtc3(:,:) = 0._wp
         ENDIF
         ! integration down to ilevel
@@ -333,21 +308,38 @@
      ENDIF
      !
-      IF( ln_timing )   CALL timing_stop('dia_hth')
-     !
    END SUBROUTINE dia_hth
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                                       Empty module
-   !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
-CONTAINS
-   SUBROUTINE dia_hth( kt )         ! Empty routine
-      IMPLICIT NONE
-      INTEGER, INTENT( in ) :: kt
-      WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
-   END SUBROUTINE dia_hth
-#endif
-
-   !!======================================================================
+
+   SUBROUTINE dia_hth_init
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hth_init  ***
+      !!
+      !! ** Purpose: Initialization for ML and thermocline depths
+      !!
+      !! ** Action : Read namelist flag to permit or not extra Ml depths and thermocline depths
+      !!---------------------------------------------------------------------------
+      INTEGER ::   ierror, ios   ! local integer
+      !!
+      NAMELIST/namhth/ ln_diahth
+      !!---------------------------------------------------------------------------
+       !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics'
+         WRITE(numout,*) '~~~~~~~~~~~~ '
+      ENDIF
+      REWIND( numnam_ref )              ! Namelist namhth in reference namelist
+      READ  ( numnam_ref, namhth, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhth in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist namhth in configuration namelist
+      READ  ( numnam_cfg, namhth, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhth in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namhth )
+
+      IF(lwp) THEN
+         WRITE(numout,*) '   Namelist  namhth :'
+         WRITE(numout,*) '      check the heat and salt budgets (T) or not (F)       ln_diahth = ', ln_diahth
+      ENDIF
+      !
+   END SUBROUTINE dia_hth_init
 END MODULE diahth
Index: /NEMO/releases/release-4.0/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/releases/release-4.0/src/OCE/nemogcm.F90	(revision 11090)
+++ /NEMO/releases/release-4.0/src/OCE/nemogcm.F90	(revision 11091)
@@ -485,4 +485,5 @@
       IF( lk_diadct    )   CALL dia_dct_init    ! Sections tranports
                            CALL dia_hsb_init    ! heat content, salt content and volume budgets
+                           CALL dia_hth_init    ! extra ML depth diagnostics, thermocline depths
                            CALL     trd_init    ! Mixed-layer/Vorticity/Integral constraints trends
                            CALL dia_obs_init    ! Initialize observational data
Index: /NEMO/releases/release-4.0/src/OCE/step.F90
===================================================================
--- /NEMO/releases/release-4.0/src/OCE/step.F90	(revision 11090)
+++ /NEMO/releases/release-4.0/src/OCE/step.F90	(revision 11091)
@@ -212,9 +212,9 @@
       ! diagnostics and outputs
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats
-      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
-      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
-      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports
-                         CALL dia_ar5 ( kstp )        ! ar5 diag
+      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
+      IF( ln_diacfl  )   CALL dia_cfl( kstp )         ! Courant number diagnostics
+      IF( ln_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
+      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
+                         CALL dia_ar5( kstp )         ! ar5 diag
       IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
                          CALL dia_wri ( kstp )        ! ocean model: outputs
