Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/pdControl/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/pdControl/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/pdControl/config.card	(revision 3136)
@@ -40,9 +40,10 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
+
 #========================================================================
 #D-- ATM -
@@ -52,7 +53,7 @@
 Restart= y
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -63,9 +64,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -76,9 +77,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -88,9 +89,9 @@
 [MBG]
 WriteFrequency="1M 1Y"
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -103,9 +104,9 @@
 Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -116,9 +117,9 @@
 Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -127,11 +128,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/piControl/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/piControl/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/IPSLCM/piControl/config.card	(revision 3136)
@@ -40,9 +40,10 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
+
 #========================================================================
 #D-- ATM -
@@ -52,7 +53,7 @@
 Restart= y
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -63,9 +64,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -76,9 +77,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -88,9 +89,9 @@
 [MBG]
 WriteFrequency="1M 1Y"
-Restart= n
+Restart= y
 ##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=1949-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM605-LR-pdCtrl-Restart
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -103,9 +104,9 @@
 Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -116,9 +117,9 @@
 Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -127,11 +128,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= n
+Restart= y
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=1859-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=O1T09V06
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=/ccc/store/cont003/dsm/p86mart/IGCM_OUT/IPSLCM6/DEVT/piControl
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZ/clim_360d/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZ/clim_360d/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZ/clim_360d/config.card	(revision 3136)
@@ -52,9 +52,9 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=2000-01-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=EXP00
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZ/DEVT/amip
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -63,9 +63,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= y
+Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/amip/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/amip/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/amip/config.card	(revision 3136)
@@ -56,9 +56,9 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=2000-01-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=EXP00
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -67,9 +67,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= y
+Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -80,10 +80,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -92,10 +93,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-#-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_360d/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_360d/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_360d/config.card	(revision 3136)
@@ -56,9 +56,9 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=2000-01-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=EXP00
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -67,9 +67,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= y
+Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -80,10 +80,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -92,10 +93,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-#-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
Index: CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_pdControl/config.card
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_pdControl/config.card	(revision 3134)
+++ CONFIG/UNIFORM/v6/IPSLCM6/EXPERIMENTS/LMDZOR/clim_pdControl/config.card	(revision 3136)
@@ -56,9 +56,9 @@
 OverRule=n
 #D- Last day of the experience used as restart for all components
-RestartDate=2000-01-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=EXP00
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -67,9 +67,9 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart= y
+Restart= n
 # Last day of the experience used as restart for this component if Restart=y
-RestartDate=2029-12-31
+RestartDate=2319-12-31
 #D- Define restart simulation name for all components
-RestartJobName=CM606-LR-pdCtrl-01
+RestartJobName=CM607-LR-pdCtrl-01
 #D- Path Server Group Login
 RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
@@ -80,10 +80,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-##-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
@@ -92,10 +93,11 @@
 WriteFrequency="1M"
 # If config_Restarts_OverRule == 'n' next 4 params are read
