!@sum contains code related to 3 layers snow model
!@auth I.Aleinov
!@ver  1.0




      MODULE SNOW_MODEL 4,2
!@sum SNOW_MODEL does column physics for 3 layers snow model
!@auth I.Aleinov
!@ver  1.0
      USE FILEMANAGER, only: openunit
      USE CONSTANT, only :
     $     rho_water => rhow,   ! (kg m-3)
     $     rho_ice => rhoi,     ! (kg m-3)
     $     lhm_kg => lhm,       ! (J kg-1)
     $     lhe_kg => lhe,       ! (J kg-1)
     $     shi_kg => shi,       ! (J kg-1 C-1)
     $     tfrz => tf,          ! (K)
     $     sigma => stbo        ! (W m-2 K-4)

      IMPLICIT NONE
      SAVE
      PRIVATE

      PUBLIC snow_adv, snow_redistr, snow_fraction, i_earth, j_earth
      public MIN_SNOW_THICKNESS, MIN_FRACT_COVER, rho_water,
     &     rho_fresh_snow

ccc physical parameters
!@var rho_fresh_snow density of fresh snow (kg m-3)
      real*8, parameter :: rho_fresh_snow =  150.d0 ! (kg m-3)
!@var lat_fusion volumetric latent heat of fusion for water  (J m-3)
      real*8, parameter :: lat_fusion = lhm_kg * rho_water ! (J m-3)
!@var lat_evap volumetric latent heat of evaporation for water  (J m-3)
      real*8, parameter :: lat_evap = lhe_kg * rho_water ! (J m-3)
!@var max_fract_water max fraction of free water in a wet snow (1)
      real*8, parameter :: max_fract_water = .055d0

ccc model parameters
!@var EPS small number (used as cutoff parameter) (1)
      real*8, parameter :: EPS = 1.d-8  ! was 1.d-12
!@var MAX_NL maximal number of snow layers (only 3 are used, really) (1)
      integer, parameter :: MAX_NL = 16
!@var TOTAL_NL maximal number of snow layers (actual number used) (1)
ccc actually MAX_NL can be set = TOTAL_NL but should test first ...
      integer, parameter :: TOTAL_NL=3
!@var MIN_SNOW_THICKNESS minimal thickness of snow (m)
ccc trying to increase MIN_SNOW_THICKNESS for stability
ccc      real*8, parameter :: MIN_SNOW_THICKNESS =  0.01d0  ! was 0.09d0
      real*8, parameter :: MIN_SNOW_THICKNESS =  0.1d0
!@var MIN_FRACT_COVER minimal snow cover (cutoff param) (1)
      real*8, parameter :: MIN_FRACT_COVER = 0.0001d0

!@var DEB_CH channel for debug output
      integer :: DEB_CH = 0
!@var i_earth, j_earth coordinate of current point (for debugging)
      integer i_earth, j_earth

      CONTAINS


      subroutine pass_water( wsn, hsn, dz, nl, 2
     &                       water_down, heat_down,
     &                       lat_fusion, max_fract_water,
     &                       rho_water, rho_fresh_snow)
!@sum  pass extra water to the lower layers
!@auth I.Aleinov
!@ver  1.0
      implicit none
      integer nl
      real*8  wsn(nl), hsn(nl), dz(nl)
      real*8 water_down, heat_down, lat_fusion, max_fract_water
      real*8 rho_water, rho_fresh_snow
      integer n
      real*8 ice, free_water, ice_old
      integer pass_water_flag

      pass_water_flag = 0
      do n=1,nl
        ice_old = min(wsn(n),-hsn(n)/lat_fusion)
        ice = ice_old
        wsn(n) = wsn(n) + water_down
        hsn(n) = hsn(n) + heat_down
        water_down = 0.d0
        heat_down = 0.d0
        if( hsn(n) .ge. 0.d0 .or. wsn(n) .le. 0.d0 ) then ! all melted
          water_down = wsn(n)
          heat_down = hsn(n)
          wsn(n) = 0.d0
          hsn(n) = 0.d0
          dz(n) = 0.d0
          pass_water_flag = 1
        else if (hsn(n).gt.-wsn(n)*lat_fusion) then !
          ice = -hsn(n)/lat_fusion
          free_water = wsn(n) - ice
