Index: /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12077)
@@ -264,32 +264,32 @@
 
           <!-- * variable related to ice shelf forcing * -->
-          <field id="isftfrz_par"     long_name="fzp temperature at ocean/isf interface"        unit="degC"     />
-          <field id="isftfrz_cav"     long_name="fzp temperature at ocean/isf interface"        unit="degC"     />
-          <field id="fwfisf_cav"      long_name="Ice shelf melting"                             unit="kg/m2/s"  />
-          <field id="fwfisf_par"      long_name="Ice shelf melting"                             unit="kg/m2/s"  />
+          <field id="isftfrz_cav"     long_name="freezing point temperature at ocean/isf interface"                unit="degC"     />
+          <field id="isftfrz_par"     long_name="freezing point temperature in the parametrization boundary layer" unit="degC"     />
+          <field id="fwfisf_cav"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
+          <field id="fwfisf_par"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
           <field id="qoceisf_cav"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
           <field id="qoceisf_par"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
           <field id="qlatisf_cav"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
           <field id="qlatisf_par"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
-          <field id="qhcisf_cav"      long_name="Ice shelf heat content flux"                   unit="W/m2"     />
-          <field id="qhcisf_par"      long_name="Ice shelf heat content flux"                   unit="W/m2"     />
-          <field id="fwfisf3d_cav"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" />
-          <field id="fwfisf3d_par"    long_name="Ice shelf melting"                             unit="kg/m2/s"  grid_ref="grid_T_3D" />
+          <field id="qhcisf_cav"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
+          <field id="qhcisf_par"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
+          <field id="fwfisf3d_cav"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
+          <field id="fwfisf3d_par"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
           <field id="qoceisf3d_cav"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
           <field id="qoceisf3d_par"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
           <field id="qlatisf3d_cav"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
           <field id="qlatisf3d_par"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
-          <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" />
-          <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux"                   unit="W/m2"     grid_ref="grid_T_3D" />
-          <field id="ttbl_cav"        long_name="temperature in the tracer sample depth      "  unit="C"        />
-          <field id="ttbl_par"        long_name="temperature in the tracer sample depth      "  unit="C"        />
-          <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting        "  unit="C"        />
-          <field id="isfthermald_par" long_name="thermal driving of ice shelf melting        "  unit="C"        />
-          <field id="isfgammat"       long_name="transfert coefficient for isf (temperature) "  unit="m/s"      />
-          <field id="isfgammas"       long_name="transfert coefficient for isf (salinity)    "  unit="m/s"      />
-          <field id="stbl"            long_name="salinity in the Losh tbl                    "  unit="PSU"      />
-          <field id="utbl"            long_name="zonal current in the Losh tbl at T point    "  unit="m/s"      />
-          <field id="vtbl"            long_name="merid current in the Losh tbl at T point    "  unit="m/s"      />
-          <field id="isfustar"        long_name="ustar at T point used in ice shelf melting  "  unit="m/s"      />
+          <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
+          <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
+          <field id="ttbl_cav"        long_name="temperature in Losch tbl"                      unit="degC"     />
+          <field id="ttbl_par"        long_name="temperature in the parametrisation boundary layer" unit="degC" />
+          <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting"          unit="degC"     />
+          <field id="isfthermald_par" long_name="thermal driving of ice shelf melting"          unit="degC"     />
+          <field id="isfgammat"       long_name="Ice shelf heat-transfert velocity"             unit="m/s"      />
+          <field id="isfgammas"       long_name="Ice shelf salt-transfert velocity"             unit="m/s"      />
+          <field id="stbl"            long_name="salinity in the Losh tbl"                      unit="1e-3"     />
+          <field id="utbl"            long_name="zonal current in the Losh tbl at T point"      unit="m/s"      />
+          <field id="vtbl"            long_name="merid current in the Losh tbl at T point"      unit="m/s"      />
+          <field id="isfustar"        long_name="ustar at T point used in ice shelf melting"    unit="m/s"      />
           <field id="qconisf"         long_name="Conductive heat flux through the ice shelf"    unit="W/m2"     />
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/namelist_ref	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/SHARED/namelist_ref	(revision 12077)
@@ -439,5 +439,7 @@
    ! ---------------- ice shelf load -------------------------------
    !
-   cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
+   cn_isfload = 'uniform'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
+      rn_isfload_T = -1.9
+      rn_isfload_S = 34.4
    !
    ! ---------------- ice shelf melt formulation -------------------------------
@@ -452,13 +454,14 @@
          cn_isfcav_mlt = '3eq'   ! ice shelf melting formulation (spe/2eq/3eq/oasis)
          !                       ! spe = fwfisf is read from a forcing field
-         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006)
-         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015)
+         !                       ! 2eq = ISOMIP  like: 2 equations formulation (Hunter et al., 2006 for a short description)
+         !                       ! 3eq = ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2016 for a short description)
          !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf
          !              !  cn_isfcav_mlt = 2eq or 3eq cases:
-         cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99)
-         !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
-         !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999)
-         rn_gammat0  = 1.4e-2    ! gammat coefficient used in blk formula
-         rn_gammas0  = 4.e-4    ! gammas coefficient used in blk formula
+         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99)
+         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0)
+         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description)
+         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description)
+         rn_gammat0  = 1.4e-2    ! gammat coefficient used in spe, vel and vel_stab gamma computation method
+         rn_gammas0  = 4.0e-4    ! gammas coefficient used in spe, vel and vel_stab gamma computation method
          !
          rn_htbl     =  30.      ! thickness of the top boundary layer    (Losh et al. 2008)
Index: /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/WED025/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/WED025/EXPREF/namelist_cfg	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/cfgs/WED025/EXPREF/namelist_cfg	(revision 12077)
@@ -211,8 +211,4 @@
 !-----------------------------------------------------------------------
    !
-   ! ---------------- ice shelf load -------------------------------
-   !
-   cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
-   !
    ! ---------------- ice shelf melt formulation -------------------------------
    !
@@ -229,5 +225,5 @@
          !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf
          !              !  cn_isfcav_mlt = 2eq or 3eq cases:
-         cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99)
+         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99)
          !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
          !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999)
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/BDY/bdyvol.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/BDY/bdyvol.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/BDY/bdyvol.F90	(revision 12077)
@@ -14,7 +14,7 @@
    USE bdy_oce        ! ocean open boundary conditions
    USE sbc_oce        ! ocean surface boundary conditions
+   USE isf_oce, ONLY : fwfisf_cav, fwfisf_par  ! ice shelf
    USE dom_oce        ! ocean space and time domain 
    USE phycst         ! physical constants
-   USE isf            ! ice shelf
    !
    USE in_out_manager ! I/O manager
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diahsb.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diahsb.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diahsb.F90	(revision 12077)
@@ -17,6 +17,6 @@
    USE phycst         ! physical constants
    USE sbc_oce        ! surface thermohaline fluxes
+   USE isf_oce        ! ice shelf fluxes
    USE sbcrnf         ! river runoff
-   USE isf            ! ice shelves
    USE domvvl         ! vertical scale factors
    USE traqsr         ! penetrative solar radiation
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diawri.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DIA/diawri.F90	(revision 12077)
@@ -26,5 +26,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers 
-   USE isf
+   USE isf_oce
    USE isfcpl
    USE dom_oce        ! ocean space and time domain
@@ -911,4 +911,5 @@
       CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
       CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height
+
       IF ( ln_isf ) THEN
          IF (ln_isfcav_mlt) THEN
@@ -918,4 +919,5 @@
             CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
          END IF
          IF (ln_isfpar_mlt) THEN
@@ -926,12 +928,8 @@
             CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
             CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
+            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
          END IF
       END IF
 
