Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90	(revision 6060)
@@ -28,5 +28,4 @@
  
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -72,11 +71,11 @@
       CASE( 'DYN' )                     ! T- and F-points
          DO jk = 1, jpk                      ! pah1 at T-point
-            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( fsdept(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
+            pah1(:,:,jk) = pahs1(:,:) * (  prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) )  ) * tmask(:,:,jk)
          END DO
          DO jk = 1, jpk                      ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zdep2 = (  fsdept(ji,jj+1,jk) + fsdept(ji+1,jj+1,jk)   &
-                     &     + fsdept(ji,jj  ,jk) + fsdept(ji+1,jj  ,jk)  ) * 0.25_wp
+                  zdep2 = (  gdept_n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk)   &
+                     &     + gdept_n(ji,jj  ,jk) + gdept_n(ji+1,jj  ,jk)  ) * 0.25_wp
                   pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * fmask(ji,jj,jk)
                END DO
@@ -89,6 +88,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  zdep1 = (  fsdept(ji,jj,jk) + fsdept(ji+1,jj,jk)  ) * 0.5_wp
-                  zdep2 = (  fsdept(ji,jj,jk) + fsdept(ji,jj+1,jk)  ) * 0.5_wp
+                  zdep1 = (  gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk)  ) * 0.5_wp
+                  zdep2 = (  gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk)  ) * 0.5_wp
                   pah1(ji,jj,jk) = pahs1(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  ) * umask(ji,jj,jk)
                   pah2(ji,jj,jk) = pahs2(ji,jj) * (  prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  ) * vmask(ji,jj,jk)
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90	(revision 6060)
@@ -52,5 +52,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 6060)
@@ -73,5 +73,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -181,16 +180,16 @@
                !                                      ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
                !                                      ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau )  )
-               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  )
+               zbu = MIN(  zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau )  )
+               zbv = MIN(  zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav )  )
                !                                      ! uslp and vslp output in zwz and zww, resp.
                zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
                zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
-               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              &
-                  &                   + zfi  * uslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   &
+               zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                                &
+                  &                   + zfi  * uslpml(ji,jj)                                                       &
+                  &                          * 0.5_wp * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk)-e3u_n(ji,jj,1) )   &
                   &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk)
-               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              &
-                  &                   + zfj  * vslpml(ji,jj)                                                     &
-                  &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   &
+               zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                                &
+                  &                   + zfj  * vslpml(ji,jj)                                                       &
+                  &                          * 0.5_wp * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk)-e3v_n(ji,jj,1) )   &
                   &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk)
 !!gm  modif to suppress omlmask.... (as in Griffies case)
@@ -198,6 +197,6 @@
 !               zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
 !               zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
-!               zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
-!               zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
+!               zci = 0.5 * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
+!               zcj = 0.5 * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
 !               zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
 !               zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
@@ -270,9 +269,9 @@
                !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
                !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai )  )
-               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  )
+               zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai )  )
+               zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj )  )
                !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
                zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0
-               zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
+               zck = gdepw_n(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )
                zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
                zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk)
@@ -281,5 +280,5 @@
 !               !                                         ! jk must be >= ML level for zfk=1. otherwise  zfk=0.
 !               zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
-!               zck = fsdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
+!               zck = gdepw(ji,jj,jk)    / MAX( hmlp(ji,jj), 10. )
 !               zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
 !               zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
@@ -441,5 +440,5 @@
                      zdks = 0._wp
                   ENDIF
-                  zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp)
+                  zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / e3w_n(ji,jj,jk+kp)
                   zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw )    ! force zdzrho >= repsln
                  END DO
@@ -451,5 +450,5 @@
          DO ji = 1, jpi
             jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1     ! MIN in case ML depth is the ocean depth
-            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)
+            z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk)
          END DO
       END DO
@@ -480,6 +479,6 @@
                      ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth)
                      zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
-                        &          - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk)
-                     ze3_e1    =  fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 
+                        &          - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj)  ) * umask(ji,jj,jk)
+                     ze3_e1    =  e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 
                      zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1  , ABS( zti_g_raw ) ), zti_g_raw )
                   ENDIF
@@ -490,6 +489,6 @@
                   ELSE
                      ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
-                        &      - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
-                     ze3_e2    =  fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj)
+                        &      - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
+                     ze3_e2    =  e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj)
                      ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2  , ABS( ztj_g_raw ) ), ztj_g_raw )
                   ENDIF
@@ -523,11 +522,11 @@
                      !
                      ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