c!!!                may be ice*max_fract_water ??
          water_down = max(0d0, free_water - ice*max_fract_water)
          wsn(n) = wsn(n) - water_down
          !
          !
          dz(n) = min( dz(n), ice*rho_water/rho_fresh_snow )
          pass_water_flag = 2
        else !
          if (wsn(n)+EPS .lt. ice_old) dz(n) = dz(n)*wsn(n)/ice_old
          dz(n) = min( dz(n), wsn(n)*rho_water/rho_fresh_snow )
        endif
        dz(n) = max( dz(n), wsn(n)*rho_water/rho_ice )
      enddo
      return
      end subroutine pass_water



      subroutine snow_fraction( 2
     &     dz, nl, prsnow, dt, fract_cover, fract_cover_new)
!@sum computes new snow fraction and returns it in fract_cover_new
!@auth I.Aleinov
!@ver  1.0
      implicit none
      integer nl
      real*8 dz(TOTAL_NL+1), prsnow, dt
      real*8 fract_cover, fract_cover_new
      real*8 dz_aver, fresh_snow

      fresh_snow = rho_water/rho_fresh_snow * prsnow*dt

      dz_aver = sum(dz(1:nl))*fract_cover + fresh_snow

      fract_cover_new = min( 1.d0, dz_aver/MIN_SNOW_THICKNESS )
      if ( fract_cover_new < MIN_FRACT_COVER ) fract_cover_new = 0.d0

      return
      end subroutine snow_fraction



      subroutine snow_redistr( 4,3
     &     dz, wsn, hsn, nl, fract_cover_ratio, tr_flux, dt )
!@sum  redistributes snow between the layers
!@auth I.Aleinov
!@ver  1.0
      implicit none
      integer nl
      real*8 dz(TOTAL_NL+1), wsn(TOTAL_NL), hsn(TOTAL_NL)
      real*8 fract_cover_ratio
      integer nlo
      real*8 dzo(TOTAL_NL+1), wsno(TOTAL_NL), hsno(TOTAL_NL)
      integer n, no
      real*8 total_dz, ddz, delta, fract
      real*8, optional :: tr_flux(0:TOTAL_NL), dt
      integer i

c$$$      if( fract_cover_new .lt. MIN_FRACT_COVER ) then
c$$$        print *, 'snow_redistr: fract_cover_new < MIN_FRACT_COVER'
c$$$        print *, 'fract_cover_new = ', fract_cover_new
c$$$        call stop_model(
c$$$     &       'snow_redistr: fract_cover_new < MIN_FRACT_COVER',255)
c$$$      endif

      if ( dz(1) == 0.d0 ) return  ! no snow - nothing to redistribute

      dzo(1:nl)  = dz(1:nl)
      wsno(1:nl) = wsn(1:nl)
      hsno(1:nl) = hsn(1:nl)
      nlo = nl

      total_dz = sum( dzo(1:nl) )

      !fract_cover_ratio = fract_cover/fract_cover_new
      total_dz = total_dz*fract_cover_ratio

      if ( total_dz .gt. MIN_SNOW_THICKNESS*1.5d0 ) then
        nl = TOTAL_NL
        dz(1) = MIN_SNOW_THICKNESS
        dz(2:nl) = (total_dz-dz(1))/(nl-1)
      else
        nl = 1
        dz(1) = total_dz
      endif

      wsn(1:TOTAL_NL) = 0.d0
      hsn(1:TOTAL_NL) = 0.d0

      no = 0
      delta = 0.d0
      do n=1,nl
        do while ( delta.lt.dz(n) .and. no.lt.nlo )
          no = no + 1
          delta = delta + dzo(no)*fract_cover_ratio
          wsn(n) = wsn(n) + wsno(no)*fract_cover_ratio
          hsn(n) = hsn(n) + hsno(no)*fract_cover_ratio
          enddo
        ddz = delta - dz(n)
        fract = ddz/(dzo(no)*fract_cover_ratio)
ccc the following is just for check
          if ( fract.lt.-EPS .or. fract.gt.1.d0+EPS ) then
            print *, 'fract= ', fract
            call stop_model('snow_redistr: internal error 3',255)
            endif
        wsn(n) = wsn(n) - fract*wsno(no)*fract_cover_ratio
        hsn(n) = hsn(n) - fract*hsno(no)*fract_cover_ratio
        if ( n.lt.nl ) then
          wsn(n+1) = wsn(n+1) + fract*wsno(no)*fract_cover_ratio
          hsn(n+1) = hsn(n+1) + fract*hsno(no)*fract_cover_ratio
          delta = ddz
        else
          if ( abs(fract).gt.EPS ) then
            print *, 'fract= ', fract
            call stop_model('snow_redistr: internal error 1',255)
            endif
          endif
        enddo

