source: trunk/NEMO/LIM_SRC_3/limthd_dh.F90 @ 869

MODULE limthd_dh
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3'                                      LIM3 sea-ice model
   !!----------------------------------------------------------------------
   !!======================================================================
   !!                       ***  MODULE limthd_dh ***
   !!                thermodynamic growth and decay of the ice 
   !!======================================================================

   !!----------------------------------------------------------------------
   !!   lim_thd_dh  : vertical accr./abl. and lateral ablation of sea ice
   !! * Modules used

   USE par_oce          ! ocean parameters
   USE phycst           ! physical constants (OCE directory) 
   USE ice_oce          ! ice variables
   USE thd_ice
   USE iceini
   USE limistate
   USE in_out_manager
   USE ice
   USE par_ice
   USE lib_mpp
      
   IMPLICIT NONE
   PRIVATE

   !! * Routine accessibility
   PUBLIC lim_thd_dh        ! called by lim_thd

   !! * Module variables
   REAL(wp)  ::           &  ! constant values
      epsi20 = 1e-20   ,  &
      epsi13 = 1e-13   ,  &
      epsi16 = 1e-16   ,  &
      zzero  = 0.e0    ,  &
      zone   = 1.e0

   !!----------------------------------------------------------------------
   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE lim_thd_dh(kideb,kiut,jl)
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dh  ***
       !!------------------------------------------------------------------
       !! ** Purpose :
       !!           This routine 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
       !! ** Steps  
       !!           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
       !!
       !! ** Arguments
       !!
       !! ** Inputs / Outputs
       !!
       !! ** External
       !!
       !! ** References : Bitz and Lipscomb, JGR 99
       !!                 Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
       !!                 Vancoppenolle, Fichefet and Bitz, GRL 2005
       !!                 Vancoppenolle et al., OM08
       !!
       !! ** History  : 
       !!   original code    01-04 (LIM)
       !!   original routine
       !!               (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium
       !!               (06/07-2005) 3D version
       !!               (03-2008)    Clean code
       !!
       !!------------------------------------------------------------------
       !! * Arguments
       INTEGER , INTENT (IN) ::  &
          kideb     ,         &  !: Start point on which the  the computation is applied
          kiut      ,         &  !: End point on which the  the computation is applied
          jl                     !: Thickness cateogry number

       !! * Local variables
       INTEGER ::             & 
          ji        ,         &  !: space index 
          jk        ,         &  !: ice layer index
          isnow     ,         &  !: switch for presence (1) or absence (0) of snow
          zji       ,         &  !: 2D corresponding indices to ji
          zjj       ,         &
          isnowic   ,         &  !: snow ice formation not
          i_ice_switch   ,    &  !: ice thickness above a certain treshold or not
          iter

       REAL(wp) ::            &
          zhsnew    ,         &  !: new snow thickness
          zihgnew   ,         &  !: switch for total ablation
          ztmelts   ,         &  !: melting point
          zhn       ,         &
          zdhcf     ,         &
          zdhbf     ,         &
          zhni      ,         &
          zhnfi     ,         &
          zihg      ,         &
          zdhnm     ,         &
          zhnnew    ,         &
          zeps = 1.0e-13,     &
          zhisn     ,         &
          zfracs    ,         &  !: fractionation coefficient for bottom salt
                                 !: entrapment
          zds       ,         &  !: increment of bottom ice salinity
          zcoeff    ,         &  !: dummy argument for snowfall partitioning
                                 !: over ice and leads
          zsm_snowice,        &  !: snow-ice salinity
          zswi1     ,         &  !: switch for computation of bottom salinity
          zswi12    ,         &  !: switch for computation of bottom salinity
          zswi2     ,         &  !: switch for computation of bottom salinity
          zgrr      ,         &  !: bottom growth rate
          zihic     ,         &  !: 
          ztform                 !: bottom formation temperature

       REAL(wp) , DIMENSION(jpij) ::  &
          zh_i      ,         &  ! ice layer thickness
          zh_s      ,         &  ! snow layer thickness
          ztfs      ,         &  ! melting point
          zhsold    ,         &  ! old snow thickness
          zqprec    ,         &  !: energy of fallen snow
          zqfont_su ,         &  ! incoming, remaining surface energy
          zqfont_bo              ! incoming, bottom energy

       REAL(wp) , DIMENSION(jpij) :: &
          z_f_surf,           &  ! surface heat for ablation
          zhgnew                 ! new ice thickness

       REAL(wp), DIMENSION(jpij) :: &
          zdh_s_mel  ,        &  ! snow melt 
          zdh_s_pre  ,        &  ! snow precipitation 
          zdh_s_sub  ,        &  ! snow sublimation
          zfsalt_melt            ! salt flux due to ice melt

       REAL(wp) , DIMENSION(jpij,jkmax) :: &
          zdeltah

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

       REAL(wp), DIMENSION(jpij,jkmax) :: &
          zqt_i_lay              !: total ice heat content

       ! Heat conservation 
       REAL(wp), DIMENSION(jpij) :: &
          zfbase,             &
          zdq_i

       INTEGER, DIMENSION(jpij) ::  &
          innermelt

       REAL(wp) :: &
          meance_dh

       INTEGER ::                   &
          num_iter_max,       &
          numce_dh

!!-----------------------------------------------------------------------------

      WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation'
      WRITE(numout,*) '~~~~~~~~~'

      zfsalt_melt(:)  = 0.0
      ftotal_fin(:)   = 0.0
      zfdt_init(:)    = 0.0
      zfdt_final(:)   = 0.0

      DO ji = kideb, kiut
         old_ht_i_b(ji) = ht_i_b(ji)
         old_ht_s_b(ji) = ht_s_b(ji)
      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.0
      zqfont_bo(:) = 0.0
      dsm_i_se_1d(:) = 0.0     
      dsm_i_si_1d(:) = 0.0     
!
!------------------------------------------------------------------------------!
!  2) Computing layer thicknesses and  snow and sea-ice enthalpies.            !
!------------------------------------------------------------------------------!
!
      ! Layer thickness
      DO ji = kideb,kiut
         zh_i(ji) = ht_i_b(ji) / nlay_i
         zh_s(ji) = ht_s_b(ji) / nlay_s
      END DO

      ! Total enthalpy of the snow
      zqt_s(:) = 0.0
      DO jk = 1, nlay_s
         DO ji = kideb,kiut
            zqt_s(ji)    =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s
         END DO
      END DO

      ! Total enthalpy of the ice
      zqt_i(:) = 0.0
      DO jk = 1, nlay_i
         DO ji = kideb,kiut
            zqt_i(ji)        =  zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
            zqt_i_lay(ji,jk) =  q_i_b(ji,jk) * ht_i_b(ji) / nlay_i
         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.0

      ! Melt of fallen snow
      DO ji = kideb, kiut
         ! tatm_ice is now in K
         zqprec(ji)     =  rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus )  
         zqfont_su(ji)  =  z_f_surf(ji) * rdt_ice
         zdeltah(ji,1)  =  MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) )
         zqfont_su(ji)  =  MAX( 0.0 , - 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)
      END DO

      ! Update total snow heat content
      zqt_s(ji)         =  MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 
      IF( lk_mpp ) CALL mpp_max(zqt_s(ji), kcom = ncomm_ice ) 

      ! 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 )
         ! Volume and mass variations of snow
         dvsbq_1d(ji)   =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji)              &
                        - zdh_s_mel(ji) )
         dvsbq_1d(ji)   =  MIN( zzero, dvsbq_1d(ji) )
         rdmsnif_1d(ji) =  rhosn*dvsbq_1d(ji)
      END DO ! ji

      !--------------------------
      ! 3.2 Surface ice ablation 
      !--------------------------
      DO ji = kideb, kiut 
         dh_i_surf(ji) =  0.0
         ! For heat conservation test
         z_f_surf(ji)  =  zqfont_su(ji) / rdt_ice ! heat conservation test
         zdq_i(ji)     =  0.0
      END DO ! ji

      DO jk = 1, nlay_i
         DO ji = kideb, kiut 
            ! melt of layer jk
            zdeltah(ji,jk)      = - zqfont_su(ji) / q_i_b(ji,jk)
            ! recompute heat available
            zqfont_su(ji)       = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) *  &
                                  q_i_b(ji,jk) 
            ! melt of layer jk cannot be higher than its thickness
            zdeltah(ji,jk)      = MAX( zdeltah(ji,jk) , - zh_i(ji) )
            ! update surface melt
            dh_i_surf(ji)       = dh_i_surf(ji) + zdeltah(ji,jk) 
            ! for energy conservation
            zdq_i(ji)           = zdq_i(ji) + zdeltah(ji,jk) *                &
                                        q_i_b(ji,jk) / rdt_ice
            ! contribution to ice-ocean salt flux 
            zji                 = MOD( npb(ji) - 1, jpi ) + 1
            zjj                 = ( npb(ji) - 1 ) / jpi + 1
            zfsalt_melt(ji)     = zfsalt_melt(ji) +                           &
                                  ( sss_io(zji,zjj) - sm_i_b(ji)   ) *        &
                                  a_i_b(ji) *                                 &
                                  MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
         END DO ! ji
      END DO ! jk

      !-------------------
      ! Conservation test
      !-------------------
      IF ( con_i ) THEN
      numce_dh = 0
      meance_dh = 0.0
      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,*) ' zji, zjj, jl :', zji, zjj, 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_io  : ', sss_io(zji,zjj)
         ENDIF

      END DO ! ji

      IF ( numce_dh .GT. 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 ! con_i

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

      DO ji = kideb, kiut
         ! if qla is positive (upwards), heat goes to the atmosphere, therefore
         ! snow sublimates, if qla is negative (downwards), snow condensates
         zdh_s_sub(ji)   =  - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice
         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  !ji

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

      DO jk = 1, nlay_s !n 
         DO ji = kideb, kiut !n
         ! 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                                                       !
!------------------------------------------------------------------------------!
!
      ! Ice basal growth / melt is given by the ratio of heat budget over basal
      ! ice heat content.  Basal heat budget is given by the difference 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 

      !-----------------------------------------------------
      ! 4.1 Basal growth - (a) salinity not varying in time 
      !-----------------------------------------------------
      IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN
         DO ji = kideb, kiut
            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
               s_i_new(ji)         =  sm_i_b(ji)
               ! Melting point in K
               ztmelts             =  - tmut * s_i_new(ji) + rtt 
               ! New ice heat content (Bitz and Lipscomb, 1999)
               ztform              =  t_i_b(ji,nlay_i)  ! t_bo_b crashes in the
                                                        ! Baltic
               q_i_b(ji,nlay_i+1)  =  rhoic * &
                                      ( cpic * ( ztmelts - ztform     )       &
                                      + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
                                                 ( ztform     - rtt ) )       &
                                      - rcp * ( ztmelts-rtt ) )
               ! Basal growth rate = - F*dt / q
               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
                                      qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
            ENDIF ! heat budget
         END DO ! ji
      ENDIF ! num_sal

      !-------------------------------------------------
      ! 4.1 Basal growth - (b) salinity varying in time 
      !-------------------------------------------------
      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
         ! the growth rate (dh_i_bott) is function of the new ice
         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice 
         ! salinity (snewice). snewice depends on dh_i_bott
         ! it converges quickly, so, no problem
         ! See Vancoppenolle et al., OM08 for more info on this

         ! Initial value (tested 1D, can be anything between 1 and 20)
         num_iter_max = 4
         s_i_new(:) = 4.0

         ! Iterative procedure
         DO iter = 1, num_iter_max
            DO ji = kideb, kiut
               IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
                  zji = MOD( npb(ji) - 1, jpi ) + 1
                  zjj = ( npb(ji) - 1 ) / jpi + 1
                  ! Melting point in K
                  ztmelts             =   - tmut * s_i_new(ji) + rtt 
                  ! New ice heat content (Bitz and Lipscomb, 1999)
                  q_i_b(ji,nlay_i+1)  =  rhoic * &
                                      ( cpic * ( ztmelts - t_bo_b(ji) )       &
                                      + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
                                                 ( t_bo_b(ji) - rtt ) )       &
                                      - rcp * ( ztmelts-rtt ) )
                  ! Bottom growth rate = - F*dt / q
                  dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) &
                                      + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
                  ! New ice salinity ( Cox and Weeks, JGR, 1988 )
                  ! zswi2  (1) if dh_i_bott/rdt .GT. 3.6e-7
                  ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
                  ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) )
                  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 = zswi1  * 0.12 +  &
                           zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + &
                           zswi2  * 0.26 /  &
                           ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
                  zds         = zfracs*sss_io(zji,zjj) - s_i_new(ji)
                  s_i_new(ji) = zfracs * sss_io(zji,zjj)
               ENDIF ! fc_bo_i
            END DO ! ji
         END DO ! iter

         ! Final values
         DO ji = kideb, kiut
            IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
               ! New ice salinity must not exceed 15 psu
               s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
               ! Metling point in K
               ztmelts     =   - tmut * s_i_new(ji) + rtt 
               ! New ice heat content (Bitz and Lipscomb, 1999)
               q_i_b(ji,nlay_i+1)  =  rhoic *                              &
                                   ( cpic * ( ztmelts - t_bo_b(ji) )       &
                                   + lfus * ( 1.0 - ( ztmelts - rtt ) /    &
                                              ( t_bo_b(ji) - rtt ) )       &
                                   - rcp * ( ztmelts-rtt ) )
               ! Basal growth rate = - F*dt / q
               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + &
                                      qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
               ! Salinity update
               ! entrapment during bottom growth
               dsm_i_se_1d(ji) =  ( s_i_new(ji)*dh_i_bott(ji) +              & 
                                   sm_i_b(ji)*ht_i_b(ji) ) /                 & 
                                   MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps )   &
                                   - sm_i_b(ji)
            ENDIF ! heat budget
         END DO ! ji
      ENDIF ! num_sal

      !----------------
      ! 4.2 Basal melt
      !----------------
      meance_dh = 0.0
      numce_dh = 0
      innermelt(:) = 0

      DO ji = kideb, kiut
         ! heat convergence at the surface > 0
         IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
                 
            s_i_new(ji)   =  s_i_b(ji,nlay_i)
            zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )

            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice ! heat conservation test
            zdq_i(ji)     =  0.0

            dh_i_bott(ji) =  0.0
         ENDIF
      END DO

      DO jk = nlay_i, 1, -1
         DO ji = kideb, kiut
            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
               ztmelts             =   - tmut * s_i_b(ji,jk) + rtt 
               IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN
                  zdeltah(ji,jk)  = - zh_i(ji)
                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
                  innermelt(ji)   = 1
               ELSE  ! normal ablation
                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * &
                                    q_i_b(ji,jk)
                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * &
                                    q_i_b(ji,jk) / rdt_ice
               ! contribution to salt flux
                  zji             = MOD( npb(ji) - 1, jpi ) + 1
                  zjj             = ( npb(ji) - 1 ) / jpi + 1
                  zfsalt_melt(ji) = zfsalt_melt(ji) +                         &
                                   ( sss_io(zji,zjj) - sm_i_b(ji)   ) *       &
                                   a_i_b(ji) * &
                                   MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
               ENDIF
            ENDIF
         END DO ! ji
      END DO ! jk

      !-------------------
      ! Conservation test
      !-------------------
      IF ( con_i ) THEN
      DO ji = kideb, kiut
         IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
            IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-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 '
                WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, 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_io  : ', sss_io(zji,zjj)
                WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
                WRITE(numout,*) ' innermelt : ', innermelt(ji)
            ENDIF
         ENDIF ! heat convergence at the surface
      END DO ! ji

      IF ( numce_dh .GT. 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 ! con_i

!
!------------------------------------------------------------------------------!
!  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)
         zdhbf = dh_i_bott(ji) 
         IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )
         ! excessive energy is sent to lateral ablation
         fsup(ji)            =  rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &
                             *  ( zdhbf - dh_i_bott(ji) ) / rdt_ice

         dh_i_bott(ji)  = zdhbf
         !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)
         !new ice thickness
         zhgnew(ji)     = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)

         ! diagnostic ( bottom ice growth )
         zji                 = MOD( npb(ji) - 1, jpi ) + 1
         zjj                 = ( npb(ji) - 1 ) / jpi + 1
         diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) &
                              / rdt_ice
         diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) &
                              / rdt_ice
         diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) &
                              / rdt_ice
      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
         !--------------------------------------------------
         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
         ! 0 if no more ice
         zhgnew(ji) =  zihgnew * zhgnew(ji) ! ice thickness is put to 0
         ! 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
         zqt_s(ji)  =  ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps )
         zdhnm      =  - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) /  &
                           MAX( zqt_s(ji) , zeps ) )
         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)

         ! Mass variations of ice and snow
         !---------------------------------
         rdmicif_1d(ji) =  rdmicif_1d(ji) + a_i_b(ji) *                        &
                           (zhgnew(ji)-ht_i_b(ji))*rhoic ! good

         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + a_i_b(ji) * &
                           (ht_s_b(ji)-zhni)*rhosn ! good too

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

      END DO

      ftotal_fin (:) = zfdt_final(:)  / rdt_ice

      !---------------------------
      ! Salt flux and heat fluxes                    
      !---------------------------
      DO ji = kideb, kiut
         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice

         ! Salt flux
         zji           = MOD( npb(ji) - 1, jpi ) + 1
         zjj           = ( npb(ji) - 1 ) / jpi + 1
         IF ( num_sal .NE. 4 ) &
         fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
                          (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
                          ( sss_io(zji,zjj) - sm_i_b(ji) ) / rdt_ice
         ! new lines
         IF ( num_sal .EQ. 4 ) &
         fseqv_1d(ji)  = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) +           &
                          (1.0 - zihgnew) * rdmicif_1d(ji) *                  &
                          ( sss_io(zji,zjj) - bulk_sal ) / rdt_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

