Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 6006)
@@ -498,8 +498,12 @@
 Contrary to the horizontal grid, the vertical grid is computed in the code and no 
 provision is made for reading it from a file. The only input file is the bathymetry 
-(in meters) (\ifile{bathy\_meter}) 
+(in meters) (\ifile{bathy\_meter}). 
 \footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the 
 \ifile{bathy\_meter} file, so that the computation of the number of wet ocean point 
 in each water column is by-passed}. 
+If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft 
+(in meters) (\ifile{isf\_draft\_meter}) is needed and all the location where the isf cavity thinnest
+ than \np{rn\_isfhmin} meters are grounded (ie masked). 
+
 After reading the bathymetry, the algorithm for vertical grid definition differs 
 between the different options:
Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 6006)
@@ -996,4 +996,39 @@
  at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div} 
 (called from \mdl{divcur}). 
+
+\section{ Ice sheet coupling}
+\label{SBC_iscpl}
+%------------------------------------------namsbc_iscpl----------------------------------------------------
+\namdisplay{namsbc_iscpl}
+%--------------------------------------------------------------------------------------------------------
+Ice sheet/ocean coupling is done through file exchange at the restart step. NEMO, at each restart step, 
+read the bathymetry and ice shelf draft variable in a netcdf file. 
+If \np{ln\_iscpl = ~true}, the isf draft is assume to be different at each restart step 
+with potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
+The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases:
+\begin{description}
+\item[Thin a cell down:]
+	T/S/ssh are unchanged and U/V in the top cell are corrected to keep the barotropic transport (bt) constant ($bt_b=bt_n$).
+\item[Enlarge  a cell:]
+	See case "Thin a cell down"
+\item[Dry a cell:]
+	mask, T/S, U/V and ssh are set to 0. Furthermore, U/V into the water column are modified to satisfy ($bt_b=bt_n$).
+\item[Wet a cell:] 
+	mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$ and U/V set to 0. If no neighbours along i,j and k, T/S/U/V and mask are set to 0.
+\item[Dry a column:]
+	mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0.
+\item[Wet a column:]
+	set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0. If no neighbour, T/S/U/V and mask set to 0.
+\end{description}
+The extrapolation is call \np{nn\_drown} times. It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps,
+ the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\
+This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe. 
+This could lead to a trend in heat/salt content and volume. In order to remove the trend and keep the conservation level as close to 0 as possible,
+ a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location. 
+Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added,
+ over a period of \np{rn\_fiscpl} time step, into the system. 
+So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\
+
+As the before and now fields are not compatible (modification of the geometry), the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$.
 %
 % ================================================================
Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namdom
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namdom	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namdom	(revision 6006)
@@ -6,4 +6,5 @@
    nn_msh      =    0      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
+   rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice
    rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
    rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namrun
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namrun	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namrun	(revision 6006)
@@ -15,4 +15,5 @@
    cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input)
    cn_ocerst_out = "restart"   !  suffix of ocean restart name (output)
+   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model (ln_iscpl = T => namsbc_iscpl)
    nn_istate   =       0   !  output the initial state (1) or not (0)
    nn_stock    =    5475   !  frequency of creation of a restart file (modulo referenced to 1)
Index: /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namsbc_iscpl
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namsbc_iscpl	(revision 6006)
+++ /branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namsbc_iscpl	(revision 6006)
@@ -0,0 +1,7 @@
+!-----------------------------------------------------------------------
+&namsbc_iscpl  !   land ice / ocean coupling option
+!-----------------------------------------------------------------------
+   nn_drown  = 10       ! number of iterations of the extrapolation loop (fill the new wet cells)
+   ln_hsb    = .false.  ! activate conservation module (conservation exact after a time of rn_fiscpl)
+   nn_fiscpl = 43800    ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency)
+/
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 6006)
@@ -200,10 +200,16 @@
 
          <!-- * variable related to ice shelf forcing * -->
-         <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
-         <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     />
-         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature)"   unit="m/s"      />
-         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)"      unit="m/s"      />
-         <field id="stbl"         long_name="salinity in the Losh tbl"                      unit="1e-3"     />
-         <field id="ttbl"         long_name="temperature in the Losh tbl"                   unit="degC"     />
+         <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  />
+         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     />
+         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "                 unit="m/s"      />
+         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "                 unit="m/s"      />
+ 	 <field id="stbl"         long_name="salinity in the Losh tbl                    "                 unit="PSU"      />
+	 <field id="ttbl"         long_name="temperature in the Losh tbl                 "                 unit="C"        />
+	 <field id="utbl"         long_name="zonal current in the Losh tbl at T point    "                 unit="m/s"      />
+	 <field id="vtbl"         long_name="merid current in the Losh tbl at T point    "                 unit="m/s"      />
+	 <field id="thermald"     long_name="thermal driving of ice shelf melting        "                 unit="C"        />
+	 <field id="tfrz"         long_name="top freezing point (used to compute melt)   "                 unit="C"        />
+	 <field id="tinsitu"      long_name="top insitu temperature (used to cmpt melt)  "                 unit="C"        />
+	 <field id="ustar"        long_name="ustar at T point used in ice shelf melting  "                 unit="m/s"      />
 
          <!-- *_oce variables available with ln_blk_clio or ln_blk_core -->
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6006)
@@ -43,4 +43,5 @@
    cn_ocerst_out = "restart"   !  suffix of ocean restart name (output)
    cn_ocerst_outdir = "."      !  directory in which to write output ocean restarts
+   ln_iscpl    = .false.   !  cavity evolution forcing or coupling to ice sheet model
    nn_istate   =       0   !  output the initial state (1) or not (0)
    ln_rst_list = .false.   !  output restarts at list of times using nn_stocklist (T) or at set frequency with nn_stock (F)
@@ -128,4 +129,5 @@
    nn_msh      =    1      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
+   rn_isfhmin  =    1.00   !  treshold (m) to discriminate grounding ice to floating ice
    rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
    rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
@@ -212,4 +214,5 @@
 !!   namsbc_rnf      river runoffs
 !!   namsbc_isf      ice shelf melting/freezing
+!!   namsbc_iscpl    coupling option between land ice model and ocean
 !!   namsbc_apr      Atmospheric Pressure
 !!   namsbc_ssr      sea surface restoring term (for T and/or S)
@@ -456,4 +459,11 @@
 /
 !-----------------------------------------------------------------------
+&namsbc_iscpl  !   land ice / ocean coupling option                     
+!-----------------------------------------------------------------------
+   nn_drown  = 10       ! number of iteration of the extrapolation loop (fill the new wet cells)
+   ln_hsb    = .false.  ! activate conservation module (conservation exact after a time of rn_fiscpl)
+   nn_fiscpl = 43800    ! (number of time step) conservation period (maybe should be fix to the coupling frequencey of restart frequency)
+/
+!-----------------------------------------------------------------------
 &namsbc_apr    !   Atmospheric pressure used as ocean forcing or in bulk
 !-----------------------------------------------------------------------
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 6006)
@@ -116,6 +116,6 @@
          &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   &
          &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   &
-         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
-      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin,   &
+         &             nn_write, ln_iscpl, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler
+      NAMELIST/namdom/ nn_bathy , rn_bathy, rn_e3zps_min, rn_e3zps_rat, nn_msh    , rn_hmin, rn_isfhmin,  &
          &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,            &
          &             rn_rdtmax, rn_rdth     , nn_baro     , nn_closea , ln_crs, &
@@ -813,15 +813,15 @@
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector loop
-            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
          END DO
          DO ji = 1, jpim1      ! NO vector opt.
-            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
+            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
                &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
          END DO
       END DO
-      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask_i, 'V', 1._wp )
-      CALL lbc_lnk( fmask_i, 'F', 1._wp )
+      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( ssvmask, 'V', 1._wp )
+      CALL lbc_lnk( ssfmask, 'F', 1._wp )
 
       ! 3. Ocean/land mask at wu-, wv- and w points 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 6006)
@@ -460,5 +460,5 @@
       ENDIF
 
-      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
+      IF( nn_timing == 1 )   CALL timing_stop('dia_fwb')
 
  9005 FORMAT(1X,A,ES24.16)
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90	(revision 6006)
@@ -135,5 +135,5 @@
       DO jk=1,nb_ana
        DO ji=1,jpmax_harmo
-          IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
+          IF (TRIM(tname(jk)) ==  Wave(ji)%cname_tide) THEN
              name(jk) = ji
              EXIT
@@ -195,6 +195,6 @@
                      ! Elevation
                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)        
-                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj)
-                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj)
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
                   END DO
                END DO
@@ -324,6 +324,6 @@
                X1= ana_amp(ji,jj,jh,1)
                X2=-ana_amp(ji,jj,jh,2)
-               out_u(ji,jj,       jh) = X1 * umask_i(ji,jj)
-               out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
+               out_u(ji,jj,       jh) = X1 * ssumask(ji,jj)
+               out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj)
             ENDDO
          ENDDO
@@ -488,5 +488,5 @@
             DO jj_sd = ji_sd, ninco
                zval2 = ABS(ztmp3(ji_sd,jj_sd))
-               IF( zval2.GE.zval1 )THEN
+               IF( zval2 >= zval1 )THEN
                   ipivot(ji_sd) = jj_sd
                   zval1         = zval2
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 6006)
@@ -46,5 +46,6 @@
    REAL(wp) ::   frc_wn_t, frc_wn_s    ! global forcing trends
    !
-   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf          , ssh_ini          !
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf 
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   surf_ini      , ssh_ini          !
    REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   ssh_hc_loc_ini, ssh_sc_loc_ini   !
    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   hc_loc_ini, sc_loc_ini, e3t_ini  !
@@ -138,4 +139,5 @@
       ! 2 -  Content variations !
       ! ------------------------ !
+      ! glob_sum_full is needed because you keep the full interior domain to compute the sum (iscpl)
       zdiff_v2 = 0._wp
       zdiff_hc = 0._wp
@@ -143,5 +145,5 @@
 
       ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum( surf(:,:) * ( sshn(:,:) - ssh_ini(:,:) ) )
+      zdiff_v1 = glob_sum_full( surf(:,:) * sshn(:,:) - surf_ini(:,:) * ssh_ini(:,:) )
 
       ! heat & salt content variation (associated with ssh)
@@ -158,18 +160,18 @@
             z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
          END IF
-         z_ssh_hc = glob_sum( z2d0 ) 
-         z_ssh_sc = glob_sum( z2d1 ) 
+         z_ssh_hc = glob_sum_full( z2d0 ) 
+         z_ssh_sc = glob_sum_full( z2d1 ) 
       ENDIF
 
       DO jk = 1, jpkm1
          ! volume variation (calculated with scale factors)
-         zdiff_v2 = zdiff_v2 + glob_sum( surf(:,:) * tmask(:,:,jk) &
-            &                           * ( fse3t_n(:,:,jk) - e3t_ini(:,:,jk) ) )
-         ! heat content variation
-         zdiff_hc = zdiff_hc + glob_sum(  surf(:,:) * tmask(:,:,jk) & 
-            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) - hc_loc_ini(:,:,jk) ) )
+         zdiff_v2 = zdiff_v2 + glob_sum_full( surf    (:,:) * fse3t_n  (:,:,jk) * tmask(:,:,jk) &
+              &                             - surf_ini(:,:) *   e3t_ini(:,:,jk) )
+         ! heat content variation 
+         zdiff_hc = zdiff_hc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) &
+              &                             - surf_ini(:,:) * hc_loc_ini(:,:,jk) )
          ! salt content variation
-         zdiff_sc = zdiff_sc + glob_sum(  surf(:,:) * tmask(:,:,jk)   &
-            &                           * ( fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) - sc_loc_ini(:,:,jk) ) )
+         zdiff_sc = zdiff_sc + glob_sum_full( surf    (:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) * tsn(:,:,jk,jp_sal) & 
+              &                             - surf_ini(:,:) * sc_loc_ini(:,:,jk) )
       ENDDO
 
@@ -191,5 +193,5 @@
       zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors)
       DO jk = 1, jpkm1
-         zvol_tot  = zvol_tot + glob_sum( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
+         zvol_tot  = zvol_tot + glob_sum_full( surf(:,:) * tmask(:,:,jk) * fse3t_n(:,:,jk) )
       END DO
 
@@ -214,6 +216,6 @@
         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu)
         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J) 
-        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3)
-        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)  
+        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9   )              ! Salt content variation (psu*km3)
+        CALL iom_put( 'bgvolssh' , zdiff_v1  * 1.e-9   )              ! volume ssh variation (km3)  
         CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3) 
         CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C) 
@@ -261,4 +263,5 @@
               CALL iom_get( numror, 'frc_wn_s', frc_wn_s )
            ENDIF
+           CALL iom_get( numror, jpdom_autoglo, 'surf_ini', surf_ini ) ! ice sheet coupling
            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini )
            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini )
@@ -273,9 +276,11 @@
           IF(lwp) WRITE(numout,*) ' dia_hsb at initial state '
           IF(lwp) WRITE(numout,*) '~~~~~~~'
-          ssh_ini(:,:) = sshn(:,:)                                       ! initial ssh
+          surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
+          ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
           DO jk = 1, jpk
-             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)                        ! initial vertical scale factors
-             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk)   ! initial heat content
-             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk)   ! initial salt content
+             ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
+             e3t_ini   (:,:,jk) = fse3t_n(:,:,jk)    * tmask(:,:,jk)                    ! initial vertical scale factors
+             hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
+             sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * fse3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
           END DO
           frc_v = 0._wp                                           ! volume       trend due to forcing
@@ -312,4 +317,5 @@
            CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s )
         ENDIF
+        CALL iom_rstput( kt, nitrst, numrow, 'surf_ini', surf_ini )      ! ice sheet coupling
         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini )
         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini )
@@ -379,6 +385,6 @@
       ! 1 - Allocate memory !
       ! ------------------- !
