source: branches/dev_r4028_CNRS_LIM3_MV2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4506

MODULE limthd_dh
   !!======================================================================
   !!                       ***  MODULE limthd_dh ***
   !!  LIM-3 :   thermodynamic growth and decay of the ice 
   !!======================================================================
   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
   !!                 ! 2005-06 (M. Vancoppenolle) 3D version 
   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice
   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes 
   !!----------------------------------------------------------------------
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3'                                      LIM3 sea-ice model
   !!----------------------------------------------------------------------
   !!   lim_thd_dh    : vertical accr./abl. and lateral ablation of sea ice
   !!----------------------------------------------------------------------
   USE par_oce        ! ocean parameters
   USE phycst         ! physical constants (OCE directory) 
   USE sbc_oce        ! Surface boundary condition: ocean fields
   USE ice            ! LIM variables
   USE par_ice        ! LIM parameters
   USE thd_ice        ! LIM thermodynamics
   USE in_out_manager ! I/O manager
   USE lib_mpp        ! MPP library
   USE wrk_nemo       ! work arrays
   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  

   IMPLICIT NONE
   PRIVATE

   PUBLIC   lim_thd_dh   ! called by lim_thd

   REAL(wp) ::   epsi20 = 1.e-20   ! constant values
   REAL(wp) ::   epsi10 = 1.e-10   !
   REAL(wp) ::   epsi13 = 1.e-13   !
   REAL(wp) ::   zzero  = 0._wp    !
   REAL(wp) ::   zone   = 1._wp    !

   !!----------------------------------------------------------------------
   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
   !! $Id$
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE lim_thd_dh( kideb, kiut, jl )
      !!------------------------------------------------------------------
      !!                ***  ROUTINE lim_thd_dh  ***
      !!
      !! ** Purpose :   determines variations of ice and snow thicknesses.
      !!
      !! ** Method  :   Ice/Snow surface melting arises from imbalance in surface fluxes
      !!              Bottom accretion/ablation arises from flux budget
      !!              Snow thickness can increase by precipitation and decrease by sublimation
      !!              If snow load excesses Archmiede limit, snow-ice is formed by
      !!              the flooding of sea-water in the snow
      !!
      !!                 1) Compute available flux of heat for surface ablation
      !!                 2) Compute snow and sea ice enthalpies
      !!                 3) Surface ablation and sublimation
      !!                 4) Bottom accretion/ablation
      !!                 5) Case of Total ablation
      !!                 6) Snow ice formation
      !!
      !! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
      !!              Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
      !!              Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let. 
      !!              Vancoppenolle et al.,2009, Ocean Modelling
      !!------------------------------------------------------------------
      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
      !! 
      INTEGER  ::   ji , jk        ! dummy loop indices
      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji
      INTEGER  ::   i_use          ! debugging flag
      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow
      INTEGER  ::   isnowic        ! snow ice formation not
      INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not
      INTEGER  ::   iter

      REAL(wp) ::   zzfmass_i, zihgnew                     ! local scalar
      REAL(wp) ::   zzfmass_s, zhsnew, ztmelts             ! local scalar
      REAL(wp) ::   zhn, zdhcf, zdhbf, zhni, zhnfi, zihg   !
      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       !
      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
      REAL(wp) ::   zs_snic  ! snow-ice salinity
      REAL(wp) ::   zswi1        ! switch for computation of bottom salinity
      REAL(wp) ::   zswi12       ! switch for computation of bottom salinity
      REAL(wp) ::   zswi2        ! switch for computation of bottom salinity
      REAL(wp) ::   zgrr         ! bottom growth rate
      REAL(wp) ::   zt_i_new     ! bottom formation temperature

      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
      REAL(wp) ::   zsstK        ! SST in Kelvin

      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! melting point
      REAL(wp), POINTER, DIMENSION(:) ::   zhsold      ! old snow thickness
      REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow
      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_su   ! incoming, remaining surface energy
      REAL(wp), POINTER, DIMENSION(:) ::   zqfont_bo   ! incoming, bottom energy
      REAL(wp), POINTER, DIMENSION(:) ::   z_f_surf    ! surface heat for ablation
      REAL(wp), POINTER, DIMENSION(:) ::   zhgnew      ! new ice thickness
      REAL(wp), POINTER, DIMENSION(:) ::   zfmass_i    ! 

      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt 
      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation 
      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation

      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah

      ! Pathological cases
      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_init   ! total incoming heat for ice melt
      REAL(wp), POINTER, DIMENSION(:) ::   zfdt_final  ! total remaing heat for ice melt
      REAL(wp), POINTER, DIMENSION(:) ::   zqt_i       ! total ice heat content
      REAL(wp), POINTER, DIMENSION(:) ::   zqt_s       ! total snow heat content
      REAL(wp), POINTER, DIMENSION(:) ::   zqt_dummy   ! dummy heat content

      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content

      ! mass and salt flux (clem)
      REAL(wp) :: zdvres, zdvsur, zdvbot, zswitch_sal
      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume...

      ! Heat conservation 
      INTEGER  ::   num_iter_max, numce_dh
      REAL(wp) ::   meance_dh
      REAL(wp) ::   zinda 
      REAL(wp), POINTER, DIMENSION(:) ::   zinnermelt
      REAL(wp), POINTER, DIMENSION(:) ::   zfbase, zdq_i
      !!------------------------------------------------------------------

      ! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
      SELECT CASE( num_sal )                       ! varying salinity or not
         CASE( 1, 3, 4 ) ;   zswitch_sal = 0       ! prescribed salinity profile
         CASE( 2 )       ;   zswitch_sal = 1       ! varying salinity profile
      END SELECT

      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i )
      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay )

      CALL wrk_alloc( jpij, zviold, zvsold ) ! clem
      
      ftotal_fin(:) = 0._wp
      zfdt_init (:) = 0._wp
      zfdt_final(:) = 0._wp

      dh_i_surf (:) = 0._wp
      dh_i_bott (:) = 0._wp
      dh_snowice(:) = 0._wp

      DO ji = kideb, kiut
         old_ht_i_b(ji) = ht_i_b(ji)
         old_ht_s_b(ji) = ht_s_b(ji)
         zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem
         zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem
      END DO
      !
      !------------------------------------------------------------------------------!
      !  1) Calculate available heat for surface ablation                            !
      !------------------------------------------------------------------------------!
      !
      DO ji = kideb, kiut
         isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  )  )
         ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt
         z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
         z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) )  )
         zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice
      END DO ! ji

      zqfont_su  (:) = 0._wp
      zqfont_bo  (:) = 0._wp
      dsm_i_se_1d(:) = 0._wp     
      dsm_i_si_1d(:) = 0._wp   
      !
      !------------------------------------------------------------------------------!
      !  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
      !------------------------------------------------------------------------------!
      !
      DO ji = kideb, kiut     ! Layer thickness
         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
      END DO
      !
      zqt_s(:) = 0._wp        ! Total enthalpy of the snow
      DO jk = 1, nlay_s
         DO ji = kideb, kiut
            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )
         END DO
      END DO
      !
      zqt_i(:) = 0._wp        ! Total enthalpy of the ice
      DO jk = 1, nlay_i
         DO ji = kideb, kiut
            zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i )
            zqt_i(ji)        =  zqt_i(ji) + zzc
            zqt_i_lay(ji,jk) =              zzc
         END DO
      END DO
      !
      !------------------------------------------------------------------------------|
      !  3) Surface ablation and sublimation                                         |
      !------------------------------------------------------------------------------|
      !
      !-------------------------
      ! 3.1 Snow precips / melt
      !-------------------------
      ! Snow accumulation in one thermodynamic time step
      ! snowfall is partitionned between leads and ice
      ! if snow fall was uniform, a fraction (1-at_i) would fall into leads
      ! but because of the winds, more snow falls on leads than on sea ice
      ! and a greater fraction (1-at_i)^beta of the total mass of snow 
      ! (beta < 1) falls in leads.
      ! In reality, beta depends on wind speed, 
      ! and should decrease with increasing wind speed but here, it is 
      ! considered as a constant. an average value is 0.66
      ! Martin Vancoppenolle, December 2006

      ! Snow fall
      DO ji = kideb, kiut
         zcoeff = ( 1.0 - ( 1.0 - at_i_b(ji) )**betas ) / at_i_b(ji) 
         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
      END DO
      zdh_s_mel(:) =  0._wp

      ! Melt of fallen snow
      DO ji = kideb, kiut
         ! tatm_ice is now in K
         zqprec   (ji)   =  rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )  
         zqfont_su(ji)   =  z_f_surf(ji) * rdt_ice
         zdeltah  (ji,1) =  MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
         zqfont_su(ji)   =  MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1)              ) * zqprec(ji)
         zdeltah  (ji,1) =  MAX( - zdh_s_pre(ji) , zdeltah(ji,1) )
         zdh_s_mel(ji)   =  zdh_s_mel(ji) + zdeltah(ji,1)
         ! heat conservation
         qt_s_in(ji,jl)  =  qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)
         zqt_s  (ji)     =  zqt_s  (ji)    + zqprec(ji) * zdh_s_pre(ji)
         zqt_s  (ji)     =  MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 
      END DO


      ! Snow melt due to surface heat imbalance
      DO jk = 1, nlay_s
         DO ji = kideb, kiut
            zdeltah  (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk)
            zqfont_su(ji)    =  MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 
            zdeltah  (ji,jk) =  MAX( zdeltah(ji,jk) , - zh_s(ji) )
            zdh_s_mel(ji)    =  zdh_s_mel(ji) + zdeltah(ji,jk)        ! resulting melt of snow    
         END DO
      END DO

      ! Apply snow melt to snow depth
      DO ji = kideb, kiut
         dh_s_tot(ji)   =  zdh_s_mel(ji) + zdh_s_pre(ji)
         ! Old and new snow depths
         zhsold(ji)     =  ht_s_b(ji)
         zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
         ! If snow is still present zhn = 1, else zhn = 0
         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  )
         ht_s_b(ji)     =  MAX( zzero , zhsnew )
         ! we recompute dh_s_tot (clem) 
         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)
         ! Volume and mass variations of snow
         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) )
         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) )
         !clem rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)
      END DO ! ji

      !--------------------------
      ! 3.2 Surface ice ablation 
      !--------------------------
      DO ji = kideb, kiut 
         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test
         zdq_i    (ji) =  0._wp
      END DO ! ji

      DO jk = 1, nlay_i
         DO ji = kideb, kiut 
            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]

            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K]

            zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]

            zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0

            zfmdt          = - zqfont_su(ji) / zdE                 ! Mass flux to the ocean [kg/m2, >0]

            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0]

            zqfont_su(ji)  = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE )
                                                                   ! Energy remaining in case of melting of the full layer [J/m2, >0]
            zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]

            dh_i_surf(ji)  = dh_i_surf(ji) + zdeltah(ji,jk)        ! Cumulate surface melt

            zfmdt          = - rhoic*zdeltah(ji,jk)                ! Recompute mass flux [kg/m2, >0]

            zQm            = MAX ( zfmdt, 0._wp ) * zEw            ! Energy of the melt water sent to the ocean [J/m2, <0]

            ! Contribution to salt flux 
            sfx_thd_1d(ji)   = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice
            ! Contribution to heat flux
            rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)
            
            ! Conservation test
            zdq_i(ji)      = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice  ! ? still don't know
         END DO
      END DO

      !                                          !-------------------
      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
         !                                       !-------------------
         numce_dh  = 0
         meance_dh = 0._wp
         DO ji = kideb, kiut
            IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN
               numce_dh  = numce_dh + 1
               meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)
            ENDIF
            IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
               WRITE(numout,*) ' ALERTE heat loss for surface melt '
               WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl
               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
               WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji)
               WRITE(numout,*) ' zdq_i        : ', zdq_i(ji)
               WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
               WRITE(numout,*) ' fc_bo_i      : ', fc_bo_i(ji)
               WRITE(numout,*) ' fbif_1d      : ', fbif_1d(ji)
               WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji)
               WRITE(numout,*) ' s_i_new      : ', s_i_new(ji)
               WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij)
            ENDIF
         END DO
         !
         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
         WRITE(numout,*) ' Error report - Category : ', jl
         WRITE(numout,*) ' ~~~~~~~~~~~~ '
         WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh
         WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh
         !
      ENDIF

      !----------------------
      ! 3.3 Snow sublimation
      !----------------------

      DO ji = kideb, kiut
         ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