-      IF ( ln_isf ) THEN
-         IF (ln_isfcav_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
-         IF (ln_isfpar_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
-      END IF
-      
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DOM/domwri.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DOM/domwri.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DOM/domwri.F90	(revision 12077)
@@ -16,5 +16,5 @@
    !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
    !!----------------------------------------------------------------------
-   USE isf             ! ice shelf
+   !
    USE dom_oce         ! ocean space and time domain
    USE phycst ,   ONLY :   rsmall
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/divhor.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/divhor.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/divhor.F90	(revision 12077)
@@ -19,10 +19,9 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
-   USE isf
-   USE isfutils
    USE dom_oce         ! ocean space and time domain
-   USE sbc_oce, ONLY : ln_rnf ! surface boundary condition: ocean
-   USE sbcrnf          ! river runoff 
-   USE isfhdiv         ! ice shelf
+   USE sbc_oce, ONLY : ln_rnf      ! river runoff
+   USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 
+   USE isf_oce, ONLY : ln_isf      ! ice shelf
+   USE isfhdiv, ONLY : isf_hdiv    ! ice shelf
 #if defined key_asminc   
    USE asminc          ! Assimilation increment
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynatf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynatf.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynatf.F90	(revision 12077)
@@ -41,6 +41,6 @@
    USE trddyn         ! trend manager: dynamics
    USE trdken         ! trend manager: kinetic energy
-   USE isf       , ONLY: ln_isf     ! ice shelf
-   USE isfdynatf , ONLY: isf_dynatf ! ice shelf 
+   USE isf_oce   , ONLY: ln_isf     ! ice shelf
+   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 
    !
    USE in_out_manager ! I/O manager
@@ -223,4 +223,5 @@
             ! PM: we could probably define a generic subroutine to do the in depth correction
             !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
+            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
             IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt )
             !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynhpg.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynhpg.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynhpg.F90	(revision 12077)
@@ -31,6 +31,6 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
-   USE isf             ! ice shelf  (risfload variable)
-   USE isfload         ! ice shelf  (isf_load routine )
+   USE isf_oce , ONLY : risfload  ! ice shelf  (risfload variable)
+   USE isfload , ONLY : isf_load  ! ice shelf  (isf_load routine )
    USE sbc_oce         ! surface variable (only for the flag with ice shelf)
    USE dom_oce         ! ocean space and time domain
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/dynspg_ts.F90	(revision 12077)
@@ -31,8 +31,7 @@
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! surface boundary condition: ocean
+   USE isf_oce         ! ice shelf variable (fwfisf)
    USE zdf_oce         ! vertical physics: variables
    USE zdfdrg          ! vertical physics: top/bottom drag coef.
-   USE isf             ! ice shelf variable (fwfisf)
-   USE isfutils
    USE sbcapr          ! surface boundary condition: atmospheric pressure
    USE dynadv    , ONLY: ln_dynadv_vec
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/sshwzv.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/sshwzv.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DYN/sshwzv.F90	(revision 12077)
@@ -19,6 +19,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
-   USE isf            ! ice shelf
-   USE isfutils
+   USE isf_oce        ! ice shelf
    USE dom_oce        ! ocean space and time domain variables 
    USE sbc_oce        ! surface boundary condition: ocean
Index: MO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isf.F90	(revision 12076)
+++ 	(revision )
@@ -1,249 +1,0 @@
-MODULE isf
-   !!======================================================================
-   !!                       ***  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
-   !
-   !-------------------------------------------------------
-   ! 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]
-   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
-   !-------------------------------------------------------
-   !
-   ! 0.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
-   !
-   ! 0.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
-   !
-   ! 0.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                      !:
-   !
-   ! 0.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: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
-   !! 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 : 
-      !!
-      !! ** Method  : 
-      !!
-      !!----------------------------------------------------------------------
-      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_alloc()
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE isf_alloc  ***
-      !!
-      !! ** Purpose : 
-      !!
-      !! ** Method  : 
-      !!
-      !!----------------------------------------------------------------------
-      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
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isf_oce.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isf_oce.F90	(revision 12077)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isf_oce.F90	(revision 12077)
@@ -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: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
+   !! 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/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcav.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcav.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcav.F90	(revision 12077)
@@ -13,5 +13,5 @@
    !!   isf_cav       : update ice shelf melting under ice shelf
    !!----------------------------------------------------------------------
-   USE isf            ! ice shelf public variables
+   USE isf_oce        ! ice shelf public variables
    !
    USE isfrst   , ONLY: isfrst_write, isfrst_read ! ice shelf restart read/write subroutine
@@ -25,5 +25,5 @@
    USE par_oce  , ONLY: jpi,jpj                ! ocean space and time domain
    USE phycst   , ONLY: grav,rau0,r1_rau0_rcp  ! physical constants
-   USE eosbn2   , ONLY: l_useCT                ! l_useCT
+   USE eosbn2   , ONLY: ln_teos10               ! use ln_teos10 or not
    !
    USE in_out_manager ! I/O manager
@@ -85,5 +85,5 @@
       !
       ! initialisation
-      IF (TRIM(cn_gammablk) == 'HJ99' ) zqoce_b (:,:) = ptsc(:,:,jp_tem) * rau0_rcp ! last time step total heat fluxes (to speed up convergence)
+      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
@@ -104,8 +104,8 @@
          ! define if we need to iterate
          SELECT CASE ( cn_gammablk )
-         CASE ( 'spe','ad15' )
+         CASE ( 'spe','vel' )
             ! no convergence needed
             lit = .FALSE.
-         CASE ( 'hj99' )
+         CASE ( 'vel_stab' )
             ! compute error between 2 iterations
             zerr = MAXVAL(ABS(zqoce(:,:) - zqoce_b(:,:)))
@@ -113,5 +113,5 @@
             ! define if iteration needed
             IF (nit >= 100) THEN              ! too much iteration
-               CALL ctl_stop( 'STOP', 'isf_cav: HJ99 gamma formulation had too many iterations ...' )
+               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.
@@ -212,5 +212,5 @@
          ! coeficient for linearisation of potential tfreez
          ! Crude approximation for pressure (but commonly used)
-         IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
+         IF ( ln_teos10 ) THEN   ! linearisation from Jourdain et al. (2017)
             risf_lamb1 =-0.0564_wp
             risf_lamb2 = 0.0773_wp
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavgam.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavgam.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavgam.F90	(revision 12077)
@@ -10,5 +10,5 @@
    !!   isfcav_gammats       : compute exchange coeficient gamma 
    !!----------------------------------------------------------------------
-   USE isf
+   USE isf_oce
    USE isfutils, ONLY: debug
    USE isftbl  , ONLY: isf_tbl
@@ -46,5 +46,5 @@
       !!
       !! ** Method     : select the gamma formulation
-      !!                 3 method available (cst, AD15 and HJ99)
+      !!                 3 method available (cst, vel and vel_stab)
       !!---------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
@@ -66,5 +66,5 @@
          ! gamma is constant (specified in namelist)
          ! nothing to do
-      CASE ('ad15', 'hj99')
+      CASE ('vel', 'vel_stab')
          ! compute velocity in tbl
          CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
@@ -90,8 +90,8 @@
          pgt(:,:) = rn_gammat0
          pgs(:,:) = rn_gammas0
-      CASE ( 'ad15' ) ! gamma is proportional to u*
-         CALL gammats_AD15 (                   zutbl, zvtbl, rCd0_top, r_ke0_top,               pgt, pgs )
-      CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_HJ99 (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pgt, pgs )
+      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 (Kmm, 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)')
@@ -116,5 +116,5 @@
    !!-----------------------------------------------------------------------------------------------------
    !