#if defined key_lim_fdd 
!(presently Activated) 
         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 
                                         * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 
         rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 
                                         * ( zhnnew       - ht_s_b(ji) )*rhosn

!        Equivalent salt flux (1) Snow-ice formation component
!        -----------------------------------------------------
         zji                 = MOD( npb(ji) - 1, jpi ) + 1
         zjj                 = ( npb(ji) - 1 ) / jpi + 1

         zsm_snowice  = ( rhoic - rhosn ) / rhoic *            &
                        sss_io(zji,zjj) 

         IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji)

         IF ( num_sal .NE. 4 ) &
         fseqv_1d(ji)   = fseqv_1d(ji)   + &
                          ( sss_io(zji,zjj) - zsm_snowice ) * &
                            a_i_b(ji)   * &
                          ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
         ! new lines
         IF ( num_sal .EQ. 4 ) &
         fseqv_1d(ji)   = fseqv_1d(ji)   + &
                          ( sss_io(zji,zjj) - bulk_sal    ) * &
                            a_i_b(ji)   * &
                          ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice

         ! entrapment during snow ice formation
         i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
         isnowic      = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * &
                        i_ice_switch
         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
             dsm_i_si_1d(ji)  = ( zsm_snowice*dh_snowice(ji) &
                            + sm_i_b(ji) * ht_i_b(ji)                          & 
                            / MAX( ht_i_b(ji) + dh_snowice(ji), zeps)          &
                            - sm_i_b(ji) ) * isnowic     

#else
         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) &
                                         * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic &
                                         + ( zhnnew - ht_s_b(ji) ) * rhosn )
#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).LE.0.0 ) a_i_b(ji) = 0.0

         ! diagnostic ( snow ice growth )
         zji                 = MOD( npb(ji) - 1, jpi ) + 1
         zjj                 = ( npb(ji) - 1 ) / jpi + 1
         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / &
                                 rdt_ice

      END DO !ji

    END SUBROUTINE lim_thd_dh
#else
   !!======================================================================
   !!                       ***  MODULE limthd_dh   ***
   !!                              no sea ice model  
   !!======================================================================
CONTAINS
   SUBROUTINE lim_thd_dh          ! Empty routine
   END SUBROUTINE lim_thd_dh
#endif
 END MODULE limthd_dh