Changeset 9863 for NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE

Timestamp:

2018-06-30T12:51:02+02:00 (6 years ago)

Author:

gm

Message:

#1911 (ENHANCE-04): simplified implementation of the Euler stepping at nit000

Location:

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE

Files:

• ASM/asminc.F90 (modified) (6 diffs)
• DIA/diacfl.F90 (modified) (4 diffs)
• DOM/dom_oce.F90 (modified) (2 diffs)
• DOM/domain.F90 (modified) (4 diffs)
• DOM/domvvl.F90 (modified) (27 diffs)
• DOM/iscplrst.F90 (modified) (1 diff)
• DOM/istate.F90 (modified) (1 diff)
• DOM/restart.F90 (moved) (moved from NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/IOM/restart.F90) (11 diffs)
• DYN/dynnxt.F90 (modified) (9 diffs)
• DYN/dynspg.F90 (modified) (1 diff)
• DYN/dynspg_ts.F90 (modified) (13 diffs)
• DYN/dynzdf.F90 (modified) (21 diffs)
• DYN/sshwzv.F90 (modified) (9 diffs)
• DYN/wet_dry.F90 (modified) (7 diffs)
• TRA/traadv.F90 (modified) (1 diff)
• TRA/traldf.F90 (modified) (3 diffs)
• TRA/traldf_iso.F90 (modified) (3 diffs)
• TRA/traldf_triad.F90 (modified) (3 diffs)
• TRA/tranpc.F90 (modified) (2 diffs)
• TRA/tranxt.F90 (modified) (4 diffs)
• TRA/traqsr.F90 (modified) (1 diff)
• TRA/trazdf.F90 (modified) (3 diffs)
• TRD/trdtra.F90 (modified) (1 diff)
• nemogcm.F90 (modified) (1 diff)
• step.F90 (modified) (2 diffs)

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/ASM/asminc.F90

@@ -491 +491 @@
       ENDIF
       !
-      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', neuler
+      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', ln_1st_euler
       !
       IF( lk_asminc ) THEN                            !==  data assimilation  ==!
@@ -579 +579 @@
          IF ( kt == nitdin_r ) THEN
             !
-            neuler = 0  ! Force Euler forward step
+            l_1st_euler = .TRUE.  ! Force Euler forward step
             !
             ! Initialize the now fields with the background + increment
@@ -677 +677 @@
          IF ( kt == nitdin_r ) THEN
             !
-            neuler = 0                    ! Force Euler forward step
+            l_1st_euler = .TRUE.          ! Force Euler forward step
             !
             ! Initialize the now fields with the background + increment
@@ -752 +752 @@
          IF ( kt == nitdin_r ) THEN
             !
-            neuler = 0                                   ! Force Euler forward step
+            l_1st_euler = .TRUE.                         ! Force Euler forward step
             !
             sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
@@ -758 +758 @@
             sshb(:,:) = sshn(:,:)                        ! Update before fields
             e3t_b(:,:,:) = e3t_n(:,:,:)
-!!gm why not e3u_b, e3v_b, gdept_b ????
+            
+!!gm BUG :   missing the update of all other scale factors (e3u e3v e3w  etc... _n and _b) 
+!!           see dom_vvl_init 
             !
             DEALLOCATE( ssh_bkg    )
@@ -894 +896 @@
          IF ( kt == nitdin_r ) THEN
             !
-            neuler = 0                    ! Force Euler forward step
+            l_1st_euler = .TRUE.             ! Force Euler forward step
             !
             ! Sea-ice : SI3 case

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DIA/diacfl.F90

@@ -55 +55 @@
       !
       INTEGER :: ji, jj, jk   ! dummy loop indices
-      REAL(wp)::   z2dt, zCu_max, zCv_max, zCw_max       ! local scalars
+      REAL(wp)::   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
@@ -62 +62 @@
       IF( ln_timing )   CALL timing_start('dia_cfl')
       !
-      !                       ! setup timestep multiplier to account for initial Eulerian timestep
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;    z2dt = rdt
-      ELSE                                        ;    z2dt = rdt * 2._wp
-      ENDIF
-      !
       !                
       DO jk = 1, jpk       ! calculate Courant numbers
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-               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
-               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction
+               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * r2dt / e1u  (ji,jj)      ! for i-direction
+               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * r2dt / e2v  (ji,jj)      ! for j-direction
+               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * r2dt / e3w_n(ji,jj,jk)   ! for k-direction
             END DO
          END DO         
@@ -120 +115 @@
          WRITE(numcfl,*) '******************************************'
          WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
-         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCu_max
          WRITE(numcfl,*) '******************************************'
          WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
-         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCv_max
          WRITE(numcfl,*) '******************************************'
          WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
-         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max
+         WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', r2dt/rCw_max
          CLOSE( numcfl ) 
          !
@@ -133 +128 @@
          WRITE(numout,*) 'dia_cfl : Maximum Courant number information for the run '
          WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', z2dt/rCu_max
-         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', z2dt/rCv_max
-         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', z2dt/rCw_max
+         WRITE(numout,*) '   Max Cu = ', rCu_max, ' at (i,j,k) = (',nCu_loc(1),nCu_loc(2),nCu_loc(3),') => dt/C = ', r2dt/rCu_max
+         WRITE(numout,*) '   Max Cv = ', rCv_max, ' at (i,j,k) = (',nCv_loc(1),nCv_loc(2),nCv_loc(3),') => dt/C = ', r2dt/rCv_max
+         WRITE(numout,*) '   Max Cw = ', rCw_max, ' at (i,j,k) = (',nCw_loc(1),nCw_loc(2),nCw_loc(3),') => dt/C = ', r2dt/rCw_max
       ENDIF
       !

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/dom_oce.F90

@@ -35 +35 @@
    REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer
    REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter
-   INTEGER , PUBLIC ::   nn_euler       !: =0 start with forward time step or not (=1)
+   LOGICAL , PUBLIC ::   ln_1st_euler   !: =0 start with forward time step or not (=1)
    LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet
    LOGICAL , PUBLIC ::   ln_crs         !: Apply grid coarsening to dynamical model output or online passive tracers
@@ -57 +57 @@
    !                                   !! old non-DOCTOR names still used in the model
    REAL(wp), PUBLIC ::   atfp           !: asselin time filter parameter
-   REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer
+   REAL(wp), PUBLIC ::   rdt            !: time step for the dynamics and tracer (=rn_rdt)
 
    !                                   !!! associated variables
-   INTEGER , PUBLIC ::   neuler         !: restart euler forward option (0=Euler)
-   REAL(wp), PUBLIC ::   r2dt           !: = 2*rdt except at nit000 (=rdt) if neuler=0
+   LOGICAL , PUBLIC ::   l_1st_euler    !: Euler 1st time-step flag (=T if ln_restart=F or ln_1st_euler=T)
+   REAL(wp), PUBLIC ::   r2dt, r1_2dt   !: = 2*rdt and 1/(2*rdt) except if l_1st_euler=T)
 
    !!----------------------------------------------------------------------

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domain.F90

@@ -288 +288 @@
       INTEGER  ::   ios   ! Local integer
       !
-      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
-         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
-         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
-         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
+      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                   &
+         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl   ,     &
+         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate   ,     &
+         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, ln_1st_euler,     &
          &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios
       NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask
@@ -323 +323 @@
          WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir )
          WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart
-         WRITE(numout,*) '      start with forward time step    nn_euler        = ', nn_euler
+         WRITE(numout,*) '      start with forward time step    ln_1st_euler    = ', ln_1st_euler
          WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl
          WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000
@@ -361 +361 @@
       nstocklist = nn_stocklist
       nwrite = nn_write
-      neuler = nn_euler
-      IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN
+      IF( ln_rstart ) THEN       ! restart : set 1st time-step scheme (LF or forward) 
+         l_1st_euler = ln_1st_euler
+      ELSE                       ! start from rest : always an Euler scheme for the 1st time-step
          IF(lwp) WRITE(numout,*)  
          IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