-   SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in
+   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in
       &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
       !!----------------------------------------------------------------------
@@ -123,6 +123,5 @@
       !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
       !!
-      !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312
-      !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
+      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
       !!---------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
@@ -146,8 +145,8 @@
       CALL iom_put('isfustar',zustar(:,:))
       !
-   END SUBROUTINE gammats_AD15
-
-   SUBROUTINE gammats_HJ99( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
-      &                                                                 pgt  , pgs    )  ! ==>> out gammats [m/s]
+   END SUBROUTINE gammats_vel
+
+   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
+      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s]
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange coefficient 
@@ -186,5 +185,5 @@
       !
       ! compute ustar
-      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + r_ke0_top ) )
+      zustar(:,:) = SQRT( pCd * ( putbl(:,:) * putbl(:,:) + pvtbl(:,:) * pvtbl(:,:) + pke2 ) )
       !
       ! output ustar
@@ -254,5 +253,5 @@
       END DO
 
-   END SUBROUTINE gammats_HJ99
+   END SUBROUTINE gammats_vel_stab
 
 END MODULE isfcavgam
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavmlt.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavmlt.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcavmlt.F90	(revision 12077)
@@ -12,5 +12,5 @@
    !!----------------------------------------------------------------------
 
-   USE isf                      ! ice shelf
+   USE isf_oce                  ! ice shelf
    USE isftbl , ONLY: isf_tbl   ! ice shelf depth average
    USE isfutils,ONLY: debug     ! debug subroutine
@@ -77,7 +77,7 @@
       !
       IF (ln_isfdebug) THEN
-         CALL debug( 'isfcav_mlt:', pqhc (:,:) )
-         CALL debug( 'isfcav_mlt:', pqoce(:,:) )
-         CALL debug( 'isfcav_mlt:', pqfwf(:,:) )
+         CALL debug( 'isfcav_mlt qhc  :', pqhc (:,:) )
+         CALL debug( 'isfcav_mlt qoce :', pqoce(:,:) )
+         CALL debug( 'isfcav_mlt qfwf :', pqfwf(:,:) )
       END IF
       !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfcpl.F90	(revision 12077)
@@ -12,5 +12,5 @@
    !!   isfrst : read/write iceshelf variables in/from restart
    !!----------------------------------------------------------------------
-   USE isf                              ! ice shelf variable
+   USE isf_oce                          ! ice shelf variable
    USE isfutils, ONLY : debug
    USE lib_mpp , ONLY: mpp_sum, mpp_max ! mpp routine
@@ -25,6 +25,6 @@
    PRIVATE
 
-   PUBLIC isfcpl_rst_write, isfcpl_init  ! iceshelf restart read and write 
-   PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons
+   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
@@ -48,4 +48,13 @@
    SUBROUTINE isfcpl_init(Kbb, Kmm, Kaa)
       !!---------------------------------------------------------------------
+      !!                   ***  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
       !!
       !!---------------------------------------------------------------------
@@ -480,5 +489,5 @@
       INTEGER  ::   jisf                              ! start, end and current position in the increment array
       INTEGER  ::   ingb, ifind                       ! 0/1 target found or need to be found 
-      INTEGER  ::   ingb, ifind, nisfl_area           ! global number of cell concerned by the wet->dry case 
+      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
       !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdiags.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdiags.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdiags.F90	(revision 12077)
@@ -1,7 +1,6 @@
 MODULE isfdiags
    !!======================================================================
-   !!                       ***  MODULE  sbcisf  ***
-   !! Surface module :  update surface ocean boundary condition under ice
-   !!                   shelf
+   !!                       ***  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
@@ -16,6 +15,5 @@
    USE in_out_manager ! I/O manager
    USE dom_oce
-   USE isf            ! ice shelf variable
-   USE isfutils
+   USE isf_oce        ! ice shelf variable
    USE iom            ! 
 
@@ -38,5 +36,6 @@
       !!                  ***  ROUTINE isf_diags_flx ***
       !!
-      !! ** Purpose : 
+      !! ** Purpose : manage the 2d and 3d flux outputs of the ice shelf module
+      !!              fwf, latent heat flux, heat content flux, oce->ice heat flux
       !!
       !!----------------------------------------------------------------------
@@ -59,9 +58,9 @@
       cvarqhc  = 'qhcisf_'//cdisf  ; cvarqhc3d  = 'qhcisf3d_'//cdisf
       !
-      ! output 2d melt rate, latent heat and heat content flux
-      CALL iom_put( TRIM(cvarqfwf), pqfwf(:,:) )   ! isf mass flux (opposite sign)
-      CALL iom_put( TRIM(cvarqoce), pqoce(:,:) )   ! isf oce to ice   flux  (cpo*gt*dt)
-      CALL iom_put( TRIM(cvarqlat), pqlat(:,:) )   ! isf oce to ice   flux  (cpo*gt*dt)
-      CALL iom_put( TRIM(cvarqhc) , pqhc (:,:) )   ! isf heat content flux  (cpo*fwf*Tfrz)
+      ! 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
@@ -77,5 +76,6 @@
       !!                  ***  ROUTINE isf_diags_2dto3d ***
       !!
-      !! ** Purpose : 
+      !! ** Purpose : compute the 3d flux outputs as they are injected into NEMO 
+      !!              (ie uniformaly spread into the top boundary layer or parametrisation layer)
       !!
       !!----------------------------------------------------------------------
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdynatf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdynatf.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfdynatf.F90	(revision 12077)
@@ -11,5 +11,5 @@
    !!-------------------------------------------------------------------------
 
-   USE isf
+   USE isf_oce
 
    USE phycst , ONLY: r1_rau0         ! physical constant
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfhdiv.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfhdiv.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfhdiv.F90	(revision 12077)
@@ -1,6 +1,18 @@
 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
+   !!----------------------------------------------------------------------
 
-   USE isf                    ! ice shelf
-   USE isfutils 
+   !!----------------------------------------------------------------------
+   !!   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
@@ -40,5 +52,8 @@
          IF ( ln_isfcpl .AND. kt /= 0 ) THEN
             !
-            ! correct divergence only for the first time step
+            ! 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(Kmm, risfcpl_vol       , phdiv)
             IF ( kt == nit000+1 ) CALL isf_hdiv_cpl(Kmm, risfcpl_vol*0.5_wp, phdiv)
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfload.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfload.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfload.F90	(revision 12077)
@@ -11,11 +11,11 @@
    !!----------------------------------------------------------------------
 
-   USE isf, ONLY: cn_isfload   ! ice shelf variables
+   USE isf_oce, ONLY: cn_isfload, rn_isfload_T, rn_isfload_S ! ice shelf variables
 
    USE dom_oce, ONLY: e3w, gdept, risfdep, mikt     ! vertical scale factor
    USE eosbn2 , ONLY: eos                           ! eos routine
 
-   USE lib_mpp, ONLY: ctl_stop ! ctl_stop routine
-   USE in_out_manager  ! 
+   USE lib_mpp, ONLY: ctl_stop                               ! ctl_stop routine
+   USE in_out_manager                                        ! 
 
    IMPLICIT NONE
@@ -45,6 +45,6 @@
       ! ice shelf cavity
       SELECT CASE ( cn_isfload )
-      CASE ( 'isomip' )
-         CALL isf_load_isomip ( Kmm, pisfload )
+      CASE ( 'uniform' )
+         CALL isf_load_uniform ( Kmm, pisfload )
       CASE DEFAULT
          CALL ctl_stop('STOP','method cn_isfload to compute ice shelf load does not exist (isomip), check your namelist')
@@ -53,5 +53,5 @@
    END SUBROUTINE isf_load
 
