Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 5845)
@@ -60,5 +60,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -67,5 +66,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -187,5 +185,5 @@
                !    closure : flux grazing is redistributed below level jpkbio
                zzoobod = tmminz * zzoo * zzoo
-               xksi(ji,jj) = xksi(ji,jj) + (1-fdbod) * zzoobod * fse3t(ji,jj,jk)
+               xksi(ji,jj) = xksi(ji,jj) + (1-fdbod) * zzoobod * e3t_n(ji,jj,jk)
                zboddet = fdbod * zzoobod
 
@@ -242,5 +240,5 @@
                 IF( ln_diatrc .OR. lk_iomput ) THEN
                   ! convert fluxes in per day
-                  ze3t = fse3t(ji,jj,jk) * 86400.
+                  ze3t = e3t_n(ji,jj,jk) * 86400.
                   zw2d(ji,jj,1)  = zw2d(ji,jj,1)  + zno3phy * ze3t
                   zw2d(ji,jj,2)  = zw2d(ji,jj,2)  + znh4phy * ze3t
@@ -363,5 +361,5 @@
                 IF( ln_diatrc .OR. lk_iomput ) THEN
                   ! convert fluxes in per day
-                  ze3t = fse3t(ji,jj,jk) * 86400.
+                  ze3t = e3t_n(ji,jj,jk) * 86400.
                   zw2d(ji,jj,1)  = zw2d(ji,jj,1)  + zno3phy * ze3t
                   zw2d(ji,jj,2)  = zw2d(ji,jj,2)  + znh4phy * ze3t
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 5845)
@@ -42,5 +42,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -95,5 +94,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
-               ze3t = 1. / fse3t(ji,jj,jk)
+               ze3t = 1. / e3t_n(ji,jj,jk)
                tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + ze3t * dmin3(ji,jj,jk) * xksi(ji,jj)
             END DO
@@ -110,5 +109,5 @@
          DO ji = fs_2, fs_jpim1
             ikt = mbkt(ji,jj) 
-            tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) + sedlam * sedpocn(ji,jj) / fse3t(ji,jj,ikt) 
+            tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) + sedlam * sedpocn(ji,jj) / e3t_n(ji,jj,ikt) 
             ! Deposition of organic matter in the sediment
             zwork = vsed * trn(ji,jj,ikt,jpdet)
@@ -121,5 +120,5 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1
-            tra(ji,jj,1,jpno3) = tra(ji,jj,1,jpno3) + zgeolpoc * cmask(ji,jj) / areacot / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpno3) = tra(ji,jj,1,jpno3) + zgeolpoc * cmask(ji,jj) / areacot / e3t_n(ji,jj,1)
          END DO
       END DO
@@ -212,6 +211,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zfluo = ( fsdepw(ji,jj,jk  ) / fsdepw(ji,jj,jpkb) )**xhr
-               zfluu = ( fsdepw(ji,jj,jk+1) / fsdepw(ji,jj,jpkb) )**xhr
+               zfluo = ( gdepw_n(ji,jj,jk  ) / gdepw_n(ji,jj,jpkb) )**xhr
+               zfluu = ( gdepw_n(ji,jj,jk+1) / gdepw_n(ji,jj,jpkb) )**xhr
                IF( zfluo.GT.1. )   zfluo = 1._wp
                zdm0(ji,jj,jk) = zfluo - zfluu
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5845)
@@ -40,6 +40,4 @@
    REAL(wp), PUBLIC ::  reddom    ! redfield ratio (C:N) for DOM
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -47,5 +45,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -105,6 +102,6 @@
                zkr  = xkr0 + xkrp * EXP( xlr * zpig )
                zkg  = xkg0 + xkgp * EXP( xlg * zpig )
-               zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * fse3t(ji,jj,jk-1) )
-               zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * fse3t(ji,jj,jk-1) )
+               zparr(ji,jj,jk) = zparr(ji,jj,jk-1) * EXP( -zkr * e3t_n(ji,jj,jk-1) )
+               zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t_n(ji,jj,jk-1) )
             END DO
         END DO
@@ -116,6 +113,6 @@
                zkr  = xkr0 + xkrp * EXP( xlr * zpig )
                zkg  = xkg0 + xkgp * EXP( xlg * zpig )
