Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/dia25h.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/dia25h.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/dia25h.F90	(revision 12178)
@@ -55,8 +55,8 @@
       REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
       READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist' )
       REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
       READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist' )
       IF(lwm) WRITE ( numond, nam_dia25h )
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diacfl.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diacfl.F90	(revision 12178)
@@ -29,9 +29,4 @@
    REAL(wp)              ::   rCu_max, rCv_max, rCw_max   ! associated run max Courant number 
 
-!!gm CAUTION: need to declare these arrays here, otherwise the calculation fails in multi-proc !
-!!gm          I don't understand why.
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
-!!gm end
-
    PUBLIC   dia_cfl       ! routine called by step.F90
    PUBLIC   dia_cfl_init  ! routine called by nemogcm
@@ -55,8 +50,8 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER                ::   ji, jj, jk                            ! dummy loop indices
-      REAL(wp)               ::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
-      INTEGER , DIMENSION(3) ::   iloc_u , iloc_v , iloc_w , iloc       ! workspace
-!!gm this does not work      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl         ! workspace
+      INTEGER                          ::   ji, jj, jk                       ! dummy loop indices
+      REAL(wp)                         ::   z2dt, zCu_max, zCv_max, zCw_max  ! local scalars
+      INTEGER , DIMENSION(3)           ::   iloc_u , iloc_v , iloc_w , iloc  ! workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zCu_cfl, zCv_cfl, zCw_cfl        ! workspace
       !!----------------------------------------------------------------------
       !
@@ -71,5 +66,5 @@
       DO jk = 1, jpk       ! calculate Courant numbers
          DO jj = 1, jpj
-            DO ji = 1, fs_jpim1   ! vector opt.
+            DO ji = 1, jpi
                zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
                zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
@@ -111,5 +106,5 @@
       !                    ! write out to file
       IF( lwp ) THEN