-   SUBROUTINE isf_load_isomip( Kmm, pisfload )
+   SUBROUTINE isf_load_uniform( Kmm, pisfload )
       !!--------------------------------------------------------------------
       !!                  ***  SUBROUTINE isf_load  ***
@@ -79,6 +79,6 @@
       znad = 1._wp                     !- To use density and not density anomaly
       !
-      !                                !- assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
-      zts_top(:,:,jp_tem) = -1.9_wp   ;   zts_top(:,:,jp_sal) = 34.4_wp
+      !                                !- 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 
@@ -113,5 +113,5 @@
       END DO
       !
-   END SUBROUTINE isf_load_isomip
+   END SUBROUTINE isf_load_uniform
 
 END MODULE isfload
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfpar.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfpar.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfpar.F90	(revision 12077)
@@ -14,11 +14,11 @@
    !!   isfpar       : compute ice shelf melt using a prametrisation of ice shelf cavities
    !!----------------------------------------------------------------------
-   USE isf            ! ice shelf
+   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 isfparmlt, ONLY: isfpar_mlt     ! ice shelf melt formulation subroutine
-   USE isfdiags , ONLY: isf_diags_flx  ! ice shelf diags subroutine
    !
    USE dom_oce  , ONLY: bathy          ! ocean space and time domain
@@ -48,4 +48,11 @@
       !!
       !! ** 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.
       !!
       !!---------------------------------------------------------------------
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfparmlt.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfparmlt.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfparmlt.F90	(revision 12077)
@@ -8,5 +8,5 @@
    !!----------------------------------------------------------------------
 
-   USE isf            ! ice shelf
+   USE isf_oce                  ! ice shelf
    USE isftbl , ONLY: isf_tbl   ! ice shelf depth average
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfstp.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfstp.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfstp.F90	(revision 12077)
@@ -14,5 +14,5 @@
    !!----------------------------------------------------------------------
    !
-   USE isf            ! isf variables
+   USE isf_oce                                      ! isf variables
    USE isfload, ONLY: isf_load                      ! ice shelf load
    USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer
@@ -187,12 +187,12 @@
             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) 
+               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 drag coef.    used (from namdrg_top) rn_Cd0      = ', r_Cdmin_top
-                  WRITE(numout,*) '         top background ke used (from namdrg_top) rn_ke0      = ', r_ke0_top
+                  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
@@ -222,6 +222,11 @@
          END IF
 
-         IF (ln_isfcav) WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload)
+         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
@@ -245,5 +250,5 @@
       IF ( l_isfoasis .AND. ln_isf ) THEN
          !
-         CALL ctl_stop( ' OASIS and ice shelf not tested' )
+         CALL ctl_stop( ' ln_ctl and ice shelf not tested' )
          !
          ! NEMO coupled to ATMO model with isf cavity need oasis method for melt computation 
@@ -277,9 +282,11 @@
       INTEGER               :: ios                  ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
-      NAMELIST/namisf/ ln_isf       ,                                                                               & 
-         &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, 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   , cn_isfdir
+      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
       !!----------------------------------------------------------------------
       !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isftbl.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isftbl.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isftbl.F90	(revision 12077)
@@ -14,5 +14,5 @@
    !!----------------------------------------------------------------------
 
-   USE isf     ! ice shelf variables
+   USE isf_oce ! ice shelf variables
 
    USE dom_oce ! vertical scale factor and depth
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfutils.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfutils.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ISF/isfutils.F90	(revision 12077)
@@ -15,5 +15,7 @@
    USE lib_fortran   , ONLY: glob_sum, glob_min, glob_max             ! compute global value
    USE par_oce       , ONLY: jpi,jpj,jpk                              ! domain size
-   USE in_out_manager, ONLY: wp, lwp, numout                          ! miscelenious
+   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
@@ -54,20 +56,39 @@
       !!                  ***  ROUTINE isf_debug2d  ***
       !!
-      !! ** Purpose : add debug print
+      !! ** Purpose : add debug print for 2d variables
       !!
       !!-------------------------- IN  -------------------------------------
-      CHARACTER(LEN=256)          , INTENT(in   ) :: cdtxt
+      CHARACTER(LEN=*)            , INTENT(in   ) :: cdtxt
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pvar
       !!--------------------------------------------------------------------
-      REAL(wp) :: zmin, zmax, zsum
+      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(:,:) )
       !
-      IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      ! 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)
       !
-      CALL FLUSH(numout)
+      ! 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
@@ -77,23 +98,43 @@
       !!                  ***  ROUTINE isf_debug3d  ***
       !!
-      !! ** Purpose : add debug print
+      !! ** Purpose : add debug print for 3d variables
       !!
       !!-------------------------- IN  -------------------------------------
-      CHARACTER(LEN=256)              , INTENT(in   ) :: cdtxt
+      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(:,:,:) )
       !
-      IF (lwp) WRITE(numout,*) TRIM(cdtxt),' (min, max, sum) : ',zmin, zmax, zsum
+      ! 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)
       !
-      CALL FLUSH(numout)
+      ! 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
-
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldfslp.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldfslp.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldfslp.F90	(revision 12077)
@@ -21,5 +21,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
-   USE isf            ! ice shelf
+   USE isf_oce        ! ice shelf
    USE dom_oce        ! ocean space and time domain
 !   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldftra.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldftra.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/LDF/ldftra.F90	(revision 12077)
@@ -664,5 +664,5 @@
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
                   zhw(ji,jj) = zhw(ji,jj) + ze3w
@@ -682,5 +682,5 @@
                   ! eddies using the isopycnal slopes calculated in ldfslp.F : 
                   ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-                  ze3w = e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
+                  ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbccpl.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbccpl.F90	(revision 12077)
@@ -27,4 +27,5 @@
    USE sbcwave         ! surface boundary condition: waves
    USE phycst          ! physical constants
+   USE isf_oce , ONLY : l_isfoasis, fwfisf_oasis ! ice shelf boundary condition
 #if defined key_si3
    USE ice            ! ice variables
@@ -36,5 +37,4 @@
    USE eosbn2         ! 
    USE sbcrnf  , ONLY : l_rnfcpl
-   USE isf     , ONLY : ln_isf, l_isfoasis, fwfisf_oasis
 #if defined key_cice
    USE ice_domain_size, only: ncat
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcfwb.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcfwb.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcfwb.F90	(revision 12077)
@@ -17,8 +17,8 @@
    USE dom_oce        ! ocean space and time domain
    USE sbc_oce        ! surface ocean boundary condition
+   USE isf_oce , ONLY : fwfisf_cav, fwfisf_par                    ! ice shelf melting contribution
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass
    USE phycst         ! physical constants
    USE sbcrnf         ! ocean runoffs
-   USE isf            ! ice shelf melting contribution
    USE sbcssr         ! Sea-Surface damping terms
    !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcrnf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcrnf.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcrnf.F90	(revision 12077)
@@ -19,5 +19,4 @@
    USE phycst         ! physical constants
    USE sbc_oce        ! surface boundary condition variables
-   USE isf            ! ice shelf
    USE eosbn2         ! Equation Of State
    USE closea         ! closed seas
@@ -127,7 +126,4 @@
                rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0
             END WHERE
-            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg
-               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rLfusisf * r1_rau0_rcp
-            END WHERE
          ELSE                                                        ! use SST as runoffs temperature
             !CEOD River is fresh water so must at least be 0 unless we consider ice
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcssm.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcssm.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/SBC/sbcssm.F90	(revision 12077)
@@ -61,10 +61,6 @@
       !
       !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zts(ji,jj,jp_tem) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm)
-            zts(ji,jj,jp_sal) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm)
-         END DO
-      END DO
+      zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
+      zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
       !
       IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        !