-         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : nn_euler is forced to 0 '   
-         neuler = 0
+         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used '   
+         l_1st_euler = .TRUE.
       ENDIF
       !                             ! control of output frequency
@@ -374 +375 @@
          nstock = nitend
       ENDIF
-      IF ( nwrite == 0 ) THEN
+      IF( nwrite == 0 ) THEN
          WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
          CALL ctl_warn( ctmp1 )

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90

@@ -56 +56 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors
+   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_a   ! baroclinic scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
@@ -76 +76 @@
       IF( ln_vvl_zstar )   dom_vvl_alloc = 0
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
-            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
+         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   &
+            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
             &      STAT = dom_vvl_alloc        )
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
@@ -103 +103 @@
       !!               - interpolate scale factors
       !!
-      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
+      !! ** Action  : - e3t_(n/b) and te3t_(n/b)
       !!              - Regrid: e3(u/v)_n
       !!                        e3(u/v)_b       
@@ -129 +129 @@
       IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
       !
-      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
+      !                    ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf
       CALL dom_vvl_rst( nit000, 'READ' )
       e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
@@ -280 +280 @@
       !!
       !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
-      !!               - tilde_e3t_a: after increment of vertical scale factor 
+      !!               - te3t_a: after increment of vertical scale factor 
       !!                              in z_tilde case
       !!               - e3(t/u/v)_a
@@ -353 +353 @@
          ! II - after z_tilde increments of vertical scale factors
          ! =======================================================
-         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
+         te3t_a(:,:,:) = 0._wp  ! te3t_a used to store tendency terms
 
          ! 1 - High frequency divergence term
@@ -359 +359 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+               te3t_a(:,:,jk) = te3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          ENDIF
@@ -371 +371 @@
          IF( ln_vvl_ztilde ) THEN
             DO jk = 1, jpk
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
+               te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk)
             END DO
          ENDIF
@@ -383 +383 @@
                DO ji = 1, fs_jpim1   ! vector opt.
                   un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
+                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj  ,jk) )
                   vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
-                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
+                     &            * ( te3t_b(ji,jj,jk) - te3t_b(ji  ,jj+1,jk) )
                   zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
                   zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
@@ -400 +400 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
+                  te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
                      &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
                      &                                            ) * r1_e1e2t(ji,jj)
@@ -414 +414 @@
          ! Leapfrog time stepping
          ! ~~~~~~~~~~~~~~~~~~~~~~
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            z2dt =  rdt
-         ELSE
-            z2dt = 2.0_wp * rdt
-         ENDIF
-         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
-         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
+         CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp )
+         te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:)
 
          ! Maximum deformation control
@@ -426 +421 @@
          ze3t(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
+            ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
          END DO
          z_tmax = MAXVAL( ze3t(:,:,:) )
@@ -446 +441 @@
             ENDIF
             IF (lwp) THEN
-               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
+               WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
                WRITE(numout, *) 'at i, j, k=', ijk_max
-               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
+               WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
                WRITE(numout, *) 'at i, j, k=', ijk_min            
-               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
+               CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
             ENDIF
          ENDIF
          ! - ML - end test
          ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
-         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
-         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
+         te3t_a(:,:,:) = MIN( te3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
+         te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
 
          !
@@ -462 +457 @@
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
+            dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk)
          END DO
          ! III - Barotropic repartition of the sea surface height over the baroclinic profile
@@ -470 +465 @@
          !        i.e. locally and not spread over the water column.
          !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
-         zht(:,:) = 0.
+         zht(:,:) = 0._wp
          DO jk = 1, jpkm1
-            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            zht(:,:)  = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
          z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -484 +479 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -494 +489 @@
             z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
             IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
-            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
+            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax
          END IF
          !
@@ -573 +568 @@
       !!               - recompute depths and water height fields
       !!
-      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
+      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step
       !!               - Recompute:
       !!                    e3(u/v)_b       
@@ -587 +582 @@
       INTEGER, INTENT( in ) ::   kt   ! time step
       !
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   zcoef        ! local scalar
+      INTEGER  ::   ji, jj, jk    ! dummy loop indices
+      REAL(wp) ::   zcoef, ze3f   ! local scalar
       !!----------------------------------------------------------------------
       !
@@ -605 +600 @@
       ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         IF( neuler == 0 .AND. kt == nit000 ) THEN
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
+         IF( l_1st_euler ) THEN
+            te3t_n(:,:,:) = te3t_a(:,:,:)
          ELSE
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
-            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     ze3f = te3t_n(ji,jj,jk)   &
+                        & + atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) )
+                     te3t_b(ji,jj,jk) = ze3f
+                     te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
          ENDIF
-         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
       gdept_b(:,:,:) = gdept_n(:,:,:)
@@ -806 +808 @@
             CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
             !
-            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
-            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
+            id1 = iom_varid( numror, 'e3t_b'      , ldstop = .FALSE. )
+            id2 = iom_varid( numror, 'e3t_n'      , ldstop = .FALSE. )
             id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
-            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. )
             !                             ! --------- !
             !                             ! all cases !
@@ -823 +825 @@
                   e3t_b(:,:,:) = e3t_0(:,:,:)
                END WHERE
-               IF( neuler == 0 ) THEN
+               IF( l_1st_euler ) THEN
                   e3t_b(:,:,:) = e3t_n(:,:,:)
                ENDIF
@@ -829 +831 @@
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
                e3t_n(:,:,:) = e3t_b(:,:,:)
-               neuler = 0
+               l_1st_euler = .TRUE.
             ELSE IF( id2 > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
                e3t_b(:,:,:) = e3t_n(:,:,:)
-               neuler = 0
+               l_1st_euler = .TRUE.
             ELSE
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
-               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
                DO jk = 1, jpk
                   e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
@@ -850 +852 @@
                END DO
                e3t_b(:,:,:) = e3t_n(:,:,:)
-               neuler = 0
+               l_1st_euler = .TRUE.
             ENDIF
             !                             ! ----------- !
@@ -862 +864 @@
                !                          ! ----------------------- !
                IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
-                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
-                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
+                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios )
+                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios )
                ELSE                            ! one at least array is missing
-                  tilde_e3t_b(:,:,:) = 0.0_wp
-                  tilde_e3t_n(:,:,:) = 0.0_wp
+                  te3t_b(:,:,:) = 0.0_wp
+                  te3t_n(:,:,:) = 0.0_wp
                ENDIF
                !                          ! ------------ !
@@ -942 +944 @@
 
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
-               tilde_e3t_b(:,:,:) = 0._wp
-               tilde_e3t_n(:,:,:) = 0._wp
+               te3t_b(:,:,:) = 0._wp
+               te3t_n(:,:,:) = 0._wp
                IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
             END IF
@@ -960 +962 @@
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
             !                                        ! ----------------------- !
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
-            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios)
+            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios)
          END IF
          !                                           ! -------------!    

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/iscplrst.F90

@@ -89 +89 @@
       END IF
       !
-      neuler = 0              ! next step is an euler time step
+      l_1st_euler = .TRUE.    ! next step is an euler time step
       !
       !                       ! set _b and _n variables equal

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/istate.F90

@@ -92 +92 @@
          !                                    ! ---------------
          numror = 0                           ! define numror = 0 -> no restart file to read
-         neuler = 0                           ! Set time-step indicator at nit000 (euler forward)
+         l_1st_euler = .TRUE.                 ! Set a Euler forward 1sr time-step
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !                                    ! Initialization of ocean to zero

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/restart.F90

@@ -8 +8 @@
    !!            2.0  !  2006-07  (S. Masson)  use IOM for restart
    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