-Restart=n
-#-- Last day of the experience used as restart for this component if Restart=y
-RestartDate=2000-01-31
-# Define restart simulation name for this component
-RestartJobName=EXP00
-RestartPath=${ARCHIVE}/IGCM_OUT/LMDZOR/DEVT/amip
+Restart= n
+# Last day of the experience used as restart for this component if Restart=y
+RestartDate=2319-12-31
+#D- Define restart simulation name for all components
+RestartJobName=CM607-LR-pdCtrl-01
+#D- Path Server Group Login
+RestartPath=${R_IN}/RESTART/IPSLCM6/DEVT/pdControl
 
 #========================================================================
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcisf.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcisf.F90	(revision 3136)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcisf.F90	(revision 3136)
@@ -0,0 +1,961 @@
+MODULE sbcisf
+   !!======================================================================
+   !!                       ***  MODULE  sbcisf  ***
+   !! Surface module :  update surface ocean boundary condition under ice
+   !!                   shelf
+   !!======================================================================
+   !! History :  3.2   !  2011-02  (C.Harris  ) Original code isf cav
+   !!            X.X   !  2006-02  (C. Wang   ) Original code bg03
+   !!            3.4   !  2013-03  (P. Mathiot) Merging + parametrization
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_isf        : update sbc under ice shelf
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE eosbn2          ! equation of state
+   USE sbc_oce         ! surface boundary condition: ocean fields
+   USE lbclnk          !
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE wrk_nemo        ! Memory allocation
+   USE timing          ! Timing
+   USE lib_fortran     ! glob_sum
+   USE zdfbfr
+   USE fldread         ! read input field at current time step
+
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
+
+   ! public in order to be able to output then 
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf
+   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
+   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence 
+   INTEGER , PUBLIC ::   nn_isfblk                   !: 
+   INTEGER , PUBLIC ::   nn_gammablk                 !:
+   LOGICAL , PUBLIC ::   ln_conserve                 !:
+   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
+   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient 
+   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
+
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl 
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
+   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
+
+   LOGICAL, PUBLIC ::   l_isfcpl = .false.       ! isf recieved from oasis
+
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
+   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
+
+!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
+   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
+   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
+   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
+   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
+   
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
+   !! $Id: sbcisf.F90 7607 2017-01-25 15:37:31Z cetlod $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+ 
+  SUBROUTINE sbc_isf(kt)
+
+    INTEGER, INTENT(in)          ::   kt         ! ocean time step
+    INTEGER                      ::   ji, jj, jk
+    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)                     ::   zhk
+    REAL(wp)                     ::   zt_frz, zpress
+    REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
+    REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
+    REAL(wp)                            :: zhisf
+
+
+      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
+
+         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkt(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt, mbkt(ji,jj)
+                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
+               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
+
+               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
+               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+            END DO
+         END DO
+
+         ! compute salf and heat flux
+         IF (nn_isf == 1) THEN
+            ! realistic ice shelf formulation
+            ! compute T/S/U/V for the top boundary layer
+            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
+            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
+            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
+            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
+            ! iom print
+            CALL iom_put('ttbl',ttbl(:,:))
+            CALL iom_put('stbl',stbl(:,:))
+            CALL iom_put('utbl',utbl(:,:))
+            CALL iom_put('vtbl',vtbl(:,:))
+            ! compute fwf and heat flux
+            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
+            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * lfusisf              ! heat        flux
+            ENDIF
+
+         ELSE IF (nn_isf == 2) THEN
+            ! Beckmann and Goosse parametrisation 
+            stbl(:,:)   = soce
+            CALL sbc_isf_bg03(kt)
+
+         ELSE IF (nn_isf == 3) THEN
+            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
+            IF( .NOT.l_isfcpl ) THEN
+               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
+               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+            ENDIF
+            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            stbl(:,:)   = soce
+
+         ELSE IF (nn_isf == 4) THEN
+            ! specified fwf and heat flux forcing beneath the ice shelf
+            IF( .NOT.l_isfcpl ) THEN
+               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
+               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
+               fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+            ENDIF
+            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
+            stbl(:,:)   = soce
+
+         END IF
+         ! compute tsc due to isf
+         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
+!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
+         
+         ! salt effect already take into account in vertical advection
+         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
+
+         ! output
+         IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf)
+         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce )
+
+         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
+         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
+ 
+         ! lbclnk
+         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
+         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
+         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
+
+         ! Diagnostics
+         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
+            !
+            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
+            !
+            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
+            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
+            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
+            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
+            !
+            DO jj = 1,jpj
+               DO ji = 1,jpi
+                  ikt = misfkt(ji,jj)
+                  ikb = misfkb(ji,jj)
+                  DO jk = ikt, ikb - 1
+                     zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj)    * zhisf
+                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf
+                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj)      * zhisf
+                  END DO
+                  jk = ikb
+                  zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
+                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * zhisf * ralpha(ji,jj) 
+                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj)
+                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * zhisf * ralpha(ji,jj)
+               END DO
+            END DO
+            !
+            CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) )
+            CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) )
+            CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) )
+            CALL iom_put( 'qhcisf'   , zqhcisf2d (:,:  ) )
+            !
+            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
+            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
+            !
+         END IF
+         ! 
+      END IF
+      !
+      !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
+              & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend
+            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
+            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
+         ELSE
+            fwfisf_b(:,:)    = fwfisf(:,:)
+            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
+         END IF
+      ENDIF
+      !
+      IF( lrst_oce ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) )
+         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
+         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
+      ENDIF
+       !
+  END SUBROUTINE sbc_isf
+
+  SUBROUTINE sbc_isf_init
+
+    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
+    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)                     ::   zhk
+    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
+    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
+    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
+    INTEGER           ::   ios           ! Local integer output status for namelist read
+
+      !
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
+                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
+      !
+      !
+         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs 
+         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
+         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namsbc_isf )
+
+
+         IF ( lwp ) WRITE(numout,*)
+         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
+         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
+         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
+         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
+         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
+         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
+         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
+         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
+         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
+         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
+         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
+            rdivisf = 1._wp
+         ELSE
+            rdivisf = 0._wp
+         END IF
+         !
+         ! Allocate public variable
+         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
+         !
+         ! initialisation
+         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
+         risf_tsc(:,:,:)  = 0._wp
+         !
+         ! define isf tbl tickness, top and bottom indice
+         IF      (nn_isf == 1) THEN
+            rhisf_tbl(:,:) = rn_hisf_tbl
+            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
+            IF( .NOT.l_isfcpl ) THEN
+               ALLOCATE( sf_rnfisf(1), STAT=ierror )
+               ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+             ENDIF
+
+            !: read effective lenght (BG03)
+            IF (nn_isf == 2) THEN
+               ! Read Data and save some integral values
+               CALL iom_open( sn_Leff_isf%clname, inum )
+               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
+               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
+               CALL iom_close(inum)
+               !
+               risfLeff = risfLeff*1000           !: convertion in m
+            END IF
+
+           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
+            CALL iom_open( sn_depmax_isf%clname, inum )
+            cvarhisf = TRIM(sn_depmax_isf%clvar)
+            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
+            CALL iom_close(inum)
+            !
+            CALL iom_open( sn_depmin_isf%clname, inum )
+            cvarzisf = TRIM(sn_depmin_isf%clvar)
+            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
+            CALL iom_close(inum)
+            !
+            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
+
+           !! compute first level of the top boundary layer
+           DO ji = 1, jpi
+              DO jj = 1, jpj
+                  jk = 2
+                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
+                  misfkt(ji,jj) = jk-1
+               END DO
+            END DO
+
+         ELSE IF ( nn_isf == 4 ) THEN
+            ! as in nn_isf == 1
+            rhisf_tbl(:,:) = rn_hisf_tbl
+            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
+            
+            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
+            IF( .NOT.l_isfcpl ) THEN
+               ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
+               ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
+               ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
+               CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
+               !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
+            ENDIF
+         END IF
+         ! save initial top boundary layer thickness         
+         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
+         ! 
+   END SUBROUTINE sbc_isf_init
+      
+
+
+  INTEGER FUNCTION sbc_isf_alloc()
+      !!----------------------------------------------------------------------
+      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
+      !!----------------------------------------------------------------------
+      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
+      IF( .NOT. ALLOCATED( qisf ) ) THEN
+         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
+               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
+               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
+               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
+               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
+               &    STAT= sbc_isf_alloc )
+         !
+         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
+         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
+         !
+      ENDIF
+  END FUNCTION
+
+  SUBROUTINE sbc_isf_bg03(kt)
+   !!==========================================================================
+   !!                 *** SUBROUTINE sbcisf_bg03  ***
+   !! add net heat and fresh water flux from ice shelf melting
+   !! into the adjacent ocean using the parameterisation by
+   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
+   !!  (hereafter BG)
+   !!==========================================================================
+   !!----------------------------------------------------------------------
+   !!   sbc_isf_bg03      : routine called from sbcmod
+   !!----------------------------------------------------------------------
+   !!
+   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
+   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
+   !!
+   !! History :
+   !!      !  06-02  (C. Wang) Original code
+   !!----------------------------------------------------------------------
+
+    INTEGER, INTENT ( in ) :: kt
+
+    INTEGER :: ji, jj, jk, jish  !temporary integer
+    INTEGER :: ijkmin
+    INTEGER :: ii, ij, ik 
+    INTEGER :: inum
+
+    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
+    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
+    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
+    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
+    
+    !!----------------------------------------------------------------------
+    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
+     !
+
+    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
+    DO ji = 1, jpi
+       DO jj = 1, jpj
+          ik = misfkt(ji,jj)
+          !! Initialize arrays to 0 (each step)
+          zt_sum = 0.e0_wp
+          IF ( ik .GT. 1 ) THEN
+    ! 3. -----------the average temperature between 200m and 600m ---------------------
+             DO jk = misfkt(ji,jj),misfkb(ji,jj)
+             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
+             ! after verif with UNESCO, wrong sign in BG eq. 2
+             ! Calculate freezing temperature
+                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
+                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
+                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
+             ENDDO
+             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
+    
+    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
+          ! For those corresponding to zonal boundary    
+             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
+                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
+             
+             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                  
+             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
+             !add to salinity trend
+          ELSE
+             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
+          END IF
+       ENDDO
+    ENDDO
+    !
+    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+  END SUBROUTINE sbc_isf_bg03
+
+   SUBROUTINE sbc_isf_cav( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_isf_cav  ***
+      !!
+      !! ** Purpose :   handle surface boundary condition under ice shelf
+      !!
+      !! ** Method  : -
+      !!
+      !! ** Action  :   utau, vtau : remain unchanged
+      !!                taum, wndm : remain unchanged
+      !!                qns        : update heat flux below ice shelf
+      !!                emp, emps  : update freshwater flux below ice shelf
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in)          ::   kt         ! ocean time step
+      !
+      LOGICAL :: ln_isomip = .true.
+      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
+      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
+      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
+      REAL(wp) ::   zlamb1, zlamb2, zlamb3
+      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
+      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
+      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
+      REAL(wp) ::   zgammat, zgammas
+      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
+      INTEGER  ::   ierror     ! return error code
+      LOGICAL  ::   lit=.TRUE.
+      INTEGER  ::   nit
+      !!---------------------------------------------------------------------
+      !
+      ! coeficient for linearisation of tfreez
+      zlamb1=-0.0575
+      zlamb2=0.0901
+      zlamb3=-7.61e-04
+      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
+      !
+      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
+
+      zcfac=0.0_wp 
+      IF (ln_conserve)  zcfac=1.0_wp
+      zpress(:,:)=0.0_wp
+      zgammat2d(:,:)=0.0_wp
+      zgammas2d(:,:)=0.0_wp
+      !
+      !
+!CDIR COLLAPSE
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! Crude approximation for pressure (but commonly used)
+            ! 1e-04 to convert from Pa to dBar
+            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
+            !
+         END DO
+      END DO
+
+! Calculate in-situ temperature (ref to surface)
+      zti(:,:)=tinsitu( ttbl, stbl, zpress )
+! Calculate freezing temperature
+      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
+
+      
+      zhtflx=0._wp ; zfwflx=0._wp
+      IF (nn_isfblk == 1) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF (mikt(ji,jj) > 1 ) THEN
+                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
+                  DO WHILE ( lit )
+! compute gamma
+                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
+! zhtflx is upward heat flux (out of ocean)
+                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
+! zwflx is upward water flux
+                     zfwflx = - zhtflx/lfusisf
+! test convergence and compute gammat
+                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
+
+                     nit = nit + 1
+                     IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+
+! save gammat and compute zhtflx_b
+                     zgammat2d(ji,jj)=zgammat
+                     zhtflx_b = zhtflx
+                  END DO
+
+                  qisf(ji,jj) = - zhtflx
+! For genuine ISOMIP protocol this should probably be something like
+                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
+               ELSE
+                  fwfisf(ji,jj) = 0._wp
+                  qisf(ji,jj)   = 0._wp
+               END IF
+            !
+            END DO
+         END DO
+
+      ELSE IF (nn_isfblk == 2 ) THEN
+
+! More complicated 3 equation thermodynamics as in MITgcm
+!CDIR COLLAPSE
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+               IF (mikt(ji,jj) > 1 ) THEN
+                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
+                  DO WHILE ( lit )
+                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
+
+                     zeps1=rcp*rau0*zgammat
+                     zeps2=lfusisf*rau0*zgammas
+                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
+                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
+                     zeps6=zeps4-zti(ji,jj)
+                     zeps7=zeps4-tsurf
+                     zaqe=zlamb1 * (zeps1 + zeps3)
+                     zaqer=0.5/zaqe
+                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
+                     zcqe=zeps2*stbl(ji,jj)
+                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
+! Presumably zdis can never be negative because gammas is very small compared to gammat
+                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
+                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
+                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
+  
+! zfwflx is upward water flux
+                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
+! zhtflx is upward heat flux (out of ocean)
+! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
+                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
+! zwflx is upward water flux
+! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
+                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
+! test convergence and compute gammat
+                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
+
+                     nit = nit + 1
+                     IF (nit .GE. 51) THEN
+                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
+                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
+                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+                     END IF
+! save gammat and compute zhtflx_b
+                     zgammat2d(ji,jj)=zgammat
+                     zgammas2d(ji,jj)=zgammas
+                     zhtflx_b = zhtflx
+
+                  END DO
+! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
+                  qisf(ji,jj) = - zhtflx 
+! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
+                  fwfisf(ji,jj) = zfwflx 
+               ELSE
+                  fwfisf(ji,jj) = 0._wp
+                  qisf(ji,jj)   = 0._wp
+               ENDIF
+               !
+            END DO
+         END DO
+      ENDIF
+      ! lbclnk
+      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
+      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
+      ! output
+      CALL iom_put('isfgammat', zgammat2d)
+      CALL iom_put('isfgammas', zgammas2d)
+      !
+      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
+
+   END SUBROUTINE sbc_isf_cav
+
+   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
+      !!----------------------------------------------------------------------
+      !! ** Purpose    : compute the coefficient echange for heat flux
+      !!
+      !! ** Method     : gamma assume constant or depends of u* and stability
+      !!
+      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
+      !!                Jenkins et al., 2010, JPO, p2298-2312
+      !!---------------------------------------------------------------------
+      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
+      INTEGER , INTENT(in)    :: ji,jj
+      LOGICAL , INTENT(inout) :: lit
+
+      INTEGER  :: ikt                 ! loop index
+      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
+      REAL(wp) :: zdku, zdkv                 ! U, V shear 
+      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number 
+      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
+      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
+      REAL(wp) :: zhmax                      ! limitation of mol
+      REAL(wp) :: zetastar                   ! stability parameter
+      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence 
+      REAL(wp) :: zcoef                      ! temporary coef
+      REAL(wp) :: zdep
+      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
+      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
+      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
+      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
+      REAL(wp), DIMENSION(2) :: zts, zab
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_gammablk == 0 ) THEN
+      !! gamma is constant (specified in namelist)
+         gt = rn_gammat0
+         gs = rn_gammas0
+         lit = .FALSE.
+      ELSE IF ( nn_gammablk == 1 ) THEN
+      !! gamma is assume to be proportional to u* 
+      !! WARNING in case of Losh 2008 tbl parametrization, 
+      !! you have to used the mean value of u in the boundary layer) 
+      !! not yet coded
+      !! Jenkins et al., 2010, JPO, p2298-2312
+         ikt = mikt(ji,jj)
+      !! Compute U and V at T points
+   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
+   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
+          zut = utbl(ji,jj)
+          zvt = vtbl(ji,jj)
+
+      !! compute ustar
+         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
+      !! Compute mean value over the TBL
+
+      !! Compute gammats
+         gt = zustar * rn_gammat0
+         gs = zustar * rn_gammas0
+         lit = .FALSE.
+      ELSE IF ( nn_gammablk == 2 ) THEN
+      !! gamma depends of stability of boundary layer
+      !! WARNING in case of Losh 2008 tbl parametrization, 
+      !! you have to used the mean value of u in the boundary layer) 
+      !! not yet coded
+      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
+      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
+               ikt = mikt(ji,jj)
+
+      !! Compute U and V at T points
+               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
+               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
+
+      !! compute ustar
+               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
+               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
+                 gt = rn_gammat0
+                 gs = rn_gammas0
+               ELSE
+      !! compute Rc number (as done in zdfric.F90)
+               zcoef = 0.5 / fse3w(ji,jj,ikt)
+               !                                            ! shear of horizontal velocity
+               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
+                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
+               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
+                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
+               !                                            ! richardson number (minimum value set to zero)
+               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
+
+      !! compute Pr and Sc number (can be improved)
+               zPr =   13.8
+               zSc = 2432.0
+
+      !! compute gamma mole
+               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
+               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
+
+      !! compute bouyancy 
+               zts(jp_tem) = ttbl(ji,jj)
+               zts(jp_sal) = stbl(ji,jj)
+               zdep        = fsdepw(ji,jj,ikt)
+               !
+               CALL eos_rab( zts, zdep, zab )
+                  !
+      !! compute length scale 
+               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      !! compute Monin Obukov Length
+               ! Maximum boundary layer depth
+               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
+               ! Compute Monin obukhov length scale at the surface and Ekman depth:
+               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
+               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+
+      !! compute eta* (stability parameter)
+               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
+
+      !! compute the sublayer thickness
+               zhnu = 5 * znu / zustar
+      !! compute gamma turb
+               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
+               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
+
+      !! compute gammats
+               gt = zustar / (zgturb + zgmolet)
+               gs = zustar / (zgturb + zgmoles)
+               END IF
+      END IF
+
+   END SUBROUTINE
+
+   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
+      !!
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer 
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
+      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
+      
+      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
+
+      REAL(wp) :: ze3, zhk
+      REAL(wp), DIMENSION(:,:), POINTER :: zikt
+
+      INTEGER :: ji,jj,jk
+      INTEGER :: ikt,ikb
+      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
+
+      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
+      CALL wrk_alloc( jpi,jpj, zikt )
+
+      ! get first and last level of tbl
+      mkt(:,:) = misfkt(:,:)
+      mkb(:,:) = misfkb(:,:)
+
+      varout(:,:)=0._wp
+      DO jj = 2,jpj
+         DO ji = 2,jpi
+            IF (ssmask(ji,jj) == 1) THEN
+               ikt = mkt(ji,jj)
+               ikb = mkb(ji,jj)
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3t_n(ji,jj,jk)
+                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
+                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
+                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
+                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
+                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
+               IF (cptin == 'T') &
+                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
+               IF (cptin == 'U') &
+                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
+               IF (cptin == 'V') &
+                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
+            END IF
+         END DO
+      END DO
+
+      CALL wrk_dealloc( jpi,jpj, mkt, mkb )      
+      CALL wrk_dealloc( jpi,jpj, zikt ) 
+
+      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
+      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
+
+   END SUBROUTINE sbc_isf_tbl
+      
+
+   SUBROUTINE sbc_isf_div( phdivn )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE sbc_isf_div  ***
+      !!       
+      !! ** Purpose :   update the horizontal divergence with the runoff inflow
+      !!
+      !! ** Method  :   
+      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the 
+      !!                          divergence and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   decreased by the runoff inflow
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   ikt, ikb 
+      INTEGER  ::   nk_isf
+      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
+      REAL(wp)     ::   zfact     ! local scalar
+      !!----------------------------------------------------------------------
+      !
+      zfact   = 0.5_wp
+      !
+      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkt(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               ! test on tmask useless ?????
+               DO jk = ikt, mbkt(ji,jj)
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
+               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
+
+               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
+               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
+            END DO
+         END DO
+      END IF ! vvl case
+      !
+      DO jj = 1,jpj
+         DO ji = 1,jpi
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
+                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
+               END DO
+               ! level partially include in ice shelf boundary layer 
+               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
+                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
+            !==   ice shelf melting mass distributed over several levels   ==!
+         END DO
+      END DO
+      !
+   END SUBROUTINE sbc_isf_div
+                        
+   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE eos_init  ***
+      !!
+      !! ** Purpose :   Compute the in-situ temperature [Celcius]
+      !!
+      !! ** Method  :   
+      !!
+      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
+      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
+!      REAL(wp) :: fsatg
+!      REAL(wp) :: pfps, pfpt, pfphp 
+      REAL(wp) :: zt, zs, zp, zh, zq, zxk
+      INTEGER  :: ji, jj            ! dummy loop indices
+      !
+      CALL wrk_alloc( jpi,jpj, pti  )
+      ! 
+      DO jj=1,jpj
+         DO ji=1,jpi
+            zh = ppress(ji,jj)
+! Theta1
+            zt = ptem(ji,jj)
+            zs = psal(ji,jj)
+            zp = 0.0
+            zxk= zh * fsatg( zs, zt, zp )
+            zt = zt + 0.5 * zxk
+            zq = zxk
+! Theta2
+            zp = zp + 0.5 * zh
+            zxk= zh*fsatg( zs, zt, zp )
+            zt = zt + 0.29289322 * ( zxk - zq )
+            zq = 0.58578644 * zxk + 0.121320344 * zq
+! Theta3
+            zxk= zh * fsatg( zs, zt, zp )
+            zt = zt + 1.707106781 * ( zxk - zq )
+            zq = 3.414213562 * zxk - 4.121320344 * zq
+! Theta4
+            zp = zp + 0.5 * zh
+            zxk= zh * fsatg( zs, zt, zp )
+            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi,jpj, pti  )
+      !
+   END FUNCTION tinsitu
+   !
+   FUNCTION fsatg( pfps, pfpt, pfphp )
+      !!----------------------------------------------------------------------
+      !!                 ***  FUNCTION fsatg  ***
+      !!
+      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
+      !!
+      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
+      !! 
+      !! ** units      :   pressure        pfphp    decibars
+      !!                   temperature     pfpt     deg celsius (ipts-68)
+      !!                   salinity        pfps     (ipss-78)
+      !!                   adiabatic       fsatg    deg. c/decibar
+      !!----------------------------------------------------------------------
+      REAL(wp) :: pfps, pfpt, pfphp 
+      REAL(wp) :: fsatg
+      !
+      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
+        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
+        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
+        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
+        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
+      !
+    END FUNCTION fsatg
+    !!======================================================================
+END MODULE sbcisf
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90	(revision 3136)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbcmod.F90	(revision 3136)
@@ -0,0 +1,517 @@
+MODULE sbcmod
+   !!======================================================================
+   !!                       ***  MODULE  sbcmod  ***
+   !! Surface module :  provide to the ocean its surface boundary condition
+   !!======================================================================
+   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
+   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
+   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
+   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
+   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
+   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting                    
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_init       : read namsbc namelist
+   !!   sbc            : surface ocean momentum, heat and freshwater boundary conditions
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE phycst           ! physical constants
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE trc_oce          ! shared ocean-passive tracers variables
+   USE sbc_ice          ! Surface boundary condition: ice fields
+   USE sbcdcy           ! surface boundary condition: diurnal cycle
+   USE sbcssm           ! surface boundary condition: sea-surface mean variables
+   USE sbcapr           ! surface boundary condition: atmospheric pressure
+   USE sbcana           ! surface boundary condition: analytical formulation
+   USE sbcflx           ! surface boundary condition: flux formulation
+   USE sbcblk_clio      ! surface boundary condition: bulk formulation : CLIO
+   USE sbcblk_core      ! surface boundary condition: bulk formulation : CORE
+   USE sbcblk_mfs       ! surface boundary condition: bulk formulation : MFS
+   USE sbcice_if        ! surface boundary condition: ice-if sea-ice model
+   USE sbcice_lim       ! surface boundary condition: LIM 3.0 sea-ice model
+   USE sbcice_lim_2     ! surface boundary condition: LIM 2.0 sea-ice model
+   USE sbcice_cice      ! surface boundary condition: CICE    sea-ice model
+   USE sbccpl           ! surface boundary condition: coupled florulation
+   USE cpl_oasis3       ! OASIS routines for coupling
+   USE sbcssr           ! surface boundary condition: sea surface restoring
+   USE sbcrnf           ! surface boundary condition: runoffs
+   USE sbcisf           ! surface boundary condition: ice shelf
+   USE sbcfwb           ! surface boundary condition: freshwater budget
+   USE closea           ! closed sea
+   USE icbstp           ! Icebergs!
+
+   USE prtctl           ! Print control                    (prt_ctl routine)
+   USE iom              ! IOM library
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! MPP library
+   USE timing           ! Timing
+   USE sbcwave          ! Wave module
+   USE bdy_par          ! Require lk_bdy
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc        ! routine called by step.F90
+   PUBLIC   sbc_init   ! routine called by opa.F90
+   
+   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
+      
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO-consortium (2011) 
+   !! $Id: sbcmod.F90 6459 2016-04-08 16:00:11Z lovato $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_init
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_init ***
+      !!
+      !! ** Purpose :   Initialisation of the ocean surface boundary computation
+      !!
+      !! ** Method  :   Read the namsbc namelist and set derived parameters
+      !!                Call init routines for all other SBC modules that have one
+      !!
+      !! ** Action  : - read namsbc parameters
+      !!              - nsbc: type of sbc
+      !!----------------------------------------------------------------------
+      INTEGER ::   icpt   ! local integer
+      !!
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   &
+         &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   &
+         &             nn_lsm    , nn_limflx , nn_components, ln_cpl
+      INTEGER  ::   ios
+      INTEGER  ::   ierr, ierr0, ierr1, ierr2, ierr3, jpm
+      LOGICAL  ::   ll_purecpl
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+
+      REWIND( numnam_ref )              ! Namelist namsbc in reference namelist : Surface boundary
+      READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namsbc )
+
+      !                          ! overwrite namelist parameter using CPP key information
+      IF( Agrif_Root() ) THEN                ! AGRIF zoom
+        IF( lk_lim2 )   nn_ice      = 2
+        IF( lk_lim3 )   nn_ice      = 3
+        IF( lk_cice )   nn_ice      = 4
+      ENDIF
+      IF( cp_cfg == 'gyre' ) THEN            ! GYRE configuration
+          ln_ana      = .TRUE.   
+          nn_ice      =   0
+      ENDIF
+
+      IF(lwp) THEN               ! Control print
+         WRITE(numout,*) '        Namelist namsbc (partly overwritten with CPP key setting)'
+         WRITE(numout,*) '           frequency update of sbc (and ice)             nn_fsbc     = ', nn_fsbc
+         WRITE(numout,*) '           Type of sbc : '
+         WRITE(numout,*) '              analytical formulation                     ln_ana      = ', ln_ana
+         WRITE(numout,*) '              flux       formulation                     ln_flx      = ', ln_flx
+         WRITE(numout,*) '              CLIO bulk  formulation                     ln_blk_clio = ', ln_blk_clio
+         WRITE(numout,*) '              CORE bulk  formulation                     ln_blk_core = ', ln_blk_core
+         WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
+         WRITE(numout,*) '              ocean-atmosphere coupled formulation       ln_cpl      = ', ln_cpl
+         WRITE(numout,*) '              forced-coupled mixed formulation           ln_mixcpl   = ', ln_mixcpl
+         WRITE(numout,*) '              OASIS coupling (with atm or sas)           lk_oasis    = ', lk_oasis
+         WRITE(numout,*) '              components of your executable              nn_components = ', nn_components
+         WRITE(numout,*) '              Multicategory heat flux formulation (LIM3) nn_limflx   = ', nn_limflx
+         WRITE(numout,*) '           Misc. options of sbc : '
+         WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
+         WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
+         WRITE(numout,*) '              ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd = ', nn_ice_embd
+         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
+         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
+         WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf
+         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr
+         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb
+         WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
+         WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm
+      ENDIF
+
+      ! LIM3 Multi-category heat flux formulation
+      SELECT CASE ( nn_limflx)
+      CASE ( -1 )
+         IF(lwp) WRITE(numout,*) '              Use of per-category fluxes (nn_limflx = -1) '
+      CASE ( 0  )
+         IF(lwp) WRITE(numout,*) '              Average per-category fluxes (nn_limflx = 0) ' 
+      CASE ( 1  )
+         IF(lwp) WRITE(numout,*) '              Average then redistribute per-category fluxes (nn_limflx = 1) '
+      CASE ( 2  )
+         IF(lwp) WRITE(numout,*) '              Redistribute a single flux over categories (nn_limflx = 2) '
+      END SELECT
+      !
+      IF ( nn_components /= jp_iam_nemo .AND. .NOT. lk_oasis )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
+      IF ( nn_components == jp_iam_opa .AND. ln_cpl )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA' )
+      IF ( nn_components == jp_iam_opa .AND. ln_mixcpl )   &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      IF ( ln_cpl .AND. .NOT. lk_oasis )    &
+         &      CALL ctl_stop( 'STOP', 'sbc_init : OASIS-coupled atmosphere model, but key_oasis3 disabled' )
+      IF( ln_mixcpl .AND. .NOT. lk_oasis )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires the cpp key key_oasis3' )
+      IF( ln_mixcpl .AND. .NOT. ln_cpl )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) requires ln_cpl = T' )
+      IF( ln_mixcpl .AND. nn_components /= jp_iam_nemo )    &
+         &      CALL ctl_stop( 'the forced-coupled mixed mode (ln_mixcpl) is not yet working with sas-opa coupling via oasis' )
+
+      !                              ! allocate sbc arrays
+      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
+
+      !                          ! Checks:
+      IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf 
+         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
+         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp
+         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp
+         rdivisf       = 0.0_wp
+      END IF
+      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero
+
+      sfx(:,:) = 0.0_wp                            ! the salt flux due to freezing/melting will be computed (i.e. will be non-zero) 
+                                                   ! only if sea-ice is present
+ 
+      fmmflx(:,:) = 0.0_wp                        ! freezing-melting array initialisation
+      
+      taum(:,:) = 0.0_wp                           ! Initialise taum for use in gls in case of reduced restart
+
+      !                                            ! restartability   
+      IF( ( nn_ice == 2 .OR. nn_ice ==3 ) .AND. .NOT.( ln_blk_clio .OR. ln_blk_core .OR. ln_cpl ) )   &
+         &   CALL ctl_stop( 'LIM sea-ice model requires a bulk formulation or coupled configuration' )
+      IF( nn_ice == 4 .AND. .NOT.( ln_blk_core .OR. ln_cpl ) )   &
+         &   CALL ctl_stop( 'CICE sea-ice model requires ln_blk_core or ln_cpl' )
+      IF( nn_ice == 4 .AND. lk_agrif )   &
+         &   CALL ctl_stop( 'CICE sea-ice model not currently available with AGRIF' )
+      IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
+         &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' )
+      IF( ( nn_ice /= 3 ) .AND. ( nn_limflx >= 0 ) )   &
+         &   WRITE(numout,*) 'The nn_limflx>=0 option has no effect if sea ice model is not LIM3'
+      IF( ( nn_ice == 3 ) .AND. ( ln_cpl ) .AND. ( ( nn_limflx == -1 ) .OR. ( nn_limflx == 1 ) ) )   &
+         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
+      IF( ( nn_ice == 3 ) .AND. ( .NOT. ln_cpl ) .AND. ( nn_limflx == 2 ) )   &
+         &   CALL ctl_stop( 'The chosen nn_limflx for LIM3 in forced mode cannot be 2' )
+
+      IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
+
+      IF( ln_dm2dc .AND. .NOT.( ln_flx .OR. ln_blk_core ) .AND. nn_components /= jp_iam_opa )   &
+         &   CALL ctl_stop( 'diurnal cycle into qsr field from daily values requires a flux or core-bulk formulation' )
+      
+      IF ( ln_wave ) THEN
+      !Activated wave module but neither drag nor stokes drift activated
+         IF ( .NOT.(ln_cdgw .OR. ln_sdw) )   THEN
+            CALL ctl_warn( 'Ask for wave coupling but nor drag coefficient (ln_cdgw=F) neither stokes drift activated (ln_sdw=F)' )
+      !drag coefficient read from wave model definable only with mfs bulk formulae and core 
+         ELSEIF (ln_cdgw .AND. .NOT.(ln_blk_mfs .OR. ln_blk_core) )       THEN       
+             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
+         ENDIF
+      ELSE
+      IF ( ln_cdgw .OR. ln_sdw  )                                         & 
+         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
+         & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
+      ENDIF 
+      !                          ! Choice of the Surface Boudary Condition (set nsbc)
+      ll_purecpl = ln_cpl .AND. .NOT. ln_mixcpl
+      !
+      icpt = 0
+      IF( ln_ana          ) THEN   ;   nsbc = jp_ana     ; icpt = icpt + 1   ;   ENDIF       ! analytical           formulation
+      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
+      IF( ln_blk_clio     ) THEN   ;   nsbc = jp_clio    ; icpt = icpt + 1   ;   ENDIF       ! CLIO bulk            formulation
+      IF( ln_blk_core     ) THEN   ;   nsbc = jp_core    ; icpt = icpt + 1   ;   ENDIF       ! CORE bulk            formulation
+      IF( ln_blk_mfs      ) THEN   ;   nsbc = jp_mfs     ; icpt = icpt + 1   ;   ENDIF       ! MFS  bulk            formulation
+      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
+      IF( cp_cfg == 'gyre') THEN   ;   nsbc = jp_gyre                        ;   ENDIF       ! GYRE analytical      formulation
+      IF( nn_components == jp_iam_opa )   &
+         &                  THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
+      IF( lk_esopa        )            nsbc = jp_esopa                                       ! esopa test, ALL formulations
+      !
+      IF( icpt /= 1 .AND. .NOT.lk_esopa ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '           E R R O R in setting the sbc, one and only one namelist/CPP key option '
+         WRITE(numout,*) '                     must be choosen. You choose ', icpt, ' option(s)'
+         WRITE(numout,*) '                     We stop'
+         nstop = nstop + 1
+      ENDIF
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( nsbc == jp_esopa   )   WRITE(numout,*) '              ESOPA test All surface boundary conditions'
+         IF( nsbc == jp_gyre    )   WRITE(numout,*) '              GYRE analytical formulation'
+         IF( nsbc == jp_ana     )   WRITE(numout,*) '              analytical formulation'
+         IF( nsbc == jp_flx     )   WRITE(numout,*) '              flux formulation'
+         IF( nsbc == jp_clio    )   WRITE(numout,*) '              CLIO bulk formulation'
+         IF( nsbc == jp_core    )   WRITE(numout,*) '              CORE bulk formulation'
+         IF( nsbc == jp_purecpl )   WRITE(numout,*) '              pure coupled formulation'
+         IF( nsbc == jp_mfs     )   WRITE(numout,*) '              MFS Bulk formulation'
+         IF( nsbc == jp_none    )   WRITE(numout,*) '              OPA coupled to SAS via oasis'
+         IF( ln_mixcpl          )   WRITE(numout,*) '              + forced-coupled mixed formulation'
+         IF( nn_components/= jp_iam_nemo )  &
+            &                       WRITE(numout,*) '              + OASIS coupled SAS'
+      ENDIF
+      !
+      IF( lk_oasis )   CALL sbc_cpl_init (nn_ice)   ! OASIS initialisation. must be done before: (1) first time step
+      !                                                     !                                            (2) the use of nn_fsbc
+
+!     nn_fsbc initialization if OPA-SAS coupling via OASIS
+!     sas model time step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
+      IF ( nn_components /= jp_iam_nemo ) THEN
+
+         IF ( nn_components == jp_iam_opa ) nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
+         IF ( nn_components == jp_iam_sas ) nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
+         !
+         IF(lwp)THEN
+            WRITE(numout,*)
+            WRITE(numout,*)"   OPA-SAS coupled via OASIS : nn_fsbc re-defined from OASIS namcouple ", nn_fsbc
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+
+      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
+          MOD( nstock             , nn_fsbc) /= 0 ) THEN 
+         WRITE(ctmp1,*) 'experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
+         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
+      ENDIF
+      !
+      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
+         &  CALL ctl_warn( 'nn_fsbc is NOT a multiple of the number of time steps in a day' )
+      !
+      IF( ln_dm2dc .AND. ( ( NINT(rday) / ( nn_fsbc * NINT(rdt) ) )  < 8 ) )   &
+         &   CALL ctl_warn( 'diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+
+                               CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
+      !
+      IF( ln_ssr           )   CALL sbc_ssr_init               ! Sea-Surface Restoring initialisation
+      !
+      IF( nn_isf   /= 0    )   CALL sbc_isf_init               ! Compute iceshelves
+
+                               CALL sbc_rnf_init               ! Runof initialisation
+      !
+      IF( nn_ice == 3      )   CALL sbc_lim_init               ! LIM3 initialisation
+
+      IF( nn_ice == 4      )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
+      
+   END SUBROUTINE sbc_init
+
+
+   SUBROUTINE sbc( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the ocean surface boundary
+      !!                condition (momentum, heat and freshwater fluxes)
+      !!
+      !! ** Method  :   blah blah  to be written ????????? 
+      !!                CAUTION : never mask the surface stress field (tke sbc)
+      !!
+      !! ** Action  : - set the ocean surface boundary condition at before and now 
+      !!                time step, i.e.  
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
+      !!              - updte the ice fraction : fr_i
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt       ! ocean time step
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('sbc')
+      !
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! qsr_b (:,:) = qsr (:,:)
+         emp_b(:,:) = emp(:,:)
+         sfx_b(:,:) = sfx(:,:)
+         IF ( ln_rnf ) THEN
+            rnf_b    (:,:  ) = rnf    (:,:  )
+            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+         ENDIF
+         IF( nn_isf /= 0  )  THEN
+            fwfisf_b  (:,:  ) = fwfisf  (:,:  )               
+            risf_tsc_b(:,:,:) = risf_tsc(:,:,:)              
+         ENDIF
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+      !
+      IF ( .NOT. lk_bdy ) then
+         IF( ln_apr_dyn ) CALL sbc_apr( kt )                ! atmospheric pressure provided at kt+0.5*nn_fsbc
+      ENDIF
+                                                         ! (caution called before sbc_ssm)
+      !
+      IF( nn_components /= jp_iam_sas )   CALL sbc_ssm( kt )   ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      !                                                        ! averaged over nf_sbc time-step
+
+      IF (ln_wave) CALL sbc_wave( kt )
+                                                   !==  sbc formulation  ==!
+                                                            
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
+      CASE( jp_gyre  )   ;   CALL sbc_gyre    ( kt )                    ! analytical formulation : GYRE configuration
+      CASE( jp_ana   )   ;   CALL sbc_ana     ( kt )                    ! analytical formulation : uniform sbc
+      CASE( jp_flx   )   ;   CALL sbc_flx     ( kt )                    ! flux formulation
+      CASE( jp_clio  )   ;   CALL sbc_blk_clio( kt )                    ! bulk formulation : CLIO for the ocean
+      CASE( jp_core  )   
+         IF( nn_components == jp_iam_sas ) &
+            &                CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: SAS receiving fields from OPA 
+                             CALL sbc_blk_core( kt )                    ! bulk formulation : CORE for the ocean
+                                                                        ! from oce: sea surface variables (sst_m, sss_m,  ssu_m,  ssv_m)
+      CASE( jp_purecpl )  ;  CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! pure coupled formulation
+                                                                        !
+      CASE( jp_mfs   )   ;   CALL sbc_blk_mfs ( kt )                    ! bulk formulation : MFS for the ocean
+      CASE( jp_none  ) 
+         IF( nn_components == jp_iam_opa ) &
+                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! OPA-SAS coupling: OPA receiving fields from SAS
+      CASE( jp_esopa )                                
+                             CALL sbc_ana     ( kt )                    ! ESOPA, test ALL the formulations
+                             CALL sbc_gyre    ( kt )                    !
+                             CALL sbc_flx     ( kt )                    !
+                             CALL sbc_blk_clio( kt )                    !
+                             CALL sbc_blk_core( kt )                    !
+                             CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   !
+      END SELECT
+
+      IF( ln_mixcpl )        CALL sbc_cpl_rcv ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing
+
+
+      !                                            !==  Misc. Options  ==!
+      
+      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
+      CASE(  1 )   ;         CALL sbc_ice_if   ( kt )                ! Ice-cover climatology ("Ice-if" model)
+      CASE(  2 )   ;         CALL sbc_ice_lim_2( kt, nsbc )          ! LIM-2 ice model
+      CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model
+      CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model
+      END SELECT                                              
+
+      IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs
+
+      IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves
+
+      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes
+ 
+      IF( ln_ssr         )   CALL sbc_ssr( kt )                   ! add SST/SSS damping term
+
+      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc )  ! control the freshwater budget
+
+      IF( nn_closea == 1 )   CALL sbc_clo( kt )                   ! treatment of closed sea in the model domain 
+      !                                                           ! (update freshwater fluxes)
+!RBbug do not understand why see ticket 667
+!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
+      CALL lbc_lnk( emp, 'T', 1. )
+      !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b  )   ! before non solar heat flux (T-point)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b  ) ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b  )    ! before     freshwater flux (T-point)
+            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
+            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b )  ! before salt flux (T-point)
+            ELSE
+               sfx_b (:,:) = sfx(:,:)
+            ENDIF
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:) 
+            vtau_b(:,:) = vtau(:,:)
+            qns_b (:,:) = qns (:,:)
+            emp_b (:,:) = emp(:,:)
+            sfx_b (:,:) = sfx(:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp  )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx  )
+      ENDIF
+
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
+         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
+         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux  
+                                                                ! (includes virtual salt flux beneath ice 
+                                                                ! in linear free surface case)
+         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux
+         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux 
+         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
+         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
+         IF( nn_ice > 0 .OR. nn_components == jp_iam_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction 
+         CALL iom_put( "taum"  , taum       )                   ! wind stress module 
+         CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
+      ENDIF
+      !
+      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at 
+      CALL iom_put( "vtau", vtau )   ! j-wind stress    each time step in sea-ice)
+      !
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 )
+         CALL prt_ctl(tab3d_1=tmask            , clinfo1=' tmask    - : ', mask1=tmask, ovlap=1, kdim=jpk )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' sst      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sss      - : ', mask1=tmask, ovlap=1, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau             , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau             , clinfo2=' vtau     - : ', mask2=vmask, ovlap=1 )
+      ENDIF
+
+      IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc')
+      !
+   END SUBROUTINE sbc
+
+
+   SUBROUTINE sbc_final
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_final  ***
+      !!
+      !! ** Purpose :   Finalize CICE (if used)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_ice == 4 )   CALL cice_sbc_final
+      !
+   END SUBROUTINE sbc_final
+
+   !!======================================================================
+END MODULE sbcmod
Index: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trasbc.F90
===================================================================
--- CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trasbc.F90	(revision 3136)
+++ CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trasbc.F90	(revision 3136)
@@ -0,0 +1,292 @@
+MODULE trasbc
+   !!==============================================================================
+   !!                       ***  MODULE  trasbc  ***
+   !! Ocean active tracers:  surface boundary condition
+   !!==============================================================================
+   !! History :  OPA  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
+   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
+   !!  NEMO      1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_sbc      : update the tracer trend at ocean surface
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE dom_oce         ! ocean space domain variables
+   USE phycst          ! physical constant
+   USE sbcmod          ! ln_rnf  
+   USE sbcrnf          ! River runoff  
+   USE sbcisf          ! Ice shelf   
+   USE traqsr          ! solar radiation penetration
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
+   !
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+   USE iom
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
+   USE eosbn2
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_sbc    ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id: trasbc.F90 6688 2016-06-13 12:50:45Z lovato $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_sbc ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_sbc  ***
+      !!                   
+      !! ** Purpose :   Compute the tracer surface boundary condition trend of
+      !!      (flux through the interface, concentration/dilution effect)
+      !!      and add it to the general trend of tracer equations.
+      !!
+      !! ** Method :
+      !!      Following Roullet and Madec (2000), the air-sea flux can be divided 
+      !!      into three effects: (1) Fext, external forcing; 
+      !!      (2) Fwi, concentration/dilution effect due to water exchanged 
+      !!         at the surface by evaporation, precipitations and runoff (E-P-R); 
+      !!      (3) Fwe, tracer carried with the water that is exchanged. 
+      !!            - salinity    : salt flux only due to freezing/melting
+      !!            sa = sa +  sfx / rau0 / e3t  for k=1
+      !!
+      !!      Fext, flux through the air-sea interface for temperature and salt: 
+      !!            - temperature : heat flux q (w/m2). If penetrative solar
+      !!         radiation q is only the non solar part of the heat flux, the
+      !!         solar part is added in traqsr.F routine.
+      !!            ta = ta + q /(rau0 rcp e3t)  for k=1
+      !!            - salinity    : no salt flux
+      !!
+      !!      The formulation for Fwb and Fwi vary according to the free 
+      !!      surface formulation (linear or variable volume). 
+      !!      * Linear free surface
+      !!            The surface freshwater flux modifies the ocean volume
+      !!         and thus the concentration of a tracer and the temperature.
+      !!         First order of the effect of surface freshwater exchange 
+      !!         for salinity, it can be neglected on temperature (especially
+      !!         as the temperature of precipitations and runoffs is usually
+      !!         unknown).
+      !!            - temperature : we assume that the temperature of both
+      !!         precipitations and runoffs is equal to the SST, thus there
+      !!         is no additional flux since in this case, the concentration
+      !!         dilution effect is balanced by the net heat flux associated
+      !!         to the freshwater exchange (Fwe+Fwi=0):
+      !!            (Tp P - Te E) + SST (P-E) = 0 when Tp=Te=SST
+      !!            - salinity    : evaporation, precipitation and runoff
+      !!         water has a zero salinity  but there is a salt flux due to 
+      !!         freezing/melting, thus:
+      !!            sa = sa + emp * sn / rau0 / e3t   for k=1
+      !!                    + sfx    / rau0 / e3t
+      !!         where emp, the surface freshwater budget (evaporation minus
+      !!         precipitation minus runoff) given in kg/m2/s is divided
+      !!         by rau0 (density of sea water) to obtain m/s.    
+      !!         Note: even though Fwe does not appear explicitly for 
+      !!         temperature in this routine, the heat carried by the water
+      !!         exchanged through the surface is part of the total heat flux
+      !!         forcing and must be taken into account in the global heat
+      !!         balance).
+      !!      * nonlinear free surface (variable volume, lk_vvl)
+      !!         contrary to the linear free surface case, Fwi is properly 
+      !!         taken into account by using the true layer thicknesses to       
+      !!         calculate tracer content and advection. There is no need to 
+      !!         deal with it in this routine.
+      !!           - temperature: Fwe=SST (P-E+R) is added to Fext.
+      !!           - salinity:  Fwe = 0, there is no surface flux of salt.
+      !!
+      !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
+      !!                with the tracer surface boundary condition 
+      !!              - send trends to trdtra module (l_trdtra=T)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices  
+      INTEGER  ::   ikt, ikb 
+      INTEGER  ::   nk_isf
+      REAL(wp) ::   zfact, z1_e3t, zdep
+      REAL(wp) ::   zalpha, zhk
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_sbc')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_sbc : TRAcer Surface Boundary Condition'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+
+      IF( l_trdtra ) THEN                    !* Save ta and sa trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+      ENDIF
+
+!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration
+      IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration
+         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns
+         qsr(:,:) = 0.e0                     ! qsr set to zero
+      ENDIF
+
+      !----------------------------------------
+      !        EMP, SFX and QNS effects
+      !----------------------------------------
+      !                                          Set before sbc tracer content fields
+      !                                          ************************************
+      IF( kt == nit000 ) THEN                      ! Set the forcing field at nit000 - 1
+         !                                         ! -----------------------------------
+         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
+              & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
+            zfact = 0.5_wp
+            CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
+            CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
+         ELSE                                         ! No restart or restart not found: Euler forward time stepping
+            zfact = 1._wp
+            sbc_tsc(:,:,:) = 0._wp
+            sbc_tsc_b(:,:,:) = 0._wp
+         ENDIF
+      ELSE                                         ! Swap of forcing fields
+         !                                         ! ----------------------
+         zfact = 0.5_wp
+         sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
+      ENDIF
+      !                                          Compute now sbc tracer content fields
+      !                                          *************************************
+
+                                                   ! Concentration dilution effect on (t,s) due to  
+                                                   ! evaporation, precipitation and qns, but not river runoff 
+                                               
+      IF( lk_vvl ) THEN                            ! Variable Volume case  ==>> heat content of mass flux is in qns
+         DO jj = 1, jpj
+            DO ji = 1, jpi 
+               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                              ! non solar heat flux
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0     * sfx(ji,jj)                              ! salt flux due to freezing/melting
+            END DO
+         END DO
+      ELSE                                         ! Constant Volume case ==>> Concentration dilution effect
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! temperature : heat flux
+               sbc_tsc(ji,jj,jp_tem) = r1_rau0_rcp * qns(ji,jj)                          &   ! non solar heat flux
+                  &                  + r1_rau0     * emp(ji,jj)  * tsn(ji,jj,1,jp_tem)       ! concent./dilut. effect
+               ! salinity    : salt flux + concent./dilut. effect (both in sfx)
+               sbc_tsc(ji,jj,jp_sal) = r1_rau0  * (  sfx(ji,jj)                          &   ! salt flux (freezing/melting)
+                  &                                + emp(ji,jj) * tsn(ji,jj,1,jp_sal) )      ! concent./dilut. effect
+            END DO
+         END DO
+         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
+         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
+      ENDIF
+      ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
+      DO jn = 1, jpts
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               z1_e3t = zfact / fse3t(ji,jj,1)
+               tsa(ji,jj,1,jn) = tsa(ji,jj,1,jn) + ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) * z1_e3t
+            END DO
+         END DO
+      END DO
+      !                                          Write in the ocean restart file
+      !                                          *******************************
+      IF( lrst_oce ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface tracer content forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) )
+         CALL iom_rstput( kt, nitrst, numrow, 'sbc_sc_b', sbc_tsc(:,:,jp_sal) )
+      ENDIF
+      !
+      !
+      !----------------------------------------
+      !       Ice Shelf effects (ISF)
+      !     tbl treated as in Losh (2008) JGR
+      !----------------------------------------
+      !
+      IF( nn_isf > 0 ) THEN
+         zfact = 0.5e0
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1
+         
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
+   
+               ! level fully include in the ice shelf boundary layer
+               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
+               ! sign - because fwf sign of evapo (rnf sign of precip)
+               DO jk = ikt, ikb - 1
+               ! compute tfreez for the temperature correction (we add water at freezing temperature)
+               ! compute trend
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj)
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
+               END DO
+   
+               ! level partially include in ice shelf boundary layer 
+               ! compute tfreez for the temperature correction (we add water at freezing temperature)
+               ! compute trend
+               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
+               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
+            END DO
+         END DO
+      END IF
+      !
+      !----------------------------------------
+      !        River Runoff effects
+      !----------------------------------------
+      !
+      IF( ln_rnf ) THEN         ! input of heat and salt due to river runoff 
+         zfact = 0.5_wp
+         DO jj = 2, jpj 
+            DO ji = fs_2, fs_jpim1
+               IF( rnf(ji,jj) /= 0._wp ) THEN
+                  zdep = zfact / h_rnf(ji,jj)
+                  DO jk = 1, nk_rnf(ji,jj)
+                                        tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   &
+                                          &               +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
+                     IF( ln_rnf_sal )   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   &
+                                          &               +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
+                  END DO
+               ENDIF
+            END DO  
+         END DO  
+      ENDIF
+
+      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*tsn(:,:,1,jp_tem) )   ! runoff term on sst
+      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*tsn(:,:,1,jp_sal) )   ! runoff term on sss
+
+      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
+      ENDIF
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_sbc')
+      !
+   END SUBROUTINE tra_sbc
+
+   !!======================================================================
+END MODULE trasbc