@@ -73,5 +69,5 @@
          ssv_m(:,:) = vv(:,:,1,Kbb)
          IF( l_useCT )  THEN    ;   sst_m(:,:) = eos_pt_from_ct( zts(:,:,jp_tem), zts(:,:,jp_sal) )
-         ELSE                    ;   sst_m(:,:) = zts(:,:,jp_tem)
+         ELSE                   ;   sst_m(:,:) = zts(:,:,jp_tem)
          ENDIF
          sss_m(:,:) = zts(:,:,jp_sal)
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/eosbn2.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/eosbn2.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/eosbn2.F90	(revision 12077)
@@ -29,5 +29,4 @@
    !!   eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass
    !!   eos_insitu_2d : Compute the in situ density for 2d fields
-   !!   bn2           : Compute the Brunt-Vaisala frequency
    !!   bn2           : compute the Brunt-Vaisala frequency
    !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traatf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traatf.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traatf.F90	(revision 12077)
@@ -29,5 +29,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbcrnf          ! river runoffs
-   USE isf             ! ice shelf melting
+   USE isf_oce         ! ice shelf melting
    USE zdf_oce         ! ocean vertical mixing
    USE domvvl          ! variable volume
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traisf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traisf.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/TRA/traisf.F90	(revision 12077)
@@ -10,14 +10,9 @@
    !!   tra_isf       : update the tracer trend at ocean surface
    !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers
-   USE dom_oce        ! ocean space domain variables
-   USE phycst         ! physical constant
-   USE eosbn2         ! Equation Of State
-   USE isf            ! Ice shelf variable
-   USE isfutils       !
-   !
-   USE in_out_manager ! I/O manager
-   USE iom            ! xIOS server
-   USE timing         ! Timing
+   USE isf_oce                                     ! Ice shelf variables
+   USE dom_oce , ONLY : e3t, r1_e1e2t            ! ocean space domain variables
+   USE isfutils, ONLY : debug                      ! debug option
+   USE timing  , ONLY : timing_start, timing_stop  ! Timing
+   USE in_out_manager                              ! I/O manager
 
    IMPLICIT NONE
@@ -65,4 +60,10 @@
       IF ( ln_isfcpl ) 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).
+         ! in accordance to this, the heat content flux due to injected water need to be added in the temperature and salt trend
+         ! at time step nit000 and nit000+1
          IF ( kt == nit000  ) CALL tra_isf_cpl(Kmm, risfcpl_tsc       , pts(:,:,:,:,Krhs))
          IF ( kt == nit000+1) CALL tra_isf_cpl(Kmm, risfcpl_tsc*0.5_wp, pts(:,:,:,:,Krhs))
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ZDF/zdfmxl.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ZDF/zdfmxl.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/ZDF/zdfmxl.F90	(revision 12077)
@@ -12,5 +12,5 @@
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
-   USE isf            ! ice shelf
+   USE isf_oce        ! ice shelf
    USE dom_oce        ! ocean space and time domain variables
    USE trc_oce  , ONLY: l_offline         ! ocean space and time domain variables
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step.F90	(revision 12077)
@@ -115,10 +115,10 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Update external forcing (tides, open boundaries, and surface boundary condition (including sea-ice)
-      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_isf     )   CALL isf_stp ( kstp, Nnn )
+      ! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( ln_tide    )   CALL sbc_tide( kstp )                        ! update tide potential
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                        ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
       IF( ln_bdy     )   CALL bdy_dta ( kstp, Nnn, kt_offset = +1 )   ! update dynamic & tracer data at open boundaries
+      IF( ln_isf     )   CALL isf_stp ( kstp, Nnn )
                          CALL sbc     ( kstp, Nbb, Nnn )              ! Sea Boundary Condition (including sea-ice)
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step_oce.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step_oce.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/step_oce.F90	(revision 12077)
@@ -22,5 +22,5 @@
    USE sbcwave         ! Wave intialisation
 
-   USE isf
+   USE isf_oce         ! ice shelf boundary condition
    USE isfstp          ! ice shelf boundary condition     (isf_stp routine)
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/EXPREF/namelist_cfg	(revision 12077)
@@ -165,5 +165,7 @@
    ! ---------------- ice shelf load -------------------------------
    !
-   cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
+   cn_isfload = 'uniform'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
+      rn_isfload_T = -1.0
+      rn_isfload_S =  34.2
    !
    ! ---------------- ice shelf melt formulation -------------------------------
@@ -182,7 +184,8 @@
          !                       ! oasis = fwfisf is given by oasis and pattern by file sn_isfcav_fwf
          !              !  cn_isfcav_mlt = 2eq or 3eq cases:
-         cn_gammablk = 'ad15'    ! scheme to compute gammat/s (spe,ad15,hj99)
-         !                       ! ad15 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010)
-         !                       ! hj99 = velocity and stability dependent Gamma    (Holland et al. 1999)
+         cn_gammablk = 'vel'     ! scheme to compute gammat/s (spe,ad15,hj99)
+         !                       ! spe      = constant transfert velocity (rn_gammat0, rn_gammas0)
+         !                       ! vel      = velocity dependent transfert velocity (u* * gammat/s) (Asay-Davis et al. 2016 for a short description)
+         !                       ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description)
          rn_gammat0  = 0.0215    ! gammat coefficient used in blk formula
          rn_gammas0  = 0.614e-3  ! gammas coefficient used in blk formula
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/dtatsd.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/dtatsd.F90	(revision 12077)
@@ -69,8 +69,8 @@
       REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
       READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtsd in reference namelist' )
       REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run
       READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtsd in configuration namelist' )
       IF(lwm) WRITE ( numond, namtsd )
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/eosbn2.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/eosbn2.F90	(revision 12077)
@@ -29,5 +29,6 @@
    !!   eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass
    !!   eos_insitu_2d : Compute the in situ density for 2d fields
-   !!   bn2           : Compute the Brunt-Vaisala frequency
+   !!   bn2           : compute the Brunt-Vaisala frequency
+   !!   eos_pt_from_ct: compute the potential temperature from the Conservative Temperature
    !!   eos_rab       : generic interface of in situ thermal/haline expansion ratio 
    !!   eos_rab_3d    : compute in situ thermal/haline expansion ratio
@@ -74,7 +75,7 @@
 
    !                               !!** Namelist nameos **
-   LOGICAL , PUBLIC ::   ln_TEOS10   ! determine if eos_pt_from_ct is used to compute sst_m
-   LOGICAL , PUBLIC ::   ln_EOS80   ! determine if eos_pt_from_ct is used to compute sst_m
-   LOGICAL , PUBLIC ::   ln_SEOS   ! determine if eos_pt_from_ct is used to compute sst_m
+   LOGICAL , PUBLIC ::   ln_TEOS10
+   LOGICAL , PUBLIC ::   ln_EOS80
+   LOGICAL , PUBLIC ::   ln_SEOS
    LOGICAL , PUBLIC ::   ln_LEOS   ! determine if linear eos is used 
 
@@ -624,5 +625,5 @@
 
 
-   SUBROUTINE rab_3d( pts, pab )
+   SUBROUTINE rab_3d( pts, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_3d  ***
@@ -634,4 +635,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts   ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab   ! thermal/haline expansion ratio
@@ -652,5 +654,5 @@
                DO ji = 1, jpi
                   !
-                  zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth
+                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
@@ -710,5 +712,5 @@
                   zt  = pts (ji,jj,jk,jp_tem) - 10._wp   ! pot. temperature anomaly (t-T0)
                   zs  = pts (ji,jj,jk,jp_sal) - 35._wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)                ! depth in meters at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters at t-point
                   ztm = tmask(ji,jj,jk)                  ! land/sea bottom mask = surf. mask
                   !
