Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/dom_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/dom_oce.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/dom_oce.F90	(revision 12166)
@@ -32,9 +32,7 @@
    LOGICAL , PUBLIC ::   ln_linssh      !: =T  linear free surface ==>> model level are fixed in time
    LOGICAL , PUBLIC ::   ln_meshmask    !: =T  create a mesh-mask file (mesh_mask.nc)
-   REAL(wp), PUBLIC ::   rn_isfhmin     !: threshold to discriminate grounded ice to floating ice
    REAL(wp), PUBLIC ::   rn_rdt         !: time step for the dynamics and tracer
    REAL(wp), PUBLIC ::   rn_atfp        !: asselin time filter parameter
    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
 
@@ -158,4 +156,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   e3t_1d  , e3w_1d   !: reference vertical scale factors at T- and W-pts (m)
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep, bathy
 
    !!----------------------------------------------------------------------
@@ -170,7 +169,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level             (ISF)
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level (ISF)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF)
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: top first wet T-, U-, V-, F-level           (ISF)
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts
@@ -274,6 +271,7 @@
          &                      hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6)  )
          !
-         !
-      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(7) )
+      ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7)  ) 
+         !
+      ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) )
          !
       ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
@@ -281,6 +279,5 @@
          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
          !
-      ALLOCATE( misfdep(jpi,jpj) , mikt(jpi,jpj) , miku(jpi,jpj) ,     &
-         &      risfdep(jpi,jpj) , mikv(jpi,jpj) , mikf(jpi,jpj) , STAT=ierr(10) )
+      ALLOCATE( 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: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domain.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domain.F90	(revision 12166)
@@ -140,11 +140,4 @@
                                        ! Read in masks to define closed seas and lakes 
       !
-      DO jj = 1, jpj                   ! depth of the iceshelves
-         DO ji = 1, jpi
-            ik = mikt(ji,jj)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ik)
-         END DO
-      END DO
-      !
       ht_0(:,:) = 0._wp  ! Reference ocean thickness
       hu_0(:,:) = 0._wp
@@ -192,9 +185,8 @@
       IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
       !
-      IF( ln_meshmask .AND. .NOT.ln_iscpl )                        CALL dom_wri     ! Create a domain file
-      IF( ln_meshmask .AND.      ln_iscpl .AND. .NOT.ln_rstart )   CALL dom_wri     ! Create a domain file
-      IF(                                       .NOT.ln_rstart )   CALL dom_ctl     ! Domain control
-      !
-      IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file
+      IF( ln_meshmask    )   CALL dom_wri       ! Create a domain file
+      IF( .NOT.ln_rstart )   CALL dom_ctl       ! Domain control
+      !
+      IF( ln_write_cfg   )   CALL cfg_write     ! create the configuration file
       !
       IF(lwp) THEN
@@ -292,6 +284,6 @@
          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
-         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios
-      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask
+         &             ln_cfmeta, ln_xios_read, nn_wxios
+      NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask
 #if defined key_netcdf4
       NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
@@ -343,5 +335,4 @@
          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
          IF( TRIM(Agrif_CFixed()) == '0' ) THEN
             WRITE(numout,*) '      READ restart for a single file using XIOS ln_xios_read =', ln_xios_read
@@ -417,5 +408,4 @@
          WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
          WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]'
          WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt
          WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domvvl.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domvvl.F90	(revision 12166)
@@ -36,4 +36,5 @@
 
    PUBLIC  dom_vvl_init       ! called by domain.F90
+   PUBLIC  dom_vvl_zgr        ! called by isfcpl.F90
    PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
    PUBLIC  dom_vvl_sf_swp     ! called by step.F90
@@ -115,21 +116,47 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      !
+      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
+      !
+      !                    ! Allocate module arrays
+      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
+      !
+      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
+      CALL dom_vvl_rst( nit000, 'READ' )
+      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      !
+      CALL dom_vvl_zgr() ! interpolation scale factor, depth and water column
+      !
+   END SUBROUTINE dom_vvl_init
+   !
+   SUBROUTINE dom_vvl_zgr()
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_init  ***
+      !!                   
+      !! ** Purpose :  Interpolation of all scale factors, 
+      !!               depths and water column heights
+      !!
+      !! ** Method  :  - interpolate scale factors
+      !!
+      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
+      !!              - Regrid: e3(u/v)_n
+      !!                        e3(u/v)_b       
+      !!                        e3w_n           
+      !!                        e3(u/v)w_b      
+      !!                        e3(u/v)w_n      
+      !!                        gdept_n, gdepw_n and gde3w_n
+      !!              - h(t/u/v)_0
+      !!              - frq_rst_e3t and frq_rst_hdv
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
+      !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
       REAL(wp)::   zcoef
       !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-      !
