Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90	(revision 12746)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/domain.F90	(revision 12747)
@@ -284,5 +284,5 @@
          &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
          &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
-         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios, ln_rstdate, ln_rst_eos
+         &             ln_cfmeta, ln_xios_read, nn_wxios, ln_rstdate, ln_rst_eos
 
       NAMELIST/namdom/ ln_linssh, rn_rdt, rn_atfp, ln_crs, ln_meshmask
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/dommsk.F90	(revision 12746)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DOM/dommsk.F90	(revision 12747)
@@ -93,5 +93,5 @@
       INTEGER  ::   ijf, ijl       !   -       -
       INTEGER  ::   iktop, ikbot   !   -       -
-      INTEGER  ::   ios, inum
+      INTEGER  ::   ios
       !!
       INTEGER  :: inum             !  logical unit for shlat2d
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 12746)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/DYN/dynnxt.F90	(revision 12747)
@@ -29,5 +29,4 @@
    USE sbc_oce        ! Surface boundary condition: ocean fields
    USE sbcrnf         ! river runoffs
-   USE sbcisf         ! ice shelf
    USE phycst         ! physical constants
    USE dynadv         ! dynamics: vector invariant versus flux form
@@ -35,4 +34,5 @@
    USE domvvl         ! variable volume
    USE bdy_oce   , ONLY: ln_bdy
+   USE isf_oce        ! ice shelf
    USE bdydta         ! ocean open boundary conditions
    USE bdydyn         ! ocean open boundary conditions
@@ -98,10 +98,10 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zfwf
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
       !
       IF( ln_timing    )   CALL timing_start('dyn_nxt')
-      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
+      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)    , zfwf(jpi,jpj))
       IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
       !
@@ -237,10 +237,26 @@
             ENDIF
 