-   !!            - -  !  2010-10  (C. Ethe, G. Madec) TRC-TRA merge (T-S in 4D)
-   !!            3.7  !  2014-01  (G. Madec) suppression of curl and hdiv from the restart
-   !!             -   !  2014-12  (G. Madec) remove KPP scheme
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   rst_opn    : open the ocean restart file
-   !!   rst_write  : write the ocean restart file
-   !!   rst_read   : read the ocean restart file
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_ice         ! only lk_si3 
-   USE phycst          ! physical constants
-   USE eosbn2          ! equation of state            (eos bn2 routine)
-   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
+   !!            - -  !  2010-10  (C. Ethe, G. Madec)  TRC-TRA merge (T-S in 4D)
+   !!            3.7  !  2014-01  (G. Madec)  suppression of curl and hdiv from the restart
+   !!             -   !  2014-12  (G. Madec)  remove KPP scheme
+   !!            4.0  !  2018-06  (G. Madec)  introduce l_1st_euler
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   rst_opn       : open the ocean restart file in write mode
+   !!   rst_write     : write the ocean restart file
+   !!   rst_read_open : open the ocean restart file in read mode
+   !!   rst_read      : read the ocean restart file
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_ice        ! only lk_si3 
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state            (eos bn2 routine)
+   USE trdmxl_oce     ! ocean active mixed layer tracers trends variables
+   USE diurnal_bulk   ! 
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE diurnal_bulk
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O module
 
    IMPLICIT NONE
@@ -34 +36 @@
    PUBLIC   rst_opn         ! routine called by step module
    PUBLIC   rst_write       ! routine called by step module
+   PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init
    PUBLIC   rst_read        ! routine called by istate module
-   PUBLIC   rst_read_open   ! routine called in rst_read and (possibly) in dom_vvl_init
 
    !! * Substitutions
@@ -144 +146 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step
       !!----------------------------------------------------------------------
-                     IF(lwxios) CALL iom_swap(      cwxios_context          )
-                     CALL iom_rstput( kt, nitrst, numrow, 'rdt'    , rdt       , ldxios = lwxios)   ! dynamics time step
-
-      IF ( .NOT. ln_diurnal_only ) THEN
-                     CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub, ldxios = lwxios        )     ! before fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb, ldxios = lwxios        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb, ldxios = lwxios      )
-                     !
-                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un, ldxios = lwxios        )     ! now fields
-                     CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn, ldxios = lwxios        )
-                     CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios )
-                     CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn, ldxios = lwxios      )
-                     CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop, ldxios = lwxios      )
-                  ! extra variable needed for the ice sheet coupling
-                  IF ( ln_iscpl ) THEN 
-                     CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask, ldxios = lwxios ) ! need to extrapolate T/S
-                     CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask, ldxios = lwxios ) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask, ldxios = lwxios ) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask, ldxios = lwxios) ! need to correct barotropic velocity
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )   ! need to compute temperature correction
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3u_n', e3u_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'e3v_n', e3v_n(:,:,:), ldxios = lwxios )   ! need to compute bt conservation
-                     CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n(:,:,:), ldxios = lwxios ) ! need to compute extrapolation if vvl
-                  END IF
-      ENDIF
-      
-      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios )  
-      IF(lwxios) CALL iom_swap(      cxios_context          )
+      IF( lwxios )   CALL iom_swap( cwxios_context )
+         
+      CALL    iom_rstput( kt, nitrst, numrow, 'rdt'    , rn_rdt           , ldxios = lwxios )   ! dynamics time step
+      !
+      IF( .NOT. ln_diurnal_only ) THEN
+         CALL iom_rstput( kt, nitrst, numrow, 'ub'     , ub               , ldxios = lwxios )   ! before fields
+         CALL iom_rstput( kt, nitrst, numrow, 'vb'     , vb               , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'tb'     , tsb(:,:,:,jp_tem), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sb'     , tsb(:,:,:,jp_sal), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb             , ldxios = lwxios )
+         !
+         CALL iom_rstput( kt, nitrst, numrow, 'un'     , un               , ldxios = lwxios )     ! now fields
+         CALL iom_rstput( kt, nitrst, numrow, 'vn'     , vn               , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'tn'     , tsn(:,:,:,jp_tem), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sn'     , tsn(:,:,:,jp_sal), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn             , ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop             , ldxios = lwxios )
+         !
+         IF( ln_iscpl ) THEN          ! extra variable needed for the ice sheet coupling
+            CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask  , ldxios = lwxios )    ! need to extrapolate T/S
+            CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask  , ldxios = lwxios )    ! need to correct barotropic velocity
+            CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask  , ldxios = lwxios )    ! need to correct barotropic velocity
+            CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask , ldxios = lwxios )    ! need to correct barotropic velocity
+            CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n  , ldxios = lwxios )    ! need to compute temperature correction
+            CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n  , ldxios = lwxios )    ! need to compute bt conservation
+            CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n  , ldxios = lwxios )    ! need to compute bt conservation
+            CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n, ldxios = lwxios )    ! need to compute extrapolation if vvl
+         ENDIF
+      ENDIF
+      !
+      IF( ln_diurnal )   CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios )  
+      IF( lwxios     )   CALL iom_swap( cxios_context )
       IF( kt == nitrst ) THEN
-         IF(.NOT.lwxios) THEN
-            CALL iom_close( numrow )     ! close the restart file (only at last time step)
-         ELSE
-            CALL iom_context_finalize(      cwxios_context          )
+         IF( lwxios ) THEN   ;   CALL iom_context_finalize( cwxios_context )
+         ELSE                ;   CALL iom_close( numrow )     ! close the restart file (only at last time step)
          ENDIF
 !!gm         IF( .NOT. lk_trdmld )   lrst_oce = .FALSE.
 !!gm  not sure what to do here   ===>>>  ask to Sebastian
          lrst_oce = .FALSE.
-            IF( ln_rst_list ) THEN
-               nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1))
-               nitrst = nstocklist( nrst_lst )
-            ENDIF
+         IF( ln_rst_list ) THEN
+            nrst_lst = MIN(nrst_lst + 1, SIZE(nstocklist,1))
+            nitrst = nstocklist( nrst_lst )
+         ENDIF
       ENDIF
       !
@@ -202 +203 @@
       !!                the file has already been opened
       !!----------------------------------------------------------------------
-      INTEGER        ::   jlibalt = jprstlib
-      LOGICAL        ::   llok
-      CHARACTER(lc)  ::   clpath   ! full path to ocean output restart file
+      INTEGER       ::   jlibalt = jprstlib
+      LOGICAL       ::   llok
+      CHARACTER(lc) ::   clpath   ! full path to ocean output restart file
       !!----------------------------------------------------------------------
       !
@@ -238 +239 @@
          ENDIF 
       ENDIF
-
+      !
    END SUBROUTINE rst_read_open
 
@@ -254 +255 @@
       REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
       !!----------------------------------------------------------------------
-
+      !
       CALL rst_read_open           ! open restart for reading (if not already opened)
 
@@ -260 +261 @@
       IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN
          CALL iom_get( numror, 'rdt', zrdt, ldxios = lrxios )
-         IF( zrdt /= rdt )   neuler = 0
+         IF( zrdt /= rn_rdt ) THEN
+            IF(lwp) WRITE( numout,*)
+            IF(lwp) WRITE( numout,*) 'rst_read:  rdt not equal to the read one'
+            IF(lwp) WRITE( numout,*)
+            IF(lwp) WRITE( numout,*) '      ==>>>   forced euler first time-step'
+            l_1st_euler =  .TRUE.
+         ENDIF
       ENDIF
 
@@ -266 +273 @@
       IF( ln_diurnal ) CALL iom_get( numror, jpdom_autoglo, 'Dsst' , x_dsst, ldxios = lrxios ) 
       IF ( ln_diurnal_only ) THEN 
