Index: branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90
===================================================================
--- branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 9849)
+++ branches/UKMO/dev_isf_flx_UKESM_r9321/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90	(revision 9855)
@@ -18,13 +18,13 @@
    USE eosbn2          ! equation of state
    USE sbc_oce         ! surface boundary condition: ocean fields
+   USE zdfbfr          !
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
+   USE fldread         ! read input field at current time step
    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
-   USE lib_fortran, ONLY: glob_sum
 
    IMPLICIT NONE
@@ -35,37 +35,35 @@
    ! 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, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2]
    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
+   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
+   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
+   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
+   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
 
    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 (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m]
    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
-
-
-   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 ?
+   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base
+
+   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K]
+   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s]
+   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3]
+   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C]
+   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg]
 
 !: 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
+   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
+   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
+   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
+   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
    
    !! * Substitutions
@@ -76,26 +74,38 @@
    !! 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, ijkmin, inum, ierror
-    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
-    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
-    REAL(wp)                     ::   rmin
-    REAL(wp)                     ::   zhk
-    REAL(wp)                     ::   zt_frz, zpress
-    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
-
-    REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
-    REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
-      !
       !!---------------------------------------------------------------------
-      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
+      !!                  ***  ROUTINE sbc_isf  ***
+      !!
+      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf 
+      !!              melting and freezing 
+      !!
+      !! ** Method  :  4 parameterizations are available according to nn_isf 
+      !!               nn_isf = 1 : Realistic ice_shelf formulation
+      !!                        2 : Beckmann & Goose parameterization
+      !!                        3 : Specified runoff in deptht (Mathiot & al. )
+      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)          ::   kt         ! ocean time step
+      INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
+      INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+      REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
+      REAL(wp)                     ::   rmin
+      REAL(wp)                     ::   zhk
+      REAL(wp)                     ::   zt_frz, zpress
+      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
+
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zfwfisf3d, zqhcisf3d, zqlatisf3d
+      REAL(wp), DIMENSION(:,:  ), POINTER :: zqhcisf2d
+      REAL(wp), DIMENSION(:,:  ), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep) 
+      !
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, &
+                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
       !
       !
@@ -121,12 +131,6 @@
          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
@@ -185,9 +189,7 @@
             
             ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
-            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
+            ALLOCATE( sf_fwfisf(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' )
          END IF
          
@@ -206,4 +208,6 @@
 
       IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
+         ! allocation
+         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  )
 
          ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
@@ -229,11 +233,11 @@
 
          ! compute salf and heat flux
-         IF (nn_isf == 1) THEN
-            ! realistic ice shelf formulation
+         SELECT CASE ( nn_isf )
+         CASE ( 1 )    ! 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')
+            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
+            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
             ! iom print
             CALL iom_put('ttbl',ttbl(:,:))
@@ -244,13 +248,11 @@
             CALL sbc_isf_cav (kt)
 
-         ELSE IF (nn_isf == 2) THEN
-            ! Beckmann and Goosse parametrisation 
+         CASE ( 2 )    ! 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)
+         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
-            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf  flux from the isf (fwfisf <0 mean melting) 
 
             IF( lk_oasis) THEN
@@ -294,12 +296,10 @@
             ENDIF
 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux
             stbl(:,:)   = soce
 
-         ELSE IF (nn_isf == 4) THEN
-            ! specified fwf and heat flux forcing beneath the ice shelf
+         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
             CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
-            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
-            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+            fwfisf(:,:) =   sf_fwfisf(1)%fnow(:,:,1)            ! fwf  flux from the isf (fwfisf <0 mean melting)
 
             IF( lk_oasis) THEN
@@ -343,17 +343,22 @@
             ENDIF
 
-            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
-            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
+            qisf(:,:)   = fwfisf(:,:) * rlfusisf              ! heat        flux
             stbl(:,:)   = soce
 
-         END IF
+         END SELECT
+
          ! 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 !
+         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
+         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
+         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj))
+            END DO
+         END DO
+         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
          
