Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8983)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8984)
@@ -154,16 +154,16 @@
    !!----------------------------------------------------------------------
    !                                     !!** ice-generic parameters namelist (nampar) **
-   INTEGER           , PUBLIC ::   jpl             !: number of ice  categories 
-   INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers 
-   INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers 
-   INTEGER           , PUBLIC ::   nn_monocat      !: virtual ITD mono-category parameterizations (1-4) or not (0)
-   LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F)
-   LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F)
-   REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere
-   REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere
-   CHARACTER(len=256), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
-   CHARACTER(len=256), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
-   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
-   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
+   INTEGER           , PUBLIC ::   jpl              !: number of ice  categories 
+   INTEGER           , PUBLIC ::   nlay_i           !: number of ice  layers 
+   INTEGER           , PUBLIC ::   nlay_s           !: number of snow layers 
+   INTEGER           , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1-4) or not (0)
+   LOGICAL           , PUBLIC ::   ln_icedyn        !: flag for ice dynamics (T) or not (F)
+   LOGICAL           , PUBLIC ::   ln_icethd        !: flag for ice thermo   (T) or not (F)
+   REAL(wp)          , PUBLIC ::   rn_amax_n        !: maximum ice concentration Northern hemisphere
+   REAL(wp)          , PUBLIC ::   rn_amax_s        !: maximum ice concentration Southern hemisphere
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_in     !: suffix of ice restart name (input)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_out    !: suffix of ice restart name (output)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir  !: ice restart input directory
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir !: ice restart output directory
 
    !                                     !!** ice-itd namelist (namitd) **
@@ -192,13 +192,16 @@
    !                                      !   =-1  Do nothing (needs N(cat) fluxes)
    !                                      !   = 0  Average N(cat) fluxes then apply the average over the N(cat) ice 
-   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice
-   !                                      !                                   using T-ice and albedo sensitivity
+   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity
    !                                      !   = 2  Redistribute a single flux over categories
-
-   INTEGER , PUBLIC            ::   nice_jules            !: Choice of jules coupling 
-   !                                                      !  Associated indices
-   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   !: no Jules coupling (ice thermodynamics forced via qsr and qns)
-   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
-   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
+   INTEGER , PUBLIC            ::   nice_jules           !: Choice of jules coupling 
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0  !: no Jules coupling (ice thermodynamics forced via qsr and qns)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1  !: emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it)
+   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2  !: active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90)
+
+   !                                     !!** ice-vertical diffusion namelist (namthd_zdf) **
+   LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964)
+   LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007)
+   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
+   REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]   
 
    !                                     !!** ice-salinity namelist (namthd_sal) **
@@ -211,5 +214,5 @@
    REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]
 
-   !                                     !!** namelist namthd_pnd
+   !                                     !!** ice-ponds namelist (namthd_pnd)
    LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012
    LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater
@@ -255,8 +258,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange   [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sni !: snow ice growth component of wfx_snw [kg.m-2.s-1]
-   ! MV MP 2016
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1]
-   ! END MV MP 2016
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice  [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice    [kg.m-2.s-1]
@@ -353,5 +354,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sz_i     !: ice salinity          [PSU]
 
-   ! MV MP 2016
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond fraction per grid cell area
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area [m]
@@ -361,5 +361,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area [m]
-   ! END MV MP 2016
 
    !!----------------------------------------------------------------------
@@ -376,6 +375,6 @@
    !! * Ice thickness distribution variables
    !!----------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max         !: Boundary of ice thickness categories in thickness space
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean        !: Mean ice thickness in catgories 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_max         !: Boundary of ice thickness categories in thickness space
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_mean        !: Mean ice thickness in catgories 
    !
    !!----------------------------------------------------------------------
@@ -385,14 +384,14 @@
    ! trp   ''         ''     ''       advection (transport of ice)
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vi   !: transport of ice volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vs   !: transport of snw volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_sv   !: transport of salt content
-   !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat     !: snw/ice heat content variation   [W/m2] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sice     !: ice salt content variation   [] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s] 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_ei   !: transport of ice enthalpy (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_es   !: transport of snw enthalpy (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_sv   !: transport of salt content
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_heat     !: snw/ice heat content variation   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_sice     !: ice salt content variation   [] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice     !: ice volume variation   [m/s] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s] 
 
    !
@@ -401,8 +400,8 @@
    !!----------------------------------------------------------------------
    ! Extra sea ice diagnostics to address the data request
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_si          !: Temperature at Snow-ice interface (K) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_si          !: Temperature at Snow-ice interface (K) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tm_si         !: mean temperature at the snow-ice interface (K) 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_fc_su    !: Surface conduction flux (W/m2)
 
    !
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8983)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8984)
@@ -251,10 +251,8 @@
             pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:)
          END DO
-         ! MV MP 2016
          IF ( ln_pnd_H12 ) THEN
-            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
-            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
+            pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:)
+            pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:)
          ENDIF
-         ! END MV MP 2016
       END DO
       !
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8983)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8984)
@@ -17,15 +17,14 @@
    !!  ice_thd_zdf      : select the appropriate routine for vertical heat diffusion calculation
    !!  ice_thd_zdf_BL99 : 
-   !!  ice_thd_enmelt   :
    !!  ice_thd_zdf_init :
    !!----------------------------------------------------------------------
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants (ocean directory) 
-   USE ice            ! sea-ice: variables
-   USE ice1D          ! sea-ice: thermodynamics variables
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants (ocean directory) 
+   USE ice             ! sea-ice: variables
+   USE icethd_zdf_BL99 ! sea-ice: vertical diffusion (Bitz and Lipscomb, 1999) 
    !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! MPP library
+   USE lib_fortran     ! fortran utilities (glob_sum + no signed zero)
 
    IMPLICIT NONE
@@ -35,19 +34,12 @@
    PUBLIC   ice_thd_zdf_init   ! called by icestp
 
+   INTEGER ::   nice_zdf       ! Choice of the type of vertical heat diffusion formulation
+   !                                 ! associated indices:
+   INTEGER, PARAMETER ::   np_BL99 = 1   ! Bitz and Lipscomb (1999)
+!! INTEGER, PARAMETER ::   np_XXXX = 2
+
    !!** namelist (namthd_zdf) **
-   LOGICAL  ::   ln_zdf_BL99      ! Heat diffusion follows Bitz and Lipscomb (1999)
-   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
-   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
-   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K]
+   LOGICAL ::   ln_zdf_BL99    ! Heat diffusion follows Bitz and Lipscomb (1999)
 