-         IF(lwp) WRITE( numout, * ) &
-         &   "rst_read:- ln_diurnal_only set, setting rhop=rau0" 
+         IF(lwp) WRITE( numout,*) 'rst_read: ln_diurnal_only set, setting rhop=rau0'
          rhop = rau0
          CALL iom_get( numror, jpdom_autoglo, 'tn'     , w3d, ldxios = lrxios ) 
@@ -274 +280 @@
       ENDIF  
       
-      IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0 ) THEN
+      IF( iom_varid( numror, 'ub', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'ub'     , ub, ldxios = lrxios                )   ! before fields
          CALL iom_get( numror, jpdom_autoglo, 'vb'     , vb, ldxios = lrxios                )
@@ -281 +287 @@
          CALL iom_get( numror, jpdom_autoglo, 'sshb'   , sshb, ldxios = lrxios              )
       ELSE
-         neuler = 0
+         l_1st_euler =  .TRUE.      ! before field not found, forced euler 1st time-step
       ENDIF
       !
@@ -295 +301 @@
       ENDIF
       !
-      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
-         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
-         ub   (:,:,:)   = un   (:,:,:)
-         vb   (:,:,:)   = vn   (:,:,:)
-         sshb (:,:)     = sshn (:,:)
-         !
-         IF( .NOT.ln_linssh ) THEN
-            DO jk = 1, jpk
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-            END DO
-         ENDIF
-         !
+      IF( l_1st_euler ) THEN              ! Euler restart
+         tsb (:,:,:,:) = tsn (:,:,:,:)          ! all before fields set to now values
+         ub  (:,:,:)   = un  (:,:,:)
+         vb  (:,:,:)   = vn  (:,:,:)
+         sshb(:,:)     = sshn(:,:)
+         IF( .NOT.ln_linssh )   e3t_b(:,:,:) = e3t_n(:,:,:)
       ENDIF
       !

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynnxt.F90

@@ -64 +64 @@
 CONTAINS
 
-   SUBROUTINE dyn_nxt ( kt )
+   SUBROUTINE dyn_nxt( kt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_nxt  ***
@@ -92 +92 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ikt          ! local integers
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef   ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, zvf          !   -      -
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
@@ -132 +132 @@
             ! so that asselin contribution is removed at the same time 
             DO jk = 1, jpkm1
-               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
-               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
+               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) ) * umask(:,:,jk)
+               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) ) * vmask(:,:,jk)
             END DO  
          ENDIF
@@ -152 +152 @@
 !!$   Do we need a call to bdy_vol here??
       !
-      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
-         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
-         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
-         !
-         !                                  ! Kinetic energy and Conversion
-         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
-         !
-         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
-            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
-            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
-            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
+      IF( l_trddyn ) THEN             !* prepare the atf trend computation + some diagnostics
+         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )    ! Kinetic energy and Conversion
+         IF( ln_dyn_trd ) THEN                                       ! 3D output: total momentum trends
+            IF( ln_dynadv_vec ) THEN 
+               zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_2dt
+               zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_2dt
+            ELSE
+               zua(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt
+               zva(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt
+            ENDIF
+            CALL iom_put( "utrd_tot", zua )                          ! total momentum trends, except the asselin filter
             CALL iom_put( "vtrd_tot", zva )
          ENDIF
-         !
-         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
-         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
-         !                                  !  computation of the asselin filter trends)
+         zua(:,:,:) = un(:,:,:)                    ! save the now velocity before the asselin filter
+         zva(:,:,:) = vn(:,:,:)                    ! (caution: the Asselin filter trends computation will be shifted by 1 timestep)
       ENDIF
 
       ! Time filter and swap of dynamics arrays
-      ! ------------------------------------------
-      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
+      ! ---------------------------------------
+      IF( l_1st_euler ) THEN        !==  Euler at 1st time-step  ==!   (swap only)
          DO jk = 1, jpkm1
             un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
             vn(:,:,jk) = va(:,:,jk)
          END DO
-         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
-!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a  
+         IF( .NOT.ln_linssh ) THEN                          ! e3._n <-- e3._a
             DO jk = 1, jpkm1
-!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
-!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
-               !
                e3t_n(:,:,jk) = e3t_a(:,:,jk)
                e3u_n(:,:,jk) = e3u_a(:,:,jk)
                e3v_n(:,:,jk) = e3v_a(:,:,jk)
             END DO
-!!gm BUG end
-         ENDIF
-                                                            ! 
-         
-      ELSE                                             !* Leap-Frog : Asselin filter and swap
+         ENDIF
+         !
+      ELSE                          !==  Leap-Frog  ==!   (Asselin filter and swap)
+         !
          !                                ! =============!
          IF( ln_linssh ) THEN             ! Fixed volume !
@@ -213 +205 @@
          ELSE                             ! Variable volume !
             !                             ! ================!
-            ! Before scale factor at t-points
-            ! (used as a now filtered scale factor until the swap)
-            ! ----------------------------------------------------
+            ! Before scale factor at t-points   (used as a now filtered scale factor until the swap)
             DO jk = 1, jpkm1
                e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
             END DO
             ! Add volume filter correction: compatibility with tracer advection scheme
-            ! => time filter + conservation correction (only at the first level)
+            !    => time filter + conservation correction (only at the first level)
             zcoef = atfp * rdt * r1_rau0
 
@@ -232 +222 @@
                            IF( jk <=  nk_rnf(ji,jj)  ) THEN
                                e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
-                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
+                                  &              * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
                            ENDIF
-                        ENDDO
-                     ENDDO
-                  ENDDO
+                        END DO
+                     END DO
+                  END DO
                ELSE
                   e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
                ENDIF
-            END IF
+            ENDIF
 
             IF ( ln_isf ) THEN   ! if ice shelf melting
@@ -253 +243 @@
                   END DO
                END DO
-            END IF
+            ENDIF
             !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
@@ -322 +312 @@
          ENDIF
          !
-      ENDIF ! neuler =/0
+      ENDIF    ! end Leap-Frog time stepping
       !
       ! Set "now" and "before" barotropic velocities for next time step:
-      ! JC: Would be more clever to swap variables than to make a full vertical
-      ! integration
+      ! JC: Would be more clever to swap variables than to make a full vertical integration
       !
       !
@@ -360 +349 @@
       ENDIF
       IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
-         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
-         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * r1_2dt
+         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * r1_2dt
          CALL trd_dyn( zua, zva, jpdyn_atf, kt )
       ENDIF

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg.F90

@@ -74 +74 @@
       !
       INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
+      REAL(wp) ::   zg_2, zintp, zgrau0r, zld   ! local scalars
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90

@@ -1 +1 @@
 MODULE dynspg_ts
-
-   !! Includes ROMS wd scheme with diagnostic outputs ; un and ua updates are commented out ! 
-
    !!======================================================================
    !!                   ***  MODULE  dynspg_ts  ***
@@ -35 +32 @@
    USE sbcisf          ! ice shelf variable (fwfisf)
    USE sbcapr          ! surface boundary condition: atmospheric pressure
-   USE dynadv    , ONLY: ln_dynadv_vec
+   USE dynadv   , ONLY : ln_dynadv_vec
    USE dynvor          ! vortivity scheme indicators
    USE phycst          ! physical constants
@@ -85 +82 @@
    REAL(wp) ::   r1_2  = 0.5_wp           !
 
+   REAL(wp) ::   r1_2dt_b, r2dt_bf               ! local scalars
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -151 +150 @@
       INTEGER  ::   ikbu, iktu, noffset   ! local integers
       INTEGER  ::   ikbv, iktv            !   -      -
-      REAL(wp) ::   r1_2dt_b, z2dt_bf               ! local scalars
       REAL(wp) ::   zx1, zx2, zu_spg, zhura, z1_hu  !   -      -
       REAL(wp) ::   zy1, zy2, zv_spg, zhvra, z1_hv  !   -      -
@@ -182 +180 @@
 !     zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
       !                                         ! reciprocal of baroclinic time step 
-      IF( kt == nit000 .AND. neuler == 0 ) THEN   ;   z2dt_bf =          rdt
-      ELSE                                        ;   z2dt_bf = 2.0_wp * rdt
-      ENDIF
-      r1_2dt_b = 1.0_wp / z2dt_bf 
+      IF( l_1st_euler ) THEN   ;   r2dt_bf =          rdt
+      ELSE                     ;   r2dt_bf = 2.0_wp * rdt
+      ENDIF
+      r1_2dt_b = 1.0_wp / r2dt_bf 
       !
       ll_init     = ln_bt_av                    ! if no time averaging, then no specific restart 
@@ -194 +192 @@
       ENDIF
       !
-      IF( kt == nit000 ) THEN                   !* initialisation
+      IF( kt == nit000 ) THEN                   !* initialisation 1st time-step
          !
          IF(lwp) WRITE(numout,*)
@@ -201 +199 @@
          IF(lwp) WRITE(numout,*)
          !
-         IF( neuler == 0 )   ll_init=.TRUE.
-         !
-         IF( ln_bt_fw .OR. neuler == 0 ) THEN
+         IF( l_1st_euler )   ll_init = .TRUE.
+         !
+         IF( ln_bt_fw .OR. l_1st_euler ) THEN
             ll_fw_start =.TRUE.
             noffset     = 0
@@ -212 +210 @@
          ! Set averaging weights and cycle length:
          CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
+         !
+      ELSEIF( kt == nit000 + 1 ) THEN           !* initialisation 2nd time-step
+         !
+         IF( .NOT.ln_bt_fw .AND. l_1st_euler ) THEN
+            ll_fw_start = .FALSE.
+            CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
+         ENDIF
          !
       ENDIF
@@ -340 +345 @@
          END SELECT
       ENDIF
-      !
-      ! If forward start at previous time step, and centered integration, 
-      ! then update averaging weights:
-      IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN
-         ll_fw_start=.FALSE.
-         CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 )
-      ENDIF
-                          
+      !                          
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
@@ -1203 +1201 @@
          zwx(:,:) = un_adv(:,:)
          zwy(:,:) = vn_adv(:,:)
