diff --git a/arch/MERCATOR/arch-X64_KARA_GCC_OMPI_DEBUG.fcm b/arch/MERCATOR/arch-X64_KARA_GCC_OMPI_DEBUG.fcm
index bf2cb6ce28411c008b39f3182a72604096bd7be0..dad25a47a46c45ec5c16edf5254ab27e98d8918f 100644
--- a/arch/MERCATOR/arch-X64_KARA_GCC_OMPI_DEBUG.fcm
+++ b/arch/MERCATOR/arch-X64_KARA_GCC_OMPI_DEBUG.fcm
@@ -49,7 +49,7 @@
 
 %CPP                 cpp -Dkey_nosignedzero
 %FC                  mpif90 -c -cpp
-%FCFLAGS             -fdefault-real-8 -O0 -g -fbacktrace -ftree-vectorize -funroll-all-loops -march=skylake-avx512 -ffree-line-length-none -fcheck=all -finit-real=nan -Wno-missing-include-dirs
+%FCFLAGS             -fdefault-real-8 -O0 -g -fbacktrace -ftree-vectorize -funroll-all-loops -march=skylake-avx512 -ffree-line-length-none -fcheck=all -finit-real=snan -Wno-missing-include-dirs -ffpe-trap=invalid,zero,overflow
 %FFLAGS              %FCFLAGS
 %LD                  mpif90
 %LDFLAGS             -lstdc++
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index 31674026d3baecba5ce200124b1c99e82d9cc903..19c4f23474597ea6e151299d54d21fb8ddc09cd7 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -1538,6 +1538,20 @@
 !-----------------------------------------------------------------------
 &namsto        ! Stochastic parametrization of EOS                      (default: OFF)
 !-----------------------------------------------------------------------
+   ln_sto_ldf  = .false.   ! stochastic lateral diffusion
+   rn_ldf_std  =  0.1      ! lateral diffusion standard deviation (in percent)
+   rn_ldf_tcor =  1440.    ! lateral diffusion correlation timescale (in timesteps)
+   ln_sto_hpg  = .false.   ! stochastic pressure gradient
+   rn_hpg_std  =  0.1      ! density gradient standard deviation (in percent)
+   rn_hpg_tcor =  1440.    ! density gradient correlation timescale (in timesteps)
+   ln_sto_pstar  = .false. ! stochastic ice strength
+   rn_pstar_std  =  0.1    ! ice strength standard deviation (in percent)
+   rn_pstar_tcor =  1440.  ! ice strength correlation timescale (in timesteps)
+   nn_pstar_ord  =  1      ! order of autoregressive processes
+   nn_pstar_flt  =  0      ! passes of Laplacian filter
+   ln_sto_trd  = .false.   ! stochastic model trend
+   rn_trd_std  =  0.1      ! trend standard deviation (in percent)
+   rn_trd_tcor =  1440.    ! trend correlation timescale (in timesteps)
    ln_sto_eos  = .false.   ! stochastic equation of state
    nn_sto_eos  = 1         ! number of independent random walks
    rn_eos_stdxy = 1.4       ! random walk horz. standard deviation (in grid points)
@@ -1546,6 +1560,14 @@
    nn_eos_ord  = 1         ! order of autoregressive processes
    nn_eos_flt  = 0         ! passes of Laplacian filter
    rn_eos_lim  = 2.0       ! limitation factor (default = 3.0)
+   ln_sto_trc  = .false.   ! stochastic tracer dynamics
+   nn_sto_trc  = 1         ! number of independent random walks
+   rn_trc_stdxy = 1.4      ! random walk horz. standard deviation (in grid points)
+   rn_trc_stdz = 0.7       ! random walk vert. standard deviation (in grid points)
+   rn_trc_tcor = 1440.     ! random walk time correlation (in timesteps)
+   nn_trc_ord  = 1         ! order of autoregressive processes
+   nn_trc_flt  = 0         ! passes of Laplacian filter
+   rn_trc_lim  = 3.0       ! limitation factor (default = 3.0)
    ln_rststo   = .false.   ! start from mean parameter (F) or from restart file (T)
    ln_rstseed  = .true.    ! read seed of RNG from restart file
    cn_storst_in  = "restart_sto" !  suffix of stochastic parameter restart file (input)