-      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), &
-         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror )
+      ALLOCATE( hc_loc_ini(jpi,jpj,jpk), sc_loc_ini(jpi,jpj,jpk), surf_ini(jpi,jpj), &
+         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror  )
       IF( ierror > 0 ) THEN
          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6006)
@@ -32,4 +32,5 @@
    REAL(wp), PUBLIC ::   rn_bathy        !: depth of flat bottom (active if nn_bathy=0; if =0 depth=jpkm1)
    REAL(wp), PUBLIC ::   rn_hmin         !: minimum ocean depth (>0) or minimum number of ocean levels (<0)
+   REAL(wp), PUBLIC ::   rn_isfhmin      !: threshold to discriminate grounded ice to floating ice
    REAL(wp), PUBLIC ::   rn_e3zps_min    !: miminum thickness for partial steps (meters)
    REAL(wp), PUBLIC ::   rn_e3zps_rat    !: minimum thickness ration for partial steps
@@ -43,4 +44,5 @@
    INTEGER , PUBLIC ::   nn_closea       !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
    INTEGER , PUBLIC ::   nn_euler        !: =0 start with forward time step or not (=1)
+   LOGICAL , PUBLIC ::   ln_iscpl       !: coupling with ice sheet
    LOGICAL , PUBLIC ::   ln_crs          !: Apply grid coarsening to dynamical model output or online passive tracers
 
@@ -251,13 +253,14 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function
 
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF)
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                       (ISF)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask                   !: surface domain T-point mask 
-
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
@@ -386,13 +389,14 @@
          &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) )
 
-      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      &
-         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), &
-         &     bmask  (jpi,jpj) ,                                                       &
-         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )
+      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       &
+         &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       & 
+         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
+         &     bmask(jpi,jpj)   ,                                                          &
+         &     mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
 
 ! (ISF) Allocation of basic array   
-      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),                   &
-         &      mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           &
-         &      mikf(jpi,jpj), ssmask(jpi,jpj), STAT=ierr(10) )
+      ALLOCATE( misfdep(jpi,jpj) , risfdep(jpi,jpj),     &
+         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           &
+         &     mikf(jpi,jpj), STAT=ierr(10) )
 
       ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     & 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6006)
@@ -112,9 +112,12 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
 
                              CALL dom_stp      ! time step
-      IF( nmsh /= 0      )   CALL dom_wri      ! Create a domain file
+      !
+      IF( nmsh /= 0 .AND. .NOT. ln_iscpl )                         CALL dom_wri      ! Create a domain file
+      IF( nmsh /= 0 .AND.       ln_iscpl .AND. .NOT. ln_rstart )   CALL dom_wri      ! Create a domain file
+      !
       IF( .NOT.ln_rstart )   CALL dom_ctl      ! Domain control
       !
@@ -135,15 +138,15 @@
       !!----------------------------------------------------------------------
       USE ioipsl
-      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               &
-                       nn_no   , cn_exp  , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,  &
-         &             nn_it000, nn_itend, nn_date0    , nn_time0     , nn_leapy  , nn_istate ,  &
-         &             nn_stock, nn_write, ln_dimgnnn  , ln_mskland   , ln_clobber, nn_chunksz,  &
-         &             nn_euler, ln_cfmeta
-      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
-         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
-         &             rn_rdtmax, rn_rdth     , nn_closea , ln_crs,    &
-         &             jphgr_msh, &
-         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m, &
-         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, &
+      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
+                       nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
+         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
+         &             nn_stock, nn_write , ln_dimgnnn  , ln_mskland   , ln_clobber, nn_chunksz,     &
+         &             nn_euler, ln_cfmeta, ln_iscpl
+      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, &
+         &             nn_acc   , rn_atfp , rn_rdt      , rn_rdtmin   ,                              &
+         &             rn_rdtmax, rn_rdth , nn_closea   , ln_crs      ,                              &
+         &             jphgr_msh,                                                                    &
+         &             ppglam0, ppgphi0, ppe1_deg, ppe2_deg, ppe1_m, ppe2_m,                         &
+         &             ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh,                  &
          &             ppa2, ppkth2, ppacr2
 #if defined key_netcdf4
@@ -193,4 +196,5 @@
          WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber
          WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz
+         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl   = ', ln_iscpl
       ENDIF
 
@@ -260,4 +264,5 @@
          WRITE(numout,*) '      min depth of the ocean    (>0) or    rn_hmin   = ', rn_hmin
          WRITE(numout,*) '      min number of ocean level (<0)       '
+         WRITE(numout,*) '      treshold to open the isf cavity   rn_isfhmin   = ', rn_isfhmin, ' (m)'
          WRITE(numout,*) '      minimum thickness of partial      rn_e3zps_min = ', rn_e3zps_min, ' (m)'
          WRITE(numout,*) '         step level                     rn_e3zps_rat = ', rn_e3zps_rat
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6006)
@@ -183,4 +183,6 @@
       ! --------------------
       tmask_i(:,:) = ssmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf
+
+      tmask_h(:,:) = 1._wp                 ! 0 on the halo and 1 elsewhere
       iif = jpreci                         ! ???
       iil = nlci - jpreci + 1
@@ -188,8 +190,8 @@
       ijl = nlcj - jprecj + 1
 
-      tmask_i( 1 :iif,   :   ) = 0._wp      ! first columns
-      tmask_i(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
-      tmask_i(   :   , 1 :ijf) = 0._wp      ! first rows
-      tmask_i(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
+      tmask_h( 1 :iif,   :   ) = 0._wp      ! first columns
+      tmask_h(iil:jpi,   :   ) = 0._wp      ! last  columns (including mpp extra columns)
+      tmask_h(   :   , 1 :ijf) = 0._wp      ! first rows
+      tmask_h(   :   ,ijl:jpj) = 0._wp      ! last  rows (including mpp extra rows)
 
       ! north fold mask
@@ -202,8 +204,11 @@
          IF( mjg(nlej) == jpjglo ) THEN                  ! only half of the nlcj-1 row
             DO ji = iif+1, iil-1
-               tmask_i(ji,nlej-1) = tmask_i(ji,nlej-1) * tpol(mig(ji))
+               tmask_h(ji,nlej-1) = tmask_h(ji,nlej-1) * tpol(mig(ji))
             END DO
          ENDIF
       ENDIF
+     
+      tmask_i(:,:) = tmask_i(:,:) * tmask_h(:,:)
+
       IF( jperio == 5 .OR. jperio == 6 ) THEN      ! F-point pivot
          tpol(     1    :jpiglo) = 0._wp
@@ -225,21 +230,21 @@
          END DO
       END DO
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
       DO jj = 1, jpjm1
          DO ji = 1, fs_jpim1   ! vector loop
-            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
          END DO
          DO ji = 1, jpim1      ! NO vector opt.
-            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
+            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
                &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
          END DO
       END DO
-      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask, 'V', 1._wp )
-      CALL lbc_lnk( fmask, 'F', 1._wp )
-      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask_i, 'V', 1._wp )
-      CALL lbc_lnk( fmask_i, 'F', 1._wp )
+      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask  , 'V', 1._wp )
+      CALL lbc_lnk( fmask  , 'F', 1._wp )
+      CALL lbc_lnk( ssumask, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( ssvmask, 'V', 1._wp )
+      CALL lbc_lnk( ssfmask, 'F', 1._wp )
 
       ! 3. Ocean/land mask at wu-, wv- and w points 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90	(revision 6006)
@@ -28,5 +28,5 @@
 CONTAINS
 
-   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid )
+   SUBROUTINE dom_ngb( plon, plat, kii, kjj, cdgrid, kkk )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dom_ngb  ***
@@ -39,6 +39,8 @@
       REAL(wp)        , INTENT(in   ) ::   plon, plat   ! longitude,latitude of the point
       INTEGER         , INTENT(  out) ::   kii, kjj     ! i-,j-index of the closes grid point
+      INTEGER         , INTENT(in   ), OPTIONAL :: kkk  ! k-index of the mask level used
       CHARACTER(len=1), INTENT(in   ) ::   cdgrid       ! grid name 'T', 'U', 'V', 'W'
       !
+      INTEGER :: ik         ! working level
       INTEGER , DIMENSION(2) ::   iloc
       REAL(wp)               ::   zlon, zmini
@@ -51,17 +53,23 @@
       !
       zmask(:,:) = 0._wp
+      ik = 1
+      IF ( PRESENT(kkk) ) ik=kkk
       SELECT CASE( cdgrid )
-      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,1)
-      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,1)
-      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,1)
-      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,1)
+      CASE( 'U' )  ; zglam(:,:) = glamu(:,:) ; zgphi(:,:) = gphiu(:,:) ; zmask(nldi:nlei,nldj:nlej) = umask(nldi:nlei,nldj:nlej,ik)
+      CASE( 'V' )  ; zglam(:,:) = glamv(:,:) ; zgphi(:,:) = gphiv(:,:) ; zmask(nldi:nlei,nldj:nlej) = vmask(nldi:nlei,nldj:nlej,ik)
+      CASE( 'F' )  ; zglam(:,:) = glamf(:,:) ; zgphi(:,:) = gphif(:,:) ; zmask(nldi:nlei,nldj:nlej) = fmask(nldi:nlei,nldj:nlej,ik)
+      CASE DEFAULT ; zglam(:,:) = glamt(:,:) ; zgphi(:,:) = gphit(:,:) ; zmask(nldi:nlei,nldj:nlej) = tmask(nldi:nlei,nldj:nlej,ik)
       END SELECT
 
-      zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
-      zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
-      IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
-      IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
+      IF (jphgr_msh /= 2 .AND. jphgr_msh /= 3) THEN
+         zlon       = MOD( plon       + 720., 360. )                                     ! plon between    0 and 360
+         zglam(:,:) = MOD( zglam(:,:) + 720., 360. )                                     ! glam between    0 and 360
+         IF( zlon > 270. )   zlon = zlon - 360.                                          ! zlon between  -90 and 270
+         IF( zlon <  90. )   WHERE( zglam(:,:) > 180. ) zglam(:,:) = zglam(:,:) - 360.   ! glam between -180 and 180
+         zglam(:,:) = zglam(:,:) - zlon
+      ELSE
+         zglam(:,:) = zglam(:,:) - plon
+      END IF
 
-      zglam(:,:) = zglam(:,:) - zlon
       zgphi(:,:) = zgphi(:,:) - plat
       zdist(:,:) = zglam(:,:) * zglam(:,:) + zgphi(:,:) * zgphi(:,:)
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 6006)
@@ -24,5 +24,5 @@
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O manager library
-   USE restart         ! ocean restart
+   USE restart, ONLY : rst_read_open    ! ocean restart
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -192,6 +192,6 @@
          hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
       END DO
-      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )
-      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )
+      hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1. - ssumask(:,:) )
+      hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1. - ssvmask(:,:) )
 
       ! Restoring frequencies for z_tilde coordinate
@@ -325,5 +325,5 @@
          ! --------------------------------------------------
          IF( ln_vvl_ztilde ) THEN
-            IF( kt .GT. nit000 ) THEN
+            IF( kt > nit000 ) THEN
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
@@ -419,5 +419,5 @@
          IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
          ! - ML - test: for the moment, stop simulation for too large e3_t variations
-         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
+         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
             IF( lk_mpp ) THEN
                CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
@@ -538,6 +538,6 @@
       END DO
       !                                        ! Inverse of the local depth
-      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
 
       CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6006)
@@ -379,5 +379,5 @@
       !!              - bathy : meter bathymetry (in meters)
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
+      INTEGER  ::   ji, jj, jk            ! dummy loop indices
       INTEGER  ::   inum                      ! temporary logical unit
       INTEGER  ::   ierror                    ! error flag
@@ -541,4 +541,11 @@
                CALL iom_close( inum )
                WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp
+
+               ! set grounded point to 0 
+               ! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
+               WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin )
+                  misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+                  mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
+               END WHERE
             END IF
             !       
@@ -578,11 +585,4 @@
       !                       
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
-         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin
-         IF ( ln_isfcav ) THEN
-            WHERE (bathy == risfdep)
-               bathy   = 0.0_wp ; risfdep = 0.0_wp
-            END WHERE
-         END IF
-         ! end patch
          IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
          ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
@@ -835,7 +835,7 @@
    END SUBROUTINE zgr_bot_level
 
-      SUBROUTINE zgr_top_level
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bot_level  ***
+   SUBROUTINE zgr_top_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_top_level  ***
       !!
       !! ** Purpose :   defines the vertical index of ocean top (mik. arrays)
@@ -963,7 +963,6 @@
       REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
       REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
-      REAL(wp) ::   zmax             ! Maximum depth
       REAL(wp) ::   zdiff            ! temporary scalar
-      REAL(wp) ::   zrefdep          ! temporary scalar
+      REAL(wp) ::   zmax             ! temporary scalar
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt
       !!---------------------------------------------------------------------
@@ -1004,6 +1003,4 @@
       END DO
 
-      IF ( ln_isfcav ) CALL zgr_isf
-
       ! Scale factors and depth at T- and W-points
       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
@@ -1013,5 +1010,750 @@
          e3w_0  (:,:,jk) = e3w_1d  (jk)
       END DO