-      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
-      !
-      !                    ! Allocate module arrays
-      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
-      !
-      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL dom_vvl_rst( nit000, 'READ' )
-      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
@@ -263,6 +290,5 @@
       ENDIF
       !
-   END SUBROUTINE dom_vvl_init
-
+   END SUBROUTINE dom_vvl_zgr
 
    SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
@@ -813,7 +839,9 @@
             id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
             id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
+            !
             !                             ! --------- !
             !                             ! all cases !
             !                             ! --------- !
+            !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
                CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
@@ -1033,5 +1061,4 @@
       !
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
-      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
       !
       IF(lwp) THEN                   ! Print the choice
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domwri.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domwri.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domwri.F90	(revision 12166)
@@ -16,4 +16,5 @@
    !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
    !!----------------------------------------------------------------------
+   !
    USE dom_oce         ! ocean space and time domain
    USE phycst ,   ONLY :   rsmall
@@ -155,10 +156,8 @@
       
       ! note that mbkt is set to 1 over land ==> use surface tmask
-      zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
+      zprt(:,:) = REAL( mbkt(:,:) , wp )
       CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
-      zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
+      zprt(:,:) = REAL( mikt(:,:) , wp )
       CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
-      zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
       !															             ! vertical mesh
       CALL iom_rstput( 0, 0, inum, 'e3t_1d', e3t_1d, ktype = jp_r8  )    !    ! scale factors
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domzgr.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domzgr.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/domzgr.F90	(revision 12166)
@@ -71,5 +71,6 @@
       INTEGER, DIMENSION(:,:), INTENT(out) ::   k_top, k_bot   ! ocean first and last level indices
       !
-      INTEGER  ::   jk                  ! dummy loop index
+      INTEGER  ::   ji,jj,jk            ! dummy loop index
+      INTEGER  ::   ikt, ikb            ! top/bot index
       INTEGER  ::   ioptio, ibat, ios   ! local integer
       REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
@@ -138,6 +139,15 @@
       !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
       CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1