#if defined key_coupled
         zdh_s_sub(ji)    =  0._wp      ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
#else
         !                              ! forced  mode: snow thickness change due to sublimation
         zdh_s_sub(ji)    =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
#endif
         dh_s_tot (ji)    =  dh_s_tot(ji) + zdh_s_sub(ji)
         zdhcf            =  ht_s_b(ji) + zdh_s_sub(ji) 
         ht_s_b   (ji)    =  MAX( zzero , zdhcf )
         ! we recompute dh_s_tot 
         dh_s_tot (ji)    =  ht_s_b(ji) - zhsold(ji)
         qt_s_in  (ji,jl) =  qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)
      END DO

      zqt_dummy(:) = 0.e0
      DO jk = 1, nlay_s
         DO ji = kideb,kiut
            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation
         END DO
      END DO

      DO jk = 1, nlay_s
         DO ji = kideb, kiut
            ! In case of disparition of the snow, we have to update the snow temperatures
            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) )  )
            t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
            q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
         END DO
      END DO

      !
      !------------------------------------------------------------------------------!
      ! 4) Basal growth / melt                                                       !
      !------------------------------------------------------------------------------!
      !
      !------------------
      ! 4.1 Basal growth 
      !------------------
      ! Basal growth is driven by heat imbalance at the ice-ocean interface,
      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux 
      ! (qlbbqb) and the turbulent ocean flux (fbif). 
      ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 

      ! If salinity varies in time, an iterative procedure is required, because
      ! the involved quantities are inter-dependent.
      ! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
      ! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
      ! -> need for an iterative procedure, which converges quickly

      IF ( num_sal == 2 ) THEN
         num_iter_max = 5
      ELSE
         num_iter_max = 1
      ENDIF

      ! Initialize dh_i_bott
      dh_i_bott(:) = 0.0e0

      ! Iterative procedure
      DO iter = 1, num_iter_max
         DO ji = kideb, kiut

            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN

               ! New bottom ice salinity (Cox & Weeks, JGR88 )
               !--- zswi1  if dh/dt < 2.0e-8
               !--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7 
               !--- zswi2  if dh/dt > 3.6e-7
               zgrr               = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
               zswi2              = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) )
               zswi12             = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
               zswi1              = 1. - zswi2 * zswi12
               zfracs             = MIN ( zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
                  &               + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  , 0.5 )

               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1

               s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
                                  + ( 1. - zswitch_sal ) * sm_i_b(ji) 
               ! New ice growth
               ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)

               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
               
               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
                  &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   &
                  &               + rcp  * ( ztmelts-rtt )          

               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)

               zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)

               dh_i_bott(ji)      = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / ( zdE * rhoic )

               !!! not sure we still need the next line... useful to keep this in memory ? check limthd_ent...
               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)

            ENDIF ! fc_bo_i
         END DO ! ji
      END DO ! iter

      ! Contribution to Energy and Salt Fluxes
      DO ji = kideb, kiut

         zEw            = rcp  * ( t_bo_b(ji) - rt0 )                   ! Specific enthalpy of seawater (J/kg, < 0)

         zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0)

         ! Contribution to energy flux to the ocean (J/m2)
         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zEw         * a_i_b(ji) * zfmdt

         ! Contribution to salt flux  ()
         sfx_thd_1d(ji) = sfx_thd_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice

      END DO

      !----------------
      ! 4.2 Basal melt
      !----------------
      meance_dh = 0._wp
      numce_dh  = 0
      zinnermelt(:) = 0._wp

      DO ji = kideb, kiut
         ! heat convergence at the surface > 0
         IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0._wp  ) THEN
            s_i_new(ji) = s_i_b(ji,nlay_i) ! is this line still useful ?
            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test
            zdq_i(ji)     =  0._wp
            dh_i_bott(ji) =  0._wp
         ENDIF
      END DO

      DO jk = nlay_i, 1, -1
         DO ji = kideb, kiut
            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN

               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K)

               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting

                  zdeltah   (ji,jk) = - zh_i(ji)  ! internal melting occurs when the internal temperature is above freezing     
                                                  ! this should normally not happen, but sometimes, heat diffusion leads to this

                  dh_i_bott (ji)    = dh_i_bott(ji) + zdeltah(ji,jk)

                  zinnermelt(ji)    = 1._wp

                  zQm               = 0.          ! Not sure which specific enthalpy we should use here (MV HC 2014)
                                                  ! If that happens, heat is probably not well counted
                                                  ! Put zero by default, but more bug analysis should be done to investigate this case

               ELSE                               !!! Basal melting

                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)

                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)

                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0)

                  zfmdt             = - zqfont_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0)

                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change

                  zqfont_bo(ji)     = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * rhoic * ( - zdE ) ! Update heat available

                  zdeltah(ji,jk)    = MAX( zdeltah(ji,jk), - zh_i(ji) ) ! Update thickness change

                  dh_i_bott(ji)     = dh_i_bott(ji) + zdeltah(ji,jk)    ! Update basal melt

                  zfmdt             = - zdeltah(ji,jk) * rhoic          ! Mass flux x time step

                  zQm               = MAX ( zfmdt , 0.0 ) * zEw         ! Heat exchanged with ocean

                  zdq_i(ji)         = zdq_i(ji) + zdeltah(ji,jk) * rhoic * zdE * r1_rdtice ! for heat conservation 

               ENDIF

               ! Contribution to salt flux
               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    &
                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice

               ! Contribution to energy flux to the ocean (J/m2)
               rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)

            ENDIF
         END DO ! ji
      END DO ! jk

      !                                          !-------------------
      IF( con_i .AND. jiindex_1d > 0 ) THEN      ! Conservation test
      !                                          !-------------------
         DO ji = kideb, kiut
            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0  ) THEN
               IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN
                  numce_dh  = numce_dh + 1
                  meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)
               ENDIF
               IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl
                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
                  WRITE(numout,*) ' zfbase    : ', zfbase(ji)
                  WRITE(numout,*) ' zdq_i     : ', zdq_i(ji)
                  WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
                  WRITE(numout,*) ' fc_bo_i   : ', fc_bo_i(ji)
                  WRITE(numout,*) ' fbif_1d   : ', fbif_1d(ji)
                  WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji)
                  WRITE(numout,*) ' s_i_new   : ', s_i_new(ji)
                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij)
                  WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
                  WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) )
               ENDIF
            ENDIF
         END DO
         IF( numce_dh > 0 )   meance_dh = meance_dh / numce_dh
         WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh
         WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh
         WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)
         !
      ENDIF

      !
      !------------------------------------------------------------------------------!
      !  5) Pathological cases                                                       !
      !------------------------------------------------------------------------------!
      !
      !----------------------------------------------
      ! 5.1 Excessive ablation in a 1-category model
      !----------------------------------------------

      DO ji = kideb, kiut
         !                     ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
         IF( jpl == 1 ) THEN   ;   zdhbf = MAX( hmelt , dh_i_bott(ji) )
         ELSE                  ;   zdhbf =              dh_i_bott(ji) 
         ENDIF
         zdvres        = zdhbf - dh_i_bott(ji)
         dh_i_bott(ji) = zdhbf
         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice
         !                     ! excessive energy is sent to lateral ablation
         zinda = MAX( 0._wp, SIGN( 1._wp , 1.0 - at_i_b(ji) - epsi10 ) )
         fsup(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi10 ) * zdvres * r1_rdtice
      END DO

      !-----------------------------------
      ! 5.2 More than available ice melts
      !-----------------------------------
      ! then heat applied minus heat content at previous time step should equal heat remaining 
      !
      DO ji = kideb, kiut
         ! Adapt the remaining energy if too much ice melts
         !--------------------------------------------------
         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) )
         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first
         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom
         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem
         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem

         ! new ice thickness (clem)
         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0
 
         !                     !since ice volume is only used for outputs, we keep it global for all categories
         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji)

         ! remaining heat
         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) ) 

         ! If snow remains, energy is used to melt snow
         zhni =  ht_s_b(ji)      ! snow depth at previous time step
         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow 

         ! energy of melting of remaining snow
         zinda = MAX( 0._wp, SIGN( 1._wp , zhni - epsi10 ) )
         zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi10 ) * zinda
         zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
         zhnfi     =  zhni + zdhnm
         zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
         ht_s_b(ji)     =  MAX( zzero , zhnfi )
         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji)
         ! we recompute dh_s_tot (clem)
         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji)

         ! Mass variations of ice and snow
         !---------------------------------
         !                                              ! mass variation of the jl category
         zzfmass_s = - a_i_b(ji) * ( zhni       - ht_s_b(ji) ) * rhosn   ! snow
         zzfmass_i =   a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic   ! ice  
         !
         zfmass_i(ji) = zzfmass_i                       ! ice variation saved to compute salt flux (see below)
         !
         !                                              ! mass variation cumulated over category
         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow 
         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice 

         ! Remaining heat to the ocean 
         !---------------------------------
         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt

         ! residual salt flux (clem)
         !--------------------------
         ! surface
         sfx_thd_1d(ji)    = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvsur * rhoic * r1_rdtice
         ! bottom
         IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting
            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvbot * rhoic * r1_rdtice
         ELSE                                                          ! growth
            sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice
         ENDIF
         !
         ! diagnostic 
         ii = MOD( npb(ji) - 1, jpi ) + 1
         ij = ( npb(ji) - 1 ) / jpi + 1
         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice
         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
      END DO

      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice

      !---------------------------
      ! heat fluxes                    
      !---------------------------
      DO ji = kideb, kiut
         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
         !
         ! Heat flux
         ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
         ! excessive total  ablation energy (focea) sent to the ocean
         qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice

         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0
         fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   &
            &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice
      END DO  ! ji

      !-------------------------------------------
      ! Correct temperature, energy and thickness
      !-------------------------------------------
      DO ji = kideb, kiut
         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
         t_su_b(ji) =  zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt
      END DO  ! ji

      DO jk = 1, nlay_i
         DO ji = kideb, kiut
            zihgnew      =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 
            t_i_b(ji,jk) =  zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt
            q_i_b(ji,jk) =  zihgnew * q_i_b(ji,jk)
         END DO
      END DO  ! ji

      DO ji = kideb, kiut
         ht_i_b(ji) = zhgnew(ji)
      END DO ! ji 
      
      !
      !------------------------------------------------------------------------------|
      !  6) Snow-Ice formation                                                       |
      !------------------------------------------------------------------------------|
      ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 
      ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
      DO ji = kideb, kiut
         !
         dh_snowice(ji) = MAX(  zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
         zhgnew(ji)     = MAX(  zhgnew(ji) , zhgnew(ji) + dh_snowice(ji)  )
         zhnnew         = MIN(  ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji)  )

         !  Changes in ice volume and ice mass.
         dvnbq_1d  (ji) =                a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) )
         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn

         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 
         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 

         ! Salinity of snow ice
         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1

         zs_snic = zswitch_sal          * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic     &
                 + ( 1. - zswitch_sal ) * sm_i_b(ji)

         ! entrapment during snow ice formation
         ! clem: new salinity difference stored (to be used in limthd_ent.F90)
         IF (  num_sal == 2  ) THEN
            i_ice_switch = MAX( 0._wp , SIGN( 1._wp , zhgnew(ji) - epsi10 ) )
            ! salinity dif due to snow-ice formation
            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch     
            ! salinity dif due to bottom growth 
            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN
               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi10 ) * i_ice_switch
            ENDIF
         ENDIF

         !  Actualize new snow and ice thickness.
         ht_s_b(ji)  = zhnnew
         ht_i_b(ji)  = zhgnew(ji)

         ! Total ablation ! new lines added to debug
         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp

         ! diagnostic ( snow ice growth )
         ii = MOD( npb(ji) - 1, jpi ) + 1 !MV HC 2014 useless
         ij =    ( npb(ji) - 1 ) / jpi + 1 ! MV HC 2014 useless
         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice

         ! Contribution to energy flux to the ocean [J/m2], >0
         zfmdt          = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0.0 )    ! <0
         zsstK          = sst_m(ii,ij) + rt0                                
         zEw            = rcp * ( zsstK - rt0 )
         zQm            = zfmdt * zEw 
         rdq_ice_1d(ji) = rdq_ice_1d(ji) + zQm * a_i_b(ji)

         ! Contribution to salt flux
         sfx_thd_1d(ji) = sfx_thd_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice 
          
         ! Contribution to mass fluxes
         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 
         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 

      END DO !ji
      !
      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i )
      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )
      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i )
      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay )
      !
      CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem
      !
   END SUBROUTINE lim_thd_dh
   
#else
   !!----------------------------------------------------------------------
   !!   Default option                               NO  LIM3 sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_thd_dh          ! Empty routine
   END SUBROUTINE lim_thd_dh
#endif

   !!======================================================================
END MODULE limthd_dh