@@ -730,5 +732,5 @@
                   zt  = pts (ji,jj,jk,jp_tem) - (-1._wp)
                   zs  = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)                 ! depth in meters at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)                 ! depth in meters at t-point
                   ztm = tmask(ji,jj,jk)                   ! land/sea bottom mask = surf. mask
                   !
@@ -757,5 +759,5 @@
 
 
-   SUBROUTINE rab_2d( pts, pdep, pab )
+   SUBROUTINE rab_2d( pts, pdep, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_2d  ***
@@ -765,4 +767,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj)         , INTENT(in   ) ::   pdep   ! depth                  [m]
@@ -891,5 +894,5 @@
 
 
-   SUBROUTINE rab_0d( pts, pdep, pab )
+   SUBROUTINE rab_0d( pts, pdep, pab, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE rab_0d  ***
@@ -899,4 +902,5 @@
       !! ** Action  : - pab     : thermal/haline expansion ratio at T-points
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpts)    , INTENT(in   ) ::   pts    ! pot. temperature & salinity
       REAL(wp),                      INTENT(in   ) ::   pdep   ! depth                  [m]
@@ -999,5 +1003,5 @@
 
 
-   SUBROUTINE bn2( pts, pab, pn2 )
+   SUBROUTINE bn2( pts, pab, pn2, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bn2  ***
@@ -1013,4 +1017,5 @@
       !!
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pts   ! pot. temperature and salinity   [Celsius,psu]
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::  pab   ! thermal/haline expansion coef.  [Celsius-1,psu-1]
@@ -1026,6 +1031,6 @@
          DO jj = 1, jpj          ! surface and bottom value set to zero one for all in istate.F90
             DO ji = 1, jpi
-               zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
-                  &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) ) 
+               zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+                  &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) ) 
                   !
                zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw 
@@ -1034,5 +1039,5 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
+                  &            / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
             END DO
          END DO
@@ -1203,5 +1208,5 @@
 
 
-   SUBROUTINE eos_pen( pts, pab_pe, ppen )
+   SUBROUTINE eos_pen( pts, pab_pe, ppen, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE eos_pen  ***
@@ -1223,4 +1228,5 @@
       !!                    pab_pe(:,:,:,jp_sal) is beta_pe
       !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   Kmm   ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts     ! pot. temperature & salinity
       REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pab_pe  ! alpha_pe and beta_pe
@@ -1242,5 +1248,5 @@
                DO ji = 1, jpi
                   !
-                  zh  = gdept_n(ji,jj,jk) * r1_Z0                                ! depth
+                  zh  = gdept(ji,jj,jk,Kmm) * r1_Z0                                ! depth
                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature
                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
@@ -1306,5 +1312,5 @@
                   zt  = pts(ji,jj,jk,jp_tem) - 10._wp  ! temperature anomaly (t-T0)
                   zs = pts (ji,jj,jk,jp_sal) - 35._wp  ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)              ! depth in meters  at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)              ! depth in meters  at t-point
                   ztm = tmask(ji,jj,jk)                ! tmask
                   zn  = 0.5_wp * zh * r1_rau0 * ztm
@@ -1326,5 +1332,5 @@
                   zt  = pts(ji,jj,jk,jp_tem) - (-1._wp)  ! temperature anomaly (t-T0)
                   zs = pts (ji,jj,jk,jp_sal) - 34.2_wp   ! abs. salinity anomaly (s-S0)
-                  zh  = gdept_n(ji,jj,jk)                ! depth in meters  at t-point
+                  zh  = gdept(ji,jj,jk,Kmm)                ! depth in meters  at t-point
                   ztm = tmask(ji,jj,jk)                  ! tmask
                   zn  = 0.5_wp * zh * r1_rau0 * ztm
@@ -1368,9 +1374,9 @@
       REWIND( numnam_ref )              ! Namelist nameos in reference namelist : equation of state
       READ  ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nameos in reference namelist' )
       !
       REWIND( numnam_cfg )              ! Namelist nameos in configuration namelist : equation of state
       READ  ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nameos in configuration namelist' )
       IF(lwm) WRITE( numond, nameos )
       !
@@ -1782,4 +1788,7 @@
          !
       CASE( np_seos )                        !==  Simplified EOS     ==!
+
+         r1_S0  = 0.875_wp/35.16504_wp   ! Used to convert CT in potential temperature when using bulk formulae (eos_pt_from_ct)
+         
          IF(lwp) THEN
             WRITE(numout,*)
Index: MO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isf.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isf.F90	(revision 12076)
+++ 	(revision )
@@ -1,264 +1,0 @@
-MODULE isf
-   !!======================================================================
-   !!                       ***  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 in_out_manager, ONLY: wp, jpi,jpj, jpk, 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_vtide        !: tidal background velocity (can be different to what is used in the 
-   REAL(wp)          , PUBLIC :: rn_htbl         !: Losch top boundary layer thickness [m]
-   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 = 0.0_wp          !: ISOMIP+ no 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
-   !-------------------------------------------------------
-   !
-   ! 0.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
-   !
-   ! 0.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
-   !
-   ! 0.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                      !:
-   !
-   ! 0.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: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
-   !! 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
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isf_oce.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 12077)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isf_oce.F90	(revision 12077)
@@ -0,0 +1,267 @@
+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_vtide        !: tidal background velocity (can be different to what is used in the 
+   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 = 0.0_wp          !: ISOMIP+ no 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: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
+   !! 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/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfcavgam.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfcavgam.F90	(revision 12077)
@@ -10,9 +10,9 @@
    !!   isfcav_gammats       : compute exchange coeficient gamma 
    !!----------------------------------------------------------------------
-   USE isf
+   USE isf_oce
    USE isfutils, ONLY: debug
    USE isftbl  , ONLY: isf_tbl
 
-   USE oce     , ONLY: un, vn, rn2         ! ocean dynamics and tracers
+   USE oce     , ONLY: uu, vv, rn2         ! ocean dynamics and tracers
    USE phycst  , ONLY: grav, vkarmn        ! physical constant
    USE eosbn2  , ONLY: eos_rab             ! equation of state
@@ -41,14 +41,15 @@
    !!-----------------------------------------------------------------------------------------------------
    !
-   SUBROUTINE isfcav_gammats( pttbl, pstbl, pqoce, pqfwf, pgt, pgs )
+   SUBROUTINE isfcav_gammats( Kmm, 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, AD15 and HJ99)
+      !!                 3 method available (cst, vel and vel_stab)
       !!---------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt  , pgs      ! gamma t and gamma s 
       !!-------------------------- IN  -------------------------------------
+      INTEGER                                     :: Kmm             ! ocean time level index
       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
@@ -57,13 +58,16 @@
       !!---------------------------------------------------------------------
       !
-      ! compute velocity in the tbl if needed
+      !==========================================
+      ! 1.: compute velocity in the tbl if needed
+      !==========================================
+      !
       SELECT CASE ( cn_gammablk )
       CASE ( 'spe'  ) 
          ! gamma is constant (specified in namelist)
          ! nothing to do
-      CASE ('ad15', 'hj99')
+      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)
+         CALL isf_tbl(Kmm, uu(:,:,:,Kmm) ,zutbl(:,:),'U', miku, rhisf_tbl_cav)
+         CALL isf_tbl(Kmm, vv(:,:,:,Kmm) ,zvtbl(:,:),'V', mikv, rhisf_tbl_cav)
          !
          ! mask velocity in tbl with ice shelf mask
@@ -78,18 +82,24 @@
       END SELECT
       ! 
-      ! compute gamma
+      !==========================================
+      ! 2.: compute gamma
+      !==========================================
+      !
       SELECT CASE ( cn_gammablk )
       CASE ( 'spe'  ) ! gamma is constant (specified in namelist)
          pgt(:,:) = rn_gammat0
          pgs(:,:) = rn_gammas0