ccc for tracers
      if ( present(tr_flux) ) then
        if ( nlo < TOTAL_NL ) wsno(nlo+1:TOTAL_NL) = 0.d0
        tr_flux(0) = 0.d0
        do i=1,TOTAL_NL
          tr_flux(i) =
     &         - (wsn(i) - wsno(i)*fract_cover_ratio)/dt
     &         + tr_flux(i-1)
        enddo

        if ( abs(tr_flux(TOTAL_NL)) > 1.d-12 )
     &       call stop_model("SNOW: snow redistr flux error",255)
      endif

      return
      end subroutine snow_redistr



      subroutine snow_adv(dz, wsn, hsn, nl, 1,4
     &    srht, trht, snht, htpr, evaporation, pr, dt,
     &    t_ground, dz_ground,
     &    water_to_ground, heat_to_ground,
     &    radiation_out, snsh_dt, evap_dt, ! fb_or_fv,
     &     tr_flux )
      implicit none
!@sum  a wrapper that calles real snow_adv (introduced for debugging)
!@auth I.Aleinov
!@ver  1.0
ccc input:
      integer nl
      real*8 srht, trht, snht, htpr, evaporation
      real*8 pr, dt, t_ground, dz_ground
      real*8 snsh_dt, evap_dt !, fb_or_fv

ccc output:
      real*8 water_to_ground, heat_to_ground
      real*8 radiation_out
ccc data arrays
      real*8 dz(TOTAL_NL+1), wsn(TOTAL_NL), hsn(TOTAL_NL)

ccc tracer variables
!@var tr_flux flux of water between snow layers (>0 is down) (m/s)
      real*8 tr_flux(0:TOTAL_NL)
      real*8 wsn_o(MAX_NL), evap_o
      integer nl_o

ccc for debug
      real*8 total_energy, total_water  !, lat_evap
      integer i, retcode
      integer, save :: counter=0

      counter = counter + 1
      !rint *, counter

ccc for tracers
      wsn_o(:) = 0.d0
      nl_o = nl
      wsn_o(1:nl) = wsn(1:nl)
      evap_o = evaporation

ccc checking if the model conserves energy (part 1) (for debugging)
      total_energy = 0.d0
      total_water = 0.d0
      do i=1,nl
        total_energy = total_energy - hsn(i)
        total_water = total_water - wsn(i)
      enddo

      call snow_adv_1(dz, wsn, hsn, nl,
     &    srht, trht, snht, htpr, evaporation, pr, dt,
     &    t_ground, dz_ground,
     &    water_to_ground, heat_to_ground,
     &    radiation_out, snsh_dt, evap_dt, retcode )

!      if (fb_or_fv .le. 0.) return
ccc checking if the model conserves energy (part 2) (for debugging)
      do i=1,nl
        total_energy = total_energy + hsn(i)
        total_water = total_water + wsn(i)
      enddo
      total_energy = total_energy -
     &    (srht+trht-snht-lat_evap*evaporation+htpr)*dt
      total_energy = total_energy + heat_to_ground*dt + radiation_out*dt
      if ( abs(total_energy) .gt. 1.d0 ) then
        print*, "total energy error",i_earth, j_earth,total_energy,
     *       heat_to_ground*dt,radiation_out*dt
        call stop_model('snow_adv: total energy error',255)
      end if
      total_water = total_water
     &     - (pr - evaporation - water_to_ground)*dt
      if ( abs(total_water)/dt .gt. 1.d-15 ) then
        print*, "water cons error",i_earth, j_earth, total_water/dt
        if ( abs(total_water)/dt .gt. 1.d-10 )
     &    call stop_model('snow_adv: water conservation error',255)
      end if

c$$$    if( fr_type .lt. 1.d-6 .and. abs(total_energy) .gt. 1.d-6 ) then
c$$$      print*, "total energy error",i_earth, j_earth,total_energy
c$$$      call abort
c$$$    endif

ccc for tracers
      !!!tr_flux(0) = tr_flux(0) + pr - evaporation
      tr_flux(0) = pr - evaporation
      do i=1,TOTAL_NL
        tr_flux(i) = -(wsn(i) - wsn_o(i))/dt
     &       + tr_flux(i-1)
      enddo