-         ! salt effect already take into account in vertical advection
-         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
+         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
+         risf_tsc(:,:,jp_sal) = 0.0_wp
 
          ! output
@@ -361,12 +366,9 @@
          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.)
+         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.)
+         CALL lbc_lnk(qisf(:,:)           ,'T',1.)
 
 !=============================================================================================================================================
@@ -405,5 +407,5 @@
 !=============================================================================================================================================
 
-         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         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
@@ -416,9 +418,12 @@
                risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
             END IF
-         ENDIF
+         END IF
          ! 
+         ! deallocation
+         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  )
       END IF
-  
+      !  
   END SUBROUTINE sbc_isf
+
 
   INTEGER FUNCTION sbc_isf_alloc()
@@ -435,79 +440,68 @@
                &    STAT= sbc_isf_alloc )
          !
-         IF( lk_mpp                  )   CALL mpp_sum ( 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 IF
   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) 
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_isf_bg03  ***
+      !!
+      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
+      !!          into the adjacent ocean
+      !!
+      !! ** Method  :   See reference
+      !!
+      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
+      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
+      !!         (hereafter BG)
+      !! History :
+      !!         06-02  (C. Wang) Original code
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT ( in ) :: kt
+      !
+      INTEGER  :: ji, jj, jk ! dummy loop index
+      INTEGER  :: ik         ! current level
+      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')
+      !
+      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 > 1 ) THEN
+               ! 1. -----------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
+                  CALL eos_fzp(stbl(ji,jj), zt_frz, gdept_n(ji,jj,jk)) 
+                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
+               END DO
+               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
+               ! 2. ------------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,jk)
              
-             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')
+               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !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
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
+      !
   END SUBROUTINE sbc_isf_bg03
 
@@ -527,163 +521,134 @@
       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
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      INTEGER  ::   nit
       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
+      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zcfac
+      REAL(wp) ::   zeps = 1.e-20_wp        
+      REAL(wp) ::   zerr
+      REAL(wp), DIMENSION(:,:), POINTER ::   ztfrz, zsfrz
+      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas 
+      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b
+      LOGICAL  ::   lit
       !!---------------------------------------------------------------------
-      !
-      ! coeficient for linearisation of tfreez
-      zlamb1=-0.0575
-      zlamb2=0.0901
-      zlamb3=-7.61e-04
+      ! coeficient for linearisation of potential tfreez
+      ! Crude approximation for pressure (but commonly used)
+      IF ( ln_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
+         zlamb1 =-0.0564_wp
+         zlamb2 = 0.0773_wp
+         zlamb3 =-7.8633e-8 * grav * rau0
+      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
+         zlamb1 =-0.0573_wp
+         zlamb2 = 0.0832_wp
+         zlamb3 =-7.53e-8 * grav * rau0
+      ENDIF
+      !
       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
+      CALL wrk_alloc( jpi,jpj, ztfrz  , zsfrz  , zgammat, zgammas  )
+      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
+
+      ! initialisation
+      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
+      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp    
+      zfwflx (:,:) = 0.0_wp     ; zsfrz(:,:) = 0.0_wp
+
+      ! compute ice shelf melting
+      nit = 1 ; lit = .TRUE.
+      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
+         SELECT CASE ( nn_isfblk )
+         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
+            ! Calculate freezing temperature
+            CALL eos_fzp( stbl(:,:), ztfrz(:,:), risfdep(:,:) )
+
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+            
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-ztfrz(ji,jj))
+                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
+               END DO
+            END DO
+
+            ! Compute heat flux and upward fresh water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
+            ! compute gammat every where (2d)
+            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
+
+            ! compute upward heat flux zhtflx and upward water flux zwflx
+            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! compute coeficient to solve the 2nd order equation
+                  zeps1 = rcp*rau0*zgammat(ji,jj)
+                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
+                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
+                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
+                  zeps6 = zeps4-ttbl(ji,jj)
+                  zeps7 = zeps4-tsurf
+                  zaqe  = zlamb1 * (zeps1 + zeps3)
+                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
+                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
+                  zcqe  = zeps2*stbl(ji,jj)
+                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
+
+                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
+                  ! compute s freeze
+                  zsfrz(ji,jj)=(-zbqe-SQRT(zdis))*zaqer
+                  IF ( zsfrz(ji,jj) < 0.0_wp ) zsfrz(ji,jj)=(-zbqe+SQRT(zdis))*zaqer
+
+                  ! compute t freeze (eq. 22)
+                  ztfrz(ji,jj)=zeps4+zlamb1*zsfrz(ji,jj)
+  
+                  ! zfwflx is upward water flux
+                  ! zhtflx is upward heat flux (out of ocean)
+                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
+                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz(ji,jj)-stbl(ji,jj)) / MAX(zsfrz(ji,jj),zeps)
+                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - ztfrz(ji,jj) ) 
+               END DO
+            END DO
+
+            ! compute heat and water flux
+            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+         END SELECT
+
+         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
+         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
+         ELSE                           
+            ! check total number of iteration
+            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
+            ELSE                 ; nit = nit + 1
+            END IF
+
+            ! compute error between 2 iterations
+            ! if needed save gammat and compute zhtflx_b for next iteration
+            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
+            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
+            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
+            END IF
+         END IF
       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 )
