Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 6075)
@@ -34,5 +34,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6075)
@@ -165,5 +165,5 @@
          ! heat content variation
          zdiff_hc = zdiff_hc + glob_sum_full(  surf(:,:) * tmask(:,:,jk)                                   & 
-            &                           * e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) )
+            &                           * e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - surf_ini(:,:) * hc_loc_ini(:,:,jk) ) 
          ! salt content variation
          zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk)                           &
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIU/cool_skin.F90	(revision 6075)
@@ -40,5 +40,4 @@
 
       !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6075)
@@ -119,6 +119,6 @@
          CALL wrk_alloc( jpi,jpj,   z1_hu_0, z1_hv_0 )
          !
-         z1_hu_0(:,:) = umask_i(:,:) / ( hu_0(:,:) + 1._wp - umask_i(:,:) )     ! _i mask due to ISF
-         z1_hv_0(:,:) = vmask_i(:,:) / ( hv_0(:,:) + 1._wp - vmask_i(:,:) )
+         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
+         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
          !
          !        before       !          now          !       after         !
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6075)
@@ -537,5 +537,4 @@
       END DO
       !                                        ! Inverse of the local depth
-<<<<<<< .working
 !!gm BUG ?  don't understand the use of umask_i here .....
       r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 6075)
@@ -32,5 +32,4 @@
    PUBLIC   iscpl_cons       
    !! * Substitutions  
-#  include "domzgr_substitute.h90"  
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -95,5 +94,5 @@
       ! 
       zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-      IF ( lk_vvl ) zdssh = 0.0_wp ! already included in the levels by definition
+      IF (.NOT. ln_linssh ) zdssh = 0.0_wp ! already included in the levels by definition
       
       DO jk = 1,jpk-1
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 6075)
@@ -30,6 +30,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) ::   hdiv_iscpl
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl
-   !! * Substitutions  
-#  include "domzgr_substitute.h90"  
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 6075)
@@ -32,5 +32,4 @@
    !!
    !! * Substitutions  
-#  include "domzgr_substitute.h90"  
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -211,5 +210,6 @@
 
 !=============================================================================
-      IF ( lk_vvl ) THEN
+!PM: Is this needed since introduction of VVL by default?
+      IF (.NOT.ln_linssh) THEN
       ! Reconstruction of all vertical scale factors at now time steps
       ! =============================================================================
@@ -264,5 +264,5 @@
       ! t-, u- and v- water column thickness
       ! ------------------------------------
-         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
+         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
          DO jk = 1, jpkm1
             hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
@@ -386,5 +386,6 @@
 
       ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
-      IF ( lk_vvl ) THEN
+      !PM: Is this IF needed since change to VVL by default
+      IF (.NOT.ln_linssh) THEN
          DO jk = 2,jpk-1
             DO jj = 1,jpj
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 6075)
@@ -228,7 +228,7 @@
             DO ji = 1, jpi
                ikt=mikt(ji,jj)
-               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
                DO jk=2,ikt-1
-                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk) &
+                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
                      &                              * (1._wp - tmask(ji,jj,jk))
                END DO
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 6075)
@@ -795,5 +795,5 @@
                             &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
                             &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
-                            &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
+                            &               +    hu_n(ji,jj)  * zu_frc(ji,jj) ) &
                             &   ) * zhura
 
@@ -801,5 +801,5 @@
                             &     + rdtbt * ( zhvst_e(ji,jj)  *    zwy(ji,jj)   &
                             &               + zhvp2_e(ji,jj)  * zv_trd(ji,jj)   &
-                            &               +      hv(ji,jj)  * zv_frc(ji,jj) ) &
+                            &               +    hv_n(ji,jj)  * zv_frc(ji,jj) ) &
                             &   ) * zhvra
                END DO
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 6075)
@@ -79,6 +79,4 @@
    !!----------------------------------------------------------------------
 