-               zparr(ji,jj,jk) = zparr(ji,jj,jk) / ( zkr * fse3t(ji,jj,jk) ) * ( 1 - EXP( -zkr * fse3t(ji,jj,jk) ) )
-               zparg(ji,jj,jk) = zparg(ji,jj,jk) / ( zkg * fse3t(ji,jj,jk) ) * ( 1 - EXP( -zkg * fse3t(ji,jj,jk) ) )
+               zparr(ji,jj,jk) = zparr(ji,jj,jk) / ( zkr * e3t_n(ji,jj,jk) ) * ( 1 - EXP( -zkr * e3t_n(ji,jj,jk) ) )
+               zparg(ji,jj,jk) = zparg(ji,jj,jk) / ( zkg * e3t_n(ji,jj,jk) ) * ( 1 - EXP( -zkg * e3t_n(ji,jj,jk) ) )
                etot (ji,jj,jk) = MAX( zparr(ji,jj,jk) + zparg(ji,jj,jk), 1.e-15 )
             END DO
@@ -138,5 +135,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            heup(ji,jj) = fsdepw(ji,jj,neln(ji,jj))
+            heup(ji,jj) = gdepw_n(ji,jj,neln(ji,jj))
          END DO
       END DO 
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zsed.F90	(revision 5845)
@@ -34,6 +34,4 @@
    REAL(wp), PUBLIC ::   xhr         ! coeff for martin''s remineralisation profile
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -102,5 +100,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ztra(ji,jj,jk)  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+               ztra(ji,jj,jk)  = - ( zwork(ji,jj,jk) - zwork(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
                tra(ji,jj,jk,jpdet) = tra(ji,jj,jk,jpdet) + ztra(ji,jj,jk) 
             END DO
@@ -111,7 +109,7 @@
          IF( iom_use( "TDETSED" ) ) THEN
             CALL wrk_alloc( jpi, jpj, zw2d )
-            zw2d(:,:) =  ztra(:,:,1) * fse3t(:,:,1) * 86400.
+            zw2d(:,:) =  ztra(:,:,1) * e3t_n(:,:,1) * 86400.
             DO jk = 2, jpkm1
-               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * fse3t(:,:,jk) * 86400.
+               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * e3t_n(:,:,jk) * 86400.
             END DO
             CALL iom_put( "TDETSED", zw2d )
@@ -121,7 +119,7 @@
          IF( ln_diatrc ) THEN 
             CALL wrk_alloc( jpi, jpj, zw2d )
-            zw2d(:,:) =  ztra(:,:,1) * fse3t(:,:,1) * 86400.
+            zw2d(:,:) =  ztra(:,:,1) * e3t_n(:,:,1) * 86400.
             DO jk = 2, jpkm1
-               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * fse3t(:,:,jk) * 86400.
+               zw2d(:,:) = zw2d(:,:) + ztra(:,:,jk) * e3t_n(:,:,jk) * 86400.
             END DO
             trc2d(:,:,jp_pcs0_2d + 7) = zw2d(:,:)
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 5845)
@@ -34,6 +34,4 @@
    PUBLIC  p4z_bio    
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -70,10 +68,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) > hmld(ji,jj) )   xdiss(ji,jj,jk) = 0.01
+!!gm  :  use nmln  and test on jk ...  less memory acces
+               IF( gdepw_n(ji,jj,jk+1) > hmld(ji,jj) )   xdiss(ji,jj,jk) = 0.01
             END DO 
          END DO
       END DO
 
-          
       CALL p4z_opt  ( kt, knt )     ! Optic: PAR in the water column
       CALL p4z_sink ( kt, knt )     ! vertical flux of particulate organic matter
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 5845)
@@ -164,6 +164,4 @@
    REAL(wp) :: devk55  = 0.3692E-3      
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -244,5 +242,4 @@
 
 
-
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
@@ -252,5 +249,5 @@
 
                ! SET PRESSION
-               zpres   = 1.025e-1 * fsdept(ji,jj,jk)
+               zpres   = 1.025e-1 * gdept_n(ji,jj,jk)
 
                ! SET ABSOLUTE TEMPERATURE
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 5845)
@@ -30,15 +30,13 @@
    PUBLIC   p4z_fechem_init ! called in trcsms_pisces.F90
 
