Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4795)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 4990)
@@ -169,10 +169,29 @@
       fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
       fsdepw_b(:,:,1) = 0.0_wp
-      DO jk = 2, jpk
-         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
-         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
-         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
-         fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk)
-         fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1)
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            DO jk = 2,mikt(ji,jj)-1
+               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
+               fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk)
+               fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk)
+            END DO
+            IF (mikt(ji,jj) .GT. 1) THEN
+               jk = mikt(ji,jj)
+               fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+               fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk)
+               fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk)
+            END IF
+            DO jk = mikt(ji,jj)+1, jpk
+               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+               fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)
+               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1)
+            END DO
+         END DO
       END DO
 
@@ -185,6 +204,6 @@
          hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
-      hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
+      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
+      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
 
       ! Restoring frequencies for z_tilde coordinate
@@ -293,5 +312,5 @@
       !                                                ! --------------------------------------------- !
 
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
          ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
@@ -313,5 +332,5 @@
             zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
          END DO
-         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
+         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
 
          ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
@@ -338,5 +357,5 @@
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          END IF
@@ -386,6 +405,6 @@
          ! d - thickness diffusion transport: boundary conditions
          !     (stored for tracer advction and continuity equation)
-         CALL lbc_lnk( un_td , 'U' , -1.)
-         CALL lbc_lnk( vn_td , 'V' , -1.)
+         CALL lbc_lnk( un_td , 'U' , -1._wp)
+         CALL lbc_lnk( vn_td , 'V' , -1._wp)
 
          ! 4 - Time stepping of baroclinic scale factors
@@ -398,5 +417,5 @@
             z2dt = 2.0_wp * rdt
          ENDIF
-         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
+         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
          tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
 
@@ -453,5 +472,5 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
             dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
@@ -463,5 +482,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)
+            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -531,6 +550,6 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
+      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
 
       CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
@@ -570,5 +589,5 @@
       !! * Local declarations
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
-      INTEGER                             :: jk       ! dummy loop indices
+      INTEGER                             :: ji,jj,jk       ! dummy loop indices
       !!----------------------------------------------------------------------
 
@@ -622,8 +641,23 @@
       fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
-      DO jk = 2, jpk
-         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
-         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
-         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            DO jk = 2,mikt(ji,jj)-1
+               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
+            END DO
+            IF (mikt(ji,jj) .GT. 1) THEN
+               jk = mikt(ji,jj)
+               fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+            END IF
+            DO jk = mikt(ji,jj)+1, jpk
+               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
+               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+            END DO
+         END DO
       END DO
       ! Local depth and Inverse of the local depth of the water column at u- and v- points
@@ -702,5 +736,5 @@
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
          ! boundary conditions
-         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
+         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
          !               ! ------------------------------------- !
@@ -720,5 +754,5 @@
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
          ! boundary conditions
-         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
+         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
          !               ! ------------------------------------- !
@@ -738,5 +772,5 @@
          IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
          ! boundary conditions
-         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
+         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
          pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
          !               ! ------------------------------------- !
@@ -815,11 +849,25 @@
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
                CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
+               ! needed to restart if land processor not computed 
+               IF(lwp) write(numout,*) 'dom_vvl_rst : fse3t_b and fse3t_n found in restart files'
+               WHERE ( tmask(:,:,:) == 0.0_wp ) 
+                  fse3t_n(:,:,:) = e3t_0(:,:,:)
+                  fse3t_b(:,:,:) = e3t_0(:,:,:)
+               END WHERE
                IF( neuler == 0 ) THEN
                   fse3t_b(:,:,:) = fse3t_n(:,:,:)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.'
+               IF(lwp) write(numout,*) 'neuler is forced to 0'
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
+               fse3t_n(:,:,:) = fse3t_b(:,:,:)
+               neuler = 0
+            ELSE IF( id2 > 0 ) THEN
                IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
                IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
+               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
                fse3t_b(:,:,:) = fse3t_n(:,:,:)
                neuler = 0
@@ -830,5 +878,5 @@
                DO jk=1,jpk
                   fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
-                      &                            / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &
+                      &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &
                       &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
                END DO
@@ -964,4 +1012,5 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
+      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
 
       IF(lwp) THEN                   ! Print the choice
