!-------------------------------------------------------------------------------
! Snow module (snow)
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
! Author : Shaochun Huang
! Version : 0.1
! Edited by: stefan.liersch@pik-potsdam.de
! Modified: 2015-02-07:
!-------------------------------------------------------------------------------
module snow

  use utilities, only : dp, path_max_length

  implicit none

  ! threshold temperature for snow fall
  real(dp), save, dimension(:), allocatable :: bsn_tsnfall
  ! threshold temperature for snow melt
  real(dp), save, dimension(:), allocatable :: bsn_tmelt
  ! snow melt rate for the degree - day method
  real(dp), save, dimension(:), allocatable :: bsn_smrate
  ! glacier melt rate for the degree - day method
  real(dp), save, dimension(:), allocatable :: bsn_gmrate
  ! snow melt, calc in snom(), mm
  real(dp), save :: sml
  ! precipitation as snow, mm
  real(dp), save :: snowVal
  ! initial snow content, mm (from .bsn)
  real(dp), save :: snow1
  real(dp), save :: prcor = 1. !
  ! threshold temperature for snow fall
  real(dp), save, dimension(:), allocatable :: tsnfall
  real(dp) :: tsnfall0
  ! threshold temperature for snow melt
  real(dp), save, dimension(:), allocatable :: tmelt
  real(dp), save :: tmelt0
  ! snow melt rate for the degree - day method (cm/dg/day)
  real(dp), save, dimension(:), allocatable :: smrate
  real(dp) :: smrate0 = 1.0
  ! glacier melt rate for the degree - day method (mm/dg/day)
  real(dp), save, dimension(:), allocatable :: gmrate
  real(dp) :: gmrate0 = 10.
  ! daily max temp. for subbasin, readcli, degree C
  real(dp), save, dimension(:), allocatable :: tmx
  ! elevation centroids, m
  real(dp), save, dimension(:), allocatable :: elev0
  ! snow water content in HYDROTOPE, mm
  real(dp), save, dimension(:, :), allocatable :: snoa
  ! accumulated snow water content in the catchment at the end of the year, mm
  real(dp), save :: snow_acc_mm
  ! accumulated snow water content in the catchment at the end of the previous year, mm
  real(dp), save :: snow_acc_mm0
  ! accumulated glacier water content in the catchment at the end of the year, mm
  real(dp), save :: glac_acc_mm
  ! accumulated glacier water content in the catchment at the end of the previous , mm
  real(dp), save :: glac_acc_mm0
  ! accumulated value of snow cover, mm / water content of snow, mm
  real(dp), save, dimension(:), allocatable :: sno


  !#############################################################################
  !
  ! VARIABLES & PARAMETERS
  !
  !#############################################################################

  logical, save :: bSnowModule = .true.

  ! ! Day when snow is converted into ice and output is written into GIS file
  ! This might be different in different regions
  integer, parameter :: gla_day_out = 273
  ! ! Daily factor used to convert accumulated snow (snoa) to glacier ice, dimensionless.
  ! !real(dp), parameter :: snow2ice = 0.001
  ! real(dp), parameter :: snow2ice = 1./365.25
  !
  ! ! Total glacier melt (gla_melt_bottom = melt from snow on top of the glacier), mm.
  ! real(dp), save :: gla_melt_total
  !
  ! mean T, °C
  real(dp), save :: tmit
  ! mean T, °C
  real(dp), save :: tx_tmp
  ! Max T, °C
  real(dp), save :: tmax
  ! mean T, °C
  real(dp), save :: tmax_tmp
  ! Min T, °C
  real(dp), save :: tmin
  ! mean T, °C
  real(dp), save :: tmin_tmp
  ! Mean precipitation, mm
  real(dp), save :: precipe
  ! elevation - based corrected HRU precipitation, mm
  real(dp), save :: precip_elev_cor
  ! snow depth, cm
  real(dp), save :: hsn0
  ! snow melt, mm
  real(dp), save :: smle
  ! glacier melt, mm
  real(dp), save :: gmle

  ! weighted subbasin precipitation, mm
  real(dp), save :: xprecip
  ! weighted subbasin elevation - based corrected precipitation, mm
  real(dp), save :: xprecip_elev_cor
  ! weighted subbasin max T, °C
  real(dp), save :: xtmax
  ! weighted subbasin mean T, °C
  real(dp), save :: xtmit
  ! weighted subbasin min T, °C
  real(dp), save :: xtmin
  ! weighted subbasin snow melt, mm
  real(dp), save :: xsml
  ! weighted subbasin snow, mm
  real(dp), save :: xsnow
  ! weighted subbasin snow fall as precpitation, mm
  real(dp), save :: psnow
  ! hydrotope snow melt from snow pack, mm
  real(dp), save :: vsn
  ! weighted subbasin snow melt from snow pack, mm
  real(dp), save :: xvsn

  !real(dp), save :: gmrate ! glacier melt rate for the degree-day method NOW DEFINED IN COMMON.F90
  ! precipitation factor
  real(dp), save :: xgrad1
  ! tempreture factor
  real(dp), save :: tgrad1
  ! holding capacity
  real(dp), save :: ulmax0 = 1.
  ! fresh snow density
  real(dp), save :: rnew = 0.08
  ! snow water density in HYDROTOPE, (g / cm3)
  real(dp), save, dimension(:, :), allocatable :: rsn
  ! volumetric content of liquid water in snow (g / cm3)
  real(dp), save, dimension(:, :), allocatable :: sul
  ! volumetric content of ice in snow (g / cm3)
  real(dp), save, dimension(:, :), allocatable :: suz
  ! glacier water equivalent, mm
  real(dp), save, dimension(:, :), allocatable :: gla


  ! for output only, otherwise local variables
  real(dp), save :: balanc, hsn, ulmax
  ! counter for GIS output
  ! for glacier output (snow module)
  integer, save :: ieapg = 0
  !#############################################################################

  !*****************************************************************************
  ! Output variable ids
  integer :: &
    glacier_weq_output_id = 0, &
    snowfall_weq_output_id = 0, &
    snow_depth_weq_output_id = 0


  namelist / SNOW_PARAMETERS / &
    tsnfall0, &
    smrate0, &
    gmrate0, &
    prcor, &
    snow1, &
    xgrad1, &
    tgrad1, &
    ulmax0, &
    rnew, &
    tmelt0, &
    bSnowModule