-         IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN
+         IF( .NOT.l_1st_euler ) THEN
             un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )
             vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )
@@ -1305 +1303 @@
       !! ** Purpose : Set time-splitting weights for temporal averaging (or not)
       !!----------------------------------------------------------------------
-      LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
-      LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
-      INTEGER, INTENT(inout) :: jpit      ! cycle length    
-      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
-                                                         zwgt2    ! Secondary weights
-      
+      LOGICAL                       , INTENT(in   ) ::   ll_av          ! temporal averaging=.true.
+      LOGICAL                       , INTENT(in   ) ::   ll_fw          ! forward time splitting =.true.
+      INTEGER                       , INTENT(inout) ::   jpit           ! cycle length    
+      REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, zwgt2   ! Primary & Secondary weights
+      !
       INTEGER ::  jic, jn, ji                      ! temporary integers
       REAL(wp) :: za1, za2
       !!----------------------------------------------------------------------
-
+      !
       zwgt1(:) = 0._wp
       zwgt2(:) = 0._wp
-
+      !
       ! Set time index when averaged value is requested
-      IF (ll_fw) THEN 
-         jic = nn_baro
-      ELSE
-         jic = 2 * nn_baro
-      ENDIF
-
-      ! Set primary weights:
-      IF (ll_av) THEN
-           ! Define simple boxcar window for primary weights 
-           ! (width = nn_baro, centered around jic)     
+      IF ( ll_fw ) THEN   ;   jic =     nn_baro
+      ELSE                ;   jic = 2 * nn_baro
+      ENDIF
+
+      !                 !==  Set primary weights  ==!
+      !
+      IF (ll_av) THEN         !* Define simple boxcar window for primary weights 
+         !                    !  (width = nn_baro, centered around jic)     
          SELECT CASE ( nn_bt_flt )
-              CASE( 0 )  ! No averaging
-                 zwgt1(jic) = 1._wp
-                 jpit = jic
-
-              CASE( 1 )  ! Boxcar, width = nn_baro
-                 DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
-                    IF (za1 < 0.5_wp) THEN
-                      zwgt1(jn) = 1._wp
-                      jpit = jn
-                    ENDIF
-                 ENDDO
-
-              CASE( 2 )  ! Boxcar, width = 2 * nn_baro
-                 DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
-                    IF (za1 < 1._wp) THEN
-                      zwgt1(jn) = 1._wp
-                      jpit = jn
-                    ENDIF
-                 ENDDO
-              CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
+         CASE( 0 )  ! No averaging
+            zwgt1(jic) = 1._wp
+            jpit = jic
+            !
+         CASE( 1 )  ! Boxcar, width = nn_baro
+            DO jn = 1, 3*nn_baro
+               za1 = ABS(float(jn-jic))/float(nn_baro) 
+               IF ( za1 < 0.5_wp ) THEN
+                  zwgt1(jn) = 1._wp
+                  jpit = jn
+               ENDIF
+            END DO
+            !
+         CASE( 2 )  ! Boxcar, width = 2 * nn_baro
+            DO jn = 1, 3*nn_baro
+               za1 = ABS(float(jn-jic))/float(nn_baro) 
+                  IF ( za1 < 1._wp ) THEN
+                     zwgt1(jn) = 1._wp
+                     jpit = jn
+                  ENDIF
+               END DO
+         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
          END SELECT
-
-      ELSE ! No time averaging
+         !
+      ELSE                    !* No time averaging
          zwgt1(jic) = 1._wp
          jpit = jic
       ENDIF
     
-      ! Set secondary weights
+      !                 !==  Set secondary weights  ==!
+      !
       DO jn = 1, jpit
-        DO ji = jn, jpit
-             zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
-        END DO
+         DO ji = jn, jpit
+            zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
+         END DO
       END DO
 
-      ! Normalize weigths:
-      za1 = 1._wp / SUM(zwgt1(1:jpit))
-      za2 = 1._wp / SUM(zwgt2(1:jpit))
+      !                 !==  Normalize weights  ==!
+      !
+      za1 = 1._wp / SUM( zwgt1(1:jpit) )
+      za2 = 1._wp / SUM( zwgt2(1:jpit) )
       DO jn = 1, jpit
-        zwgt1(jn) = zwgt1(jn) * za1
-        zwgt2(jn) = zwgt2(jn) * za2
+         zwgt1(jn) = zwgt1(jn) * za1
+         zwgt2(jn) = zwgt2(jn) * za2
       END DO
       !
@@ -1539 +1536 @@
       !
       !                             ! read restart when needed
+!!gm what's happen when starting with an euler time-step BUT not from rest ?
+!!   this case correspond to a restart with only now time-step available...
       CALL ts_rst( nit000, 'READ' )
       !
@@ -1548 +1547 @@
          CALL iom_set_rstw_var_active('vn_bf')
          !
-         IF (.NOT.ln_bt_av) THEN
+         IF ( .NOT.ln_bt_av ) THEN
             CALL iom_set_rstw_var_active('sshbb_e')
             CALL iom_set_rstw_var_active('ubb_e')

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynzdf.F90

@@ -11 +11 @@
    !!----------------------------------------------------------------------
    !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
+   !!       zdf_trd   : diagnose the zdf velocity trends and the KE dissipation trend 
+!!gm                        ==>>> zdf_trd currently not used
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -26 +28 @@
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
+   USE iom             ! IOM library
    USE prtctl         ! Print control
    USE timing         ! Timing
@@ -67 +70 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   ji, jj, jk         ! dummy loop indices
-      INTEGER  ::   iku, ikv           ! local integers
-      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
-      REAL(wp) ::   zzws, ze3va        !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace 
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
+      INTEGER  ::   iku, ikv              ! local integers
+      REAL(wp) ::   zzwi, ze3ua, z2dt_2   ! local scalars
+      REAL(wp) ::   zzws, ze3va           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwd, zws   ! 3D workspace 
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv    !  -      -
       !!---------------------------------------------------------------------
       !