+      !
+      CALL iom_put('isftfrz'  , ztfrz * (1._wp - tmask(:,:,1)) * ssmask(:,:))
+      CALL iom_put('isfsfrz'  , zsfrz * (1._wp - tmask(:,:,1)) * ssmask(:,:))
+      CALL iom_put('isfgammat', zgammat)
+      CALL iom_put('isfgammas', zgammas)
+      ! 
+      CALL wrk_dealloc( jpi,jpj, ztfrz  , zgammat, zgammas  )
+      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b )
       !
       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 )
+   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
       !!----------------------------------------------------------------------
       !! ** Purpose    : compute the coefficient echange for heat flux
@@ -694,10 +659,10 @@
       !!                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), DIMENSION(:,:), INTENT(out) :: pgt, pgs
+      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf
+      !
+      INTEGER  :: ikt                        
+      INTEGER  :: ji, jj                     ! loop index
+      REAL(wp), DIMENSION(:,:), POINTER :: 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 
@@ -709,169 +674,212 @@
       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) :: zeps = 1.0e-20_wp    
+      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
+      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
       REAL(wp), DIMENSION(2) :: zts, zab
       !!---------------------------------------------------------------------
-      !
-      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)
+      CALL wrk_alloc( jpi,jpj, zustar )
+      !
+      SELECT CASE ( nn_gammablk )
+      CASE ( 0 ) ! gamma is constant (specified in namelist)
+         !! ISOMIP formulation (Hunter et al, 2006)
+         pgt(:,:) = rn_gammat0
+         pgs(:,:) = rn_gammas0
+
+      CASE ( 1 ) ! gamma is assume to be proportional to u*
+         !! Jenkins et al., 2010, JPO, p2298-2312
+         !! Adopted by Asay-Davis et al. (2015)
+
+         !! compute ustar (eq. 24)
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! Compute gammats
+         pgt(:,:) = zustar(:,:) * rn_gammat0
+         pgs(:,:) = zustar(:,:) * rn_gammas0
+      
+      CASE ( 2 ) ! gamma depends of stability of boundary layer
+         !! 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)
+         !! compute ustar
+         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) )
+
+         !! compute Pr and Sc number (can be improved)
+         zPr =   13.8_wp
+         zSc = 2432.0_wp
+
+         !! compute gamma mole
+         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
+         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
+
+         !! compute gamma
+         DO ji=2,jpi
+            DO jj=2,jpj
                ikt = mikt(ji,jj)
 
-      !! 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
+               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think
+                  pgt = rn_gammat0
+                  pgs = rn_gammas0
                ELSE