-   !! * Shared module variables
-   LOGICAL          ::  ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
-   LOGICAL          ::  ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
-   REAL(wp), PUBLIC ::  xlam1        !: scavenging rate of Iron 
-   REAL(wp), PUBLIC ::  xlamdust     !: scavenging rate of Iron by dust 
-   REAL(wp), PUBLIC ::  ligand       !: ligand concentration in the ocean 
-
+   LOGICAL          ::   ln_fechem    !: boolean for complex iron chemistry following Tagliabue and voelker
+   LOGICAL          ::   ln_ligvar    !: boolean for variable ligand concentration following Tagliabue and voelker
+   REAL(wp), PUBLIC ::   xlam1        !: scavenging rate of Iron 
+   REAL(wp), PUBLIC ::   xlamdust     !: scavenging rate of Iron by dust 
+   REAL(wp), PUBLIC ::   ligand       !: ligand concentration in the ocean 
+
+!!gm Not DOCTOR norm !!!
    REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -61,8 +59,8 @@
       !!                    and one particulate form (ln_fechem)
       !!---------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
       !
       INTEGER  ::   ji, jj, jk, jic
+      CHARACTER (len=25) :: charout
       REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
       REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll
@@ -79,16 +77,14 @@
       REAL(wp) :: ztfe, zoxy
       REAL(wp) :: zstep
-      CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_fechem')
       !
-      ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig )
+      CALL wrk_alloc( jpi,jpj,jpk,   zFe3, zFeL1, zTL1, ztotlig )
       zFe3 (:,:,:) = 0.
       zFeL1(:,:,:) = 0.
       zTL1 (:,:,:) = 0.
       IF( ln_fechem ) THEN
-         CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+         CALL wrk_alloc( jpi,jpj,jpk,   zFe2, zFeL2, zTL2, zFeP )
          zFe2 (:,:,:) = 0.
          zFeL2(:,:,:) = 0.
@@ -253,5 +249,6 @@
                zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
                zlamfac = MIN( 1.  , zlamfac )
-               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
+!!gm very small BUG :  it is unlikely but possible that gdept_n = 0  .....
+               zdep    = MIN( 1., 1000. / gdept_n(ji,jj,jk) )
                zlam1b  = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
                zcoag   = zfeequi * zlam1b * zstep + 1E-4 * ( 1. - zlamfac ) * zdep * zstep * trb(ji,jj,jk,jpfer)
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 5845)
@@ -59,6 +59,4 @@
    REAL(wp) ::  xconv  = 0.01_wp / 3600._wp !: coefficients for conversion 
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -182,5 +180,5 @@
             oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! compute the trend
-            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / e3t_n(ji,jj,1)
 
             ! Compute O2 flux 
@@ -188,5 +186,5 @@
             zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
             zoflx(ji,jj) = zfld16 - zflu16
-            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / e3t_n(ji,jj,1)
          END DO
       END DO
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5845)
@@ -51,6 +51,4 @@
    REAL(wp), DIMENSION(3,61), PUBLIC ::   xkrgb   !: tabulated attenuation coefficients for RGB absorption
    
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -101,7 +99,7 @@
                irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
                !                                                         
-               ekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk)
-               ekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk)
-               ekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk)
+               ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t_n(ji,jj,jk)
+               ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t_n(ji,jj,jk)
+               ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t_n(ji,jj,jk)
             END DO
          END DO
@@ -162,5 +160,5 @@
                  neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
                  !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
-                 heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
+                 heup(ji,jj) = gdepw_n(ji,jj,jk+1)     ! Euphotic layer depth
               ENDIF
            END DO
@@ -179,10 +177,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * fse3t(ji,jj,jk) ! remineralisation
-                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * fse3t(ji,jj,jk) ! production
-                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
-                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
-                  zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
+               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t_n(ji,jj,jk) ! remineralisation
+                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * e3t_n(ji,jj,jk) ! production
+                  zdepmoy(ji,jj) = zdepmoy(ji,jj) +                       e3t_n(ji,jj,jk)
                ENDIF
             END DO
@@ -196,5 +194,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                   z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
                   emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
@@ -260,5 +258,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * xsi0r )
+                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t_n(ji,jj,jk-1) * xsi0r )
                   pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) )
                   pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) )
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5845)
@@ -54,6 +54,4 @@
    REAL(wp) :: texcret2               !: 1 - excret2        
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -277,5 +275,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+               IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                   zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * zmixnano(ji,jj)
                   zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
