Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 9630)
@@ -11,5 +11,5 @@
 ============================================================================================================
     -->
-   <field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" detect_missing_value="true"> <!-- time step automaticaly defined -->
+   <field_definition level="1" prec="4" operation="average" enabled=".TRUE." default_value="1.e20" detect_missing_value="false"> <!-- time step automaticaly defined -->
 
     <!-- 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 9630)
@@ -494,4 +494,11 @@
                           !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate
                           ! 2 = velocity and stability dependent Gamma    Holland et al. 1999
+/
+!-----------------------------------------------------------------------
+&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)
 /
 !-----------------------------------------------------------------------
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 9630)
@@ -43,4 +43,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
 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 9630)
@@ -138,5 +138,5 @@
          &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , ln_rstdate, 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
+         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler, ln_iscpl
       NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   &
          &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  &
@@ -193,4 +193,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
 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 9630)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 9630)
@@ -0,0 +1,327 @@
+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 =   e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) + e12t(ji  ,jjm1) * tmask(ji  ,jjm1,jk) &
+!                       &    + e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) + e12t(jip1,jj  ) * tmask(jip1,jj  ,jk)
+!
+!                     IF ( zsum /= 0._wp ) THEN
+!                        zjip1_ratio   = e12t(jip1,jj  ) * tmask(jip1,jj  ,jk) / zsum
+!                        zjim1_ratio   = e12t(jim1,jj  ) * tmask(jim1,jj  ,jk) / zsum
+!                        zjjp1_ratio   = e12t(ji  ,jjp1) * tmask(ji  ,jjp1,jk) / zsum
+!                        zjjm1_ratio   = e12t(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
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*) 'adding the divergence correction contribution at nit000' 
+
+      !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)
+      !         phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + rhdivdiff(ji,jj,jk)
+      !      END DO
+      !   END DO
+      !END DO
+      !
+      phdivn(:,:,:) =  phdivn(:,:,:) + rhdivdiff(:,:,:)
+
+   END SUBROUTINE iscpl_div
+
+END MODULE iscplhsb
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 9630)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 9630)
@@ -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/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 9630)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 9630)
@@ -0,0 +1,483 @@
+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, jtop
+      !!
+      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb, zhdivdiff, zfp_ui, zfm_ui, zfp_vi, zfm_vi
+      REAL(wp):: zdz, zdzm1, zdzp1, zvol
+      REAL(wp),DIMENSION(2)::ztrdtsb
+      !!
+      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, zuflux, zvflux
+      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, zuflux, zvflux                       ) 
+      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,nn_drown    ! 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 - umask_i(:,:) ) * umask_i(:,:)
+         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+
+      END IF
+
+!=============================================================================
+! compute velocity
+
+      ub(:,:,:)=un(:,:,:)
+      vb(:,:,:)=vn(:,:,:)
+
+      un(:,:,:)=ub(:,:,:)*umask(:,:,:)
+      vn(:,:,:)=vb(:,:,:)*vmask(:,:,:)
+    
+
+      !-----initialise the divergence and tracer correction------
+      rhdivdiff(:,:,:) = 0._wp
+      rhdivdiff_trac(:,:,:,:) = 0._wp      
+
+      !-----the change of flux in each velocity cell 
+      zuflux(:,:,:) = (pe3u_b(:,:,:)*pumask_b(:,:,:) - fse3u_n(:,:,:)*umask(:,:,:))*ub(:,:,:)
+      zvflux(:,:,:) = (pe3v_b(:,:,:)*pvmask_b(:,:,:) - fse3v_n(:,:,:)*vmask(:,:,:))*vb(:,:,:)
+     
+
+      ! compute divergence and tracer correction (whether the horizontal velocity cell is closed or open)    
+      DO jj = 1,jpj
+         DO ji = 1,jpim1   ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
+            !--------u-cell is closed-----------  
+            DO jk  = 1,jpk
+                 !------the contribution to the west side of u-cell
+                 IF(tmask(ji,jj,jk) == 1) THEN   ! the west cell is still active 
+                    zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
+                    rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol 
+                    rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zuflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e2u(ji,jj)/zvol
+                 ELSEIF(mikt(ji,jj) > 1) THEN  ! the west cell becomes dry, put the waterflux to the top of the west column
+                    jtop=mikt(ji,jj)
+                    zvol=e1e2t(ji,jj)*fse3t(ji,jj,jtop)
+                    rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
+                    rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2)-zuflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e2u(ji,jj)/zvol
+                 !ELSE ! for the closing T-column : nothing to be done 
+                 ENDIF
+
+                 !------the contribution to the east side of u-cell
+                 IF(tmask(ji+1,jj,jk) == 1) THEN  ! the east cell is still active
+                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jk)
+                    rhdivdiff(ji+1,jj,jk) = rhdivdiff(ji+1,jj,jk) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
+                    rhdivdiff_trac(ji+1,jj,jk,1:2) = rhdivdiff_trac(ji+1,jj,jk,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jk,1:2)*e2u(ji,jj)/zvol
+                 ELSEIF(mikt(ji+1,jj) > 1) THEN    ! the east cell becomes dry, put the waterflux to the top of the east column
+                    jtop = mikt(ji+1,jj)
+                    zvol = e1e2t(ji+1,jj)*fse3t(ji+1,jj,jtop)
+                    rhdivdiff(ji+1,jj,jtop) = rhdivdiff(ji+1,jj,jtop) - zuflux(ji,jj,jk)*e2u(ji,jj)/zvol
+                    rhdivdiff_trac(ji+1,jj,jtop,1:2) = rhdivdiff_trac(ji+1,jj,jtop,1:2) + zuflux(ji,jj,jk)*tsb(ji+1,jj,jtop,1:2)*e2u(ji,jj)/zvol
+                 !ELSE  ! for the closing T-column : nothing to be done
+                 ENDIF
+            ENDDO
+         END DO
+      END DO
+
+
+      DO jj = 1,jpjm1  ! jpim1 rather than jpi for the correct lbc_lnk synchronisation
+         DO ji = 1,jpi   
+            !-----------v-cell is closed-------------
+            DO jk  = 1,jpk
+                 !------the contribution to the south side of v-cell
+                 IF(tmask(ji,jj,jk) == 1) THEN  ! the south cell is still active
+                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jk)
+                   rhdivdiff(ji,jj,jk) = rhdivdiff(ji,jj,jk) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
+                   rhdivdiff_trac(ji,jj,jk,1:2) =  rhdivdiff_trac(ji,jj,jk,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jk,1:2)*e1v(ji,jj)/zvol
+                 ELSEIF(mikt(ji,jj) > 1) THEN  !the south cell becomes dry, put the waterflux to the top of the south column
+                   jtop = mikt(ji,jj)
+                   zvol = e1e2t(ji,jj)*fse3t(ji,jj,jtop)
+                   rhdivdiff(ji,jj,jtop) = rhdivdiff(ji,jj,jtop) + zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
+                   rhdivdiff_trac(ji,jj,jtop,1:2) = rhdivdiff_trac(ji,jj,jtop,1:2) - zvflux(ji,jj,jk)*tsb(ji,jj,jtop,1:2)*e1v(ji,jj)/zvol
+                 ENDIF
+
+                 !------the contribution to the north side of v-cell
+                 IF(tmask(ji,jj+1,jk) == 1) THEN ! the north cell is still active
+                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jk)
+                    rhdivdiff(ji,jj+1,jk) = rhdivdiff(ji,jj+1,jk) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
+                    rhdivdiff_trac(ji,jj+1,jk,1:2) = rhdivdiff_trac(ji,jj+1,jk,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jk,1:2)*e1v(ji,jj)/zvol
+                 ELSEIF(mikt(ji,jj+1) > 1) THEN  ! the north cell becomes dry, put the waterflux to the top of the north column
+                    jtop = mikt(ji,jj+1)
+                    zvol = e1e2t(ji,jj+1)*fse3t(ji,jj+1,jtop)
+                    rhdivdiff(ji,jj+1,jtop) = rhdivdiff(ji,jj+1,jtop) - zvflux(ji,jj,jk)*e1v(ji,jj)/zvol
+                    rhdivdiff_trac(ji,jj+1,jtop,1:2) = rhdivdiff_trac(ji,jj+1,jtop,1:2) + zvflux(ji,jj,jk)*tsb(ji,jj+1,jtop,1:2)*e1v(ji,jj)/zvol
+                 ENDIF
+            ENDDO
+
+         END DO
+      END DO 
+      
+
+      CALL lbc_lnk(rhdivdiff,   'T',1._wp)
+      CALL lbc_lnk(rhdivdiff_trac(:,:,:,1),   'T',1._wp)
+      CALL lbc_lnk(rhdivdiff_trac(:,:,:,2),   'T',1._wp)
+
+
+      ! 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,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
+!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
+!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- 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
+      ! 
+      ! 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, zuflux, zvflux    ) 
+      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/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 9630)
@@ -37,4 +37,5 @@
    USE eosbn2          ! equation of state            (eos bn2 routine)
    USE domvvl          ! varying vertical mesh
