source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90 @ 8626

MODULE icethd_dh
   !!======================================================================
   !!                       ***  MODULE icethd_dh ***
   !!   seaice : thermodynamic growth and melt 
   !!======================================================================
   !! 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 wfx_snw & wfx_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'                                       ESIM sea-ice model
   !!----------------------------------------------------------------------
   !!   ice_thd_dh        : vertical sea-ice growth and melt
   !!   ice_thd_snwblow   : distribute snow fall between ice and ocean
  !!----------------------------------------------------------------------
   USE dom_oce        ! ocean space and time domain
   USE phycst         ! physical constants
   USE ice            ! sea-ice: variables
   USE ice1D          ! sea-ice: thermodynamics variables
   USE icethd_sal     ! sea-ice: salinity profiles
   !
   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_dh        ! called by ice_thd
   PUBLIC   ice_thd_snwblow   ! called in sbcblk/sbcclio/sbccpl and here

   INTERFACE ice_thd_snwblow
      MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
   END INTERFACE

   !!----------------------------------------------------------------------
   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
   !! $Id: icethd_dh.F90 8420 2017-08-08 12:18:46Z clem $
   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE ice_thd_dh
      !!------------------------------------------------------------------
      !!                ***  ROUTINE ice_thd_dh  ***
      !!
      !! ** Purpose :   compute ice and snow thickness changes due to growing/melting
      !!
      !! ** 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  ::   ji, jk       ! dummy loop indices
      INTEGER  ::   iter         ! local integer

      REAL(wp) ::   ztmelts      ! local scalar
      REAL(wp) ::   zdum       
      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
      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) ::   z1_rho       ! 1/(rhosn+rau0-rhoic)

      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), DIMENSION(jpij) ::   zqprec      ! energy of fallen snow                       (J.m-3)
      REAL(wp), DIMENSION(jpij) ::   zq_su       ! heat for surface ablation                   (J.m-2)
      REAL(wp), DIMENSION(jpij) ::   zq_bo       ! heat for bottom ablation                    (J.m-2)
      REAL(wp), DIMENSION(jpij) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
      REAL(wp), DIMENSION(jpij) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
      REAL(wp), DIMENSION(jpij) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)

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

      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
      REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
      INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting 

      REAL(wp), DIMENSION(jpij) ::   zeh_i       ! total ice heat content  (J.m-2)
      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing

      REAL(wp) ::   zswitch_sal

      INTEGER  ::   num_iter_max      ! Heat conservation 
      !!------------------------------------------------------------------

      ! Discriminate between time varying salinity and constant
      SELECT CASE( nn_icesal )                  ! varying salinity or not
         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
      END SELECT

      ! initialize layer thicknesses and enthalpies
      h_i_old (1:npti,0:nlay_i+1) = 0._wp
      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
      DO jk = 1, nlay_i
         DO ji = 1, npti
            h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
            eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
         END DO
      END DO
      !
      !------------------------------------------------------------------------------!
      !  1) Calculate available heat for surface and bottom ablation                 !
      !------------------------------------------------------------------------------!
      !
      DO ji = 1, npti
         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 

         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
      END DO

      !------------------------------------------------------------------------------!
      ! If snow temperature is above freezing point, then snow melts 
      ! (should not happen but sometimes it does)
      !------------------------------------------------------------------------------!
      DO ji = 1, npti
         IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
            ! Contribution to heat flux to the ocean [W.m-2], < 0  
            hfx_res_1d(ji) = hfx_res_1d(ji) + e_s_1d(ji,1) * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
            ! Contribution to mass flux
            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji)
            ! updates
            h_s_1d(ji)    = 0._wp
            e_s_1d (ji,1) = 0._wp
            t_s_1d (ji,1) = rt0
         END IF
      END DO

      !------------------------------------------------------------!
      !  2) Computing layer thicknesses and enthalpies.            !
      !------------------------------------------------------------!
      DO jk = 1, nlay_i
         DO ji = 1, npti
            zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
            zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
         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

      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing

      zdeltah(1:npti,:) = 0._wp
      DO ji = 1, npti
         !-----------
         ! Snow fall
         !-----------
         ! thickness change
         zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji)
         ! enthalpy of the precip (>0, J.m-3)
         zqprec   (ji) = - qprec_ice_1d(ji)   
         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
         ! heat flux from snow precip (>0, W.m-2)
         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
         ! mass flux, <0
         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice

         !---------------------
         ! Melt of falling snow
         !---------------------
         ! thickness change
         rswitch        = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
         zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
         zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 
         ! heat used to melt snow (W.m-2, >0)
         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
         ! snow melting only = water into the ocean (then without snow precip), >0
         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
         ! updates available heat + precipitations after melting
         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
         zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1)

         ! update thickness
         h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
      END DO

      ! If heat still available (zq_su > 0), then melt more snow
      zdeltah(1:npti,:) = 0._wp
      zdh_s_mel(1:npti) = 0._wp
      DO jk = 1, nlay_s
         DO ji = 1, npti
            ! thickness change
            rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) 
            rswitch          = rswitch * ( MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) ) 
            zdeltah  (ji,jk) = - rswitch * zq_su(ji) / MAX( e_s_1d(ji,jk), epsi20 )
            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - h_s_1d(ji) ) ! bound melting
            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
            ! heat used to melt snow(W.m-2, >0)
            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
            ! snow melting only = water into the ocean (then without snow precip)
            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk)
            ! updates available heat + thickness
            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
         END DO
      END DO

      !------------------------------
      ! 3.2 Sublimation (part1: snow) 
      !------------------------------
      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
      ! clem comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
      zdeltah(1:npti,:) = 0._wp
      DO ji = 1, npti
         zdh_s_sub(ji)  = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
         ! remaining evap in kg.m-2 (used for ice melting later on)
         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
         zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
         hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1)  &
            &                              ) * a_i_1d(ji) * r1_rdtice
         ! Mass flux by sublimation
         wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice

         ! new snow thickness
         h_s_1d(ji)    =  MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
         ! update precipitations after sublimation and correct sublimation
         zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
         zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
      END DO
      
      ! --- Update snow diags --- !
      DO ji = 1, npti
         dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
      END DO

      !-------------------------------------------
      ! 3.3 Update temperature, energy
      !-------------------------------------------
      ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
      DO jk = 1, nlay_s
         DO ji = 1,npti
            rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
            e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *           &
              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
              &              ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
         END DO
      END DO

      !--------------------------
      ! 3.4 Surface ice ablation 
      !--------------------------
      zdeltah(1:npti,:) = 0._wp ! important
      DO jk = 1, nlay_i
         DO ji = 1, npti
            ztmelts           = - tmut * sz_i_1d(ji,jk)          ! Melting point of layer k [C]
            
            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting

               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
               zdE            = 0._wp                                 ! Specific enthalpy difference   (J/kg, <0)
                                                                      ! set up at 0 since no energy is needed to melt water...(it is already melted)
               zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )          ! internal melting occurs when the internal temperature is above freezing     
                                                                      ! this should normally not happen, but sometimes, heat diffusion leads to this
               zfmdt          = - zdeltah(ji,jk) * rhoic              ! Mass flux x time step > 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]

               ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 
               hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
               
               ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
               sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
               
               ! Contribution to mass flux
               wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice

            ELSE                                !!! Surface melting
               
               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
               
               zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
               
               zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
               
               zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
               
               zq_su(ji)      = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
               
               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            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
               
               ! Contribution to salt flux >0 (clem: using s_i_1d and not sz_i_1d(jk) is ok)
               sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
               
               ! Contribution to heat flux [W.m-2], < 0
               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
               
               ! Total heat flux used in this process [W.m-2], > 0  
               hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
               
               ! Contribution to mass flux
               wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
               
            END IF
            ! ----------------------
            ! Sublimation part2: ice
            ! ----------------------
            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted.
            !                          It must be corrected at some point)
            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_rdtice
            ! Heat flux [W.m-2], < 0
            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
            ! Mass flux > 0
            wfx_ice_sub_1d(ji) =  wfx_ice_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice

            ! update remaining mass flux
            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
            
            ! record which layers have disappeared (for bottom melting) 
            !    => icount=0 : no layer has vanished
            !    => icount=5 : 5 layers have vanished
            rswitch       = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 
            icount(ji,jk) = NINT( rswitch )
            zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
                        
            ! update heat content (J.m-2) and layer thickness
            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
            h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
         END DO
      END DO
      ! update ice thickness
      DO ji = 1, npti
         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
      END DO

      ! remaining "potential" evap is sent to ocean
      DO ji = 1, npti
         wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
      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 
      ! (fhld) and the turbulent ocean flux (fhtur). 
      ! fc_bo_i is positive downwards. fhtur and fhld 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

      num_iter_max = 1
      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
      
      ! Iterative procedure
      DO ji = 1, npti
         IF(  zf_tt(ji) < 0._wp  ) THEN
            DO iter = 1, num_iter_max

               ! 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 , epsi10 ) )
               zswi2              = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
               zswi12             = MAX( 0._wp , SIGN( 1._wp , 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 )

               s_i_new(ji)        = zswitch_sal * zfracs * sss_1d(ji)  &  ! New ice salinity
                                  + ( 1. - zswitch_sal ) * s_i_1d(ji) 
               ! New ice growth
               ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C)

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

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

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

               dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
               
            END DO
            ! Contribution to Energy and Salt Fluxes                                    
            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)

            ztmelts        = - tmut * s_i_new(ji)                ! New ice melting point (C)
            
            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
            
            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
               &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
               &               + rcp  * ztmelts          
            
            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
            
            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
            
            ! Contribution to heat flux to the ocean [W.m-2], >0  
            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice

            ! Total heat flux used in this process [W.m-2], <0  
            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
            
            ! Contribution to salt flux, <0
            sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice

            ! Contribution to mass flux, <0
            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice

            ! update heat content (J.m-2) and layer thickness
            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * (-zEi * rhoic)
            h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)

         ENDIF

      END DO

      !----------------
      ! 4.2 Basal melt
      !----------------
      zdeltah(1:npti,:) = 0._wp ! important
      DO jk = nlay_i, 1, -1
         DO ji = 1, npti
            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting 

               ztmelts = - tmut * sz_i_1d(ji,jk)  ! Melting point of layer jk (C)

               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting

                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
                  zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
                                                                    ! set up at 0 since no energy is needed to melt water...(it is already melted)
                  zdeltah   (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) )  ! 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)

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

                  ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 
                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice

                  ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
                  sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
                                    
                  ! Contribution to mass flux
                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice

                  ! update heat content (J.m-2) and layer thickness
                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)

               ELSE                               !!! Basal melting

                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
                  zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0)
                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)

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

                  zdeltah(ji,jk)  = - zfmdt * r1_rhoic         ! Gross thickness change

                  zdeltah(ji,jk)  = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
                  
                  zq_bo(ji)       = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors

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

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

                  zQm             = zfmdt * zEw         ! Heat exchanged with ocean

                  ! Contribution to heat flux to the ocean [W.m-2], <0  
                  hfx_thd_1d(ji)  = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice

                  ! Contribution to salt flux (clem: using s_i_1d and not sz_i_1d(jk) is ok)
                  sfx_bom_1d(ji)  = sfx_bom_1d(ji) - rhoic *  a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_rdtice
                  
                  ! Total heat flux used in this process [W.m-2], >0  
                  hfx_bom_1d(ji)  = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
                  
                  ! Contribution to mass flux
                  wfx_bom_1d(ji)  =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice

                  ! update heat content (J.m-2) and layer thickness
                  eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
               ENDIF
           
            ENDIF
         END DO
      END DO

      !-------------------------------------------
      ! Update temperature, energy
      !-------------------------------------------
      DO ji = 1, npti
         h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) )
      END DO  

      !-------------------------------------------
      ! 5. What to do with remaining energy
      !-------------------------------------------
      ! If heat still available for melting and snow remains, then melt more snow
      !-------------------------------------------
      zdeltah(1:npti,:) = 0._wp ! important
      DO ji = 1, npti
         zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
         rswitch         = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
         zdeltah  (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1)
        
         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
         ! heat used to melt snow
         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
         ! Contribution to mass flux
         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
         !    
         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
         hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice

         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
      END DO

      !
      !------------------------------------------------------------------------------|
      !  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
      z1_rho = 1._wp / ( rhosn+rau0-rhoic )
      DO ji = 1, npti
         !
         dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho )

         h_i_1d(ji)    = h_i_1d(ji) + dh_snowice(ji)
         h_s_1d(ji)    = h_s_1d(ji) - dh_snowice(ji)

         ! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
         zfmdt          = ( rhosn - rhoic ) * dh_snowice(ji)    ! <0
         zEw            = rcp * sst_1d(ji)
         zQm            = zfmdt * zEw 
         
         ! Contribution to heat flux
         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 

         ! Contribution to salt flux
         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 

         ! Case constant salinity in time: virtual salt flux to keep salinity constant
         IF( nn_icesal /= 2 )  THEN
            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean 
         ENDIF

         ! Contribution to mass flux
         ! All snow is thrown in the ocean, and seawater is taken to replace the volume
         wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
         wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice

         ! update heat content (J.m-2) and layer thickness
         eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
         h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
         
      END DO

      !
      !-------------------------------------------
      ! Update temperature, energy
      !-------------------------------------------
      DO ji = 1, npti
         rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
         t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
      END DO

      DO jk = 1, nlay_s
         DO ji = 1,npti
            ! mask enthalpy
            rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - h_s_1d(ji) )  )
            e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
            ! recalculate t_s_1d from e_s_1d
            t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhosn * r1_cpic + lfus * r1_cpic )
         END DO
      END DO

      ! --- ensure that a_i = 0 where h_i = 0 ---
      WHERE( h_i_1d(1:npti) == 0._wp )   a_i_1d(1:npti) = 0._wp
      !
   END SUBROUTINE ice_thd_dh


   !!--------------------------------------------------------------------------
   !! INTERFACE ice_thd_snwblow
   !! ** Purpose :   Compute distribution of precip over the ice
   !!--------------------------------------------------------------------------
!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
   SUBROUTINE ice_thd_snwblow_2d( pin, pout )
      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
      REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
      pout = ( 1._wp - ( pin )**rn_blow_s )
   END SUBROUTINE ice_thd_snwblow_2d

   SUBROUTINE ice_thd_snwblow_1d( pin, pout )
      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
      pout = ( 1._wp - ( pin )**rn_blow_s )
   END SUBROUTINE ice_thd_snwblow_1d

#else
   !!----------------------------------------------------------------------
   !!   Default option                                NO ESIM sea-ice model
   !!----------------------------------------------------------------------
#endif

   !!======================================================================
END MODULE icethd_dh