+      
+      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
+      IF ( ln_isfcav ) CALL zgr_isf
+
+      ! Scale factors and depth at T- and W-points
+      IF ( .NOT. ln_isfcav ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ik = mbathy(ji,jj)
+               IF( ik > 0 ) THEN               ! ocean point only
+                  ! max ocean level case
+                  IF( ik == jpkm1 ) THEN
+                     zdepwp = bathy(ji,jj)
+                     ze3tp  = bathy(ji,jj) - gdepw_1d(ik)
+                     ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
+                     e3t_0(ji,jj,ik  ) = ze3tp
+                     e3t_0(ji,jj,ik+1) = ze3tp
+                     e3w_0(ji,jj,ik  ) = ze3wp
+                     e3w_0(ji,jj,ik+1) = ze3tp
+                     gdepw_0(ji,jj,ik+1) = zdepwp
+                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
+                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
+                     !
+                  ELSE                         ! standard case
+                     IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
+                     ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
+                     ENDIF
+   !gm Bug?  check the gdepw_1d
+                     !       ... on ik
+                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
+                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
+                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
+                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
+                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
+                     e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
+                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                     !       ... on ik+1
+                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
+                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+         !
+         it = 0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ik = mbathy(ji,jj)
+               IF( ik > 0 ) THEN               ! ocean point only
+                  e3tp (ji,jj) = e3t_0(ji,jj,ik)
+                  e3wp (ji,jj) = e3w_0(ji,jj,ik)
+                  ! test
+                  zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  )
+                  IF( zdiff <= 0._wp .AND. lwp ) THEN 
+                     it = it + 1
+                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
+                     WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+                     WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
+                     WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  )
+                  ENDIF
+               ENDIF
+            END DO
+         END DO
+      END IF
+      !
+      ! Scale factors and depth at U-, V-, UW and VW-points
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3u_0 (:,:,jk) = e3t_1d(jk)
+         e3v_0 (:,:,jk) = e3t_1d(jk)
+         e3uw_0(:,:,jk) = e3w_1d(jk)
+         e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+
+      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
+               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
+               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
+               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      IF ( ln_isfcav ) THEN
+      ! (ISF) define e3uw (adapted for 2 cells in the water column)
+         DO jj = 2, jpjm1 
+            DO ji = 2, fs_jpim1   ! vector opt. 
+               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
+               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
+               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
+               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
+               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
+                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
+            END DO
+         END DO
+      END IF
+
+      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
+      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
+      !
+
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
+         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
+         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
+      END DO
+      
+      ! Scale factor at F-point
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
+      END DO
+!!gm  bug ? :  must be a do loop with mj0,mj1
       ! 
+      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
+      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
+      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
+      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
+      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
+
+      ! Control of the sign
+      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
+      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
+      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
+      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
+     
+      ! Compute gdep3w_0 (vertical sum of e3w)
+      IF ( ln_isfcav ) THEN ! if cavity
+         WHERE (misfdep == 0) misfdep = 1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
+               DO jk = 2, misfdep(ji,jj)
+                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+               IF (misfdep(ji,jj) >= 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
+               DO jk = misfdep(ji,jj) + 1, jpk
+                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+      ELSE ! no cavity
+         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
+         DO jk = 2, jpk
+            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
+         END DO
+      END IF
+      !                                               ! ================= !
+      IF(lwp .AND. ll_print) THEN                     !   Control print   !
+         !                                            ! ================= !
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ik = MAX( mbathy(ji,jj), 1 )
+               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
+               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
+               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
+               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
+               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
+               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+      ENDIF  
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
+      !
+   END SUBROUTINE zgr_zps
+
+   SUBROUTINE zgr_isf
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_isf  ***
+      !!   
+      !! ** Purpose :   check the bathymetry in levels
+      !!   
+      !! ** Method  :   THe water column have to contained at least 2 cells
+      !!                Bathymetry and isfdraft are modified (dig/close) to respect
+      !!                this criterion.
+      !!   
+      !! ** Action  : - test compatibility between isfdraft and bathy 
+      !!              - bathy and isfdraft are modified
+      !!----------------------------------------------------------------------
+      !!   
+      INTEGER  ::   ji, jj, jl, jk       ! dummy loop indices
+      INTEGER  ::   ik, it               ! temporary integers
+      INTEGER  ::   icompt, ibtest       ! (ISF)
+      INTEGER  ::   ibtestim1, ibtestip1 ! (ISF)
+      INTEGER  ::   ibtestjm1, ibtestjp1 ! (ISF)
+      REAL(wp) ::   zdepth           ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zmax             ! Maximum and minimum depth
+      REAL(wp) ::   zbathydiff       ! isf temporary scalar
+      REAL(wp) ::   zrisfdepdiff     ! isf temporary scalar
+      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
+      REAL(wp) ::   zdepwp           ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zdiff            ! temporary scalar
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
+      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
+      !
+      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
+      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
+
+      ! (ISF) compute misfdep
+      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
+      ELSEWHERE                      ;                         misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
+      END WHERE  
+
+      ! Compute misfdep for ocean points (i.e. first wet level) 
+      ! find the first ocean level such that the first level thickness 
+      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 
+      ! e3t_0 is the reference level thickness 
+      DO jk = 2, jpkm1 
+         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
+         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
+      END DO 
+      WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )
+         risfdep(:,:) = 0. ; misfdep(:,:) = 1
+      END WHERE
+
+! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
+      icompt = 0 
+! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
+      DO jl = 1, 10     
+         ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
+         WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin)
+            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
+            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
+         END WHERE
+         WHERE (mbathy(:,:) <= 0) 
+            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
+            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
+         ENDWHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
+            misfdep( 1 ,:) = misfdep(jpim1,:)            ! local domain is cyclic east-west 
+            misfdep(jpi,:) = misfdep(  2  ,:) 
+            mbathy( 1 ,:)  = mbathy(jpim1,:)             ! local domain is cyclic east-west
+            mbathy(jpi,:)  = mbathy(  2  ,:)
+         ENDIF
+
+         ! split last cell if possible (only where water column is 2 cell or less)
+         ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
+         IF ( .NOT. ln_iscpl) THEN
+            DO jk = jpkm1, 1, -1
+               zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+               WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
+                  mbathy(:,:) = jk
+                  bathy(:,:)  = zmax
+               END WHERE
+            END DO
+         END IF
+ 
+         ! split top cell if possible (only where water column is 2 cell or less)
+         DO jk = 2, jpkm1
+            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
+            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
+               misfdep(:,:) = jk
+               risfdep(:,:) = zmax
+            END WHERE
+         END DO
+
+ 
+ ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! find the minimum change option:
+               ! test bathy
+               IF (risfdep(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN
+                     zbathydiff  =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
+                         &            + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
+                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+ 
+                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
+                        IF (zbathydiff <= zrisfdepdiff) THEN
+                           bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
+                           mbathy(ji,jj)= mbathy(ji,jj) + 1
+                        ELSE
+                           risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
+                           misfdep(ji,jj) = misfdep(ji,jj) - 1
+                        END IF
+                     ENDIF
+                  ELSE
+                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN
+                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
+                        misfdep(ji,jj) = misfdep(ji,jj) - 1
+                     END IF
+                  END IF
+               END IF
+            END DO
+         END DO
+ 
+         ! At least 2 levels for water thickness at T, U, and V point.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! find the minimum change option:
+               ! test bathy
+               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj)   - ( gdepw_1d(mbathy (ji,jj)+1) &
+                         &                             + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj)  ) & 
+                         &                             - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)= misfdep(ji,jj) - 1
+                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)= misfdep(ji,jj) - 1
+                     risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+ ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(bathy(ji,jj    ) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji,jj) = mbathy(ji,jj) + 1
+                        bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
+                        risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
+                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
+                           &                             + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji,jj  ) - ( gdepw_1d(misfdep(ji,jj  )  ) &
+                           &                             - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
+                        bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
+                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
+                  END IF
+               ENDIF
+            END DO
+         END DO
+ 
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF 
+ 
+ ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                  zbathydiff  =ABS(  bathy(ji  ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
+                       &                              + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
+                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
+                       &                              - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
+                  IF (zbathydiff <= zrisfdepdiff) THEN
+                     mbathy(ji,jj) = mbathy(ji,jj) + 1
+                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
+                  ELSE
+                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
+                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
+                  END IF
+                  ELSE
+                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
+                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF 
+ 
+ ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN
+                  IF ( .NOT. ln_iscpl ) THEN 
+                     zbathydiff  =ABS(  bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
+                          &                              + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
+                     zrisfdepdiff=ABS(risfdep(ji  ,jj) - ( gdepw_1d(misfdep(ji  ,jj)  ) &
+                          &                              - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
+                     IF (zbathydiff <= zrisfdepdiff) THEN
+                        mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
+                        bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
+                     ELSE
+                        misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
+                        risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
+                     END IF
+                  ELSE
+                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
+                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
+                  ENDIF
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+      END DO
+      ! end dig bathy/ice shelf to be compatible
+      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
+      DO jl = 1,20
+ 
+ ! remove single point "bay" on isf coast line in the ice shelf draft'
+         DO jk = 2, jpk
+            WHERE (misfdep==0) misfdep=jpk
+            zmask=0._wp
+            WHERE (misfdep <= jk) zmask=1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF (misfdep(ji,jj) == jk) THEN
+                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
+                     IF (ibtest <= 1) THEN
+                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         WHERE (misfdep==jpk)
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
+         END WHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ 
+ ! remove single point "bay" on bathy coast line beneath an ice shelf'
+         DO jk = jpk,1,-1
+            zmask=0._wp
+            WHERE (mbathy >= jk ) zmask=1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN
+                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
+                     IF (ibtest <= 1) THEN
+                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
+                        IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         WHERE (mbathy==0)
+             misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
+         END WHERE
+         IF( lk_mpp ) THEN
+            zbathy(:,:)  = FLOAT( misfdep(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            misfdep(:,:) = INT( zbathy(:,:) )
+
+            CALL lbc_lnk( risfdep,'T', 1. )
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:)  = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:)  = INT( zbathy(:,:) )
+         ENDIF
+ 
+ ! fill hole in ice shelf
+         zmisfdep = misfdep
+         zrisfdep = risfdep
+         WHERE (zmisfdep <= 1._wp) zmisfdep=jpk
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
+               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
+               IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj  ) ) ibtestim1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj  ) ) ibtestip1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj-1) ) ibtestjm1 = jpk
+               IF( zmisfdep(ji,jj) >= mbathy(ji  ,jj+1) ) ibtestjp1 = jpk
+               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
+               IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
+               END IF
+               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN
+                  misfdep(ji,jj) = ibtest
+                  risfdep(ji,jj) = gdepw_1d(ibtest)
+               ENDIF
+            ENDDO
+         ENDDO
+ 
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ !
+ !! fill hole in bathymetry
+         zmbathy (:,:)=mbathy (:,:)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
+               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
+               IF( zmbathy(ji,jj) <  misfdep(ji-1,jj  ) ) ibtestim1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji+1,jj  ) ) ibtestip1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj-1) ) ibtestjm1 = 0
+               IF( zmbathy(ji,jj) <  misfdep(ji  ,jj+1) ) ibtestjp1 = 0
+               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
+               IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
+               END IF
+               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN
+                  mbathy(ji,jj) = ibtest
+                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
+               ENDIF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie U point water column less than 2 cells), mask U
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
+                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy,  'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep, 'T', 1. ) 
+            CALL lbc_lnk( bathy,   'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy,  'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie U point water column less than 2 cells), mask U
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
+                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie V point water column less than 2 cells), mask V
+         DO jj = 1, jpjm1
+            DO ji = 1, jpi
+               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
+                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check (ie V point water column less than 2 cells), mask V
+         DO jj = 1, jpjm1
+            DO ji = 1, jpi
+               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
+                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
+               END IF
+            END DO
+         END DO
+         IF( lk_mpp ) THEN 
+            zbathy(:,:)  = FLOAT( misfdep(:,:) ) 
+            CALL lbc_lnk( zbathy, 'T', 1. ) 
+            misfdep(:,:) = INT( zbathy(:,:) ) 
+
+            CALL lbc_lnk( risfdep,'T', 1. ) 
+            CALL lbc_lnk( bathy,  'T', 1. )
+
+            zbathy(:,:) = FLOAT( mbathy(:,:) )
+            CALL lbc_lnk( zbathy, 'T', 1. )
+            mbathy(:,:) = INT( zbathy(:,:) )
+         ENDIF 
+ ! if not compatible after all check, mask T
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
+                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
+               END IF
+            END DO
+         END DO
+ 
+         WHERE (mbathy(:,:) == 1)
+            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
+         END WHERE
+      END DO 
+! end check compatibility ice shelf/bathy
+      ! remove very shallow ice shelf (less than ~ 10m if 75L)
+      IF( icompt == 0 ) THEN 
+         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
+      ELSE 
+         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
+      ENDIF 
+
+      ! compute scale factor and depth at T- and W- points
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1035,4 +1777,5 @@
                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
                   ENDIF
+      !            gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
 !gm Bug?  check the gdepw_1d
                   !       ... on ik
@@ -1040,12 +1783,9 @@
                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
-                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
-                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
-                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   &
-                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
+                  e3t_0  (ji,jj,ik  ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik  )
+                  e3w_0  (ji,jj,ik  ) = gdept_0(ji,jj,ik  ) - gdept_1d(ik-1)
                   !       ... on ik+1
                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
-                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
                ENDIF
             ENDIF
@@ -1073,680 +1813,50 @@
       END DO
       !
-      IF ( ln_isfcav ) THEN
       ! (ISF) Definition of e3t, u, v, w for ISF case
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               ik = misfdep(ji,jj) 
-               IF( ik > 1 ) THEN               ! ice shelf point only 
-                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
-                  gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            ik = misfdep(ji,jj) 
+            IF( ik > 1 ) THEN               ! ice shelf point only 
+               IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik) 
+               gdepw_0(ji,jj,ik) = risfdep(ji,jj) 
 !gm Bug?  check the gdepw_0 
-               !       ... on ik 
-                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
-                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
-                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
-                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
-                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
-
-                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
-                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
-                  ENDIF 
-               !       ... on ik / ik-1 
-                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 
-                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
+            !       ... on ik 
+               gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   & 
+                  &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   & 
+                  &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      ) 
+               e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 
+               e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
+
+               IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column) 
+                  e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 
+               ENDIF 
+            !       ... on ik / ik-1 
+               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) 
+               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 
-                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+               gdept_0(ji,jj,ik-1) = gdept_1d(ik-1)
+            ENDIF 
+         END DO 
+      END DO 
+   
+      it = 0 
+      DO jj = 1, jpj 
+         DO ji = 1, jpi 
+            ik = misfdep(ji,jj) 
+            IF( ik > 1 ) THEN               ! ice shelf point only 
+               e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
+               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
+            ! test 
+               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
+               IF( zdiff <= 0. .AND. lwp ) THEN  
+                  it = it + 1 
+                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
+                  WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
+                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
+                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
                ENDIF 
