Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5839)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 5845)
@@ -73,5 +73,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -181,6 +180,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,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) )
@@ -188,9 +187,9 @@
                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) )   &
+                  &                          * 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) )   &
+                  &                          * 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)
@@ -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