-   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation
-   !                                           ! associated indices:
-   INTEGER, PARAMETER ::   np_BL99 = 1         ! Bitz and Lipscomb (1999)
-
-   INTEGER, PARAMETER ::   np_zdf_jules_OFF = 0   !   compute all temperatures from qsr and qns
-   INTEGER, PARAMETER ::   np_zdf_jules_SND = 1   !   compute conductive heat flux and surface temperature from qsr and qns
-   INTEGER, PARAMETER ::   np_zdf_jules_RCV = 2   !   compute snow and inner ice temperatures from qcnd
-   
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2017)
@@ -67,13 +59,21 @@
       SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver
       !
-      !                             !-------------      
-      CASE( np_BL99 )               ! BL99 solver
-         !                          !-------------
+      !                             !-------------!      
+      CASE( np_BL99 )               ! BL99 solver !
+         !                          !-------------!
          SELECT CASE( nice_jules )
-         CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF )   ! No Jules coupler     
-         CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND )   ! Jules coupler is emulated (send/recieve)
-                                       CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )
-         CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV )   ! Jules coupler is active (receive only)
+         !                         ! No Jules coupler ==> default option
+         CASE( np_jules_OFF    )   ;   CALL ice_thd_zdf_BL99 ( np_jules_OFF    )
+         !
+         !                         ! Jules coupler is emulated => 1st call to get the needed fields (conduction...)
+         !                                                        2nd call to use these fields to calculate heat diffusion   
+         CASE( np_jules_EMULE  )   ;   CALL ice_thd_zdf_BL99 ( np_jules_EMULE  )
+                                       CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE )
+         !
+         !                         ! Jules coupler is active ==> Met Office default option
+         CASE( np_jules_ACTIVE )   ;   CALL ice_thd_zdf_BL99 ( np_jules_ACTIVE )
+         !
          END SELECT
+         !
       END SELECT
       