-      CASE ( 'ad15' ) ! gamma is proportional to u*
-         CALL gammats_AD15 (              zutbl, zvtbl, rCd0_top, rn_vtide**2,               pgt, pgs )
-      CASE ( 'hj99' ) ! gamma depends of stability of boundary layer and u*
-         CALL gammats_HJ99 (pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top  , pqoce, pqfwf, pgt, pgs )
+      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 (Kmm, 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
       !
-      ! ouput exchange coeficient and tbl velocity
+      !==========================================
+      ! 3.: output and debug
+      !==========================================
+      !
       CALL iom_put('isfgammat', pgt(:,:))
       CALL iom_put('isfgammas', pgs(:,:))
@@ -106,5 +116,5 @@
    !!-----------------------------------------------------------------------------------------------------
    !
-   SUBROUTINE gammats_AD15(putbl, pvtbl, pCd, pke2, &   ! <<== in
+   SUBROUTINE gammats_vel( putbl, pvtbl, pCd, pke2, &   ! <<== in
       &                                  pgt, pgs   )   ! ==>> out gammats [m/s]
       !!----------------------------------------------------------------------
@@ -113,6 +123,5 @@
       !! ** Method     : gamma is velocity dependent ( gt= gt0 * Ustar )
       !!
-      !! ** Reference  : Jenkins et al., 2010, JPO, p2298-2312
-      !!                 Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
+      !! ** Reference  : Asay-Davis et al., Geosci. Model Dev., 9, 2471-2497, 2016
       !!---------------------------------------------------------------------
       !!-------------------------- OUT -------------------------------------
@@ -136,8 +145,8 @@
       CALL iom_put('isfustar',zustar(:,:))
       !
-   END SUBROUTINE gammats_AD15
-
-   SUBROUTINE gammats_HJ99( pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
-      &                                                            pgt  , pgs    )  ! ==>> out gammats [m/s]
+   END SUBROUTINE gammats_vel
+
+   SUBROUTINE gammats_vel_stab( Kmm, pttbl, pstbl, putbl, pvtbl, pCd, pke2, pqoce, pqfwf, &  ! <<== in
+      &                                                                     pgt  , pgs    )  ! ==>> out gammats [m/s]
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange coefficient 
@@ -150,4 +159,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: pgt, pgs     ! gammat and gammas
       !!-------------------------- IN  -------------------------------------
+      INTEGER                                     :: Kmm            ! ocean time level index
       REAL(wp),                     INTENT(in   ) :: pke2           ! background velocity squared
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pqoce, pqfwf   ! surface heat flux and fwf flux
@@ -199,10 +209,10 @@
                ! 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)
+               zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm)
                !                                            ! 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)  )
+               zdku = zcoef * (  uu(ji-1,jj  ,ikt  ,Kmm) + uu(ji,jj,ikt  ,Kmm)  &
+                  &             -uu(ji-1,jj  ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm)  )
+               zdkv = zcoef * (  vv(ji  ,jj-1,ikt  ,Kmm) + vv(ji,jj,ikt  ,Kmm)  &
+                  &             -vv(ji  ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm)  )
                !                                            ! richardson number (minimum value set to zero)
                zRc = MAX(rn2(ji,jj,ikt+1), 0._wp) / MAX( zdku*zdku + zdkv*zdkv, zeps )
@@ -211,7 +221,7 @@
                zts(jp_tem) = pttbl(ji,jj)
                zts(jp_sal) = pstbl(ji,jj)
-               zdep        = gdepw_n(ji,jj,ikt)
-               !
-               CALL eos_rab( zts, zdep, zab )
+               zdep        = gdepw(ji,jj,ikt,Kmm)
+               !
+               CALL eos_rab( zts, zdep, zab, Kmm )
                !
                ! compute length scale (Eq ??)
@@ -220,5 +230,5 @@
                ! 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
+               zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp
                !
                ! Compute Monin obukhov length scale at the surface and Ekman depth: (Eq ??)
@@ -243,5 +253,5 @@
       END DO
 
-   END SUBROUTINE gammats_HJ99
+   END SUBROUTINE gammats_vel_stab
 
 END MODULE isfcavgam
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfstp.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/isfstp.F90	(revision 12077)
@@ -14,5 +14,5 @@
    !!----------------------------------------------------------------------
    !
-   USE isf            ! isf variables
+   USE isf_oce                                      ! isf variables
    USE isfload, ONLY: isf_load                      ! ice shelf load
    USE isftbl , ONLY: isf_tbl_lvl                   ! ice shelf boundary layer
@@ -21,9 +21,10 @@
    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 dom_oce, ONLY: ht, e3t, 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
@@ -42,5 +43,5 @@
 CONTAINS
  
-  SUBROUTINE isf_stp( kt )
+  SUBROUTINE isf_stp( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isf_stp  ***
@@ -58,11 +59,16 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('isf')
       !
+      !=======================================================================
+      ! 1.: compute melt and associated heat fluxes in the ice shelf cavities
+      !=======================================================================
+      !
       IF ( ln_isfcav_mlt ) THEN
          !
-         ! before time step 
+         ! 1.1: before time step 
          IF ( kt /= nit000 ) THEN 
             risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
@@ -70,16 +76,20 @@
          END IF
          !
-         ! compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
+         ! 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)
-         !
-         ! compute ice shelf melt
-         CALL isf_cav( kt, risf_cav_tsc, fwfisf_cav)
+         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav)
+         !
+         ! 1.3: compute ice shelf melt
+         CALL isf_cav( kt, Kmm, 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
          !
-         ! before time step 
+         ! 2.1: before time step 
          IF ( kt /= nit000 ) THEN 
             risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
@@ -87,20 +97,19 @@
          END IF
          !
-         ! 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
+         ! 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)
-         !
-         ! compute ice shelf melt
-         CALL isf_par( kt, risf_par_tsc, fwfisf_par)
-         !
-      END IF
-      !
-      IF ( ln_isfcpl ) THEN
-         ! after step nit000 + 2 we do not need anymore the risfcpl_ arrays
-         IF ( kt == nit000 + 2 ) CALL isf_dealloc_cpl()
-
-         IF ( lrst_oce ) CALL isfcpl_rst_write(kt)
-      END IF
+         CALL isf_tbl_lvl(ht, e3t(:,:,:,Kmm), misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par)
+         !
+         ! 2.3: compute ice shelf melt
+         CALL isf_par( kt, Kmm, 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, Kmm)
       !
       IF( ln_timing )   CALL timing_stop('isf')
@@ -108,5 +117,5 @@
    END SUBROUTINE isf_stp
 
-   SUBROUTINE isf_init
+   SUBROUTINE isf_init(Kbb, Kmm, Kaa)
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE isfstp_init  ***
@@ -122,4 +131,5 @@
       !!              - call cav/param/isfcpl init routine
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa      ! ocean time level indices
       !
       ! constrain: l_isfoasis need to be known
@@ -135,5 +145,5 @@
       !
       ! compute ice shelf load
-      IF ( ln_isfcav ) CALL isf_load( risfload )
+      IF ( ln_isfcav ) CALL isf_load( Kmm, risfload )
       !
       ! terminate routine now if no ice shelf melt formulation specify
@@ -150,5 +160,5 @@
          !---------------------------------------------------------------------------------------------------------------------
          ! initialisation ice sheet coupling
-         IF ( ln_isfcpl ) CALL isfcpl_init()
+         IF( ln_isfcpl ) CALL isfcpl_init(Kbb, Kmm, Kaa)
          !
       END IF