-            END DO 
+            ENDIF 
          END DO 
-      ! 
-         it = 0 
-         DO jj = 1, jpj 
-            DO ji = 1, jpi 
-               ik = misfdep(ji,jj) 
-               IF( ik > 1 ) THEN               ! ice shelf point only 
-                  e3tp (ji,jj) = e3t_0(ji,jj,ik  ) 
-                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 
-               ! test 
-                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  ) 
-                  IF( zdiff <= 0. .AND. lwp ) THEN  
-                     it = it + 1 
-                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
-                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 
-                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
-                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj) 
-                  ENDIF 
-               ENDIF 
-            END DO 
-         END DO 
-      END IF
-      ! END (ISF)
-
-      ! Scale factors and depth at U-, V-, UW and VW-points
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3u_0 (:,:,jk) = e3t_1d(jk)
-         e3v_0 (:,:,jk) = e3t_1d(jk)
-         e3uw_0(:,:,jk) = e3w_1d(jk)
-         e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
-               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
-               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
-               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
-            END DO
-         END DO
-      END DO
-      IF ( ln_isfcav ) THEN
-      ! (ISF) define e3uw (adapted for 2 cells in the water column)
-         DO jj = 2, jpjm1 
-            DO ji = 2, fs_jpim1   ! vector opt. 
-               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
-               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
-               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
-                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
-               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
-               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
-               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
-                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
-            END DO
-         END DO
-      END IF
-
-      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions
-      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp )
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk)
-         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk)
-         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk)
-      END DO
-      
-      ! Scale factor at F-point
-      DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions
-      !
-      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk)
-      END DO
-!!gm  bug ? :  must be a do loop with mj0,mj1
-      ! 
-      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
-      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 
-      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 
-      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 
-      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 
-
-      ! Control of the sign
-      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' )
-      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' )
-      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' )
-      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' )
-     
-      ! Compute gdep3w_0 (vertical sum of e3w)
-      IF ( ln_isfcav ) THEN ! if cavity
-         WHERE (misfdep == 0) misfdep = 1
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
-               DO jk = 2, misfdep(ji,jj)
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-               END DO
-               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
-               DO jk = misfdep(ji,jj) + 1, jpk
-                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 
-               END DO
-            END DO
-         END DO
-      ELSE ! no cavity
-         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
-         DO jk = 2, jpk
-            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
-         END DO
-      END IF
-      !                                               ! ================= !
-      IF(lwp .AND. ll_print) THEN                     !   Control print   !
-         !                                            ! ================= !
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX( mbathy(ji,jj), 1 )
-               zprt(ji,jj,1) = e3t_0   (ji,jj,ik)
-               zprt(ji,jj,2) = e3w_0   (ji,jj,ik)
-               zprt(ji,jj,3) = e3u_0   (ji,jj,ik)
-               zprt(ji,jj,4) = e3v_0   (ji,jj,ik)
-               zprt(ji,jj,5) = e3f_0   (ji,jj,ik)
-               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-      ENDIF  
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zprt )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps')
-      !
-   END SUBROUTINE zgr_zps
-
-   SUBROUTINE zgr_isf
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_isf  ***
-      !!   
-      !! ** Purpose :   check the bathymetry in levels
-      !!   
-      !! ** Method  :   THe water column have to contained at least 2 cells
-      !!                Bathymetry and isfdraft are modified (dig/close) to respect
-      !!                this criterion.
-      !!                 
-      !!   
-      !! ** Action  : - test compatibility between isfdraft and bathy 
-      !!              - bathy and isfdraft are modified
-      !!----------------------------------------------------------------------
-      !!   
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  ::   ik, it           ! temporary integers
-      INTEGER  ::   id, jd, nprocd
-      INTEGER  ::   icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1   ! (ISF)
-      LOGICAL  ::   ll_print         ! Allow  control print for debugging
-      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
-      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
-      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
-      REAL(wp) ::   zdiff            ! temporary scalar
-      REAL(wp) ::   zrefdep          ! temporary scalar
-      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH)
-      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH)
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_isf')
-      !
-      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep)
-      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy )
-
-
-      ! (ISF) compute misfdep
-      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1  
-      ELSEWHERE                      ;                          misfdep(:,:) = 2   ! iceshelf : initialize misfdep to second level 
-      END WHERE  
-
-      ! Compute misfdep for ocean points (i.e. first wet level) 
-      ! find the first ocean level such that the first level thickness 
-      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 
-      ! e3t_0 is the reference level thickness 
-      DO jk = 2, jpkm1 
-         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
-         WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth )   misfdep(:,:) = jk+1 
       END DO 
-      WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT. 0._wp)
-         risfdep(:,:) = 0. ; misfdep(:,:) = 1
-      END WHERE
- 
-! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
-      icompt = 0 
-! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
-      DO jl = 1, 10     
-         WHERE (bathy(:,:) .EQ. risfdep(:,:) )
-            misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
-            mbathy (:,:) = 0 ; bathy  (:,:) = 0._wp
-         END WHERE
-         WHERE (mbathy(:,:) <= 0) 
-            misfdep(:,:) = 0; risfdep(:,:) = 0._wp 
-            mbathy (:,:) = 0; bathy  (:,:) = 0._wp
-         ENDWHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
-         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
-            misfdep( 1 ,:) = misfdep(jpim1,:)           ! local domain is cyclic east-west 
-            misfdep(jpi,:) = misfdep(  2  ,:) 
-         ENDIF
-
-         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
-            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west
-            mbathy(jpi,:) = mbathy(  2  ,:)
-         ENDIF
-
-         ! split last cell if possible (only where water column is 2 cell or less)
-         DO jk = jpkm1, 1, -1
-            zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-            WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
-               mbathy(:,:) = jk
-               bathy(:,:)  = zmax
-            END WHERE
-         END DO
- 
-         ! split top cell if possible (only where water column is 2 cell or less)
-         DO jk = 2, jpkm1
-            zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
-            WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
-               misfdep(:,:) = jk
-               risfdep(:,:) = zmax
-            END WHERE
-         END DO
-
- 
- ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! find the minimum change option:
-               ! test bathy
-               IF (risfdep(ji,jj) .GT. 1) THEN
-               zbathydiff =ABS(bathy(ji,jj)   - (gdepw_1d(mbathy (ji,jj)+1) &
-                 &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-               zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) &
-                 &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
- 
-                  IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT. misfdep(ji,jj)) THEN
-                     IF (zbathydiff .LE. zrisfdepdiff) THEN
-                        bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
-                        mbathy(ji,jj)= mbathy(ji,jj) + 1
-                     ELSE
-                        risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
-                        misfdep(ji,jj) = misfdep(ji,jj) - 1
-                     END IF
-                  END IF
-               END IF
-            END DO
-         END DO
- 
-          ! At least 2 levels for water thickness at T, U, and V point.
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               ! find the minimum change option:
-               ! test bathy
-               IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff  =ABS(bathy(ji,jj)  - (gdepw_1d(mbathy (ji,jj)+1)&
-                    & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  )  &
-                    & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)= misfdep(ji,jj) - 1
-                     risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
- ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(bathy(ji,jj    ) - (gdepw_1d(mbathy (ji,jj)+1)   &
-                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj  )+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj+1) - (gdepw_1d(misfdep(ji,jj+1))  &
-                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj  )) &
-                   &    + MIN( e3zps_min, e3t_1d(mbathy(ji,jj   )+1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj+1)  = misfdep(ji,jj+1) - 1
-                     risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) &
-                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- ! point V misdep(ji,jj) EQ mbathy(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji,jj+1) - (gdepw_1d(mbathy (ji,jj+1)+1) &
-                   &   + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji,jj  ) - (gdepw_1d(misfdep(ji,jj  )  )  &
-                   &   - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
-                     bathy  (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1)  ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)   = misfdep(ji,jj) - 1
-                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji,jj  )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj  )  )*e3zps_rat )
-                  END IF
-               ENDIF
-            END DO
-         END DO
- 
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- 
- ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji  ,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
-                    &   + MIN( e3zps_min, e3t_1d(mbathy (ji  ,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji+1,jj) - (gdepw_1d(misfdep(ji+1,jj)) &
-                    &  - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji,jj) = mbathy(ji,jj) + 1
-                     bathy(ji,jj)  = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
-                     risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
-                  END IF
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- 
- ! point U misfdep(ji,jj) EQ bathy(ji,jj+1) 
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GT. 1) THEN
-                  zbathydiff =ABS(  bathy(ji+1,jj) - (gdepw_1d(mbathy (ji+1,jj)+1) &
-                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
-                  zrisfdepdiff=ABS(risfdep(ji  ,jj) - (gdepw_1d(misfdep(ji  ,jj)  )  &
-                      &  - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)-1)*e3zps_rat )))
-                  IF (zbathydiff .LE. zrisfdepdiff) THEN
-                     mbathy(ji+1,jj)  = mbathy (ji+1,jj) + 1
-                     bathy (ji+1,jj)  = gdepw_1d(mbathy (ji+1,jj)  )  &
-                      &   + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
-                  ELSE
-                     misfdep(ji,jj)   = misfdep(ji  ,jj) - 1
-                     risfdep(ji,jj)   = gdepw_1d(misfdep(ji  ,jj)+1) &
-                      &   - MIN( e3zps_min, e3t_1d(misfdep(ji  ,jj)   )*e3zps_rat )
-                  END IF
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
-      END DO
-      ! end dig bathy/ice shelf to be compatible
-      ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
-      DO jl = 1,20
- 
- ! remove single point "bay" on isf coast line in the ice shelf draft'
-         DO jk = 2, jpk
-            WHERE (misfdep==0) misfdep=jpk
-            zmask=0
-            WHERE (misfdep .LE. jk) zmask=1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (misfdep(ji,jj) .EQ. jk) THEN
-                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
-                        risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) misfdep(ji,jj) = jpk
-                     END IF
-                  END IF
-               END DO
-            END DO
-         END DO
-         WHERE (misfdep==jpk)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
-         END WHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- 
- ! remove single point "bay" on bathy coast line beneath an ice shelf'
-         DO jk = jpk,1,-1
-            zmask=0
-            WHERE (mbathy .GE. jk ) zmask=1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (mbathy(ji,jj) .EQ. jk .AND. misfdep(ji,jj) .GE. 2) THEN
-                     ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
-                     IF (ibtest .LE. 1) THEN
-                        bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
-                        IF (misfdep(ji,jj) .GT. mbathy(ji,jj)) mbathy(ji,jj) = 0
-                     END IF
-                  END IF
-               END DO
-            END DO
-         END DO
-         WHERE (mbathy==0)
-             misfdep=0 ; risfdep=0. ; mbathy=0 ; bathy=0.
-         END WHERE
-         IF( lk_mpp ) THEN
-            zbathy(:,:) = FLOAT( misfdep(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            misfdep(:,:) = INT( zbathy(:,:) )
-            CALL lbc_lnk( risfdep, 'T', 1. )
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF
- 
- ! fill hole in ice shelf
-         zmisfdep = misfdep
-         zrisfdep = risfdep
-         WHERE (zmisfdep .LE. 1) zmisfdep=jpk
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ibtestim1 = zmisfdep(ji-1,jj  ) ; ibtestip1 = zmisfdep(ji+1,jj  )
-               ibtestjm1 = zmisfdep(ji  ,jj-1) ; ibtestjp1 = zmisfdep(ji  ,jj+1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji-1,jj  ) ) ibtestim1 = jpk!MAX(0, mbathy(ji-1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji+1,jj  ) ) ibtestip1 = jpk!MAX(0, mbathy(ji+1,jj  ) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj-1) ) ibtestjm1 = jpk!MAX(0, mbathy(ji  ,jj-1) - 1)
-               IF( zmisfdep(ji,jj) .GE. mbathy(ji  ,jj+1) ) ibtestjp1 = jpk!MAX(0, mbathy(ji  ,jj+1) - 1)
-               ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == jpk .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
-               END IF
-               IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) .GE. 2) THEN
-                  misfdep(ji,jj) = ibtest
-                  risfdep(ji,jj) = gdepw_1d(ibtest)
-               ENDIF
-            ENDDO
-         ENDDO
- 
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- !
- !! fill hole in bathymetry
-         zmbathy (:,:)=mbathy (:,:)
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ibtestim1 = zmbathy(ji-1,jj  ) ; ibtestip1 = zmbathy(ji+1,jj  )
-               ibtestjm1 = zmbathy(ji  ,jj-1) ; ibtestjp1 = zmbathy(ji  ,jj+1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji-1,jj  ) ) ibtestim1 = 0!MIN(jpk-1, misfdep(ji-1,jj  ) + 1)
-               IF( zmbathy(ji,jj) .LT. misfdep(ji+1,jj  ) ) ibtestip1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj-1) ) ibtestjm1 = 0
-               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0
-               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
-               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
-               END IF
-               IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) .GE. 2) THEN
-                  mbathy(ji,jj) = ibtest
-                  bathy(ji,jj)  = gdepw_1d(ibtest+1) 
-               ENDIF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie U point water column less than 2 cells), mask U
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
-                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie U point water column less than 2 cells), mask U
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-               IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji+1,jj) .GE. 1) THEN
-                  mbathy(ji+1,jj)  = mbathy(ji+1,jj) - 1;   bathy(ji+1,jj)   = gdepw_1d(mbathy(ji+1,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie V point water column less than 2 cells), mask V
-         DO jj = 1, jpjm1
-            DO ji = 1, jpi
-               IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
-                  mbathy(ji,jj)  = mbathy(ji,jj) - 1 ; bathy(ji,jj)   = gdepw_1d(mbathy(ji,jj)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check (ie V point water column less than 2 cells), mask V
-         DO jj = 1, jpjm1
-            DO ji = 1, jpi
-               IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) .GE. 1 .AND. mbathy(ji,jj+1) .GE. 1) THEN
-                  mbathy(ji,jj+1)  = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1)   = gdepw_1d(mbathy(ji,jj+1)+1) ;
-               END IF
-            END DO
-         END DO
-         IF( lk_mpp ) THEN 
-            zbathy(:,:) = FLOAT( misfdep(:,:) ) 
-            CALL lbc_lnk( zbathy, 'T', 1. ) 
-            misfdep(:,:) = INT( zbathy(:,:) ) 
-            CALL lbc_lnk( risfdep, 'T', 1. ) 
-            CALL lbc_lnk( bathy, 'T', 1. )
-            zbathy(:,:) = FLOAT( mbathy(:,:) )
-            CALL lbc_lnk( zbathy, 'T', 1. )
-            mbathy(:,:) = INT( zbathy(:,:) )
-         ENDIF 
- ! if not compatible after all check, mask T
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
-                  misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj)  = 0 ; bathy(ji,jj)   = 0._wp ;
-               END IF
-            END DO
-         END DO
- 
-         WHERE (mbathy(:,:) == 1)
-            mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
-         END WHERE
-      END DO 
-! end check compatibility ice shelf/bathy
-      ! remove very shallow ice shelf (less than ~ 10m if 75L)
-      WHERE (misfdep(:,:) <= 5)
-         misfdep = 1; risfdep = 0.0_wp;
-      END WHERE
-
-      IF( icompt == 0 ) THEN 
-         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry' 
-      ELSE 
-         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 
-      ENDIF 
 
       CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep )
@@ -1755,5 +1865,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('zgr_isf')
       
-   END SUBROUTINE
+   END SUBROUTINE zgr_isf
 
    SUBROUTINE zgr_sco
@@ -2084,11 +2194,11 @@
       fsde3w(:,:,:) = gdep3w_0(:,:,:)
       !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      where (e3t_0   (:,:,:) == 0.0)  e3t_0(:,:,:)  = 1.0_wp
+      where (e3u_0   (:,:,:) == 0.0)  e3u_0(:,:,:)  = 1.0_wp
+      where (e3v_0   (:,:,:) == 0.0)  e3v_0(:,:,:)  = 1.0_wp
+      where (e3f_0   (:,:,:) == 0.0)  e3f_0(:,:,:)  = 1.0_wp
+      where (e3w_0   (:,:,:) == 0.0)  e3w_0(:,:,:)  = 1.0_wp
+      where (e3uw_0  (:,:,:) == 0.0)  e3uw_0(:,:,:) = 1.0_wp
+      where (e3vw_0  (:,:,:) == 0.0)  e3vw_0(:,:,:) = 1.0_wp
 
 #if defined key_agrif
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 6006)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 6006)
@@ -0,0 +1,322 @@
+MODULE iscplhsb
+   !!======================================================================
+   !!                       ***  MODULE  iscplhsb***
+   !! Ocean forcing: ice sheet/ocean coupling (conservation)
+   !!=====================================================================
+   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   iscpl_alloc    : variable allocation
+   !!   iscpl_hsb      : compute and store the input of heat/salt/volume 
+   !!                    into the system due to the coupling process
+   !!   iscpl_div      : correction of divergence to keep volume conservation
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE domwri          ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! surface boundary condition variables
+   USE oce             ! global tra/dyn variable
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! MPP library
+   USE wrk_nemo        ! Memory allocation
+   USE lbclnk          !
+   USE domngb          !
+   USE iscplini
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   iscpl_div   
+   PUBLIC   iscpl_cons       
+   !! * Substitutions  
+#  include "domzgr_substitute.h90"  
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE iscpl_cons(ptmask_b, psmask_b, pe3t_b, pts_flx, pvol_flx, prdt_iscpl)
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_cons  ***
+      !! 
+      !! ** Purpose :   compute input into the system during the coupling step
+      !!                compute the correction term
+      !!                compute where the correction have to be applied
+      !! 
+      !! ** Method  :   compute tsn*e3t-tsb*e3tb and e3t-e3t_b
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b    !! mask before
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b      !! scale factor before
+      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b    !! mask before
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pts_flx     !! corrective flux to have tracer conservation
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(out) :: pvol_flx    !! corrective flux to have volume conservation
+      REAL(wp),                     INTENT(in ) :: prdt_iscpl  !! coupling period 
+      !!
+      INTEGER :: ji, jj, jk                                    !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1
+      !!
+      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
+      REAL(wp):: r1_rdtiscpl
+      REAL(wp):: zjip1_ratio  , zjim1_ratio  , zjjp1_ratio  , zjjm1_ratio
+      !!
+      REAL(wp):: zde3t, zdtem, zdsal
+      REAL(wp), DIMENSION(:,:), POINTER :: zdssh
+      !!
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zlon, zlat
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
+      INTEGER , DIMENSION(:), ALLOCATABLE :: ixpts, iypts, izpts, inpts
+      INTEGER :: jpts, npts
+
+      CALL wrk_alloc(jpi,jpj, zdssh )
+
+      ! get imbalance (volume heat and salt)
+      ! initialisation difference
+      zde3t = 0.0_wp; zdsal = 0.0_wp ; zdtem = 0.0_wp
+
+      ! initialisation correction term
+      pvol_flx(:,:,:  ) = 0.0_wp
+      pts_flx (:,:,:,:) = 0.0_wp
+      
+      r1_rdtiscpl = 1._wp / prdt_iscpl 
+
+      ! mask tsn and tsb 
+      tsb(:,:,:,jp_tem)=tsb(:,:,:,jp_tem)*ptmask_b(:,:,:); tsn(:,:,:,jp_tem)=tsn(:,:,:,jp_tem)*tmask(:,:,:);
+      tsb(:,:,:,jp_sal)=tsb(:,:,:,jp_sal)*ptmask_b(:,:,:); tsn(:,:,:,jp_sal)=tsn(:,:,:,jp_sal)*tmask(:,:,:);
+
+      !==============================================================================
+      ! diagnose the heat, salt and volume input and compute the correction variable
+      !==============================================================================
+
+      ! 
+      zdssh(:,:) = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+      IF ( lk_vvl ) zdssh = 0.0_wp ! already included in the levels by definition
+      
+      DO jk = 1,jpk-1
+         DO jj = 2,jpj-1
+            DO ji = fs_2,fs_jpim1
+               IF (tmask_h(ji,jj) == 1._wp) THEN
+
+                  ! volume differences
+                  zde3t = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk)
+
+                  ! heat diff
+                  zdtem = tsn(ji,jj,jk,jp_tem) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                        - tsb(ji,jj,jk,jp_tem) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
+                  ! salt diff
+                  zdsal = tsn(ji,jj,jk,jp_sal) * fse3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                        - tsb(ji,jj,jk,jp_sal) * pe3t_b (ji,jj,jk) * ptmask_b(ji,jj,jk)
+               
+                  ! shh changes
+                  IF ( ptmask_b(ji,jj,jk) == 1._wp .OR. tmask(ji,jj,jk) == 1._wp ) THEN 
+                     zde3t = zde3t + zdssh(ji,jj) ! zdssh = 0 if vvl
+                     zdssh(ji,jj) = 0._wp
+                  END IF
+
+                  ! volume, heat and salt differences in each cell 
+                  pvol_flx(ji,jj,jk)       =   pvol_flx(ji,jj,jk)        + zde3t * r1_rdtiscpl
+                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * r1_rdtiscpl 
+                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * r1_rdtiscpl
+
+                  ! case where we close a cell: check if the neighbour cells are wet 
+                  IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
+
+                     jip1=ji+1 ; jim1=ji-1 ; jjp1=jj+1 ; jjm1=jj-1 ;
+
+                     zsum =   e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
+                       &    + e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk)
+
+                     IF ( zsum /= 0._wp ) THEN
+                        zjip1_ratio   = e1e2t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
+                        zjim1_ratio   = e1e2t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
+                        zjjp1_ratio   = e1e2t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
+                        zjjm1_ratio   = e1e2t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) / zsum
+
+                        pvol_flx(ji  ,jjp1,jk       ) = pvol_flx(ji  ,jjp1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjp1_ratio
+                        pvol_flx(ji  ,jjm1,jk       ) = pvol_flx(ji  ,jjm1,jk       ) + pvol_flx(ji,jj,jk       ) * zjjm1_ratio
+                        pvol_flx(jip1,jj  ,jk       ) = pvol_flx(jip1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjip1_ratio
+                        pvol_flx(jim1,jj  ,jk       ) = pvol_flx(jim1,jj  ,jk       ) + pvol_flx(ji,jj,jk       ) * zjim1_ratio
+                        pts_flx (ji  ,jjp1,jk,jp_sal) = pts_flx (ji  ,jjp1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjp1_ratio
+                        pts_flx (ji  ,jjm1,jk,jp_sal) = pts_flx (ji  ,jjm1,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjjm1_ratio
+                        pts_flx (jip1,jj  ,jk,jp_sal) = pts_flx (jip1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjip1_ratio
+                        pts_flx (jim1,jj  ,jk,jp_sal) = pts_flx (jim1,jj  ,jk,jp_sal) + pts_flx (ji,jj,jk,jp_sal) * zjim1_ratio
+                        pts_flx (ji  ,jjp1,jk,jp_tem) = pts_flx (ji  ,jjp1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjp1_ratio
+                        pts_flx (ji  ,jjm1,jk,jp_tem) = pts_flx (ji  ,jjm1,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjjm1_ratio
+                        pts_flx (jip1,jj  ,jk,jp_tem) = pts_flx (jip1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjip1_ratio
+                        pts_flx (jim1,jj  ,jk,jp_tem) = pts_flx (jim1,jj  ,jk,jp_tem) + pts_flx (ji,jj,jk,jp_tem) * zjim1_ratio
+
+                        ! set to 0 the cell we distributed over neigbourg cells
+                        pvol_flx(ji,jj,jk       ) = 0._wp
+                        pts_flx (ji,jj,jk,jp_sal) = 0._wp
+                        pts_flx (ji,jj,jk,jp_tem) = 0._wp
+
+                     ELSE IF (zsum == 0._wp ) THEN
+                        ! case where we close a cell and no adjacent cell open
+                        ! check if the cell beneath is wet
+                        IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
+                           pvol_flx(ji,jj,jk+1)       =  pvol_flx(ji,jj,jk+1)        + pvol_flx(ji,jj,jk)
+                           pts_flx (ji,jj,jk+1,jp_sal)=  pts_flx (ji,jj,jk+1,jp_sal) + pts_flx (ji,jj,jk,jp_sal)
+                           pts_flx (ji,jj,jk+1,jp_tem)=  pts_flx (ji,jj,jk+1,jp_tem) + pts_flx (ji,jj,jk,jp_tem)
+
+                           ! set to 0 the cell we distributed over neigbourg cells
+                           pvol_flx(ji,jj,jk       ) = 0._wp
+                           pts_flx (ji,jj,jk,jp_sal) = 0._wp
+                           pts_flx (ji,jj,jk,jp_tem) = 0._wp
+                        ELSE
+                        ! case no adjacent cell on the horizontal and on the vertical
+                           IF ( lwp ) THEN   ! JMM : cAution this warning may occur on any mpp subdomain but numout is only
+                                             ! open for narea== 1 (lwp=T)
+                           WRITE(numout,*) 'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
+                           WRITE(numout,*) '                     ',mig(ji),' ',mjg(jj),' ',jk
+                           WRITE(numout,*) '                     ',ji,' ',jj,' ',jk,' ',narea
+                           WRITE(numout,*) ' we are now looking for the closest wet cell on the horizontal '
+                           ENDIF
+                        ! We deal with these points later.
+                        END IF
+                     END IF
+                  END IF
+               END IF
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1.)
+      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1.)
+      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1.)
+
+      ! if no neighbour wet cell in case of 2close a cell", need to find the nearest wet point 
+      ! allocation and initialisation of the list of problematic point
+      ALLOCATE(inpts(jpnij))
+      inpts(:)=0
+
+      ! fill narea location with the number of problematic point
+      DO jk = 1,jpk-1
+         DO jj = 2,jpj-1
+            DO ji = fs_2,fs_jpim1
+               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
+                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
+                  inpts(narea) = inpts(narea) + 1 
+               END IF
+            END DO
+         END DO
+      END DO
+
+      ! build array of total problematic point on each cpu (share to each cpu)
+      CALL mpp_max(inpts,jpnij) 
+
+      ! size of the new variable
+      npts  = SUM(inpts)    
+      
+      ! allocation of the coordinates, correction, index vector for the problematic points 
+      ALLOCATE(ixpts(npts), iypts(npts), izpts(npts), zcorr_vol(npts), zcorr_sal(npts), zcorr_tem(npts), zlon(npts), zlat(npts))
+      ixpts(:) = -9999 ; iypts(:) = -9999 ; izpts(:) = -9999 ; zlon(:) = -1.0e20_wp ; zlat(:) = -1.0e20_wp
+      zcorr_vol(:) = -1.0e20_wp
+      zcorr_sal(:) = -1.0e20_wp
+      zcorr_tem(:) = -1.0e20_wp
+
+      ! fill new variable
+      jpts = SUM(inpts(1:narea-1))
+      DO jk = 1,jpk-1
+         DO jj = 2,jpj-1
+            DO ji = fs_2,fs_jpim1
+               IF (     ptmask_b(ji,jj,jk) == 1._wp .AND. tmask(ji,jj,jk+1)  == 0._wp .AND. tmask_h(ji,jj) == 1._wp  &
+                  .AND. SUM(tmask(ji-1:ji+1,jj,jk)) + SUM(tmask(ji,jj-1:jj+1,jk)) == 0._wp) THEN
+                  jpts = jpts + 1  ! positioning in the inpts vector for the area narea
+                  ixpts(jpts) = ji           ; iypts(jpts) = jj ; izpts(jpts) = jk
+                  zlon (jpts) = glamt(ji,jj) ; zlat (jpts) = gphit(ji,jj)
+                  zcorr_vol(jpts) = pvol_flx(ji,jj,jk)
+                  zcorr_sal(jpts) = pts_flx (ji,jj,jk,jp_sal)
+                  zcorr_tem(jpts) = pts_flx (ji,jj,jk,jp_tem)
+
+                  ! set flx to 0 (safer)
+                  pvol_flx(ji,jj,jk       ) = 0.0_wp
+                  pts_flx (ji,jj,jk,jp_sal) = 0.0_wp
+                  pts_flx (ji,jj,jk,jp_tem) = 0.0_wp
+               END IF
+            END DO
+         END DO
+      END DO
+
+      ! build array of total problematic point on each cpu (share to each cpu)
+      ! point coordinates
+      CALL mpp_max(zlat ,npts)
+      CALL mpp_max(zlon ,npts)
+      CALL mpp_max(izpts,npts)
+
+      ! correction values 
+      CALL mpp_max(zcorr_vol,npts)
+      CALL mpp_max(zcorr_sal,npts)
+      CALL mpp_max(zcorr_tem,npts)
+
+      ! put correction term in the closest cell          
+      DO jpts = 1,npts
+         CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
+         DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
+            DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
+               jk = izpts(jpts)
+
+               IF (tmask_h(ji,jj) == 1._wp) THEN
+                  ! correct the vol_flx in the closest cell
+                  pvol_flx(ji,jj,jk)        =  pvol_flx(ji,jj,jk       ) + zcorr_vol(jpts)
+                  pts_flx (ji,jj,jk,jp_sal) =  pts_flx (ji,jj,jk,jp_sal) + zcorr_sal(jpts)
+                  pts_flx (ji,jj,jk,jp_tem) =  pts_flx (ji,jj,jk,jp_tem) + zcorr_tem(jpts)
+
+                  ! set correction to 0
+                  zcorr_vol(jpts) = 0.0_wp
+                  zcorr_sal(jpts) = 0.0_wp
+                  zcorr_tem(jpts) = 0.0_wp
+               END IF
+            END DO
+         END DO
+      END DO
+
+      ! deallocate variables 
+      DEALLOCATE(inpts)
+      DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
+    
+      ! add contribution store on the hallo (lbclnk remove one of the contribution)
+      pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
+      pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
+      pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
+
+      ! compute sum over the halo and set it to 0.
+      CALL lbc_sum(pvol_flx(:,:,:       ),'T',1._wp)
+      CALL lbc_sum(pts_flx (:,:,:,jp_sal),'T',1._wp)
+      CALL lbc_sum(pts_flx (:,:,:,jp_tem),'T',1._wp)
+
+      ! deallocate variables
+      CALL wrk_dealloc(jpi,jpj, zdssh ) 
+
+   END SUBROUTINE iscpl_cons
+
+   SUBROUTINE iscpl_div( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE iscpl_div  ***
+      !!
+      !! ** Purpose :   update the horizontal divergenc
+      !!
+      !! ** Method  :
+      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   increase by the iscpl correction term
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE iscpl_div
+
+END MODULE iscplhsb
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 6006)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 6006)
@@ -0,0 +1,87 @@
+MODULE iscplini
+   !!======================================================================
+   !!                       ***  MODULE  sbciscpl***
+   !! Ocean forcing:  river runoff
+   !!=====================================================================
+   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   iscpl_init     : initialisation routine (namelist)
+   !!   iscpl_alloc    : allocation of correction variables
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE oce             ! global tra/dyn variable
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! MPP library
+   USE in_out_manager  ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   iscpl_init      
+   PUBLIC   iscpl_alloc 
+   !!                                                      !!* namsbc_iscpl namelist *
+   LOGICAL , PUBLIC                                        ::   ln_hsb
+   INTEGER , PUBLIC                                        ::   nn_fiscpl, nstp_iscpl
+   INTEGER , PUBLIC                                        ::   nn_drown
+   REAL(wp), PUBLIC                                        ::   rdt_iscpl
+   !!                                                      !!* namsbc_iscpl namelist *
+   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)
+   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION iscpl_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE sbc_iscpl_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc )
+         !
+      IF( lk_mpp          )   CALL mpp_sum ( iscpl_alloc )
+      IF( iscpl_alloc > 0 )   CALL ctl_warn('iscpl_alloc: allocation of arrays failed')
+   END FUNCTION iscpl_alloc
+
+   SUBROUTINE iscpl_init()
+      INTEGER ::   ios           ! Local integer output status for namelist read
+      NAMELIST/namsbc_iscpl/nn_fiscpl,ln_hsb,nn_drown
+      !!----------------------------------------------------------------------
+      !                                   ! ============
+      !                                   !   Namelist
+      !                                   ! ============
+      !
+      nn_fiscpl = 0
+      ln_hsb    = .FALSE.
+      REWIND( numnam_ref )              ! Namelist namsbc_iscpl in reference namelist : Ice sheet coupling
+      READ  ( numnam_ref, namsbc_iscpl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_iscpl in configuration namelist : Ice Sheet coupling
+      READ  ( numnam_cfg, namsbc_iscpl, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_iscpl in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc_iscpl )
+      !
+      nstp_iscpl=MIN(nn_fiscpl, nitend-nit000+1) ! the coupling period have to be less or egal than the total number of time step
+      rdt_iscpl = nstp_iscpl * rn_rdt
+      !
+      IF (lwp) THEN
+         WRITE(numout,*) 'iscpl_rst:'
+         WRITE(numout,*) '~~~~~~~~~'
+         WRITE(numout,*) ' coupling     flag (ln_iscpl )            = ', ln_iscpl
+         WRITE(numout,*) ' conservation flag (ln_hsb   )            = ', ln_hsb
+         WRITE(numout,*) ' nb of stp for cons (rn_fiscpl)           = ', nstp_iscpl
+         IF (nstp_iscpl .NE. nn_fiscpl) WRITE(numout,*) 'W A R N I N G: nb of stp for cons has been modified &
+            &                                           (larger than run length)'
+         WRITE(numout,*) ' coupling time step                       = ', rdt_iscpl
+         WRITE(numout,*) ' number of call of the extrapolation loop = ', nn_drown
+      END IF
+
+   END SUBROUTINE iscpl_init
+
+END MODULE iscplini
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 6006)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 6006)
@@ -0,0 +1,442 @@
+MODULE iscplrst
+   !!======================================================================
+   !!                       ***  MODULE  iscplrst***
+   !! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
+   !!=====================================================================
+   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   iscpl_stp          : step management
+   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE domwri          ! ocean space and time domain
+   USE domvvl, ONLY : dom_vvl_interpol
+   USE phycst          ! physical constants
+   USE sbc_oce         ! surface boundary condition variables
+   USE oce             ! global tra/dyn variable
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! MPP library
+   USE wrk_nemo        ! Memory allocation
+   USE lbclnk          ! communication
+   USE iscplini        ! ice sheet coupling: initialisation
+   USE iscplhsb        ! ice sheet coupling: conservation
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   iscpl_stp          ! step management 
+   !!
+   !! * Substitutions  
+#  include "domzgr_substitute.h90"  
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE iscpl_stp
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_stp  ***
+      !! 
+      !! ** Purpose : compute initialisation
+      !!              compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)   
+      !!              compute correction term if needed
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER  ::   inum0
+      REAL(wp), DIMENSION(:,:  ), POINTER :: zsmask_b
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmask_b, zumask_b, zvmask_b
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ze3t_b  , ze3u_b  , ze3v_b  
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zdepw_b
+      CHARACTER(20) :: cfile
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! mask before
+      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! e3   before
+      CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
+      CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
+
+
+      !! get restart variable
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b   ) ! need to extrapolate T/S
+      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b   ) ! need to correct barotropic velocity
+      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b   ) ! need to correct barotropic velocity
+      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b   ) ! need to correct barotropic velocity
+      CALL iom_get( numror, jpdom_autoglo, 'fse3t_n' , ze3t_b(:,:,:) )  ! need to compute temperature correction
+      CALL iom_get( numror, jpdom_autoglo, 'fse3u_n' , ze3u_b(:,:,:) )  ! need to correct barotropic velocity
+      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n' , ze3v_b(:,:,:) )  ! need to correct barotropic velocity
+      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
+
+      !! read namelist
+      CALL iscpl_init()
+
+      !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
+      CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
+
+      !! compute correction if conservation needed
+      IF ( ln_hsb ) THEN
+         IF( iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
+         CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
+      END IF
+         
+      !! print mesh/mask
+      IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
+
+      IF ( ln_hsb ) THEN
+         cfile='correction'
+         cfile = TRIM( cfile )
+         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
+         CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
+         CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
+         CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
+         CALL iom_close ( inum0 )
+      END IF
+
+      CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b )  
+      CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b  ,ze3u_b  ,ze3v_b   )  
+      CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b                    )
+      CALL wrk_dealloc(jpi,jpj,     zsmask_b                   )
+
+      !! next step is an euler time step
+      neuler = 0
+
+      !! set _b and _n variables equal
+      tsb (:,:,:,:) = tsn (:,:,:,:)
+      ub  (:,:,:  ) = un  (:,:,:  )
+      vb  (:,:,:  ) = vn  (:,:,:  )
+      sshb(:,:    ) = sshn(:,:)
+
+      !! set _b and _n vertical scale factor equal
+      fse3t_b (:,:,:) = fse3t_n (:,:,:)
+      fse3u_b (:,:,:) = fse3u_n (:,:,:)
+      fse3v_b (:,:,:) = fse3v_n (:,:,:)
+
+      IF ( lk_vvl ) THEN
+         fse3uw_b(:,:,:) = fse3uw_n(:,:,:)
+         fse3vw_b(:,:,:) = fse3vw_n(:,:,:)
+         fsdept_b(:,:,:) = fsdept_n(:,:,:)
+         fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+         hu_b (:,:) = hu(:,:)
+         hv_b (:,:) = hv(:,:)
+         hur_b(:,:) = hur(:,:)
+         hvr_b(:,:) = hvr(:,:)
+      END IF
+      !
+   END SUBROUTINE iscpl_stp
+   
+   SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_rst_interpol  ***
+      !! 
+      !! ** Purpose :   compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves 
+      !!                compute 2d fields of heat, salt and volume correction
+      !! 
+      !! ** Method  :   tn, sn : extrapolation from neigbourg cells
+      !!                un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b    !! mask before
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pe3t_b  , pe3u_b  , pe3v_b      !! scale factor before
+      REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
+      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
+      !!
+      INTEGER :: ji, jj, jk, iz          !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
+      !!
+      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
+      REAL(wp):: zdz, zdzm1, zdzp1
+      !!
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, zhu1, zhv1
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
+      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
+      !!----------------------------------------------------------------------
+
+      !! allocate variables
+      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
+      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d        ) 
+      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
+      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
+      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
+      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, zhu1, zhv1             ) 
+
+      !! mask value to be sure
+      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
+      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
+      
+      ! compute wmask
+      zwmaskn(:,:,1) = tmask   (:,:,1)
+      zwmaskb(:,:,1) = ptmask_b(:,:,1)
+      DO jk = 2,jpk
+         zwmaskn(:,:,jk) =  tmask  (:,:,jk) *  tmask  (:,:,jk-1)
+         zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
+      END DO
+           
+      ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
+      sshb (:,:)=sshn(:,:)
+      zssh0(:,:)=sshn(:,:)
+      zsmask0(:,:) = psmask_b(:,:)
+      zsmask1(:,:) = psmask_b(:,:)
+      DO iz = 1,10    ! need to be tuned (configuration dependent) (OK for ISOMIP+)
+         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
+         DO jj = 2,jpj-1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               jip1=ji+1; jim1=ji-1;
+               jjp1=jj+1; jjm1=jj-1;
+               summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
+               IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
+                  sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj)     &
+                  &           + zssh0(jim1,jj)*zsmask0(jim1,jj)     &
+                  &           + zssh0(ji,jjp1)*zsmask0(ji,jjp1)     &
+                  &           + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
+                  zsmask1(ji,jj)=1._wp
+               END IF
+            END DO
+         END DO
+         CALL lbc_lnk(sshn,'T',1._wp)
+         CALL lbc_lnk(zsmask1,'T',1._wp)
+         zssh0   = sshn
+         zsmask0 = zsmask1
+      END DO
+      sshn(:,:) = sshn(:,:) * ssmask(:,:)
+
+!=============================================================================
+      IF ( lk_vvl ) THEN
+      ! Reconstruction of all vertical scale factors at now time steps
+      ! =============================================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+         DO jk = 1,jpk
+            DO jj=1,jpj
+               DO ji=1,jpi
+                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
+                     fse3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+
+      ! t- and w- points depth
+      ! ----------------------
+         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+         fsdepw_n(:,:,1) = 0.0_wp
+         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - 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) > 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
+
+      ! t-, u- and v- water column thickness
+      ! ------------------------------------
+         ht(:,:) = 0._wp ; hu(:,:) = 0._wp ; hv(:,:) = 0._wp
+         DO jk = 1, jpkm1
+            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         !                                        ! Inverse of the local depth
+         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
+
+      END IF
+
+!=============================================================================
+! compute velocity
+! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
+      ub(:,:,:)=un(:,:,:)
+      vb(:,:,:)=vn(:,:,:)
+      DO jk = 1,jpk
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/fse3u_n(ji,jj,jk)*umask(ji,jj,jk);
+               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/fse3v_n(ji,jj,jk)*vmask(ji,jj,jk);
+            END DO
+         END DO
+      END DO
+
+! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
+! compute barotropic velocity now and after 
+      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
+      zbub(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
+      zbvb(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
+      zbun(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
+      zbvn(:,:)   = SUM(ztrp,DIM=3)
+
+      ! new water column
+      zhu1=0.0_wp ;
+      zhv1=0.0_wp ;
+      DO jk  = 1,jpk
+        zhu1(:,:) = zhu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+        zhv1(:,:) = zhv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      
+      ! compute correction      
+      zucorr = 0._wp
+      zvcorr = 0._wp
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
+               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
+            END IF
+            IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
+               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
+            END IF
+         END DO
+      END DO 
+      
+      ! update velocity
+      DO jk  = 1,jpk
+         un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
+         vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
+      END DO
+
+!=============================================================================
+      ! compute temp and salt
+      ! compute new tn and sn if we open a new cell
+      tsb (:,:,:,:) = tsn(:,:,:,:)
+      zts0(:,:,:,:) = tsn(:,:,:,:)
+      ztmask1(:,:,:) = ptmask_b(:,:,:)
+      ztmask0(:,:,:) = ptmask_b(:,:,:)
+      DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
+          DO jk = 1,jpk-1
+              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
+              DO jj = 2,jpj-1
+                 DO ji = fs_2,fs_jpim1
+                      jip1=ji+1; jim1=ji-1;
+                      jjp1=jj+1; jjm1=jj-1;
+                      summsk= (ztmask0(jip1,jj  ,jk)+ztmask0(jim1,jj  ,jk)+ztmask0(ji  ,jjp1,jk)+ztmask0(ji  ,jjm1,jk))
+                      IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
+                      !! horizontal basic extrapolation
+                         tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1)*ztmask0(jip1,jj  ,jk) &
+                         &                +zts0(jim1,jj  ,jk,1)*ztmask0(jim1,jj  ,jk) &
+                         &                +zts0(ji  ,jjp1,jk,1)*ztmask0(ji  ,jjp1,jk) &
+                         &                +zts0(ji  ,jjm1,jk,1)*ztmask0(ji  ,jjm1,jk) ) / summsk
+                         tsn(ji,jj,jk,2)=( zts0(jip1,jj  ,jk,2)*ztmask0(jip1,jj  ,jk) &
+                         &                +zts0(jim1,jj  ,jk,2)*ztmask0(jim1,jj  ,jk) &
+                         &                +zts0(ji  ,jjp1,jk,2)*ztmask0(ji  ,jjp1,jk) &
+                         &                +zts0(ji  ,jjm1,jk,2)*ztmask0(ji  ,jjm1,jk) ) / summsk
+                         ztmask1(ji,jj,jk)=1
+                      ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
+                      !! vertical extrapolation if horizontal extrapolation failed
+                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
+                         IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
+                            tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                            &                +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
+                            tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                            &                +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
+                            ztmask1(ji,jj,jk)=1._wp
+                         END IF
+                      END IF
+                  END DO
+              END DO
+          END DO
+          
+          CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
+          CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
+          CALL lbc_lnk(ztmask1,     'T',1._wp)
+
+          ! update
+          zts0(:,:,:,:) = tsn(:,:,:,:)
+          ztmask0 = ztmask1
+
+      END DO
+
+      ! mask new tsn field
+      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
+      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
+
+      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+      IF ( lk_vvl ) THEN
+         DO jk = 2,jpk-1
+            DO jj = 1,jpj
+               DO ji = 1,jpi
+                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
+                     !compute weight
+                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
+                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
+                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  )  - fsdepw_n(ji,jj,jk  ))
+                     IF (zdz .LT. 0._wp) THEN 
+                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
+                     END IF
+                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
+                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
+                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
+                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
+                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
+                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
+                  END IF
+               END DO
+            END DO
+         END DO               
+      END IF
+
+      ! closed pool
+      ! -----------------------------------------------------------------------------------------
+      ! case we open a cell but no neigbour cells available to get an estimate of T and S
+      WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp) 
+         tsn(:,:,:,2) = -99._wp  ! Special value for closed pool (checking purpose in output.init)
+         tmask(:,:,:) = 0._wp    ! set mask to 0 to run
+         umask(:,:,:) = 0._wp
+         vmask(:,:,:) = 0._wp
+      END WHERE
+      
+      ! set mbkt and mikt to 1 in thiese location
+      WHERE (SUM(tmask,dim=3) == 0)
+         mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
+         mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
+      END WHERE
+      ! -------------------------------------------------------------------------------------------
+      ! compute new tn and sn if we close cell 
+      ! nothing to do
+      ! 
+      ! deallocation tmp arrays
+      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
+      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
+      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
+      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
+      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
+      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
+
+   END SUBROUTINE iscpl_rst_interpol
+
+END MODULE iscplrst
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6006)
@@ -35,4 +35,5 @@
    USE dtauvd          ! data: U & V current             (dta_uvd routine)
    USE domvvl          ! varying vertical mesh