@@ -86 +89 @@
          ENDIF
       ENDIF
-      !                             !* set time step
-      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
-      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
-      ENDIF
+      !
+      z2dt_2 = r2dt * 0.5_wp        !* =rdt except in 1st Euler time step where it is equal to rdt/2
+      !
       !
       !                             !* explicit top/bottom drag case
@@ -111 +113 @@
       ELSE                                      ! applied on thickness weighted velocity
          DO jk = 1, jpkm1
-            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
-               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
-            va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
-               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
+            ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + r2dt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)
+            va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + r2dt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -133 +133 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -144 +144 @@
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
-                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
+                  va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
                END DO
             END DO
@@ -153 +153 @@
       !              !==  Vertical diffusion on u  ==!
       !
-      !                    !* Matrix construction
-      zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+      SELECT CASE( nldf_dyn )    !* Matrix construction
+      !
+      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
-                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zzwi = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) + akzu(ji,jj,jk  ) )   &
+                     &            / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - r2dt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) )   &
+                     &            / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -171 +170 @@
             END DO
          END DO
-      CASE DEFAULT               ! iso-level lateral mixing
+      CASE DEFAULT                  ! iso-level lateral mixing
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1 
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zzwi = - z2dt_2 * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi
                   zws(ji,jj,jk) = zzws
@@ -186 +185 @@
       END SELECT
       !
-      DO jj = 2, jpjm1     !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1           !* Surface boundary conditions
+         DO ji = fs_2, fs_jpim1     ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -204 +203 @@
                iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
                ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
-               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
+               zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
             END DO
          END DO
@@ -213 +212 @@
                   iku = miku(ji,jj)       ! ocean top level at u- and v-points 
                   ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
-                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
                END DO
             END DO
@@ -245 +244 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
-            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
+            ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rau0 ) * umask(ji,jj,1)
          END DO
       END DO
@@ -272 +270 @@
       !              !==  Vertical diffusion on v  ==!
       !
-      !                       !* Matrix construction
-      zdt = r2dt * 0.5
-      SELECT CASE( nldf_dyn )
-      CASE( np_lap_i )           ! rotated lateral mixing: add its vertical mixing (akzu)
+      !                      
+      SELECT CASE( nldf_dyn )    !* Matrix construction
+      CASE( np_lap_i )              ! rotated lateral mixing: add its vertical mixing (akzu)
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
-                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zzwi = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) + akzv(ji,jj,jk  ) )   &
+                     &            / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - r2dt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) )   &
+                     &            / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -290 +287 @@
             END DO
          END DO
-      CASE DEFAULT               ! iso-level lateral mixing
+      CASE DEFAULT                  ! iso-level lateral mixing
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1   
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-                  zzwi = - zdt * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
-                  zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zzwi = - z2dt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
                   zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
                   zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
@@ -305 +302 @@
       END SELECT
       !
-      DO jj = 2, jpjm1        !* Surface boundary conditions
-         DO ji = fs_2, fs_jpim1   ! vector opt.
+      DO jj = 2, jpjm1           !* Surface boundary conditions
+         DO ji = fs_2, fs_jpim1     ! vector opt.
             zwi(ji,jj,1) = 0._wp
             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
@@ -322 +319 @@
                ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
                ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
-               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
+               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
             END DO
          END DO
@@ -330 +327 @@
                   ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
                   ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
-                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
                END DO
             END DO
@@ -362 +359 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.          
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
-            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
+            va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rau0 ) * vmask(ji,jj,1)
          END DO
       END DO
@@ -387 +383 @@
       END DO
       !
-      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
+      IF( l_trddyn ) THEN                        ! save the vertical diffusive trends for further diagnostics
+         IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
+            ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_2dt - ztrdu(:,:,:)
+            ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_2dt - ztrdv(:,:,:)
+         ELSE                                      ! applied on thickness weighted velocity
+            ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ztrdu(:,:,:)
+            ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ztrdv(:,:,:)
+         ENDIF
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
          DEALLOCATE( ztrdu, ztrdv ) 
@@ -401 +402 @@
    END SUBROUTINE dyn_zdf
 
+!!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing
+
+   SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_trd  ***
+      !!
+      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
+      !!              together with the Leap-Frog time stepping using an 
+      !!              implicit scheme.
+      !!
+      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
+      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
+      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
+      !!               - update the after velocity with the implicit vertical mixing.
+      !!      This requires to solver the following system: 
+      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
+      !!      with the following surface/top/bottom boundary condition:
+      !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
+      !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
+      !!
+      !! ** Action :   (ua,va)   after velocity 
+      !!---------------------------------------------------------------------
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   ptrdu, ptrdv   ! 3D work arrays use for zdf trends diag
+      INTEGER , INTENT(in   )                         ::   kt             ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      REAL(wp) ::   zzz              ! local scalar
+      REAL(wp) ::   zavmu, zavmum1   !   -      -
+      REAL(wp) ::   zavmv, zavmvm1   !   -      -
+      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   z2d    !  -      -
+      !!---------------------------------------------------------------------
+      !
+      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )   ! apply lateral boundary condition on (ua,va)
+      !
+      !
+      !                 !==  momentum trend due to vertical diffusion  ==!
+      !
+      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
+         ptrdu(:,:,:) = (              ua(:,:,:) -              ub(:,:,:) )                * r1_2dt - ptrdu(:,:,:)
+         ptrdv(:,:,:) = (              va(:,:,:) -              vb(:,:,:) )                * r1_2dt - ptrdv(:,:,:)
+      ELSE                                      ! applied on thickness weighted velocity
+         ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_2dt - ptrdu(:,:,:)
+         ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_2dt - ptrdv(:,:,:)
+      ENDIF
+      CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt )
+      !
+      !
+      !                 !==  KE dissipation trend due to vertical diffusion  ==!
+      !
+      IF( iom_use( 'dispkevfo' ) ) THEN   ! ocean kinetic energy dissipation per unit area
+         !                                ! due to v friction (v=vertical) 
+         !                                ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points)
+         !                                ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of 
+         !                                ! now by before shears, i.e. the source term of TKE (local positivity is not ensured).
+         !                                ! Note also that now e3 scale factors are used as after one are not computed !
+         !
+         CALL wrk_alloc(jpi,jpj,   z2d )
+         z2d(:,:) = 0._wp
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  zavmu   = 0.5 * ( avm(ji+1,jj,jk) + avm(ji  ,jj,jk) )
+                  zavmum1 = 0.5 * ( avm(ji  ,jj,jk) + avm(ji-1,jj,jk) )
+                  zavmv   = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj  ,jk) )
+                  zavmvm1 = 0.5 * ( avm(ji,jj  ,jk) + avm(ji,jj-1,jk) )
+               
+                  z2d(ji,jj) = z2d(ji,jj)  +  (                                                                                  &
+                     &   zavmu   * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) )**2 / e3uw_n(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
+                     & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
+                     & + zavmv   * ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) )**2 / e3vw_n(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
+                     & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
+                     &                        )
+!!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step  ===>>> use it ?
+!!                                                                                     No since in zdfshé only kz tke (or gls) is used
+!!
+!!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of
+!!gm     now by before shears, i.e. the source term of TKE (local positivity is not ensured).
+!!       CAUTION: requires to compute e3uw_a and e3vw_a !!!
+!                  z2d(ji,jj) = z2d(ji,jj)  + (                                                                                 &
+!                     &    avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) / e3uw_n(ji  ,jj,jk)                        &
+!                     &                     * ( ua(ji  ,jj,jk-1) - ua(ji  ,jj,jk) ) / e3uw_a(ji  ,jj,jk) * wumask(ji  ,jj,jk)   &
+!                     &  + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk)                        &
+!                     &                       ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk)   &
+!                     &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) / e3vw_n(ji,jj  ,jk)                        &
+!                     &                       ( va(ji,jj  ,jk-1) - va(ji,jj  ,jk) ) / e3vw_a(ji,jj  ,jk) * wvmask(ji,jj  ,jk)   &
+!                     &  + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk)                        &
+!                     &                       ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk)   &
+!                     &                       )
+!!gm end
+               END DO
+            END DO
+         END DO
+         zzz= - 0.5_wp* rau0           ! caution sign minus here
+         z2d(:,:) = zzz * z2d(:,:) 
+         CALL lbc_lnk( z2d,'T', 1. )
+         CALL iom_put( 'dispkevfo', z2d )
+      ENDIF
+      !
+   END SUBROUTINE zdf_trd
+   
+!!gm end not used
+
    !!==============================================================================
 END MODULE dynzdf

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/sshwzv.F90