-         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(2x,i6,3x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
          WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
@@ -172,8 +167,4 @@
       rCw_max = 0._wp
       !
-!!gm required to work
-      ALLOCATE ( zCu_cfl(jpi,jpj,jpk), zCv_cfl(jpi,jpj,jpk), zCw_cfl(jpi,jpj,jpk) )
-!!gm end
-      !      
    END SUBROUTINE dia_cfl_init
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diadct.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diadct.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diadct.F90	(revision 12178)
@@ -11,8 +11,6 @@
    !!            3.4  ! 09/2011 (C Bricaud)
    !!----------------------------------------------------------------------
-#if defined key_diadct
-   !!----------------------------------------------------------------------
-   !!   'key_diadct' :
-   !!----------------------------------------------------------------------
+   !! does not work with agrif
+#if ! defined key_agrif
    !!----------------------------------------------------------------------
    !!   dia_dct      :  Compute the transport through a sec.
@@ -42,16 +40,11 @@
 
    PUBLIC   dia_dct      ! routine called by step.F90
-   PUBLIC   dia_dct_init ! routine called by opa.F90
-   PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90 
-   PRIVATE  readsec
-   PRIVATE  removepoints
-   PRIVATE  transport
-   PRIVATE  dia_dct_wri
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
-
-   INTEGER :: nn_dct        ! Frequency of computation
-   INTEGER :: nn_dctwri     ! Frequency of output
-   INTEGER :: nn_secdebug   ! Number of the section to debug
+   PUBLIC   dia_dct_init ! routine called by nemogcm.F90
+
+   !                         !!** namelist variables **
+   LOGICAL, PUBLIC ::   ln_diadct     !: Calculate transport thru a section or not
+   INTEGER         ::   nn_dct        !  Frequency of computation
+   INTEGER         ::   nn_dctwri     !  Frequency of output
+   INTEGER         ::   nn_secdebug   !  Number of the section to debug
    
    INTEGER, PARAMETER :: nb_class_max  = 10
@@ -104,19 +97,16 @@
 CONTAINS
  
-  INTEGER FUNCTION diadct_alloc() 
-     !!---------------------------------------------------------------------- 
-     !!                   ***  FUNCTION diadct_alloc  *** 
-     !!---------------------------------------------------------------------- 
-     INTEGER :: ierr(2) 
-     !!---------------------------------------------------------------------- 
-
-     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 
-     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) ) 
-
-     diadct_alloc = MAXVAL( ierr ) 
-     IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
- 
-  END FUNCTION diadct_alloc 
-
+   INTEGER FUNCTION diadct_alloc() 
+      !!---------------------------------------------------------------------- 
+      !!                   ***  FUNCTION diadct_alloc  *** 
+      !!---------------------------------------------------------------------- 
+
+      ALLOCATE( transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), &
+         &      transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=diadct_alloc ) 
+
+      CALL mpp_sum( 'diadct', diadct_alloc ) 
+      IF( diadct_alloc /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' ) 
+
+   END FUNCTION diadct_alloc
 
    SUBROUTINE dia_dct_init
@@ -130,15 +120,15 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
+      NAMELIST/nam_diadct/ln_diadct, nn_dct, nn_dctwri, nn_secdebug
       !!---------------------------------------------------------------------
 
-     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
-     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
-901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
-
-     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
-     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
-902  IF( ios >  0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
-     IF(lwm) WRITE ( numond, namdct )
+     REWIND( numnam_ref )              ! Namelist nam_diadct in reference namelist : Diagnostic: transport through sections
+     READ  ( numnam_ref, nam_diadct, IOSTAT = ios, ERR = 901)
+901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadct in reference namelist' )
+
+     REWIND( numnam_cfg )              ! Namelist nam_diadct in configuration namelist : Diagnostic: transport through sections
+     READ  ( numnam_cfg, nam_diadct, IOSTAT = ios, ERR = 902 )
+902  IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diadct in configuration namelist' )
+     IF(lwm) WRITE ( numond, nam_diadct )
 
      IF( lwp ) THEN
@@ -146,6 +136,7 @@
         WRITE(numout,*) "diadct_init: compute transports through sections "
         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
-        WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
-        WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
+        WRITE(numout,*) "       Calculate transport thru sections: ln_diadct = ", ln_diadct
+        WRITE(numout,*) "       Frequency of computation:          nn_dct    = ", nn_dct
+        WRITE(numout,*) "       Frequency of write:                nn_dctwri = ", nn_dctwri
 
         IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
@@ -155,23 +146,29 @@
         ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
         ENDIF
-
+     ENDIF
+
+     IF( ln_diadct ) THEN
+        ! control
         IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
-          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
-
+           &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
+
+        ! allocate dia_dct arrays
+        IF( diadct_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'diadct_alloc: failed to allocate arrays' )
+
+        !Read section_ijglobal.diadct
+        CALL readsec
+
+        !open output file
+        IF( lwm ) THEN
+           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+        ENDIF
+
+        ! Initialise arrays to zero 
+        transports_3d(:,:,:,:)=0.0 
+        transports_2d(:,:,:)  =0.0 
+        !
      ENDIF
-
-     !Read section_ijglobal.diadct
-     CALL readsec
-
-     !open output file
-     IF( lwm ) THEN
-        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-        CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-     ENDIF
-
-     ! Initialise arrays to zero 
-     transports_3d(:,:,:,:)=0.0 
-     transports_2d(:,:,:)  =0.0 
      !
   END SUBROUTINE dia_dct_init
@@ -1241,21 +1238,16 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option :                                       Dummy module
+   !!   Dummy module                                             
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
-   PUBLIC 
-   !! $Id$
+   LOGICAL, PUBLIC ::   ln_diadct = .FALSE.
 CONTAINS
-
-   SUBROUTINE dia_dct_init          ! Dummy routine
+   SUBROUTINE dia_dct_init
       IMPLICIT NONE
-      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?'
    END SUBROUTINE dia_dct_init
-
-   SUBROUTINE dia_dct( kt )         ! Dummy routine
+   SUBROUTINE dia_dct( kt )
       IMPLICIT NONE
-      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
-      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
+      INTEGER, INTENT(in) ::   kt
    END SUBROUTINE dia_dct
+   !
 #endif
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaharm.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaharm.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaharm.F90	(revision 12178)
@@ -5,8 +5,4 @@
    !!======================================================================
    !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code
-   !!----------------------------------------------------------------------
-#if defined key_diaharm
-   !!----------------------------------------------------------------------
-   !!   'key_diaharm'
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
@@ -26,6 +22,4 @@
    IMPLICIT NONE
    PRIVATE
-
-   LOGICAL, PUBLIC, PARAMETER :: lk_diaharm  = .TRUE.
    
    INTEGER, PARAMETER :: jpincomax    = 2.*jpmax_harmo
@@ -33,8 +27,9 @@
 
    !                         !!** namelist variables **
-   INTEGER ::   nit000_han    ! First time step used for harmonic analysis
-   INTEGER ::   nitend_han    ! Last time step used for harmonic analysis
-   INTEGER ::   nstep_han     ! Time step frequency for harmonic analysis
-   INTEGER ::   nb_ana        ! Number of harmonics to analyse
+   LOGICAL, PUBLIC ::   ln_diaharm    ! Choose tidal harmonic output or not
+   INTEGER         ::   nit000_han    ! First time step used for harmonic analysis
+   INTEGER         ::   nitend_han    ! Last time step used for harmonic analysis
+   INTEGER         ::   nstep_han     ! Time step frequency for harmonic analysis
+   INTEGER         ::   nb_ana        ! Number of harmonics to analyse
 
    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name
@@ -53,5 +48,6 @@
    CHARACTER (LEN=4), DIMENSION(jpmax_harmo) ::   tname   ! Names of tidal constituents ('M2', 'K1',...)
 
-   PUBLIC   dia_harm   ! routine called by step.F90
+   PUBLIC   dia_harm        ! routine called by step.F90
+   PUBLIC   dia_harm_init   ! routine called by nemogcm.F90
 
    !!----------------------------------------------------------------------
@@ -71,8 +67,8 @@
       !!
       !!--------------------------------------------------------------------
-      INTEGER :: jh, nhan, jk, ji
+      INTEGER ::   jh, nhan, ji
       INTEGER ::   ios                 ! Local integer output status for namelist read
 
-      NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
+      NAMELIST/nam_diaharm/ ln_diaharm, nit000_han, nitend_han, nstep_han, tname
       !!----------------------------------------------------------------------
 
@@ -83,79 +79,84 @@
       ENDIF
       !
-      IF( .NOT. ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
-      !
-      CALL tide_init_Wave
-      !
       REWIND( numnam_ref )              ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
       READ  ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_diaharm in reference namelist' )
       REWIND( numnam_cfg )              ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
       READ  ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist' )
       IF(lwm) WRITE ( numond, nam_diaharm )
       !
       IF(lwp) THEN
-         WRITE(numout,*) 'First time step used for analysis:  nit000_han= ', nit000_han
-         WRITE(numout,*) 'Last  time step used for analysis:  nitend_han= ', nitend_han
-         WRITE(numout,*) 'Time step frequency for harmonic analysis:  nstep_han= ', nstep_han
+         WRITE(numout,*) 'Tidal diagnostics = ', ln_diaharm
+         WRITE(numout,*) '   First time step used for analysis:         nit000_han= ', nit000_han
+         WRITE(numout,*) '   Last  time step used for analysis:         nitend_han= ', nitend_han
+         WRITE(numout,*) '   Time step frequency for harmonic analysis: nstep_han = ', nstep_han
       ENDIF
 
-      ! Basic checks on harmonic analysis time window:
-      ! ----------------------------------------------
-      IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
-         &                                       ' restart capability not implemented' )
-      IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
-         &                                       'restart capability not implemented' )
-
-      IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
-         &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
-
-      nb_ana = 0
-      DO jk=1,jpmax_harmo
-         DO ji=1,jpmax_harmo
-            IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
-               nb_ana=nb_ana+1
-            ENDIF
-         END DO
-      END DO
-      !
-      IF(lwp) THEN
-         WRITE(numout,*) '        Namelist nam_diaharm'
-         WRITE(numout,*) '        nb_ana    = ', nb_ana
-         CALL flush(numout)
+      IF( ln_diaharm .AND. .NOT.ln_tide )   CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
+
+      IF( ln_diaharm ) THEN
+
+         CALL tide_init_Wave
+         !
+         ! Basic checks on harmonic analysis time window:
+         ! ----------------------------------------------
+         IF( nit000 > nit000_han )   CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000',   &
+            &                                       ' restart capability not implemented' )
+         IF( nitend < nitend_han )   CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend',   &
+            &                                       'restart capability not implemented' )
+
+         IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 )   &
+            &                        CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
+         !
+         nb_ana = 0
+         DO jh=1,jpmax_harmo
+            DO ji=1,jpmax_harmo
+               IF(TRIM(tname(jh)) == Wave(ji)%cname_tide) THEN
+                  nb_ana=nb_ana+1
+               ENDIF
+            END DO
+         END DO
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) '        Namelist nam_diaharm'
+            WRITE(numout,*) '        nb_ana    = ', nb_ana
+            CALL flush(numout)
+         ENDIF
+         !
+         IF (nb_ana > jpmax_harmo) THEN
+            WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
+            WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
+            CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
+         ENDIF
+
+         ALLOCATE(name    (nb_ana))
+         DO jh=1,nb_ana
+            DO ji=1,jpmax_harmo
+               IF (TRIM(tname(jh)) ==  Wave(ji)%cname_tide) THEN
+                  name(jh) = ji
+                  EXIT
+               END IF
+            END DO
+         END DO
+
+         ! Initialize frequency array:
+         ! ---------------------------
+         ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
+
+         CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
+
+         IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
+
+         DO jh = 1, nb_ana
+            IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
+         END DO
+
+         ! Initialize temporary arrays:
+         ! ----------------------------
+         ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
+         ana_temp(:,:,:,:) = 0._wp
+
       ENDIF