ccc checking if preserve water
      if ( abs( tr_flux(TOTAL_NL) - water_to_ground ) > 1.d-15 ) then
        if ( DEB_CH == 0 )
     $       call openunit("snow_debug", DEB_CH, .false., .false.)
        write(DEB_CH,*) "snow_adv: H2O error "
     &       , abs( tr_flux(TOTAL_NL) - water_to_ground )
     &       , tr_flux(TOTAL_NL) , water_to_ground
     &       , nl, nl_o, counter, evaporation, evap_o
        !call abort
        !stop 'snow_adv: H2O error'
      endif

      return
      end subroutine snow_adv



      subroutine snow_adv_1(dz, wsn, hsn, nl, 1,17
     &    srht, trht, snht, htpr, evaporation, pr, dt,
     &    t_ground, dz_ground,
     &    water_to_ground, heat_to_ground,
     &    radiation_out, snsh_dt, evap_dt, retcode )
      implicit none
!@sum main program that does column snow physics
!@auth I.Aleinov
!@ver  1.0
ccc input:
      integer nl
      real*8 srht, trht, snht, htpr, evaporation
      real*8 pr, dt, t_ground, dz_ground
      real*8 snsh_dt, evap_dt

ccc constants: (now defined as global params)
      real*8 k_ground, c_ground
ccc output:
      real*8 water_to_ground, heat_to_ground
      real*8 radiation_out
      integer retcode

ccc main parameters: layer thickness, water equivalent, heat content
      real*8 dz(TOTAL_NL+1), wsn(TOTAL_NL), hsn(TOTAL_NL)
ccc!!! I wonder if the following arrays should have dim (nl) instead
ccc    of (MAX_NL) to force the allocation on a stack (for OpenMP) ?
      real*8 tsn(MAX_NL+1), csn(MAX_NL), ksn(MAX_NL+1), isn(MAX_NL)
      real*8 wsn_o(MAX_NL), hsn_o(MAX_NL), dz_o(MAX_NL+1)

      real*8 water_down, heat_down, fresh_snow, rho_snow, flux_in
      real*8 mass_layer, mass_above, scale_rho
      real*8 flux_in_deriv, flux_corr
      real*8 delta_tsn_impl ! shft of temperature due to implicit method
      integer n, nl_o
      real*8 dz_he(MAX_NL+1)

ccc!!!  check if lat_evap shoud be replaced by heat of sublimation

ccc !!! ground properties should be passed as formal parameters !!!
      k_ground =        3.4d0    !    
      c_ground =        1.d5     !  

ccc will need initial values if all snow melted and for tracers
      wsn_o(1:nl) = wsn(1:nl)
      hsn_o(1:nl) = hsn(1:nl)
      dz_o(1:nl) = dz(1:nl)
      nl_o = nl

      !!! just to deceive the optimizer
      if ( dz_o(1) > 1.d6 ) call stop_model("snow_adv_1: what?!!",255)

ccc just in case, may need reasonable tsn(1) if all melted
      tsn(1) = 0.d0

ccc fluxes in/out
      heat_to_ground = 0.d0
      water_to_ground = 0.d0


      call check_rho_snow( dz, wsn, nl )



ccc compute amount of fresh snow
ccc !!! insert evaporation into computation of the amount of fresh snow
ccc use fresh_snow only to change the depth of the first layer
      fresh_snow = rho_water/rho_fresh_snow
     &                * min(pr*dt-evaporation*dt, -htpr*dt/lat_fusion)
      if(fresh_snow > 0.d0) then
        dz(1) = dz(1) + fresh_snow
      else
        if ( wsn(1) < EPS ) goto 1000
        !goto 1000
      endif

ccc water and heat to pass through the snow pack
      water_down = (pr - evaporation)*dt
      heat_down = htpr*dt
ccc will include heat of evaporation later
c!!      heat_down = heat_down - lat_evap*evaporation*dt

!#ifdef 
!      call check_rho_snow( dz, wsn, nl )
!#endif

      call pass_water( wsn, hsn, dz, nl, water_down, heat_down,
     &                       lat_fusion, max_fract_water,
     &                       rho_water, rho_fresh_snow)
      heat_to_ground = heat_to_ground + heat_down/dt
      water_to_ground = water_to_ground + water_down/dt

      if ( sum( wsn(1:nl) ) < EPS ) goto 1000


      call check_rho_snow( dz, wsn, nl )


      call snow_redistr(dz, wsn, hsn, nl, 1.d0)
      dz(nl+1) = dz_ground

