Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/dia25h.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/dia25h.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/dia25h.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diacfl.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diacfl.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diacfl.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diadct.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diadct.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diadct.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaharm.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaharm.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaharm.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diahsb.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diahsb.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diahsb.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaptr.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaptr.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diaptr.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diatmb.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diatmb.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diatmb.F90	(revision 11573)
@@ -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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diawri.F90	(revision 11238)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/OCE/DIA/diawri.F90	(revision 11573)
@@ -210,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
@@ -220,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.
@@ -426,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
@@ -442,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  ==!
@@ -450,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
       ! -----------------
@@ -459,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
@@ -465,5 +471,5 @@
       clop = "ave("//TRIM(clop)//")"
 #endif
-      zout = nwrite * rdt
+      zout = nn_write * rdt
       zmax = ( nitend - nit000 + 1 ) * rdt
 
@@ -496,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
@@ -504,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
@@ -540,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
@@ -553,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
@@ -566,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
@@ -657,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 )
@@ -665,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 
@@ -750,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,*) '~~~~~~ '
@@ -814,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
@@ -842,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.
@@ -903,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