-      !
-      IF (nb_ana > jpmax_harmo) THEN
-         WRITE(ctmp1,*) ' nb_ana must be lower than jpmax_harmo'
-         WRITE(ctmp2,*) ' jpmax_harmo= ', jpmax_harmo
-         CALL ctl_stop( 'dia_harm_init', ctmp1, ctmp2 )
-      ENDIF
-
-      ALLOCATE(name    (nb_ana))
-      DO jk=1,nb_ana
-       DO ji=1,jpmax_harmo
-          IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN
-             name(jk) = ji
-             EXIT
-          END IF
-       END DO
-      END DO
-
-      ! Initialize frequency array:
-      ! ---------------------------
-      ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
-
-      CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
-
-      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency '
-
-      DO jh = 1, nb_ana
-        IF(lwp) WRITE(numout,*) '                    : ',tname(jh),' ',ana_freq(jh)
-      END DO
-
-      ! Initialize temporary arrays:
-      ! ----------------------------
-      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
-      ana_temp(:,:,:,:) = 0._wp
 
    END SUBROUTINE dia_harm_init
@@ -177,6 +178,4 @@
       !!--------------------------------------------------------------------
       IF( ln_timing )   CALL timing_start('dia_harm')
-      !
-      IF( kt == nit000 )   CALL dia_harm_init
       !
       IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