ccc compute spec. heat and thermal conductivity
      do n=1,nl
        rho_snow = wsn(n)*rho_water/dz(n)
        !csn(n) = (1.9d6) * rho_snow / rho_ice
        csn(n) = shi_kg * rho_snow
        ksn(n) = (3.22d-6) * rho_snow**2
      enddo
      ksn(nl+1) = k_ground

ccc compute temperature of the layers (and amount of ice)
      do n=1,nl
        if ( hsn(n) .gt. 0.d0 ) then
          print *, 'OOPS, No snow in layer ', n
          call stop_model('snow_adv_1: empty snow layer found (1)',255)
        else if ( hsn(n) .gt. -wsn(n)*lat_fusion ) then
          tsn(n) = 0.d0
          isn(n) = -hsn(n)/lat_fusion
        else
          tsn(n) = (hsn(n)+wsn(n)*lat_fusion)/(csn(n)*dz(n))
          isn(n) = wsn(n)
        endif
      enddo
      tsn(nl+1) = t_ground

c!!! this is for debugging
      if(tsn(1).lt.-120.d0) call stop_model("SNOW:tsn<-120",255)

ccc compute incomming heat flux (from atm.)
ccc include all fluxes except htpr (which is already included)
      flux_in =
     &      srht
     &      +trht - sigma*(tsn(1)+tfrz)**4
     &      -lat_evap*evaporation
     &      -snht

      flux_in_deriv =
     &     -4.d0*sigma*(tsn(1)+tfrz)**3 - lat_evap*evap_dt - snsh_dt

      radiation_out = sigma*(tsn(1)+tfrz)**4

ccc this is a hack to keep heat_eq stable: force dz >= 0.1
      !dz_he(1:nl+1) = dz(1:nl+1)
      !dz_he(1) = max( dz_he(1), MIN_SNOW_THICKNESS )

ccc solve heat equation
      call heat_eq(
     &     dz, tsn, hsn, csn, ksn, nl,
     &     flux_in, flux_in_deriv, flux_corr, dt )

c!!! this is for debugging
      if(tsn(1).lt.-120.d0) call stop_model("SNOW:he:tsn<-120",255)

      heat_to_ground = heat_to_ground + flux_in

ccc now redistribute extra flux proportionallly to input derivatives
ccc snht_out_cor, delta_tsn_impl
      delta_tsn_impl = flux_corr / flux_in_deriv
      radiation_out = radiation_out -
     &     (-4.d0*sigma*(tsn(1)+tfrz)**3)*delta_tsn_impl
      snht = snht + snsh_dt * delta_tsn_impl
      evaporation = evaporation + evap_dt * delta_tsn_impl

ccc and now remove (add) water due to extra evaporation.
c!!! this may make wsn(1) negative, the only way I see now to prevent it
c!!! is to keep minimal thickness of snow big enough

      water_down = - evap_dt * delta_tsn_impl * dt ! ??


      call check_rho_snow( dz, wsn, nl )


ccc pass extra water down
ccc      water_down = 0.d0
      heat_down = 0.d0
      call pass_water( wsn, hsn, dz, nl, water_down, heat_down,
     &                       lat_fusion, max_fract_water,
     &                       rho_water, rho_fresh_snow)
      heat_to_ground = heat_to_ground + heat_down/dt
      water_to_ground = water_to_ground + water_down/dt

      if ( sum( wsn(1:nl) ) < EPS ) goto 1000

      call snow_redistr(dz, wsn, hsn, nl, 1.d0)
      dz(nl+1) = dz_ground

ccc compute temperature of the layers
      do n=1,nl
        if ( hsn(n) .gt. 0.d0 ) then
          print *, 'OOPS, No snow in layer ', n
          call stop_model('snow_adv_1: empty snow layer found (2)',255)
        else if ( hsn(n) .gt. -wsn(n)*lat_fusion ) then
          tsn(n) = 0.d0
        else
          tsn(n) = (hsn(n)+wsn(n)*lat_fusion)/(csn(n)*dz(n))
        endif
      enddo
      tsn(nl+1) = t_ground