+   USE iscplrst        ! ice sheet coupling
    !
    USE in_out_manager  ! I/O manager
@@ -85,6 +86,7 @@
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
-         CALL rst_read                           ! Read the restart file
-         CALL day_init                           ! model calendar (using both namelist and restart infos)
+         CALL rst_read                        ! Read the restart file
+         IF (ln_iscpl)       CALL iscpl_stp   ! extraloate restart to wet and dry
+         CALL day_init                        ! model calendar (using both namelist and restart infos)
       ELSE
          !                                    ! Start from rest
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90	(revision 6006)
@@ -23,4 +23,6 @@
    USE sbcrnf          ! river runoff 
    USE sbcisf          ! ice shelf
+   USE iscplhsb        ! ice sheet / ocean coupling
+   USE iscplini        ! ice sheet / ocean coupling
    !
    USE in_out_manager  ! I/O manager
@@ -93,4 +95,6 @@
       IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field)
       !
+      IF( ln_iscpl  .AND. ln_hsb     )   CALL iscpl_div( hdivn )   !==  ice sheet  ==!   (update hdivn field)
+      !
       CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change)
       !
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 6006)
@@ -95,5 +95,4 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   iku, ikv     ! local integers
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
@@ -320,6 +319,6 @@
             hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
          END DO