@@ -81,863 +81,4 @@
    
    
-   SUBROUTINE ice_thd_zdf_BL99( k_jules )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
-      !!
-      !! ** Purpose :   computes the time evolution of snow and sea-ice temperature
-      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm
-      !!
-      !! ** Method  :   solves the heat equation diffusion with a Neumann boundary
-      !!              condition at the surface and a Dirichlet one at the bottom. 
-      !!              Solar radiation is partially absorbed into the ice.
-      !!              The specific heat and thermal conductivities depend on ice 
-      !!              salinity and temperature to take into account brine pocket   
-      !!              melting. The numerical scheme is an iterative Crank-Nicolson
-      !!              on a non-uniform multilayer grid in the ice and snow system.
-      !!
-      !!           The successive steps of this routine are
-      !!           1.  initialization of ice-snow layers thicknesses
-      !!           2.  Internal absorbed and transmitted radiation
-      !!           Then iterative procedure begins
-      !!           3.  Thermal conductivity
-      !!           4.  Kappa factors
-      !!           5.  specific heat in the ice
-      !!           6.  eta factors
-      !!           7.  surface flux computation
-      !!           8.  tridiagonal system terms
-      !!           9.  solving the tridiagonal system with Gauss elimination
-      !!           Iterative procedure ends according to a criterion on evolution
-      !!           of temperature
-      !!           10. Fluxes at the interfaces
-      !!
-      !! ** Inputs / Ouputs : (global commons)
-      !!           surface temperature : t_su_1d
-      !!           ice/snow temperatures   : t_i_1d, t_s_1d
-      !!           ice salinities          : sz_i_1d
-      !!           number of layers in the ice/snow: nlay_i, nlay_s
-      !!           total ice/snow thickness : h_i_1d, h_s_1d
-      !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND)
-      !
-      INTEGER ::   ji, jk         ! spatial loop index
-      INTEGER ::   jm             ! current reference number of equation
-      INTEGER ::   jm_mint, jm_maxt
-      INTEGER ::   iconv          ! number of iterations in iterative procedure
-      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
-      !
-      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
-      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
-      !
-      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
-      REAL(wp) ::   zg1       =  2._wp        !
-      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
-      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
-      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
-      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
-      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
-      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
-      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
-      REAL(wp) ::   zcpi                      ! Ice specific heat
-      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
-      REAL(wp) ::   zfac                      ! dummy factor
-      !
-      REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow
-      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
-      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
-      REAL(wp), DIMENSION(jpij) ::   zh_s, z1_h_s ! snow layer thickness
-      REAL(wp), DIMENSION(jpij) ::   zqns_ice_b   ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(jpij) ::   zfnet        ! surface flux function
-      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function
-      !
-      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
-      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
-      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
-      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
-      !
-      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor
-      !
-      ! Mono-category
-      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
-      REAL(wp) :: zhe           ! dummy factor
-      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
-      !!------------------------------------------------------------------     
-
-      ! --- diag error on heat diffusion - PART 1 --- !
-      DO ji = 1, npti
-         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
-            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
-      END DO
-
-      !------------------
-      ! 1) Initialization
-      !------------------
-      DO ji = 1, npti
-         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
-         ! layer thickness
-         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
-         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
-      END DO
-      !
-      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
-      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
-      END WHERE
-      !
-      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
-      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
-      END WHERE
-      !
-      ! Store initial temperatures and non solar heat fluxes
-      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
-         !
-         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1
-         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value
-         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux 
-         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
-         !
-         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
-         !
-      ENDIF   
-      !
-      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
-      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
-
-      !-------------
-      ! 2) Radiation
-      !-------------
-      ! --- Transmission/absorption of solar radiation in the ice --- !
-!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value
-!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1
-!
-!     DO ji = 1, npti
-!
-!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )     
-!
-!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness
-!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque
-!
-!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation 
-!
-!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji)
-!              zftrice(ji) = qsr_ice_tr_1d(ji)
-!              i0(ji) = zfrqsr_tr_i
-!
-!     END DO
-
-      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)
-      DO jk = 1, nlay_s
-         DO ji = 1, npti
-            !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
-            !                             ! radiation absorbed by the layer-th snow layer
-            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
-         END DO
-      END DO
-      !
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )
-      DO jk = 1, nlay_i 
-         DO ji = 1, npti
-            !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
-            !                             ! radiation absorbed by the layer-th ice layer
-            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
-         END DO
-      END DO
-      !
-      ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
-      !
-      iconv    =  0          ! number of iterations
-      zdti_max =  1000._wp   ! maximal value of error on all points
-      !
-      !                                                          !============================!
-      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
-         !                                                       !============================!
-         iconv = iconv + 1
-         !
-         ztib(1:npti,:) = t_i_1d(1:npti,:)
-         ztsb(1:npti,:) = t_s_1d(1:npti,:)
-         !
-         !--------------------------------
-         ! 3) Sea ice thermal conductivity
-         !--------------------------------
-         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
-            !
-            DO ji = 1, npti
-               ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
-                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
-               END DO
-            END DO
-            !
-         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
-            !
-            DO ji = 1, npti
-               ztcond_i(ji,0)      = rcdic + 0.09_wp * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
-                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
-                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
-                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
-               END DO
-            END DO
-            !
-         ENDIF
-         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )        
-         !
-         !--- G(he) : enhancement of thermal conductivity in mono-category case
-         ! Computation of effective thermal conductivity G(h)
-         ! Used in mono-category case only to simulate an ITD implicitly
-         ! Fichefet and Morales Maqueda, JGR 1997
-         zghe(1:npti) = 1._wp
-         !
-         SELECT CASE ( nn_monocat )
-         !
-         CASE ( 1 , 3 )
-            !
-            zepsilon = 0.1_wp
-            DO ji = 1, npti
-               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity
-               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )     ! Effective thickness he (zhe)
-               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
-                  zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he)
-               ENDIF
-            END DO
-            !
-         END SELECT
-         !
-         !-----------------
-         ! 4) kappa factors
-         !-----------------
-         !--- Snow
-         DO jk = 0, nlay_s-1
-            DO ji = 1, npti
-               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
-            END DO
-         END DO
-         DO ji = 1, npti   ! Snow-ice interface
-            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
-            IF( zfac > epsi10 ) THEN
-               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
-            ELSE
-               zkappa_s(ji,nlay_s) = 0._wp
-            ENDIF
-         END DO
-
-         !--- Ice
-         DO jk = 0, nlay_i
-            DO ji = 1, npti
-               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
-            END DO
-         END DO
-         DO ji = 1, npti   ! Snow-ice interface
-            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
-         END DO
-         !
-         !--------------------------------------
-         ! 5) Sea ice specific heat, eta factors
-         !--------------------------------------
-         DO jk = 1, nlay_i
-            DO ji = 1, npti
-               zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
-            END DO
-         END DO
-
-         DO jk = 1, nlay_s
-            DO ji = 1, npti
-               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
-            END DO
-         END DO
-         
-         !
-         !----------------------------------------!
-         !                                        !
-         !       JULES IF (OFF or SND MODE)       !
-         !                                        !
-         !----------------------------------------!
-         !
-         
-         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode
-         
-            !
-            ! In OFF mode the original BL99 temperature computation is used
-            ! (with qsr_ice, qns_ice and dqns_ice as inputs)
-            !
-            ! In SND mode, the computation is required to compute the conduction fluxes
-            !
-         
-            !----------------------------
-            ! 6) surface flux computation
-            !----------------------------
-
-            ! update of the non solar flux according to the update in T_su
-            DO ji = 1, npti
-               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
-            END DO
-
-            DO ji = 1, npti
-               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar
-            END DO
-            !
-            !----------------------------
-            ! 7) tridiagonal system terms
-            !----------------------------
-            !!layer denotes the number of the layer in the snow or in the ice
-            !!jm denotes the reference number of the equation in the tridiagonal
-            !!system, terms of tridiagonal system are indexed as following :
-            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
-            !!ice interior terms (top equation has the same form as the others)
-            ztrid   (1:npti,:,:) = 0._wp
-            zindterm(1:npti,:)   = 0._wp
-            zindtbis(1:npti,:)   = 0._wp
-            zdiagbis(1:npti,:)   = 0._wp
-
-            DO jm = nlay_s + 2, nlay_s + nlay_i 
-            DO ji = 1, npti
-               jk = jm - nlay_s - 1
-               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
-               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
-               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-            END DO
-            END DO
-
-            jm =  nlay_s + nlay_i + 1
-            DO ji = 1, npti
-            !!ice bottom term
-            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
-            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
-            ztrid(ji,jm,3)  = 0.0
-            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
-            END DO
-
-
-            DO ji = 1, npti
-            !                               !---------------------!
-            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells !
-               !                            !---------------------!
-               ! snow interior terms (bottom equation has the same form as the others)
-               DO jm = 3, nlay_s + 1
-                  jk = jm - 1
-                  ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-                  ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                  ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                  zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
-               END DO
-
-               ! case of only one layer in the ice (ice equation is altered)
-               IF ( nlay_i == 1 ) THEN
-                  ztrid(ji,nlay_s+2,3)  = 0.0
-                  zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
-               ENDIF
-
-               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
-
-                  jm_min(ji) = 1
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! surface equation
-                  ztrid(ji,1,1)  = 0.0
-                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
-                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
-
-                  ! first layer of snow equation
-                  ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
-                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
-
-               ELSE                            !--  case 2 : surface is melting
-                  !
-                  jm_min(ji) = 2
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! first layer of snow equation
-                  ztrid(ji,2,1)  = 0.0
-                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
-               ENDIF
-               !                            !---------------------!
-            ELSE                            ! cells without snow  !
-               !                            !---------------------!
-               !
-               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
-                  !
-                  jm_min(ji) = nlay_s + 1
-                  jm_max(ji) = nlay_i + nlay_s + 1
-
-                  ! surface equation	
-                  ztrid(ji,jm_min(ji),1)  = 0.0
-                  ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
-                  ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
-                  zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
-
-                  ! first layer of ice equation
-                  ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
-                  ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
-                  ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)  
-                  zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
-
-                  ! case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,jm_min(ji),1)    = 0.0
-                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
-                     ztrid(ji,jm_min(ji),3)    =  zkappa_i(ji,0) * 2.0
-                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
-                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,jm_min(ji)+1,3)  = 0.0
-                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
-                  ENDIF
-
-               ELSE                            !--  case 2 : surface is melting
-
-                  jm_min(ji)    =  nlay_s + 2
-                  jm_max(ji)    =  nlay_i + nlay_s + 1
-
-                  ! first layer of ice equation
-                  ztrid(ji,jm_min(ji),1)  = 0.0
-                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
-                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
-                     &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
-
-                  ! case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,jm_min(ji),1)  = 0.0
-                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,jm_min(ji),3)  = 0.0
-                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
-                        &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
-                  ENDIF
-
-               ENDIF
-            ENDIF
-            !
-            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
-            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
-            !
-            END DO
-            !
-            !------------------------------
-            ! 8) tridiagonal system solving
-            !------------------------------
-            ! Solve the tridiagonal system with Gauss elimination method.
-            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-            jm_maxt = 0
-            jm_mint = nlay_i+5
-            DO ji = 1, npti
-               jm_mint = MIN(jm_min(ji),jm_mint)
-               jm_maxt = MAX(jm_max(ji),jm_maxt)
-            END DO
-
-            DO jk = jm_mint+1, jm_maxt
-               DO ji = 1, npti
-                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
-                  zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
-                  zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
-               END DO
-            END DO
-
-            ! ice temperatures
-            DO ji = 1, npti
-               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-            END DO
-
-            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-               DO ji = 1, npti
-                  jk = jm - nlay_s - 1
-                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
-               END DO
-            END DO
-
-            DO ji = 1, npti
-               ! snow temperatures      
-               IF( h_s_1d(ji) > 0._wp ) THEN
-                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
-               ENDIF
-               ! surface temperature
-               ztsub(ji) = t_su_1d(ji)
-               IF( t_su_1d(ji) < rt0 ) THEN
-                  t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                     &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
-               ENDIF
-            END DO
-            !
-            !--------------------------------------------------------------
-            ! 9) Has the scheme converged ?, end of the iterative procedure
-            !--------------------------------------------------------------
-            ! check that nowhere it has started to melt
-            ! zdti_max is a measure of error, it has to be under zdti_bnd
-            zdti_max = 0._wp
-            DO ji = 1, npti
-               t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
-               zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
-            END DO
-
-            DO jk = 1, nlay_s
-               DO ji = 1, npti
-                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
-                  zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
-               END DO
-            END DO
-
-            DO jk = 1, nlay_i
-               DO ji = 1, npti
-                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
-                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
-                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
-               END DO
-            END DO
-
-            ! Compute spatial maximum over all errors
-            ! note that this could be optimized substantially by iterating only the non-converging points
-            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
-         !
-         !----------------------------------------!
-         !                                        !
-         !       JULES IF (RCV MODE)              !
-         !                                        !
-         !----------------------------------------!
-         !
-
-         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode
-         
-            !
-            ! In RCV mode, we use a modified BL99 solver 
-            ! with conduction flux (qcn_ice) as forcing term
-            !
-            !----------------------------
-            ! 7) tridiagonal system terms
-            !----------------------------
-            !!layer denotes the number of the layer in the snow or in the ice
-            !!jm denotes the reference number of the equation in the tridiagonal
-            !!system, terms of tridiagonal system are indexed as following :
-            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
-            !!ice interior terms (top equation has the same form as the others)
-            ztrid   (1:npti,:,:) = 0._wp
-            zindterm(1:npti,:)   = 0._wp
-            zindtbis(1:npti,:)   = 0._wp
-            zdiagbis(1:npti,:)   = 0._wp
-
-            DO jm = nlay_s + 2, nlay_s + nlay_i 
-               DO ji = 1, npti
-                  jk = jm - nlay_s - 1
-                  ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
-                  ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-                  ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
-                  zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-               END DO
-            ENDDO
-
-            jm =  nlay_s + nlay_i + 1
-            DO ji = 1, npti
-               !!ice bottom term
-               ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
-               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
-               ztrid(ji,jm,3)  = 0.0
-               zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-                  &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
-            ENDDO
-
-
-            DO ji = 1, npti
-               !                              !---------------------!
-               IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
-                  !                           !---------------------!
-                  ! snow interior terms (bottom equation has the same form as the others)
-                  DO jm = 3, nlay_s + 1
-                     jk = jm - 1
-                     ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-                     ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                     ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                     zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
-                  END DO
-                  
-                  ! case of only one layer in the ice (ice equation is altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,nlay_s+2,3)  = 0.0
-                     zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
-                  ENDIF
-                  
-                  jm_min(ji) = 2
-                  jm_max(ji) = nlay_i + nlay_s + 1
-                  
-                  ! first layer of snow equation
-                  ztrid(ji,2,1)  = 0.0
-                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1)
-                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
-                  
-                  !                            !---------------------!
-               ELSE                            ! cells without snow  !
-                  !                            !---------------------!
-                  jm_min(ji)    =  nlay_s + 2
-                  jm_max(ji)    =  nlay_i + nlay_s + 1
-                  
-                  ! first layer of ice equation
-                  ztrid(ji,jm_min(ji),1)  = 0.0
-                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)  
-                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
-                  
-                  ! case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,jm_min(ji),1)  = 0.0
-                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)
-                     ztrid(ji,jm_min(ji),3)  = 0.0
-                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)  &
-                        &                    + qcn_ice_1d(ji) )
-                  ENDIF
-                  
-               ENDIF
-               !
-               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
-               zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
-               !
-            END DO
-            !
-            !------------------------------
-            ! 8) tridiagonal system solving
-            !------------------------------
-            ! Solve the tridiagonal system with Gauss elimination method.
-            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-            jm_maxt = 0
-            jm_mint = nlay_i+5
-            DO ji = 1, npti
-               jm_mint = MIN(jm_min(ji),jm_mint)
-               jm_maxt = MAX(jm_max(ji),jm_maxt)
-            END DO
-            
-            DO jk = jm_mint+1, jm_maxt
-               DO ji = 1, npti
-                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
-                  zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
-                  zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
-               END DO
-            END DO
-            
-            ! ice temperatures
-           DO ji = 1, npti
-               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-            END DO
-
-            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-               DO ji = 1, npti
-                  jk = jm - nlay_s - 1
-                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
-               END DO
-            END DO
-            
-            ! snow temperatures      
-            DO ji = 1, npti
-               IF( h_s_1d(ji) > 0._wp ) THEN
-                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
-               ENDIF
-            END DO
-            !
-            !--------------------------------------------------------------
-            ! 9) Has the scheme converged ?, end of the iterative procedure
-            !--------------------------------------------------------------
-            ! check that nowhere it has started to melt
-            ! zdti_max is a measure of error, it has to be under zdti_bnd
-            zdti_max = 0._wp
-
-            DO jk = 1, nlay_s
-               DO ji = 1, npti
-                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
-                  zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
-               END DO
-            END DO
-            
-            DO jk = 1, nlay_i
-               DO ji = 1, npti
-                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
-                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
-                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
-               END DO
-            END DO
-
-            ! Compute spatial maximum over all errors
-            ! note that this could be optimized substantially by iterating only the non-converging points
-            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice )
-
-         ENDIF ! k_jules
-         
-      END DO  ! End of the do while iterative procedure
-      
-      IF( ln_icectl .AND. lwp ) THEN
-         WRITE(numout,*) ' zdti_max : ', zdti_max
-         WRITE(numout,*) ' iconv    : ', iconv
-      ENDIF
-      
-      !
-      !-----------------------------
-      ! 10) Fluxes at the interfaces
-      !-----------------------------
-      !
-      ! --- update conduction fluxes
-      !
-      DO ji = 1, npti
-         !                                ! surface ice conduction flux
-         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
-            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
-         !                                ! bottom ice conduction flux
-         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
-      END DO
-      
-      !
-      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
-      !
-      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  ! OFF or SND mode
-         DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
-         END DO
-      ELSE        ! RCV mode
-         DO ji = 1, npti
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji) 
-         END DO
-      ENDIF
-      
-      !
-      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
-      !
-
-      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF
-         
-         CALL ice_thd_enmelt       
-         
-         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-         DO ji = 1, npti
-            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
-               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
-            
-            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN
-               
-               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-                  zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
-               ELSE                          ! case T_su = 0degC
-                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
-               ENDIF
-               
-            ELSE ! RCV CASE
-            
-               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
-            
-            ENDIF
-            !
-            ! total heat sink to be sent to the ocean
-            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
-            !
-            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
-            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
-            !
-         END DO
-         !
-         ! --- SIMIP diagnostics
-         !
-         DO ji = 1, npti
-            !--- Conduction fluxes (positive downwards)
-            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
-   
-            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
-               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
-            ELSE
-               t_si_1d(ji) = t_su_1d(ji)
-            ENDIF
-         END DO
-         !
-      ENDIF
-      !
-      !---------------------------------------------------------------------------------------
-      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
-      !---------------------------------------------------------------------------------------
-      ! effective conductivity and 1st layer temperature (needed by Met Office)
-      DO ji = 1, npti
-         IF( h_s_1d(ji) > 0.1_wp ) THEN 
-            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
-         ELSE
-            IF( h_i_1d(ji) > 0.1_wp ) THEN
-               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
-            ELSE
-               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp
-            ENDIF
-         ENDIF
-         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
-      END DO
-      !
-      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode
-         ! Restore temperatures to their initial values
-         t_s_1d    (1:npti,:) = ztsold(1:npti,:)
-         t_i_1d    (1:npti,:) = ztiold(1:npti,:)
-         qcn_ice_1d(1:npti)   = fc_su (1:npti)
-      ENDIF
-      !
-   END SUBROUTINE ice_thd_zdf_BL99
-
-
-   SUBROUTINE ice_thd_enmelt
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_thd_enmelt *** 
-      !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
-      !!
-      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts  ! local scalar 
-      !!-------------------------------------------------------------------
-      !
-      DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = 1, npti
-            ztmelts      = - tmut  * sz_i_1d(ji,jk)
-            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
-                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
-            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
-               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
-               &                    - rcp  *   ztmelts )
-         END DO
-      END DO
-      DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = 1, npti
-            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
-         END DO
-      END DO
-      !
-   END SUBROUTINE ice_thd_enmelt
-
-
    SUBROUTINE ice_thd_zdf_init
       !!-----------------------------------------------------------------------