-   !! * Substitutions 
-#  include "domzgr_substitute.h90"
 CONTAINS
 
@@ -529,5 +527,5 @@
       !! * Modules used
       USE dom_oce, ONLY : &             ! Ocean space and time domain variables
-         & gdept_n       
+         & gdept_n,       &      
          & gdept_1d      
       USE phycst, ONLY : &              ! Physical constants
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90	(revision 6075)
@@ -56,6 +56,4 @@
    !!----------------------------------------------------------------------
 
-   !! * Substitutions 
-#  include "domzgr_substitute.h90" 
 CONTAINS
 
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6075)
@@ -43,6 +43,4 @@
    !!----------------------------------------------------------------------
 
-!! * Substitutions
-#  include "domzgr_substitute.h90"  
 
 CONTAINS
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6075)
@@ -150,5 +150,5 @@
          DO jj = 1,jpj
             DO ji = 1,jpi
-               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj))
+               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
             END DO
          END DO
@@ -299,5 +299,5 @@
             DO jj = 1, jpj
                 ik = 2
-                DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
+                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
                 misfkt(ji,jj) = ik-1
             END DO
@@ -324,16 +324,16 @@
             ikb = misfkt(ji,jj)
             ! thickness of boundary layer at least the top level thickness
-            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
+            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
 
             ! determine the deepest level influenced by the boundary layer
             DO jk = ikt+1, mbkt(ji,jj)
-               IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
-            END DO
-            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
+               IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+            END DO
+            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
             misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
             r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
 
-            zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
-            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
+            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
          END DO
       END DO
@@ -608,5 +608,5 @@
                ELSE
                   !! compute Rc number (as done in zdfric.F90)
-                  zcoef = 0.5_wp / fse3w(ji,jj,ikt)
+                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
                   !                                            ! shear of horizontal velocity
                   zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
@@ -620,5 +620,5 @@
                   zts(jp_tem) = ttbl(ji,jj)
                   zts(jp_sal) = stbl(ji,jj)
-                  zdep        = fsdepw(ji,jj,ikt)
+                  zdep        = gdepw_n(ji,jj,ikt)
                   !
                   CALL eos_rab( zts, zdep, zab )
@@ -627,10 +627,10 @@
                   zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-      !! compute Monin Obukov Length
-               ! Maximum boundary layer depth
-               zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
-               ! Compute Monin obukhov length scale at the surface and Ekman depth:
-               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
-               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+                  !! compute Monin Obukov Length
+                  ! Maximum boundary layer depth
+                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
+                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
+                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
 
                   !! compute eta* (stability parameter)
@@ -686,20 +686,20 @@
                ikt = miku(ji,jj) ; ikb = miku(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3u_n(ji,jj,ikt))
 
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbku(ji,jj)
-                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = fse3u_n(ji,jj,jk)
+                  ze3 = e3u_n(ji,jj,jk)
                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
@@ -717,20 +717,20 @@
                ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
                ! thickness of boundary layer at least the top level thickness
-               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
 
                ! determine the deepest level influenced by the boundary layer
                DO jk = ikt+1, mbkv(ji,jj)
-                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
-               END DO
-               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
 
                ! level fully include in the ice shelf boundary layer
                DO jk = ikt, ikb - 1
-                  ze3 = fse3v_n(ji,jj,jk)
+                  ze3 = e3v_n(ji,jj,jk)
                   pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
                END DO
 
                ! level partially include in ice shelf boundary layer 
-               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
                pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
             END DO
Index: /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 6074)
+++ /branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 6075)
@@ -419,10 +419,10 @@
                ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
                !
-               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku)    ! i-direction: case 1
-               ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku)    ! -     -      case 2
-               ENDIF
-
-               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv)    ! j-direction: case 1
-               ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv)    ! -     -      case 2
+               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
+               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
+               ENDIF
+
+               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
+               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
                ENDIF
 