contains

  subroutine snow_initialise(mb, meap, subbasin_input_file_id, mstruc)
    use input, only: get_config_fid
    use output, only: output_register_hydrotope_var, output_register_subbasin_var
    integer, intent(in) :: mb, meap, subbasin_input_file_id
    integer, dimension(:, :, :), intent(in) :: mstruc

    read(get_config_fid(), SNOW_PARAMETERS)
    glacier_weq_output_id = output_register_hydrotope_var("glacier_weq", .False.)
    snowfall_weq_output_id = output_register_subbasin_var("snowfall_weq")
    snow_depth_weq_output_id = output_register_hydrotope_var("snow_depth_weq")
    call snow_allocate(mb, meap)
    call snow_read_input(subbasin_input_file_id)

    if (bSnowModule) then
      allocate(rsn(mb, meap))
      rsn = 0.
      allocate(sul(mb, meap))
      sul = 0.
      allocate(suz(mb, meap))
      suz = 0.
      allocate(gla(mb, meap))
      gla = 0.
      gla = mstruc(:, :, 6) ! initialise glacier depth [mm]
      snoa = snow1
    end if
  end subroutine snow_initialise


  subroutine snow_read_input(subbasin_input_file_id)
    use input, only : read_real_column
    integer, intent(in) :: subbasin_input_file_id

    call read_real_column(subbasin_input_file_id, "sno", sno, snow1)
    call read_real_column(subbasin_input_file_id, "elev0", elev0, 0.0_dp)

  end subroutine snow_read_input


  subroutine snow_allocate(mb, meap)
    integer, intent(in) :: mb, meap
    allocate(elev0(mb))
    allocate(gmrate(mb))
    allocate(smrate(mb))
    allocate(sno(mb))
    allocate(snoa(mb, meap))
    allocate(tmelt(mb))
    allocate(tmx(mb))
    allocate(tsnfall(mb))
    elev0 = 0.
    gmrate = gmrate0
    smrate = smrate0
    sno = 0.
    snoa = 0.
    tmelt = tmelt0
    tmx = 0.
    tsnfall = tsnfall0
  end subroutine snow_allocate

  subroutine snow_allocate_subcatch(n_subcatch)
    integer, intent(in) :: n_subcatch
    allocate(bsn_tsnfall(n_subcatch))
    bsn_tsnfall = 0.
    allocate(bsn_tmelt(n_subcatch))
    bsn_tmelt = 0.
    allocate(bsn_smrate(n_subcatch))
    bsn_smrate = 0.
    allocate(bsn_gmrate(n_subcatch))
    bsn_gmrate = 0.
  end subroutine snow_allocate_subcatch

  subroutine dealloc_snow
   deallocate(sno)
   deallocate(snoa)
   deallocate(tmx)
  end subroutine dealloc_snow


  subroutine snow_initialise_subbasin
    xprecip = 0.
    xprecip_elev_cor = 0.
    xtmax = 0.
    xtmit = 0.
    xtmin = 0.
    xsml = 0.
    xsnow = 0.
    psnow = 0.
    xvsn = 0.
  end subroutine snow_initialise_subbasin
  !------------------------------------------------------------------------------

  !------------------------------------------------------------------------------
  subroutine snow_process(j, jea, ida, mstruc, precip, tmn, tx)
    use output, only: output_store_hydrotope_value

    integer, intent(in) :: ida
    integer, dimension(:, :, :), intent(in) :: mstruc
    real(dp), intent(in) :: precip
    real(dp), dimension(:), intent(in) :: tmn
    real(dp), dimension(:), intent(in) :: tx
    integer, intent(in) :: j, jea
    real(dp) :: elevh

    smle = 0.
    snowVal = 0.
    gmle = 0.
    elevh = mstruc(j, jea, 5) ! hydrotope elevation
    vsn = 0.

    ! elevation-based correction of precipitation and air temperature
    if (precip > 0.) then
      precipe = precip * (1. + xgrad1 * (elevh - elev0(j)))
    else
      precipe = 0.
    end if
    precipe = max(0.0_dp, precipe)
    precip_elev_cor = precipe

    tmit = tx(j) + tgrad1 * (elevh - elev0(j))
    tmax = tmx(j) + tgrad1 * (elevh - elev0(j))
    tmin = tmn(j) + tgrad1 * (elevh - elev0(j))

    ! rsn = snow water density in HYDROTOPE, (g/cm3)
    ! snoa = snow water content, mm
    if (rsn(j, jea) > 0. ) then
      hsn0 = snoa(j, jea) * .1 / rsn(j, jea) ! hsn0 = snow depth, cm
    else
      hsn0 = 0.
    end if

    if (tmit <= tsnfall(j) ) then
      snoa(j, jea) = snoa(j, jea) + precipe * prcor! snow ! add precipitation to snow water content, mm
      snowVal = precipe * prcor ! actual snow fall, mm
      precipe = 0.
    end if

    if (snoa(j, jea) > 0. ) then
      call snow_melting(j, jea, ida)
    else
      if (gla(j, jea) > 0. ) call snow_glacier_melt(j, jea)
    end if

    if (ida == gla_day_out) &
      call output_store_hydrotope_value(glacier_weq_output_id, j, jea, gla(j, jea))

  end subroutine snow_process
  !------------------------------------------------------------------------------

  !------------------------------------------------------------------------------
  subroutine snow_melting(j, jea, ida)
    integer, intent(in) :: ida
    integer, intent(in) :: j, jea
    real(dp) :: ph, smli, sml0, ulmax, tsn, sfreez, svz, svl, snoa0
    real(dp) :: precipe0, smr

    sfreez = 0.
    tsn = 0.
    ph = 0.
    svz = 0.
    svl = 0.
    smli = 0.
    sml0 = 0.
    hsn = 0.
    vsn = 0.
    !    gla_melt_total = 0.
    !#### snow pack depth, hsn, hsn0 (last step), cm
    ! snow, mm
    ! hsn0, snow depth, cm
    ! hsn, cm
    hsn = hsn0 + snowVal * .1 / rnew ! actual snow depth, cm
    snoa0 = hsn0 * rsn(j, jea) * 10. ! previous day snow water content, mm
    precipe0 = precipe ! current day precipitation

    !#### different snowmeltrate for forest and opensite (cm/degree)
    smr = smrate(j) ! cm / °C

    !#### snow temperature
    if (tmit < 0.) tsn = tmit
    !#### calculate ulmax (holding capacity of snowpack)
    ulmax = ulmax0 ! 0 or 1
    if (rsn(j, jea) /= 0. .AND. ulmax0 == 1.) ulmax = - 0.11 * rsn(j, jea) + 0.11
    if (rsn(j, jea) == 0. .AND. ulmax0 == 1.) ulmax = 0.
    ulmax = max(0.0_dp, ulmax)

    !#### calculate intensity of snow melt (cm/day)
    if (tmit > tmelt(j) ) then
      smli = smr * (tmit - tmelt(j)) ! cm (cm / °C * °C )
    end if

    !#### calculate intensity of freezing of liquid water in snowpack (cm/day)
    if (tmit < 0.) sfreez = smr * sqrt(- tmit) ! cm / day
    if (sfreez >= (sul(j, jea) * hsn0) ) sfreez = sul(j, jea) * hsn0 ! is actually g / cm2, similar to cm / day?

    !#### calculate the compression sml0 (cm/day)
    !#### the snowmelt water percolates into pack, hits the ice, and freezes until
    !#### the temperature of the pack is 0.
    if (suz(j, jea) > 0.) then
      if (hsn0 > 0.) then
          ph = .15 * rsn(j, jea) * hsn0 ** 2. / exp(- 0.08 * tsn + 21. * rsn(j, jea))
      end if
      sml0 = smli / (suz(j, jea) * 0.917) + ph
      if (sml0 > hsn) sml0 = hsn
      hsn = hsn - sml0
    else
      smli = 0.
      hsn = snowVal * 0.1 / rnew
    end if

    if (hsn .ne. 0.) then ! SL: equals or lower equal?
      !#### recalculate suz (volumetric content of ice in snow g/cm3)
      svz = 0.917 * suz(j, jea) * hsn0 + snowVal * 0.1 + sfreez
      if (smli .gt. svz) smli = svz
      suz(j, jea) = (svz - smli) / (hsn * 0.917)

      !#### recalculate sul (volumetric content of liquid water in snow g/cm3, svl cm)
      svl = sul(j, jea) * hsn0 + precipe * 0.1 + smli
      if (sfreez > svl) sfreez = svl
      sul(j, jea) = (svl - sfreez) / hsn

      !#### water outflow from snowpack
      if (sul(j, jea) > ulmax ) then
        vsn = hsn * 10. * (sul(j, jea) - ulmax)
        if (vsn > precipe0) smle = vsn - precipe0
        sul(j, jea) = ulmax
      end if
    else
      sul(j, jea) = 0.
      suz(j, jea) = 0.
      vsn = rsn(j, jea) * hsn0 * 10. + precipe0
      smle = rsn(j, jea) * hsn0 * 10.
    end if

    !#### recalculate rsn
    rsn(j, jea) = sul(j, jea) + 0.917 * suz(j, jea)
    if (hsn == 0.) rsn(j, jea) = 0.

    !#### recalculate the total precipitation from snowpack
    precipe = vsn
    snoa(j, jea) = hsn * 10. * rsn(j, jea)
    balanc = snowVal + precipe0 + (snoa0 - snoa(j, jea)) - vsn ! SL: balanc is never used hereafter

    !#### conversion from snow to ice
    !   if(ida.eq.273) then
    if (ida == gla_day_out) then
      ! SL: Maybe one should not convert all the snow into ice at day gla_day_out
      !     but every day a fraction only (snow2ice) !
      !if ( snoa(j, jea) > 2000. ) then
      gla(j, jea) = gla(j, jea) + snoa(j, jea) ! * snow2ice
      snoa(j, jea) = 0. !snoa(j, jea) * (1. - snow2ice)
      hsn = 0. !snoa(j, jea) * .1 / rsn(j, jea)
      rsn(j, jea) = 0.
    end if
  !end if

  end subroutine snow_melting
  !------------------------------------------------------------------------------

  !------------------------------------------------------------------------------
  subroutine snow_glacier_melt(j, jea)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES DAILY ICE MELT
    !              WHEN AVERAGE AIR TEMP EXCEEDS 0 DEGREES C
    !**** CALLED IN: SUBBASIN
    !     snow melt rate in SWAT = 4.57: sml = 4.57 * tmax
    !     snow melt rate in HBV = 3.2
    !     current version in SWIM: sml = 3.2 * tx(j
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      precipe = precipitation in the hydrotop, mm
    !      gmle = ice melt, mm
    !      gla(j, jea) = glacier water equivalent, mm
    !      tmax = daily max temp, read in readcli, degree C
    !      tmit = daily average temp, read in readcli, degree C
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      gmr = ice melt rate, coef for the degree-day method, mm/degree
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    integer, intent(in) :: j, jea

    if (tmit > tmelt(j) ) then
      gmle = gmrate(j) * (tmit - tmelt(j))
      if (gmle > gla(j, jea)) gmle = gla(j, jea)
      gla(j, jea) = gla(j, jea) - gmle
      precipe = precipe + gmle
    else
      gmle = 0.
    end if

  end subroutine snow_glacier_melt

  subroutine snow_degree_day_melting(j, jea)
    !**** PURPOSE: THIS SUBROUTINE COMPUTES DAILY SNOW MELT
    !              WHEN AVERAGE AIR TEMP EXCEEDS 0 DEGREES C
    !**** CALLED IN: HYDROTOPE
    !     snow melt rate in SWAT = 4.57: sml = 4.57 * tmx(j)
    !     snow melt rate in HBV = 3.2
    !     current version in SWIM: sml = 3.2 * tx(j
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    !     PARAMETERS & VARIABLES
    !
    !      >>>>> COMMON PARAMETERS & VARIABLES
    !      sml = snow melt, mm
    !      smrate = snow melt rate for the degree-day method
    !      sno(j) = water content of snow, mm
    !      tmelt = threshold temperature for snow melt
    !      tmx(j) = daily max temp, read in readcli, degree C
    !      tx(j) = daily average temp, read in readcli, degree C
    !      >>>>>

    !      >>>>> STATIC PARAMETERS
    !      smr = snow melt rate, coef for the degree-day method, mm/degree
    !      >>>>>
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    !**** Include common parameters
    !     T.Vetter modification: smrate from .bsn file; tmelt instead of 0;
    !     T.Vetter modification: sml = smr * (tmx(j)-tmelt)

    integer j, jea
    real(dp) smr

    smr = smrate(j)
    sml = 0.
    !if (tmx(j).gt.tmelt) then
    sml = smr * (tmx(j) - tmelt(j))
    if (sml .gt. snoa(j, jea)) sml = snoa(j, jea)
    snoa(j, jea) = snoa(j, jea) - sml
    !end if
    return
  end subroutine snow_degree_day_melting



!------------------------------------------------------------------------------
end module snow
!------------------------------------------------------------------------------