ccc repack the layers
      mass_above = 0.d0
      do n=1,nl
        if( dz(n) .gt. EPS ) then
          mass_layer = wsn(n)*rho_water
          mass_above = mass_above + .5d0 * mass_layer
          scale_rho = .5d-7 * 9.8d0 * mass_above
     &     * exp(14.643d0 - 4000.d0/(tsn(n)+tfrz)
     &                  - .02d0 * mass_layer/dz(n))
     &     * dt
          scale_rho = 1.d0 + scale_rho
          dz(n) = dz(n) / scale_rho
          if( dz(n).lt.mass_layer/rho_ice )
     &                   dz(n) = mass_layer/rho_ice
          mass_above = mass_above + .5d0 * mass_layer
        endif
      enddo


      call check_rho_snow( dz, wsn, nl )


c!!! this is for debugging
      if(tsn(1).lt.-120.d0) then
        print*,"tsn error",i_earth, j_earth,1,tsn(1:nl)
        call stop_model('snow_adv_1: tsn error',255)
      end if

      if ( radiation_out > 500.d0 ) call stop_model("SNOW:T>0",255)

      retcode = 0
      return

 1000 continue ! all melted
ccc!!! may create water_to_ground<0 if evaporation is not properly
ccc    limited. Check later.
      water_to_ground = sum( wsn_o(1:nl_o) )/dt + pr - evaporation
      heat_to_ground  = sum( hsn_o(1:nl_o) )/dt + htpr
     &     - lat_evap*evaporation
     &     - snht + srht + trht - sigma*(tsn(1)+tfrz)**4
      radiation_out = sigma*(tsn(1)+tfrz)**4
      if ( radiation_out > 500.d0 ) call stop_model("SNOW:T>0",255)
      wsn(1:nl) = 0.d0
      hsn(1:nl) = 0.d0
      nl = 1
      return
      end subroutine snow_adv_1



      subroutine heat_eq( 1,1
     &     dz, tsn, hsn, csn, ksn, nl,
     &     flux_in, flux_in_deriv, flux_corr, dt)
      implicit none
!@sum solves heat transport equation for snow layers
!@auth I.Aleinov
!@ver  1.0
!@alg This solver is currently using a half-implicit scheme.
!@+   Implicitness is controlled by the parameters:
!@+       gamma - for upper boundary
!@+       eta(MAX_NL+1) - for the rest of the boundaries
!@+   In general we are trying to use half-implicit method (gamma = .5)
!@+   at the surface. But this may introduce a systematic error when
!@+   temperature of the snow is close to 0 C.
!@+   When option DO_EXPLIC_0 is enabled the program checks for possible
!@+   error and if necessary recomputes the solution with reduced gamma,
!@+   i.e. in more explicit way.
!@+   DO_EXPLIC_0 is enabled by default. One may try to turn it off if
!@+   there are serious stability problems.
!@calls :TRIDIAG
      integer nl
      real*8 dz(nl+1), tsn(nl+1), hsn(nl)
      real*8 csn(nl), ksn(nl+1)
ccc n+1 corresponds to the first layer of the ground
      real*8 flux_in, flux_in_deriv, flux_corr, dt, syst_flux_err

      real*8 eta(MAX_NL+1), tnew(MAX_NL), gamma
ccc arrays for coefficients:
      real*8 a(MAX_NL), b(MAX_NL), c(MAX_NL), f(MAX_NL)
      real*8 left, right, dt_to_cdz
      integer n, iter
      integer :: itermax=2

ccc setting implicit/explicit choice for each layer
ccc eta=1 - implicit, eta=0 - explicit
ccc for now using explicit method for layers with T=0.d0

      do n=1,nl
         eta(n) = .5d0
ccc        eta(n) = 1.d0
c!!        if( tsn(n) .ge. 0.d0 ) eta(n) = 0.d0
        enddo
      eta(nl+1) = 0.d0
ccc      gamma = .5d0
c!! I am trying to make it nearly expicit for a test !!!
      gamma = .5d0

!!! testing: trying to improve stability
ccc switching to explicit method for very thin snow
      if ( dz(1) < MIN_SNOW_THICKNESS * .5d0 ) then
        eta(1) = 1.d0
        gamma = 1.d0
      endif

ccc In general we are trying to use half-implicit method (gamma = .5d0)
ccc on the interface. But this introduces a systematic error when
ccc tnew(1) > 0. So in the case DO_EXPLIC_0 we move to explicit method
ccc (i.e. reduce gamma)