@@ -321,5 +319,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
+                  IF( gdepw_n(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                      zprnch(ji,jj,jk) = zprnch(ji,jj,jk) * zmixnano(ji,jj)
                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
@@ -462,5 +460,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-               zw2d(:,:) = zw2d(:,:) + zprorca (:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
+               zw2d(:,:) = zw2d(:,:) + zprorca (:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
              ENDDO
              CALL iom_put( "INTPPPHY" , zw2d )
@@ -468,5 +466,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated  primary produc. by diatom
+                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated  primary produc. by diatom
              ENDDO
              CALL iom_put( "INTPPPHY2" , zw2d )
@@ -475,5 +473,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + ( zprorca(:,:,jk) + zprorcad(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
+                zw2d(:,:) = zw2d(:,:) + ( zprorca(:,:,jk) + zprorcad(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
              ENDDO
              CALL iom_put( "INTPP" , zw2d )
@@ -482,5 +480,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + ( zpronew(:,:,jk) + zpronewd(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
+                zw2d(:,:) = zw2d(:,:) + ( zpronew(:,:,jk) + zpronewd(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
              ENDDO
              CALL iom_put( "INTPNEW" , zw2d )
@@ -489,5 +487,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + ( zprofen(:,:,jk) + zprofed(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert integr. bfe prod
+                zw2d(:,:) = zw2d(:,:) + ( zprofen(:,:,jk) + zprofed(:,:,jk) ) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk) ! vert integr. bfe prod
              ENDDO
             CALL iom_put( "INTPBFE" , zw2d )
@@ -496,5 +494,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * zysopt(:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert integr. bsi prod
+                zw2d(:,:) = zw2d(:,:) + zprorcad(:,:,jk) * zysopt(:,:,jk) * e3t_n(:,:,jk) * zfact * tmask(:,:,jk)  ! vert integr. bsi prod
              ENDDO
              CALL iom_put( "INTPBSI" , zw2d )
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 5845)
@@ -50,6 +50,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitnh4   !: -    -    -    -   -
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -103,9 +101,9 @@
             DO ji = 1, jpi
                zdep = MAX( hmld(ji,jj), heup(ji,jj) )
-               IF( fsdept(ji,jj,jk) < zdep ) THEN
+               IF( gdept_n(ji,jj,jk) < zdep ) THEN
                   zdepbac(ji,jj,jk) = MIN( 0.7 * ( trb(ji,jj,jk,jpzoo) + 2.* trb(ji,jj,jk,jpmes) ), 4.e-6 )
                   ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
                ELSE
-                  zdepmin = MIN( 1., zdep / fsdept(ji,jj,jk) )
+                  zdepmin = MIN( 1., zdep / gdept_n(ji,jj,jk) )
                   zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
                   zdepprod(ji,jj,jk) = zdepmin**0.273
@@ -283,5 +281,5 @@
                ! ----------------------------------------------------------
                zdep     = MAX( hmld(ji,jj), heup(ji,jj) ) 
-               zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
+               zdep     = MAX( 0., gdept_n(ji,jj,jk) - zdep )
                ztem     = MAX( tsn(ji,jj,1,jp_tem), 0. )
                zfactdep = xsilab * EXP(-( xsiremlab - xsirem ) * znusil2 * zdep / wsbio2 ) * ztem / ( ztem + 10. )
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5845)
@@ -25,5 +25,4 @@
    PUBLIC   p4z_sbc_init   
 
-   !! * Shared module variables
    LOGICAL , PUBLIC  :: ln_dust     !: boolean for dust input from the atmosphere
    LOGICAL , PUBLIC  :: ln_solub    !: boolean for variable solubility of atmospheric iron
@@ -45,5 +44,4 @@
    LOGICAL , PUBLIC  :: ll_sbc
 
-   !! * Module variables
    LOGICAL  ::  ll_solub
 
@@ -80,7 +78,5 @@
    REAL(wp), PUBLIC :: rivdininput, rivdipinput, rivdsiinput
 
-
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -89,5 +85,4 @@
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -163,5 +158,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * fse3t(ji,jj,1) + rtrn )
+                  nitdep(ji,jj) = sf_ndepo(1)%fnow(ji,jj,1) / rno3 / ( 14E6 * ryyss * e3t_n(ji,jj,1) + rtrn )
                END DO
             END DO
@@ -267,5 +262,5 @@
       IF( lk_offline ) THEN
         nk_rnf(:,:) = 1
-        h_rnf (:,:) = fsdept(:,:,1)
+        h_rnf (:,:) = gdept_n(:,:,1)
       ENDIF
 
@@ -456,5 +451,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zexpide   = MIN( 8.,( fsdept(ji,jj,jk) / 500. )**(-1.5) )
+                  zexpide   = MIN( 8.,( gdept_n(ji,jj,jk) / 500. )**(-1.5) )
                   zdenitide = -0.9543 + 0.7662 * LOG( zexpide ) - 0.235 * LOG( zexpide )**2
                   zcmask(ji,jj,jk) = zcmask(ji,jj,jk) * MIN( 1., EXP( zdenitide ) / 0.5 )
@@ -466,5 +461,5 @@
          ironsed(:,:,jpk) = 0._wp
          DO jk = 1, jpkm1
-            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( fse3t(:,:,jk) * rday )
+            ironsed(:,:,jk) = sedfeinput * zcmask(:,:,jk) / ( e3t_n(:,:,jk) * rday )
          END DO
          DEALLOCATE( zcmask)
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5845)
@@ -32,12 +32,8 @@
    PUBLIC   p4z_sed_alloc
  
-
-   !! * Module variables
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
    REAL(wp) :: r1_rday                  !: inverse of rday
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -100,5 +96,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zdep    = rfact2 / fse3t(ji,jj,1)
+               zdep    = rfact2 / e3t_n(ji,jj,1)
                zwflux  = fmmflx(ji,jj) / 1000._wp
                zfminus = MIN( 0._wp, -zwflux ) * trb(ji,jj,1,jpfer) * zdep
@@ -111,5 +107,5 @@
          ! 
          IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironice" ) )   &
-            &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
+            &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
          !
          CALL wrk_dealloc( jpi, jpj, zironice )
@@ -125,15 +121,15 @@
          !                                              ! Iron and Si deposition at the surface
          IF( ln_solub ) THEN
-            zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
+            zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
          ELSE
-            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
+            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
          ENDIF
-         zsidep(:,:) = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 28.1 
-         zpdep (:,:) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / fse3t(:,:,1) / 31. / po4r 
+         zsidep(:,:) = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
+         zpdep (:,:) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
          !                                              ! Iron solubilization of particles in the water column
          !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j
          zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 270. * rday )
          DO jk = 2, jpkm1
-            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -fsdept(:,:,jk) / 540. )
+            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / 540. )
          END DO
          !                                              ! Iron solubilization of particles in the water column
@@ -145,5 +141,5 @@
             IF( knt == nrdttrc ) THEN
                 IF( iom_use( "Irondep" ) )   &
-                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
+                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
                 IF( iom_use( "pdust" ) )   &
                 &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface
@@ -151,5 +147,5 @@
          ELSE                                    
             IF( ln_diatrc )  &
-              &  trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+              &  trc2d(:,:,jp_pcs0_2d + 11) = zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1)
          ENDIF
          CALL wrk_dealloc( jpi, jpj,      zpdep, zsidep )
@@ -206,5 +202,5 @@
          DO ji = 1, jpi
             ikt  = mbkt(ji,jj)
-            zdep = fse3t(ji,jj,ikt) / xstep
+            zdep = e3t_n(ji,jj,ikt) / xstep
             zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
             zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) )
@@ -230,5 +226,5 @@
               zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) )
               zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
-              zdep  = LOG10( fsdepw(ji,jj,ikt+1) )
+              zdep  = LOG10( gdepw_n(ji,jj,ikt+1) )
               zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
               &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
@@ -279,5 +275,5 @@
          DO ji = 1, jpi
             ikt  = mbkt(ji,jj)
-            zdep = xstep / fse3t(ji,jj,ikt) 
+            zdep = xstep / e3t_n(ji,jj,ikt) 
             zws4 = zwsbio4(ji,jj) * zdep
             zwsc = zwscal (ji,jj) * zdep
@@ -305,5 +301,5 @@
          DO ji = 1, jpi
             ikt  = mbkt(ji,jj)
-            zdep = xstep / fse3t(ji,jj,ikt) 
+            zdep = xstep / e3t_n(ji,jj,ikt) 
             zws4 = zwsbio4(ji,jj) * zdep
             zws3 = zwsbio3(ji,jj) * zdep
@@ -336,5 +332,5 @@
             tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) )
             tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt
-            sdenit(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt)
+            sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
 #endif
          END DO
@@ -388,5 +384,5 @@
                zwork1(:,:) = 0.
                DO jk = 1, jpkm1
-                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * zfact * fse3t(:,:,jk) * tmask(:,:,jk)
+                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
                ENDDO
                CALL iom_put( "INTNFIX" , zwork1 ) 
@@ -395,5 +391,5 @@
       ELSE
          IF( ln_diatrc )  &
-            &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+            &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1)
       ENDIF
       !
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 5845)
@@ -65,6 +65,4 @@
 #endif
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -108,5 +106,5 @@
             DO ji = 1,jpi
                zmax  = MAX( heup(ji,jj), hmld(ji,jj) )
-               zfact = MAX( 0., fsdepw(ji,jj,jk+1) - zmax ) / 5000._wp
+               zfact = MAX( 0., gdepw_n(ji,jj,jk+1) - zmax ) / 5000._wp
                wsbio4(ji,jj,jk) = wsbio2 + ( 200.- wsbio2 ) * zfact
             END DO
@@ -137,5 +135,5 @@
              DO ji = 1, jpi
                 IF( tmask(ji,jj,jk) == 1) THEN
-                   zwsmax =  0.5 * fse3t(ji,jj,jk) / xstep
+                   zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep
                    iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
                    iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