@@ -68 +68 @@
       INTEGER, INTENT(in) ::   kt   ! time step
       ! 
-      INTEGER  ::   jk            ! dummy loop indice
-      REAL(wp) ::   z2dt, zcoef   ! local scalars
+      INTEGER  ::   jk         ! dummy loop indice
+      REAL(wp) ::   z1_2rau0   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
       !!----------------------------------------------------------------------
@@ -81 +81 @@
       ENDIF
       !
-      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
-      zcoef = 0.5_wp * r1_rau0
+      z1_2rau0 = 0.5_wp * r1_rau0
 
       !                                           !------------------------------!
       !                                           !   After Sea Surface Height   !
       !                                           !------------------------------!
-      IF(ln_wd_il) THEN
-         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
-      ENDIF
+
+      IF(ln_wd_il)   CALL wad_lmt( sshb, z1_2rau0 * (emp_b(:,:) + emp(:,:)), r2dt )
 
       CALL div_hor( kt )                               ! Horizontal divergence
@@ -102 +99 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+      ssha(:,:) = (  sshb(:,:) - r2dt * ( z1_2rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
       !
 #if defined key_agrif
@@ -143 +140 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp) ::   z1_2dt       ! local scalars
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
       !!----------------------------------------------------------------------
@@ -159 +155 @@
       !                                           !     Now Vertical Velocity    !
       !                                           !------------------------------!
-      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
       !
       IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
@@ -180 +174 @@
             ! computation of w
             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
-               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
+               &                         + r1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
          END DO
          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
@@ -188 +182 @@
             ! computation of w
             wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
-               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
+               &                         + r1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -243 +237 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
-      !              !==  Euler time-stepping: no filter, just swap  ==!
-      IF ( neuler == 0 .AND. kt == nit000 ) THEN
-         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
-         !
-      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
-         !                                                  ! before <-- now filtered
+      !              
+      IF ( l_1st_euler ) THEN    !==  Euler time-stepping  ==!   no filter, just swap
+         !
+         sshn(:,:) = ssha(:,:)               ! now    <-- after  (before already = now)
+         !
+      ELSE                       !==  Leap-Frog time-stepping  ==!   Asselin filter + swap
+         !
+         !                                   ! before <-- now filtered
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
-         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
+         IF( .NOT.ln_linssh ) THEN           ! before <-- with forcing removed
             zcoef = atfp * rdt * r1_rau0
             sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
@@ -256 +252 @@
                &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
          ENDIF
-         sshn(:,:) = ssha(:,:)                              ! now <-- after
+         sshn(:,:) = ssha(:,:)               ! now <-- after
       ENDIF
       !

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/wet_dry.F90

@@ -117 +117 @@
 
 
-   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
+   SUBROUTINE wad_lmt( sshb1, sshemp, p2dt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
@@ -129 +129 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1        !!gm DOCTOR names: should start with p !
       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   sshemp
-      REAL(wp)                , INTENT(in   ) ::   z2dt
+      REAL(wp)                , INTENT(in   ) ::   p2dt
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
@@ -220 +220 @@
                   &   + MIN( zflxv1(ji,jj) , 0._wp ) - MAX( zflxv1(ji,  jj-1) , 0._wp) 
                !
-               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
-               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)
+               zdep1 = (zzflxp + zzflxn) * p2dt / ztmp
+               zdep2 = ht_0(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - p2dt * sshemp(ji,jj)
                !
                IF( zdep1 > zdep2 ) THEN
                   wdmask(ji, jj) = 0._wp
-                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
-                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                  zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * p2dt ) / ( zflxp(ji,jj) * p2dt )
+                  !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * p2dt ) / ( zzflxp * p2dt )
                   ! flag if the limiter has been used but stop flagging if the only
                   ! changes have zeroed the coefficient since further iterations will
@@ -270 +270 @@
 
 
-   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
+   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, pdt )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE wad_lmt  ***
@@ -280 +280 @@
       !! ** Action  : - calculate flux limiter and W/D flag
       !!----------------------------------------------------------------------
-      REAL(wp)                , INTENT(in   ) ::   rdtbt    ! ocean time-step index
+      REAL(wp)                , INTENT(in   ) ::   pdt    ! external mode time-step [s]
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc  
       !
       INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
       INTEGER  ::   jflag               ! local integer
-      REAL(wp) ::   z2dt
       REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
       REAL(wp) ::   zzflxp, zzflxn      ! local scalars
@@ -298 +297 @@
       jflag  = 0
       zdepwd = 50._wp   ! maximum depth that ocean cells can have W/D processes
-      !
-      z2dt = rdtbt   
       !
       zflxp(:,:)   = 0._wp
@@ -347 +344 @@
                   &   + min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
           
-               zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
-               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - z2dt * zssh_frc(ji,jj)
+               zdep1 = (zzflxp + zzflxn) * pdt / ztmp
+               zdep2 = ht_0(ji,jj) + sshn_e(ji,jj) - rn_wdmin1 - pdt * zssh_frc(ji,jj)
           
                IF(zdep1 > zdep2) THEN
-                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
-                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zzflxp * z2dt )
+                 zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * pdt ) / ( zflxp(ji,jj) * pdt )
+                 !zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * pdt ) / ( zzflxp * pdt )
                  ! flag if the limiter has been used but stop flagging if the only
                  ! changes have zeroed the coefficient since further iterations will

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traadv.F90

@@ -92 +92 @@
       IF( ln_timing )   CALL timing_start('tra_adv')
       !
-      !                                          ! set time step
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =         rdt   ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp * rdt   ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-      !
       !                                         !==  effective transport  ==!
       zun(:,:,jpk) = 0._wp

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf.F90

@@ -55 +55 @@
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !!
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd   ! 4D workspace
       !!----------------------------------------------------------------------
       !
@@ -61 +61 @@
       !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ALLOCATE( ztrd(jpi,jpj,jpk,jpts) ) 
+         ztrd(:,:,:,:) = tsa(:,:,:,:) 
       ENDIF
       !
@@ -78 +77 @@
       !
       IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
-         DEALLOCATE( ztrdt, ztrds ) 
+         ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:)
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrd(:,:,:,jp_tem) )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrd(:,:,:,jp_sal) )
+         DEALLOCATE( ztrd ) 
       ENDIF
       !                                        !* print mean trends (used for debugging)

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf_iso.F90

@@ -108 +108 @@
       REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
       REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
-      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign                 !   -      -
       REAL(wp), DIMENSION(jpi,jpj)     ::   zdkt, zdk1t, z2d
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zftu, zftv, ztfw 
@@ -127 +127 @@
       IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
          &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
-      !
-      !                                            ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
-      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
-      ENDIF
-      z1_2dt = 1._wp / z2dt
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -191 +185 @@
                      DO ji = 1, fs_jpim1
                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
-                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
-                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                        zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt
                      END DO
                   END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traldf_triad.F90

@@ -85 +85 @@
       INTEGER  ::  ip,jp,kp         ! dummy loop indices
       INTEGER  ::  ierr            ! local integer
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      -
-      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      -
+      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
+      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign         !   -      -
       !
       REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
@@ -110 +110 @@
       l_hst = .FALSE.
       l_ptr = .FALSE.
-      IF( cdtype == 'TRA' .AND. ln_diaptr )                                                 l_ptr = .TRUE. 
-      IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
-         &                        iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")  ) )   l_hst = .TRUE.
-      !
-      !                                                        ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdt      ! at nit000   (Euler)
-      ELSE                                        ;   z2dt = 2.* rdt      !             (Leapfrog)
+      IF( cdtype == 'TRA' ) THEN
+         IF ( ln_diaptr )                                                 l_ptr = .TRUE. 
+         IF ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+            & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr")      )   l_hst = .TRUE.
       ENDIF
-      z1_2dt = 1._wp / z2dt
       !
       IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
@@ -202 +198 @@
                      DO ji = 1, fs_jpim1
                         ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk)
-                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
-                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                        zcoef0 = r2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * r1_2dt
                      END DO
                   END DO

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/tranpc.F90

@@ -65 +65 @@
       LOGICAL  ::   l_bottom_reached, l_column_treated
       REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
-      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
+      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw
       REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0)
       REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point...
@@ -301 +301 @@
          !
          IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
-            z1_r2dt = 1._wp / (2._wp * rdt)
-            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * z1_r2dt
-            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * z1_r2dt
+            ztrdt(:,:,:) = ( tsa(:,:,:,jp_tem) - ztrdt(:,:,:) ) * r1_2dt
+            ztrds(:,:,:) = ( tsa(:,:,:,jp_sal) - ztrds(:,:,:) ) * r1_2dt
             CALL trd_tra( kt, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, 'TRA', jp_sal, jptra_npc, ztrds )

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/tranxt.F90

@@ -111 +111 @@
       IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
  
-      ! set time step size (Euler/Leapfrog)
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
-      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
-      ENDIF
-
-      ! trends computation initialisation
-      IF( l_trdtra )   THEN                    
+      IF( l_trdtra )   THEN               ! trends computation initialisation
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,jpk) = 0._wp
@@ -142 +136 @@
       ENDIF
 
-      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
+      IF( l_1st_euler ) THEN        ! Euler time-stepping at first time-step (only swap)
          DO jn = 1, jpts
             DO jk = 1, jpkm1
@@ -156 +150 @@
          END IF
          !
-      ELSE                                            ! Leap-Frog + Asselin filter time stepping
+      ELSE                          ! Leap-Frog + Asselin filter time stepping
          !
          IF( ln_linssh ) THEN   ;   CALL tra_nxt_fix( kt, nit000,      'TRA', tsb, tsn, tsa, jpts )  ! linear free surface 
@@ -163 +157 @@
          ENDIF
          !
-         CALL lbc_lnk_multi( tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., &
-                  &          tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., &
-                  &          tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  )
+!!gm
+!         CALL lbc_lnk_multi( tsb(:,:,:,jp_tem), 'T', 1., tsb(:,:,:,jp_sal), 'T', 1., &
+!            &                tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., &
+!            &                tsa(:,:,:,jp_tem), 'T', 1., tsa(:,:,:,jp_sal), 'T', 1.  )
+!!gm
+         CALL lbc_lnk_multi( tsb, 'T', 1., tsn, 'T', 1., tsa, 'T', 1.  )
+!!gm
          !
       ENDIF     

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/traqsr.F90

@@ -133 +133 @@
       !                         !-----------------------------------!
       IF( kt == nit000 ) THEN          !==  1st time step  ==!
-!!gm case neuler  not taken into account....
-         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0 ) THEN    ! read in restart
+         IF( ln_rstart .AND. iom_varid( numror, 'qsr_hc_b', ldstop = .FALSE. ) > 0  .AND. .NOT.l_1st_euler ) THEN    ! read in restart
             IF(lwp) WRITE(numout,*) '          nit000-1 qsr tracer content forcing field read in the restart file'
             z1_2 = 0.5_wp

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRA/trazdf.F90

@@ -52 +52 @@
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
-      INTEGER  ::   jk   ! Dummy loop indices
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
+      INTEGER  ::   jk, jts   ! Dummy loop indices
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   ztrd   ! 4D workspace
       !!---------------------------------------------------------------------
       !
@@ -64 +64 @@
       ENDIF
       !
-      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =      rdt   ! at nit000, =   rdt (restarting with Euler time stepping)
-      ELSEIF( kt <= nit000 + 1           ) THEN   ;   r2dt = 2. * rdt   ! otherwise, = 2 rdt (leapfrog)
-      ENDIF
-      !
       IF( l_trdtra )   THEN                  !* Save ta and sa trends
-         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ALLOCATE( ztrd(jpi,jpj,jpk,jpts) )
+         ztrd(:,:,:,:) = tsa(:,:,:,:)
       ENDIF
       !
@@ -85 +80 @@
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &
-               &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &
-              &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)
+         DO jts = 1, jpts
+            DO jk = 1, jpkm1
+               ztrd(:,:,jk,jts) = ( ( tsa(:,:,jk,jts)*e3t_a(:,:,jk) - tsb(:,:,jk,jts)*e3t_b(:,:,jk) ) / (e3t_n(:,:,jk)*r2dt) )   &
+                  &            - ztrd(:,:,jk,jts)
+            END DO
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
-         CALL lbc_lnk_multi( ztrdt, 'T', 1. , ztrds, 'T', 1. )
+         CALL lbc_lnk( ztrd, 'T', 1. )
 !!gm
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
-         DEALLOCATE( ztrdt , ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrd(:,:,:,jp_tem) )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrd(:,:,:,jp_sal) )
+         DEALLOCATE( ztrd )
       ENDIF
       !                                          ! print mean trends (used for debugging)

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/TRD/trdtra.F90

@@ -237 +237 @@
       INTEGER                   , INTENT(in   ) ::   kt      ! time step
       !!----------------------------------------------------------------------
-
-      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdt (restart with Euler time stepping)
-      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdt (leapfrog)
-      ENDIF
 
       !                   ! 3D output of tracers trends using IOM interface

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/nemogcm.F90

@@ -151 +151 @@
       !                            !==   time stepping   ==!
       !                            !-----------------------!
+      !
+      !                                               !== set the model time-step  ==!
+      !
+      IF( l_1st_euler ) THEN   ;   r2dt =         rn_rdt   ;   l_1st_euler = .TRUE.    ! start or restart with Euler 1st time-step
+      ELSE                     ;   r2dt = 2._wp * rn_rdt   ;   l_1st_euler = .FALSE.   ! restart with leapfrog
+      ENDIF
+      r1_2dt = 1._wp / r2dt
+      ! NB: if l_1st_euler=T, r2dt will be set to 2*rdt at the end of the 1st time-step (in step.F90)
+      !     Done here (not in domain.F90) as in ASM initialization an Euler 1st time step can be forced
+      !
+      !
       istp = nit000
       !

NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/step.F90

@@ -34 +34 @@
 
    !!----------------------------------------------------------------------
-   !!   stp             : OPA system time-stepping
-   !!----------------------------------------------------------------------
-   USE step_oce         ! time stepping definition modules
+   !!   stp           : OPA system time-stepping
+   !!----------------------------------------------------------------------
+   USE step_oce       ! time stepping definition modules
    !
-   USE iom              ! xIOs server
+   USE iom            ! xIOs server
 
    IMPLICIT NONE
@@ -323 +323 @@
 #endif
       !
+      IF( l_1st_euler ) THEN
+         r2dt   = 2._wp * rn_rdt                            ! recover Leap-frog time-step
+         r1_2dt = 1._wp / r2dt
+         l_1st_euler = .FALSE.
+      ENDIF
+      !
       IF( ln_timing )   CALL timing_stop('stp')
       !