-      !! compute Rc number (as done in zdfric.F90)
-               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 Rc number (as done in zdfric.F90)
+                  zcoef = 0.5_wp / 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) / MAX( zdku*zdku + zdkv*zdkv, zeps )
+
+                  !! compute bouyancy 
+                  zts(jp_tem) = ttbl(ji,jj)
+                  zts(jp_sal) = stbl(ji,jj)
+                  zdep        = fsdepw(ji,jj,ikt)
                   !
-      !! 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)
+                  CALL eos_rab( zts, zdep, zab )
+                  !
+                  !! compute length scale 
+                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+                  !! compute Monin Obukov Length
+                  ! Maximum boundary layer depth
+                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp
+                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
+                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
+                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
+
+                  !! compute eta* (stability parameter)
+                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp)))
+
+                  !! compute the sublayer thickness
+                  zhnu = 5 * znu / zustar(ji,jj)
+
+                  !! compute gamma turb
+                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
+                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
+
+                  !! compute gammats
+                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
+                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
                END IF
-      END IF
-
-   END SUBROUTINE
-
-   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
+            END DO
+         END DO
+         CALL lbc_lnk(pgt(:,:),'T',1.)
+         CALL lbc_lnk(pgs(:,:),'T',1.)
+      END SELECT
+      CALL wrk_dealloc( jpi,jpj, zustar )
+      !
+   END SUBROUTINE sbc_isf_gammats
+
+   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
       !!----------------------------------------------------------------------
       !!                  ***  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
+      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin
+      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout
+      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out
+      !
+      REAL(wp) :: ze3, zhk
+      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl
+
+      INTEGER :: ji, jj, jk                  ! loop index
+      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl
+      !!----------------------------------------------------------------------
+      ! allocation
+      CALL wrk_alloc( jpi,jpj, zhisf_tbl)
       
-      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)
+      ! initialisation
+      pvarout(:,:)=0._wp
+   
+      SELECT CASE ( cd_ptin )
+      CASE ( 'U' ) ! compute U in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbku(ji,jj)
+                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3u_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+      
+      CASE ( 'V' ) ! compute V in the top boundary layer at T- point 
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
+               ! thickness of boundary layer at least the top level thickness
+               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt))
+
+               ! determine the deepest level influenced by the boundary layer
+               DO jk = ikt+1, mbkv(ji,jj)
+                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
+               END DO
+               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
+
+               ! level fully include in the ice shelf boundary layer
+               DO jk = ikt, ikb - 1
+                  ze3 = fse3v_n(ji,jj,jk)
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
+               END DO
+
+               ! level partially include in ice shelf boundary layer 
+               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+         DO jj = 2,jpj
+            DO ji = 2,jpi
+               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
+            END DO
+         END DO
+         CALL lbc_lnk(pvarout,'T',-1.)
+
+      CASE ( 'T' ) ! compute T in the top boundary layer at T- point 
+         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
                   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
+                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,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.)
-
+               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
+            END DO
+         END DO
+      END SELECT
+
+      ! mask mean tbl value
+      pvarout(:,:) = pvarout(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
+
+      ! deallocation
+      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )      
+      !
    END SUBROUTINE sbc_isf_tbl
       
@@ -889,16 +897,15 @@
       !! ** Action  :   phdivn   decreased by the runoff inflow
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
-      !!
+      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
+      REAL(wp) ::   zhk
+      REAL(wp) ::   zfact     ! local scalar
       !!----------------------------------------------------------------------
       !
       zfact   = 0.5_wp
       !
-      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
+      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water
          DO jj = 1,jpj
             DO ji = 1,jpi
@@ -909,7 +916,6 @@
 
                ! 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
+               DO jk = ikt+1, mbkt(ji,jj)
+                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  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.
@@ -917,10 +923,11 @@
                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
+               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
-      !
+      END IF 
+      !
+      !==   ice shelf melting distributed over several levels   ==!
       DO jj = 1,jpj
          DO ji = 1,jpi
@@ -930,89 +937,13 @@
                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
+                    &              * 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   ==!
+                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
          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