-         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
-         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
+         hur_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
+         hvr_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
       ENDIF
       !
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 6006)
@@ -142,8 +142,20 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , sshn      )
                      CALL iom_rstput( kt, nitrst, numrow, 'rhop'   , rhop      )
+
+                  ! extra variable needed for the ice sheet coupling
+                  IF ( ln_iscpl ) THEN 
+                     CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask     ) ! need to extrapolate T/S
+                     CALL iom_rstput( kt, nitrst, numrow, 'umask'  , umask     ) ! need to correct barotropic velocity
+                     CALL iom_rstput( kt, nitrst, numrow, 'vmask'  , vmask     ) ! need to correct barotropic velocity
+                     CALL iom_rstput( kt, nitrst, numrow, 'smask'  , ssmask    ) ! need to correct barotropic velocity
+                     CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )   ! need to compute temperature correction
+                     CALL iom_rstput( kt, nitrst, numrow, 'fse3u_n', fse3u_n(:,:,:) )   ! need to compute bt conservation
+                     CALL iom_rstput( kt, nitrst, numrow, 'fse3v_n', fse3v_n(:,:,:) )   ! need to compute bt conservation
+                     CALL iom_rstput( kt, nitrst, numrow, 'fsdepw_n', fsdepw_n(:,:,:) ) ! need to compute extrapolation if vvl
+                  END IF
       ENDIF
       
       IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst    )  
-      
+
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -209,5 +221,4 @@
       REAL(wp) ::   zrdt, zrdttra1
       INTEGER  ::   jk
-      LOGICAL  ::   llok
       !!----------------------------------------------------------------------
 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 6006)
@@ -17,4 +17,5 @@
    !!----------------------------------------------------------------------
    !!   lbc_lnk      : generic interface for mpp_lnk_3d and mpp_lnk_2d routines defined in lib_mpp
+   !!   lbc_sum      : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d routines defined in lib_mpp
    !!   lbc_lnk_e    : generic interface for mpp_lnk_2d_e routine defined in lib_mpp
    !!   lbc_bdy_lnk  : generic interface for mpp_lnk_bdy_2d and mpp_lnk_bdy_3d routines defined in lib_mpp
@@ -31,4 +32,9 @@
    END INTERFACE
 
+!JMM interface not defined if not key_mpp_mpi : likely do not compile without this CPP key !!!!
+   INTERFACE lbc_sum
+      MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d
+   END INTERFACE
+
    INTERFACE lbc_bdy_lnk
       MODULE PROCEDURE mpp_lnk_bdy_2d, mpp_lnk_bdy_3d
@@ -45,4 +51,5 @@
    PUBLIC lbc_lnk       ! ocean lateral boundary conditions
    PUBLIC lbc_lnk_multi ! modified ocean lateral boundary conditions
+   PUBLIC lbc_sum
    PUBLIC lbc_lnk_e
    PUBLIC lbc_bdy_lnk   ! ocean lateral BDY boundary conditions
@@ -59,8 +66,11 @@
    !!   Default option                              shared memory computing
    !!----------------------------------------------------------------------
-   !!   lbc_lnk      : generic interface for lbc_lnk_3d and lbc_lnk_2d
-   !!   lbc_lnk_3d   : set the lateral boundary condition on a 3D variable on ocean mesh
-   !!   lbc_lnk_2d   : set the lateral boundary condition on a 2D variable on ocean mesh
-   !!   lbc_bdy_lnk  : set the lateral BDY boundary condition
+   !!   lbc_sum       : generic interface for mpp_lnk_sum_3d and mpp_lnk_sum_2d 
+   !!   lbc_lnk_sum_3d: compute sum over the halos on a 3D variable on ocean mesh
+   !!   lbc_lnk_sum_3d: compute sum over the halos on a 2D variable on ocean mesh
+   !!   lbc_lnk       : generic interface for lbc_lnk_3d and lbc_lnk_2d
+   !!   lbc_lnk_3d    : set the lateral boundary condition on a 3D variable on ocean mesh
+   !!   lbc_lnk_2d    : set the lateral boundary condition on a 2D variable on ocean mesh
+   !!   lbc_bdy_lnk   : set the lateral BDY boundary condition
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers   
@@ -74,4 +84,8 @@
    INTERFACE lbc_lnk
       MODULE PROCEDURE lbc_lnk_3d_gather, lbc_lnk_3d, lbc_lnk_2d
+   END INTERFACE
+
+   INTERFACE lbc_sum
+      MODULE PROCEDURE mpp_lnk_sum_3d, mpp_lnk_sum_2d
    END INTERFACE
 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 6006)
@@ -72,4 +72,5 @@
    PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e
    PUBLIC   mpp_lnk_2d_9 
+   PUBLIC   mpp_lnk_sum_3d, mpp_lnk_sum_2d
    PUBLIC   mppscatter, mppgather
    PUBLIC   mpp_ini_ice, mpp_ini_znl
@@ -1402,4 +1403,347 @@
    END SUBROUTINE mpp_lnk_2d_e
 
+   SUBROUTINE mpp_lnk_sum_3d( ptab, cd_type, psgn, cd_mpp, pval )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_sum_3d  ***
+      !!
+      !! ** Purpose :   Message passing manadgement (sum the overlap region)
+      !!
+      !! ** Method  :   Use mppsend and mpprecv function for passing mask
+      !!      between processors following neighboring subdomains.
+      !!            domain parameters
+      !!                    nlci   : first dimension of the local subdomain
+      !!                    nlcj   : second dimension of the local subdomain
+      !!                    nbondi : mark for "east-west local boundary"
+      !!                    nbondj : mark for "north-south local boundary"
+      !!                    noea   : number for local neighboring processors
+      !!                    nowe   : number for local neighboring processors
+      !!                    noso   : number for local neighboring processors
+      !!                    nono   : number for local neighboring processors
+      !!
+      !! ** Action  :   ptab with update value at its periphery
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
+      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
+      !                                                             ! = T , U , V , F , W points
+      REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
+      !                                                             ! =  1. , the sign is kept
+      CHARACTER(len=3), OPTIONAL      , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
+      REAL(wp)        , OPTIONAL      , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
+      !!
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
+      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
+      REAL(wp) ::   zland
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
+      !
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east
+
+      !!----------------------------------------------------------------------
+      
+      ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   &
+         &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  )
+
+      !
+      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
+      ELSE                         ;   zland = 0.e0      ! zero by default
+      ENDIF
+
+      ! 1. standard boundary treatment
+      ! ------------------------------
+      ! 2. East and west directions exchange
+      ! ------------------------------------
+      ! we play with the neigbours AND the row number because of the periodicity
+      !
+      SELECT CASE ( nbondi )      ! Read lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+      iihom = nlci-jpreci
+         DO jl = 1, jpreci
+            zt3ew(:,jl,:,1) = ptab(jl      ,:,:) ; ptab(jl      ,:,:) = 0.0_wp
+            zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp 
+         END DO
+      END SELECT
+      !
+      !                           ! Migrations
+      imigr = jpreci * jpj * jpk
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      END SELECT
+      !
+      !                           ! Write lateral conditions
+      iihom = nlci-nreci
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, jpreci
+            ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jpreci
+            ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2)
+            ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jpreci
+            ptab(jpreci+jl,:,:) = ptab(jpreci+jl,:,:) + zt3we(:,jl,:,2)
+         END DO
+      END SELECT
+
+
+      ! 3. North and south directions
+      ! -----------------------------
+      ! always closed : we play only with the neigbours
+      !
+      IF( nbondj /= 2 ) THEN      ! Read lateral conditions
+         ijhom = nlcj-jprecj
+         DO jl = 1, jprecj
+            zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp
+            zt3ns(:,jl,:,1) = ptab(:,jl      ,:) ; ptab(:,jl      ,:) = 0.0_wp
+         END DO
+      ENDIF
+      !
+      !                           ! Migrations
+      imigr = jprecj * jpi * jpk
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      CASE ( 0 )
+         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      END SELECT
+      !
+      !                           ! Write lateral conditions
+      ijhom = nlcj-nrecj
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, jprecj
+            ptab(:,ijhom+jl,:) = ptab(:,ijhom+jl,:) + zt3ns(:,jl,:,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jprecj
+            ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl,:,2)
+            ptab(:,ijhom +jl,:) = ptab(:,ijhom +jl,:) + zt3ns(:,jl,:,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jprecj
+            ptab(:,jprecj+jl,:) = ptab(:,jprecj+jl,:) + zt3sn(:,jl   ,:,2)
+         END DO
+      END SELECT
+
+
+      ! 4. north fold treatment
+      ! -----------------------
+      !
+      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+         !
+         SELECT CASE ( jpni )
+         CASE ( 1 )     ;   CALL lbc_nfd      ( ptab, cd_type, psgn )   ! only 1 northern proc, no mpp
+         CASE DEFAULT   ;   CALL mpp_lbc_north( ptab, cd_type, psgn )   ! for all northern procs.
+         END SELECT
+         !
+      ENDIF
+      !
+      DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we )
+      !
+   END SUBROUTINE mpp_lnk_sum_3d
+
+   SUBROUTINE mpp_lnk_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_sum_2d  ***
+      !!
+      !! ** Purpose :   Message passing manadgement for 2d array (sum the overlap region)
+      !!
+      !! ** Method  :   Use mppsend and mpprecv function for passing mask
+      !!      between processors following neighboring subdomains.
+      !!            domain parameters
+      !!                    nlci   : first dimension of the local subdomain
+      !!                    nlcj   : second dimension of the local subdomain
+      !!                    nbondi : mark for "east-west local boundary"
+      !!                    nbondj : mark for "north-south local boundary"
+      !!                    noea   : number for local neighboring processors
+      !!                    nowe   : number for local neighboring processors
+      !!                    noso   : number for local neighboring processors
+      !!                    nono   : number for local neighboring processors
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d     ! 2D array on which the boundary condition is applied
+      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
+      !                                                         ! = T , U , V , F , W and I points
+      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
+      !                                                         ! =  1. , the sign is kept
+      CHARACTER(len=3), OPTIONAL  , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
+      REAL(wp)        , OPTIONAL  , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
+      !!
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
+      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
+      REAL(wp) ::   zland
+      INTEGER, DIMENSION(MPI_STATUS_SIZE) ::   ml_stat   ! for key_mpi_isend
+      !
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east
+
+      !!----------------------------------------------------------------------
+
+      ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2),  &
+         &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   )
+
+      !
+      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
+      ELSE                         ;   zland = 0.e0      ! zero by default
+      ENDIF
+
+      ! 1. standard boundary treatment
+      ! ------------------------------
+      ! 2. East and west directions exchange
+      ! ------------------------------------
+      ! we play with the neigbours AND the row number because of the periodicity
+      !
+      SELECT CASE ( nbondi )      ! Read lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+         iihom = nlci - jpreci
+         DO jl = 1, jpreci
+            zt2ew(:,jl,1) = pt2d(jl       ,:) ; pt2d(jl       ,:) = 0.0_wp
+            zt2we(:,jl,1) = pt2d(iihom +jl,:) ; pt2d(iihom +jl,:) = 0.0_wp
+         END DO
+      END SELECT
+      !
+      !                           ! Migrations
+      imigr = jpreci * jpj
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      !                           ! Write lateral conditions
+      iihom = nlci-nreci
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, jpreci
+            pt2d(iihom+jl,:) = pt2d(iihom+jl,:) + zt2ew(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jpreci
+            pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2)
+            pt2d(iihom +jl,:) = pt2d(iihom +jl,:) + zt2ew(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jpreci
+            pt2d(jpreci+jl,:) = pt2d(jpreci+jl,:) + zt2we(:,jl,2)
+         END DO
+      END SELECT
+
+
+      ! 3. North and south directions
+      ! -----------------------------
+      ! always closed : we play only with the neigbours
+      !
+      IF( nbondj /= 2 ) THEN      ! Read lateral conditions
+         ijhom = nlcj - jprecj
+         DO jl = 1, jprecj
+            zt2sn(:,jl,1) = pt2d(:,ijhom +jl) ; pt2d(:,ijhom +jl) = 0.0_wp
+            zt2ns(:,jl,1) = pt2d(:,jl       ) ; pt2d(:,jl       ) = 0.0_wp
+         END DO
+      ENDIF
+      !
+      !                           ! Migrations
+      imigr = jprecj * jpi
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      END SELECT
+      !
+      !                           ! Write lateral conditions
+      ijhom = nlcj-nrecj
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, jprecj
+            pt2d(:,ijhom+jl) = pt2d(:,ijhom+jl) + zt2ns(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jprecj
+            pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2)
+            pt2d(:,ijhom +jl) = pt2d(:,ijhom +jl) + zt2ns(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jprecj
+            pt2d(:,jprecj+jl) = pt2d(:,jprecj+jl) + zt2sn(:,jl,2)
+         END DO
+      END SELECT
+
+
+      ! 4. north fold treatment
+      ! -----------------------
+      !
+      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+         !
+         SELECT CASE ( jpni )
+         CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d, cd_type, psgn )   ! only 1 northern proc, no mpp
+         CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d, cd_type, psgn )   ! for all northern procs.
+         END SELECT
+         !
+      ENDIF
+      !
+      DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
+      !
+   END SUBROUTINE mpp_lnk_sum_2d
 
    SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req )
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 6006)
@@ -136,5 +136,5 @@
 
       imask(:,:)=1
-      WHERE ( zdta(:,:) - zdtaisf(:,:) <= 0. ) imask = 0
+      WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0
 
       !  1. Dimension arrays for subdomains
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 6006)
@@ -119,4 +119,6 @@
          IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
          IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+         IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0  
+         IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0  
          IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
          IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
@@ -230,11 +232,11 @@
             CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
             CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
-            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
-            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
+            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
+            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
             ! iom print
             CALL iom_put('ttbl',ttbl(:,:))
             CALL iom_put('stbl',stbl(:,:))
-            CALL iom_put('utbl',utbl(:,:))
-            CALL iom_put('vtbl',vtbl(:,:))
+            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
+            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
             ! compute fwf and heat flux
             CALL sbc_isf_cav (kt)
@@ -356,6 +358,4 @@
     IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
      !
-
-    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
     DO ji = 1, jpi
        DO jj = 1, jpj
@@ -366,9 +366,6 @@
     ! 3. -----------the average temperature between 200m and 600m ---------------------
              DO jk = misfkt(ji,jj),misfkb(ji,jj)
-             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
-             ! after verif with UNESCO, wrong sign in BG eq. 2
              ! Calculate freezing temperature
-                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
-                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
+                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik)) 
                 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
              ENDDO
@@ -418,5 +415,5 @@
       REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
       REAL(wp) ::   zgammat, zgammas
-      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
+      REAL(wp) ::   zeps = 1.e-20_wp        !==   Local constant initialization   ==!
       INTEGER  ::   ji, jj     ! dummy loop indices
       INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
@@ -503,10 +500,10 @@
                      zeps1=rcp*rau0*zgammat
                      zeps2=lfusisf*rau0*zgammas
-                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
+                     zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps)
                      zeps4=zlamb2+zlamb3*risfdep(ji,jj)
                      zeps6=zeps4-zti(ji,jj)
                      zeps7=zeps4-tsurf
                      zaqe=zlamb1 * (zeps1 + zeps3)
-                     zaqer=0.5/zaqe
+                     zaqer=0.5/MIN(zaqe,-zeps)
                      zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
                      zcqe=zeps2*stbl(ji,jj)
@@ -517,6 +514,6 @@
                      zfrz(ji,jj)=zeps4+zlamb1*zsfrz
   
-! zfwflx is upward water flux
-                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
+! zfwflx is upward water flux (MAX usefull if kappa = 0
+                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) )
 ! zhtflx is upward heat flux (out of ocean)
 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90	(revision 6006)
@@ -100,5 +100,4 @@
       !
       CALL wrk_alloc( jpi,jpj,jpk,jpts,   zts_dta )
-      !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi,jpj,jpk,jpts,   ztrdts ) 
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 6006)
@@ -22,4 +22,5 @@
    USE sbcrnf          ! River runoff  
    USE sbcisf          ! Ice shelf   
+   USE iscplini        ! Ice sheet coupling
    USE traqsr          ! solar radiation penetration
    USE trd_oce         ! trends: ocean variables
@@ -117,5 +118,4 @@
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices  
       INTEGER  ::   ikt, ikb 
-      INTEGER  ::   nk_isf
       REAL(wp) ::   zfact, z1_e3t, zdep
       REAL(wp) ::   zalpha, zhk
@@ -219,5 +219,5 @@
       !
       IF( nn_isf > 0 ) THEN
-         zfact = 0.5e0
+         zfact = 0.5_wp
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1
@@ -230,5 +230,4 @@
                ! sign - because fwf sign of evapo (rnf sign of precip)
                DO jk = ikt, ikb - 1
-               ! compute tfreez for the temperature correction (we add water at freezing temperature)
                ! compute trend
                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
@@ -239,5 +238,4 @@
    
                ! level partially include in ice shelf boundary layer 
-               ! compute tfreez for the temperature correction (we add water at freezing temperature)
                ! compute trend
                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
@@ -279,5 +277,23 @@
       ENDIF
  
-      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
+      !----------------------------------------
+      !        Ice Sheet coupling imbalance correction to have conservation
+      !----------------------------------------
+      !
+      IF( ln_iscpl .AND. ln_hsb) THEN         ! input of heat and salt due to river runoff 
+         DO jk = 1,jpk
+            DO jj = 2, jpj 
+               DO ji = fs_2, fs_jpim1
+                  zdep = 1._wp / fse3t_n(ji,jj,jk) 
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) - htsc_iscpl(ji,jj,jk,jp_tem)                       &
+                      &                                         * zdep
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) - htsc_iscpl(ji,jj,jk,jp_sal)                       &
+                      &                                         * zdep  
+               END DO  
+            END DO  
+         END DO
+      ENDIF
+
+      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
Index: /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90
===================================================================
--- /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 6005)
+++ /branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/lib_fortran.F90	(revision 6006)
@@ -24,6 +24,7 @@
    PRIVATE
 
-   PUBLIC   glob_sum   ! used in many places
-   PUBLIC   DDPDD      ! also used in closea module
+   PUBLIC   glob_sum      ! used in many places (masked with tmask_i)
+   PUBLIC   glob_sum_full ! used in many places (masked with tmask_h, ie omly over the halos)
+   PUBLIC   DDPDD         ! also used in closea module
    PUBLIC   glob_min, glob_max
 #if defined key_nosignedzero
@@ -34,4 +35,7 @@
       MODULE PROCEDURE glob_sum_1d, glob_sum_2d, glob_sum_3d, &
          &             glob_sum_2d_a, glob_sum_3d_a
+   END INTERFACE
+   INTERFACE glob_sum_full
+      MODULE PROCEDURE glob_sum_full_2d, glob_sum_full_3d
    END INTERFACE
    INTERFACE glob_min
@@ -156,4 +160,44 @@
       !
    END FUNCTION glob_sum_3d_a
+
+   FUNCTION glob_sum_full_2d( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_full_2d ***
+      !!
+      !! ** Purpose : perform a sum in calling DDPDD routine (nomask)
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_sum_full_2d   ! global sum
+      !!
+      !!-----------------------------------------------------------------------
+      !
+      glob_sum_full_2d = SUM( ptab(:,:) * tmask_h(:,:) )
+      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_2d )
+      !
+   END FUNCTION glob_sum_full_2d
+
+   FUNCTION glob_sum_full_3d( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_full_3d ***
+      !!
+      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine (nomask)
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_sum_full_3d   ! global sum
+      !!
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   ijpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      ijpk = SIZE(ptab,3)
+      !
+      glob_sum_full_3d = 0.e0
+      DO jk = 1, ijpk
+         glob_sum_full_3d = glob_sum_full_3d + SUM( ptab(:,:,jk) * tmask_h(:,:) )
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( glob_sum_full_3d )
+      !
+   END FUNCTION glob_sum_full_3d
+
 
 #else  
@@ -314,4 +358,65 @@
    END FUNCTION glob_sum_3d_a   
 
+   FUNCTION glob_sum_full_2d( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_full_2d ***
+      !!
+      !! ** Purpose : perform a sum in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:) ::   ptab
+      REAL(wp)                             ::   glob_sum_full_2d   ! global sum (nomask)
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj   ! dummy loop indices
+      !!-----------------------------------------------------------------------
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jj = 1, jpj
+         DO ji =1, jpi
+         ztmp =  ptab(ji,jj) * tmask_h(ji,jj) 
+         CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
+      glob_sum_full_2d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_sum_full_2d
+
+   FUNCTION glob_sum_full_3d( ptab )
+      !!----------------------------------------------------------------------
+      !!                  ***  FUNCTION  glob_sum_full_3d ***
+      !!
+      !! ** Purpose : perform a sum on a 3D array in calling DDPDD routine
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(:,:,:) ::   ptab
+      REAL(wp)                               ::   glob_sum_full_3d   ! global sum (nomask)
+      !!
+      COMPLEX(wp)::   ctmp
+      REAL(wp)   ::   ztmp
+      INTEGER    ::   ji, jj, jk   ! dummy loop indices
+      INTEGER    ::   ijpk ! local variables: size of ptab
+      !!-----------------------------------------------------------------------
+      !
+      ijpk = SIZE(ptab,3)
+      !
+      ztmp = 0.e0
+      ctmp = CMPLX( 0.e0, 0.e0, wp )
+      DO jk = 1, ijpk
+         DO jj = 1, jpj
+            DO ji =1, jpi
+            ztmp =  ptab(ji,jj,jk) * tmask_h(ji,jj)
+            CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+            END DO
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( ctmp )   ! sum over the global domain
+      glob_sum_full_3d = REAL(ctmp,wp)
+      !
+   END FUNCTION glob_sum_full_3d
+
+
+
 #endif
 