-                     zti_coord = znot_thru_surface * ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj)
-                     ztj_coord = znot_thru_surface * ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj)     ! unmasked
+                     zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj  ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj)
+                     ztj_coord = znot_thru_surface * ( gdept_n(ji  ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj)     ! unmasked
                      zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces
                      ztj_g_raw = ztj_raw - ztj_coord
                      ! additional limit required in bilaplacian case
-                     ze3_e1    = fse3w(ji+ip,jj   ,jk+kp) * r1_e1u(ji,jj)
-                     ze3_e2    = fse3w(ji   ,jj+jp,jk+kp) * r1_e2v(ji,jj)
+                     ze3_e1    = e3w_n(ji+ip,jj   ,jk+kp) * r1_e1u(ji,jj)
+                     ze3_e2    = e3w_n(ji   ,jj+jp,jk+kp) * r1_e2v(ji,jj)
                      ! NB: hard coded factor 5 (can be a namelist parameter...)
                      zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw )
@@ -542,8 +541,8 @@
                      zti_g_lim =          ( zfacti   * zti_g_lim                       &
                         &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
-                        &                           * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
+                        &                           * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
                      ztj_g_lim =          ( zfactj   * ztj_g_lim                       &
                         &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
-                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
+                        &                           * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
                      !
                      triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim
@@ -577,8 +576,8 @@
                      triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw
                      !
-                     zbu  = e1e2u(ji   ,jj   ) * fse3u(ji   ,jj   ,jk   )
-                     zbv  = e1e2v(ji   ,jj   ) * fse3v(ji   ,jj   ,jk   )
-                     zbti = e1e2t(ji+ip,jj   ) * fse3w(ji+ip,jj   ,jk+kp)
-                     zbtj = e1e2t(ji   ,jj+jp) * fse3w(ji   ,jj+jp,jk+kp)
+                     zbu  = e1e2u(ji   ,jj   ) * e3u_n(ji   ,jj   ,jk   )
+                     zbv  = e1e2v(ji   ,jj   ) * e3v_n(ji   ,jj   ,jk   )
+                     zbti = e1e2t(ji+ip,jj   ) * e3w_n(ji+ip,jj   ,jk+kp)
+                     zbtj = e1e2t(ji   ,jj+jp) * e3w_n(ji   ,jj+jp,jk+kp)
                      !
                      wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim      ! masked
@@ -682,6 +681,6 @@
             !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
             !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau )  )
-            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav )  )
+            zbu = MIN(  zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau )  )
+            zbv = MIN(  zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav )  )
             !                        !- Slope at u- & v-points (uslpml, vslpml)
             uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
@@ -705,6 +704,6 @@
             !                        !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
             !                           kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
-            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai )  )
-            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  )
+            zbi = MIN(  zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai )  )
+            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj )  )
             !                        !- i- & j-slope at w-points (wslpiml, wslpjml)
             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
@@ -767,5 +766,5 @@
 
          !!gm I no longer understand this.....
-!!gm         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
+!!gm         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (.NOT.ln_linssh .AND. ln_rstart) ) THEN
 !            IF(lwp)   WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
 !
@@ -775,12 +774,12 @@
 !
 !            ! set the slope of diffusion to the slope of s-surfaces
-!            !      ( c a u t i o n : minus sign as fsdep has positive value )
+!            !      ( c a u t i o n : minus sign as dep has positive value )
 !            DO jk = 1, jpk
 !               DO jj = 2, jpjm1
 !                  DO ji = fs_2, fs_jpim1   ! vector opt.
-!                     uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
-!                     vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
-!                     wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
-!                     wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                     uslp (ji,jj,jk) = - ( gdept_n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+!                     vslp (ji,jj,jk) = - ( gdept_n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+!                     wslpi(ji,jj,jk) = - ( gdepw_n(ji+1,jj,jk) - gdepw_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
+!                     wslpj(ji,jj,jk) = - ( gdepw_n(ji,jj+1,jk) - gdepw_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
 !                  END DO
 !               END DO
Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 6060)
@@ -81,5 +81,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -515,9 +514,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
                   zhw(ji,jj) = zhw(ji,jj) + ze3w
@@ -533,9 +532,9 @@
                   ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
                   zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
-                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk)
                   ! Compute elements required for the inverse time scale of baroclinic
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
@@ -711,10 +710,10 @@
       !
       DO jk = 1, jpkm1                                         ! e2u e3u u_eiv = -dk[psi_uw]
-         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * fse3u(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) )
       END DO
       CALL iom_put( "uoce_eiv", zw3d )
       !
       DO jk = 1, jpkm1                                         ! e1v e3v v_eiv = -dk[psi_vw]
-         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * fse3v(:,:,jk) )
+         zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) )
       END DO
       CALL iom_put( "voce_eiv", zw3d )
