Index: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90	(revision 9254)
+++ branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/domvvl.F90	(revision 9255)
@@ -65,5 +65,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO-Consortium (2015) 
-   !! $Id: domvvl.F90 9065 2017-12-14 17:00:50Z jchanut $
+   !! $Id: domvvl.F90 9190 2018-01-06 14:18:23Z gm $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -119,6 +119,4 @@
       REAL(wp)::   zcoef
       !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('dom_vvl_init')
       !
       IF(lwp) WRITE(numout,*)
@@ -240,6 +238,4 @@
       ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('dom_vvl_init')
-      !
    END SUBROUTINE dom_vvl_init
 
@@ -389,6 +385,5 @@
          !                       ! d - thickness diffusion transport: boundary conditions
          !                             (stored for tracer advction and continuity equation)
-         CALL lbc_lnk( un_td , 'U' , -1._wp)
-         CALL lbc_lnk( vn_td , 'V' , -1._wp)
+         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
 
          ! 4 - Time stepping of baroclinic scale factors
@@ -675,6 +670,4 @@
       !!----------------------------------------------------------------------
       !
-      IF( ln_timing )   CALL timing_start('dom_vvl_interpol')
-      !
       IF(ln_wd_il) THEN
         zlnwd = 1.0_wp
@@ -762,6 +755,4 @@
       END SELECT
       !
-      IF( ln_timing )   CALL timing_stop('dom_vvl_interpol')
-      !
    END SUBROUTINE dom_vvl_interpol
 
@@ -785,6 +776,4 @@
       INTEGER ::   id1, id2, id3, id4, id5     ! local integers
       !!----------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('dom_vvl_rst')
       !
       IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
@@ -914,16 +903,17 @@
             ELSE
                !
-               ! Just to read set ssh in fact, called latter once vertical grid
-               ! is set up:
-               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
+               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n
+               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )  
+               ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn)
                !
                DO jk=1,jpk
                   e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
-                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
+                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points
                END DO
                e3t_n(:,:,:) = e3t_b(:,:,:)
-!                sshn(:,:)=0._wp
-!                e3t_n(:,:,:)=e3t_0(:,:,:)
-!                e3t_b(:,:,:)=e3t_0(:,:,:)
+!!$                sshn(:,:)=0._wp
+!!$                e3t_n(:,:,:)=e3t_0(:,:,:)
+!!$                e3t_b(:,:,:)=e3t_0(:,:,:)
                !
             END IF           ! end of ll_wd edits
@@ -958,6 +948,4 @@
       ENDIF
       !
-      IF( ln_timing )   CALL timing_stop('dom_vvl_rst')
-      !
    END SUBROUTINE dom_vvl_rst
 
@@ -979,6 +967,5 @@
       REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist : 
       READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
-      !
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
       REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
@@ -990,31 +977,25 @@
          WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
-         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
-         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
-         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
-         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
+         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
+         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
+         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
+         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
+         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
          WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
-         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
-         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
-         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
-         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
-         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
-         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
+         WRITE(numout,*) '      !'
+         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
+         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
          IF( ln_vvl_ztilde_as_zstar ) THEN
-            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
-            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
-            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
-            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
-            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
-            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
+            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
+            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
+            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
+            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
+            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
+            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
          ELSE
-            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
-            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
-            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
-            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
+            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
+            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
          ENDIF
-         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
-         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
+         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
       ENDIF
       !
@@ -1030,15 +1011,12 @@
       IF(lwp) THEN                   ! Print the choice
          WRITE(numout,*)
-         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
-         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
-         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
-         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
-         ! - ML - Option not developed yet
-         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
-         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
+         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
+         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
+         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
+         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
       ENDIF
       !
 #if defined key_agrif
-      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
+      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
 #endif
       !
Index: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/usrdef_istate.F90
===================================================================
--- branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/usrdef_istate.F90	(revision 9254)
+++ branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/VORTEX/MY_SRC/usrdef_istate.F90	(revision 9255)
@@ -110,6 +110,6 @@
             DO jk=1, jpk
                zdu = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji+1,jj,jk))
-               zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
                IF (zdu < zH) THEN
+                  zf = (zH-1._wp-zdu+EXP(zdu-zH)) / (zH-1._wp+EXP(-zH))
                   pu(ji,jj,jk) = (za * zf * zy * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji+1,jj,jk)
                ELSE
@@ -126,6 +126,6 @@
             DO jk=1, jpk
                zdv = 0.5_wp * (pdept(ji  ,jj,jk) + pdept(ji,jj+1,jk))
-               zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
                IF (zdv < zH) THEN
+                  zf = (zH-1._wp-zdv+EXP(zdv-zH)) / (zH-1._wp+EXP(-zH))
                   pv(ji,jj,jk) = -(za * zf * zx * EXP(-(zx**2+zy**2)/zlambda**2)) * ptmask(ji,jj,jk) * ptmask(ji,jj+1,jk)
                ELSE