@@ -422,5 +421,5 @@
       INTEGER, INTENT(in) ::   init 
       !
-      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
+      INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
       REAL(wp)                        :: zval1, zval2, zx1
       REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2
@@ -434,11 +433,11 @@
          ztmp3(:,:) = 0._wp
          !
-         DO jk1_sd = 1, nsparse
-            DO jk2_sd = 1, nsparse
-               nisparse(jk2_sd) = nisparse(jk2_sd)
-               njsparse(jk2_sd) = njsparse(jk2_sd)
-               IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
-                  ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  &
-                     &                                     + valuesparse(jk1_sd)*valuesparse(jk2_sd)
+         DO jh1_sd = 1, nsparse
+            DO jh2_sd = 1, nsparse
+               nisparse(jh2_sd) = nisparse(jh2_sd)
+               njsparse(jh2_sd) = njsparse(jh2_sd)
+               IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
+                  ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd))  &
+                     &                                     + valuesparse(jh1_sd)*valuesparse(jh2_sd)
                ENDIF
             END DO
@@ -515,16 +514,4 @@
    END SUBROUTINE SUR_DETERMINE
 
-#else
-   !!----------------------------------------------------------------------
-   !!   Default case :   Empty module
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
-CONTAINS
-   SUBROUTINE dia_harm ( kt )     ! Empty routine
-      INTEGER, INTENT( IN ) :: kt  
-      WRITE(*,*) 'dia_harm: you should not have seen this print'
-   END SUBROUTINE dia_harm
-#endif
-
    !!======================================================================
 END MODULE diaharm
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahsb.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahsb.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahsb.F90	(revision 12178)
@@ -362,8 +362,8 @@
       REWIND( numnam_ref )              ! Namelist namhsb in reference namelist
       READ  ( numnam_ref, namhsb, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namhsb in reference namelist' )
       REWIND( numnam_cfg )              ! Namelist namhsb in configuration namelist
       READ  ( numnam_cfg, namhsb, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namhsb in configuration namelist' )
       IF(lwm) WRITE( numond, namhsb )
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahth.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahth.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diahth.F90	(revision 12178)
@@ -5,9 +5,13 @@
    !!======================================================================
    !! 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.
    !!----------------------------------------------------------------------
    !!   dia_hth      : Compute varius diagnostics associated with the mixed layer
@@ -20,4 +24,5 @@
    USE lib_mpp         ! MPP library
    USE iom             ! I/O library
+   USE timing          ! preformance summary
 
    IMPLICIT NONE
@@ -25,7 +30,13 @@
 
    PUBLIC   dia_hth       ! routine called by step.F90
-   PUBLIC   dia_hth_init  ! routine called by nemogcm.F90
-
-   LOGICAL, PUBLIC ::   ll_diahth   !: Compute further diagnostics of ML and thermocline depth
+   PUBLIC   dia_hth_alloc ! routine called by nemogcm.F90
+
+   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, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hth    !: depth of the max vertical temperature gradient [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd20   !: depth of 20 C isotherm                         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hd28   !: depth of 28 C isotherm                         [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htc3   !: heat content of first 300 m                    [W]
 
    !!----------------------------------------------------------------------
@@ -36,299 +47,306 @@
 CONTAINS
 
+   FUNCTION dia_hth_alloc()
+      !!---------------------------------------------------------------------
+      INTEGER :: dia_hth_alloc
+      !!---------------------------------------------------------------------
+      !
+      ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd28(jpi,jpj), htc3(jpi,jpj), STAT=dia_hth_alloc )
+      !
+      CALL mpp_sum ( 'diahth', dia_hth_alloc )
+      IF(dia_hth_alloc /= 0)   CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' )
+      !
+   END FUNCTION dia_hth_alloc
+
 
    SUBROUTINE dia_hth( kt )
-     !!---------------------------------------------------------------------
-     !!                  ***  ROUTINE dia_hth  ***
-     !!
-     !! ** 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 :
-     !!-------------------------------------------------------------------
-     INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-     !!
-     INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
-     INTEGER                          ::   iid, ilevel           ! temporary integers
-     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
-   ! 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
-        ! ------------------------------------------------------------- !
-        ! thermocline depth: strongest vertical gradient of temperature !
-        ! turbocline depth (mixing layer depth): avt = zavt5            !
-        ! MLD: rho = rho(1) + zrho3                                     !
-        ! 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
-              DO ji = 1, jpi
-                 !
-                 zzdep = gdepw_n(ji,jj,jk)
-                 zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / 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   ;   zhth    (ji,jj) = zzdep                ! max and depth of dT/dz
-                 ENDIF
-
-                 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
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_hth  ***
+      !!
+      !! ** 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 : 
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!
+      INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
+      INTEGER                          ::   iid, ilevel           ! temporary integers
+      INTEGER, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   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), 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
+      !!----------------------------------------------------------------------
+      IF( ln_timing )   CALL timing_start('dia_hth')
+
+      IF( kt == nit000 ) THEN
+         !                                      ! allocate dia_hth array
+         IF( dia_hth_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_hth : unable to allocate standard arrays' )
+
+         IF(.NOT. ALLOCATED(ik20) ) THEN
+            ALLOCATE(ik20(jpi,jpj), ik28(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)
+            CALL mpp_sum('diahth', 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
+
+      ! 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) 
+            hth     (ji,jj) = zztmp
+            zabs2   (ji,jj) = zztmp
+            ztm2    (ji,jj) = zztmp
+            zrho10_3(ji,jj) = zztmp
+            zpycn   (ji,jj) = zztmp
         END DO
-
-        IF (iom_use("mlddzt")) CALL iom_put( "mlddzt", zhth*tmask(:,:,1) )            ! depth of the thermocline
-        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
-        ENDIF
-     ENDIF
-
-     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
-           ! 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 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
-                         &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
-                         &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
-                    zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
-                         &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
-                    zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
-                    zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
-                    zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
-                    zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
-                 ELSE
-                    zdelr(ji,jj) = 0._wp
-                 ENDIF
-              END DO
-           END DO
-        ELSE
-           zdelr(:,:) = 0._wp
-        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 = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
-                 !
-                 zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - 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
-
-        IF (iom_use("mld_dt02")) CALL iom_put( "mld_dt02", zabs2*tmask(:,:,1)        )   ! MLD abs(delta t) - 0.2
-        IF (iom_use("topthdep")) CALL iom_put( "topthdep", ztm2*tmask(:,:,1)         )   ! T(10) - 0.2
-        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)
-     ENDIF
-
-     IF(iom_use("20d") .OR. iom_use("28d")) THEN
-        ! ----------------------------------- !
-        ! 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 = tsn(ji,jj,jk,jp_tem)
-                 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
-              !
-              zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
-              !
-              iid = ik20(ji,jj)
-              IF( iid /= 1 ) THEN
-                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                      &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
-                      &  / ( 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
-                 zhd20(ji,jj) = 0._wp
-              ENDIF
-              !
-              iid = ik28(ji,jj)
-              IF( iid /= 1 ) THEN
-                 zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                      &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                      &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
-                      &  / (  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
-                 zhd28(ji,jj) = 0._wp
-              ENDIF
-
-           END DO
-        END DO
-        CALL iom_put( "20d", zhd20 )   ! depth of the 20 isotherm
-        CALL iom_put( "28d", zhd28 )   ! depth of the 28 isotherm
-     ENDIF
-
-     ! ----------------------------- !
-     !  Heat content of first 300 m  !
-     ! ----------------------------- !
-     IF (iom_use("hc300")) THEN
-        ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
-        ilevel   = 0
-        zthick_0 = 0._wp
-        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
-        ENDIF
-        ! integration down to ilevel
-        DO jk = 1, ilevel
-           zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
-           zhtc3  (:,:) = zhtc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
-        END DO
-        ! deepest layer
-        zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              zhtc3(ji,jj) = zhtc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
-                   &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
-           END DO
-        END DO
-        ! from temperature to heat contain
-        zcoef = rau0 * rcp
-        zhtc3(:,:) = zcoef * zhtc3(:,:)
-        CALL iom_put( "hc300", zhtc3*tmask(:,:,1) )      ! first 300m heat content
-     ENDIF
-     !
+      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
+      
+      ! 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 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
+                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
+                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
+               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
+                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
+               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
+               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
+               zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
+               zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
+            ELSE
+               zdelr(ji,jj) = 0._wp
+            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 = gdepw_n(ji,jj,jk)
+               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / 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( 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
+      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 = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
+               !
+               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - 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_dt02", 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( "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 = tsn(ji,jj,jk,jp_tem)
+               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
+            !
+            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
+            !
+            iid = ik20(ji,jj)
+            IF( iid /= 1 ) THEN 
+               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
+                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
+                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
+                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
+               hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
+            ELSE 
+               hd20(ji,jj) = 0._wp
+            ENDIF
+            !
+            iid = ik28(ji,jj)
+            IF( iid /= 1 ) THEN 
+               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
+                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
+                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
+                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
+               hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
+            ELSE 
+               hd28(ji,jj) = 0._wp
+            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
+
+      ! ----------------------------- !
+      !  Heat content of first 300 m  !
+      ! ----------------------------- !
+
+      ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...)
+      ilevel   = 0
+      zthick_0 = 0._wp
+      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(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
+      ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
+      ENDIF
+      ! integration down to ilevel
+      DO jk = 1, ilevel
+         zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
+         htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
+      END DO
+      ! deepest layer
+      zthick(:,:) = 300. - zthick(:,:)   !   remaining thickness to reach 300m
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
+               &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
+         END DO
+      END DO
+      ! from temperature to heat contain
+      zcoef = rau0 * rcp
+      htc3(:,:) = zcoef * htc3(:,:)
+      CALL iom_put( "hc300", htc3 )      ! first 300m heat content
+      !
+      IF( ln_timing )   CALL timing_stop('dia_hth')
+      !
    END SUBROUTINE dia_hth
 
-
-   SUBROUTINE dia_hth_init
-      !!---------------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_hth_init  ***
-      !!
-      !! ** Purpose: Initialization for ML and thermocline depths
-      !!
-      !! ** Action : If any upper ocean diagnostic required by xml file, set in dia_hth
-      !!---------------------------------------------------------------------------
-       !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dia_hth_init : heat and salt budgets diagnostics'
-         WRITE(numout,*) '~~~~~~~~~~~~ '
-      ENDIF
-      ll_diahth = 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").OR. &
-           & iom_use("20d") .OR. iom_use("28d") .OR. iom_use("hc300")
-      IF(lwp) THEN
-         WRITE(numout,*) '   output upper ocean diagnostics (T) or not (F)       ll_diahth = ', ll_diahth
-      ENDIF
-      !
-   END SUBROUTINE dia_hth_init
+#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
+
+   !!======================================================================
 END MODULE diahth
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaptr.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diaptr.F90	(revision 12178)
@@ -393,9 +393,9 @@
       REWIND( numnam_ref )              ! Namelist namptr in reference namelist : Poleward transport
       READ  ( numnam_ref, namptr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist', lwp )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namptr in reference namelist' )
 
       REWIND( numnam_cfg )              ! Namelist namptr in configuration namelist : Poleward transport
       READ  ( numnam_cfg, namptr, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist', lwp )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namptr in configuration namelist' )
       IF(lwm) WRITE ( numond, namptr )
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diatmb.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diatmb.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diatmb.F90	(revision 12178)
@@ -43,9 +43,9 @@
       REWIND( numnam_ref )              ! Read Namelist nam_diatmb in reference namelist : TMB diagnostics
       READ  ( numnam_ref, nam_diatmb, IOSTAT=ios, ERR= 901 )
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist', lwp )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diatmb in reference namelist' )
  
       REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
       READ  ( numnam_cfg, nam_diatmb, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist', lwp )
+902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_diatmb in configuration namelist' )
       IF(lwm) WRITE ( numond, nam_diatmb )
 
Index: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90	(revision 11143)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90	(revision 12178)
@@ -43,5 +43,4 @@
    USE zdfdrg         ! ocean vertical physics: top/bottom friction
    USE zdfmxl         ! mixed layer
-   USE zdfosm         ! mixed layer
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -211,4 +210,6 @@
       ENDIF
 
+      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
+      !
       CALL iom_put( "woce", wn )                   ! vertical velocity
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
@@ -221,4 +222,6 @@
          IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
       ENDIF
+      !
+      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
 
       CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
@@ -427,5 +430,5 @@
       !!      define all the NETCDF files and fields
       !!      At each time step call histdef to compute the mean if ncessary
-      !!      Each nwrite time step, output the instantaneous or mean fields
+      !!      Each nn_write time step, output the instantaneous or mean fields
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
@@ -443,6 +446,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
       !!----------------------------------------------------------------------
-      ! 
-      IF( ln_timing )   CALL timing_start('dia_wri')
       !
       IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
@@ -451,4 +452,8 @@
       ENDIF
       !
+      IF( nn_write == -1 )   RETURN   ! we will never do any output
+      ! 
+      IF( ln_timing )   CALL timing_start('dia_wri')
+      !
       ! 0. Initialisation
       ! -----------------
@@ -460,5 +465,5 @@
       clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
 #if defined key_diainstant
-      zsto = nwrite * rdt
+      zsto = nn_write * rdt
       clop = "inst("//TRIM(clop)//")"
 #else
@@ -466,5 +471,5 @@
       clop = "ave("//TRIM(clop)//")"
 #endif
-      zout = nwrite * rdt
+      zout = nn_write * rdt
       zmax = ( nitend - nit000 + 1 ) * rdt
 
@@ -497,5 +502,5 @@
          ! WRITE root name in date.file for use by postpro
          IF(lwp) THEN
-            CALL dia_nam( clhstnam, nwrite,' ' )
+            CALL dia_nam( clhstnam, nn_write,' ' )
             CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
             WRITE(inum,*) clhstnam
@@ -505,5 +510,5 @@
          ! Define the T grid FILE ( nid_T )
 
-         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
+         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
          IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
@@ -541,5 +546,5 @@
          ! Define the U grid FILE ( nid_U )
 
-         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
+         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
          IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
          CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
@@ -554,5 +559,5 @@
          ! Define the V grid FILE ( nid_V )
 
-         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
+         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
          IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
          CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
@@ -567,5 +572,5 @@
          ! Define the W grid FILE ( nid_W )
 
-         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
+         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
          IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
          CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
@@ -658,5 +663,5 @@
          ENDIF
 
-         IF( .NOT. ln_cpl ) THEN
+         IF( ln_ssr ) THEN
             CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -666,14 +671,5 @@
                &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          ENDIF
-
-         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
-            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
-               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         ENDIF
-         
+       
          clmx ="l_max(only(x))"    ! max index on a period
 !         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX 
@@ -751,5 +747,5 @@
       ! donne le nombre d'elements, et ndex la liste des indices a sortir
 
-      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN 
+      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN 
          WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
          WRITE(numout,*) '~~~~~~ '
@@ -815,14 +811,8 @@
       ENDIF
 
-      IF( .NOT. ln_cpl ) THEN
+      IF( ln_ssr ) THEN
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
-         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
-      ENDIF
-      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
-         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
-         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
@@ -843,5 +833,9 @@
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
-      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
+      IF( ln_zad_Aimp ) THEN
+         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current
+      ELSE
+         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current
+      ENDIF
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
       CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
@@ -904,5 +898,9 @@
       CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
       CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
-      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
+      IF( ln_zad_Aimp ) THEN
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
+      ELSE
+         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
+      ENDIF
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
@@ -928,22 +926,4 @@
          CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
       ENDIF
-
-      IF( ln_zdfosm ) THEN
-         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)   )    ! now boundary-layer depth
-         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)    )    ! now mixed-layer depth
-         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask       )    ! w-level diffusion
-         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask       )    ! now w-level viscosity
-         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask       )    ! non-local t forcing
-         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask       )    ! non-local s forcing
-         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*wmask       )    ! non-local u forcing
-         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamu*wmask       )    ! non-local v forcing
-      ENDIF
-     !     ! CALL histwrite( id_i, "zla", kt, zla*tmask(:,:,1)   , jpi*jpj, idex)         ! now Langmuir #
-     !     ! CALL histwrite( id_i, "zvstr", kt, zvstr*tmask(:,:,1)   , jpi*jpj, idex)     ! now mixed velocity scale
-     !     ! CALL histwrite( id_i, "zustar", kt, zustar*tmask(:,:,1)   , jpi*jpj, idex)   ! now friction velocity scale
-     !     ! CALL histwrite( id_i, "zwstrl", kt, zwstrl*tmask(:,:,1)   , jpi*jpj, idex)   ! now Langmuir velocity scale
-     !     ! CALL histwrite( id_i, "zwstrc", kt, zwstrc*tmask(:,:,1)   , jpi*jpj, idex)   ! now convective velocity scale
-     !     ! CALL histwrite( id_i, "zwb_ent", kt, zwb_ent*tmask(:,:,1)   , jpi*jpj, idex) ! now upward turb buoyancy entrainment flux
-     !     ! CALL histwrite( id_i, "zdb_bl", kt, zdb_bl*tmask(:,:,1)   , jpi*jpj, idex)   ! now db at ml base
  
 #if defined key_si3