-            IF ( ln_isf ) THEN
-               DO jk = 1, jpkm1
-                  e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_b(:,:) - fwfisf(:,:) ) * tmask(:,:,jk) & 
-                                &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
-               END DO
-            ENDIF
+            IF ( ln_isf ) THEN 
+               IF ( ln_isfpar_mlt ) THEN
+                  DO jk = 1, jpkm1
+                     e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_cav_b(:,:) - fwfisf_cav(:,:) ) * tmask(:,:,jk) & 
+                                   &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+                  END DO
+               END IF
+               IF ( ln_isfpar_mlt ) THEN
+                  DO jk = 1, jpkm1
+                     e3t_b(:,:,jk) = e3t_b(:,:,jk) - zcoef * ( fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * tmask(:,:,jk) & 
+                                   &                       * e3t_n(:,:,jk) /  ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+                  END DO
+               END IF
+               IF ( ln_isfcpl .AND. ln_rstart .AND. kt == nit000+1 ) THEN
+                  ! risfcpl_vol in m3/s
+                  zfwf(:,:) = SUM(risfcpl_vol,DIM=3)
+                  DO jk = 1, jpkm1
+                     e3t_b(:,:,jk) =   e3t_b(:,:,jk) - zcoef * zfwf(:,:) * r1_e1e2t(:,:) * r1_rau0 * tmask(:,:,jk) &
+                                   &                         * e3t_n(:,:,jk) / ( ht_n(:,:) + 1._wp - ssmask(:,:) )
+                  END DO
+               END IF
+            END IF
             !
             IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
@@ -357,5 +373,5 @@
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
+      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve, zfwf )
       IF( l_trddyn     )   DEALLOCATE( zua, zva )
       IF( ln_timing    )   CALL timing_stop('dyn_nxt')
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isf_oce.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isf_oce.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isf_oce.F90	(revision 12747)
@@ -0,0 +1,266 @@
+MODULE isf_oce
+   !!======================================================================
+   !!                       ***  MODULE  sbcisf  ***
+   !! Surface module :  compute iceshelf melt and heat flux
+   !!======================================================================
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!            4.1  !  2019-09  (P. Mathiot) Split param/explicit ice shelf and re-organisation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isf          : define and allocate ice shelf variables
+   !!----------------------------------------------------------------------
+
+   USE par_oce       , ONLY: jpi, jpj, jpk
+   USE in_out_manager, ONLY: wp, jpts ! I/O manager
+   USE lib_mpp       , ONLY: ctl_stop, mpp_sum      ! MPP library
+   USE fldread        ! read input fields
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC   isf_alloc, isf_alloc_par, isf_alloc_cav, isf_alloc_cpl, isf_dealloc_cpl
+   !
+   !-------------------------------------------------------
+   ! 0 :              namelist parameter
+   !-------------------------------------------------------
+   !
+   ! 0.1 -------- ice shelf cavity parameter --------------
+   CHARACTER(LEN=256), PUBLIC :: cn_isfdir
+   LOGICAL           , PUBLIC :: ln_isf
+   LOGICAL           , PUBLIC :: ln_isfdebug
+   !
+   ! 0.2 -------- ice shelf cavity opened namelist parameter -------------
+   LOGICAL           , PUBLIC :: ln_isfcav_mlt   !: logical for the use of ice shelf parametrisation
+   REAL(wp)          , PUBLIC :: rn_gammat0      !: temperature exchange coeficient    []
+   REAL(wp)          , PUBLIC :: rn_gammas0      !: salinity    exchange coeficient    []
+   REAL(wp)          , PUBLIC :: rn_htbl         !: Losch top boundary layer thickness [m]
+   REAL(wp)          , PUBLIC :: rn_isfload_T    !: 
+   REAL(wp)          , PUBLIC :: rn_isfload_S    !: 
+   CHARACTER(LEN=256), PUBLIC :: cn_gammablk     !: gamma formulation
+   CHARACTER(LEN=256), PUBLIC :: cn_isfcav_mlt   !: melt formulation (cavity/param)
+   CHARACTER(LEN=256), PUBLIC :: cn_isfload      !: ice shelf load computation method
+   TYPE(FLD_N)       , PUBLIC :: sn_isfcav_fwf   !: information about the isf melting file to be read
+   !
+   ! 0.3 -------- ice shelf cavity parametrised namelist parameter -------------
+   LOGICAL           , PUBLIC :: ln_isfpar_mlt   !: logical for the computation of melt inside the cavity
+   CHARACTER(LEN=256), PUBLIC :: cn_isfpar_mlt   !: melt formulation (cavity/param)
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_fwf   !: information about the isf melting file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmax  !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_zmin  !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC :: sn_isfpar_Leff  !: information about the effective length     file to be read
+   !
+   ! 0.4 -------- coupling namelist parameter -------------
+   LOGICAL, PUBLIC :: ln_isfcpl      !:
+   LOGICAL, PUBLIC :: ln_isfcpl_cons !:
+   INTEGER, PUBLIC :: nn_drown       !:
+   !
+   !-------------------------------------------------------
+   ! 1 :              ice shelf parameter
+   !-------------------------------------------------------
+   !
+   REAL(wp), PARAMETER, PUBLIC :: rLfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
+   REAL(wp), PARAMETER, PUBLIC :: rcpisf = 2000.0_wp       !: specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PARAMETER, PUBLIC :: rkappa = 1.54e-6_wp      !: heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PARAMETER, PUBLIC :: rhoisf = 920.0_wp        !: volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PARAMETER, PUBLIC :: rtsurf = -20.0           !: surface temperature                    [C]
+   !
+   !-------------------------------------------------------
+   ! 2 :              ice shelf global variables
+   !-------------------------------------------------------
+   !
+   ! 2.1 -------- ice shelf cavity parameter --------------
+   LOGICAL , PUBLIC            :: l_isfoasis
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   risfload                    !: ice shelf load
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    ::   fwfisf_oasis
+   !
+   ! 2.2 -------- ice shelf cavity melt namelist parameter -------------
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: mskisf_cav                    !:
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: misfkt_cav   , misfkb_cav     !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf_tbl_cav, rfrac_tbl_cav  !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: fwfisf_cav   , fwfisf_cav_b   !: before and now net fwf from the ice shelf        [kg/m2/s]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     :: sf_isfcav_fwf                 !:
+   !
+   REAL(wp) , PUBLIC                                      :: risf_lamb1, risf_lamb2, risf_lamb3  ! freezing point linearization coeficient
+   !
+   ! 2.3 -------- ice shelf param. melt namelist parameter -------------
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: mskisf_par                    !:
+   INTEGER  , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: misfkt_par   , misfkb_par     !:
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf_tbl_par, rfrac_tbl_par  !: 
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: fwfisf_par   , fwfisf_par_b   !: before and now net fwf from the ice shelf        [kg/m2/s]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_par_tsc , risf_par_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     :: sf_isfpar_fwf                 !:
+   !
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: rhisf0_tbl_par                !: thickness of tbl (initial value)  [m]
+   REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: risfLeff                      !:
+   !
+   ! 2.4 -------- coupling namelist parameter -------------
+   INTEGER , PUBLIC                                        ::   nstp_iscpl   !:
+   REAL(wp), PUBLIC                                        ::   rdt_iscpl    !: 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   risfcpl_ssh, risfcpl_cons_ssh, risfcpl_cons_ssh_b               !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   risfcpl_vol, risfcpl_cons_vol, risfcpl_cons_vol_b  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   risfcpl_tsc, risfcpl_cons_tsc, risfcpl_cons_tsc_b  !:
+   !
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE isf_alloc_par()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc_par  ***
+      !!
+      !! ** Purpose : 
+      !!
+      !! ** Method  : 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0       ! set to zero if no array to be allocated
+      !
+      ALLOCATE(risfLeff(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(misfkt_par(jpi,jpj), misfkb_par(jpi,jpj), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rfrac_tbl_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rhisf_tbl_par(jpi,jpj), rhisf0_tbl_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_par(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
+   END SUBROUTINE isf_alloc_par
+
+   SUBROUTINE isf_alloc_cav()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc_cav  ***
+      !!
+      !! ** Purpose : 
+      !!
+      !! ** Method  : 
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0       ! set to zero if no array to be allocated
+      !
+      ALLOCATE(misfkt_cav(jpi,jpj), misfkb_cav(jpi,jpj), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rfrac_tbl_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( rhisf_tbl_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
+   END SUBROUTINE isf_alloc_cav
+
+   SUBROUTINE isf_alloc_cpl()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc_cpl  ***
+      !!
+      !! ** Purpose : allocate array use for the ice sheet coupling
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0
+      !
+      ALLOCATE( risfcpl_ssh(jpi,jpj), risfcpl_tsc(jpi,jpj,jpk,jpts), risfcpl_vol(jpi,jpj,jpk), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      risfcpl_tsc(:,:,:,:) = 0.0 ; risfcpl_vol(:,:,:) = 0.0 ; risfcpl_ssh(:,:) = 0.0
+
+      IF ( ln_isfcpl_cons) THEN
+         ALLOCATE( risfcpl_cons_tsc(jpi,jpj,jpk,jpts) , risfcpl_cons_vol(jpi,jpj,jpk) ,risfcpl_cons_ssh(jpi,jpj), STAT=ialloc )
+         ierr = ierr + ialloc
+         !
+         risfcpl_cons_tsc(:,:,:,:) = 0.0 ; risfcpl_cons_vol(:,:,:) = 0.0 ; risfcpl_cons_ssh(:,:) = 0.0
+         !
+      END IF
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop('STOP','isfcpl: failed to allocate arrays.')
+      !
+   END SUBROUTINE isf_alloc_cpl
+
+   SUBROUTINE isf_dealloc_cpl()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_dealloc_cpl  ***
+      !!
+      !! ** Purpose : de-allocate useless public 3d array used for ice sheet coupling
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      ierr = 0
+      !
+      DEALLOCATE( risfcpl_ssh, risfcpl_tsc, risfcpl_vol, STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop('STOP','isfcpl: failed to deallocate arrays.')
+      !
+   END SUBROUTINE isf_dealloc_cpl
+
+   SUBROUTINE isf_alloc()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_alloc  ***
+      !!
+      !! ** Purpose : allocate array used for the ice shelf cavity (cav and par)
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr, ialloc
+      !!----------------------------------------------------------------------
+      !
+      ierr = 0       ! set to zero if no array to be allocated
+      !
+      ALLOCATE(fwfisf_par(jpi,jpj)  , fwfisf_par_b(jpi,jpj), &
+         &     fwfisf_cav(jpi,jpj)  , fwfisf_cav_b(jpi,jpj), &
+         &     fwfisf_oasis(jpi,jpj),            STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risf_par_tsc(jpi,jpj,jpts), risf_par_tsc_b(jpi,jpj,jpts), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risf_cav_tsc(jpi,jpj,jpts), risf_cav_tsc_b(jpi,jpj,jpts), STAT=ialloc )
+      ierr = ierr + ialloc
+      !
+      ALLOCATE(risfload(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      ALLOCATE( mskisf_cav(jpi,jpj), STAT=ialloc)
+      ierr = ierr + ialloc
+      !
+      CALL mpp_sum ( 'isf', ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'isf: failed to allocate arrays.' )
+      !
+      ! initalisation of fwf and tsc array to 0
+      risfload(:,:)       = 0.0_wp
+      fwfisf_oasis(:,:)   = 0.0_wp
+      fwfisf_par(:,:)     = 0.0_wp    ; fwfisf_par_b(:,:)     = 0.0_wp
+      fwfisf_cav(:,:)     = 0.0_wp    ; fwfisf_cav_b(:,:)     = 0.0_wp
+      risf_cav_tsc(:,:,:) = 0.0_wp    ; risf_cav_tsc_b(:,:,:) = 0.0_wp
+      risf_par_tsc(:,:,:) = 0.0_wp    ; risf_par_tsc_b(:,:,:) = 0.0_wp
+      !
+
+   END SUBROUTINE isf_alloc
+
+END MODULE isf_oce
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcav.F90	(revision 12747)
@@ -0,0 +1,232 @@
+MODULE isfcav
+   !!======================================================================
+   !!                       ***  MODULE  isfcav  ***
+   !! Ice shelf cavity module :  update ice shelf melting under ice
+   !!                            shelf
+   !!======================================================================
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!            4.1  !  2019-09  (P. Mathiot) Split ice shelf cavity and ice shelf parametrisation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isf_cav       : update ice shelf melting under ice shelf
+   !!----------------------------------------------------------------------
+   USE isf_oce        ! ice shelf public variables
+   !
+   USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
+   USE isfutils , ONLY: debug          ! ice shelf debug subroutine
+   USE isftbl   , ONLY: isf_tbl        ! ice shelf top boundary layer properties subroutine
+   USE isfcavmlt, ONLY: isfcav_mlt     ! ice shelf melt formulation subroutine
+   USE isfcavgam, ONLY: isfcav_gammats ! ice shelf melt exchange coeficient subroutine
+   USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
+   !
+   USE oce      , ONLY: tsn                    ! ocean tracers
+   USE par_oce  , ONLY: jpi,jpj                ! ocean space and time domain
+   USE phycst   , ONLY: grav,rau0,r1_rau0_rcp  ! physical constants
+   USE eosbn2   , ONLY: ln_teos10               ! use ln_teos10 or not
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE fldread        ! read input field at current time step
+   USE lbclnk         ! lbclnk
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC   isf_cav, isf_cav_init ! routine called in isfmlt
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE isf_cav( kt, ptsc, pqfwf )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE isf_cav  ***
+      !!
+      !! ** Purpose :   handle surface boundary condition under ice shelf
+      !!
+      !! ** Method  :   based on Mathiot et al. (2017)
+      !!
+      !! ** Action  :   - compute geometry of the Losch top bournary layer (see Losch et al. 2008)
+      !!                - depending on the chooses option
+      !!                   - compute temperature/salt in the tbl
+      !!                   - compute exchange coeficient
+      !!                   - compute heat and fwf fluxes
+      !!                   - output
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pqfwf  ! ice shelf melt (>0 out)
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc   ! T & S ice shelf cavity contents
+      !!-------------------------- IN  --------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      LOGICAL :: lit
+      INTEGER :: nit
+      REAL(wp) :: zerr
+      REAL(wp), DIMENSION(jpi,jpj) :: zqlat, zqoce, zqhc, zqh  ! heat fluxes
+      REAL(wp), DIMENSION(jpi,jpj) :: zqoce_b                  !
+      REAL(wp), DIMENSION(jpi,jpj) :: zgammat, zgammas         ! exchange coeficient
+      REAL(wp), DIMENSION(jpi,jpj) :: zttbl, zstbl             ! temp. and sal. in top boundary layer
+      !!---------------------------------------------------------------------
+      !
+      ! compute T/S/U/V for the top boundary layer
+      CALL isf_tbl(tsn(:,:,:,jp_tem), zttbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
+      CALL isf_tbl(tsn(:,:,:,jp_sal), zstbl(:,:),'T', misfkt_cav, rhisf_tbl_cav, misfkb_cav, rfrac_tbl_cav )
+      ! 
+      ! output T/S/U/V for the top boundary layer
+      CALL iom_put('ttbl_cav',zttbl(:,:) * mskisf_cav(:,:))
+      CALL iom_put('stbl'    ,zstbl(:,:) * mskisf_cav(:,:))
+      !
+      ! initialisation
+      IF (TRIM(cn_gammablk) == 'vel_stab' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rau0_rcp ! last time step total heat fluxes (to speed up convergence)
+      !
+      ! compute ice shelf melting
+      nit = 1 ; lit = .TRUE.
+      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
+         !
+         ! compute gammat everywhere (2d)
+         ! useless if melt specified
+         IF ( TRIM(cn_isfcav_mlt) .NE. 'spe' ) THEN
+            CALL isfcav_gammats( zttbl, zstbl, zqoce  , pqfwf,  &
+               &                               zgammat, zgammas )
+         END IF
+         !   
+         ! compute tfrz, latent heat and melt (2d)
+         CALL isfcav_mlt(kt, zgammat, zgammas, zttbl, zstbl, &
+            &                         zqhc   , zqoce, pqfwf  )
+         !
+         ! define if we need to iterate
+         SELECT CASE ( cn_gammablk )
+         CASE ( 'spe','vel' )
+            ! no convergence needed
+            lit = .FALSE.
+         CASE ( 'vel_stab' )
+            ! compute error between 2 iterations
+            zerr = MAXVAL(ABS(zqoce(:,:) - zqoce_b(:,:)))
+            !
+            ! define if iteration needed
+            IF (nit >= 100) THEN              ! too much iteration
+               CALL ctl_stop( 'STOP', 'isf_cav: vel_stab gamma formulation had too many iterations ...' )
+            ELSE IF ( zerr <= 0.01_wp ) THEN  ! convergence is achieve
+               lit = .FALSE.
+            ELSE                              ! converge is not yet achieve
+               nit = nit + 1
+               zqoce_b(:,:) = zqoce(:,:)
+            END IF
+         END SELECT
+         !
+      END DO
+      !
+      ! compute heat and water flux ( > 0 out )
+      pqfwf(:,:) = pqfwf(:,:) * mskisf_cav(:,:)
+      zqoce(:,:) = zqoce(:,:) * mskisf_cav(:,:)
+      zqhc (:,:) = zqhc(:,:)  * mskisf_cav(:,:)
+      !
+      ! compute heat content flux ( > 0 out )
+      zqlat(:,:) = - pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
+      !
+      ! total heat flux ( >0 out )
+      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
+      !
+      ! lbclnk on melt
+      CALL lbc_lnk_multi( 'isfmlt', zqh, 'T', 1., pqfwf, 'T', 1.)
+      !
+      ! output fluxes
+      CALL isf_diags_flx( misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
+      !
+      ! set temperature content
+      ptsc(:,:,jp_tem) = - zqh(:,:) * r1_rau0_rcp
+      !
+      ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+      IF (lrst_oce) CALL isfrst_write(kt, 'cav', ptsc, pqfwf)
+      !
+      IF ( ln_isfdebug ) THEN
+         CALL debug('isf_cav: ptsc T',ptsc(:,:,1))
+         CALL debug('isf_cav: ptsc S',ptsc(:,:,2))
+         CALL debug('isf_cav: pqfwf fwf',pqfwf(:,:))
+      END IF
+      !
+   END SUBROUTINE isf_cav
+
+   SUBROUTINE isf_cav_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_cav_init ***
+      !!
+      !! ** Purpose : initialisation of variable needed to compute melt under an ice shelf
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER :: ierr
+      !!---------------------------------------------------------------------
+      !
+      !==============
+      ! 0: allocation
+      !==============
+      !
+      CALL isf_alloc_cav()
+      !
+      !==================
+      ! 1: initialisation
+      !==================
+      !
+      ! top and bottom level of the 'top boundary layer'
+      misfkt_cav(:,:)    = mikt(:,:) ; misfkb_cav(:,:)    = 1
+      !
+      ! thickness of 'tbl' and fraction of bottom cell affected by 'tbl'
+      rhisf_tbl_cav(:,:) = 0.0_wp    ; rfrac_tbl_cav(:,:) = 0.0_wp
+      !
+      ! cavity mask
+      mskisf_cav(:,:)    = (1._wp - tmask(:,:,1)) * ssmask(:,:)
+      !
+      !================
+      ! 2: read restart
+      !================
+      !
+      ! read cav variable from restart
+      IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
+      !
+      !==========================================
+      ! 3: specific allocation and initialisation (depending of scheme choice)
+      !==========================================
+      !
+      SELECT CASE ( TRIM(cn_isfcav_mlt) )
+      CASE( 'spe' )
+
+         ALLOCATE( sf_isfcav_fwf(1), STAT=ierr )
+         ALLOCATE( sf_isfcav_fwf(1)%fnow(jpi,jpj,1), sf_isfcav_fwf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_isfcav_fwf, (/ sn_isfcav_fwf /), cn_isfdir, 'isf_cav_init', 'read fresh water flux isf data', 'namisf' )
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '  ==>> The ice shelf melt inside the cavity is read from forcing files'
+
+      CASE( '2eq' )
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '  ==>> The original ISOMIP melt formulation is used to compute melt under the ice shelves'
+
+      CASE( '3eq' )
+         ! coeficient for linearisation of potential tfreez
+         ! Crude approximation for pressure (but commonly used)
+         IF ( ln_teos10 ) THEN   ! linearisation from Jourdain et al. (2017)
+            risf_lamb1 =-0.0564_wp
+            risf_lamb2 = 0.0773_wp
+            risf_lamb3 =-7.8633e-8 * grav * rau0
+         ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
+            risf_lamb1 =-0.0573_wp
+            risf_lamb2 = 0.0832_wp
+            risf_lamb3 =-7.5300e-8 * grav * rau0
+         ENDIF
+
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '  ==>> The 3 equations melt formulation is used to compute melt under the ice shelves'
+
+      CASE DEFAULT
+         CALL ctl_stop(' cn_isfcav_mlt method unknown (spe, 2eq, 3eq), check namelist')
+      END SELECT
+      !
+   END SUBROUTINE isf_cav_init
+
+END MODULE isfcav
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavgam.F90	(revision 12747)
@@ -0,0 +1,255 @@
+MODULE isfcavgam
+   !!======================================================================
+   !!                       ***  MODULE  isfgammats  ***
+   !! Ice shelf gamma module :  compute exchange coeficient at the ice/ocean interface
+   !!======================================================================
+   !! History :  4.1  !  (P. Mathiot) original
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfcav_gammats       : compute exchange coeficient gamma 
+   !!----------------------------------------------------------------------
+   USE isf_oce
+   USE isfutils, ONLY: debug
+   USE isftbl  , ONLY: isf_tbl
+
+   USE oce     , ONLY: un, vn, rn2         ! ocean dynamics and tracers
+   USE phycst  , ONLY: grav, vkarmn        ! physical constant
+   USE eosbn2  , ONLY: eos_rab             ! equation of state
+   USE zdfdrg  , ONLY: rCd0_top, r_ke0_top ! vertical physics: top/bottom drag coef.
+   USE iom     , ONLY: iom_put             !
+   USE lib_mpp , ONLY: ctl_stop
+
+   USE dom_oce        ! ocean space and time domain
+   USE in_out_manager ! I/O manager
+   !
+   IMPLICIT NONE
+   !
+   PRIVATE
+   !
+   PUBLIC   isfcav_gammats
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+   !
+   !!-----------------------------------------------------------------------------------------------------
+   !!                              PUBLIC SUBROUTINES 
+   !!-----------------------------------------------------------------------------------------------------
+   !
+   SUBROUTINE isfcav_gammats( pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange for heat and fwf flux
+      !!
+      !! ** Method     : select the gamma formulation
+      !!                 3 method available (cst, vel and vel_stab)
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt  , pgs      ! gamma t and gamma s 
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf    ! isf heat and fwf
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl    ! top boundary layer tracer
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)                :: zutbl, zvtbl    ! top boundary layer velocity
+      !!---------------------------------------------------------------------
+      !
+      !==========================================
+      ! 1.: compute velocity in the tbl if needed
+      !==========================================
+      !
+      SELECT CASE ( cn_gammablk )
+      CASE ( 'spe'  ) 
+         ! gamma is constant (specified in namelist)
+         ! nothing to do
+      CASE ('vel', 'vel_stab')
+         ! compute velocity in tbl
+         CALL isf_tbl(un(:,:,:) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
+         CALL isf_tbl(vn(:,:,:) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav)
+         !
+         ! mask velocity in tbl with ice shelf mask
+         zutbl(:,:) = zutbl(:,:) * mskisf_cav(:,:)
+         zvtbl(:,:) = zvtbl(:,:) * mskisf_cav(:,:)
+         !
+         ! output
+         CALL iom_put('utbl',zutbl(:,:))
+         CALL iom_put('vtbl',zvtbl(:,:))
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
+      END SELECT
+      ! 
+      !==========================================
+      ! 2.: compute gamma
+      !==========================================
+      !
+      SELECT CASE ( cn_gammablk )
+      CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+      CASE ( 'vel' ) ! gamma is proportional to u*
+         CALL gammats_vel      (              zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
+      CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
+         CALL gammats_vel_stab (pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
+      END SELECT
+      !
+      !==========================================
+      ! 3.: output and debug
+      !==========================================
+      !
+      CALL iom_put('isfgammat', pgt(:,:))
+      CALL iom_put('isfgammas', pgs(:,:))
+      !
+      IF (ln_isfdebug) THEN
+         CALL debug( 'isfcav_gam pgt:', pgt(:,:) )
+         CALL debug( 'isfcav_gam pgs:', pgs(:,:) )
+      END IF
+      !
+   END SUBROUTINE isfcav_gammats
+   !
+   !!-----------------------------------------------------------------------------------------------------
+   !!                              PRIVATE SUBROUTINES 
+   !!-----------------------------------------------------------------------------------------------------
+   !
+   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in
+      &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange coefficient 
+      !!
+      !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
+      !!
+      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas [m/s]
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl ! velocity in the losch top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd          ! drag coefficient
+      REAL(wp),                     INTENT(in   ) :: pke2         ! background velocity
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar
+      !!---------------------------------------------------------------------
+      !
+      ! compute ustar (AD15 eq. 27)
+      zustar(:,:) = SQRT( pCd(:,:) * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) ) * mskisf_cav(:,:)
+      !
+      ! Compute gammats
+      pgt(:,:) = zustar(:,:) * rn_gammat0
+      pgs(:,:) = zustar(:,:) * rn_gammas0
+      !
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
+      !
+   END SUBROUTINE gammats_vel
+
+   SUBROUTINE gammats_vel_stab( pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
+      &                                                                pgt  , pgs    )  ! ==>> out gammats [m/s]
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange coefficient 
+      !!
+      !! ** Method     : gamma is velocity dependent and stability dependent
+      !!
+      !! ** Reference  : Holland and Jenkins, 1999, JPO, p1787-1800
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pCd            ! drag coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: putbl, pvtbl   ! velocity in the losch top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl   ! tracer   in the losch top boundary layer
+      !!---------------------------------------------------------------------
+      INTEGER  :: ji, jj                     ! loop index
+      INTEGER  :: ikt                        ! local integer
+      REAL(wp) :: zdku, zdkv                 ! U, V shear 
+      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
+      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
+      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
+      REAL(wp) :: zhmax                      ! limitation of mol
+      REAL(wp) :: zetastar                   ! stability parameter
+      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence 
+      REAL(wp) :: zcoef                      ! temporary coef
+      REAL(wp) :: zdep
+      REAL(wp) :: zeps = 1.0e-20_wp    
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
+      REAL(wp), DIMENSION(2) :: zts, zab
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
+      !!---------------------------------------------------------------------
+      !
+      ! compute ustar
+      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )
+      !
+      ! output ustar
+      CALL iom_put('isfustar',zustar(:,:))
+      !
+      ! compute Pr and Sc number (eq ??)
+      zPr =   13.8_wp
+      zSc = 2432.0_wp
+      !
+      ! compute gamma mole (eq ??)
+      zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+      zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+      !
+      ! compute gamma
+      DO ji = 2, jpi
+         DO jj = 2, jpj
+            ikt = mikt(ji,jj)
+
+            IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
+               pgt = rn_gammat0
+               pgs = rn_gammas0
+            ELSE
+               ! compute Rc number (as done in zdfric.F90)
+!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
+               zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
+               !                                            ! shear of horizontal velocity
+               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
+                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
+                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+               !                                            ! richardson number (minimum value set to zero)
+               zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
+
+               ! compute bouyancy 
+               zts(jp_tem) = pttbl(ji,jj)
+               zts(jp_sal) = pstbl(ji,jj)
+               zdep        = gdepw_n(ji,jj,ikt)
+               !
+               CALL eos_rab( zts, zdep, zab )
+               !
+               ! compute length scale (Eq ??)
+               zbuofdep = grav * ( zab(jp_tem) * pqoce(ji,jj) - zab(jp_sal) * pqfwf(ji,jj) )
+               !
+               ! compute Monin Obukov Length
+               ! Maximum boundary layer depth (Eq ??)
+               zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
+               !
+               ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
+               zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+               !
+               ! compute eta* (stability parameter) (Eq ??)
+               zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
+               !
+               ! compute the sublayer thickness (Eq ??)
+               zhnu = 5 * znu / zustar(ji,jj)
+               !
+               ! compute gamma turb (Eq ??)
+               zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
+               &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+               !
+               ! compute gammats
+               pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+               pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
+            END IF
+         END DO
+      END DO
+
+   END SUBROUTINE gammats_vel_stab
+
+END MODULE isfcavgam
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcavmlt.F90	(revision 12747)
@@ -0,0 +1,311 @@
+MODULE isfcavmlt
+   !!======================================================================
+   !!                       ***  MODULE  isfcavmlt  ***
+   !! ice shelf module :  update surface ocean boundary condition under ice
+   !!                   shelves
+   !!======================================================================
+   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfcav_mlt    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
+   !!----------------------------------------------------------------------
+
+   USE isf_oce                  ! ice shelf
+   USE isftbl , ONLY: isf_tbl   ! ice shelf depth average
+   USE isfutils,ONLY: debug     ! debug subroutine
+
+   USE dom_oce                            ! ocean space and time domain
+   USE oce    , ONLY: tsn                 ! ocean dynamics and tracers
+   USE phycst , ONLY: rcp, rau0, rau0_rcp ! physical constants
+   USE eosbn2 , ONLY: eos_fzp             ! equation of state
+
+   USE in_out_manager              ! I/O manager
+   USE iom        , ONLY: iom_put  ! I/O library
+   USE fldread    , ONLY: fld_read, FLD, FLD_N !
+   USE lib_fortran, ONLY: glob_sum !
+   USE lib_mpp    , ONLY: ctl_stop !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   isfcav_mlt
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfcav_mlt(kt, pgt, pgs , pttbl, pstbl, &
+      &                           pqhc, pqoce, pqfwf  )
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt  ***
+      !!
+      !! ** Purpose    : compute or read ice shelf fwf/heat fluxes in the ice shelf cavity
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER, INTENT(in) :: kt
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs    ! gamma t and gamma s
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! top boundary layer tracer
+      !!---------------------------------------------------------------------
+      !
+      ! compute latent heat and melt (2d)
+      SELECT CASE ( cn_isfcav_mlt )
+      CASE ( 'spe' )   ! ice shelf melt specified (read input file, and heat fluxes derived from
+         CALL isfcav_mlt_spe( kt, pstbl,               &
+            &                  pqhc, pqoce, pqfwf  )
+      CASE ( '2eq' )   !  ISOMIP  formulation (2 equations) for volume flux (Hunter et al., 2006)
+         CALL isfcav_mlt_2eq( pgt, pttbl, pstbl,       &
+            &                  pqhc , pqoce, pqfwf )
+      CASE ( '3eq' )   ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+         CALL isfcav_mlt_3eq( pgt, pgs , pttbl, pstbl, &
+            &                  pqhc, pqoce, pqfwf  )
+      CASE ( 'oasis' ) ! fwf pass trough oasis
+         CALL isfcav_mlt_oasis( kt, pstbl,              &
+            &                   pqhc, pqoce, pqfwf  )
+      CASE DEFAULT
+         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfcav (should not see this)')
+      END SELECT
+      !
+      IF (ln_isfdebug) THEN
+         CALL debug( 'isfcav_mlt:', pqhc (:,:) )
+         CALL debug( 'isfcav_mlt:', pqoce(:,:) )
+         CALL debug( 'isfcav_mlt:', pqfwf(:,:) )
+      END IF
+      !
+   END SUBROUTINE isfcav_mlt
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfcav_mlt_spe(kt, pstbl,          &  ! <<== in
+      &                      pqhc , pqoce, pqfwf )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_spe  ***
+      !!
+      !! ** Purpose    : - read ice shelf melt from forcing file
+      !!                 - compute ocea-ice heat flux (assuming it is equal to latent heat)
+      !!                 - compute heat content flux
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
+      !!--------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
+      !!--------------------------------------------------------------------
+      !
+      ! Compute freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! read input file
+      CALL fld_read ( kt, 1, sf_isfcav_fwf )
+      !
+      ! define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)      ! fwf                ( >0 out)
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
+      !
+      ! output freezing point at the interface
+      CALL iom_put('isftfrz_cav', ztfrz )
+      !
+   END SUBROUTINE isfcav_mlt_spe
+
+   SUBROUTINE isfcav_mlt_2eq(pgt , pttbl, pstbl, &  ! <<== in
+      &                      pqhc, pqoce, pqfwf  )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_2eq  ***
+      !!
+      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using ISOMIP formulation (Hunter et al., 2006)
+      !!
+      !! ** Method     : The ice shelf melt latent heat is defined as being equal to the ocean/ice heat flux.
+      !!                 From this we can derived the fwf, ocean/ice heat flux and the heat content flux as being :
+      !!                   qfwf  = Gammat * Rau0 * Cp * ( Tw - Tfrz ) / Lf 
+      !!                   qhoce = qlat
+      !!                   qhc   = qfwf * Cp * Tfrz
+      !!
+      !! ** Reference  : Hunter,  J.  R.:  Specification  for  test  models  of  ice  shelf  cavities,  
+      !!                 Tech.  Rep.  June,  Antarctic  Climate  &  Ecosystems  Cooperative  Research  Centre,  available  at:  
+      !!                 http://staff.acecrc.org.au/~bkgalton/ISOMIP/test_cavities.pdf (last access: 21 July 2016), 2006.
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! hean content, ocean-ice heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt           ! temperature exchange coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl  ! temperature and salinity in top boundary layer
+      !!--------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
+      !!--------------------------------------------------------------------
+      !
+      ! Calculate freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! thermal driving
+      zthd (:,:) = ( pttbl(:,:) - ztfrz(:,:) ) * mskisf_cav(:,:)
+      !
+      ! compute ocean-ice heat flux and then derive fwf assuming that ocean heat flux equal latent heat
+      pqfwf(:,:) = - pgt(:,:) * rau0_rcp * zthd(:,:) / rLfusisf    ! fresh water flux  ( > 0 out )
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf                         ! ocea-ice flux     ( > 0 out )
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp                 ! heat content flux ( > 0 out )
+      !
+      ! output thermal driving and freezinpoint at the ice shelf interface
+      CALL iom_put('isfthermald_cav', zthd )
+      CALL iom_put('isftfrz_cav'    , ztfrz )
+      !
+   END SUBROUTINE isfcav_mlt_2eq
+
+   SUBROUTINE isfcav_mlt_3eq(pgt, pgs , pttbl, pstbl, &  ! <<== in
+      &                           pqhc, pqoce, pqfwf  )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!
+      !!                          ***  ROUTINE isfcav_mlt_3eq  ***
+      !!
+      !! ** Purpose    : Compute ice shelf fwf/heqt fluxes using the 3 equation formulation 
+      !!
+      !! ** Method     : The melt rate is determined considering the heat balance, the salt balance
+      !!                 at the phase change interface and a linearisation of the equation of state.
+      !!
+      !! ** Reference  : - Holland, D. M. and Jenkins, A.,
+      !!                   Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf,
+      !!                   J. Phys. Oceanogr., 29, 1999.
+      !!                 - Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, 
+      !!                   R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, 
+      !!                   P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.:
+      !!                   Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: 
+      !!                   MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), 
+      !!                   Geosci. Model Dev., 9, 2471-2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016. 
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pgt  , pgs          ! heat/salt exchange coeficient
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pttbl, pstbl        ! mean temperature and salinity in top boundary layer
+      !!--------------------------------------------------------------------
+      REAL(wp) :: zeps1,zeps2,zeps3,zeps4,zeps6,zeps7       ! dummy local scalar for quadratic equation resolution
+      REAL(wp) :: zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac     ! dummy local scalar for quadratic equation resolution
+      REAL(wp) :: zeps = 1.e-20
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz         ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj) :: zqcon         ! conductive flux through the ice shelf
+      REAL(wp), DIMENSION(jpi,jpj) :: zthd          ! thermal driving
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      !!--------------------------------------------------------------------
+      !
+      ! compute upward heat flux zhtflx and upward water flux zwflx
+      ! Resolution of a 3d equation from equation 24, 25 and 26 (note conduction through the ice has been added to Eq 24)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            !
+            ! compute coeficient to solve the 2nd order equation
+            zeps1 = rau0_rcp * pgt(ji,jj)
+            zeps2 = rLfusisf * rau0 * pgs(ji,jj)
+            zeps3 = rhoisf * rcpisf * rkappa / MAX(risfdep(ji,jj),zeps)
+            zeps4 = risf_lamb2 + risf_lamb3 * risfdep(ji,jj)
+            zeps6 = zeps4 - pttbl(ji,jj)
+            zeps7 = zeps4 - rtsurf
+            !
+            ! solve the 2nd order equation to find zsfrz
+            zaqe  = risf_lamb1 * (zeps1 + zeps3)
+            zaqer = 0.5_wp / MIN(zaqe,-zeps)
+            zbqe  = zeps1 * zeps6 + zeps3 * zeps7 - zeps2
+            zcqe  = zeps2 * pstbl(ji,jj)
+            zdis  = zbqe * zbqe - 4.0_wp * zaqe * zcqe               
+            !
+            ! Presumably zdis can never be negative because gammas is very small compared to gammat
+            zsfrz=(-zbqe - SQRT(zdis)) * zaqer
+            IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe + SQRT(zdis)) * zaqer  ! check this if this if is needed
+            !
+            ! compute t freeze (eq. 25)
+            ztfrz(ji,jj) = zeps4 + risf_lamb1 * zsfrz
+            !
+            ! thermal driving
+            zthd(ji,jj) = ( pttbl(ji,jj) - ztfrz(ji,jj) )
+            !
+            ! compute the upward water and heat flux (eq. 24 and eq. 26)
+            pqfwf(ji,jj) = rau0     * pgs(ji,jj) * ( zsfrz - pstbl(ji,jj) ) / MAX(zsfrz,zeps) ! fresh water flux    (> 0 out)
+            pqoce(ji,jj) = rau0_rcp * pgt(ji,jj) * zthd (ji,jj)                               ! ocean-ice heat flux (> 0 out)
+            pqhc (ji,jj) = rcp      * pqfwf(ji,jj) * ztfrz(ji,jj)                             ! heat content   flux (> 0 out)
+            !
+            zqcon(ji,jj) = zeps3 * ( ztfrz(ji,jj) - rtsurf )
+            !
+         END DO
+      END DO
+      !
+      ! output conductive heat flux through the ice
+      CALL iom_put('qconisf', zqcon(:,:) * mskisf_cav(:,:) )
+      !
+      ! output thermal driving and freezing point at the interface
+      CALL iom_put('isfthermald_cav', zthd (:,:) * mskisf_cav(:,:) )
+      CALL iom_put('isftfrz_cav'    , ztfrz(:,:) * mskisf_cav(:,:) )
+      !
+   END SUBROUTINE isfcav_mlt_3eq
+
+   SUBROUTINE isfcav_mlt_oasis(kt, pstbl,          &  ! <<== in
+      &                        pqhc , pqoce, pqfwf )  ! ==>> out
+      !!----------------------------------------------------------------------
+      !!                          ***  ROUTINE isfcav_mlt_oasis  ***
+      !!
+      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
+      !!
+      !! ** Purpose    : - read ice shelf melt from forcing file => pattern
+      !!                 - total amount of fwf is given by sbccpl (fwfisf_cpl)
+      !!                 - scale fwf and compute heat fluxes
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pstbl               ! salinity in tbl
+      !!--------------------------------------------------------------------
+      REAL(wp)                     :: zfwf_fld, zfwf_oasis               ! total fwf in the forcing fields (pattern) and from the oasis interface (amount)
+      REAL(wp), DIMENSION(jpi,jpj) :: ztfrz                              ! tbl freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj) :: zfwf                               ! 2d fwf map after scaling
+      !!--------------------------------------------------------------------
+      !
+      ! Calculate freezing temperature
+      CALL eos_fzp( pstbl(:,:), ztfrz(:,:), risfdep(:,:) )
+      !
+      ! read input file
+      CALL fld_read ( kt, 1, sf_isfcav_fwf )
+      !
+      ! ice shelf 2d map
+      zfwf(:,:) = - sf_isfcav_fwf(1)%fnow(:,:,1)
+      !
+      ! compute glob sum from input file
+      ! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
+      zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
+      !
+      ! compute glob sum from atm->oce ice shelf fwf
+      ! (PM) should consider delay sum as in fwb (1 time step offset if I well understood)
+      zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
+      !
+      ! scale fwf
+      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
+      ! 
+      ! define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out)
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out)
+      !
+      CALL iom_put('isftfrz_cav', ztfrz * mskisf_cav(:,:) )
+      !
+   END SUBROUTINE isfcav_mlt_oasis
+
+END MODULE isfcavmlt
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfcpl.F90	(revision 12747)
@@ -0,0 +1,747 @@
+MODULE isfcpl
+   !!======================================================================
+   !!                       ***  MODULE  isfcpl  ***
+   !!
+   !! iceshelf coupling module : module managing the coupling between NEMO and an ice sheet model
+   !!
+   !!======================================================================
+   !! History :  4.1  !  2019-07  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfrst : read/write iceshelf variables in/from restart
+   !!----------------------------------------------------------------------
+   USE isf_oce                          ! ice shelf variable
+   USE isfutils, ONLY : debug
+   USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
+   USE domvvl  , ONLY: dom_vvl_zgr      ! vertical scale factor interpolation
+   !
+   USE oce            ! ocean dynamics and tracers
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   !
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isfcpl_rst_write, isfcpl_init                    ! iceshelf restart read and write 
+   PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons  ! iceshelf correction for ssh, tra, dyn and conservation 
+
+   TYPE isfcons
+      INTEGER :: ii     ! i global
+      INTEGER :: jj     ! j global
+      INTEGER :: kk     ! k level
+      REAL(wp):: dvol   ! volume increment
+      REAL(wp):: dsal   ! salt increment
+      REAL(wp):: dtem   ! heat increment
+      REAL(wp):: lon    ! lon
+      REAL(wp):: lat    ! lat
+      INTEGER :: ngb    ! 0/1 (valid location or not (ie on halo or no neigbourg))
+   END TYPE
+   !
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE isfcpl_init()
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE iscpl_init  ***
+      !! 
+      !! ** Purpose : correct ocean state for new wet cell and horizontal divergence 
+      !!              correction for the dynamical adjustement
+      !!
+      !! ** Action : - compute ssh on new wet cell
+      !!             - compute T/S on new wet cell
+      !!             - compute horizontal divergence correction as a volume flux
+      !!             - compute the T/S/vol correction increment to keep trend to 0
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER :: id
+      !!----------------------------------------------------------------------
+      !
+      ! start on an euler time step
+      neuler = 0
+      ! 
+      ! allocation and initialisation to 0
+      CALL isf_alloc_cpl()
+      !
+      ! check presence of variable needed for coupling
+      ! iom_varid return 0 if not found
+      id = 1
+      id = id * iom_varid(numror, 'ssmask', ldstop = .false.)
+      id = id * iom_varid(numror, 'tmask' , ldstop = .false.)
+      id = id * iom_varid(numror, 'e3t_n' , ldstop = .false.)
+      id = id * iom_varid(numror, 'e3u_n' , ldstop = .false.)
+      id = id * iom_varid(numror, 'e3v_n' , ldstop = .false.)
+      IF(lwp) WRITE(numout,*) ' isfcpl_init:', id
+      IF (id == 0) THEN
+         IF(lwp) WRITE(numout,*) ' isfcpl_init: restart variables for ice sheet coupling are missing, skip coupling for this leg ' 
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
+         IF(lwp) WRITE(numout,*) ''
+      ELSE
+         ! extrapolation ssh
+         CALL isfcpl_ssh()
+         !
+         ! extrapolation tracer properties
+         CALL isfcpl_tra()
+         !
+         ! correction of the horizontal divergence and associated temp. and salt content flux
+         ! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
+         !           - decide how to manage thickness level change in conservation
+         CALL isfcpl_vol()
+         !
+         ! apply the 'conservation' method
+         IF ( ln_isfcpl_cons ) CALL isfcpl_cons()
+         !
+      END IF
+      !
+      ! mask velocity properly (mask used in restart not compatible with new mask)
+      un(:,:,:) = un(:,:,:) * umask(:,:,:)
+      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+      !
+      ! all before fields set to now values
+      tsb  (:,:,:,:) = tsn  (:,:,:,:)
+      ub   (:,:,:)   = un   (:,:,:)
+      vb   (:,:,:)   = vn   (:,:,:)
+      sshb (:,:)     = sshn (:,:)
+      e3t_b(:,:,:)   = e3t_n(:,:,:)
+ 
+      ! prepare writing restart
+      IF( lwxios ) THEN
+         CALL iom_set_rstw_var_active('ssmask')
+         CALL iom_set_rstw_var_active('tmask')
+         CALL iom_set_rstw_var_active('e3t_n')
+         CALL iom_set_rstw_var_active('e3u_n')
+         CALL iom_set_rstw_var_active('e3v_n')
+      END IF
+      !
+   END SUBROUTINE isfcpl_init
+   ! 
+   SUBROUTINE isfcpl_rst_write(kt)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE iscpl_rst_write  ***
+      !! 
+      !! ** Purpose : write icesheet coupling variables in restart
+      !!
+      !!-------------------------- IN  --------------------------------------
+      INTEGER, INTENT(in) :: kt
+      !!----------------------------------------------------------------------
+      !
+      IF( lwxios ) CALL iom_swap( cwxios_context )
+      CALL iom_rstput( kt, nitrst, numrow, 'tmask'  , tmask , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask, ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3t_n'  , e3t_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3u_n'  , e3u_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'e3v_n'  , e3v_n , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', gdepw_n , ldxios = lwxios )
+      IF( lwxios ) CALL iom_swap( cxios_context )
+      !
+   END SUBROUTINE isfcpl_rst_write
+
+   SUBROUTINE isfcpl_ssh()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_ssh  ***
+      !! 
+      !! ** Purpose :   basic guess of ssh in new wet cell
+      !! 
+      !! ** Method  :   basic extrapolation from neigbourg cells
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: ji, jj, jd, jk      !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1
+      !!
+      REAL(wp):: zsummsk
+      REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
+      !!----------------------------------------------------------------------
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'ssmask'  , zssmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+
+      ! compute new ssh if we open a full water column 
+      ! rude average of the closest neigbourgs (e1e2t not taking into account)
+      !
+      zssh(:,:)     = sshn(:,:)
+      zssmask0(:,:) = zssmask_b(:,:)
+      !
+      DO jd = 1, nn_drown
+         !
+         zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
+         DO jj = 2,jpj-1
+            DO ji = 2,jpi-1
+               jip1=ji+1; jim1=ji-1;
+               jjp1=jj+1; jjm1=jj-1;
+               !
+               zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
+               !
+               IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
+                  sshn(ji,jj)=( zssh(jip1,jj)*zssmask0(jip1,jj)     &
+                  &           + zssh(jim1,jj)*zssmask0(jim1,jj)     &
+                  &           + zssh(ji,jjp1)*zssmask0(ji,jjp1)     &
+                  &           + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
+                  zssmask_b(ji,jj) = 1._wp
+               ENDIF
+            END DO
+         END DO
+         !
+         zssh(:,:) = sshn(:,:)
+         zssmask0(:,:) = zssmask_b(:,:)
+         !
+         CALL lbc_lnk_multi( 'iscplrst', zssh, 'T', 1., zssmask0, 'T', 1. )
+         !
+      END DO
+      !
+      ! update sshn
+      sshn(:,:) = zssh(:,:) * ssmask(:,:)
+      !
+      sshb(:,:) = sshn(:,:)
+      !
+      IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',sshn(:,:))
+      !
+      ! recompute the vertical scale factor, depth and water thickness
+      IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from sshn (new wet cell)'
+      IF(lwp) write(numout,*) '~~~~~~~~~~~'
+      DO jk = 1, jpk
+         e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+             &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+             &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
+      END DO
+      e3t_b(:,:,:) = e3t_n(:,:,:)
+      CALL dom_vvl_zgr()
+      !
+   END SUBROUTINE isfcpl_ssh
+
+   SUBROUTINE isfcpl_tra()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_tra  ***
+      !! 
+      !! ** Purpose :   compute new tn, sn in case of evolving geometry of ice shelves 
+      !! 
+      !! ** Method  :   tn, sn : basic extrapolation from neigbourg cells
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
+      !REAL(wp), DIMENSION(:,:,:  ), INTENT(in ) :: pdepw_b                         !! depth w before
+      !!
+      INTEGER :: ji, jj, jk, jd          !! loop index
+      INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
+      !!
+      REAL(wp):: zsummsk
+      REAL(wp):: zdz, zdzm1, zdzp1
+      !!
+      REAL(wp), DIMENSION(jpi,jpj)          :: zdmask 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask0, zwmaskn
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: ztmask1, zwmaskb, ztmp3d
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
+      !!----------------------------------------------------------------------
+      ! 
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+      !CALL iom_get( numror, jpdom_autoglo, 'wmask'  , zwmask_b, ldxios = lrxios   ) ! need to extrapolate T/S
+      !CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:), ldxios = lrxios ) ! need to interpol vertical profile (vvl)
+      !
+      ! 
+      ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+      !PM: Is this IF needed since change to VVL by default
+      !bugged : to be corrected (PM)
+      ! back up original t/s/mask
+      !tsb (:,:,:,:) = tsn(:,:,:,:)
+      !
+     ! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
+
+!      IF (.NOT.ln_linssh) THEN
+!         DO jk = 2,jpk-1
+!            DO jj = 1,jpj
+!               DO ji = 1,jpi
+!                  IF (wmask(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ztmask_b(ji,jj,1)==0._wp) ) THEN
+!
+!                     !compute weight
+!                     zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw_n(ji,jj,jk+1))
+!                     zdzm1 = MAX(0._wp,gdepw_n(ji,jj,jk  ) - pdepw_b(ji,jj,jk  ))
+!                     zdz   = e3t_n(ji,jj,jk) - zdzp1 - zdzm1 ! if isf : e3t = gdepw_n(ji,jj,jk+1)- gdepw_n(ji,jj,jk)
+!
+!                     IF (zdz .LT. 0._wp) THEN
+!                        CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
+!                     END IF
+!
+!                     tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
+!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_tem) &
+!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
+!
+!                     tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
+!                        &                   + zdz  *tsb(ji,jj,jk  ,jp_sal) &
+!                        &                   + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
+!
+!                  END IF
+!               END DO
+!            END DO
+!         END DO
+!      END IF
+
+      zts0(:,:,:,:)  = tsn(:,:,:,:)
+      ztmask0(:,:,:) = ztmask_b(:,:,:)
+      ztmask1(:,:,:) = ztmask_b(:,:,:)
+      !
+      ! iterate the extrapolation processes nn_drown times
+      DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
+         DO jk = 1,jpk-1
+            !
+            ! define new wet cell
+            zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
+            !
+            DO jj = 2,jpj-1
+               DO ji = 2,jpi-1
+                  jip1=ji+1; jim1=ji-1;
+                  jjp1=jj+1; jjm1=jj-1;
+                  !
+                  ! check if a wet neigbourg cell is present
+                  zsummsk = ztmask0(jip1,jj  ,jk) + ztmask0(jim1,jj  ,jk) &
+                          + ztmask0(ji  ,jjp1,jk) + ztmask0(ji  ,jjm1,jk)
+                  !
+                  ! if neigbourg wet cell available at the same level
+                  IF ( zdmask(ji,jj) == 1._wp  .AND. zsummsk /= 0._wp ) THEN
+                     !
+                     ! horizontal basic extrapolation
+                     tsn(ji,jj,jk,1)=( zts0(jip1,jj  ,jk,1) * ztmask0(jip1,jj  ,jk) &
+                     &               + zts0(jim1,jj  ,jk,1) * ztmask0(jim1,jj  ,jk) &
+                     &               + zts0(ji  ,jjp1,jk,1) * ztmask0(ji  ,jjp1,jk) &
+                     &               + zts0(ji  ,jjm1,jk,1) * ztmask0(ji  ,jjm1,jk) ) / zsummsk
+                     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) ) / zsummsk
+                     !
+                     ! update mask for next pass
+                     ztmask1(ji,jj,jk)=1
+                     !
+                  ! in case no neigbourg wet cell available at the same level
+                  ! check if a wet cell is available below
+                  ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
+                     !
+                     ! vertical extrapolation if horizontal extrapolation failed
+                     jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
+                     !
+                     ! check if a wet neigbourg cell is present
+                     zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
+                     IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
+                        tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1)     &
+                        &               + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
+                        tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1)     &
+                        &               + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
+                        !
+                        ! update mask for next pass
+                        ztmask1(ji,jj,jk)=1._wp
+                     END IF
+                  END IF
+               END DO
+            END DO
+         END DO
+         !
+         ! update temperature and salinity and mask
+         zts0(:,:,:,:)  = tsn(:,:,:,:)
+         ztmask0(:,:,:) = ztmask1(:,:,:)
+         !
+         CALL lbc_lnk_multi( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1., zts0(:,:,:,jp_sal), 'T', 1., ztmask0, 'T', 1.)
+         !
+      END DO  ! nn_drown
+      !
+      ! mask new tsn field
+      tsn(:,:,:,jp_tem) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
+      tsn(:,:,:,jp_sal) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
+      !
+      ! sanity check
+      ! -----------------------------------------------------------------------------------------
+      ! case we open a cell but no neigbour cells available to get an estimate of T and S
+      DO jk = 1,jpk-1
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               IF (tmask(ji,jj,jk) == 1._wp .AND. tsn(ji,jj,jk,2) == 0._wp)              &
+                  &   CALL ctl_stop('STOP', 'failing to fill all new weet cell,     &
+                  &                          try increase nn_drown or activate XXXX &
+                  &                         in your domain cfg computation'         )
+            END DO
+         END DO
+      END DO
+      ! 
+   END SUBROUTINE isfcpl_tra
+
+   SUBROUTINE isfcpl_vol()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_vol  ***
+      !! 
+      !! ** Purpose : compute the correction of the local divergence to apply  
+      !!              during the first time step after the coupling.
+      !!
+      !! ** Method  : - compute horizontal vol div. before/after coupling
+      !!              - compute vertical input
+      !!              - compute correction
+      !!                
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER :: ji, jj, jk 
+      INTEGER :: ikb, ikt
+      !!
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln  ! vol flux div.         before/after coupling
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b  ! vertical scale factor before/after coupling
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b        ! mask                  before       coupling
+      !!----------------------------------------------------------------------
+      !
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b, ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'e3u_n'  , ze3u_b  , ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'e3v_n'  , ze3v_b  , ldxios = lrxios )
+      !
+      ! 1.0: compute horizontal volume flux divergence difference before-after coupling
+      !
+      DO jk = 1, jpk                                 ! Horizontal slab
+         ! 1.1: get volume flux before coupling (>0 out)
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zqvolb(ji,jj,jk) =  (   e2u(ji,jj) * ze3u_b(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * ze3u_b(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                  &                  + e1v(ji,jj) * ze3v_b(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * ze3v_b(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+                  &                * ztmask_b(ji,jj,jk)
+            END DO
+         ENDDO
+         !
+         ! 1.2: get volume flux after coupling (>0 out)
+         ! properly mask velocity 
+         ! (velocity are still mask with old mask at this stage)
+         un(:,:,jk) = un(:,:,jk) * umask(:,:,jk)
+         vn(:,:,jk) = vn(:,:,jk) * vmask(:,:,jk)
+         ! compute volume flux divergence after coupling
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zqvoln(ji,jj,jk) = (   e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  ) * e3u_n(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
+                  &                 + e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1) * e3v_n(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  ) &
+                  &               * tmask(ji,jj,jk)
+            END DO
+         ENDDO
+         !
+         ! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
+         !      correction to add is _b - _n
+         risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
+      END DO
+      !
+      ! 2.0: include the contribution of the vertical velocity in the volume flux correction
+      !
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            !
+            ikt = mikt(ji,jj)
+            IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
+               risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1))  ! test sign
+            ENDIF
+            !
+         END DO
+      ENDDO
+      !
+      CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1. )
+      !
+      ! 3.0: set total correction (div, tra, ssh)
+      !
+      ! 3.1: mask volume flux divergence correction
+      risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
+      !
+      ! 3.2: get 3d tra increment to apply at the first time step
+      ! temperature and salt content flux computed using local tsn 
+      ! (very simple advection scheme)
+      ! (>0 out)
+      risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_tem)
+      risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * tsn(:,:,:,jp_sal)
+      !
+      ! 3.3: ssh correction (for dynspg_ts)
+      risfcpl_ssh(:,:) = 0.0
+      DO jk = 1,jpk
+         risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
+      END DO
+
+   END SUBROUTINE isfcpl_vol
+
+   SUBROUTINE isfcpl_cons()
+      !!---------------------------------------------------------------------- 
+      !!                   ***  ROUTINE iscpl_cons  ***
+      !! 
+      !! ** Purpose :   compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
+      !!                removed or added during the coupling processes (wet or dry new cell)
+      !! 
+      !! ** Method  :   - compare volume/heat/salt before and after
+      !!                - look for the closest wet cells (share amoung neigbourgs if there are)
+      !!                - build the correction increment to applied at each time step
+      !!                
+      !!----------------------------------------------------------------------
+      !
+      TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
+      !
+      INTEGER  ::   ji   , jj  , jk  , jproc          ! loop index
+      INTEGER  ::   jip1 , jim1, jjp1, jjm1           ! dummy indices
+      INTEGER  ::   iig  , ijg, ik                    ! dummy indices
+      INTEGER  ::   jisf                              ! start, end and current position in the increment array
+      INTEGER  ::   ingb, ifind                       ! 0/1 target found or need to be found 
+      INTEGER  ::   nisfl_area                        ! global number of cell concerned by the wet->dry case 
+      INTEGER, DIMENSION(jpnij) :: nisfl              ! local  number of cell concerned by the wet->dry case
+      !
+      REAL(wp) ::   z1_sum, z1_rdtiscpl
+      REAL(wp) ::   zdtem, zdsal, zdvol, zratio       ! tem, sal, vol increment
+      REAL(wp) ::   zlon , zlat                       ! target location  
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b    ! mask before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b      ! scale factor before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b      ! scale factor before
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b      ! scale factor before
+      !!----------------------------------------------------------------------
+
+      !==============================================================================
+      ! 1.0: initialisation
+      !==============================================================================
+
+      ! get restart variable
+      CALL iom_get( numror, jpdom_autoglo, 'tmask'  , ztmask_b(:,:,:), ldxios = lrxios   ) ! need to extrapolate T/S
+      CALL iom_get( numror, jpdom_autoglo, 'e3t_n'  , ze3t_b(:,:,:)  , ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'tn'     , zt_b(:,:,:)    , ldxios = lrxios )
+      CALL iom_get( numror, jpdom_autoglo, 'sn'     , zs_b(:,:,:)    , ldxios = lrxios )
+
+      ! compute run length
+      nstp_iscpl  = nitend - nit000 + 1
+      rdt_iscpl   = nstp_iscpl * rn_rdt
+      z1_rdtiscpl = 1._wp / rdt_iscpl 
+
+      IF (lwp) WRITE(numout,*) '            nb of stp for cons  = ', nstp_iscpl
+      IF (lwp) WRITE(numout,*) '            coupling time step  = ', rdt_iscpl
+
+      ! initialisation correction
+      risfcpl_cons_vol = 0.0
+      risfcpl_cons_ssh = 0.0
+      risfcpl_cons_tsc = 0.0
+
+      !==============================================================================
+      ! 2.0: diagnose the heat, salt and volume input and compute the correction variable
+      !      for case where we wet a cell or cell still wet (no change in cell status)
+      !==============================================================================
+
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+
+               ! volume diff
+               zdvol = e3t_n(ji,jj,jk) * tmask(ji,jj,jk) - ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+
+               ! heat diff
+               zdtem = tsn (ji,jj,jk,jp_tem) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                     - zt_b(ji,jj,jk)        * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+
+               ! salt diff
+               zdsal = tsn(ji,jj,jk,jp_sal) *  e3t_n(ji,jj,jk) *  tmask  (ji,jj,jk)   &
+                     - zs_b(ji,jj,jk)       * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
+            
+               ! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
+               ! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
+               risfcpl_cons_vol(ji,jj,jk)        = (   zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk)        ) * z1_rdtiscpl
+               risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
+               risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
+
+            END DO
+         END DO
+      END DO
+      !
+      !==============================================================================
+      ! 3.0: diagnose the heat, salt and volume input and compute the correction variable
+      !      for case where we close a cell
+      !==============================================================================
+      !
+      ! compute the total number of point receiving a correction increment for each processor
+      ! local
+      nisfl(:)=0
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+               jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
+               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      ! global 
+      CALL mpp_sum('isfcpl',nisfl  )
+      !
+      ! allocate list of point receiving correction
+      ALLOCATE(zisfpts(nisfl(narea)))
+      !
+      zisfpts(:) = isfcons(0,0,0,-HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), -HUGE(1.0), 0)
+      !
+      ! start computing the correction and fill zisfpts
+      ! local
+      jisf = 0
+      DO jk = 1,jpk-1
+         DO jj = nldj,nlej
+            DO ji = nldi,nlei
+               IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
+
+                  jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
+
+                  zdvol = risfcpl_cons_vol(ji,jj,jk       )
+                  zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
+                  zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
+
+                  IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
+                     ! spread correction amoung neigbourg wet cells (horizontal direction first)
+                     ! as it is a rude correction corner and lateral cell have the same weight
+                     !
+                     z1_sum =  1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
+                     !
+                     ! lateral cells
+                     IF (tmask(jip1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jj  ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj  , jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(ji  ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(ji  ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji  , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     !
+                     ! corner  cells
+                     IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
+                     !
+                  ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
+                     ! spread correction amoung neigbourg wet cells (vertical direction)
+                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk+1, zdvol, zdsal, zdtem, 1., 0)
+                  ELSE
+                     ! need to find where to put correction in later on
+                     CALL update_isfpts(zisfpts, jisf, ji  , jj  , jk  , zdvol, zdsal, zdtem, 1., 1)
+                  END IF
+               END IF
+            END DO
+         END DO
+      END DO
+      !
+      ! share data among all processes because for some point we need to find the closest wet point (could be on other process)
+      DO jproc=1,jpnij
+         ! 
+         ! share total number of isf point treated for proc jproc
+         IF (jproc==narea) THEN
+            nisfl_area=nisfl(jproc)
+         ELSE
+            nisfl_area=0
+         END IF
+         CALL mpp_max('isfcpl',nisfl_area)
+         !
+         DO jisf = 1,nisfl_area
+            !
+            IF (jproc==narea) THEN
+               ! indices (conversion to global indices and sharing)
+               iig = zisfpts(jisf)%ii       ; ijg = zisfpts(jisf)%jj       ; ik = zisfpts(jisf)%kk
+               !
+               ! data
+               zdvol = zisfpts(jisf)%dvol   ; zdsal = zisfpts(jisf)%dsal   ; zdtem = zisfpts(jisf)%dtem
+               !
+               ! location
+               zlat = zisfpts(jisf)%lat     ; zlon = zisfpts(jisf)%lon
+               !
+               ! find flag
+               ingb = zisfpts(jisf)%ngb
+            ELSE
+               iig  =0   ; ijg  =0   ; ik   =0  
+               zdvol=-HUGE(1.0) ; zdsal=-HUGE(1.0) ; zdtem=-HUGE(1.0)
+               zlat =-HUGE(1.0) ; zlon =-HUGE(1.0)   
+               ingb = 0
+            END IF
+            !
+            ! share data (need synchronisation of data as get_correction call a global com)
+            CALL mpp_max('isfcpl',iig)   ; CALL mpp_max('isfcpl',ijg)   ; CALL mpp_max('isfcpl',ik)
+            CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
+            CALL mpp_max('isfcpl',zlat)  ; CALL mpp_max('isfcpl',zlon)
+            CALL mpp_max('isfcpl',ingb)
+            !
+            ! fill the 3d correction array
+            CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
+         END DO
+      END DO
+      !
+      !==============================================================================
+      ! 4.0: finalisation and compute ssh equivalent of the volume correction
+      !==============================================================================
+      !
+      ! mask (>0 out)
+      risfcpl_cons_vol(:,:,:       ) = risfcpl_cons_vol(:,:,:       ) * tmask(:,:,:)
+      risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
+      risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
+      !
+      ! add lbclnk
+      CALL lbc_lnk_multi( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1., risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1., &
+         &                            risfcpl_cons_vol(:,:,:)       , 'T', 1.)
+      !
+      ! ssh correction (for dynspg_ts)
+      DO jk = 1,jpk
+         risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
+      END DO
+      risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
+      !
+   END SUBROUTINE isfcpl_cons
+   !
+   SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE update_isfpts  ***
+      !!
+      !! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
+      !!
+      !! ** Action  : update the list of point
+      !!
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
+      INTEGER,                     INTENT(inout) :: kpts
+      !!----------------------------------------------------------------------
+      INTEGER,      INTENT(in   )           :: ki, kj, kk                  !    target location (kfind=0) 
+      !                                                                    ! or source location (kfind=1)
+      INTEGER,      INTENT(in   ), OPTIONAL :: kfind                       ! 0  target cell already found
+      !                                                                    ! 1  target to be determined
+      REAL(wp),     INTENT(in   )           :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment 
+      !                                                                    ! and ratio in case increment span over multiple cells.
+      !!----------------------------------------------------------------------
+      INTEGER :: ifind
+      !!----------------------------------------------------------------------
+      ! 
+      ! increment position
+      kpts = kpts + 1
+      !
+      ! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
+      IF ( PRESENT(kfind) ) THEN
+         ifind = kfind
+      ELSE
+         ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
+      END IF
+      !
+      ! update isfpts structure
+      sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
+      !
+   END SUBROUTINE update_isfpts
+   !
+   SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE get_correction  ***
+      !!
+      !! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
+      !!             - Scale the correction depending of pratio (case where multiple wet neigbourgs)
+      !!             - Fill the correction array
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER , INTENT(in) :: ki, kj, kk, kfind        ! target point indices
+      REAL(wp), INTENT(in) :: plon, plat               ! target point lon/lat
+      REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
+      !!----------------------------------------------------------------------
+      INTEGER :: jj, ji, iig, ijg
+      !!----------------------------------------------------------------------
+      !
+      ! define global indice of correction location
+      iig = ki ; ijg = kj
+      IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
+      !
+      ! fill the correction array
+      DO jj = mj0(ijg),mj1(ijg)
+         DO ji = mi0(iig),mi1(iig)
+            ! correct the vol_flx and corresponding heat/salt flx in the closest cell
+            risfcpl_cons_vol(ji,jj,kk)        =  risfcpl_cons_vol(ji,jj,kk       ) + pvolinc
+            risfcpl_cons_tsc(ji,jj,kk,jp_sal) =  risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
+            risfcpl_cons_tsc(ji,jj,kk,jp_tem) =  risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
+         END DO
+      END DO
+
+   END SUBROUTINE get_correction
+
+END MODULE isfcpl
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfdiags.F90	(revision 12747)
@@ -0,0 +1,113 @@
+MODULE isfdiags
+   !!======================================================================
+   !!                       ***  MODULE  isfdiags  ***
+   !! ice shelf diagnostics module :  manage the 2d and 3d flux outputs from the ice shelf module
+   !!======================================================================
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_isf       : update sbc under ice shelf
+   !!----------------------------------------------------------------------
+
+   USE in_out_manager ! I/O manager
+   USE dom_oce
+   USE isf_oce        ! ice shelf variable
+   USE iom            ! 
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC   isf_diags_flx
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE isf_diags_flx(ktop, kbot, phtbl, pfrac, cdisf, pqfwf, pqoce, pqlat, pqhc)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_diags_flx ***
+      !!
+      !! ** Purpose : manage the 2d and 3d flux outputs of the ice shelf module
+      !!              fwf, latent heat flux, heat content flux, oce->ice heat flux
+      !!
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      !!-------------------------- IN  -------------------------------------
+      INTEGER , DIMENSION(jpi,jpj), INTENT(in) :: ktop , kbot               ! top and bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac              ! thickness of the tbl and fraction of last cell affected by the tbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqfwf, pqoce, pqlat, pqhc ! 2d var to map in 3d
+      CHARACTER(LEN=3), INTENT(in) :: cdisf                                 ! parametrisation or interactive melt
+      !!---------------------------------------------------------------------
+      CHARACTER(LEN=256) :: cvarqfwf  , cvarqoce  , cvarqlat  , cvarqhc
+      CHARACTER(LEN=256) :: cvarqfwf3d, cvarqoce3d, cvarqlat3d, cvarqhc3d
+      !!---------------------------------------------------------------------
+      !
+      ! output melt
+      cvarqfwf = 'fwfisf_'//cdisf  ; cvarqfwf3d = 'fwfisf3d_'//cdisf
+      cvarqoce = 'qoceisf_'//cdisf ; cvarqoce3d = 'qoceisf3d_'//cdisf
+      cvarqlat = 'qlatisf_'//cdisf ; cvarqlat3d = 'qlatisf3d_'//cdisf 
+      cvarqhc  = 'qhcisf_'//cdisf  ; cvarqhc3d  = 'qhcisf3d_'//cdisf
+      !
+      ! output 2d melt rate, latent heat and heat content flux from the injected water
+      CALL iom_put( TRIM(cvarqfwf), pqfwf(:,:) )   ! mass         flux ( >0 out )
+      CALL iom_put( TRIM(cvarqoce), pqoce(:,:) )   ! oce to ice   flux ( >0 out )
+      CALL iom_put( TRIM(cvarqlat), pqlat(:,:) )   ! latent heat  flux ( >0 out )
+      CALL iom_put( TRIM(cvarqhc) , pqhc (:,:) )   ! heat content flux ( >0 out )
+      !
+      ! output 3d Diagnostics
+      IF ( iom_use( TRIM(cvarqfwf3d) ) ) CALL isf_diags_2dto3d( ktop, kbot, phtbl, pfrac, TRIM(cvarqfwf3d) , pqfwf(:,:))
+      IF ( iom_use( TRIM(cvarqoce3d) ) ) CALL isf_diags_2dto3d( ktop, kbot, phtbl, pfrac, TRIM(cvarqoce3d) , pqoce(:,:))
+      IF ( iom_use( TRIM(cvarqlat3d) ) ) CALL isf_diags_2dto3d( ktop, kbot, phtbl, pfrac, TRIM(cvarqlat3d) , pqoce(:,:))
+      IF ( iom_use( TRIM(cvarqhc3d)  ) ) CALL isf_diags_2dto3d( ktop, kbot, phtbl, pfrac, TRIM(cvarqhc3d)  , pqhc (:,:))
+      !
+   END SUBROUTINE
+
+   SUBROUTINE isf_diags_2dto3d(ktop, kbot, phtbl, pfrac, cdvar, pvar2d)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_diags_2dto3d ***
+      !!
+      !! ** Purpose : compute the 3d flux outputs as they are injected into NEMO 
+      !!              (ie uniformaly spread into the top boundary layer or parametrisation layer)
+      !!
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      !!-------------------------- IN  -------------------------------------
+      INTEGER , DIMENSION(jpi,jpj), INTENT(in) :: ktop , kbot   ! top and bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phtbl, pfrac  ! thickness of the tbl and fraction of last cell affected by the tbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pvar2d        ! 2d var to map in 3d
+      CHARACTER(LEN=256), INTENT(in) :: cdvar
+      !!---------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk                       ! loop indices
+      INTEGER  :: ikt, ikb                         ! top and bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj)     :: zvar2d   !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvar3d   ! 3d var to output
+      !!---------------------------------------------------------------------
+      !
+      ! compute 3d output
+      zvar2d(:,:) = pvar2d(:,:) / phtbl(:,:)
+      zvar3d(:,:,:) = 0._wp
+      !
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = ktop(ji,jj)
+            ikb = kbot(ji,jj)
+            DO jk = ikt, ikb - 1
+               zvar3d(ji,jj,jk) = zvar2d(ji,jj) * e3t_n(ji,jj,jk)
+            END DO
+            zvar3d(ji,jj,ikb) = zvar2d(ji,jj) * e3t_n(ji,jj,ikb) * pfrac(ji,jj)
+         END DO
+      END DO
+      !
+      CALL iom_put( TRIM(cdvar) , zvar3d(:,:,:))
+      !
+   END SUBROUTINE isf_diags_2dto3d
+
+END MODULE isfdiags
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfhdiv.F90	(revision 12747)
@@ -0,0 +1,139 @@
+MODULE isfhdiv
+   !!======================================================================
+   !!                       ***  MODULE  isfhdiv  ***
+   !! ice shelf horizontal divergence module :  update the horizontal divergence
+   !!                   with the ice shelf melt and coupling correction
+   !!======================================================================
+   !! History :  4.0  !  2019-09  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isf_hdiv    : update the horizontal divergence with the ice shelf 
+   !!                 melt and coupling correction
+   !!----------------------------------------------------------------------
+
+   USE isf_oce                ! ice shelf
+
+   USE dom_oce                ! time and space domain
+   USE phycst , ONLY: r1_rau0 ! physical constant
+   USE in_out_manager         !
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_hdiv
+
+CONTAINS
+
+   SUBROUTINE isf_hdiv( kt, phdiv )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_hdiv  ***
+      !!       
+      !! ** Purpose :   update the horizontal divergence with the ice shelf contribution
+      !!                (parametrisation, explicit, ice sheet coupling conservation
+      !!                 increment)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdiv   ! horizontal divergence
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
+      !
+      IF ( ln_isf ) THEN
+         !
+         ! ice shelf cavity contribution
+         IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, phdiv)
+         !
+         ! ice shelf parametrisation contribution
+         IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, phdiv)
+         !
+         ! ice sheet coupling contribution 
+         IF ( ln_isfcpl .AND. kt /= 0 ) THEN
+            !
+            ! Dynamical stability at start up after change in under ice shelf cavity geometry is achieve by correcting the divergence.
+            ! This is achieved by applying a volume flux in order to keep the horizontal divergence after remapping 
+            ! the same as at the end of the latest time step. So correction need to be apply at nit000 (euler time step) and
+            ! half of it at nit000+1 (leap frog time step).
+            IF ( kt == nit000   ) CALL isf_hdiv_cpl(risfcpl_vol       , phdiv)
+            IF ( kt == nit000+1 ) CALL isf_hdiv_cpl(risfcpl_vol*0.5_wp, phdiv) ! risfcpl_vol*0.5_wp = 0.5 * (risfcpl_vol_b + risfcpl_vol_n)
+            !
+            ! correct divergence every time step to remove any trend due to coupling
+            ! conservation option
+            IF ( ln_isfcpl_cons ) CALL isf_hdiv_cpl(risfcpl_cons_vol, phdiv)
+            !
+         END IF
+         !
+      END IF
+      !
+   END SUBROUTINE isf_hdiv
+
+   SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE sbc_isf_div  ***
+      !!       
+      !! ** Purpose :   update the horizontal divergence with the ice shelf inflow
+      !!
+      !! ** Method  :   pfwf is positive (outflow) and expressed as kg/m2/s
+      !!                increase the divergence
+      !!
+      !! ** Action  :   phdivn   increased by the ice shelf outflow
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
+      !!----------------------------------------------------------------------
+      INTEGER , DIMENSION(jpi,jpj), INTENT(in   ) :: ktop , kbot
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfrac, phtbl
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pfwf , pfwf_b
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikt, ikb 
+      REAL(wp), DIMENSION(jpi,jpj) :: zhdiv
+      !!----------------------------------------------------------------------
+      !
+      !==   fwf distributed over several levels   ==!
+      !
+      ! compute integrated divergence correction
+      zhdiv(:,:) = 0.5_wp * ( pfwf(:,:) + pfwf_b(:,:) ) * r1_rau0 / phtbl(:,:)
+      !
+      ! update divergence at each level affected by ice shelf top boundary layer
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            ikt = ktop(ji,jj)
+            ikb = kbot(ji,jj)
+            ! level fully include in the ice shelf boundary layer
+            DO jk = ikt, ikb - 1
+               phdiv(ji,jj,jk) = phdiv(ji,jj,jk) + zhdiv(ji,jj)
+            END DO
+            ! level partially include in ice shelf boundary layer 
+            phdiv(ji,jj,ikb) = phdiv(ji,jj,ikb) + zhdiv(ji,jj) * pfrac(ji,jj)
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_hdiv_mlt
+
+   SUBROUTINE isf_hdiv_cpl(pqvol, phdiv)
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_hdiv_cpl  ***
+      !!       
+      !! ** Purpose :   update the horizontal divergence with the ice shelf 
+      !!                coupling conservation increment
+      !!
+      !! ** Method  :   pqvol is positive (outflow) and expressed as m3/s
+      !!                increase the divergence
+      !!
+      !! ** Action  :   phdivn   increased by the ice shelf outflow
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pqvol
+      !!----------------------------------------------------------------------
+      INTEGER :: jk
+      !!----------------------------------------------------------------------
+      !
+      DO jk=1,jpk 
+         phdiv(:,:,jk) =  phdiv(:,:,jk) + pqvol(:,:,jk) * r1_e1e2t(:,:) / e3t_n(:,:,jk)
+      END DO
+      !
+   END SUBROUTINE isf_hdiv_cpl
+
+END MODULE isfhdiv
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfload.F90	(revision 12747)
@@ -0,0 +1,115 @@
+MODULE isfload
+   !!======================================================================
+   !!                       ***  MODULE  isfload  ***
+   !! isfload module :  compute ice shelf load (needed for the hpg)
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfload      : compute ice shelf load
+   !!----------------------------------------------------------------------
+
+   USE isf_oce, ONLY: cn_isfload, rn_isfload_T, rn_isfload_S ! ice shelf variables
+
+   USE dom_oce, ONLY: e3w_n, gdept_n, risfdep, mikt          ! vertical scale factor
+   USE eosbn2 , ONLY: eos                                    ! eos routine
+
+   USE lib_mpp, ONLY: ctl_stop                               ! ctl_stop routine
+   USE in_out_manager                                        ! 
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_load
+
+CONTAINS
+
+   SUBROUTINE isf_load ( pisfload )
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_load  ***
+      !!
+      !! ** Purpose : compute the ice shelf load
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
+      !!-------------------------- IN  -------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      ! quality test: ice shelf in a stratify/uniform ocean should not drive any flow.
+      !               the smaller the residual flow is, the better it is.
+      !
+      ! ice shelf cavity
+      SELECT CASE ( cn_isfload )
+      CASE ( 'uniform' )
+         CALL isf_load_uniform ( pisfload )
+      CASE DEFAULT
+         CALL ctl_stop('STOP','method cn_isfload to compute ice shelf load does not exist (isomip), check your namelist')
+      END SELECT
+      !
+   END SUBROUTINE isf_load
+
+   SUBROUTINE isf_load_uniform( pisfload )
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_load  ***
+      !!
+      !! ** Purpose : compute the ice shelf load
+      !!
+      !! ** Method  : The ice shelf is assumed to be in hydro static equilibrium
+      !!              in water at -1.9 C and 34.4 PSU. Weight of the ice shelf is
+      !!              integrated from top to bottom.
+      !!
+      !!--------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pisfload
+      !!-------------------------- IN  -------------------------------------
+      !!--------------------------------------------------------------------
+      INTEGER  :: ji, jj, jk
+      INTEGER  :: ikt
+      REAL(wp)                          :: znad        ! 
+      REAL(wp), DIMENSION(jpi,jpj)      :: zrhdtop_isf ! water density    displaced by the ice shelf (at the interface)
+      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts_top     ! water properties displaced by the ice shelf   
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: zrhd        ! water density    displaced by the ice shelf
+      !!----------------------------------------------------------------------
+      !
+      znad = 1._wp                     !- To use density and not density anomaly
+      !
+      !                                !- assume water displaced by the ice shelf is at T=rn_isfload_T and S=rn_isfload_S (rude)
+      zts_top(:,:,jp_tem) = rn_isfload_T   ;   zts_top(:,:,jp_sal) = rn_isfload_S
+      !
+      DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
+         CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )
+      END DO
+      !
+      !                                !- compute rhd at the ice/oce interface (ice shelf side)
+      CALL eos( zts_top , risfdep, zrhdtop_isf )
+      !
+      !                                !- Surface value + ice shelf gradient
+      pisfload(:,:) = 0._wp                       ! compute pressure due to ice shelf load 
+      DO jj = 1, jpj                         ! (used to compute hpgi/j for all the level from 1 to miku/v)
+         DO ji = 1, jpi                      ! divided by 2 later
+            ikt = mikt(ji,jj)
+            !
+            IF ( ikt > 1 ) THEN
+               !
+               ! top layer of the ice shelf
+               pisfload(ji,jj) = pisfload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1)
+               !
+               ! core layers of the ice shelf
+               DO jk = 2, ikt-1
+                  pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk)
+               END DO
+               !
+               ! deepest part of the ice shelf (between deepest T point and ice/ocean interface
+               pisfload(ji,jj) = pisfload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
+                  &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )
+               !
+            END IF
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_load_uniform
+
+END MODULE isfload
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90	(revision 12747)
@@ -0,0 +1,185 @@
+MODULE isfpar
+   !!======================================================================
+   !!                       ***  MODULE  isfpar  ***
+   !! ice shelf module :  update ocean boundary condition under ice
+   !!                   shelf
+   !!======================================================================
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!            4.1  !  2019-09  (P. Mathiot) Restructuration
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfpar       : compute ice shelf melt using a prametrisation of ice shelf cavities
+   !!----------------------------------------------------------------------
+   USE isf_oce        ! ice shelf
+   !
+   USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
+   USE isftbl   , ONLY: isf_tbl_ktop, isf_tbl_lvl ! ice shelf top boundary layer properties subroutine
+   USE isfparmlt, ONLY: isfpar_mlt                ! ice shelf melt formulation subroutine
+   USE isfdiags , ONLY: isf_diags_flx             ! ice shelf diags subroutine
+   USE isfutils , ONLY: debug, read_2dcstdta      ! ice shelf debug subroutine
+   !
+   USE dom_oce  , ONLY: bathy          ! ocean space and time domain
+   USE par_oce  , ONLY: jpi,jpj        ! ocean space and time domain
+   USE phycst   , ONLY: r1_rau0_rcp    ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE fldread        ! read input field at current time step
+   USE lbclnk         ! lbc_lnk
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   isf_par, isf_par_init
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+ 
+   SUBROUTINE isf_par( kt, ptsc, pqfwf )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE isf_par ***      
+      !!
+      !! ** Purpose : compute the heat and fresh water due to ice shelf melting/freezing using a parametrisation 
+      !!
+      !! ** Comment : in isf_par and all its call tree, 
+      !!              'tbl' means parametrisation layer (ie how the far field temperature/salinity is computed) 
+      !!              instead of in a proper top boundary layer as at the ice shelf ocean interface
+      !!              as the action to compute the properties of the tbl or the parametrisation layer are the same,
+      !!              (ie average T/S over a specific depth (can be across multiple levels))
+      !!              the name tbl was kept.
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT --------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pqfwf
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc
+      !!-------------------------- IN  --------------------------------------
+      INTEGER, INTENT(in) ::   kt                                           ! ocean time step
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
+      !!---------------------------------------------------------------------
+      !
+      ! compute heat content, latent heat and melt fluxes (2d)
+      CALL isfpar_mlt( kt, zqhc, zqoce, pqfwf  )
+      !
+      ! compute heat and water flux ( > 0 out )
+      pqfwf(:,:) = pqfwf(:,:) * mskisf_par(:,:)
+      zqoce(:,:) = zqoce(:,:) * mskisf_par(:,:)
+      zqhc (:,:) = zqhc(:,:)  * mskisf_par(:,:)
+      !
+      ! compute heat content flux ( > 0 out )
+      zqlat(:,:) = pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
+      !
+      ! total heat flux ( > 0 out )
+      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
+      !
+      ! lbclnk on melt and heat fluxes
+      CALL lbc_lnk_multi( 'isfmlt', zqh, 'T', 1., pqfwf, 'T', 1.)
+      !
+      ! output fluxes
+      CALL isf_diags_flx( misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc)
+      !
+      ! set temperature content
+      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rau0_rcp
+      !
+      ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+      IF (lrst_oce) CALL isfrst_write(kt, 'par', ptsc, pqfwf)
+      !
+      IF ( ln_isfdebug ) THEN
+         CALL debug('isf_par: ptsc T',ptsc(:,:,1))
+         CALL debug('isf_par: ptsc S',ptsc(:,:,2))
+         CALL debug('isf_par: pqfwf fwf',pqfwf(:,:))
+      END IF
+      !
+   END SUBROUTINE isf_par
+
+   SUBROUTINE isf_par_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_par_init  ***
+      !!
+      !! ** Purpose : initialisation of the variable needed for the parametrisation of ice shelf melt
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER               :: ierr
+      REAL(wp), DIMENSION(jpi,jpj) :: ztblmax, ztblmin
+      !!----------------------------------------------------------------------
+      !
+      ! allocation
+      CALL isf_alloc_par()
+      !
+      ! initialisation
+      misfkt_par(:,:)     = 1         ; misfkb_par(:,:)       = 1         
+      rhisf_tbl_par(:,:)  = 1e-20     ; rfrac_tbl_par(:,:)    = 0.0_wp
+      !
+      ! define isf tbl tickness, top and bottom indice
+      CALL read_2dcstdta(TRIM(sn_isfpar_zmax%clname), TRIM(sn_isfpar_zmax%clvar), ztblmax)
+      CALL read_2dcstdta(TRIM(sn_isfpar_zmin%clname), TRIM(sn_isfpar_zmin%clvar), ztblmin)
+      !
+      ! mask ice shelf parametrisation location
+      ztblmax(:,:) = ztblmax(:,:) * ssmask(:,:)
+      ztblmin(:,:) = ztblmin(:,:) * ssmask(:,:)
+      !
+      ! if param used under an ice shelf overwrite ztblmin by the ice shelf draft
+      WHERE ( risfdep > 0._wp .AND. ztblmin > 0._wp )
+         ztblmin(:,:) = risfdep(:,:)
+      END WHERE
+      !
+      ! ensure ztblmax <= bathy
+      WHERE ( ztblmax(:,:) > bathy(:,:) )
+         ztblmax(:,:) = bathy(:,:)
+      END WHERE
+      !
+      ! compute ktop and update ztblmin to gdepw_0(misfkt_par) 
+      CALL isf_tbl_ktop(ztblmin, misfkt_par) !   out: misfkt_par
+      !                                      ! inout: ztblmin
+      !
+      ! initial tbl thickness
+      rhisf0_tbl_par(:,:) = ztblmax(:,:) - ztblmin(:,:)
+      !
+      ! define iceshelf parametrisation mask
+      mskisf_par = 0
+      WHERE ( rhisf0_tbl_par(:,:) > 0._wp )
+         mskisf_par(:,:) = 1._wp
+      END WHERE
+      !
+      ! read par variable from restart
+      IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
+      !
+      SELECT CASE ( TRIM(cn_isfpar_mlt) )
+         !
+      CASE ( 'spe' )
+         !
+         ALLOCATE( sf_isfpar_fwf(1), STAT=ierr )
+         ALLOCATE( sf_isfpar_fwf(1)%fnow(jpi,jpj,1), sf_isfpar_fwf(1)%fdta(jpi,jpj,1,2) )
+         CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_isfdir, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '      ==>>>   ice melt read from forcing field (cn_isfmlt_par = spe)'
+         !
+      CASE ( 'bg03' )
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (cn_isfmlt_par = bg03)'
+         !
+         ! read effective length
+         CALL read_2dcstdta(TRIM(sn_isfpar_Leff%clname), TRIM(sn_isfpar_Leff%clvar), risfLeff)
+         risfLeff = risfLeff*1000.0_wp           !: convertion in m
+         !
+      CASE ( 'oasis' )
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '      ==>>>    isf melt provided by OASIS (cn_isfmlt_par = oasis)'
+         !
+      CASE DEFAULT
+         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
+      END SELECT
+      !
+   END SUBROUTINE isf_par_init
+
+END MODULE isfpar
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfparmlt.F90	(revision 12747)
@@ -0,0 +1,220 @@
+MODULE isfparmlt
+   !!======================================================================
+   !!                       ***  MODULE  isfparmlt  ***
+   !! Ice shelf parametrisation module :  update surface ocean boundary condition under ice
+   !!                   shelf using an ice shelf melt parametrisation
+   !!======================================================================
+   !! History :  4.0  ! original code
+   !!----------------------------------------------------------------------
+
+   USE isf_oce                  ! ice shelf
+   USE isftbl , ONLY: isf_tbl   ! ice shelf depth average
+
+   USE dom_oce                  ! ocean space and time domain
+   USE oce    , ONLY: tsn       ! ocean dynamics and tracers
+   USE phycst , ONLY: rcp, rau0 ! physical constants
+   USE eosbn2 , ONLY: eos_fzp   ! equation of state
+
+   USE in_out_manager              ! I/O manager
+   USE iom        , ONLY: iom_put  ! I/O library
+   USE fldread    , ONLY: fld_read, FLD, FLD_N !
+   USE lib_fortran, ONLY: glob_sum !
+   USE lib_mpp    , ONLY: ctl_stop !
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC  isfpar_mlt 
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PUBLIC SUBROUTINE ----------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+  SUBROUTINE isfpar_mlt( kt, pqhc, pqoce, pqfwf )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  2 parameterizations are available according
+      !!                        1 : Specified melt flux
+      !!                        2 : Beckmann & Goose parameterization
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqfwf, pqoce, pqhc  ! fresh water, ice-ocean heat and heat content fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      ! Choose among the available ice shelf parametrisation
+      SELECT CASE ( cn_isfpar_mlt )
+      CASE ( 'spe' )    ! specified runoff in depth (Mathiot et al., 2017 in preparation)
+         CALL isfpar_mlt_spe(kt, pqhc, pqoce, pqfwf)
+      CASE ( 'bg03' )    ! Beckmann and Goosse parametrisation 
+         CALL isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
+      CASE ( 'oasis' )
+         CALL isfpar_mlt_oasis( kt, pqhc, pqoce, pqfwf)
+      CASE DEFAULT
+         CALL ctl_stop('STOP', 'unknown isf melt formulation : cn_isfpar (should not see this)')
+      END SELECT
+      !
+   END SUBROUTINE isfpar_mlt
+
+! -------------------------------------------------------------------------------------------------------
+! -------------------------------- PRIVATE SUBROUTINE ---------------------------------------------------
+! -------------------------------------------------------------------------------------------------------
+
+   SUBROUTINE isfpar_mlt_spe(kt, pqhc, pqoce, pqfwf)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt_spe  ***
+      !!
+      !! ** Purpose : prescribed ice shelf melting in case ice shelf cavities are closed.
+      !!              data read into a forcing files.
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf, pqoce  ! fresh water and ice-ocean heat fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  INTENT(in) :: kt
+      !!--------------------------------------------------------------------
+      INTEGER :: jk
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
+      REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz
+      !!--------------------------------------------------------------------
+      !
+      ! 0. ------------Read specified runoff
+      CALL fld_read ( kt, 1, sf_isfpar_fwf   )
+      !
+      ! compute ptfrz
+      ! 1. ------------Mean freezing point
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+      !
+      pqfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)      ! fresh water flux from the isf (fwfisf <0 mean melting) 
+      pqoce(:,:) =   pqfwf(:,:) * rLfusisf             ! ocean/ice shelf flux assume to be equal to latent heat flux
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
+      !
+      CALL iom_put('isftfrz_par', ztfrz )
+      !
+   END SUBROUTINE isfpar_mlt_spe
+
+   SUBROUTINE isfpar_mlt_bg03(kt, pqhc, pqoce, pqfwf)
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt_bg03  ***
+      !!
+      !! ** Purpose : compute an estimate of ice shelf melting and 
+      !!              latent, ocean-ice and heat content heat fluxes
+      !!              in case cavities are closed based on the far fields T and S properties. 
+      !!
+      !! ** Method  : The ice shelf melt is computed as proportional to the differences between the 
+      !!              mean temperature and mean freezing point in front of the ice shelf averaged 
+      !!              over the ice shelf min ice shelf draft and max ice shelf draft and the freezing point
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!                interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!----------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pqhc, pqfwf, pqoce  ! fresh water and ice-ocean heat fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  INTENT(in) :: kt
+      !!--------------------------------------------------------------------
+      INTEGER :: jk
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztfrz3d        ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztfrz          ! freezing point
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztavg          ! temperature avg
+      !!----------------------------------------------------------------------
+      !
+      ! 0. ------------Mean freezing point
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+      !
+      ! 1. ------------Mean temperature
+      CALL isf_tbl(tsn(:,:,jk,jp_tem), ztavg, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+      !
+      ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
+      pqoce(:,:) =   rau0 * rcp * rn_gammat0 * risfLeff(:,:) * e1t(:,:) * ( ztavg(:,:) - ztfrz(:,:) ) * r1_e1e2t(:,:)
+      pqfwf(:,:) = - pqoce(:,:) / rLfusisf             ! derived from the latent heat flux
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux 
+      !
+      ! 3. ------------BG03 output
+      ! output ttbl
+      CALL iom_put('ttbl_par', ztavg(:,:) * mskisf_par(:,:) )
+      !
+      ! output thermal driving
+      CALL iom_put('isfthermald_par',( ztfrz(:,:) - ztavg(:,:) ) * mskisf_par(:,:))
+      !
+      ! output freezing point used to define the thermal driving and heat content fluxes
+      CALL iom_put('isftfrz_par', ztfrz )
+      !
+   END SUBROUTINE isfpar_mlt_bg03
+
+   SUBROUTINE isfpar_mlt_oasis(kt, pqhc , pqoce, pqfwf )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE isfpar_mlt_oasis  ***
+      !!
+      !! ** Purpose    : scale the fwf read from input file by the total amount received by the sbccpl interface
+      !!
+      !! ** Purpose    : - read ice shelf melt from forcing file and scale it by the input file total amount => pattern
+      !!                 - compute total amount of fwf given by sbccpl (fwfisf_oasis)
+      !!                 - scale fwf and compute heat fluxes
+      !!
+      !!---------------------------------------------------------------------
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pqhc, pqoce, pqfwf  ! heat content, latent heat and fwf fluxes
+      !!-------------------------- IN  -------------------------------------
+      INTEGER                     , INTENT(in   ) :: kt                  ! current time step
+      !!--------------------------------------------------------------------
+      INTEGER                           :: jk                            ! loop index
+      REAL(wp)                          :: zfwf_fld, zfwf_oasis          ! total fwf in the forcing fields (pattern) and from the cpl interface (amount)
+      REAL(wp), DIMENSION(jpi,jpj)      :: ztfrz                         ! tbl freezing temperature
+      REAL(wp), DIMENSION(jpi,jpj)      :: zfwf                          ! 2d fwf map after scaling
+      REAL(wp), DIMENSION(jpi,jpj,jpk)  :: ztfrz3d
+      !!--------------------------------------------------------------------
+      !
+      ! 0. ------------Read specified runoff
+      CALL fld_read ( kt, 1, sf_isfpar_fwf   )
+      !
+      ! 1. ------------Mean freezing point (needed for heat content flux)
+      DO jk = 1,jpk
+         CALL eos_fzp(tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), gdept_n(:,:,jk))
+      END DO
+      CALL isf_tbl(ztfrz3d, ztfrz, 'T', misfkt_par, rhisf_tbl_par, misfkb_par, rfrac_tbl_par )
+      !
+      ! 2. ------------Scale isf melt pattern with total amount from oasis
+      ! ice shelf 2d map
+      zfwf(:,:) = - sf_isfpar_fwf(1)%fnow(:,:,1)
+      !
+      ! compute glob sum from input file
+      ! (PM) should we consider delay sum as in fwb ? (it will offset by 1 time step if I understood well)
+      zfwf_fld = glob_sum('isfcav_mlt', e1e2t(:,:) * zfwf(:,:))
+      !
+      ! compute glob sum from atm->oce ice shelf fwf
+      ! (PM) should we consider delay sum as in fwb ?
+      zfwf_oasis = glob_sum('isfcav_mlt', e1e2t(:,:) * fwfisf_oasis(:,:))
+      !
+      ! scale fwf
+      zfwf(:,:) = zfwf(:,:) * zfwf_oasis / zfwf_fld
+      ! 
+      ! 3. -----------Define fwf and qoce
+      ! ocean heat flux is assume to be equal to the latent heat
+      pqfwf(:,:) =   zfwf(:,:)                         ! fwf                ( >0 out )
+      pqoce(:,:) = - pqfwf(:,:) * rLfusisf             ! ocean heat flux    ( >0 out ) (assumed to be the latent heat flux)
+      pqhc (:,:) =   pqfwf(:,:) * ztfrz(:,:) * rcp     ! heat content flux  ( >0 out )
+      !
+      CALL iom_put('isftfrz_par', ztfrz )
+      !
+   END SUBROUTINE isfpar_mlt_oasis
+
+END MODULE isfparmlt
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfrst.F90	(revision 12747)
@@ -0,0 +1,105 @@
+MODULE isfrst
+   !!======================================================================
+   !!                       ***  MODULE  isfrst  ***
+   !! iceshelf restart module :read/write iceshelf variables from/in restart
+   !!======================================================================
+   !! History :  4.1  !  2019-07  (P. Mathiot) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfrst : read/write iceshelf variables in/from restart
+   !!----------------------------------------------------------------------
+   !
+   USE par_oce, ONLY: jpi,jpj,jpk,jpts ! time and space domain
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   !
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isfrst_read, isfrst_write ! iceshelf restart read and write 
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+   ! 
+   SUBROUTINE isfrst_read(cdisf, ptsc, pfwf, ptsc_b, pfwf_b )
+      !!---------------------------------------------------------------------
+      !!
+      !!   isfrst_read : read iceshelf variables from restart
+      !!
+      !!-------------------------- OUT --------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(  out) :: pfwf_b
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(  out) :: ptsc_b
+      !!-------------------------- IN  --------------------------------------
+      CHARACTER(LEN=3)                 , INTENT(in   ) :: cdisf
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
+      !!----------------------------------------------------------------------
+      !
+      ! define variable name
+      cfwf_b = 'fwfisf_'//TRIM(cdisf)//'_b'
+      chc_b  = 'isf_hc_'//TRIM(cdisf)//'_b'
+      csc_b  = 'isf_sc_'//TRIM(cdisf)//'_b'
+      !
+      ! read restart
+      IF( iom_varid( numror, cfwf_b, ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
+         CALL iom_get( numror, jpdom_autoglo, cfwf_b, pfwf_b(:,:)        , ldxios = lrxios )   ! before ice shelf melt
+         CALL iom_get( numror, jpdom_autoglo, chc_b , ptsc_b (:,:,jp_tem), ldxios = lrxios )   ! before ice shelf heat flux
+         CALL iom_get( numror, jpdom_autoglo, csc_b , ptsc_b (:,:,jp_sal), ldxios = lrxios )   ! before ice shelf heat flux
+      ELSE
+         pfwf_b(:,:)   = pfwf(:,:)
+         ptsc_b(:,:,:) = ptsc(:,:,:)
+      ENDIF
+      !
+      IF( lwxios ) THEN
+         CALL iom_set_rstw_var_active(TRIM(chc_b ))
+         CALL iom_set_rstw_var_active(TRIM(csc_b ))
+         CALL iom_set_rstw_var_active(TRIM(cfwf_b))
+      ENDIF
+
+   END SUBROUTINE isfrst_read
+   ! 
+   SUBROUTINE isfrst_write(kt, cdisf, ptsc, pfwf )
+      !!---------------------------------------------------------------------
+      !!
+      !!   isfrst_write : write iceshelf variables in restart
+      !!
+      !!-------------------------- IN  --------------------------------------
+      INTEGER                          , INTENT(in   ) :: kt
+      CHARACTER(LEN=3)                 , INTENT(in   ) :: cdisf
+      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(in   ) :: pfwf
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) :: ptsc
+      !!---------------------------------------------------------------------
+      CHARACTER(LEN=256) :: cfwf_b, chc_b, csc_b
+      !!---------------------------------------------------------------------
+      !
+      ! ocean output print
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'isf : isf fwf and heat fluxes written in ocean restart file ',   &
+         &                    'at it= ', kt,' date= ', ndastp
+      IF(lwp) WRITE(numout,*) '~~~~'
+      !
+      ! define variable name
+      cfwf_b = 'fwfisf_'//TRIM(cdisf)//'_b'
+      chc_b  = 'isf_hc_'//TRIM(cdisf)//'_b'
+      csc_b  = 'isf_sc_'//TRIM(cdisf)//'_b'
+      !
+      ! write restart variable
+      IF( lwxios ) CALL iom_swap( cwxios_context )
+      CALL iom_rstput( kt, nitrst, numrow, cfwf_b, pfwf(:,:)       , ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, chc_b , ptsc(:,:,jp_tem), ldxios = lwxios )
+      CALL iom_rstput( kt, nitrst, numrow, csc_b , ptsc(:,:,jp_sal), ldxios = lwxios )
+      IF( lwxios ) CALL iom_swap( cxios_context )
+      !
+   END SUBROUTINE isfrst_write
+   !
+END MODULE isfrst
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfstp.F90	(revision 12747)
@@ -0,0 +1,303 @@
+MODULE isfstp
+   !!======================================================================
+   !!                       ***  MODULE  isfstp  ***
+   !! Surface module :  compute iceshelf load, melt and heat flux
+   !!======================================================================
+   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
+   !!            4.1  !  2019-09  (P. Mathiot) Split param/explicit ice shelf and re-organisation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfstp       : compute iceshelf melt and heat flux
+   !!----------------------------------------------------------------------
+   !
+   USE isf_oce                                      ! isf variables
+   USE isfload, ONLY: isf_load                      ! ice shelf load
+   USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer
+   USE isfpar , ONLY: isf_par, isf_par_init         ! ice shelf parametrisation
+   USE isfcav , ONLY: isf_cav, isf_cav_init         ! ice shelf cavity
+   USE isfcpl , ONLY: isfcpl_rst_write, isfcpl_init ! isf variables
+
+   USE dom_oce, ONLY: ht_n, e3t_n, ln_isfcav, ln_linssh ! ocean space and time domain
+   USE domvvl,  ONLY: ln_vvl_zstar                      ! zstar logical
+   USE zdfdrg,  ONLY: r_Cdmin_top, r_ke0_top            ! vertical physics: top/bottom drag coef.
+   !
+   USE lib_mpp, ONLY: ctl_stop, ctl_nam
+   USE fldread, ONLY: FLD, FLD_N
+   USE in_out_manager ! I/O manager
+   USE timing
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC   isf_stp, isf_init, isf_nam  ! routine called in sbcmod and divhor
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+ 
+  SUBROUTINE isf_stp( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_stp  ***
+      !!
+      !! ** Purpose : compute total heat flux and total fwf due to ice shelf melt
+      !!
+      !! ** Method  : For each case (parametrisation or explicity cavity) :
+      !!              - define the before fields
+      !!              - compute top boundary layer properties
+      !!                (in case of parametrisation, this is the 
+      !!                 depth range model array used to compute mean far fields properties)
+      !!              - compute fluxes
+      !!              - write restart variables
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('isf')
+      !
+      !=======================================================================
+      ! 1.: compute melt and associated heat fluxes in the ice shelf cavities
+      !=======================================================================
+      !
+      IF ( ln_isfcav_mlt ) THEN
+         !
+         ! 1.1: before time step 
+         IF ( kt /= nit000 ) THEN 
+            risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
+            fwfisf_cav_b(:,:)      = fwfisf_cav(:,:)
+         END IF
+         !
+         ! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
+         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         !
+         ! 1.3: compute ice shelf melt
+         CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav)
+         !
+      END IF
+      ! 
+      !=================================================================================
+      ! 2.: compute melt and associated heat fluxes for not resolved ice shelf cavities
+      !=================================================================================
+      !
+      IF ( ln_isfpar_mlt ) THEN
+         !
+         ! 2.1: before time step 
+         IF ( kt /= nit000 ) THEN 
+            risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
+            fwfisf_par_b  (:,:)   = fwfisf_par  (:,:)
+         END IF
+         !
+         ! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         ! by simplicity, we assume the top level where param applied do not change with time (done in init part)
+         rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
+         CALL isf_tbl_lvl(ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         !
+         ! 2.3: compute ice shelf melt
+         CALL isf_par( kt, risf_par_tsc, fwfisf_par)
+         !
+      END IF
+      !
+      !==================================================================================
+      ! 3.: output specific restart variable in case of coupling with an ice sheet model
+      !==================================================================================
+      !
+      IF ( ln_isfcpl .AND. lrst_oce ) CALL isfcpl_rst_write(kt)
+      !
+      IF( ln_timing )   CALL timing_stop('isf')
+      !
+   END SUBROUTINE isf_stp
+
+   SUBROUTINE isf_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isfstp_init  ***
+      !!
+      !! ** Purpose :   Initialisation of the ice shelf public variables
+      !!
+      !! ** Method  :   Read the namisf namelist, check option compatibility and set derived parameters
+      !!
+      !! ** Action  : - read namisf parameters
+      !!              - allocate memory
+      !!              - output print
+      !!              - ckeck option compatibility
+      !!              - call cav/param/isfcpl init routine
+      !!----------------------------------------------------------------------
+      !
+      ! constrain: l_isfoasis need to be known
+      !
+      ! Read namelist
+      CALL isf_nam()
+      !
+      ! Allocate public array
+      CALL isf_alloc()
+      !
+      ! check option compatibility
+      CALL isf_ctl()
+      !
+      ! compute ice shelf load
+      IF ( ln_isfcav ) CALL isf_load( risfload )
+      !
+      ! terminate routine now if no ice shelf melt formulation specify
+      IF ( ln_isf ) THEN
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation melt in the cavity
+         IF ( ln_isfcav_mlt ) CALL isf_cav_init()
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation parametrised melt
+         IF ( ln_isfpar_mlt ) CALL isf_par_init()
+         !
+         !---------------------------------------------------------------------------------------------------------------------
+         ! initialisation ice sheet coupling
+         IF ( ln_isfcpl ) CALL isfcpl_init()
+         !
+      END IF
+         
+  END SUBROUTINE isf_init
+
+  SUBROUTINE isf_ctl()
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_ctl  ***
+      !!
+      !! ** Purpose :   output print and option compatibility check
+      !!
+      !!----------------------------------------------------------------------
+      IF (lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'isf_init : ice shelf initialisation'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namisf :'
+         !
+         WRITE(numout,*) '   ice shelf cavity (open or parametrised)  ln_isf = ', ln_isf
+         WRITE(numout,*)
+         !
+         IF ( ln_isf ) THEN
+            WRITE(numout,*) '      Add debug print in isf module           ln_isfdebug     = ', ln_isfdebug
+            WRITE(numout,*)
+            WRITE(numout,*) '      melt inside the cavity                  ln_isfcav_mlt   = ', ln_isfcav_mlt
+            IF ( ln_isfcav_mlt) THEN
+               WRITE(numout,*) '         melt formulation                         cn_isfcav_mlt= ', TRIM(cn_isfcav_mlt)
+               WRITE(numout,*) '         thickness of the top boundary layer      rn_htbl      = ', rn_htbl
+               WRITE(numout,*) '         gamma formulation                        cn_gammablk  = ', TRIM(cn_gammablk) 
+               IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN 
+                  WRITE(numout,*) '         gammat coefficient                       rn_gammat0   = ', rn_gammat0  
+                  WRITE(numout,*) '         gammas coefficient                       rn_gammas0   = ', rn_gammas0  
+                  WRITE(numout,*) '         top background ke used (from namdrg_top) rn_ke0       = ', r_ke0_top
+                  WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0       = ', r_Cdmin_top
+               END IF
+            END IF
+            WRITE(numout,*) ''
+            !
+            WRITE(numout,*) '      ice shelf melt parametrisation          ln_isfpar_mlt   = ', ln_isfpar_mlt
+            IF ( ln_isfpar_mlt ) THEN
+               WRITE(numout,*) '         isf parametrisation formulation         cn_isfpar_mlt   = ', TRIM(cn_isfpar_mlt)
+            END IF
+            WRITE(numout,*) ''
+            !
+            WRITE(numout,*) '      Coupling to an ice sheet model          ln_isfcpl       = ', ln_isfcpl
+            IF ( ln_isfcpl ) THEN
+               WRITE(numout,*) '         conservation activated ln_isfcpl_cons     = ', ln_isfcpl_cons
+               WRITE(numout,*) '         number of call of the extrapolation loop  = ', nn_drown
+            ENDIF
+            WRITE(numout,*) ''
+            !
+         ELSE
+            !
+            IF ( ln_isfcav ) THEN
+               WRITE(numout,*) ''
+               WRITE(numout,*) '   W A R N I N G: ice shelf cavities are open BUT no melt will be computed or read from file !'
+               WRITE(numout,*) ''
+            END IF
+            !
+         END IF
+
+         IF (ln_isfcav) THEN
+            WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload)
+            WRITE(numout,*) '         Temperature used to compute the ice shelf load            = ', rn_isfload_T
+            WRITE(numout,*) '         Salinity    used to compute the ice shelf load            = ', rn_isfload_S
+         END IF
+         WRITE(numout,*) ''
+         FLUSH(numout)
+
+      END IF
+      !
+
+      !---------------------------------------------------------------------------------------------------------------------
+      ! sanity check  ! issue ln_isfcav not yet known as well as l_isfoasis  => move this call in isf_stp ?
+      ! melt in the cavity without cavity
+      IF ( ln_isfcav_mlt .AND. (.NOT. ln_isfcav) ) &
+          &   CALL ctl_stop('ice shelf melt in the cavity activated (ln_isfcav_mlt) but no cavity detected in domcfg (ln_isfcav), STOP' )
+      !
+      ! ice sheet coupling without cavity
+      IF ( ln_isfcpl .AND. (.NOT. ln_isfcav) ) &
+         &   CALL ctl_stop('coupling with an ice sheet model detected (ln_isfcpl) but no cavity detected in domcfg (ln_isfcav), STOP' )
+      !
+      IF ( ln_isfcpl .AND. ln_isfcpl_cons .AND. ln_linssh ) &
+         &   CALL ctl_stop( 'The coupling between NEMO and an ice sheet model with the conservation option does not work with the linssh option' )
+      !
+      IF ( l_isfoasis .AND. .NOT. ln_isf ) CALL ctl_stop( ' OASIS send ice shelf fluxes to NEMO but NEMO does not have the isf module activated' )
+      !
+      IF ( l_isfoasis .AND. ln_isf ) THEN
+         !
+         CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
+         !
+         ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
+         IF ( ln_isfcav_mlt .AND. TRIM(cn_isfcav_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis is the only option availble if fwf send by oasis' )
+         IF ( ln_isfpar_mlt .AND. TRIM(cn_isfpar_mlt) /= 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis is the only option availble if fwf send by oasis' )
+         !
+         ! oasis melt computation not tested (coded but not tested)
+         IF ( ln_isfcav_mlt .OR. ln_isfpar_mlt ) THEN
+            IF ( TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfcav_mlt = oasis not tested' )
+            IF ( TRIM(cn_isfpar_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis not tested' )
+         END IF
+         !
+         ! oasis melt computation with cavity open and cavity parametrised (not coded)
+         IF ( ln_isfcav_mlt .AND. ln_isfpar_mlt ) THEN
+            IF ( TRIM(cn_isfpar_mlt) == 'oasis' .AND. TRIM(cn_isfcav_mlt) == 'oasis' ) CALL ctl_stop( 'cn_isfpar_mlt = oasis and cn_isfcav_mlt = oasis not coded' )
+         END IF
+         !
+         ! compatibility ice shelf and vvl
+         IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
+         !
+      END IF
+   END SUBROUTINE isf_ctl
+   !
+   SUBROUTINE isf_nam
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_nam  ***
+      !!
+      !! ** Purpose :   Read ice shelf namelist cfg and ref
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER               :: ios                  ! Local integer output status for namelist read
+      !!----------------------------------------------------------------------
+      NAMELIST/namisf/ ln_isf        ,                                                           & 
+         &             cn_gammablk   , rn_gammat0    , rn_gammas0    , rn_htbl, sn_isfcav_fwf,   &
+         &             ln_isfcav_mlt , cn_isfcav_mlt , sn_isfcav_fwf ,                           &
+         &             ln_isfpar_mlt , cn_isfpar_mlt , sn_isfpar_fwf ,                           &
+         &             sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff,                           &
+         &             ln_isfcpl     , nn_drown      , ln_isfcpl_cons, ln_isfdebug,              &
+         &             cn_isfload    , rn_isfload_T  , rn_isfload_S  , cn_isfdir
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+      READ  ( numnam_ref, namisf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namisf in reference namelist' )
+      !
+      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+      READ  ( numnam_cfg, namisf, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' )
+      IF(lwm) WRITE ( numond, namisf )
+
+   END SUBROUTINE isf_nam
+   !!
+   !!======================================================================
+END MODULE isfstp
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isftbl.F90	(revision 12747)
@@ -0,0 +1,274 @@
+MODULE isftbl
+   !!======================================================================
+   !!                       ***  MODULE  isftbl  ***
+   !! isftbl module :  compute properties of top boundary layer
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isftbl       : routine to compute : 
+   !!                  - geometry of the ice shelf tbl (isf_tbl_lvl, isftbl_ktop, isftbl_kbot)
+   !!                    (top and bottom level, thickness and fraction of deepest level affected)
+   !!                  - tbl averaged properties (isf_tbl, isf_tbl_avg)
+   !!----------------------------------------------------------------------
+
+   USE isf_oce ! ice shelf variables
+
+   USE dom_oce ! vertical scale factor and depth
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
+
+CONTAINS
+
+   SUBROUTINE isf_tbl( pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
+      !!--------------------------------------------------------------------
+      !!                  ***  SUBROUTINE isf_tbl  ***
+      !!
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
+      !!
+      !! ** Method : Average properties over a specific thickness
+      !!
+      !! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
+      !!                https://doi.org/10.1029/2007JC004368 , 2008
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(len=1)                      , INTENT(in   ) :: cd_ptin       ! point of variable in/out
+      REAL(wp), DIMENSION(jpi,jpj,jpk)      , INTENT(in   ) :: pvarin        ! 3d variable to average over the tbl
+      INTEGER,  DIMENSION(jpi,jpj)          , INTENT(in   ) :: ktop          ! top level
+      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in   ) :: phtbl         ! tbl thickness
+      !!-------------------------- IN OPTIONAL -----------------------------
+      INTEGER,  DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: kbot          ! bottom level
+      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: pfrac         ! fraction of bottom cell affected by tbl
+      !!--------------------------------------------------------------------
+      INTEGER ::   ji, jj                     ! loop index
+      INTEGER , DIMENSION(jpi,jpj) :: ikbot   ! bottom level of the tbl
+      REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
+      REAL(wp), DIMENSION(jpi,jpj) :: zhtbl   ! thickness of the tbl
+      REAL(wp), DIMENSION(jpi,jpj) :: zfrac   ! thickness of the tbl
+      !!--------------------------------------------------------------------
+      ! 
+      SELECT CASE ( cd_ptin )
+      CASE ( 'U' )
+         !
+         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
+         zhtbl = phtbl
+         !
+         ! compute tbl lvl and thickness
+         CALL isf_tbl_lvl( hu_n, e3u_n, ktop, ikbot, zhtbl, zfrac )
+         !
+         ! compute tbl property at U point
+         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u_n, pvarin, zvarout )
+         !
+         ! compute tbl property at T point
+         pvarout(1,:) = 0._wp
+         DO jj = 1, jpj
+            DO ji = 2, jpi
+               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
+            END DO
+         END DO
+         ! lbclnk not needed as a final communication is done after the computation of fwf
+         ! 
+      CASE ( 'V' )
+         !
+         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
+         zhtbl = phtbl
+         !
+         ! compute tbl lvl and thickness
+         CALL isf_tbl_lvl( hv_n, e3v_n, ktop, ikbot, zhtbl, zfrac )
+         !
+         ! compute tbl property at V point
+         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v_n, pvarin, zvarout )
+         !
+         ! pvarout is an averaging of wet point
+         pvarout(:,1) = 0._wp
+         DO jj = 2, jpj
+            DO ji = 1, jpi
+               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
+            END DO
+         END DO
+         ! lbclnk not needed as a final communication is done after the computation of fwf
+         !
+      CASE ( 'T' )
+         !
+         ! compute tbl property at T point
+         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t_n, pvarin, pvarout )
+         !
+      END SELECT
+      !
+   END SUBROUTINE isf_tbl
+
+   SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_avg  ***
+      !!
+      !! ** Purpose : compute mean property in the boundary layer
+      !!
+      !! ** Method  : Depth average is made between the top level ktop and the bottom level kbot
+      !!              over a thickness phtbl. The bottom level is partially counted (pfrac).
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
+      !!-------------------------- IN  -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop, kbot   ! top and bottom level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3          ! vertical scale factor
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvarin       ! tbl property to average between ktop, kbot over phtbl
+      !!--------------------------------------------------------------------
+      INTEGER  :: ji,jj,jk                    ! loop indices
+      INTEGER  :: ikt, ikb                    ! top and bottom levels
+      !!--------------------------------------------------------------------
+      !
+      ! compute tbl top.bottom level and thickness
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            !
+            ! tbl top/bottom indices initialisation
+            ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
+            !
+            ! level fully include in the ice shelf boundary layer
+            pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj)
+            !
+            ! level partially include in ice shelf boundary layer 
+            pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj)
+            !
+         END DO
+      END DO
+
+   END SUBROUTINE isf_tbl_avg
+
+   SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_lvl  ***
+      !!
+      !! ** Purpose : - compute bottom level off the top boundary layer
+      !!              - thickness of the top boundary layer
+      !!              - fraction of the bottom level affected by the tbl
+      !!
+      !!-------------------------- OUT --------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pfrac  ! fraction of bottom level in the tbl
+      !!-------------------------- IN  --------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phw    ! water column thickness
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
+      !!-------------------------- INOUT ------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) :: phtbl  ! top boundary layer thickness
+      !!---------------------------------------------------------------------
+      INTEGER :: ji,jj,jk
+      INTEGER :: ikt, ikb
+      !!---------------------------------------------------------------------
+      !
+      ! get htbl
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            !
+            ! tbl top/bottom indices initialisation
+            ikt = ktop(ji,jj)
+            !
+            ! limit the tbl to water thickness.
+            phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) )
+            !
+            ! thickness of boundary layer must be at least the top level thickness
+            phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) )
+            !
+         END DO
+      END DO
+      !
+      ! get ktbl
+      CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
+      !
+      ! get pfrac
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            !
+            ! tbl top/bottom indices initialisation
+            ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
+            !
+            ! proportion of the bottom cell included in ice shelf boundary layer 
+            pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb)
+            !
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_tbl_lvl
+   !
+   SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_bot  ***
+      !!
+      !! ** Purpose : compute bottom level of the isf top boundary layer
+      !!
+      !!-------------------------- OUT -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl  ! top boundary layer thickness
+      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
+      !!--------------------------------------------------------------------
+      INTEGER :: ji, jj
+      INTEGER :: ikt, ikb
+      !!--------------------------------------------------------------------
+      !
+      ! phtbl need to be bounded by water column thickness before
+      ! test: if htbl = water column thickness, should return mbathy
+      ! test: if htbl = 0 should return ktop (phtbl cap to e3t(ji,jj,1))
+      !
+      ! get ktbl
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+            !
+            ! determine the deepest level influenced by the boundary layer
+            ikt = ktop(ji,jj)
+            ikb = ikt
+            DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ;  ikb = ikb + 1 ;  END DO
+            kbot(ji,jj) = ikb - 1
+            !
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_tbl_kbot
+      !
+   SUBROUTINE isf_tbl_ktop(pdep, ktop)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_tbl_top  ***
+      !!
+      !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
+      !!
+      !!-------------------------- OUT -------------------------------------
+      INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
+      !!-------------------------- IN  -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdep        ! top depth of the parametrisation influence
+      !!--------------------------------------------------------------------
+      INTEGER :: ji,jj
+      INTEGER :: ikt
+      !!--------------------------------------------------------------------
+      !
+      ! if we need to recompute the top level at every time stepcompute top level (z*, z~) 
+      ! in case of weak ht_n variation we can assume the top level of htbl to be constant
+      ! => only done using gdepw_0
+      ! be sure pdep is already correctly bounded
+      ! test: this routine run on isfdraft should return mikt
+      ! test: this routine run with pdep = 0 should return 1
+      !
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+            ! comput ktop
+            ikt = 2
+            DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ;  ikt = ikt + 1 ;  END DO
+            ktop(ji,jj) = ikt - 1
+            !
+            ! update pdep
+            pdep(ji,jj) = gdepw_0(ji,jj,ktop(ji,jj))
+         END DO
+      END DO
+      !
+   END SUBROUTINE isf_tbl_ktop
+
+END MODULE isftbl
Index: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 12747)
+++ NEMO/branches/UKMO/NEMO_4.0.2_GO8_package_ENHANCE-02_ISF_nemo/src/OCE/ISF/isfutils.F90	(revision 12747)
@@ -0,0 +1,140 @@
+MODULE isfutils
+   !!======================================================================
+   !!                       ***  MODULE  isfutils  ***
+   !! istutils module : miscelenious useful routines
+   !!======================================================================
+   !! History :  4.1  !  2019-09  (P. Mathiot) original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   isfutils       : - read_2dcstdta to read a constant input file with iom_get
+   !!                    - debug to print array sum, min, max in ocean.output
+   !!----------------------------------------------------------------------
+
+   USE iom           , ONLY: iom_open, iom_get, iom_close, jpdom_data ! read input file
+   USE lib_fortran   , ONLY: glob_sum, glob_min, glob_max             ! compute global value
+   USE par_oce       , ONLY: jpi,jpj,jpk                              ! domain size
+   USE dom_oce       , ONLY: nldi, nlei, nldj, nlej                   ! local domain
+   USE in_out_manager, ONLY: i8, wp, lwp, numout                          ! miscelenious
+   USE lib_mpp
+
+   IMPLICIT NONE
+
+   PRIVATE
+
+   INTERFACE debug
+      MODULE PROCEDURE debug2d, debug3d
+   END INTERFACE debug
+
+   PUBLIC read_2dcstdta, debug
+
+CONTAINS
+
+   SUBROUTINE read_2dcstdta(cdfile, cdvar, pvar)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE read_2dcstdta  ***
+      !!
+      !! ** Purpose : read input file
+      !!
+      !!-------------------------- OUT -------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pvar          ! output variable
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(len=256)          , INTENT(in   ) :: cdfile   ! input file name
+      CHARACTER(len=34)           , INTENT(in   ) :: cdvar    ! variable name
+      !!--------------------------------------------------------------------
+      INTEGER :: inum
+      !!--------------------------------------------------------------------
+
+      CALL iom_open( TRIM(cdfile), inum )
+      CALL iom_get( inum, jpdom_data, TRIM(cdvar), pvar)
+      CALL iom_close(inum)
+
+   END SUBROUTINE read_2dcstdta
+
+   SUBROUTINE debug2d(cdtxt,pvar)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_debug2d  ***
+      !!
+      !! ** Purpose : add debug print for 2d variables
+      !!
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(LEN=*)          , INTENT(in   ) :: cdtxt
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pvar
+      !!--------------------------------------------------------------------
+      REAL(wp)    :: zmin, zmax, zsum
+      INTEGER(i8) :: imodd, ip
+      INTEGER     :: itmps,imods, ji, jj, jk
+      !!--------------------------------------------------------------------
+      !
+      ! global min/max/sum to check data range and NaN
+      zsum = glob_sum( 'debug', pvar(:,:) )
+      zmin = glob_min( 'debug', pvar(:,:) )
+      zmax = glob_max( 'debug', pvar(:,:) )
+      !
+      ! basic check sum to check reproducibility
+      ! TRANSFER function find out the integer corresponding to pvar(i,j) bit pattern
+      ! MOD allow us to keep only the latest digits during the sum
+      ! imod is not choosen to be very large as at the end there is a classic mpp_sum
+      imodd=65521 ! highest prime number < 2**16 with i8 type
+      imods=65521 ! highest prime number < 2**16 with default integer for mpp_sum subroutine
+      itmps=0
+      DO jj=nldj,nlej
+         DO ji=nldi,nlei
+            itmps = MOD(itmps + MOD(TRANSFER(pvar(ji,jj), ip),imodd), imods)
+         END DO
+      END DO
+      CALL mpp_sum('debug',itmps)
+      !
+      ! print out
+      IF (lwp) THEN
+         WRITE(numout,*) TRIM(cdtxt),' (min, max, sum, tag) : ',zmin, zmax, zsum, itmps
+         CALL FLUSH(numout)
+      END IF
+      !
+   END SUBROUTINE debug2d
+
+   SUBROUTINE debug3d(cdtxt,pvar)
+      !!--------------------------------------------------------------------
+      !!                  ***  ROUTINE isf_debug3d  ***
+      !!
+      !! ** Purpose : add debug print for 3d variables
+      !!
+      !!-------------------------- IN  -------------------------------------
+      CHARACTER(LEN=*)              , INTENT(in   ) :: cdtxt
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvar
+      !!--------------------------------------------------------------------
+      REAL(wp) :: zmin, zmax, zsum
+      INTEGER(i8) :: imodd, ip
+      INTEGER     :: itmps,imods, ji, jj, jk
+      !!--------------------------------------------------------------------
+      !
+      ! global min/max/sum to check data range and NaN
+      zsum = glob_sum( 'debug', pvar(:,:,:) )
+      zmin = glob_min( 'debug', pvar(:,:,:) )
+      zmax = glob_max( 'debug', pvar(:,:,:) )
+      !
+      ! basic check sum to check reproducibility
+      ! TRANSFER function find out the integer corresponding to pvar(i,j) bit pattern
+      ! MOD allow us to keep only the latest digits during the sum
+      ! imod is not choosen to be very large as at the end there is a classic mpp_sum
+      imodd=65521 ! highest prime number < 2**16 with i8 type
+      imods=65521 ! highest prime number < 2**16 with default integer for mpp_sum subroutine
+      itmps=0
+      DO jk=1,jpk
+         DO jj=nldj,nlej
+            DO ji=nldi,nlei
+               itmps = MOD(itmps + MOD(TRANSFER(pvar(ji,jj,jk), ip),imodd), imods)
+            END DO
+         END DO
+      END DO
+      CALL mpp_sum('debug',itmps)
+      !
+      ! print out
+      IF (lwp) THEN
+         WRITE(numout,*) TRIM(cdtxt),' (min, max, sum, tag) : ',zmin, zmax, zsum, itmps
+         CALL FLUSH(numout)
+      END IF
+      !
+   END SUBROUTINE debug3d
+
+END MODULE isfutils