ccc DO_EXPLIC_0 enables correction of the error which is introduced
ccc by implicit scheme when tsn(1) == 0 C
ccc (of course one has to use preprocessing to do it ...)
ccc it is enabled by default



      do iter=1,itermax

ccc equation for upper snow layer
        n = 1
        dt_to_cdz = dt/(csn(n)*dz(n))
        right = 2.d0*dt_to_cdz/( dz(n)/ksn(n) + dz(n+1)/ksn(n+1) )
        a(n) = 1.d0 + right*eta(n) - dt_to_cdz*flux_in_deriv*gamma
        b(n) = 0.d0
        c(n) = -right*eta(n+1)
        f(n) = tsn(n)*( 1.d0 - right*(1.d0-eta(n))
     &                  - dt_to_cdz*flux_in_deriv*gamma )
     &         + tsn(n+1)*right*(1.d0-eta(n+1))
     &         + dt_to_cdz * flux_in

ccc all other snow layers including the lower one
        do n=2,nl

          dt_to_cdz = dt/(csn(n)*dz(n))
          right = 2.d0*dt_to_cdz/( dz(n)/ksn(n) + dz(n+1)/ksn(n+1) )
          left  = 2.d0*dt_to_cdz/( dz(n)/ksn(n) + dz(n-1)/ksn(n-1) )

          a(n) = 1.d0 + ( left + right )*eta(n)
          b(n) = -left*eta(n-1)
          c(n) = -right*eta(n+1)

          f(n) = tsn(n)*( 1.d0 - ( left + right )*(1.d0-eta(n)) )
     &         + tsn(n-1)*left*(1.d0-eta(n-1))
     &         + tsn(n+1)*right*(1.d0-eta(n+1))
          enddo

c        call sweep3diag(b, c, a, f, tnew, nl)
          call TRIDIAG(b,a,c,f,tnew,nl)

ccc flux_corr is the energy wich should be returned to the atmosphere
        flux_corr = flux_in_deriv*( tnew(1) - tsn(1) )*gamma


        syst_flux_err = flux_in_deriv*( tnew(1) - 0.d0 )*gamma
        ! if( iter/=2 .and. tnew(1)>0.d0 .and. flux_in_deriv<0.d0 ) then
        ! back to 77 :-L
        if ( iter.ne.itermax .and.
     &         tnew(1).gt.0.d0 .and. flux_in_deriv.lt.0.d0 ) then
          syst_flux_err = flux_in_deriv*( tnew(1) - 0.d0 )*gamma
          gamma = (1.d0 - syst_flux_err/flux_corr) *gamma
        else
          ! exit
          ! back to 77 :-L
          goto 77
        endif
      enddo
 77   continue


      do n=1,nl
        hsn(n) = hsn(n) + (tnew(n)-tsn(n))*csn(n)*dz(n)
        enddo

ccc flux to the ground :
      n = nl
      flux_in = - ( tsn(n+1)-tnew(n)*eta(n)-tsn(n)*(1.d0-eta(n)) ) *
     &           2.d0/( dz(n)/ksn(n) + dz(n+1)/ksn(n+1) )

      return
      end subroutine heat_eq


      subroutine check_rho_snow( dz, wsn, nl ) 4,2
      integer nl
      real*8 dz(nl+1), wsn(nl)
      integer n
      do n=1,nl
        if( wsn(n) < EPS ) cycle
        if ( wsn(n)/dz(n)*rho_water+EPS < rho_fresh_snow) then
          print *,"rho_snow < rho_fresh_snow at i,j=",i_earth, j_earth
          print *,"n, wsn, rho_sn, rho_fresh= ", n,wsn(n),
     *         wsn(n)/dz(n)*rho_water ,rho_fresh_snow
          call stop_model('check_rho_snow: rho < rho_fresh_snow',255)
        end if
      enddo

      do n=1,nl
        if( wsn(n) < EPS ) cycle
        if ( wsn(n)/dz(n)*rho_water-EPS > rho_ice) then
          print *,"rho_snow > rho_ice at i,j=",i_earth, j_earth
          print *,"n, wsn, rho_sn, rho_ice= ", n,wsn(n),
     *         wsn(n)/dz(n)*rho_water ,rho_ice
          call stop_model('check_rho_snow: rho > rho_ice',255)
        end if
      enddo


      end subroutine check_rho_snow

      END MODULE SNOW_MODEL