@@ -976,16 +117,14 @@
          WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
          WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-     ENDIF
-
+      ENDIF
       !
       IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
          CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' )
       ENDIF
-
       !                             !== set the choice of ice vertical thermodynamic formulation ==!
       ioptio = 0 
-      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties)
-      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF
-      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' )
+      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF   ! BL99 thermodynamics (linear liquidus + constant thermal properties)
+!!    IF( ln_zdf_XXXX ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_XXXX   ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' )
       !
    END SUBROUTINE ice_thd_zdf_init
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf_bl99.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf_bl99.F90	(revision 8984)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf_bl99.F90	(revision 8984)
@@ -0,0 +1,844 @@
+MODULE icethd_zdf_BL99
+   !!======================================================================
+   !!                       ***  MODULE icethd_zdf_BL99 ***
+   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 
+   !!======================================================================
+   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
+   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
+   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
+   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       ESIM sea-ice model
+   !!----------------------------------------------------------------------
+   !!  ice_thd_zdf_BL99 : vertical diffusion computation
+   !!----------------------------------------------------------------------
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants (ocean directory) 
+   USE ice            ! sea-ice: variables
+   USE ice1D          ! sea-ice: thermodynamics variables
+   USE icevar         ! sea-ice: operations
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_thd_zdf_BL99   ! called by icethd_zdf
+
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
+   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_thd_zdf_BL99( k_jules )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_thd_zdf_BL99  ***
+      !!
+      !! ** Purpose : computes the time evolution of snow and sea-ice temperature
+      !!              profiles, using the original Bitz and Lipscomb (1999) algorithm
+      !!
+      !! ** Method  : solves the heat equation diffusion with a Neumann boundary
+      !!              condition at the surface and a Dirichlet one at the bottom. 
+      !!              Solar radiation is partially absorbed into the ice.
+      !!              The specific heat and thermal conductivities depend on ice 
+      !!              salinity and temperature to take into account brine pocket   
+      !!              melting. The numerical scheme is an iterative Crank-Nicolson
+      !!              on a non-uniform multilayer grid in the ice and snow system.
+      !!
+      !!           The successive steps of this routine are
+      !!           1.  initialization of ice-snow layers thicknesses
+      !!           2.  Internal absorbed and transmitted radiation
+      !!           Then iterative procedure begins
+      !!           3.  Thermal conductivity
+      !!           4.  Kappa factors
+      !!           5.  specific heat in the ice
+      !!           6.  eta factors
+      !!           7.  surface flux computation
+      !!           8.  tridiagonal system terms
+      !!           9.  solving the tridiagonal system with Gauss elimination
+      !!           Iterative procedure ends according to a criterion on evolution
+      !!           of temperature
+      !!           10. Fluxes at the interfaces
+      !!
+      !! ** Inputs / Ouputs : (global commons)
+      !!           surface temperature              : t_su_1d
+      !!           ice/snow temperatures            : t_i_1d, t_s_1d
+      !!           ice salinities                   : sz_i_1d
+      !!           number of layers in the ice/snow : nlay_i, nlay_s
+      !!           total ice/snow thickness         : h_i_1d, h_s_1d
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE)
+      !
+      INTEGER ::   ji, jk         ! spatial loop index
+      INTEGER ::   jm             ! current reference number of equation
+      INTEGER ::   jm_mint, jm_maxt
+      INTEGER ::   iconv          ! number of iterations in iterative procedure
+      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
+      !
+      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
+      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
+      !
+      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
+      REAL(wp) ::   zg1       =  2._wp        !
+      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
+      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
+      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
+      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
+      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
+      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
+      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
+      REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
+      REAL(wp) ::   zcpi                      ! Ice specific heat
+      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
+      REAL(wp) ::   zfac                      ! dummy factor
+      !
+      REAL(wp), DIMENSION(jpij) ::   isnow        ! switch for presence (1) or absence (0) of snow
+      REAL(wp), DIMENSION(jpij) ::   ztsub        ! surface temperature at previous iteration
+      REAL(wp), DIMENSION(jpij) ::   zh_i, z1_h_i ! ice layer thickness
+      REAL(wp), DIMENSION(jpij) ::   zh_s, z1_h_s ! snow layer thickness
+      REAL(wp), DIMENSION(jpij) ::   zqns_ice_b   ! solar radiation absorbed at the surface
+      REAL(wp), DIMENSION(jpij) ::   zfnet        ! surface flux function
+      REAL(wp), DIMENSION(jpij) ::   zdqns_ice_b  ! derivative of the surface flux function
+      !
+      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
+      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
+      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
+      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
+      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
+      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
+      !
+      ! Mono-category
+      REAL(wp) ::   zepsilon   ! determines thres. above which computation of G(h) is done
+      REAL(wp) ::   zhe        ! dummy factor
+      REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity
+      !!------------------------------------------------------------------     
+
+      ! --- diag error on heat diffusion - PART 1 --- !
+      DO ji = 1, npti
+         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
+            &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
+      END DO
+
+      !------------------
+      ! 1) Initialization
+      !------------------
+      DO ji = 1, npti
+         isnow(ji) = 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
+         ! layer thickness
+         zh_i(ji) = h_i_1d(ji) * r1_nlay_i
+         zh_s(ji) = h_s_1d(ji) * r1_nlay_s
+      END DO
+      !
+      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
+      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
+      END WHERE
+      !
+      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
+      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
+      END WHERE
+      !
+      ! Store initial temperatures and non solar heat fluxes
+      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
+         !
+         ztsub      (1:npti) = t_su_1d(1:npti)                          ! surface temperature at iteration n-1
+         ztsuold    (1:npti) = t_su_1d(1:npti)                          ! surface temperature initial value
+         t_su_1d    (1:npti) = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not
+         zdqns_ice_b(1:npti) = dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux 
+         zqns_ice_b (1:npti) = qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value
+         !
+      ENDIF
+      !
+      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
+      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
+
+      !-------------
+      ! 2) Radiation
+      !-------------
+      ! --- Transmission/absorption of solar radiation in the ice --- !
+      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti)
+      DO jk = 1, nlay_s
+         DO ji = 1, npti
+            !                             ! radiation transmitted below the layer-th snow layer
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
+            !                             ! radiation absorbed by the layer-th snow layer
+            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
+         END DO
+      END DO
+      !
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) )
+      DO jk = 1, nlay_i 
+         DO ji = 1, npti
+            !                             ! radiation transmitted below the layer-th ice layer
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
+            !                             ! radiation absorbed by the layer-th ice layer
+            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
+         END DO
+      END DO
+      !
+      ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
+      !
+      iconv    = 0          ! number of iterations
+      zdti_max = 1000._wp   ! maximal value of error on all points
+      !
+      !                                                          !============================!
+      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
+         !                                                       !============================!
+         iconv = iconv + 1
+         !
+         ztib(1:npti,:) = t_i_1d(1:npti,:)
+         ztsb(1:npti,:) = t_s_1d(1:npti,:)
+         !
+         !--------------------------------
+         ! 3) Sea ice thermal conductivity
+         !--------------------------------
+         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
+            !
+            DO ji = 1, npti
+               ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
+               ztcond_i(ji,nlay_i) = rcdic + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
+            END DO
+            DO jk = 1, nlay_i-1
+               DO ji = 1, npti
+                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
+                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+               END DO
+            END DO
+            !
+         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
+            !
+            DO ji = 1, npti
+               ztcond_i(ji,0)      = rcdic + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
+                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
+               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+            END DO
+            DO jk = 1, nlay_i-1
+               DO ji = 1, npti
+                  ztcond_i(ji,jk) = rcdic + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
+                     &                     MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
+                     &                    - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
+               END DO
+            END DO
+            !
+         ENDIF
+         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )        
+         !
+         !--- G(he) : enhancement of thermal conductivity in mono-category case
+         ! Computation of effective thermal conductivity G(h)
+         ! Used in mono-category case only to simulate an ITD implicitly
+         ! Fichefet and Morales Maqueda, JGR 1997
+         zghe(1:npti) = 1._wp
+         !
+         SELECT CASE ( nn_monocat )
+         !
+         CASE ( 1 , 3 )
+            !
+            zepsilon = 0.1_wp
+            DO ji = 1, npti
+               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
+               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
+               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
+                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
+            END DO
+            !
+         END SELECT
+         !
+         !-----------------
+         ! 4) kappa factors
+         !-----------------
+         !--- Snow
+         DO jk = 0, nlay_s-1
+            DO ji = 1, npti
+               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
+            END DO
+         END DO
+         DO ji = 1, npti   ! Snow-ice interface
+            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
+            IF( zfac > epsi10 ) THEN
+               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
+            ELSE
+               zkappa_s(ji,nlay_s) = 0._wp
+            ENDIF
+         END DO
+
+         !--- Ice
+         DO jk = 0, nlay_i
+            DO ji = 1, npti
+               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
+            END DO
+         END DO
+         DO ji = 1, npti   ! Snow-ice interface
+            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+         END DO
+         !
+         !--------------------------------------
+         ! 5) Sea ice specific heat, eta factors
+         !--------------------------------------
+         DO jk = 1, nlay_i
+            DO ji = 1, npti
+               zcpi = cpic + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
+               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+            END DO
+         END DO
+
+         DO jk = 1, nlay_s
+            DO ji = 1, npti
+               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
+            END DO
+         END DO
+         !
+         !----------------------------------------!
+         !                                        !
+         !   JULES COUPLING IS OFF OR EMULATED    !
+         !                                        !
+         !----------------------------------------!
+         !
+         IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
+            !
+            ! ==> The original BL99 temperature computation is used
+            !       (with qsr_ice, qns_ice and dqns_ice as inputs)
+            !
+            !----------------------------
+            ! 6) surface flux computation
+            !----------------------------
+            ! update of the non solar flux according to the update in T_su
+            DO ji = 1, npti
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
+            END DO
+
+            DO ji = 1, npti
+               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
+            END DO
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            ! layer denotes the number of the layer in the snow or in the ice
+            ! jm denotes the reference number of the equation in the tridiagonal
+            ! system, terms of tridiagonal system are indexed as following :
+            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            ! ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
+                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
+                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+               END DO
+            END DO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
+               ! ice bottom term
+               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
+               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
+               ztrid   (ji,jm,3) = 0._wp
+               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
+            END DO
+
+            DO ji = 1, npti
+               !                               !---------------------!
+               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
+                  !                            !---------------------!
+                  ! snow interior terms (bottom equation has the same form as the others)
+                  DO jm = 3, nlay_s + 1
+                     jk = jm - 1
+                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
+                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
+                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+                  END DO
+                  
+                  ! case of only one layer in the ice (ice equation is altered)
+                  IF( nlay_i == 1 ) THEN
+                     ztrid   (ji,nlay_s+2,3) = 0._wp
+                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+                  ENDIF
+                  
+                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
+                     
+                     jm_min(ji) = 1
+                     jm_max(ji) = nlay_i + nlay_s + 1
+                     
+                     ! surface equation
+                     ztrid   (ji,1,1) = 0._wp
+                     ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
+                     ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0)
+                     zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+                     
+                     ! first layer of snow equation
+                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
+                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1)
+                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
+                     
+                  ELSE                            !--  case 2 : surface is melting
+                     !
+                     jm_min(ji) = 2
+                     jm_max(ji) = nlay_i + nlay_s + 1
+                     
+                     ! first layer of snow equation
+                     ztrid   (ji,2,1) = 0._wp
+                     ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                     ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
+                     zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+                  ENDIF
+                  !                            !---------------------!
+               ELSE                            ! cells without snow  !
+                  !                            !---------------------!
+                  !
+                  IF( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
+                     !
+                     jm_min(ji) = nlay_s + 1
+                     jm_max(ji) = nlay_i + nlay_s + 1
+                     
+                     ! surface equation	
+                     ztrid   (ji,jm_min(ji),1) = 0._wp
+                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1    
+                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
+                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+                     
+                     ! first layer of ice equation
+                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
+                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
+                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)  
+                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
+                     
+                     ! case of only one layer in the ice (surface & ice equations are altered)
+                     IF( nlay_i == 1 ) THEN
+                        ztrid   (ji,jm_min(ji),1)   = 0._wp
+                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
+                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
+                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
+                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
+                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
+                        zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))
+                     ENDIF
+                     
+                  ELSE                            !--  case 2 : surface is melting
+                     
+                     jm_min(ji) = nlay_s + 2
+                     jm_max(ji) = nlay_i + nlay_s + 1
+                     
+                     ! first layer of ice equation
+                     ztrid   (ji,jm_min(ji),1) = 0._wp
+                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
+                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
+                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 
+                     
+                     ! case of only one layer in the ice (surface & ice equations are altered)
+                     IF( nlay_i == 1 ) THEN
+                        ztrid   (ji,jm_min(ji),1) = 0._wp
+                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
+                        ztrid   (ji,jm_min(ji),3) = 0._wp
+                        zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) &
+                           &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp
+                     ENDIF
+                     
+                  ENDIF
+               ENDIF
+               !
+               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
+               !
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
+               jm_mint = MIN(jm_min(ji),jm_mint)
+               jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+
+            DO jk = jm_mint+1, jm_maxt
+               DO ji = 1, npti
+                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
+                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
+                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
+               END DO
+            END DO
+
+            ! ice temperatures
+            DO ji = 1, npti
+               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+               END DO
+            END DO
+
+            DO ji = 1, npti
+               ! snow temperatures      
+               IF( h_s_1d(ji) > 0._wp ) THEN
+                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
+               ENDIF
+               ! surface temperature
+               ztsub(ji) = t_su_1d(ji)
+               IF( t_su_1d(ji) < rt0 ) THEN
+                  t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
+                     &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+               ENDIF
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+            DO ji = 1, npti
+               t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
+               zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
+            END DO
+
+            DO jk = 1, nlay_s
+               DO ji = 1, npti
+                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+                  zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+               END DO
+            END DO
+
+            DO jk = 1, nlay_i
+               DO ji = 1, npti
+                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+               END DO
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
+         !
+         !----------------------------------------!
+         !                                        !
+         !      JULES COUPLING IS ACTIVE          !
+         !                                        !
+         !----------------------------------------!
+         !
+         ELSEIF( k_jules == np_jules_ACTIVE ) THEN
+            !
+            ! ==> we use a modified BL99 solver with conduction flux (qcn_ice) as forcing term
+            !
+            !----------------------------
+            ! 7) tridiagonal system terms
+            !----------------------------
+            ! layer denotes the number of the layer in the snow or in the ice
+            ! jm denotes the reference number of the equation in the tridiagonal
+            ! system, terms of tridiagonal system are indexed as following :
+            ! 1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+            ! ice interior terms (top equation has the same form as the others)
+            ztrid   (1:npti,:,:) = 0._wp
+            zindterm(1:npti,:)   = 0._wp
+            zindtbis(1:npti,:)   = 0._wp
+            zdiagbis(1:npti,:)   = 0._wp
+
+            DO jm = nlay_s + 2, nlay_s + nlay_i 
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
+                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
+                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+               END DO
+            ENDDO
+
+            jm =  nlay_s + nlay_i + 1
+            DO ji = 1, npti
+               ! ice bottom term
+               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
+               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
+               ztrid   (ji,jm,3) = 0._wp
+               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
+            ENDDO
+
+            DO ji = 1, npti
+               !                               !---------------------!
+               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
+                  !                            !---------------------!
+                  ! snow interior terms (bottom equation has the same form as the others)
+                  DO jm = 3, nlay_s + 1
+                     jk = jm - 1
+                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
+                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
+                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+                  END DO
+                  
+                  ! case of only one layer in the ice (ice equation is altered)
+                  IF ( nlay_i == 1 ) THEN
+                     ztrid   (ji,nlay_s+2,3) = 0._wp
+                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+                  ENDIF
+                  
+                  jm_min(ji) = 2
+                  jm_max(ji) = nlay_i + nlay_s + 1
+                  
+                  ! first layer of snow equation
+                  ztrid   (ji,2,1) = 0._wp
+                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
+                  ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1) 
+                  zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) ) 
+                  
+                  !                            !---------------------!
+               ELSE                            ! cells without snow  !
+                  !                            !---------------------!
+                  jm_min(ji) = nlay_s + 2
+                  jm_max(ji) = nlay_i + nlay_s + 1
+                  
+                  ! first layer of ice equation
+                  ztrid   (ji,jm_min(ji),1) = 0._wp
+                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)  
+                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
+                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+                  
+                  ! case of only one layer in the ice (surface & ice equations are altered)
+                  IF( nlay_i == 1 ) THEN
+                     ztrid   (ji,jm_min(ji),1) = 0._wp
+                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
+                     ztrid   (ji,jm_min(ji),3) = 0._wp
+                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
+                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
+                  ENDIF
+                  
+               ENDIF
+               !
+               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
+               !
+            END DO
+            !
+            !------------------------------
+            ! 8) tridiagonal system solving
+            !------------------------------
+            ! Solve the tridiagonal system with Gauss elimination method.
+            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+            jm_maxt = 0
+            jm_mint = nlay_i+5
+            DO ji = 1, npti
+               jm_mint = MIN(jm_min(ji),jm_mint)
+               jm_maxt = MAX(jm_max(ji),jm_maxt)
+            END DO
+            
+            DO jk = jm_mint+1, jm_maxt
+               DO ji = 1, npti
+                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
+                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
+                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
+               END DO
+            END DO
+            
+            ! ice temperatures
+           DO ji = 1, npti
+               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+            END DO
+
+            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+               DO ji = 1, npti
+                  jk = jm - nlay_s - 1
+                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+               END DO
+            END DO
+            
+            ! snow temperatures      
+            DO ji = 1, npti
+               IF( h_s_1d(ji) > 0._wp ) THEN
+                  t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1)
+               ENDIF
+            END DO
+            !
+            !--------------------------------------------------------------
+            ! 9) Has the scheme converged?, end of the iterative procedure
+            !--------------------------------------------------------------
+            ! check that nowhere it has started to melt
+            ! zdti_max is a measure of error, it has to be under zdti_bnd
+            zdti_max = 0._wp
+
+            DO jk = 1, nlay_s
+               DO ji = 1, npti
+                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+                  zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+               END DO
+            END DO
+            
+            DO jk = 1, nlay_i
+               DO ji = 1, npti
+                  ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
+                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+               END DO
+            END DO
+
+            ! Compute spatial maximum over all errors
+            ! note that this could be optimized substantially by iterating only the non-converging points
+            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice )
+
+         ENDIF ! k_jules
+         
+      END DO  ! End of the do while iterative procedure
+      
+      IF( ln_icectl .AND. lwp ) THEN
+         WRITE(numout,*) ' zdti_max : ', zdti_max
+         WRITE(numout,*) ' iconv    : ', iconv
+      ENDIF
+      
+      !
+      !-----------------------------
+      ! 10) Fluxes at the interfaces
+      !-----------------------------
+      !
+      ! --- update conduction fluxes
+      !
+      DO ji = 1, npti
+         !                                ! surface ice conduction flux
+         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
+            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
+         !                                ! bottom ice conduction flux
+         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * zg1 * ( t_bo_1d(ji) - t_i_1d(ji,nlay_i) )
+      END DO
+      
+      !
+      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- !
+      !
+      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_EMULE ) THEN
+         !
+         DO ji = 1, npti
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 
+         END DO
+         !
+      ELSEIF( k_jules == np_jules_ACTIVE ) THEN
+         !
+         DO ji = 1, npti
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji)      - qcn_ice_1d(ji) ) * a_i_1d(ji) 
+         END DO
+         !
+      ENDIF
+      
+      !
+      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- !
+      !
+      IF( k_jules == np_jules_OFF .OR. k_jules == np_jules_ACTIVE ) THEN
+         
+         CALL ice_var_enthalpy       
+         
+         ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
+         DO ji = 1, npti
+            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
+               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
+            
+            IF( k_jules == np_jules_OFF ) THEN
+               
+               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
+                  zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
+               ELSE                          ! case T_su = 0degC
+                  zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice )*a_i_1d(ji)
+               ENDIF
+               
+            ELSEIF( k_jules == np_jules_ACTIVE ) THEN
+            
+               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+            
+            ENDIF
+            !
+            ! total heat sink to be sent to the ocean
+            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+            !
+            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2   
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
+            !
+         END DO
+         !
+         ! --- SIMIP diagnostics
+         !
+         DO ji = 1, npti
+            !--- Conduction fluxes (positive downwards)
+            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su  (ji) * a_i_1d(ji) / at_i_1d(ji)
+   
+            !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
+               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
+                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
+            ELSE
+               t_si_1d(ji) = t_su_1d(ji)
+            ENDIF
+         END DO
+         !
+      ENDIF
+      !
+      !---------------------------------------------------------------------------------------
+      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes
+      !---------------------------------------------------------------------------------------
+      ! effective conductivity and 1st layer temperature (needed by Met Office)
+      DO ji = 1, npti
+         IF( h_s_1d(ji) > 0.1_wp ) THEN 
+            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
+         ELSE
+            IF( h_i_1d(ji) > 0.1_wp ) THEN
+               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
+            ELSE
+               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) / 0.1_wp
+            ENDIF
+         ENDIF
+         t1_ice_1d(ji) = isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1)
+      END DO
+      !
+      IF( k_jules == np_jules_EMULE ) THEN
+         ! Restore temperatures to their initial values
+         t_s_1d    (1:npti,:) = ztsold(1:npti,:)
+         t_i_1d    (1:npti,:) = ztiold(1:npti,:)
+         qcn_ice_1d(1:npti)   = fc_su (1:npti)
+      ENDIF
+      !
+   END SUBROUTINE ice_thd_zdf_BL99
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option       Dummy Module             No ESIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE icethd_zdf_BL99
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8983)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8984)
@@ -47,4 +47,5 @@
    !!   ice_var_itd2      : convert N-cat to jpl-cat
    !!   ice_var_bv        : brine volume
+   !!   ice_var_enthalpy  : compute ice and snow enthalpies from temperature
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -70,4 +71,5 @@
    PUBLIC   ice_var_itd2
    PUBLIC   ice_var_bv           
+   PUBLIC   ice_var_enthalpy           
 
    !!----------------------------------------------------------------------
@@ -396,5 +398,5 @@
          !                                      ! Slope of the linear profile 
          WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
-         ELSEWHERE                           ;   z_slope_s(1:npti) = 0._wp
+         ELSEWHERE                          ;   z_slope_s(1:npti) = 0._wp
          END WHERE
          
@@ -822,4 +824,34 @@
 
 
+   SUBROUTINE ice_var_enthalpy
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_enthalpy *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = 1, npti
+            ztmelts      = - tmut  * sz_i_1d(ji,jk)
+            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
+                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
+            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
+               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
+               &                    - rcp  *   ztmelts )
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = 1, npti
+            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE ice_var_enthalpy
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8983)
+++ /branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8984)
@@ -46,7 +46,6 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_lim3
-   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s
    USE icethd_dh      ! for CALL ice_thd_snwblow
-   USE icethd_zdf,ONLY:   rn_cnd_s  ! for blk_ice_qcn
 #endif
    USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