@@ -177,12 +187,12 @@
             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) 
+               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_vtide**2 = ', rn_vtide**2
-                  WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0      = ', r_Cdmin_top
+                  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_vtide**2  = ', rn_vtide**2
+                  WRITE(numout,*) '         top drag coef.    used (from namdrg_top) rn_Cd0       = ', r_Cdmin_top
                END IF
             END IF
@@ -212,6 +222,11 @@
          END IF
 
-         IF (ln_isfcav) WRITE(numout,*) '      Ice shelf load method                   cn_isfload        = ', TRIM(cn_isfload)
+         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
@@ -267,18 +282,20 @@
       INTEGER               :: ios                  ! Local integer output status for namelist read
       !!----------------------------------------------------------------------
-      NAMELIST/namisf/ ln_isf       ,                                                                               & 
-         &             ln_isfcav_mlt, cn_isfcav_mlt, cn_gammablk, rn_gammat0, rn_gammas0, rn_htbl, 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, rn_vtide,                              &
-         &             cn_isfload   , cn_isfdir
+      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, rn_vtide,    &
+         &             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', lwp )
+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', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namisf in configuration namelist' )
       IF(lwm) WRITE ( numond, namisf )
 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/istate.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/istate.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/istate.F90	(revision 12077)
@@ -50,5 +50,5 @@
 CONTAINS
 
-   SUBROUTINE istate_init
+   SUBROUTINE istate_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
@@ -56,4 +56,6 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
+      !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
 !!gm see comment further down
@@ -75,14 +77,14 @@
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
+      ts   (:,:,:,:,Kaa) = 0._wp                               ! set one for all to 0 at level jpk
       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
-      ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization
-      va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization    
+      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
+      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization    
 #endif
 
       IF( ln_rstart ) THEN                    ! Restart from a file
          !                                    ! -------------------
-         CALL rst_read                        ! Read the restart file
+         CALL rst_read( Kbb, Kmm )            ! Read the restart file
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !
@@ -95,9 +97,9 @@
          !
          IF( ln_tsd_init ) THEN               
-            CALL dta_tsd( nit000, 'ini', tsb )       ! read 3D T and S data at nit000
+            CALL dta_tsd( nit000, 'ini', ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
             !
-            sshb(:,:)   = 0._wp               ! set the ocean at rest
+            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
             IF( ll_wd ) THEN
-               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
+               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
                !
                ! Apply minimum wetdepth criterion
@@ -105,26 +107,26 @@
                DO jj = 1,jpj
                   DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
-                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
                      ENDIF
                   END DO
                END DO 
             ENDIF 
-            ub  (:,:,:) = 0._wp
-            vb  (:,:,:) = 0._wp  
+            uu  (:,:,:,Kbb) = 0._wp
+            vv  (:,:,:,Kbb) = 0._wp  
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
+            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
-         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-         sshn (:,:)     = sshb(:,:)   
-         un   (:,:,:)   = ub  (:,:,:)
-         vn   (:,:,:)   = vb  (:,:,:)
-         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-         CALL div_hor( 0 )                    ! compute interior hdivn value  
-!!gm                                    hdivn(:,:,:) = 0._wp
+         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
+         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
+         CALL div_hor( 0, Kbb, Kmm )         ! compute interior hdiv value  
+!!gm                                    hdiv(:,:,:) = 0._wp
 
 !!gm POTENTIAL BUG :
-!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
+!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
 !!             as well as gdept and gdepw....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
@@ -136,14 +138,14 @@
 !            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
 !            CALL dta_uvd( nit000, zuvd )
-!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
+!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
 !            DEALLOCATE( zuvd )
 !         ENDIF
          !
 !!gm This is to be changed !!!!
-!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
 !         IF( .NOT.ln_linssh ) THEN
 !            DO jk = 1, jpk
-!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
 !            END DO
 !         ENDIF
@@ -155,25 +157,25 @@
       ! Do it whatever the free surface method, these arrays being eventually used
       !
-      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
+      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp
+      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp
       !
-!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
+!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
-      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
       !
-      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
+      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
+      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
       !
    END SUBROUTINE istate_init
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/sbcfwb.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/sbcfwb.F90	(revision 12077)
@@ -17,10 +17,9 @@
    USE dom_oce        ! ocean space and time domain
    USE sbc_oce        ! surface ocean boundary condition
+   USE isf_oce        ! ice shelf melting contribution
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b, snwice_fmass
    USE phycst         ! physical constants
    USE sbcrnf         ! ocean runoffs
-   USE isf            ! ice shelf melting contribution
    USE sbcssr         ! Sea-Surface damping terms
-   USE tradmp         !
    !
    USE in_out_manager ! I/O manager
@@ -49,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
+   SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_fwb  ***
@@ -66,4 +65,5 @@
       INTEGER, INTENT( in ) ::   kn_fsbc  ! 
       INTEGER, INTENT( in ) ::   kn_fwb   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm      ! ocean time level index
       !
       INTEGER  ::   inum, ikty, iyear     ! local integers
@@ -157,5 +157,5 @@
             a_fwb_b = a_fwb                           ! mean sea level taking into account the ice+snow
                                                       ! sum over the global domain
-            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( sshn(:,:) + snwice_mass(:,:) * r1_rau0 ) )
+            a_fwb   = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rau0 ) )
             a_fwb   = a_fwb * 1.e+3 / ( area * rday * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/tradmp.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/tradmp.F90	(revision 12077)
@@ -72,5 +72,5 @@
 
 
-   SUBROUTINE tra_dmp( kt )
+   SUBROUTINE tra_dmp( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_dmp  ***
@@ -90,5 +90,7 @@
       !! ** Action  : - tsa: tracer trends updated with the damping trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
@@ -101,5 +103,5 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) 
+         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -113,5 +115,6 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - tsb(ji,jj,jk,jn) )
+                     pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
+                        &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
                   END DO
                END DO
@@ -124,8 +127,8 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -137,9 +140,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) )
-                     tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
+                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -150,12 +153,12 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = tsa(:,:,:,:) - ztrdts(:,:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
+         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
          DEALLOCATE( ztrdts ) 
       ENDIF
       !                           ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_dmp')
@@ -179,9 +182,9 @@
       REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation
       READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp )
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_dmp in reference namelist' )
       !
       REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation
       READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist' )
       IF(lwm) WRITE ( numond, namtra_dmp )
       !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/usrdef_sbc.F90	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP+/MY_SRC/usrdef_sbc.F90	(revision 12077)
@@ -41,5 +41,5 @@
 CONTAINS
 
-   SUBROUTINE usrdef_sbc_oce( kt )
+   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE usr_def_sbc  ***
@@ -56,4 +56,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
       !!---------------------------------------------------------------------
       !
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP/EXPREF/namelist_cfg	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/ISOMIP/EXPREF/namelist_cfg	(revision 12077)
@@ -150,5 +150,4 @@
    ! ---------------- ice shelf load -------------------------------
    !
-   cn_isfload = 'isomip'      ! scheme to compute ice shelf load (ln_isfcav = .true. in domain_cfg.nc)
    !
    ! ---------------- ice shelf melt formulation -------------------------------
Index: /NEMO/branches/2019/UKMO_MERGE_2019/tests/demo_cfgs.txt
===================================================================
--- /NEMO/branches/2019/UKMO_MERGE_2019/tests/demo_cfgs.txt	(revision 12076)
+++ /NEMO/branches/2019/UKMO_MERGE_2019/tests/demo_cfgs.txt	(revision 12077)
@@ -1,4 +1,5 @@
 CANAL OCE
 ISOMIP OCE
+ISOMIP+ OCE
 LOCK_EXCHANGE OCE
 OVERFLOW OCE
