Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplhsb.F90	(revision 5790)
@@ -0,0 +1,346 @@
+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"  
+   !!----------------------------------------------------------------------
+   !! 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_tiscpl
+      REAL(wp):: zjip1_ratio, zjim1_ratio, zjjp1_ratio, zjjm1_ratio
+      !!
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zde3t
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmp3d
+      !
+      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zlon, zlat
+      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
+      INTEGER , DIMENSION(:    ), ALLOCATABLE :: ixpts, iypts, izpts, vnpts
+      INTEGER :: jpts, npts
+
+      CALL wrk_alloc(jpi,jpj,jpk,   ztmp3d )
+      CALL wrk_alloc(jpi,jpj,       zde3t  )
+      CALL wrk_alloc(jpi,jpj,       zssh0  )
+
+    ! get unbalance (volume heat and salt)
+    ! initialisation
+       zde3t   (:,:)     = 0.0_wp
+       pvol_flx(:,:,:  ) = 0.0_wp
+       pts_flx (:,:,:,:) = 0.0_wp
+
+       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
+
+       ! mask tsn and tsb (should be useless)
+       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 non conservation of heat, salt and volume 
+       r1_tiscpl = 1._wp / (prdt_iscpl * rn_rdt) 
+       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+       IF ( lk_vvl ) zssh0 = 0.0_wp
+       DO jk = 1,jpk-1
+          DO ji = 2,jpi-1
+             DO jj = 2,jpj-1
+                ! volume differences
+                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk);
+                
+                ! shh changes
+                IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 
+                   zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj)
+                   zssh0(ji,jj) = 0._wp
+                END IF
+
+                ! ocean cell now
+                ! case where we open, enlarge or thin a cell :
+                pvol_flx(ji,jj,jk)       =                          zde3t(ji,jj) * r1_tiscpl
+                pts_flx (ji,jj,jk,jp_sal)=   tsn(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
+                pts_flx (ji,jj,jk,jp_tem)=   tsn(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
+             END DO
+          END DO
+       END DO
+       ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
+       PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
+       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
+
+       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
+       IF ( lk_vvl ) zssh0 = 0.0_wp
+       DO jk = 1,jpk-1
+          DO ji = 2,jpi-1
+             DO jj = 2,jpj-1
+                ! volume differences
+                zde3t(ji,jj) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) - pe3t_b(ji,jj,jk) * ptmask_b(ji,jj,jk);
+                
+                ! shh changes
+                IF ( ptmask_b(ji,jj,jk) == 1 .OR. tmask(ji,jj,jk) == 1 ) THEN 
+                   zde3t(ji,jj) = zde3t(ji,jj) + zssh0(ji,jj)
+                   zssh0(ji,jj) = 0._wp
+                END IF
+
+                ! ocean cell before and mask cell now
+                IF ( tmask(ji,jj,jk) == 0._wp .AND. ptmask_b(ji,jj,jk) == 1._wp ) THEN
+                   ! case where we close a cell and adjacent cell open
+                   pvol_flx(ji,jj,jk)       =                         zde3t(ji,jj) * r1_tiscpl
+                   pts_flx (ji,jj,jk,jp_sal)=  tsb(ji,jj,jk,jp_sal) * zde3t(ji,jj) * r1_tiscpl 
+                   pts_flx (ji,jj,jk,jp_tem)=  tsb(ji,jj,jk,jp_tem) * zde3t(ji,jj) * r1_tiscpl
+
+                   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 .NE. 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 .EQ. 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) .EQ. 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
+                         PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
+                         PRINT *,        '                     ',mig(ji),' ',mjg(jj),' ',jk
+                         PRINT *,        '                     ',ji,' ',jj,' ',jk,' ',narea
+                         PRINT *,        ' we are now looking for the closest wet cell on the horizontal '
+                      ! We deal with this points later.
+                      END IF
+                   END IF
+                END IF
+             END DO
+          END DO
+       END DO
+
+       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
+       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+       IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
+
+       ! allocation and initialisation of the list of problematic point
+       ALLOCATE(vnpts(jpnij))
+       vnpts(:)=0
+
+       ! fill narea location with the number of problematic point
+       DO jk = 1,jpk-1
+          DO ji = 2,jpi-1
+             DO jj = 2,jpj-1
+                IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
+                &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN
+                   vnpts(narea) = vnpts(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(vnpts,jpnij) 
+
+       ! size of the new variable
+       npts  = SUM(vnpts)    
+       
+       ! 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 ; zlat(:) = -1.0e20
+       zcorr_vol(:) = 0.0_wp
+       zcorr_sal(:) = 0.0_wp
+       zcorr_tem(:) = 0.0_wp
+
+       ! fill new variable
+       jpts = SUM(vnpts(1:narea-1))
+       DO jk = 1,jpk-1
+          DO ji = 2,jpi-1
+             DO jj = 2,jpj-1
+                IF (          ptmask_b(ji,jj,jk)      == 1 .AND. SUM(tmask(ji-1:ji+1,jj,jk)) == 0 &
+                &   .AND. SUM(tmask(ji,jj-1:jj+1,jk)) == 0 .AND.     tmask(ji,jj,jk+1)       == 0 ) THEN
+                   jpts = jpts + 1  ! positioning in the vnpts vector for the area narea
+                   PRINT *, 'corrected point ', narea, ji, jj, jk, jpts
+                   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
+                   PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt
+                END IF
+             END DO
+          END DO
+       END DO
+
+       ! build array of total problematic point on each cpu (share to each cpu)
+       CALL mpp_max(zlat ,npts)
+       CALL mpp_max(zlon ,npts)
+       CALL mpp_max(izpts,npts)
+
+       ! put correction term in the closest cell          
+       PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
+       DO jpts = 1,npts
+          CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
+          PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
+          DO jj = mj0(iypts(jpts)),mj1(iypts(jpts))
+             DO ji = mi0(ixpts(jpts)),mi1(ixpts(jpts))
+                jk = izpts(jpts)
+                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)
+             END DO
+          END DO
+       END DO
+       ! deallocate variables 
+       DEALLOCATE(vnpts)
+       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_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
+       pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
+
+       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.)
+
+       ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
+       zsumn = glob_sum     ( fse3t_n(:,:,:) * tmask  (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
+       ztmp3d(:,:,:) = 0.0
+       ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+       zsumb = glob_sum_full(ztmp3d) 
+       zsum  = glob_sum     ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
+       IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
+       ! CHECK salt
+       zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
+       zsumb = glob_sum( tsb(:,:,:,jp_sal) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
+       zsum  = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt)
+       IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
+       ! CHECK heat
+       zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
+       zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
+       zsum  = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt)
+       IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
+    !! 
+    CALL wrk_dealloc(jpi,jpj,jpk,   ztmp3d ) 
+    CALL wrk_dealloc(jpi,jpj,       zde3t  ) 
+    CALL wrk_dealloc(jpi,jpj,       zssh0  ) 
+   END SUBROUTINE iscpl_cons
+
+   SUBROUTINE iscpl_div( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE iscpl_div  ***
+      !!
+      !! ** Purpose :   update the horizontal divergenc
+      !!
+      !! ** Method  :
+      !!                CAUTION : iscpl is positive (inflow) and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   increase by the iscpl correction term
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+   END SUBROUTINE iscpl_div
+
+END MODULE iscplhsb
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplini.F90	(revision 5790)
@@ -0,0 +1,82 @@
+MODULE iscplini
+   !!======================================================================
+   !!                       ***  MODULE  sbciscpl***
+   !! Ocean forcing:  river runoff
+   !!=====================================================================
+   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   iscpl_rst          : restart correction    in case of coupling with ice sheet
+   !!   iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE oce             ! global tra/dyn variable
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! MPP library
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   iscpl_init      
+   PUBLIC   iscpl_alloc 
+   !!                                                     !!* namsbc_iscpl namelist *
+   LOGICAL , PUBLIC                                        ::   ln_hsb
+   REAL(wp), PUBLIC                                        ::   rn_fiscpl
+   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/rn_fiscpl,ln_hsb
+      !!----------------------------------------------------------------------
+      !                                   ! ============
+      !                                   !   Namelist
+      !                                   ! ============
+      !
+      rn_fiscpl = 2480.
+      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 )
+      !
+      rdt_iscpl=MAX(rn_fiscpl, nitend-nit000+1.0) ! the coupling period have to be less or egal than the total number of time step
+      !
+      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) = ', rdt_iscpl
+      END IF
+
+   END SUBROUTINE iscpl_init
+
+END MODULE iscplini
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 5790)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90	(revision 5790)
@@ -0,0 +1,488 @@
+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_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
+   USE domngb
+   USE sbc_ice, ONLY : lk_lim3
+   USE iscplini
+   USE iscplhsb
+
+   IMPLICIT NONE
+   PRIVATE
+   
+   PUBLIC   iscpl_stp       
+   PUBLIC   iscpl_rst_interpol ! routine to wet and dry
+   !!
+   !! * 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
+
+   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
+      !! 
+      !!----------------------------------------------------------------------
+      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
+!!----------------------------------------------------------------------
+      INTEGER  ::   inum0
+      CHARACTER(20) :: cfile
+
+      CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variables for evoving isf cavity
+      CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! need of this variables for interpolation)
+      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 compute volume      correction  ????
+      CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume      correction  ????
+      CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume      correction  ????
+
+      !! 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_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) )
+         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                   )
+ 
+      neuler = 0              ! Euler restart (neuler=0)
+
+      tsb    (:,:,:,:) = tsn    (:,:,:,:)  ! all before fields set to now values
+      ub     (:,:,:)   = un     (:,:,:)
+      vb     (:,:,:)   = vn     (:,:,:)
+      rotb   (:,:,:)   = rotn   (:,:,:)
+      hdivb  (:,:,:)   = hdivn  (:,:,:)
+      fse3t_b(:,:,:)   = fse3t_n(:,:,:)
+      sshb   (:,:)     = sshn   (:,:)
+      !
+   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
+      REAL(wp), DIMENSION(:,:    ), INTENT(in ) :: psmask_b                        !! mask before
+      !!
+      INTEGER :: ji, jj, jk, iz          !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
+      !!
+      REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
+      REAL(wp):: zdz, zdzm1, zdzp1
+      !!
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, hu1, hv1
+      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
+      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
+      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
+      !!
+      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
+      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d         ) 
+      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
+      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
+      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
+      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
+      !! 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 ssh     
+      zsmask0(:,:) = psmask_b(:,:)
+      zsmask1(:,:) = psmask_b(:,:)
+      
+      ! compute new ssh if we open a full water column (average of the closest neigbourgs)  
+      sshb (:,:)=sshn(:,:)
+      zssh0(:,:)=sshn(:,:)
+      DO iz = 1,10    ! need to be tuned (configuration dependent)
+         zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
+         DO ji = 2,jpi-1
+            DO jj = 2,jpj-1
+               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 .AND. summsk .NE. 0) THEN
+                  PRINT *, 'add ssh to 1 cell',ji,jj,narea
+                  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
+               END IF
+            END DO
+         END DO
+         CALL lbc_lnk(sshn,'T',1.)
+         CALL lbc_lnk(zsmask1,'T',1.)
+         zssh0   = sshn
+         zsmask0 = zsmask1
+      END DO
+      sshn(:,:) = sshn(:,:) * ssmask(:,:)
+
+!=============================================================================
+! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full
+      ztmp3d(:,:,:) = 0.0
+      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+      zsum = glob_sum_full(ztmp3d) 
+      IF (lwp) PRINT *, 'total volume correction 00 = ',zsum
+      IF ( lk_vvl ) THEN
+! compute fse3t_n
+         DO jk = 1,jpk
+            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) )
+         END DO
+         ztmp3d(:,:,:) = 0.0
+         ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+         &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+         zsum = glob_sum_full(ztmp3d) 
+         IF (lwp) PRINT *, 'total volume correction 01 = ',zsum
+
+! compute fse3u/v ... (call interpolation vvl)
+      ! Reconstruction of all vertical scale factors at now and before time steps
+      ! =============================================================================
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+         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) .GT. 1) THEN
+                  jk = mikt(ji,jj)
+                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
+                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
+                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+               END IF
+               DO jk = mikt(ji,jj)+1, jpk
+                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
+                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
+               END DO
+            END DO
+         END DO
+
+         hu(:,:) = 0._wp                          ! Ocean depth at U-points
+         hv(:,:) = 0._wp                          ! Ocean depth at V-points
+         ht(:,:) = 0._wp                          ! Ocean depth at T-points
+         DO jk = 1, jpkm1
+            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+         END DO
+         !                                        ! Inverse of the local depth
+         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
+         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
+
+      END IF
+      ztmp3d(:,:,:) = 0.0
+      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+      zsum = glob_sum_full(ztmp3d) 
+      IF (lwp) PRINT *, 'total volume correction 02 = ',zsum
+
+!=============================================================================
+
+! 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 ji = 1,jpi
+            DO jj = 1,jpj
+               un(ji,jj,jk) = ub(ji,jj,jk)*pe3u_b(ji,jj,jk)*pumask_b(ji,jj,jk)/fse3u_n(ji,jj,jk)*umask(ji,jj,jk);
+               vn(ji,jj,jk) = vb(ji,jj,jk)*pe3v_b(ji,jj,jk)*pvmask_b(ji,jj,jk)/fse3v_n(ji,jj,jk)*vmask(ji,jj,jk);
+            END DO
+         END DO
+      END DO
+! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
+! compute barotropic velocity now and after 
+      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
+      zbub(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
+      zbvb(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
+      zbun(:,:)   = SUM(ztrp,DIM=3)
+      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
+      zbvn(:,:)   = SUM(ztrp,DIM=3)
+! Already know ????????
+      hu1=0.0_wp ;
+      hv1=0.0_wp ;
+      DO jk  = 1,jpk
+        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+      END DO
+      zucorr = 0._wp
+      zvcorr = 0._wp
+      DO ji = 1,jpi
+         DO jj = 1,jpj
+            IF (zbun(ji,jj) .NE. zbub(ji,jj) .AND. hu1(ji,jj) .NE. 0._wp ) THEN
+               zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/hu1(ji,jj)
+            END IF
+            IF (zbvn(ji,jj) .NE. zbvb(ji,jj) .AND. hv1(ji,jj) .NE. 0._wp ) THEN
+               zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/hv1(ji,jj)
+            END IF
+         END DO
+      END DO 
+      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,10
+          DO jk = 1,jpk-1
+              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
+              DO ji = 2,jpi-1
+                  DO jj = 2,jpj-1
+                      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))
+                      !! horizontal extrapolation
+                      IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
+                         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
+                      END IF
+                      !! vertical extrapolation if horizontal extra failed
+                      IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN
+                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+                         summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
+                         IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0 ) 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
+                         END IF
+                      END IF
+                      !! horizontal corner extrapolation if horizontal and vertical failed
+                      IF (zdmask(ji,jj)==1 .AND. summsk==0) THEN
+                         jip1=ji+1; jim1=ji-1;
+                         jjp1=jj+1; jjm1=jj-1;
+                         summsk=(ztmask0(jip1,jjp1,jk)+ztmask0(jim1,jjm1,jk)+ztmask0(jim1,jjp1,jk)+ztmask0(jip1,jjm1,jk))
+                         IF (zdmask(ji,jj)==1 .AND. summsk .NE. 0) THEN
+                            tsn(ji,jj,jk,1)=( zts0(jip1,jjp1,jk,1)*ztmask0(jip1,jjp1,jk)     &
+                            &                +zts0(jim1,jjm1,jk,1)*ztmask0(jim1,jjm1,jk)     &
+                            &                +zts0(jim1,jjp1,jk,1)*ztmask0(jim1,jjp1,jk)     &
+                            &                +zts0(jip1,jjm1,jk,1)*ztmask0(jip1,jjm1,jk))/summsk
+                            tsn(ji,jj,jk,2)=( zts0(jip1,jjp1,jk,2)*ztmask0(jip1,jjp1,jk)     &
+                            &                +zts0(jim1,jjm1,jk,2)*ztmask0(jim1,jjm1,jk)     &
+                            &                +zts0(jim1,jjp1,jk,2)*ztmask0(jim1,jjp1,jk)     &
+                            &                +zts0(jip1,jjm1,jk,2)*ztmask0(jip1,jjm1,jk))/summsk
+                            ztmask1(ji,jj,jk)=1
+                         END IF
+                      END IF
+                  END DO
+              END DO
+          END DO
+          CALL lbc_lnk(tsn(:,:,:,1),'T',1.)
+          CALL lbc_lnk(tsn(:,:,:,2),'T',1.)
+          CALL lbc_lnk(ztmask1,'T',1.)
+          zts0(:,:,:,:) = tsn(:,:,:,:)
+          ztmask0 = ztmask1
+      END DO
+      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
+      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
+
+      ! CHECK heat
+      zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
+      zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
+      IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb
+
+      ! compute new T/S (interpolation) because of vvl
+      IF ( lk_vvl .AND. .FALSE. ) THEN
+      !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.0_wp ) THEN
+                     !compute weight
+                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
+                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
+                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  ))
+                     IF (zdz .LT. 0.0_wp) THEN 
+                        PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
+                        PRINT *, 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
+                        PRINT *, 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
+                        PRINT *, 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
+                        PRINT *, 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
+                        PRINT *, 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
+                        PRINT *, 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
+                        PRINT *, 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
+                        PRINT *, 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
+                        PRINT *, 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
+                        STOP
+                     END IF
+                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
+                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
+                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
+                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
+                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
+                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
+                  END IF
+               END DO
+            END DO
+         END DO               
+      END IF
+
+      ! CHECK heat
+      ztmp3d(:,:,:) = 0.0
+      ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
+      &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
+      zsum = glob_sum_full(ztmp3d) 
+      IF (lwp) PRINT *, 'total volume correction 03 = ',zsum
+
+      WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp) 
+         tsn(:,:,:,2)=  -99._wp
+         tmask(:,:,:) = 0.0
+         umask(:,:,:) = 0.0
+         vmask(:,:,:) = 0.0
+      END WHERE
+
+      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
+
+      !! 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
+      !! 
+      CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
+      CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
+      CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
+      CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
+      CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
+      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
+      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , hu1 , hv1            ) 
+
+   END SUBROUTINE iscpl_rst_interpol
+
+END MODULE iscplrst
Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbciscpl.F90	(revision 5779)
+++ 	(revision )
@@ -1,811 +1,0 @@
-MODULE sbc_iscpl
-   !!======================================================================
-   !!                       ***  MODULE  sbciscpl***
-   !! Ocean forcing:  river runoff
-   !!=====================================================================
-   !! History :  NEMO  ! 2015-01 P. Mathiot: original 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   sbc_iscpl_alloc: variable allocation
-   !!   rst_iscpl      : restart correction in case of coupling with ice sheet
-   !!   sbc_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 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
-   USE domngb
-   USE iom
-   USE sbc_ice, ONLY : lk_lim3
-
-   IMPLICIT NONE
-   PRIVATE
-   
-   PUBLIC   sbc_iscpl_div   
-   PUBLIC   rst_iscpl       
-   PUBLIC   rst_iscpl_interpol ! routine to wet and dry
-   PUBLIC   rst_iscpl_alloc 
-   !!                                                     !!* namsbc_iscpl namelist *
-   LOGICAL , PUBLIC                                        ::   ln_hfb
-   REAL(wp), PUBLIC                                        ::   rn_fiscpl
-   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 rst_iscpl_alloc()
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE sbc_iscpl_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( htsc_iscpl(jpi,jpj,jpk,jpts) , hdiv_iscpl(jpi,jpj,jpk) , STAT=rst_iscpl_alloc )
-         !
-      IF( lk_mpp              )   CALL mpp_sum ( rst_iscpl_alloc )
-      IF( rst_iscpl_alloc > 0 )   CALL ctl_warn('rst_iscpl_alloc: allocation of arrays failed')
-   END FUNCTION rst_iscpl_alloc
- 
-   SUBROUTINE rst_iscpl
-      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
-      INTEGER  ::   jk
-      LOGICAL  ::   llok
-!!----------------------------------------------------------------------
-      INTEGER  ::   inum0
-      CHARACTER(20) :: cfile
-
-
-      WRITE(numout,*) 'rst_read ln_iscpl : ', ln_iscpl
-      IF ( ln_iscpl ) THEN
-         CALL wrk_alloc(jpi,jpj,jpk, ztmask_b, zumask_b, zvmask_b) ! need this new variable for evoving isf cavity
-         CALL wrk_alloc(jpi,jpj,jpk, ze3t_b  , ze3u_b  , ze3v_b  ) ! need of this variable afterwards (only used for interpolation)
-         CALL wrk_alloc(jpi,jpj,jpk, zdepw_b )
-         CALL wrk_alloc(jpi,jpj,     zsmask_b                    )
-
-         IF( rst_iscpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
-
-         WRITE(numout,*) 'rst_read : modify restart to fit new geometry'
-         WRITE(numout,*) '~~~~~~~~'
-
-         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 compute volume      correction  ????
-         CALL iom_get( numror, jpdom_autoglo, 'fse3v_n', ze3v_b(:,:,:) ) ! need to compute volume      correction  ????
-         CALL iom_get( numror, jpdom_autoglo, 'fsdepw_n', zdepw_b(:,:,:) ) ! need to compute volume      correction  ????
-         !!  ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
-         CALL rst_iscpl_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b, htsc_iscpl, hdiv_iscpl )
-         IF( nmsh /= 0 .AND. ln_iscpl )   CALL dom_wri      ! Create a domain file
-
-         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_rstput( 0, 0, inum0, 'e3t_diff', (e3t_0(:,:,:)*tmask(:,:,:) - ze3t_b(:,:,:)*ztmask_b(:,:,:))*tmask(:,:,:) )
-         CALL iom_close ( inum0 )
-
-         CALL wrk_dealloc(jpi,jpj,jpk, ztmask_b,zumask_b,zvmask_b ) ! no need of this variable afterwards (only used for interpolation)
-         CALL wrk_dealloc(jpi,jpj,jpk, ze3t_b  ,ze3u_b  ,ze3v_b   ) ! no need of this variable afterwards (only used for interpolation)
-         CALL wrk_dealloc(jpi,jpj,jpk, zdepw_b )
-         CALL wrk_dealloc(jpi,jpj,     zsmask_b                     )
-      END IF
-      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
-         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values
-         ub   (:,:,:)   = un   (:,:,:)
-         vb   (:,:,:)   = vn   (:,:,:)
-         rotb (:,:,:)   = rotn (:,:,:)
-         hdivb(:,:,:)   = hdivn(:,:,:)
-         sshb (:,:)     = sshn (:,:)
-
-         IF( lk_vvl ) THEN
-            DO jk = 1, jpk
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            END DO
-         ENDIF
-
-         IF( lk_lim3 .AND. .NOT. lk_vvl ) THEN
-            DO jk = 1, jpk
-               fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
-            END DO
-         ENDIF
-
-      ENDIF
-      !
-      IF( lk_lim3 ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev' , fraqsr_1lev )
-      ENDIF
-
-   END SUBROUTINE rst_iscpl
-
-   SUBROUTINE rst_iscpl_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b, pts_flx, pvol_flx)
-      !!---------------------------------------------------------------------- 
-      !!                   ***  ROUTINE rst_iscpl  ***
-      !! 
-      !! ** 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
-      !!----------------------------------------------------------------------
-      !!######## ~ cp/paste pyton script for this routine
-      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
-      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
-      !!
-      INTEGER :: ji, jj, ii, ij, jk, iz          !! loop index
-      INTEGER :: iip1, iim1, ijp1, ijm1, jkp1, jkm1
-      INTEGER ::   ierror, inum  ! temporary integer
-      INTEGER ::   ios           ! Local integer output status for namelist read
-      !!
-      REAL(wp):: summsk, r1_tiscpl, zsum, zsum1, zarea, zsumn, zsumb
-      REAL(wp):: ziip1_ratio, ziim1_ratio, zijp1_ratio, zijm1_ratio
-      REAL(wp):: zdz, zdzm1, zdzp1
-      !!
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zdmask , zdsmask, zvcorr, zucorr, zde3t
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zbub   , zbvb   , zbun  , zbvn
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zssh0  , zssh1, hu1, hv1
-      REAL(wp), DIMENSION(:,:    ), POINTER :: zsmask0, zsmask1
-      REAL(wp), DIMENSION(:,:,:  ), POINTER :: ztmask0, ztmask1, ztrp
-      REAL(wp), DIMENSION(:,:,:  ), POINTER :: zwmaskn, zwmaskb, ztmp3d
-      REAL(wp), DIMENSION(:,:,:,:), POINTER :: zts0
-      !
-      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zlon, zlat 
-      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zcorr_vol, zcorr_tem, zcorr_sal
-      INTEGER , DIMENSION(:    ), ALLOCATABLE :: ixpts, iypts, izpts, vnpts
-      INTEGER :: jpts, npts
-      !
-      NAMELIST/namsbc_iscpl/rn_fiscpl,ln_hfb 
-      !!----------------------------------------------------------------------
-      !                                   ! ============
-      !                                   !   Namelist
-      !                                   ! ============
-      !
-      rn_fiscpl = 2480
-      ln_hfb    = .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 )
-      !
-      PRINT *, 'rn_iscpl .. =',rn_fiscpl, ln_hfb
-
-      CALL wrk_alloc(jpi,jpj,jpk,2, zts0                                   )
-      CALL wrk_alloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp, ztmp3d         ) 
-      CALL wrk_alloc(jpi,jpj,jpk,   zwmaskn, zwmaskb                       ) 
-      CALL wrk_alloc(jpi,jpj,       zsmask0, zsmask1                       ) 
-      CALL wrk_alloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
-      CALL wrk_alloc(jpi,jpj,       zbub   , zbvb    , zbun , zbvn         ) 
-      CALL wrk_alloc(jpi,jpj,       zssh0  , zssh1, hu1, hv1               ) 
-! 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 ssh     
-      zsmask0(:,:) = psmask_b(:,:)
-      zsmask1(:,:) = psmask_b(:,:)
-      
-      !    compute new ssh if we open a full water column (mean of the closest neigbourgs)  
-      sshb (:,:)=sshn(:,:)
-      zssh0(:,:)=sshn(:,:)
-      DO iz = 1,10    ! need to be tuned (configuration dependent)
-          zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
-          DO ii = 2,jpi-1
-              DO ij = 2,jpj-1
-                  iip1=ii+1; iim1=ii-1;
-                  ijp1=ij+1; ijm1=ij-1;
-                  summsk=(zsmask0(iip1,ij)+zsmask0(iim1,ij)+zsmask0(ii,ijp1)+zsmask0(ii,ijm1))
-                  IF (zdsmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
-                      PRINT *, 'add ssh to 1 cell',ii,ij,narea
-                      sshn(ii,ij)=(zssh0(iip1,ij)*zsmask0(iip1,ij)     &
-                      &           +zssh0(iim1,ij)*zsmask0(iim1,ij)     &
-                      &           +zssh0(ii,ijp1)*zsmask0(ii,ijp1)     &
-                      &           +zssh0(ii,ijm1)*zsmask0(ii,ijm1))/summsk
-                      zsmask1(ii,ij)=1
-                  END IF
-              END DO
-          END DO
-          CALL lbc_lnk(sshn,'T',1.)
-          CALL lbc_lnk(zsmask1,'T',1.)
-          zssh0   = sshn
-          zsmask0 = zsmask1
-      END DO
-      sshn(:,:) = sshn(:,:) * ssmask(:,:)
-
-!=============================================================================
-! WARNING we cannot used glob_sum for before time variable (because glob_sum use tmask_i). Need to mask the halo and glob_sum_full
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsum = glob_sum_full(ztmp3d) 
-       IF (lwp) PRINT *, 'total volume correction 00 = ',zsum
-      IF ( lk_vvl ) THEN
-! compute fse3t_n
-         DO jk = 1,jpk
-            fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1 + sshn(:,:) / ( ht_0(:,:) + 1. - ssmask(:,:) ) * tmask(:,:,jk) )
-         END DO
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsum = glob_sum_full(ztmp3d) 
-       IF (lwp) PRINT *, 'total volume correction 01 = ',zsum
-
-! compute fse3u/v ... (call interpolation vvl)
-      ! Reconstruction of all vertical scale factors at now and before time steps
-      ! =============================================================================
-      ! Horizontal scale factor interpolations
-      ! --------------------------------------
-         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) .GT. 1) THEN
-                  jk = mikt(ji,jj)
-                  fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk)
-                  fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
-                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END IF
-               DO jk = mikt(ji,jj)+1, jpk
-                  fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)
-                  fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
-                  fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj)
-               END DO
-            END DO
-         END DO
-
-         hu(:,:) = 0._wp                          ! Ocean depth at U-points
-         hv(:,:) = 0._wp                          ! Ocean depth at V-points
-         ht(:,:) = 0._wp                          ! Ocean depth at T-points
-         DO jk = 1, jpkm1
-            hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-            hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-            ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         !                                        ! Inverse of the local depth
-         hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
-         hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
-
-      END IF
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsum = glob_sum_full(ztmp3d) 
-       IF (lwp) PRINT *, 'total volume correction 02 = ',zsum
-
-!=============================================================================
-
-! 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 ii = 1,jpi
-            DO ij = 1,jpj
-               un(ii,ij,jk) = ub(ii,ij,jk)*pe3u_b(ii,ij,jk)*pumask_b(ii,ij,jk)/fse3u_n(ii,ij,jk)*umask(ii,ij,jk);
-               vn(ii,ij,jk) = vb(ii,ij,jk)*pe3v_b(ii,ij,jk)*pvmask_b(ii,ij,jk)/fse3v_n(ii,ij,jk)*vmask(ii,ij,jk);
-               !IF ( umask(ii,ij,1) == 1 .AND. pumask(ii,ij,1) == 0 ) un(ii,ij,jk) = ub(ii,ij,jk)
-               !IF ( vmask(ii,ij,1) == 1 .AND. pvmask(ii,ij,1) == 0 ) vn(ii,ij,jk) = vb(ii,ij,jk)
-            END DO
-         END DO
-      END DO
-! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
-      ! compute barotropic velocity now and after 
-      ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:); 
-      zbub(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:); 
-      zbvb(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = un(:,:,:)*fse3u_n(:,:,:); 
-      zbun(:,:)   = SUM(ztrp,DIM=3)
-      ztrp(:,:,:) = vn(:,:,:)*fse3v_n(:,:,:); 
-      zbvn(:,:)   = SUM(ztrp,DIM=3)
-! Already know ????????
-      hu1=0.0_wp ;
-      hv1=0.0_wp ;
-      DO jk  = 1,jpk
-        hu1(:,:) = hu1(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-        hv1(:,:) = hv1(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      zucorr = 0._wp
-      zvcorr = 0._wp
-      DO ii = 1,jpi
-         DO ij = 1,jpj
-            IF (zbun(ii,ij) .NE. zbub(ii,ij) .AND. hu1(ii,ij) .NE. 0._wp ) THEN
-               zucorr(ii,ij) = (zbun(ii,ij) - zbub(ii,ij))/hu1(ii,ij)
-            END IF
-            IF (zbvn(ii,ij) .NE. zbvb(ii,ij) .AND. hv1(ii,ij) .NE. 0._wp ) THEN
-               zvcorr(ii,ij) = (zbvn(ii,ij) - zbvb(ii,ij))/hv1(ii,ij)
-            END IF
-         END DO
-      END DO 
-      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,10
-          DO jk = 1,jpk-1
-              zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
-              DO ii = 2,jpi-1
-                  DO ij = 2,jpj-1
-                      iip1=ii+1; iim1=ii-1;
-                      ijp1=ij+1; ijm1=ij-1;
-                      summsk=(ztmask0(iip1,ij,jk)+ztmask0(iim1,ij,jk)+ztmask0(ii,ijp1,jk)+ztmask0(ii,ijm1,jk))
-                      !! horizontal extrapolation
-                      IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
-                         tsn(ii,ij,jk,1)=( zts0(iip1,ij,jk,1)*ztmask0(iip1,ij,jk)     &
-                         &                +zts0(iim1,ij,jk,1)*ztmask0(iim1,ij,jk)     &
-                         &                +zts0(ii,ijp1,jk,1)*ztmask0(ii,ijp1,jk)     &
-                         &                +zts0(ii,ijm1,jk,1)*ztmask0(ii,ijm1,jk))/summsk
-                         tsn(ii,ij,jk,2)=( zts0(iip1,ij,jk,2)*ztmask0(iip1,ij,jk)     &
-                         &                +zts0(iim1,ij,jk,2)*ztmask0(iim1,ij,jk)     &
-                         &                +zts0(ii,ijp1,jk,2)*ztmask0(ii,ijp1,jk)     &
-                         &                +zts0(ii,ijm1,jk,2)*ztmask0(ii,ijm1,jk))/summsk
-                         ztmask1(ii,ij,jk)=1
-                      END IF
-                      !! vertical extrapolation if horizontal extra failed
-                      IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
-                         jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
-                         summsk=(ztmask0(ii,ij,jkm1)+ztmask0(ii,ij,jkp1))
-                         IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0 ) THEN
-                            tsn(ii,ij,jk,1)=( zts0(ii,ij,jkp1,1)*ztmask0(ii,ij,jkp1)     &
-                            &                +zts0(ii,ij,jkm1,1)*ztmask0(ii,ij,jkm1))/summsk
-                            tsn(ii,ij,jk,2)=( zts0(ii,ij,jkp1,2)*ztmask0(ii,ij,jkp1)     &
-                            &                +zts0(ii,ij,jkm1,2)*ztmask0(ii,ij,jkm1))/summsk
-                            ztmask1(ii,ij,jk)=1
-                         END IF
-                      END IF
-                      !! horizontal corner extrapolation if horizontal and vertical failed
-                      IF (zdmask(ii,ij)==1 .AND. summsk==0) THEN
-                         iip1=ii+1; iim1=ii-1;
-                         ijp1=ij+1; ijm1=ij-1;
-                         summsk=(ztmask0(iip1,ijp1,jk)+ztmask0(iim1,ijm1,jk)+ztmask0(iim1,ijp1,jk)+ztmask0(iip1,ijm1,jk))
-                         IF (zdmask(ii,ij)==1 .AND. summsk .NE. 0) THEN
-                            tsn(ii,ij,jk,1)=( zts0(iip1,ijp1,jk,1)*ztmask0(iip1,ijp1,jk)     &
-                            &                +zts0(iim1,ijm1,jk,1)*ztmask0(iim1,ijm1,jk)     &
-                            &                +zts0(iim1,ijp1,jk,1)*ztmask0(iim1,ijp1,jk)     &
-                            &                +zts0(iip1,ijm1,jk,1)*ztmask0(iip1,ijm1,jk))/summsk
-                            tsn(ii,ij,jk,2)=( zts0(iip1,ijp1,jk,2)*ztmask0(iip1,ijp1,jk)     &
-                            &                +zts0(iim1,ijm1,jk,2)*ztmask0(iim1,ijm1,jk)     &
-                            &                +zts0(iim1,ijp1,jk,2)*ztmask0(iim1,ijp1,jk)     &
-                            &                +zts0(iip1,ijm1,jk,2)*ztmask0(iip1,ijm1,jk))/summsk
-                            ztmask1(ii,ij,jk)=1
-                         END IF
-                      END IF
-                  END DO
-              END DO
-          END DO
-          CALL lbc_lnk(tsn(:,:,:,1),'T',1.)
-          CALL lbc_lnk(tsn(:,:,:,2),'T',1.)
-          CALL lbc_lnk(ztmask1,'T',1.)
-          zts0(:,:,:,:) = tsn(:,:,:,:)
-          ztmask0 = ztmask1
-      END DO
-      tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
-      tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
-
-      ! CHECK heat
-      zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
-      zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * zwmaskb(:,:,:) * zwmaskn(:,:,:))
-      IF (lwp) PRINT *, 'CHECK tsn = ',zsumn, zsumb
-
-      ! compute new T/S (interpolation) because of vvl
-      IF ( lk_vvl .AND. .FALSE. ) THEN
-      !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.0_wp ) THEN
-                     !compute weight
-                     zdzp1 = MAX(0._wp,fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
-                     zdz   =           fsdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk  ) 
-                     zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk  ) - fsdepw_n(ji,jj,jk  ))
-                     IF (zdz .LT. 0.0_wp) THEN 
-                        PRINT *, 'ERROR dz n ', ji,jj,jk,zdz,fsdepw_n(ji,jj,jk+1),fsdepw_n(ji,jj,jk),fsdepw_n(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ',fse3t_n (ji,jj,jk+1),fse3t_n (ji,jj,jk),fse3t_n (ji,jj,jk-1), sshn(ji,jj)
-                        PRINT *, 'ERROR dz b             = ',pdepw_b(ji,jj,jk+1),pdepw_b(ji,jj,jk),pdepw_b(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ',fse3t_b (ji,jj,jk+1),fse3t_b (ji,jj,jk),fse3t_b (ji,jj,jk-1), sshb(ji,jj)
-                        PRINT *, 'ERROR dz 0             = ',  e3t_0 (ji,jj,jk+1),  e3t_0 (ji,jj,jk),  e3t_0 (ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ',  tmask (ji,jj,jk+1),  tmask (ji,jj,jk),  tmask (ji,jj,jk-1)
-                        PRINT *, 'ERROR dz n             = ', zwmaskn(ji,jj,jk+1), zwmaskn(ji,jj,jk), zwmaskn(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', ptmask_b(ji,jj,jk+1), ptmask_b(ji,jj,jk), ptmask_b(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', zwmaskb(ji,jj,jk+1), zwmaskb(ji,jj,jk), zwmaskb(ji,jj,jk-1)
-                        PRINT *, 'ERROR dz b             = ', gdepw_0(ji,jj,jk+1), gdepw_0(ji,jj,jk), gdepw_0(ji,jj,jk-1)
-                        STOP
-                     END IF
-                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem)      )/fse3t_n(ji,jj,jk)
-                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal)      &
-                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal)      &
-                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal)      )/fse3t_n(ji,jj,jk)
-                  END IF
-               END DO
-            END DO
-         END DO               
-      END IF
-
-      ! CHECK heat
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = fse3t_n(2:jpi-1,2:jpj-1,:) * tmask(2:jpi-1,2:jpj-1,:) &
-       &                         - pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsum = glob_sum_full(ztmp3d) 
-       IF (lwp) PRINT *, 'total volume correction 03 = ',zsum
-
-      WHERE (tmask(:,:,:) == 1.0 .AND. tsn(:,:,:,2) == 0._wp) 
-         tsn(:,:,:,2)=  -99._wp
-         tmask(:,:,:) = 0.0
-         umask(:,:,:) = 0.0
-         vmask(:,:,:) = 0.0
-      END WHERE
-
-      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
-
-!=============================================================================
-! put this stuff in a subroutine ????
-      !IF ( ln_hfb ) CALL scb_iscpl_cons(pvol_flx, pts_flx, pe3t_b, zssh0)
-    ! initialisation
-    zde3t   (:,:)     = 0.0_wp
-    pvol_flx(:,:,:  ) = 0.0_wp
-    pts_flx (:,:,:,:) = 0.0_wp
-    IF ( ln_hfb ) THEN
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 0 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 0 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 0 = ',zsum
-
-       ! mask tsn and tsb (should be useless)
-       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(:,:,:);
-       ! why tsn = tsb ????? (bug ???????)
-       ! diagnose non conservation of heat, salt and volume 
-       r1_tiscpl = 1._wp / (rn_fiscpl * rn_rdt) 
-       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-       IF ( lk_vvl ) zssh0 = 0.0_wp
-       DO jk = 1,jpk-1
-          DO ii = 2,jpi-1
-             DO ij = 2,jpj-1
-                ! volume differences
-                zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
-                
-                ! shh changes
-                IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN 
-                   zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
-                   zssh0(ii,ij) = 0._wp
-                END IF
-
-                ! ocean cell now
-!                IF (tmask(ii,ij,jk) == 1._wp) THEN
-                ! case where we open, enlarge or thin a cell :
-                   pvol_flx(ii,ij,jk)       =                          zde3t(ii,ij) * r1_tiscpl
-                   pts_flx (ii,ij,jk,jp_sal)=   tsn(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl 
-                   pts_flx (ii,ij,jk,jp_tem)=   tsn(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
-!                END IF
-             END DO
-          END DO
-       END DO
-       ! glob_sum_full because with glob summ some data can be masked. WARNING the halo have to be set at 0
-       PRINT *, 'test ', narea, SUM(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt, SUM(pvol_flx(2:jpi-1,2:jpj-1,:)) * rn_fiscpl * rn_rdt
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 1 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 1 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 1 = ',zsum
-
-       zssh0(:,:)        = sshn(:,:) * ssmask(:,:) - sshb(:,:) * psmask_b(:,:)
-       IF ( lk_vvl ) zssh0 = 0.0_wp
-       DO jk = 1,jpk-1
-          DO ii = 2,jpi-1
-             DO ij = 2,jpj-1
-                ! volume differences
-                zde3t(ii,ij) = fse3t_n(ii,ij,jk) * tmask(ii,ij,jk) - pe3t_b(ii,ij,jk) * ptmask_b(ii,ij,jk);
-                
-                ! shh changes
-                !IF ( ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) .AND. zssh0(ii,ij) .NE. 0._wp) THEN 
-                IF ( ptmask_b(ii,ij,jk) == 1 .OR. tmask(ii,ij,jk) == 1 ) THEN 
-                   zde3t(ii,ij) = zde3t(ii,ij) + zssh0(ii,ij)
-                   zssh0(ii,ij) = 0._wp
-                END IF
-
-                ! ocean cell before and mask cell now
-                IF ( tmask(ii,ij,jk) == 0._wp .AND. ptmask_b(ii,ij,jk) == 1._wp ) THEN
-                   ! case where we close a cell and adjacent cell open
-                   pvol_flx(ii,ij,jk)       =                         zde3t(ii,ij) * r1_tiscpl
-                   pts_flx (ii,ij,jk,jp_sal)=  tsb(ii,ij,jk,jp_sal) * zde3t(ii,ij) * r1_tiscpl 
-                   pts_flx (ii,ij,jk,jp_tem)=  tsb(ii,ij,jk,jp_tem) * zde3t(ii,ij) * r1_tiscpl
-
-                   iip1=ii+1 ; iim1=ii-1 ; ijp1=ij+1 ; ijm1=ij-1 ;
-
-                   zsum =   e12t(ii  ,ijp1) * tmask(ii  ,ijp1,jk) + e12t(ii  ,ijm1) * tmask(ii  ,ijm1,jk) &
-                     &    + e12t(iim1,ij  ) * tmask(iim1,ij  ,jk) + e12t(iip1,ij  ) * tmask(iip1,ij  ,jk) 
-
-                   IF ( zsum .NE. 0._wp ) THEN
-                      ziip1_ratio = e12t(iip1,ij  ) * tmask(iip1,ij  ,jk) / zsum
-                      ziim1_ratio = e12t(iim1,ij  ) * tmask(iim1,ij  ,jk) / zsum
-                      zijp1_ratio = e12t(ii  ,ijp1) * tmask(ii  ,ijp1,jk) / zsum
-                      zijm1_ratio = e12t(ii  ,ijm1) * tmask(ii  ,ijm1,jk) / zsum
-
-                      pvol_flx(ii  ,ijp1,jk       ) = pvol_flx(ii  ,ijp1,jk       ) + pvol_flx(ii,ij,jk       ) * zijp1_ratio
-                      pvol_flx(ii  ,ijm1,jk       ) = pvol_flx(ii  ,ijm1,jk       ) + pvol_flx(ii,ij,jk       ) * zijm1_ratio
-                      pvol_flx(iip1,ij  ,jk       ) = pvol_flx(iip1,ij  ,jk       ) + pvol_flx(ii,ij,jk       ) * ziip1_ratio
-                      pvol_flx(iim1,ij  ,jk       ) = pvol_flx(iim1,ij  ,jk       ) + pvol_flx(ii,ij,jk       ) * ziim1_ratio
-                      pts_flx (ii  ,ijp1,jk,jp_sal) = pts_flx (ii  ,ijp1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijp1_ratio
-                      pts_flx (ii  ,ijm1,jk,jp_sal) = pts_flx (ii  ,ijm1,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * zijm1_ratio
-                      pts_flx (iip1,ij  ,jk,jp_sal) = pts_flx (iip1,ij  ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziip1_ratio
-                      pts_flx (iim1,ij  ,jk,jp_sal) = pts_flx (iim1,ij  ,jk,jp_sal) + pts_flx (ii,ij,jk,jp_sal) * ziim1_ratio
-                      pts_flx (ii  ,ijp1,jk,jp_tem) = pts_flx (ii  ,ijp1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijp1_ratio
-                      pts_flx (ii  ,ijm1,jk,jp_tem) = pts_flx (ii  ,ijm1,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * zijm1_ratio
-                      pts_flx (iip1,ij  ,jk,jp_tem) = pts_flx (iip1,ij  ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziip1_ratio
-                      pts_flx (iim1,ij  ,jk,jp_tem) = pts_flx (iim1,ij  ,jk,jp_tem) + pts_flx (ii,ij,jk,jp_tem) * ziim1_ratio
-
-                      ! set to 0 the cell we distributed over neigbourg cells
-                      pvol_flx(ii,ij,jk       ) = 0._wp
-                      pts_flx (ii,ij,jk,jp_sal) = 0._wp
-                      pts_flx (ii,ij,jk,jp_tem) = 0._wp
-                   ELSE IF (zsum .EQ. 0._wp ) THEN
-                      ! case where we close a cell and no adjacent cell open
-                      IF ( tmask(ii,ij,jk+1) .EQ. 1._wp ) THEN
-                         pvol_flx(ii,ij,jk+1)       =  pvol_flx(ii,ij,jk+1)        + pvol_flx(ii,ij,jk)
-                         pts_flx (ii,ij,jk+1,jp_sal)=  pts_flx (ii,ij,jk+1,jp_sal) + pts_flx (ii,ij,jk,jp_sal)
-                         pts_flx (ii,ij,jk+1,jp_tem)=  pts_flx (ii,ij,jk+1,jp_tem) + pts_flx (ii,ij,jk,jp_tem)
-
-                         ! set to 0 the cell we distributed over neigbourg cells
-                         pvol_flx(ii,ij,jk       ) = 0._wp
-                         pts_flx (ii,ij,jk,jp_sal) = 0._wp
-                         pts_flx (ii,ij,jk,jp_tem) = 0._wp
-                      ELSE
-                      ! case no adjacent cell on the horizontal and on the vertical
-                         PRINT *,        'W A R N I N G iscpl: no adjacent cell on the vertical and horizontal'
-                         PRINT *,        '                     ',mig(ii),' ',mjg(ij),' ',jk
-                         PRINT *,        '                     ',ii,' ',ij,' ',jk,' ',narea
-                         PRINT *,        '                     conservation broken '
-! We deal with this point later.
-                      END IF
-                   END IF
-                END IF
-             END DO
-          END DO
-       END DO
-
-       zsum = glob_sum_full(pvol_flx(:,:,:)      ) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total volume correction 2 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total heat correction 2 = ',zsum
-       zsum = glob_sum_full(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       IF (lwp) PRINT *, 'total salt correction 2 = ',zsum
-
-          ! allocation and initialisation of the list of problematic point
-          ALLOCATE(vnpts(jpnij))
-          vnpts(:)=0
-
-          ! fill narea location with the number of problematic point
-          DO jk = 1,jpk-1
-             DO ii = 2,jpi-1
-                DO ij = 2,jpj-1
-                   IF (          ptmask_b(ii,ij,jk)      == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
-                   &   .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND.     tmask(ii,ij,jk+1)       == 0 ) THEN
-                      vnpts(narea) = vnpts(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(vnpts,jpnij) 
-
-          ! size of the new variable
-          npts  = SUM(vnpts)    
-         
-          ! allocation of the coordiantes, 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 ; zlat(:) = -1.0e20
-          zcorr_vol(:) = 0.0_wp
-          zcorr_sal(:) = 0.0_wp
-          zcorr_tem(:) = 0.0_wp
-
-          ! fill new variable
-          jpts = SUM(vnpts(1:narea-1))
-          DO jk = 1,jpk-1
-             DO ii = 2,jpi-1
-                DO ij = 2,jpj-1
-                   IF (          ptmask_b(ii,ij,jk)      == 1 .AND. SUM(tmask(ii-1:ii+1,ij,jk)) == 0 &
-                   &   .AND. SUM(tmask(ii,ij-1:ij+1,jk)) == 0 .AND.     tmask(ii,ij,jk+1)       == 0 ) THEN
-                      jpts = jpts + 1
-                      PRINT *, 'corrected point ', narea, ii, ij, jk, jpts
-                      ixpts(jpts) = ii           ; iypts(jpts) = ij ; izpts(jpts) = jk
-                      zlon (jpts) = glamt(ii,ij) ; zlat (jpts) = gphit(ii,ij)
-                      zcorr_vol(jpts) = pvol_flx(ii,ij,jk)
-                      zcorr_sal(jpts) = pts_flx (ii,ij,jk,jp_sal)
-                      zcorr_tem(jpts) = pts_flx (ii,ij,jk,jp_tem)
-                      ! set flx to 0 (safer)
-                      pvol_flx(ii,ij,jk       ) = 0.0_wp
-                      pts_flx (ii,ij,jk,jp_sal) = 0.0_wp
-                      pts_flx (ii,ij,jk,jp_tem) = 0.0_wp
-                      PRINT *, zcorr_vol(jpts)*rn_fiscpl*rn_rdt, zcorr_sal(jpts)*rn_fiscpl*rn_rdt, zcorr_tem(jpts)*rn_fiscpl*rn_rdt
-                   END IF
-                END DO
-             END DO
-          END DO
-
-          ! build array of total problematic point on each cpu (share to each cpu)
-          CALL mpp_max(zlat ,npts)
-          CALL mpp_max(zlon ,npts)
-          CALL mpp_max(izpts,npts)
-
-          ! put correction term in the closest cell          
-          PRINT *, 'corrected point1 ', narea, zlon, zlat, izpts
-          DO jpts = 1,npts
-             CALL dom_ngb(zlon(jpts), zlat(jpts), ixpts(jpts), iypts(jpts),'T', izpts(jpts))
-             PRINT *, 'corrected point2 ', narea, jpts, ixpts(jpts), iypts(jpts), izpts(jpts)
-             DO ij = mj0(iypts(jpts)),mj1(iypts(jpts))
-                DO ii = mi0(ixpts(jpts)),mi1(ixpts(jpts))
-                   jk = izpts(jpts)
-                   pvol_flx(ii,ij,jk)        =  pvol_flx(ii,ij,jk       ) + zcorr_vol(jpts)
-                   pts_flx (ii,ij,jk,jp_sal) =  pts_flx (ii,ij,jk,jp_sal) + zcorr_sal(jpts)
-                   pts_flx (ii,ij,jk,jp_tem) =  pts_flx (ii,ij,jk,jp_tem) + zcorr_tem(jpts)
-                END DO
-             END DO
-          END DO
-          ! deallocate variable 
-          DEALLOCATE(vnpts)
-          DEALLOCATE(ixpts, iypts, izpts, zcorr_vol, zcorr_sal, zcorr_tem, zlon, zlat)
-     
-       ! add contribution store on the allo (lbclnk remove one of the contribution)
-       pvol_flx(:,:,:       ) = pvol_flx(:,:,:       ) * tmask(:,:,:)
-       pts_flx (:,:,:,jp_tem) = pts_flx (:,:,:,jp_tem) * tmask(:,:,:)
-       pts_flx (:,:,:,jp_sal) = pts_flx (:,:,:,jp_sal) * tmask(:,:,:)
-
-       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.)
-
-       ! CHECK vol !!!!!!!!! warning tmask_i wrong if deals with before value, so glob_sum wrong for before value!!!!
-       zsumn = glob_sum     ( fse3t_n(:,:,:) * tmask  (:,:,:)) - glob_sum(pvol_flx(:,:,:)) * rn_fiscpl * rn_rdt
-       ztmp3d(:,:,:) = 0.0
-       ztmp3d(2:jpi-1,2:jpj-1,:) = pe3t_b(2:jpi-1,2:jpj-1,:) * ptmask_b(2:jpi-1,2:jpj-1,:)
-       zsumb = glob_sum_full(ztmp3d) 
-       zsum  = glob_sum     ( pvol_flx(:,:,:) * rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK vol = ',zsumn, zsumb, zsumn - zsumb, zsum
-       ! CHECK salt
-       zsumn = glob_sum( tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_sal)) * rn_fiscpl * rn_rdt
-       zsumb = glob_sum( tsb(:,:,:,jp_sal) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
-       zsum  = glob_sum( pts_flx(:,:,:,jp_sal)*rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK salt = ',zsumn, zsumb, zsumn - zsumb, zsum
-       ! CHECK heat
-       zsumn = glob_sum( tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) *  tmask  (:,:,:)) - glob_sum(pts_flx(:,:,:,jp_tem)) * rn_fiscpl * rn_rdt
-       zsumb = glob_sum( tsb(:,:,:,jp_tem) *  pe3t_b(:,:,:) * ptmask_b(:,:,:))
-       zsum  = glob_sum( pts_flx(:,:,:,jp_tem)*rn_fiscpl * rn_rdt)
-       IF (lwp) PRINT *, 'CHECK heat = ',zsumn, zsumb, zsumn - zsumb, zsum
-          
-    END IF
-!=============================================================================
-    !! now T/S/U are compute in agreement with the new grid.
-    !! need to interpolate on vvl grid
-    !! 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(:,:,:)
-    !! set depth_b = depth_n
-       fsdept_b(:,:,:) = fsdept_n(:,:,:)
-       fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
-    !! set hu_b
-       hu_b(:,:) = hu(:,:)
-       hv_b(:,:) = hv(:,:)
-       hur_b(:,:) = hur(:,:)
-       hvr_b(:,:) = hvr(:,:)
-    END IF
-    !! 
-    CALL wrk_dealloc(jpi,jpj,jpk,2, zts0                                   )
-    CALL wrk_dealloc(jpi,jpj,jpk,   ztmask0, ztmask1 , ztrp                ) 
-    CALL wrk_dealloc(jpi,jpj,jpk,   zwmaskn, zwmaskb , ztmp3d              ) 
-    CALL wrk_dealloc(jpi,jpj,       zsmask0, zsmask1                       ) 
-    CALL wrk_dealloc(jpi,jpj,       zdmask , zdsmask, zvcorr, zucorr, zde3t) 
-    CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
-    CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , hu1 , hv1            ) 
-   END SUBROUTINE rst_iscpl_interpol
-
-   SUBROUTINE sbc_iscpl_div( phdivn )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_iscpl_div  ***
-      !!
-      !! ** Purpose :   update the horizontal divergence with the iscpl imbalance
-      !!
-      !! ** Method  :
-      !!                CAUTION : iscpl is positive (inflow) increase the
-      !!                          divergence and expressed in m/s
-      !!
-      !! ** Action  :   phdivn   decreased by the runoff inflow
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + hdiv_iscpl(ji,jj,jk) / fse3t_n(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-   END SUBROUTINE sbc_iscpl_div
-
-END MODULE sbc_iscpl