@@ -156,5 +154,5 @@
             DO ji = 1, jpi
                IF( tmask(ji,jj,jk) == 1 ) THEN
-                 zwsmax = 0.5 * fse3t(ji,jj,jk) / xstep
+                 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
                  wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
                  wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
@@ -700,5 +698,5 @@
          zl = zmin
          zr = zmax
-         wmax = 0.5 * fse3t(1,1,jk) * rday * float(niter1max) / rfact2
+         wmax = 0.5 * e3t_n(1,1,jk) * rday * float(niter1max) / rfact2
          zdiv = xkr_zeta + xkr_eta - xkr_eta * zl
          znum = zl - 1.
@@ -844,5 +842,5 @@
             DO jj = 1, jpj      
                DO ji = 1, jpi    
-                  zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
                   zew   = zwsink2(ji,jj,jk+1)
                   psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
@@ -858,5 +856,5 @@
             DO jj = 1,jpj
                DO ji = 1, jpi
-                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
                   trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
                END DO
@@ -869,5 +867,5 @@
          DO jj = 1,jpj
             DO ji = 1, jpi
-               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
+               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
                ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
             END DO
Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 5845)
@@ -21,7 +21,9 @@
    PUBLIC trc_wri_pisces 
 
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcnam.F90 5836 2015-10-26 14:49:40Z cetlod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -57,5 +59,5 @@
          zdic(:,:) = 0.
          DO jk = 1, jpkm1
-            zdic(:,:) = zdic(:,:) + trn(:,:,jk,jpdic) * fse3t(:,:,jk) * tmask(:,:,jk) * 12.
+            zdic(:,:) = zdic(:,:) + trn(:,:,jk,jpdic) * e3t_n(:,:,jk) * tmask(:,:,jk) * 12.
          ENDDO
          CALL iom_put( 'INTDIC', zdic )     
@@ -64,5 +66,5 @@
       IF( iom_use( "O2MIN" ) .OR. iom_use ( "ZO2MIN" ) ) THEN  ! Oxygen minimum concentration and depth 
          zo2min   (:,:) = trn(:,:,1,jpoxy) * tmask(:,:,1)
-         zdepo2min(:,:) = fsdepw(:,:,1)    * tmask(:,:,1)
+         zdepo2min(:,:) = gdepw_n(:,:,1)   * tmask(:,:,1)
          DO jk = 2, jpkm1
             DO jj = 1, jpj
@@ -71,5 +73,5 @@
                      IF( trn(ji,jj,jk,jpoxy) < zo2min(ji,jj) ) then
                         zo2min   (ji,jj) = trn(ji,jj,jk,jpoxy)
-                        zdepo2min(ji,jj) = fsdepw(ji,jj,jk)
+                        zdepo2min(ji,jj) = gdepw_n(ji,jj,jk)
                      ENDIF
                   ENDIF