+   USE iscplrst        ! ice sheet coupling
    USE dynspg_oce      ! pressure gradient schemes
    USE dynspg_flt      ! filtered free surface
@@ -91,4 +92,5 @@
          !                                    ! -------------------
          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
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 9630)
@@ -35,4 +35,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iscplhsb
 
    IMPLICIT NONE
@@ -327,4 +328,8 @@
       !                                                ! ===============
 
+      IF( ln_iscpl .AND. kt == nit000)   THEN 
+         IF(lwp) WRITE(numout,*) 'for nit000, update hdivn after calling div_cur with the divergence correction'
+                          CALL iscpl_div( hdivn ) !
+      ENDIF
       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field)
       IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 9630)
@@ -169,4 +169,15 @@
                      ENDIF
                      ENDIF
+                  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
+
                      IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
 
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90	(revision 9630)
@@ -220,4 +220,12 @@
       !
       !
+
+      !---start the tendency with tracer source/sink from divergence correction
+      IF ( ln_iscpl .AND. kt == nit000) THEN        
+         IF (lwp) WRITE(numout,*) 'add tendency tsa due to re-meshing at kt = nit000'
+         tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) + rhdivdiff_trac(:,:,:,jp_tem)*tmask(:,:,:)
+         tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) + rhdivdiff_trac(:,:,:,jp_sal)*tmask(:,:,:)
+      ENDIF
+
       !----------------------------------------
       !       Ice Shelf effects (ISF)
Index: /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 9629)
+++ /branches/UKMO/dev_isf_divg_corr_GO6_package_r9385/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 9630)
@@ -71,4 +71,8 @@
    !! Energy budget of the leads (open water embedded in sea ice)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraqsr_1lev        !: fraction of solar net radiation absorbed in the first ocean level [-]
+
+  !!------array used for divergence correction
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  rhdivdiff
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  rhdivdiff_trac
  
    !! Arrays used in coupling when MEDUSA is present. These arrays need to be declared
@@ -129,5 +133,5 @@
          &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     &
          &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     &
-         &     riceload(jpi,jpj),                             STAT=ierr(2) )
+         &     riceload(jpi,jpj),  rhdivdiff(jpi,jpj,jpk), rhdivdiff_trac(jpi,jpj,jpk,2) ,  STAT=ierr(2) )
          !
       ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) )