diff --git a/src/ICE/iceitd.F90 b/src/ICE/iceitd.F90
index 1c0df3661a7cceadd32b1e22fe710b9a1ae5db81..f3e662748ae988a2ee5388c3085fb297b05b9efa 100644
--- a/src/ICE/iceitd.F90
+++ b/src/ICE/iceitd.F90
@@ -709,6 +709,8 @@ CONTAINS
       NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin, rn_himax
       !!------------------------------------------------------------------
       !
+      rn_catbnd(:) =  0._wp ! Circumvent possible initialization by compiler
+                            ! to prevent from errors when writing output
       READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namitd in reference namelist' )
       READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
diff --git a/src/NST/agrif_oce_interp.F90 b/src/NST/agrif_oce_interp.F90
index 4d4535e4a3a4834a507532f30f5bd751c4e12e68..6cd1e9ff2a3c5787c1f617b75c8b20bb3adb5fd7 100644
--- a/src/NST/agrif_oce_interp.F90
+++ b/src/NST/agrif_oce_interp.F90
@@ -918,7 +918,7 @@ CONTAINS
          IF( l_ini_child )   Kmm_a = Kbb_a  
 
          DO jn = 1,jpts
-            DO jk=k1,k2
+            DO jk=k1,k2-1
                DO jj=j1,j2
                  DO ji=i1,i2
                        ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
@@ -934,7 +934,7 @@ CONTAINS
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3w(ji,jj,k1,Kmm_a)
-                  DO jk=k1+1,k2
+                  DO jk=k1+1,k2-1
                      ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
                         & ( ptab(ji,jj,jk-1,jpts+1) + e3w(ji,jj,jk,Kmm_a) ) 
                   END DO
@@ -1047,7 +1047,7 @@ CONTAINS
             ENDIF
 
             DO jn =1, jpts
-               ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)
+               ts(i1:i2,j1:j2,1:jpkm1,jn,Krhs_a) = ptab(i1:i2,j1:j2,1:jpkm1,jn)*tmask(i1:i2,j1:j2,1:jpkm1)
             END DO
          ENDIF
 
@@ -1121,7 +1121,7 @@ CONTAINS
          item = Kmm_a
          IF( l_ini_child )   Kmm_a = Kbb_a     
 
-         DO jk=1,jpk
+         DO jk=k1,k2-1
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)*umask(ji,jj,jk)) 
@@ -1215,7 +1215,7 @@ CONTAINS
          item = Kmm_a
          IF( l_ini_child )   Kmm_a = Kbb_a     
        
-         DO jk=k1,k2
+         DO jk=k1,k2-1
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)*vmask(ji,jj,jk))
diff --git a/src/NST/agrif_oce_sponge.F90 b/src/NST/agrif_oce_sponge.F90
index 0e728c4c29bad106535c1d03cb5246747e4f8cd7..5585854cee2220f51c3bd89f5e74ac30967a3db0 100644
--- a/src/NST/agrif_oce_sponge.F90
+++ b/src/NST/agrif_oce_sponge.F90
@@ -398,7 +398,7 @@ CONTAINS
       !
       IF( before ) THEN
          DO jn = 1, jpts
-            DO jk=k1,k2
+            DO jk=k1,k2-1
                DO jj=j1,j2
                   DO ji=i1,i2
                      ! JC: masking is mandatory here: before tracer field seems 
@@ -416,7 +416,7 @@ CONTAINS
             DO jj=j1,j2
                DO ji=i1,i2
                   tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3w(ji,jj,k1,Kbb_a)
-                  DO jk=k1+1,k2
+                  DO jk=k1+1,k2-1
                      tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
                         & ( tabres(ji,jj,jk-1,jpts+1) + e3w(ji,jj,jk,Kbb_a) )
                   END DO
@@ -627,7 +627,7 @@ CONTAINS
       !!---------------------------------------------    
       !
       IF( before ) THEN
-         DO jk=k1,k2
+         DO jk=k1,k2-1
             DO jj=j1,j2
                DO ji=i1,i2
                   tabres(ji,jj,jk,m1) = e2u(ji,jj) * e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) * umask(ji,jj,jk)
@@ -672,7 +672,7 @@ CONTAINS
             END DO
          ENDIF
          !
-         ubdiff(i1:i2,j1:j2,1:jpk) = (uu(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*umask(i1:i2,j1:j2,1:jpk)
+         ubdiff(i1:i2,j1:j2,1:jpkm1) = (uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpkm1))*umask(i1:i2,j1:j2,1:jpkm1)
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !                                             ! ===============
@@ -774,7 +774,7 @@ CONTAINS
       !!--------------------------------------------- 
       
       IF( before ) THEN 
-         DO jk=k1,k2
+         DO jk=k1,k2-1
             DO jj=j1,j2
                DO ji=i1,i2
                   tabres(ji,jj,jk,m1) = e1v(ji,jj) * e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) * vmask(ji,jj,jk)
@@ -820,7 +820,7 @@ CONTAINS
             END DO
          ENDIF
          !
-         vbdiff(i1:i2,j1:j2,1:jpk) = (vv(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*vmask(i1:i2,j1:j2,1:jpk)
+         vbdiff(i1:i2,j1:j2,1:jpkm1) = (vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpkm1))*vmask(i1:i2,j1:j2,1:jpkm1)
          !
          DO jk = 1, jpkm1                                 ! Horizontal slab
             !                                             ! ===============
diff --git a/src/NST/agrif_top_interp.F90 b/src/NST/agrif_top_interp.F90
index ad7770e4372e5d17aa598ef5c82561164e93ae66..e2314a077cae5495c04fdf8711e2ae8c11c8e449 100644
--- a/src/NST/agrif_top_interp.F90
+++ b/src/NST/agrif_top_interp.F90
@@ -77,7 +77,7 @@ CONTAINS
          IF( l_ini_child )   Kmm_a = Kbb_a  
 
          DO jn = 1,jptra
-            DO jk=k1,k2
+            DO jk=k1,k2-1
                DO jj=j1,j2
                  DO ji=i1,i2
                        ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)
@@ -92,7 +92,7 @@ CONTAINS
             DO jj=j1,j2
                DO ji=i1,i2
                   ptab(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3w(ji,jj,k1,Kmm_a)
-                  DO jk=k1+1,k2
+                  DO jk=k1+1,k2-1
                      ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
                         & ( ptab(ji,jj,jk-1,jptra+1) + e3w(ji,jj,jk,Kmm_a) )
                   END DO
@@ -206,7 +206,7 @@ CONTAINS
             ENDIF
 
             DO jn=1, jptra
-                tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
+                tr(i1:i2,j1:j2,1:jpkm1,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpkm1,jn)*tmask(i1:i2,j1:j2,1:jpkm1) 
             END DO
          ENDIF
 
diff --git a/src/NST/agrif_top_sponge.F90 b/src/NST/agrif_top_sponge.F90
index 0b29cb59c2d35f57282dc3db7134cac5aaa4429e..e40c7063ca9bdbef336bc0a0c6d641f825e1c1fa 100644
--- a/src/NST/agrif_top_sponge.F90
+++ b/src/NST/agrif_top_sponge.F90
@@ -86,7 +86,7 @@ CONTAINS
       !
       IF( before ) THEN
          DO jn = 1, jptra
-            DO jk=k1,k2
+            DO jk=k1,k2-1
                DO jj=j1,j2
                   DO ji=i1,i2
                      tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a) * tmask(ji,jj,jk) 
@@ -102,7 +102,7 @@ CONTAINS
             DO jj=j1,j2
                DO ji=i1,i2
                   tabres(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3w(ji,jj,k1,Kbb_a)
-                  DO jk=k1+1,k2
+                  DO jk=k1+1,k2-1
                      tabres(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
                         & ( tabres(ji,jj,jk-1,jptra+1) + e3w(ji,jj,jk,Kbb_a) )
                   END DO
diff --git a/src/OCE/ZDF/zdfiwm.F90 b/src/OCE/ZDF/zdfiwm.F90
index 30a2dfd71a5ece2e41e15b323ac6f09b165ee8ec..4df01326275a4a1bd4fc640b8ad7367383f76370 100644
--- a/src/OCE/ZDF/zdfiwm.F90
+++ b/src/OCE/ZDF/zdfiwm.F90
@@ -145,7 +145,9 @@ CONTAINS
       !!----------------------------------------------------------------------
       !
       !                       !* Initialize appropriately certain variables
-      zav_ratio(:,:,:) = 1._wp * wmask(:,:,:)  ! important to set it to 1 here 
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         zav_ratio(ji,jj,jk) = 1._wp * wmask(ji,jj,jk)  ! important to set it to 1 here 
+      END_3D
       IF( iom_use("emix_iwm") )                         zemx_iwm (:,:,:) = 0._wp
       IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl )   zav_wave (:,:,:) = 0._wp
       !
diff --git a/src/OCE/timing.F90 b/src/OCE/timing.F90
index a7e3ca0e5c0469b0ac976444667ed4ee67266a74..c63003d8b66e5ffb901c471502339be6d98c2eeb 100644
--- a/src/OCE/timing.F90
+++ b/src/OCE/timing.F90
@@ -387,15 +387,16 @@ CONTAINS
          WRITE(numtime,*) '    warning: includes restarts writing time if output before nitend... '
          WRITE(numtime,*) ' '
          DO ji = 1, jpnij
+            zperc = 0._wp ; zsypd = 0._wp
             ztot = SUM( timing_glob(4*ji-3:4*ji-1) )
             WRITE(numtime,'(A28,F11.6,            A34,I8)') 'Computing       time : ',timing_glob(4*ji-3), ' on MPI rank : ', ji
-            IF ( ztot /= 0. ) zperc = timing_glob(4*ji-2) / ztot * 100.
+            IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-2) / ztot * 100.
             WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting lbc_lnk time : ',timing_glob(4*ji-2)   &
                &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
-            IF ( ztot /= 0. ) zperc = timing_glob(4*ji-1) / ztot * 100.
+            IF ( ztot /= 0._wp ) zperc = timing_glob(4*ji-1) / ztot * 100.
             WRITE(numtime,'(A28,F11.6,A2, F4.1,A3,A25,I8)') 'Waiting  global time : ',timing_glob(4*ji-1)   &
                &                                                         , ' (',      zperc,' %)',   ' on MPI rank : ', ji
-            zsypd = rn_Dt * REAL(nitend-nit000-1, wp) / (timing_glob(4*ji) * 365.)
+            IF ( timing_glob(4*ji) /= 0._wp ) zsypd = rn_Dt * REAL(nitend-nit000-1, wp) / (timing_glob(4*ji) * 365.)
             WRITE(numtime,'(A28,F11.6,A7,F10.3,A2,A15,I8)') 'Total           time : ',timing_glob(4*ji  )   &
                &                                                         , ' (SYPD: ', zsypd, ')',   ' on MPI rank : ', ji
          END DO
diff --git a/tools/DOMAINcfg/namelist_cfg b/tools/DOMAINcfg/namelist_cfg
index 8615fede68c3a0a58103a2e4f177baf01b9037c8..cf41b1a1d3950f2ddb4629c1b9370050c2e81b63 100644
--- a/tools/DOMAINcfg/namelist_cfg
+++ b/tools/DOMAINcfg/namelist_cfg
@@ -54,6 +54,8 @@
    !                       !      ===>>> will become the only possibility in v4.0
    !                       ! =F : e3 analytical derivative of depth function
    !                       !      only there for backward compatibility test with v3.6
+      !                      ! if ln_e3_dep = T
+      ln_dept_mid = .false.  ! =T : set T points in the middle of cells
    !                       !
    cp_cfg      =  "orca"   !  name of the configuration
    jp_cfg      =       2   !  resolution of the configuration