-      
-
+      !
+      !                                ! ice shelf draft and bathymetry
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = mikt(ji,jj)
+            ikb = mbkt(ji,jj)
+            bathy  (ji,jj) = gdepw_0(ji,jj,ikb+1)
+            risfdep(ji,jj) = gdepw_0(ji,jj,ikt  )
+         END DO
+      END DO
+      !
       !                                ! deepest/shallowest W level Above/Below ~10m
 !!gm BUG in s-coordinate this does not work!
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplhsb.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplhsb.F90	(revision 12146)
+++ 	(revision )
@@ -1,324 +1,0 @@
-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 oce             ! global tra/dyn variable
-   USE dom_oce         ! ocean space and time domain
-   USE domwri          ! ocean space and time domain
-   USE domngb          ! 
-   USE phycst          ! physical constants
-   USE sbc_oce         ! surface boundary condition variables
-   USE iscplini        ! 
-   !
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE lbclnk          !
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_div   
-   PUBLIC   iscpl_cons       
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-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 , zsumn, zjip1_ratio  , zjim1_ratio, zdtem, zde3t, z1_rdtiscpl
-      REAL(wp) ::   zarea , zsum1, zsumb, zjjp1_ratio  , zjjm1_ratio, zdsal
-      REAL(wp), DIMENSION(jpi,jpj)        ::   zdssh   ! workspace
-      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
-      !!----------------------------------------------------------------------
-
-      ! get imbalance (volume heat and salt)
-      ! initialisation difference
-      zde3t = 0._wp   ;   zdsal = 0._wp   ;   zdtem = 0._wp
-
-      ! initialisation correction term
-      pvol_flx(:,:,:  ) = 0._wp
-      pts_flx (:,:,:,:) = 0._wp
-      
-      z1_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 (.NOT. ln_linssh ) 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 = e3t_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) * e3t_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) * e3t_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 * z1_rdtiscpl
-                  pts_flx (ji,jj,jk,jp_sal)=   pts_flx (ji,jj,jk,jp_sal) + zdsal * z1_rdtiscpl 
-                  pts_flx (ji,jj,jk,jp_tem)=   pts_flx (ji,jj,jk,jp_tem) + zdtem * z1_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
-
-!!gm  ERROR !!!!
-!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
-!
-!      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.)
-      CALL ctl_stop( 'STOP', ' iscpl_cons:   please modify this MODULE !' )
-!!gm end
-      ! 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('iscplhsb', 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('iscplhsb', zlat ,npts)
-      CALL mpp_max('iscplhsb', zlon ,npts)
-      CALL mpp_max('iscplhsb', izpts,npts)
-
-      ! correction values 
-      CALL mpp_max('iscplhsb', zcorr_vol,npts)
-      CALL mpp_max('iscplhsb', zcorr_sal,npts)
-      CALL mpp_max('iscplhsb', 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(:,:,:)
-
-!!gm  ERROR !!!!
-!!    juste use tmask_i  or in case of ISF smask_i (to be created to compute the sum without halos)
-!
-!      ! 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)
-!!gm end
-
-      !
-   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) / e3t_n(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE iscpl_div
-
-END MODULE iscplhsb
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplini.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplini.F90	(revision 12146)
+++ 	(revision )
@@ -1,90 +1,0 @@
-MODULE iscplini
-   !!======================================================================
-   !!                       ***  MODULE  sbciscpl  ***
-   !! Ocean forcing:  ?????
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   iscpl_init     : initialisation routine (namelist)
-   !!   iscpl_alloc    : allocation of correction variables
-   !!----------------------------------------------------------------------
-   USE oce             ! global tra/dyn variable
-   USE dom_oce         ! ocean space and time domain
-   !
-   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    !:
-   INTEGER , PUBLIC ::   nn_drown     !:
-   
-   INTEGER , PUBLIC ::   nstp_iscpl   !:
-   REAL(wp), PUBLIC ::   rdt_iscpl    !: 
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   hdiv_iscpl   !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   htsc_iscpl   !:
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   INTEGER FUNCTION iscpl_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE sbc_iscpl_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=iscpl_alloc )
-         !
-      CALL mpp_sum ( 'iscplini', 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
-      !!----------------------------------------------------------------------
-      !
-      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' )
-      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' )
-      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
-      ENDIF
-      !
-   END SUBROUTINE iscpl_init
-
-   !!======================================================================
-END MODULE iscplini
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplrst.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/iscplrst.F90	(revision 12146)
+++ 	(revision )
@@ -1,407 +1,0 @@
-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 oce             ! global tra/dyn variable
-   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 iscplini        ! ice sheet coupling: initialisation
-   USE iscplhsb        ! ice sheet coupling: conservation
-   !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O module
-   USE lib_mpp         ! MPP library
-   USE lib_fortran     ! MPP library
-   USE lbclnk          ! communication
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   iscpl_stp          ! step management 
-   !!
-   !! * Substitutions  
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!----------------------------------------------------------------------
-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(jpi,jpj)     ::   zsmask_b
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztmask_b, zumask_b, zvmask_b
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t_b  , ze3u_b  , ze3v_b  
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepw_b
-      CHARACTER(20) :: cfile
-      !!----------------------------------------------------------------------
-      !
-      !                       ! get restart variable
-      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
-      CALL iom_get( numror, jpdom_autoglo, 'umask'  , zumask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'vmask'  , zvmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'smask'  , zsmask_b, ldxios = lrxios   ) ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:), ldxios = lrxios )  ! need to compute temperature correction
-      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b(:,:,:), ldxios = lrxios )  ! need to correct barotropic velocity
-      CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
-      !
-      CALL iscpl_init()       ! read namelist
-      !                       ! 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 )
-      !
-      IF ( ln_hsb ) THEN      ! compute correction if conservation needed
-         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
-         
-      !                       ! create  a domain file
-      IF( ln_meshmask .AND. ln_iscpl )   CALL dom_wri
-      !
-      IF ( ln_hsb ) THEN
-         cfile='correction'
-         cfile = TRIM( cfile )
-         CALL iom_open  ( cfile, inum0, ldwrt = .TRUE. )
-         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
-      !
-      neuler = 0              ! next step is an euler time step
-      !
-      !                       ! set _b and _n variables equal
-      tsb (:,:,:,:) = tsn (:,:,:,:)
-      ub  (:,:,:)   = un  (:,:,:)
-      vb  (:,:,:)   = vn  (:,:,:)
-      sshb(:,:)     = sshn(:,:)
-      !
-      !                       ! set _b and _n vertical scale factor equal
-      e3t_b (:,:,:) = e3t_n (:,:,:)
-      e3u_b (:,:,:) = e3u_n (:,:,:)
-      e3v_b (:,:,:) = e3v_n (:,:,:)
-      !
-      e3uw_b (:,:,:) = e3uw_n (:,:,:)
-      e3vw_b (:,:,:) = e3vw_n (:,:,:)
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-      hu_b   (:,:)   = hu_n   (:,:)
-      hv_b   (:,:)   = hv_n   (:,:)
-      r1_hu_b(:,:)   = r1_hu_n(:,:)
-      r1_hv_b(:,:)   = r1_hv_n(:,:)
-      !
-   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(jpi,jpj)          :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
-      REAL(wp), DIMENSION(jpi,jpj)          :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn, ztrp
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
-      !!----------------------------------------------------------------------
-      !
-      !                 ! 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
-               ENDIF
-            END DO
-         END DO
-         CALL lbc_lnk_multi( 'iscplrst', sshn, 'T', 1., zsmask1, 'T', 1. )
-         zssh0   = sshn
-         zsmask0 = zsmask1
-      END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
-
-!=============================================================================
-!PM: Is this needed since introduction of VVL by default?
-      IF ( .NOT.ln_linssh ) THEN
-      ! Reconstruction of all vertical scale factors at now time steps
-      ! =============================================================================
-      ! 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
-                     e3t_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( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-
-         ! Vertical scale factor interpolations
-         ! ------------------------------------
-         CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-         CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-         CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-         
-         ! t- and w- points depth
-         ! ----------------------
-         gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-         gdepw_n(:,:,1) = 0.0_wp
-         gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-         DO jj = 1,jpj
-            DO ji = 1,jpi
-               DO jk = 2,mikt(ji,jj)-1
-                  gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
-               END DO
-               IF (mikt(ji,jj) > 1) THEN
-                  jk = mikt(ji,jj)
-                  gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END IF
-               DO jk = mikt(ji,jj)+1, jpk
-                  gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
-                  gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-                  gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk  ) - sshn (ji,jj)
-               END DO
-            END DO
-         END DO
-
-      ! t-, u- and v- water column thickness
-      ! ------------------------------------
-         ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
-         DO jk = 1, jpkm1
-            hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-            hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
-            ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !                                        ! Inverse of the local depth
-         r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
-         r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 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)/e3u_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)/e3v_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(:,:,:)*e3u_n(:,:,:); 
-      zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:); 
-      zbvn(:,:)   = SUM(ztrp,DIM=3)
-
-      ! new water column
-      zhu1=0.0_wp ;
-      zhv1=0.0_wp ;
-      DO jk  = 1,jpk
-        zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-        zhv1(:,:) = zhv1(:,:) + e3v_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_multi( 'iscplrst', tsn(:,:,:,jp_tem), 'T', 1., tsn(:,:,:,jp_sal), 'T', 1., ztmask1, 'T', 1.)
-
-          ! 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
-      !PM: Is this IF needed since change to VVL by default
-      IF (.NOT.ln_linssh) THEN
-         DO jk = 2,jpk-1
-            DO jj = 1,jpj
-               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,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
-                     zdz   =           gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
-                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - gdepw_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)      )/e3t_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)      )/e3t_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
-      ! 
-   END SUBROUTINE iscpl_rst_interpol
-
-   !!======================================================================
-END MODULE iscplrst
Index: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/istate.F90	(revision 12146)
+++ NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DOM/istate.F90	(revision 12166)
@@ -28,5 +28,4 @@
    USE dtauvd         ! data: U & V current             (dta_uvd routine)
    USE domvvl          ! varying vertical mesh
-   USE iscplrst        ! ice sheet coupling
    USE wet_dry         ! wetting and drying (needed for wad_istate)
    USE usrdef_istate   ! User defined initial state
@@ -86,5 +85,4 @@
          !                                    ! -------------------
          CALL rst_read                        ! Read the restart file
-         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !
