suews_phys_snow.f95

Go to the documentation of this file.

MODULE snow_module
   USE evap_module, ONLY: cal_evap
   USE allocatearray, ONLY: nsurf, pavsurf, bldgsurf, conifsurf, bsoilsurf, watersurf, ncolumnsdataoutsnow
 
   IMPLICIT NONE
 
CONTAINS
   !This subroutine makes snow related calculations at the model time step. Needed for the
   !available energy in LUMPS and SUEWS. Made by LJ in Dec 2012
   !SUBROUTINES:
   !  MeltHeat - Calculation of snow related energy processes
   !  SnowCalc - Calculation of snow and soil storages
   !  Evap_SUEWS_Snow - Calculation of evaporation from the SnowPack
   !  snowRem - Removal of snow my snow clearing
   !  SnowDepletionCurve - Calculation of snow fractions
   !Last modified
   !  TS 22 Sep 2019 - split `SnowUpdate` into `albedo` and `density` parts so they can be better integrated with `SUEWS_cal_Qn` and `Snow_cal_MeltHeat`
   !  TS 17 Sep 2017 - added wrapper `Snow_cal_MeltHeat` for `SUEWS_driver`
   !  TS 04 Sep 2017 - added `veg_fr_snow` to update VegFractions with snow effect included
   !  TS 31 Aug 2017 - fixed the incomplete explicit interfaces
   !  LJ 24 Aug 2017 - added explicit interfaces
   !  LJ 3 May 2016  - Changed so that not all surface water freezes in 5-min timestep.
   !                    Re-organization of the snow routine due to this change
   !                    Calculation of albedo moved from MeltHeat to SnowCalc
   !  LJ 27 Jan 2016  - Tabs removed, cleaning of the code
   !  HCW 08 Dec 2015 - Added check for no Paved surfaces
   !  LJ 14 July 2015 - Code fixed to work with tstep.
   !  HCW 06 Mar 2015 - Unused variable 'i' removed.
   !  LJ Jan 2015     - Change the calculation from hourly timestep to timestep defined by nsh
   !  LJ May 2013     - Calculation of the energy balance for the SnowPack was modified
   !                        to use qn1_ind_snow(StoreDrainPrm)
   !=======================================================================================
   ! SUBROUTINE Snow_cal_MeltHeat( &
   !    tstep, tau_r, SnowDensMax, & !input
   !    lvS_J_kg, lv_J_kg, RadMeltFact, TempMeltFact, SnowAlbMax, &
   !    SnowDensMin, Temp_C, Precip, PrecipLimit, PrecipLimitAlb, &
   !    nsh_real, sfr_surf, Tsurf_ind, Tsurf_ind_snow, state_surf, qn1_ind_snow, &
   !    kup_ind_snow, SnowWater, deltaQi, albedo_snow, &
   !    SnowPack_in, SnowFrac_in, SnowAlb_in, SnowDens_in, SnowfallCum_in, & !input
   !    SnowPack_out, SnowFrac_out, SnowAlb_out, SnowDens_out, SnowfallCum_out, & !output
   !    mwh, Qm, QmFreez, QmRain, & ! output
   !    snowCalcSwitch, Qm_melt, Qm_freezState, Qm_rain, FreezMelt, &
   !    FreezState, FreezStateVol, rainOnSnow, SnowDepth, mw_ind, &
   !    dataOutLineSnow) !output
 
   !    IMPLICIT NONE
   !    ! INTEGER, PARAMETER::nsurf = 7
   !    ! INTEGER, PARAMETER::PavSurf = 1
   !    ! INTEGER, PARAMETER::BldgSurf = 2
   !    ! INTEGER, PARAMETER::WaterSurf = 7
   !    INTEGER, PARAMETER :: ncolumnsDataOutSnow_notime = ncolumnsDataOutSnow - 5
   !    REAL(KIND(1D0)), PARAMETER :: waterDens = 999.8395 !Density of water in 0 cel deg
 
   !    !These are input to the module
   !    ! INTEGER, INTENT(in) :: SnowUse
   !    INTEGER, INTENT(in) :: tstep
   !    ! INTEGER,INTENT(in)::bldgsurf
   !    ! INTEGER,INTENT(in)::nsurf
   !    ! INTEGER,INTENT(in)::PavSurf
   !    ! INTEGER,INTENT(in)::WaterSurf
 
   !    REAL(KIND(1D0)), INTENT(in) :: lvS_J_kg
   !    REAL(KIND(1D0)), INTENT(in) :: lv_J_kg
   !    REAL(KIND(1D0)), INTENT(in) :: RadMeltFact
   !    REAL(KIND(1D0)), INTENT(in) :: TempMeltFact
   !    REAL(KIND(1D0)), INTENT(in) :: SnowAlbMax
   !    REAL(KIND(1D0)), INTENT(in) :: SnowDensMax
   !    REAL(KIND(1D0)), INTENT(in) :: SnowDensMin
   !    REAL(KIND(1D0)), INTENT(in) :: Temp_C
   !    REAL(KIND(1D0)), INTENT(in) :: Precip
   !    REAL(KIND(1D0)), INTENT(in) :: PrecipLimit
   !    REAL(KIND(1D0)), INTENT(in) :: PrecipLimitAlb
   !    REAL(KIND(1D0)), INTENT(in) :: nsh_real
   !    REAL(KIND(1D0)), INTENT(in) :: albedo_snow
   !    REAL(KIND(1D0)), INTENT(in) :: tau_r
   !    ! REAL(KIND(1d0)),INTENT(in)::waterdens
 
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Tsurf_ind
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Tsurf_ind_snow
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: state_surf
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: qn1_ind_snow
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: kup_ind_snow
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowWater
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: deltaQi
 
   !    !Input and output as this is updated in this subroutine
   !    REAL(KIND(1D0)), INTENT(in) :: SnowAlb_in
   !    REAL(KIND(1D0)), INTENT(in) :: SnowfallCum_in
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowPack_in
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowFrac_in
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowDens_in
   !    REAL(KIND(1D0)), INTENT(out) :: SnowAlb_out
   !    REAL(KIND(1D0)), INTENT(out) :: SnowfallCum_out
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowPack_out
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowFrac_out
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowDens_out
 
   !    REAL(KIND(1D0)) :: SnowAlb
   !    REAL(KIND(1D0)) :: SnowfallCum
   !    REAL(KIND(1D0)), DIMENSION(nsurf) :: SnowPack
   !    REAL(KIND(1D0)), DIMENSION(nsurf) :: SnowFrac
   !    REAL(KIND(1D0)), DIMENSION(nsurf) :: SnowDens
 
   !    !Output:
   !    REAL(KIND(1D0)), INTENT(out) :: mwh
   !    REAL(KIND(1D0)) :: fwh
   !    REAL(KIND(1D0)), INTENT(out) :: Qm
   !    REAL(KIND(1D0)), INTENT(out) :: QmFreez
   !    REAL(KIND(1D0)), INTENT(out) :: QmRain
 
   !    ! REAL(KIND(1D0)), INTENT(out) :: veg_fr
 
   !    INTEGER, DIMENSION(nsurf), INTENT(out) :: snowCalcSwitch
 
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_melt
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_freezState
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_rain
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezMelt
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezState
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezStateVol
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: rainOnSnow
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowDepth
   !    REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: mw_ind
 
   !    REAL(KIND(1D0)), DIMENSION(ncolumnsDataOutSnow_notime), INTENT(out) :: dataOutLineSnow
 
   !    SnowPack = SnowPack_in
   !    SnowFrac = SnowFrac_in
   !    SnowAlb = SnowAlb_in
   !    SnowDens = SnowDens_in
   !    SnowfallCum = SnowfallCum_in
 
   !    ! IF (SnowUse == 1) THEN
   !    SnowDens = update_snow_dens( &
   !               tstep, SnowFrac, SnowDens, &
   !               tau_r, SnowDensMax, SnowDensMin)
 
   !    CALL MeltHeat( &
   !       lvS_J_kg, lv_J_kg, tstep*1D0, RadMeltFact, TempMeltFact, & !input
   !       SnowAlbMax, SnowDensMin, Temp_C, Precip, PrecipLimit, PrecipLimitAlb, &
   !       nsh_real, sfr_surf, Tsurf_ind, state_surf, qn1_ind_snow, &
   !       SnowWater, deltaQi, &
   !       SnowPack, SnowFrac, SnowAlb, SnowDens, SnowfallCum, & !inout
   !       mwh, fwh, Qm, QmFreez, QmRain, snowCalcSwitch, & !output
   !       Qm_melt, Qm_freezState, Qm_rain, FreezMelt, FreezState, FreezStateVol, &
   !       rainOnSnow, SnowDepth, mw_ind)
 
   !    ! ELSE ! no snow calculation
   !    !    mwh = 0
   !    !    fwh = 0
   !    !    Qm = 0
   !    !    QmFreez = 0
   !    !    QmRain = 0
   !    !    SnowfallCum = 0
   !    !    snowCalcSwitch = 0
   !    !    Qm_melt = 0
   !    !    Qm_freezState = 0
   !    !    Qm_rain = 0
   !    !    FreezMelt = 0
   !    !    FreezState = 0
   !    !    FreezStateVol = 0
   !    !    rainOnSnow = 0
   !    !    SnowDepth = 0
   !    !    mw_ind = 0
   !    !    SnowFrac = 0
 
   !    ! END IF
 
   !    ! ! update veg_fr
   !    ! CALL veg_fr_snow( &
   !    !    sfr_surf, SnowFrac, & !input
   !    !    veg_fr) !output
 
   !    SnowPack_out = SnowPack
   !    SnowFrac_out = SnowFrac
   !    SnowAlb_out = SnowAlb
   !    SnowfallCum_out = SnowfallCum
   !    SnowDens_out = SnowDens
 
   !    ! pack output into one line
   !    dataOutLineSnow = [ &
   !                      SnowPack_out(1:nsurf), mw_ind(1:nsurf), Qm_melt(1:nsurf), & !26
   !                      Qm_rain(1:nsurf), Qm_freezState(1:nsurf), SnowFrac_out(1:(nsurf - 1)), & !46
   !                      rainOnSnow(1:nsurf), & !53
   !                      qn1_ind_snow(1:nsurf), kup_ind_snow(1:nsurf), freezMelt(1:nsurf), & !74
   !                      SnowWater(1:nsurf), SnowDens_out(1:nsurf), & !88
   !                      snowDepth(1:nsurf), Tsurf_ind_snow(1:nsurf), &
   !                      albedo_snow]
   !    ! dataOutLineSnow=set_nan(dataOutLineSnow)
 
   ! END SUBROUTINE Snow_cal_MeltHeat
 
   SUBROUTINE meltheat( &
      lvS_J_kg, & !input
      lv_J_kg, &
      tstep_real, &
      RadMeltFact, &
      TempMeltFact, &
      SnowAlbMax, &
      SnowDensMin, &
      Temp_C, &
      Precip, &
      PrecipLimit, &
      PrecipLimitAlb, &
      nsh_real, &
      sfr_surf, &
      Tsurf_ind, &
      state_id, &
      qn1_ind_snow, &
      SnowWater, &
      deltaQi, &
      SnowPack, & !inoout
      SnowFrac, &
      SnowAlb, &
      SnowDens, &
      SnowfallCum, &
      mwh, & !output
      fwh, &
      Qm, &
      QmFreez, &
      QmRain, &
      snowCalcSwitch, &
      Qm_melt, &
      Qm_freezState, &
      Qm_rain, &
      FreezMelt, &
      FreezState, &
      FreezStateVol, &
      rainOnSnow, &
      SnowDepth, &
      mw_ind)
 
      IMPLICIT NONE
 
      !These are input to the module
      ! INTEGER, INTENT(in)::bldgsurf
      ! INTEGER, INTENT(in)::nsurf
      ! INTEGER, INTENT(in)::PavSurf
      ! INTEGER, INTENT(in)::WaterSurf
      REAL(KIND(1D0)), PARAMETER :: waterDens = 999.8395 !Density of water in 0 cel deg
 
      REAL(KIND(1D0)), INTENT(in) :: lvS_J_kg
      REAL(KIND(1D0)), INTENT(in) :: lv_J_kg
      REAL(KIND(1D0)), INTENT(in) :: tstep_real
      REAL(KIND(1D0)), INTENT(in) :: RadMeltFact
      REAL(KIND(1D0)), INTENT(in) :: TempMeltFact
      REAL(KIND(1D0)), INTENT(in) :: SnowAlbMax
      REAL(KIND(1D0)), INTENT(in) :: SnowDensMin
      REAL(KIND(1D0)), INTENT(in) :: Temp_C
      REAL(KIND(1D0)), INTENT(in) :: Precip
      REAL(KIND(1D0)), INTENT(in) :: PrecipLimit
      REAL(KIND(1D0)), INTENT(in) :: PrecipLimitAlb
      REAL(KIND(1D0)), INTENT(in) :: nsh_real
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Tsurf_ind
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: state_id
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: qn1_ind_snow
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowWater
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: deltaQi
 
      !Input and output as this is updated in this subroutine
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowPack
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowFrac
      REAL(KIND(1D0)), INTENT(inout) :: SnowAlb
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowDens
      REAL(KIND(1D0)), INTENT(inout) :: SnowfallCum
 
      !Output:
      REAL(KIND(1D0)), INTENT(out) :: mwh
      REAL(KIND(1D0)), INTENT(out) :: fwh
      REAL(KIND(1D0)), INTENT(out) :: Qm
      REAL(KIND(1D0)), INTENT(out) :: QmFreez
      REAL(KIND(1D0)), INTENT(out) :: QmRain
 
      !Output, dimension nsurf
      INTEGER, DIMENSION(nsurf), INTENT(out) :: snowCalcSwitch
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_melt
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_freezState
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: Qm_rain
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezMelt
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezState
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: FreezStateVol
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: rainOnSnow
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowDepth
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: mw_ind
 
      ! local variables:
      REAL(KIND(1D0)) :: AdjMeltFact
      REAL(KIND(1D0)) :: Watfreeze
 
      REAL(KIND(1D0)), PARAMETER :: cw = 4190 !,ci=2090   !Specific heat capacity of water
 
      INTEGER :: is, xx
 
      !Initialize snow variables
      mwh = 0 !Initialize snow melt and heat related to snowmelt
      fwh = 0
      qm = 0
      qmfreez = 0
      qmrain = 0
      snowcalcswitch = 0
      qm_melt = 0
      qm_freezstate = 0
      qm_rain = 0
      freezmelt = 0
      freezstate = 0
      freezstatevol = 0
      rainonsnow = 0
      snowdepth = 0
      mw_ind = 0
 
      !===dummy calculations===
      ! xx = bldgsurf
      ! xx = PavSurf
      !===dummy calculations end===
 
      !=========================================================================================
      DO is = 1, nsurf !Go each surface type through
         IF (sfr_surf(is) /= 0) THEN !If surface type existing,
 
            IF (snowpack(is) > 0) THEN !If SnowPack existing, calculate meltwater related water flows
 
               snowdepth(is) = (snowpack(is)/1000)*waterdens/snowdens(is) !Snow depth in m
 
               !Calculate meltwater related water flows with hourly degree-day method.
 
               !These are for snow melting
               IF (temp_c >= 0) THEN
                  IF (qn1_ind_snow(is) < 0) THEN
                     mw_ind(is) = tempmeltfact*temp_c !(mm C−1 h−1)*(C) = in mm h-1
                  ELSE
                     mw_ind(is) = radmeltfact*(qn1_ind_snow(is)) !(mm m2 W−1 h−1)*(W m-2)= mm h-1 ??
                  END IF
 
               ELSE !Freezing equations
                  adjmeltfact = 1 !Relationship between the temperature melt and freezing factors
                  mw_ind(is) = tempmeltfact*temp_c*adjmeltfact ! in mm h-1
               END IF
               !Previous equation give the hourly values, divide these with the timestep number
               mw_ind(is) = mw_ind(is)/nsh_real
 
               IF (mw_ind(is) > snowpack(is)) mw_ind(is) = snowpack(is) !Limited by the previous timestep SnowPack
 
               !-----------------------------------------------------
               ! Heat consumed to snowmelt/refreezing within Tstep.
               ! Converted from mm nsh-1 to mm nsh-1 and to m s-1
               qm_melt(is) = waterdens*((mw_ind(is)/tstep_real)/1000)*(lvs_j_kg - lv_j_kg)
 
               !If melt is negative this means freezing water in the SnowPack
               IF (mw_ind(is) < 0) THEN
 
                  freezmelt(is) = -mw_ind(is) !Save this to variable FreezMelt
                  mw_ind(is) = 0
 
                  !Freezing water cannot exceed meltwater store
                  IF (freezmelt(is) > snowwater(is)) freezmelt(is) = snowwater(is)
 
                  !Recalculate melt related energy
                  qm_melt(is) = waterdens*((-freezmelt(is)/tstep_real)/1000)*(lvs_j_kg - lv_j_kg)
               END IF
 
               !-----------------------------------------------------
               ! If air temperature is above zero, precipitation causes advective heat to the
               ! SnowPack. Eq (23) in Sun et al., 1999
               ! Calculation done at resolution of the model timestep
               IF (temp_c >= preciplimit .AND. precip > 0) THEN
                  qm_rain(is) = waterdens*cw*(temp_c - preciplimit)*(precip*0.001/tstep_real) !in W m-2
                  IF (qm_rain(is) < 0) THEN !Can only be positive
                     qm_rain(is) = 0
                  ELSE
                     rainonsnow(is) = precip !Save information on the rain on snow event
                  END IF
               END IF
 
            END IF !End if SnowPack
 
            !=================================================================
 
            !Freeze surface water state_id if cold enough.
            IF (tsurf_ind(is) < 0 .AND. state_id(is) > 0) THEN
 
               snowcalcswitch(is) = 1 !If water on ground this forms ice and snow calculations are made
 
               !Other surfaces than water treated first
               IF (is /= watersurf) THEN
 
                  !FreezState(is) = state_id(is)
                  !Previously all state_id could freeze in 5-min timestep. Now we calculate how much water
                  !can freeze in a timestep based on the same temperature freezing fraction.
                  freezstate(is) = -tempmeltfact*tsurf_ind(is)/nsh_real
 
                  !The amount of freezing water cannot be greater than the surface state_id
                  IF (freezstate(is) > state_id(is)) freezstate(is) = state_id(is)
 
                  IF (snowpack(is) == 0 .OR. snowfrac(is) == 0) THEN !SnowPack forms
                     freezstatevol(is) = freezstate(is)
                  ELSE !There is snow already on ground
                     freezstatevol(is) = freezstate(is)*(1 - snowfrac(is))/snowfrac(is)
                  END IF
 
                  ! If the amount of freezing water is very small and there is state_id left to the ground
                  ! no freezing of water will take place
                  IF (freezstatevol(is) < 0.00000000001 .AND. freezstate(is) < state_id(is)) THEN
                     freezstate(is) = 0
                     freezstatevol(is) = 0
                  END IF
 
                  !Calculate the heat exchange in W m-2
                  qm_freezstate(is) = -waterdens*(freezstate(is)/tstep_real/1000)*(lvs_j_kg - lv_j_kg)
 
                  !Water surface separately
               ELSE
                  !Calculate average value how much water can freeze above the water areas
                  !Equation is -hA(T-T0) = rhoV(Cp+dT +Lf) in 5-min timestep
                  !h=convective heat trasnfer,A, area of water,rwo water density,V volume, dT temperature difference
                  !before and end of the 5-min period. dT equals zero, h=100 and when multiplied with Area, the equation
                  !simplyfies to the this. LJ 14 July 2015
                  watfreeze = 100*(0 - temp_c)/(waterdens*(lvs_j_kg - lv_j_kg))
                  freezstate(is) = watfreeze
                  qm_freezstate(is) = -waterdens*(watfreeze/tstep_real/1000)*(lvs_j_kg - lv_j_kg)
               END IF
 
            END IF
 
            !======================================================================
            ! Define if any snowmelt calculations are made: SnowPack existing,
            ! freezing occuring on ground or from precip
            IF (is /= watersurf) THEN
               IF (snowpack(is) > 0 .OR. (precip > 0 .AND. tsurf_ind(is) < 0)) THEN
                  snowcalcswitch(is) = 1
               END IF
            ELSE !Water surface separately
               IF (snowpack(watersurf) > 0 .OR. freezstate(watersurf) > 0) THEN
                  snowcalcswitch(watersurf) = 1
               END IF
            END IF
 
            !Update snow density of each surface
            IF (precip > 0 .AND. tsurf_ind(is) < 0 .AND. snowpack(is) > 0) THEN
               snowdens(is) = snowdens(is)*snowpack(is)/(snowpack(is) + precip) + snowdensmin*precip/(snowpack(is) + precip)
            END IF
 
            !Weighted variables for the whole area
            mwh = mwh + mw_ind(is)*sfr_surf(is)*snowfrac(is) !Snowmelt
            fwh = fwh + freezmelt(is)*sfr_surf(is)*snowfrac(is) !Freezing water
            qm = qm + qm_melt(is)*sfr_surf(is)*snowfrac(is) !Energy consumed to the melt/freezing.
            qmrain = qmrain + qm_rain(is)*sfr_surf(is)*snowfrac(is) !Rain on snow
            qmfreez = qmfreez + deltaqi(is)*sfr_surf(is)*snowfrac(is) + qm_freezstate(is)*sfr_surf(is)*(1 - snowfrac(is)) !Freezing water
         END IF
 
      END DO !End surface type
 
      !Update snow albedo to its maximum value if precipitation exists
      IF (precip > 0 .AND. sum(snowpack) > 0 .AND. temp_c < 0) THEN
 
         snowfallcum = snowfallcum + precip
 
         IF (snowfallcum > preciplimitalb) THEN
 
            snowalb = snowalbmax
            snowfallcum = 0
         END IF
      ELSE
 
         snowfallcum = 0
      END IF
 
   END SUBROUTINE meltheat
 
   !===============================================================================================
   !===============================================================================================
   SUBROUTINE snowcalc( &
      tstep, imin, it, dectime, is, & !input
      snowCalcSwitch, & !input
      EvapMethod, CRWmin, CRWmax, nsh_real, lvS_J_kg, avdens, &
      avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, avcp, sIce_hPa, &
      PervFraction, vegfraction, addimpervious, &
      vpd_hPa, qn_e, s_hPa, ResistSurf, RA, rb, tlv, snowdensmin, SnowProf_24hr, precip, &
      PipeCapacity, RunoffToWater, &
      addVeg, SnowLimPaved, SnowLimBldg, FlowChange, drain, &
      WetThresh, stateOld, mw_ind, SoilStoreCap, rainonsnow, &
      freezmelt, freezstate, freezstatevol, &
      Qm_Melt, Qm_rain, Tsurf_ind, sfr_surf, dayofWeek_id, StoreDrainPrm, SnowPackLimit, &
      AddWater, addwaterrunoff, &
      soilstore_id, SnowPack, SurplusEvap, & !inout
      SnowFrac, SnowWater, iceFrac, SnowDens, &
      runoffAGimpervious, runoffAGveg, surplusWaterBody, &
      ev_tot, qe_tot, runoff_tot, surf_chang_tot, chSnow_tot, & ! output
      rss_surf, & ! output
      runoff_snowfree, chang, changSnow, SnowToSurf, state_id, ev_snow, &
      SnowRemoval, swe, &
      runoffPipes, mwstore, runoffwaterbody)
 
      !Calculation of snow and water balance on 5 min timestep. Treats snowfree and snow covered
      !areas separately. Weighting is taken into account in the overall values.
      !Last modified:
      !  LJ in 6 May 2015 - Modified to run with timestep
      !  HCW 06 Mar 2015 - Unused variable 'i' removed.
      !  HCW 26 Jan 2015 - Added weekday/weekend option for snow clearing profiles
      !  LJ in 24 May 2013
      !========================================================================
      USE waterdist_module, ONLY: updateflood
 
      IMPLICIT NONE
      ! INTEGER, PARAMETER::nsurf = 7! number of surface types
      ! INTEGER, PARAMETER::PavSurf = 1  !New surface classes: Grass = 5th/7 surfaces
      ! INTEGER, PARAMETER::BldgSurf = 2  !New surface classes: Grass = 5th/7 surfaces
      ! INTEGER, PARAMETER::ConifSurf = 3  !New surface classes: Grass = 5th/7 surfaces
      ! INTEGER,PARAMETER::DecidSurf = 4  !New surface classes: Grass = 5th/7 surfaces
      ! INTEGER,PARAMETER::GrassSurf = 5
      ! INTEGER, PARAMETER::BSoilSurf = 6!New surface classes: Grass = 5th/7 surfaces
      ! INTEGER, PARAMETER::WaterSurf = 7
 
      INTEGER, PARAMETER :: snowfractionchoice = 2 ! this PARAMETER is used all through the model
      REAL(KIND(1D0)), PARAMETER :: waterDens = 999.8395 !Density of water in 0 cel deg
 
      ! INTEGER,INTENT(in)::id
      ! INTEGER,INTENT(in)::nsurf
      INTEGER, INTENT(in) :: tstep
      INTEGER, INTENT(in) :: imin
      INTEGER, INTENT(in) :: it
      INTEGER, INTENT(in) :: is
 
      INTEGER, DIMENSION(nsurf), INTENT(in) :: snowCalcSwitch
 
      ! INTEGER,INTENT(in)::ConifSurf
      ! INTEGER,INTENT(in)::BSoilSurf
      ! INTEGER,INTENT(in)::BldgSurf
      ! INTEGER,INTENT(in)::PavSurf
      ! INTEGER,INTENT(in)::WaterSurf
      INTEGER, INTENT(in) :: EvapMethod !Evaporation calculated according to Rutter (1) or Shuttleworth (2)
      INTEGER, DIMENSION(3), INTENT(in) :: DayofWeek_id
 
      REAL(KIND(1D0)), INTENT(in) :: dectime
      REAL(KIND(1D0)), INTENT(in) :: CRWmin
      REAL(KIND(1D0)), INTENT(in) :: CRWmax
      REAL(KIND(1D0)), INTENT(in) :: nsh_real
      REAL(KIND(1D0)), INTENT(in) :: lvS_J_kg
      ! REAL(KIND(1d0)), INTENT(in)::lv_j_kg
      REAL(KIND(1D0)), INTENT(in) :: avdens
      REAL(KIND(1D0)), INTENT(in) :: vpd_hPa ! vapour pressure deficit [hPa]
      REAL(KIND(1D0)), INTENT(in) :: qn_e !net available energy for evaporation
      REAL(KIND(1D0)), INTENT(in) :: avRh
      REAL(KIND(1D0)), INTENT(in) :: Press_hPa
      REAL(KIND(1D0)), INTENT(in) :: Temp_C
      REAL(KIND(1D0)), INTENT(in) :: RAsnow
      REAL(KIND(1D0)), INTENT(in) :: psyc_hPa
      REAL(KIND(1D0)), INTENT(in) :: avcp
      REAL(KIND(1D0)), INTENT(in) :: sIce_hPa
      REAL(KIND(1D0)), INTENT(in) :: PervFraction
      REAL(KIND(1D0)), INTENT(in) :: vegfraction
      REAL(KIND(1D0)), INTENT(in) :: addimpervious
      ! REAL(KIND(1d0)),INTENT(in)::numPM
      REAL(KIND(1D0)), INTENT(in) :: s_hPa
      REAL(KIND(1D0)), INTENT(in) :: ResistSurf
      ! REAL(KIND(1d0)),INTENT(in)::sp
      REAL(KIND(1D0)), INTENT(in) :: RA
      REAL(KIND(1D0)), INTENT(in) :: rb
      REAL(KIND(1D0)), INTENT(in) :: tlv
      REAL(KIND(1D0)), INTENT(in) :: snowdensmin
      REAL(KIND(1D0)), INTENT(in) :: precip
      REAL(KIND(1D0)), INTENT(in) :: PipeCapacity
      REAL(KIND(1D0)), INTENT(in) :: RunoffToWater
      REAL(KIND(1D0)), INTENT(in) :: addVeg
      REAL(KIND(1D0)), INTENT(in) :: SnowLimPaved
      REAL(KIND(1D0)), INTENT(in) :: SnowLimBldg
      REAL(KIND(1D0)), INTENT(in) :: FlowChange
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: drain
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: WetThresh
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: stateOld
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: mw_ind
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SoilStoreCap
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: rainonsnow
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: freezmelt
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: freezstate
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: freezstatevol
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Qm_Melt
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Qm_rain
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: Tsurf_ind
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowPackLimit
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: AddWater
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: addwaterrunoff
      REAL(KIND(1D0)), DIMENSION(6, nsurf), INTENT(in) :: StoreDrainPrm
      REAL(KIND(1D0)), DIMENSION(0:23, 2), INTENT(in) :: SnowProf_24hr
 
      !Updated status: input and output
      REAL(KIND(1D0)), INTENT(inout) :: runoffAGveg
      REAL(KIND(1D0)), INTENT(inout) :: runoffAGimpervious
      REAL(KIND(1D0)), INTENT(inout) :: surplusWaterBody
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: soilstore_id
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowPack
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowFrac
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowWater
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: iceFrac
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(inout) :: SnowDens
      REAL(KIND(1D0)), DIMENSION(2), INTENT(inout) :: SurplusEvap
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: rss_surf
      REAL(KIND(1D0)), DIMENSION(nsurf) :: runoffSnow_surf !Initialize for runoff caused by snowmelting
      REAL(KIND(1D0)), DIMENSION(nsurf) :: runoff_snowfree
      ! REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: runoffSoil
      REAL(KIND(1D0)), DIMENSION(nsurf) :: chang
      REAL(KIND(1D0)), DIMENSION(nsurf) :: changSnow
      REAL(KIND(1D0)), DIMENSION(nsurf) :: SnowToSurf
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: state_id
      REAL(KIND(1D0)), DIMENSION(nsurf) :: SnowDepth
      REAL(KIND(1D0)), DIMENSION(nsurf) :: ev_snow
      REAL(KIND(1D0)), DIMENSION(2), INTENT(out) :: SnowRemoval
 
      REAL(KIND(1D0)), INTENT(out) :: swe
      REAL(KIND(1D0)) :: ev_snowfree
      REAL(KIND(1D0)), INTENT(out) :: ev_tot
      REAL(KIND(1D0)), INTENT(out) :: chSnow_tot
      REAL(KIND(1D0)), INTENT(out) :: qe_tot
      REAL(KIND(1D0)), INTENT(out) :: runoff_tot
      REAL(KIND(1D0)), INTENT(out) :: surf_chang_tot
      REAL(KIND(1D0)), INTENT(out) :: runoffPipes
      REAL(KIND(1D0)), INTENT(out) :: mwstore
      REAL(KIND(1D0)), INTENT(out) :: runoffwaterbody
 
      REAL(KIND(1D0)) :: qe
 
      ! REAL(KIND(1d0))::Evap_SUEWS_Snow
      REAL(KIND(1D0)) :: MeltExcess !Excess melt water that needs to leave SnowPack
      REAL(KIND(1D0)) :: snowTotInit
      REAL(KIND(1D0)) :: EvPart
      REAL(KIND(1D0)) :: runoffTest
      REAL(KIND(1D0)) :: snowFracFresh1 !Snow fraction for newly formed SnowPack
      REAL(KIND(1D0)) :: snowFracFresh2 !Snow fraction for newly formed SnowPack from state_id only
      REAL(KIND(1D0)) :: snowFracOld
      REAL(KIND(1D0)) :: WaterHoldCapFrac
      REAL(KIND(1D0)) :: FWC !Water holding capacity of snow in mm
      REAL(KIND(1D0)) :: tlv_sub
      ! REAL(KIND(1d0)):: SnowDepletionCurve
 
      INTEGER :: iu !1=weekday OR 2=weekend
      REAL(KIND(1D0)), PARAMETER :: IPThreshold_mmhr = 10 !Threshold for intense precipitation [mm hr-1]
 
      REAL(KIND(1D0)), DIMENSION(7) :: capStore ! current storage capacity [mm]
 
      !========================================================================
      !Initialize variables for the calculation of water storages and evaporation
      ev_tot = 0
      qe_tot = 0
      runoff_tot = 0
      surf_chang_tot = 0
 
      !swe = 0
      chsnow_tot = 0
      !runoffPipes = 0
      !mwstore = 0
      !runoffwaterbody = 0
      rss_surf = 0
      state_id = stateold
      snowdepth = 0
      ev_snow = 0
      !SnowRemoval = 0
      tlv_sub = 0
 
      ! Use weekday or weekend snow clearing profile
      iu = 1 !Set to 1=weekday
      IF (dayofweek_id(1) == 1 .OR. dayofweek_id(1) == 7) iu = 2 !Set to 2=weekend
 
      !write(*,*) is
      runoffsnow_surf(is) = 0 !Initialize for runoff caused by snowmelting
      runoff_snowfree(is) = 0
      ! runoffSoil(is) = 0
      chang(is) = 0
      changsnow(is) = 0
      runofftest = 0
      snowtosurf(is) = 0
      evpart = 0
      ev_snowfree = 0
      snowfracfresh1 = 0
      snowfracfresh2 = 0
      snowfracold = 0
 
      !Initial SnowPack + meltwater in it
      snowtotinit = snowpack(is) + snowwater(is)
 
      !Calculate water holding capacity (Jin et al. 1999)
      IF (snowdens(is) >= 200) THEN
         waterholdcapfrac = crwmin
      ELSE
         waterholdcapfrac = crwmin + (crwmax - crwmin)*(200 - snowdens(is))/200
      END IF
 
      !======================================================================
      ! Calculate evaporation from SnowPack and snow free surfaces (in mm)
      ! IF (SnowFrac(is)<1) CALL Evap_SUEWS !ev and qe for snow free surface out
      capstore(is) = storedrainprm(6, is)
      IF (snowfrac(is) < 1) CALL cal_evap( &
         evapmethod, state_id(is), wetthresh(is), capstore(is), & !input
         vpd_hpa, avdens, avcp, qn_e, s_hpa, psyc_hpa, resistsurf, ra, rb, tlv, &
         rss_surf(is), ev_snowfree, qe) !output
 
      IF (snowfrac(is) > 0 .AND. snowcalcswitch(is) > 0) THEN
         CALL evap_suews_snow(qm_melt(is), qm_rain(is), lvs_j_kg, avdens, avrh, press_hpa, temp_c, rasnow, &
                              psyc_hpa, tstep, avcp, sice_hpa, dectime, ev_snow(is), tlv_sub)
      END IF
 
      !If not enough water for evaporation in impervious surfaces,
      !evaporation is taken from pervious surfaces
      IF (is > 2) THEN
         IF (pervfraction /= 0) THEN
            evpart = (surplusevap(pavsurf)*sfr_surf(pavsurf) + surplusevap(bldgsurf)*sfr_surf(bldgsurf))/pervfraction
         END IF
      END IF
 
      !============================================================================
      !Water surface is treated separately
      IF (is == watersurf .AND. sfr_surf(watersurf) > 0) GO TO 606
 
      !The calculations are divided into 2 main parts
      ! 1) Surface is fully covered with snow at the beginning of the time step
      ! 2) Surface is not fully covered with snow but rather part is snow free OR
      !    surface not orginally covered with snow, but the snow forms at the current timestep
 
      !1)------------------------------------------------------------------
      !  ------------------------------------------------------------------
      IF (snowcalcswitch(is) > 0) THEN
         IF (snowpack(is) > 0 .AND. snowfrac(is) == 1) THEN
 
            ev_snow(is) = ev_snow(is) + evpart !Evaporation surplus
 
            !(Snowfall per interval+freezing of melt water and surface state_id) - (meltwater+evaporation from SnowPack)
            changsnow(is) = (precip + freezmelt(is)) - (mw_ind(is) + ev_snow(is)) !Calculate change in SnowPack (in mm)
 
            !If rain on snow event, add this water to SnowWater
            IF (rainonsnow(is) > 0) THEN
               changsnow(is) = changsnow(is) - precip
               snowwater(is) = snowwater(is) + rainonsnow(is)
            END IF
 
            snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
 
            !---------If SnowPack exists after the state_id calculations
            IF (snowpack(is) > 0) THEN
 
               !Add melted water to meltstore and freeze water according to freezMelt(is)
               snowwater(is) = snowwater(is) + mw_ind(is) - freezmelt(is)
 
               !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
               fwc = waterholdcapfrac*snowpack(is)
 
               !If FWC is exceeded, excess meltwater (MeltExcess) will leave from the SnowPack
               IF (snowwater(is) >= fwc) THEN
                  meltexcess = 0 !Initialize the excess meltwater
                  meltexcess = snowwater(is) - fwc !Calculate the exceess water
                  snowwater(is) = fwc !Update the SnowWater to the maximum it can hold
                  runoffsnow_surf(is) = runoffsnow_surf(is) + meltexcess
               END IF
 
               !At the end of the hour calculate possible snow removal
               IF (snowprof_24hr(it, iu) == 1 .AND. is < 3 .AND. (imin == (nsh_real - 1)/nsh_real*60)) &
                  CALL snow_removal( &
                  is, &
                  snowfrac, sfr_surf, &
                  snowpack, snowremoval, &
                  snowlimpaved, snowlimbldg)
               !----------If SnowPack is negative, it melts at this timestep
            ELSEIF (snowpack(is) < 0) THEN
 
               !If freezing meltwater inside this timestep, remove it from the SnowWater
               snowwater(is) = snowwater(is) - freezmelt(is) + mw_ind(is) + snowpack(is)
               snowpack(is) = 0.0 !Set the snow pack and snow
               snowfracold = 1
               snowfrac(is) = 0
               snowdens(is) = 0
 
               IF (snowwater(is) < 0) THEN !Not enough water in the meltwater store,
                  ev_snow(is) = ev_snow(is) + snowwater(is) !QUESTION: evaporation from snow is decreased?
                  IF (ev_snow(is) < 0) ev_snow(is) = 0
                  changsnow(is) = changsnow(is) + snowwater(is)
                  snowwater(is) = 0
               ELSE
                  chang(is) = snowwater(is) !Meltwater goes to surface state_id as no snow exists anymore
                  state_id(is) = state_id(is) + chang(is)
                  snowwater(is) = 0
               END IF
            END IF !SnowPack negative or positive
 
            !2)------Surface not fully covered with snow-------------------------------------------
            !  ------------------------------------------------------------------------------------
         ELSEIF (snowfrac(is) < 1) THEN
 
            !Snow calculations: SnowPack can either exist or form at the current timestep
            IF (snowpack(is) > 0) THEN
               ev_snow(is) = ev_snow(is) + evpart !Evaporation surplus
 
               !----SnowPack water balance for the whole surface area. In reality snow depth = SnowPack/SnowFrac(is)
               !(Snowfall per interval+freezing of melt water and surface state_id) - (meltwater+evaporation from SnowPack)
               changsnow(is) = (precip + freezmelt(is) + freezstatevol(is)) - (mw_ind(is) + ev_snow(is)) !Calculate change in SnowPack (in mm)
 
               !If rain on snow event, add this water to SnowWater
               IF (rainonsnow(is) > 0) THEN
                  changsnow(is) = changsnow(is) - precip
                  snowwater(is) = snowwater(is) + rainonsnow(is)
               END IF
               snowpack(is) = snowpack(is) + changsnow(is)
 
               !The fraction of snow will update when:
               !a) Surface state_id is dry but precipitation occurs =1
               !b) There is both precipitation and all surface state_id freezes =1
               !c) No precipitation but all state_id freezes at a single timestep =2
               !d) Part of the surface freezes
               IF (precip > 0 .AND. freezstate(is) == state_id(is)) THEN !both a) and b)
                  snowfracfresh1 = 1
               ELSEIF (precip == 0 .AND. freezstate(is) > 0 .AND. freezstate(is) == state_id(is)) THEN
                  snowfracfresh1 = 1
 
                  !snowFracFresh1=SnowDepletionCurve(is,SnowPack(is),SnowPackLimit(is))
                  !if (snowFracFresh1<0.001) snowFracFresh1=0.001
               ELSEIF (freezstate(is) > 0 .AND. freezstate(is) < state_id(is)) THEN !This if not all water freezes
                  snowfracfresh1 = 0.95 !Now this fraction set to something close to one. Should be improved in the future at some point
                  !if (is==1)then
                  ! write(*,*) id,it,imin,SnowFrac(is),FreezState(is),state_id(is)
                  ! pause
                  !endif
               END IF
 
               !SnowPack can also form at the current timestep (2). If this forms purely from snowfall or/and all water at surface freezes,
               !the whole surface will be covered with snow. If there is water on ground this snowfall can immediately melt
               !and in this case the snow fraction is not necessarily 1 but its information is saved to snowFracFresh that
               !is taken into account in snow fraction after calculation of state_id.
            ELSEIF (snowpack(is) == 0 .AND. tsurf_ind(is) < 0) THEN
 
               !The fraction of snow will get a value of 1 (ie full snow cover):
               !Surface state_id is dry but precipitation occurs, no precipitation but all state_id freezes at a single timestep,
               !There is both precipitation and all surface state_id freezes
               IF ((precip > 0 .AND. state_id(is) == 0) .OR. (precip == 0 .AND. freezstate(is) == state_id(is)) .OR. &
                   (precip > 0 .AND. freezstate(is) == state_id(is))) THEN
 
                  !ev=ev+EvPart
                  changsnow(is) = precip + freezstatevol(is)
                  snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
 
                  snowfracfresh1 = 1
                  icefrac(is) = freezstate(is)/(freezstate(is) + precip)
                  snowdens(is) = snowdensmin
               END IF
 
               IF (freezstate(is) > 0 .AND. freezstate(is) < state_id(is)) THEN
 
                  changsnow(is) = precip + freezstatevol(is)
                  snowpack(is) = snowpack(is) + changsnow(is) !Update SnowPack
                  snowfracfresh2 = 0.95 !Now this fraction set to something close to one. Should be improved in the future at some point
 
                  !snowFracFresh2=SnowDepletionCurve(is,SnowPack(is),SnowPackLimit(is))
                  !if (snowFracFresh2<0.001) snowFracFresh2=0.001
                  icefrac(is) = 1
                  snowdens(is) = snowdensmin
                  !write(*,*) 2,is,id,it,imin,SnowFrac(is),FreezState(is),state_id(is),state_id(is)+Precip
                  !pause
 
               END IF
            END IF
 
            !---------If SnowPack exists after the state_id calculations
            IF (snowpack(is) > 0) THEN
 
               !Add melted water to meltstore and freeze water according to freezMelt(is)
               snowwater(is) = snowwater(is) + mw_ind(is) - freezmelt(is)
 
               !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
               fwc = waterholdcapfrac*snowpack(is)
 
               !If FWC is exceeded, excess meltwater (MeltExcess) will leave from the SnowPack
               IF (snowwater(is) >= fwc) THEN
                  meltexcess = 0 !Initialize the excess meltwater
                  meltexcess = snowwater(is) - fwc !Calculate the exceess water
                  snowwater(is) = fwc !Update the SnowWater to the maximum it can hold
 
                  !If the fraction of snow is greater than 0.8 or if the surface is is buildings,
                  !the excess water will directly go to runoff. Otherwise it will flow to the
                  !snow free area via SnowToSurf(is)
                  IF ((snowfrac(is) > 0.9 .AND. is /= bldgsurf) .OR. (is == bldgsurf)) THEN
                     runoffsnow_surf(is) = runoffsnow_surf(is) + meltexcess
                  ELSE
                     snowtosurf(is) = snowtosurf(is) + meltexcess*snowfrac(is)/(1 - snowfrac(is))
                  END IF
               END IF
 
               !At the end of the hour calculate possible snow removal
               IF (snowprof_24hr(it, iu) == 1 .AND. is < 3 .AND. (imin == (nsh_real - 1)/nsh_real*60)) &
                  CALL snow_removal( &
                  is, &
                  snowfrac, sfr_surf, &
                  snowpack, snowremoval, &
                  snowlimpaved, snowlimbldg)
 
               !----------If SnowPack is negative, it melts at this timestep
            ELSEIF (snowpack(is) < 0) THEN
 
               !If freezing meltwater inside this timestep, remove it from the SnowWater
               snowwater(is) = snowwater(is) - freezmelt(is) + mw_ind(is) + snowpack(is)
 
               snowpack(is) = 0.0 !Set the snow pack and snow
               snowfracfresh1 = 0
               snowfracfresh2 = 0
               snowdens(is) = 0
 
               IF (snowwater(is) < 0) THEN !Not enough water in the meltwater store,
                  ev_snow(is) = ev_snow(is) + snowwater(is) !QUESTION: evaporation from snow is decreased.?
                  IF (ev_snow(is) < 0) ev_snow(is) = 0
                  changsnow(is) = changsnow(is) + snowwater(is)
                  snowwater(is) = 0
               ELSE
                  snowtosurf(is) = snowtosurf(is) + snowwater(is)*snowfrac(is)/(1 - snowfrac(is))
                  snowwater(is) = 0
               END IF
            END IF !SnowPack negative or positive
 
            !--------
            !Next the snow free surface (3). Calculations only done if snowfraction is smaller than 1
            IF ((is == pavsurf .OR. is == bldgsurf) .AND. snowfrac(is) < 1) THEN !Impervious surfaces (paved, buildings)
 
               !Surface store update. If precipitation is greater than the threshold, the exceeding water
               !goes directly to runoff
               IF (precip > ipthreshold_mmhr/nsh_real) THEN
                  !runoff = runoff + (precipitation+water from the snow surface+water from other surfaces-the thereshold limit)
                  runoff_snowfree(is) = runoff_snowfree(is) + (precip + snowtosurf(is) + addwater(is) - ipthreshold_mmhr/nsh_real)
                  chang(is) = ipthreshold_mmhr/nsh_real - (drain(is) + ev_snowfree + freezstate(is))
               ELSE
                  !Add precip and water from other surfaces and remove drainage, evap and freezing of state_id
                  chang(is) = precip + snowtosurf(is) + addwater(is) - (drain(is) + ev_snowfree + freezstate(is))
               END IF
 
               state_id(is) = state_id(is) + chang(is) !Change in state_id (for whole surface area areasfr(is))
 
               !Add water from impervious grids
               ! Check sfr_surf/=0 added HCW 08 Dec 2015
               IF (is == pavsurf .AND. sfr_surf(pavsurf) > 0) state_id(is) = state_id(is) + (addimpervious)/sfr_surf(pavsurf)
 
               runoff_snowfree(is) = runoff_snowfree(is) + drain(is)*addwaterrunoff(is) !Drainage (not flowing to other surfaces) goes to runoff
 
               IF (state_id(is) < 0.0) THEN !Surface state_id cannot be negative
                  surplusevap(is) = abs(state_id(is)) !take evaporation from other surfaces in mm
                  ev_snowfree = ev_snowfree - surplusevap(is)
                  state_id(is) = 0.0
               END IF
 
            ELSEIF (is >= 3 .AND. snowfrac(is) < 1) THEN ! Pervious surfaces (conif, decid, grass unirr, grass irr)
 
               ev_snowfree = ev_snowfree + evpart
 
               !Change in water stores
               IF (vegfraction > 0) THEN
                  IF (precip + addveg*(sfr_surf(is)/vegfraction) > (ipthreshold_mmhr/nsh_real)) THEN !if 5min precipitation is larger than 10 mm
                     runoff_snowfree(is) = runoff_snowfree(is) + (precip + addveg*(sfr_surf(is)/vegfraction) + &
                                                                  snowtosurf(is) + addwater(is) - (ipthreshold_mmhr/nsh_real))
                     chang(is) = (ipthreshold_mmhr/nsh_real) - (drain(is) + ev_snowfree + freezstate(is))
                  ELSE
                     chang(is) = precip + addveg*(sfr_surf(is)/vegfraction) + snowtosurf(is) + &
                                 addwater(is) - (drain(is) + ev_snowfree + freezstate(is))
                  END IF
               ELSE
                  chang(is) = precip + snowtosurf(is) + addwater(is) - (drain(is) + ev_snowfree + freezstate(is))
               END IF
 
               state_id(is) = state_id(is) + chang(is)
 
               !Add water in soil store only if ground is not frozen
               IF (temp_c > 0) THEN
                  soilstore_id(is) = soilstore_id(is) + drain(is)*addwaterrunoff(is)*(1 - snowfrac(is))
               ELSE
                  runoff_snowfree(is) = runoff_snowfree(is) + drain(is)*addwaterrunoff(is)
               END IF
 
               !If state_id of the surface is negative, remove water from soilstore
               IF (state_id(is) < 0.0) THEN
 
                  IF ((soilstore_id(is) + state_id(is)) >= 0 .AND. temp_c > 0) THEN !If water in soilstore, water is removed
 
                     soilstore_id(is) = soilstore_id(is) + state_id(is)*(1 - snowfrac(is))
                     state_id(is) = 0.0
 
                  ELSE !If not water in the soilstore evaporation does not occur
                     chang(is) = chang(is) + state_id(is)
                     ev_snowfree = ev_snowfree + state_id(is)
                     state_id(is) = 0.0
                  END IF
               END IF !state_id is negative
 
               !If soilstorage is full at this point, excess will go to surface runoff
               IF (soilstore_id(is) > soilstorecap(is)) THEN
                  runofftest = runofftest + (soilstore_id(is) - soilstorecap(is))
                  soilstore_id(is) = soilstorecap(is)
               ELSEIF (soilstore_id(is) < 0) THEN
                  soilstore_id(is) = 0
               END IF
 
            END IF !Surface type
 
         END IF !Surface fraction
         !-------------------------------------------------------------------------------------------------------------------
 
         !Calculate change in SnowPack and state_id for the respective surface areas
         !Here the case where not all surface state_id freezes is handled
         IF (snowfracfresh2 > 0) THEN
            surf_chang_tot = (state_id(is) - stateold(is))*sfr_surf(is)*(1 - snowfrac(is)) &
                             - precip*sfr_surf(is)*(1 - snowfracfresh2)
            chsnow_tot = ((snowpack(is) + snowwater(is)) - snowtotinit)*sfr_surf(is)*(1 - snowfrac(is)) &
                         - precip*sfr_surf(is)*snowfracfresh2
         ELSE
            surf_chang_tot = (state_id(is) - stateold(is))*sfr_surf(is)*(1 - snowfrac(is))
            chsnow_tot = ((snowpack(is) + snowwater(is)) - snowtotinit)*sfr_surf(is)*max(snowfrac(is), snowfracold)
         END IF
 
         !===Update snow depth, weighted SWE, and Mwstore
         IF (snowdens(is) /= 0) THEN
            snowdepth(is) = snowpack(is)*waterdens/snowdens(is)
         END IF
 
         ! Calculate overall snow water equivalent
         swe = swe + snowpack(is)*sfr_surf(is)*max(snowfrac(is), snowfracold)
         mwstore = mwstore + snowwater(is)*sfr_surf(is)*max(snowfrac(is), snowfracold)
 
         !if (id==6.and.it==13.and.imin==20) then!
         !if (id==85.and.it==3.and.imin==10) then!
         ! if (id==92.and.it==21.and.imin==35) then!
         !  write(*,*)  ((SnowPack(is)+SnowWater(is))-snowTotInit)*sfr_surf(is)*(1-SnowFrac(is)),&
         !              runoff(is)*sfr_surf(is)*(1-SnowFrac(is)),&
         !              ev*sfr_surf(is)*(1-SnowFrac(is)),&
         !              (state_id(is)-stateOld(is))*sfr_surf(is)*(1-SnowFrac(is)),Precip*sfr_surf(is)
         !  write(*,*)  changSnow(is),runoff(is),ev,chang(is),runoffTest,FreezState(is) !changSnow(is)-freezMelt(is)
         !  write(*,*)  is,Precip,runoff_per_tstep,ev_per_tstep,surf_chang_per_tstep,chSnow_per_interval
         !  write(*,*)  is,Precip-runoff_per_tstep-ev_per_tstep,surf_chang_per_tstep+chSnow_per_interval
         !  write(*,*)  is,SnowFrac(is),sfr_surf(is),sfr_surf(is)*ev_snow(is)
         !  pause
         ! endif
 
         !Only now update the new snow fractions both in the case that snow existing already on ground
         !and snow forms at the current timestep
         IF (snowfracfresh1 > 0) snowfrac(is) = snowfracfresh1
         IF (snowfracfresh2 > 0) snowfrac(is) = snowfracfresh2
 
         !Calculate new snow fraction here.
         !Tässä ongelmana että snow fraction muuttuu vain kun on sulamisvettä ja on vika tunti.
         !Tämä ei juuri koskaan toteudu johtuen lämpötilan vuorokausisyklistä
         !Kokeile tässä ajaa kahdella tavalla 1) ei tarvita Mw:tä
         !                                    2) päivitys voi tapahtua millon vain
         !if (SnowFractionChoice==2.and.imin==(nsh_real-1)/nsh_real*60) then
         IF (snowfractionchoice == 2) THEN
            IF (snowpack(is) > 0 .AND. mw_ind(is) > 0) THEN
               snowfrac(is) = snowdepletioncurve(is, snowpack(is), snowpacklimit(is))
               IF (snowfrac(is) < 0.001) snowfrac(is) = 0.001 !The snow fraction minimum is 1% of the surface
            ELSEIF (snowpack(is) == 0) THEN
               snowfrac(is) = 0
            END IF
         END IF
      END IF !end snowCalcSwitch
 
      !Add evaporation to total
      IF (is == bldgsurf .OR. is == pavsurf) THEN
         ev_tot = ev_snowfree*sfr_surf(is)*(1 - snowfrac(is)) + ev_snow(is)*sfr_surf(is)*max(snowfrac(is), snowfracold)
         qe_tot = ev_snow(is)*tlv_sub*sfr_surf(is)*snowfrac(is) + ev_snowfree*tlv*sfr_surf(is)*(1 - snowfrac(is))
      ELSE
         ev_tot = ev_snowfree*sfr_surf(is)*(1 - snowfrac(is)) + ev_snow(is)*sfr_surf(is)*max(snowfrac(is), snowfracold)
         qe_tot = ev_snow(is)*tlv_sub*sfr_surf(is)*max(snowfrac(is), snowfracold) + ev_snowfree*tlv*sfr_surf(is)*(1 - snowfrac(is))
      END IF
 
      !========RUNOFF=======================
 
      !Add runoff to pipes
      runoffpipes = runoffpipes &
                    + runoffsnow_surf(is)*sfr_surf(is)*max(snowfrac(is), snowfracold) &
                    + runoff_snowfree(is)*sfr_surf(is)*(1 - snowfrac(is)) &
                    + runofftest*sfr_surf(is)
      CALL updateflood( &
         is, runoff_snowfree, & ! input:
         sfr_surf, pipecapacity, runofftowater, &
         runoffagimpervious, surpluswaterbody, runoffagveg, runoffpipes) ! inout:
 
      runoff_tot = runoffsnow_surf(is)*sfr_surf(is)*max(snowfrac(is), snowfracold) &
                   + runoff_snowfree(is)*sfr_surf(is)*(1 - snowfrac(is)) &
                   + runofftest*sfr_surf(is)
 
      ! !===Update snow depth, weighted SWE, and Mwstore
      ! IF (SnowDens(is) /= 0) THEN
      !    SnowDepth(is) = SnowPack(is)*waterDens/SnowDens(is)
      ! END IF
 
      ! ! Calculate overall snow water equivalent
      ! swe = swe + SnowPack(is)*sfr_surf(is)*MAX(SnowFrac(is), snowfracOld)
      ! MwStore = MwStore + SnowWater(is)*sfr_surf(is)*MAX(SnowFrac(is), snowfracOld)
 
      ! !if (id==6.and.it==13.and.imin==20) then!
      ! !if (id==85.and.it==3.and.imin==10) then!
      ! ! if (id==92.and.it==21.and.imin==35) then!
      ! !  write(*,*)  ((SnowPack(is)+SnowWater(is))-snowTotInit)*sfr_surf(is)*(1-SnowFrac(is)),&
      ! !              runoff(is)*sfr_surf(is)*(1-SnowFrac(is)),&
      ! !              ev*sfr_surf(is)*(1-SnowFrac(is)),&
      ! !              (state_id(is)-stateOld(is))*sfr_surf(is)*(1-SnowFrac(is)),Precip*sfr_surf(is)
      ! !  write(*,*)  changSnow(is),runoff(is),ev,chang(is),runoffTest,FreezState(is) !changSnow(is)-freezMelt(is)
      ! !  write(*,*)  is,Precip,runoff_per_tstep,ev_per_tstep,surf_chang_per_tstep,chSnow_per_interval
      ! !  write(*,*)  is,Precip-runoff_per_tstep-ev_per_tstep,surf_chang_per_tstep+chSnow_per_interval
      ! !  write(*,*)  is,SnowFrac(is),sfr_surf(is),sfr_surf(is)*ev_snow(is)
      ! !  pause
      ! ! endif
 
      ! !Only now update the new snow fractions both in the case that snow existing already on ground
      ! !and snow forms at the current timestep
      ! IF (snowFracFresh1 > 0) SnowFrac(is) = snowFracFresh1
      ! IF (snowFracFresh2 > 0) SnowFrac(is) = snowFracFresh2
 
      ! !Calculate new snow fraction here.
      ! !Tässä ongelmana että snow fraction muuttuu vain kun on sulamisvettä ja on vika tunti.
      ! !Tämä ei juuri koskaan toteudu johtuen lämpötilan vuorokausisyklistä
      ! !Kokeile tässä ajaa kahdella tavalla 1) ei tarvita Mw:tä
      ! !                                    2) päivitys voi tapahtua millon vain
      ! !if (SnowFractionChoice==2.and.imin==(nsh_real-1)/nsh_real*60) then
      ! IF (SnowFractionChoice == 2) THEN
      !    IF (SnowPack(is) > 0 .AND. mw_ind(is) > 0) THEN
      !       SnowFrac(is) = SnowDepletionCurve(is, SnowPack(is), SnowPackLimit(is))
      !       IF (SnowFrac(is) < 0.001) SnowFrac(is) = 0.001 !The snow fraction minimum is 1% of the surface
      !    ELSEIF (SnowPack(is) == 0) THEN
      !       SnowFrac(is) = 0
      !    END IF
      ! END IF
 
      RETURN
 
      !==========================================================================
      !WATERBODY is treated separately as state_id always below ice if ice existing
      !Calculate change in SnowPack
606   changsnow(watersurf) = (precip + freezmelt(watersurf) + freezstate(watersurf)) - &
                             (mw_ind(watersurf) + ev_snow(watersurf))
 
      snowpack(watersurf) = snowpack(watersurf) + changsnow(watersurf) !Update SnowPack
      state_id(watersurf) = state_id(watersurf) + flowchange - freezstate(watersurf) !Update state_id below ice
 
      !If SnowPack exists
      IF (snowpack(watersurf) > 0) THEN
 
         !Add melted water to meltstore and freeze water according to freezMelt(is)
         snowwater(watersurf) = snowwater(watersurf) + mw_ind(watersurf) - freezmelt(watersurf)
 
         !Calculate water holding capacity (FWC: Valeo and Ho, 2004) of the SnowPack
         fwc = waterholdcapfrac*snowpack(watersurf)
 
         !If FWC is exceeded, add meltwater to state_id
         IF (snowwater(watersurf) >= fwc .AND. temp_c >= 0) THEN
            state_id(watersurf) = state_id(watersurf) + (snowwater(watersurf) - fwc)
            snowwater(watersurf) = fwc
         END IF
 
         !If SnowPack is negative, it melts at this timestep
      ELSEIF (snowpack(is) < 0) THEN
 
         !Add water to the meltwater store
         !If freezing meltwater inside this hour, remove it from the SnowWater
         snowwater(watersurf) = snowwater(watersurf) - freezmelt(watersurf) &
                                + mw_ind(watersurf)
 
         state_id(watersurf) = state_id(watersurf) + snowwater(watersurf) + snowpack(watersurf) !Add meltwater to state_id
         snowpack(watersurf) = 0
         IF (state_id(watersurf) < 0) ev_snow(watersurf) = ev_snow(watersurf) + state_id(watersurf)
 
      END IF !SnowPack negative or positive
 
      !Check water state_id separately
      IF (state_id(watersurf) > storedrainprm(5, watersurf)) THEN
         runoff_snowfree(watersurf) = runoff_snowfree(watersurf) + (state_id(watersurf) - storedrainprm(5, watersurf))
         state_id(watersurf) = storedrainprm(5, watersurf)
         runoffwaterbody = runoffwaterbody + runoff_snowfree(watersurf)*sfr_surf(watersurf)
      ELSE
         state_id(watersurf) = state_id(watersurf) + surpluswaterbody
 
         IF (state_id(watersurf) > storedrainprm(5, watersurf)) THEN
            runoffwaterbody = runoffwaterbody + (state_id(watersurf) - storedrainprm(5, watersurf))*sfr_surf(watersurf)
            state_id(watersurf) = storedrainprm(5, watersurf)
         END IF
      END IF
 
      !Change state_id of snow and surface
      chsnow_tot = ((snowpack(watersurf) + snowwater(watersurf)) - snowtotinit)*sfr_surf(watersurf)
      !ch_per_interval=ch_per_interval+(state_id(WaterSurf)-stateOld(WaterSurf))*sfr_surf(WaterSurf)
      surf_chang_tot = (state_id(watersurf) - stateold(watersurf))*sfr_surf(watersurf)
 
      !Evaporation
      ev_tot = ev_snowfree*sfr_surf(watersurf) + ev_snow(watersurf)*sfr_surf(watersurf)
      qe_tot = ev_snow(watersurf)*tlv_sub*sfr_surf(watersurf) + ev_snowfree*tlv*sfr_surf(watersurf)
      runoff_tot = runoff_snowfree(is) !The total runoff from the area
 
      IF (snowpack(watersurf) > 0) THEN !Fraction only 1 or 0
         snowfrac(watersurf) = 1
      ELSE
         snowfrac(watersurf) = 0
      END IF
 
   END SUBROUTINE snowcalc
 
   !==========================================================================
   !==========================================================================
   !Calculates evaporation from snow surface (ev_snow).
   !Last update: LJ/Jan 2019 Function changed to subroutine. tlv_sub added to output
   SUBROUTINE evap_suews_snow(Qm, QP, lvS_J_kg, avdens, avRh, Press_hPa, Temp_C, RAsnow, psyc_hPa, &
                              tstep, avcp, sIce_hPa, dectime, ev_snow, tlv_sub)
 
      USE meteo, ONLY: sat_vap_pressice
 
      IMPLICIT NONE
 
      !INPUT
      REAL(KIND(1D0)), INTENT(in) :: Qm !melt heat,
      REAL(KIND(1D0)), INTENT(in) :: QP !advect. heat
      REAL(KIND(1D0)), INTENT(in) :: lvS_J_kg !latent heat of sublimation
      REAL(KIND(1D0)), INTENT(in) :: avdens !air density
      REAL(KIND(1D0)), INTENT(in) :: avRh !relative humidity
      REAL(KIND(1D0)), INTENT(in) :: Press_hPa !air pressure
      REAL(KIND(1D0)), INTENT(in) :: Temp_C !air temperature
      REAL(KIND(1D0)), INTENT(in) :: RAsnow !aerodyn res snow
      REAL(KIND(1D0)), INTENT(in) :: psyc_hPa !psychometric constant
      REAL(KIND(1D0)), INTENT(in) :: avcp !spec. heat,
      REAL(KIND(1D0)), INTENT(in) :: sIce_hPa !satured curve on snow
      REAL(KIND(1D0)), INTENT(in) :: dectime
      INTEGER, INTENT(in) :: tstep
 
      REAL(KIND(1D0)), INTENT(out) :: ev_snow !Evaporation
      REAL(KIND(1D0)), INTENT(out) :: tlv_sub !Latent heat for sublimation
 
      !OTHER VARIABLES
      REAL(KIND(1D0)) :: e_snow, & !PM equation obe line
                         sae_snow, & !s * (Available energy)
                         qe_snow, & !Latent heat flux
                         vdrcIce, & !Vapour pressure deficit
                         esIce_hPa, & !Saturation vapor pressure over ice
                         EaIce_hPa, & !Vapour pressure
                         tstep_real !timestep as real
 
      INTEGER :: from = 1
      !-----------------------------------------------------
 
      tstep_real = real(tstep, kind(1d0))
 
      sae_snow = sice_hpa*(qp - qm) !Calculate the driving parameter in calculation of evaporation. Järvi et al. (2015)
 
      esice_hpa = sat_vap_pressice(temp_c, press_hpa, from, dectime) !Saturation vapor pressure over ice
      eaice_hpa = avrh/100*esice_hpa !Vapour pressure of water
      vdrcice = (esice_hpa - eaice_hpa)*avdens*avcp !Vapour pressure deficit
      tlv_sub = lvs_j_kg/tstep_real !Latent heat for sublimation
      e_snow = sae_snow + vdrcice/rasnow !PM equation
      qe_snow = e_snow/(sice_hpa + psyc_hpa) !Latent heat (W/m^2)
      ev_snow = qe_snow/tlv_sub !Evaporation (in mm)
 
      RETURN
 
   END SUBROUTINE evap_suews_snow
 
   !==========================================================================
   !==========================================================================
   ! Calculates mechanical removal of snow from roofs ans roads
   SUBROUTINE snow_removal( &
      is, &
      SnowFrac, sfr_surf, &
      SnowPack, SnowRemoval, &
      SnowLimPaved, SnowLimBldg)
 
      IMPLICIT NONE
      INTEGER, INTENT(in) :: is
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowFrac, sfr_surf
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowPack
      REAL(KIND(1D0)), DIMENSION(2), INTENT(out) :: SnowRemoval
      REAL(KIND(1D0)), INTENT(in) :: SnowLimPaved, SnowLimBldg
      !write(*,*) is, SnowPack(is),SnowLimPaved,SnowLimBldg
 
      IF (is == pavsurf) THEN
         IF (snowpack(pavsurf) > snowlimpaved) THEN
            snowremoval(pavsurf) = (snowpack(pavsurf) - snowlimpaved)*sfr_surf(pavsurf)*snowfrac(pavsurf)
            snowpack(pavsurf) = snowlimpaved
            !SnowPack(PavSurf)=SnowPack(PavSurf)/SnowFrac(PavSurf)
         END IF
      END IF
      IF (is == bldgsurf) THEN
         IF (snowpack(bldgsurf) > snowlimbldg) THEN
            snowremoval(2) = (snowpack(bldgsurf) - snowlimbldg)*sfr_surf(bldgsurf)*snowfrac(bldgsurf)
            snowpack(bldgsurf) = snowlimbldg
            !SnowPack(BldgSurf)=SnowPack(BldgSurf)/SnowFrac(BldgSurf)
         END IF
      END IF
      !write(*,*) is, SnowPack(is),SnowLimPaved,SnowLimBldg
      !pause
   END SUBROUTINE snow_removal
 
   !----------------------------------------------------------------------------
   !----------------------------------------------------------------------------
   FUNCTION snowdepletioncurve(is, swe, sweD) RESULT(asc)
      !This function calculates surface coverage of snow according to the
      !depletion curves in Valeo and Ho (2004).
      !INPUT: is   Surface type number
      !       swe  Snow water content
      !       sweD Limit for
 
      USE allocatearray
 
      IMPLICIT NONE
 
      INTEGER :: is
      REAL(kind(1d0)) :: asc, swed, swe
 
      ! initialisation
      asc = 1
      !Impervious surface
      IF (is == pavsurf) THEN
 
         IF (swe <= swed) THEN !Snow water equivalent below threshold
            asc = ((swe/swed))**2
         ELSE
            asc = 1
         END IF
 
         !Bldgs surface
      ELSEIF (is == bldgsurf) THEN
 
         IF (swe <= swed) THEN
            IF ((swe/swed) < 0.9) THEN
               asc = (swe/swed)*0.5
            ELSE
               asc = (swe/swed)**8
            END IF
         ELSE
            asc = 1
         END IF
      ELSEIF (is == watersurf) THEN
         IF (swe > 0) asc = 1
 
         !Vegetion surfaces
      ELSE
         IF (swe <= swed) THEN
 
            asc = 1 - ((1/3.1416)*acos(2*(swe/swed) - 1))**1.7
         ELSE
            asc = 1
         END IF
 
      END IF
 
      !asc=real(int(10000.*asc))/10000  !4 decimal precision
 
      RETURN
   END FUNCTION snowdepletioncurve
 
   SUBROUTINE veg_fr_snow( &
      sfr_surf, SnowFrac, & !input
      veg_fr) !output
 
      IMPLICIT NONE
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: sfr_surf
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowFrac
 
      REAL(KIND(1D0)), INTENT(out) :: veg_fr
 
      veg_fr = dot_product(sfr_surf(3:7), 1 - snowfrac(3:7))
 
   END SUBROUTINE veg_fr_snow
 
   !In this subroutine the snow properties are updated. The aging functions work for hourly air
   !temperature (dt=1). As the code timestep is typically smaller than one hour, the new albedo
   !and density are calculated using the timestep air temperature and then these are scaled to timestep
   !according to NSH.  Made by LJ in sprint 2014
   !Last update:
   ! TS 17 Sep 2017 - Improve the explicit interface
   ! LJ 7 July 2015 - Changed to work with shorter timestep: defined by tstep. Cleaning of the code.
   !
   !========================================================================
 
   SUBROUTINE snowupdate( &
      tstep, & !input
      Temp_C, &
      tau_a, &
      tau_f, &
      tau_r, &
      SnowDensMax, &
      SnowDensMin, &
      SnowAlbMax, &
      SnowAlbMin, &
      SnowPack_prev, SnowAlb_prev, SnowDens_prev, &
      SnowAlb_next, SnowDens_next & ! output
      )
 
      IMPLICIT NONE
 
      ! INTEGER, INTENT(in)::nsurf
      INTEGER, INTENT(in) :: tstep
 
      REAL(KIND(1D0)), INTENT(in) :: Temp_C !Air temperature
      REAL(KIND(1D0)), INTENT(in) :: tau_a
      REAL(KIND(1D0)), INTENT(in) :: tau_f
      REAL(KIND(1D0)), INTENT(in) :: tau_r
      REAL(KIND(1D0)), INTENT(in) :: SnowDensMax
      REAL(KIND(1D0)), INTENT(in) :: SnowDensMin
      REAL(KIND(1D0)), INTENT(in) :: SnowAlbMax
      REAL(KIND(1D0)), INTENT(in) :: SnowAlbMin
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowPack_prev
 
      REAL(KIND(1D0)), INTENT(in) :: SnowAlb_prev
      REAL(KIND(1D0)), INTENT(out) :: SnowAlb_next
 
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(in) :: SnowDens_prev
      REAL(KIND(1D0)), DIMENSION(nsurf), INTENT(out) :: SnowDens_next
 
      ! INTEGER::is
      REAL(KIND(1D0)) :: alb_change !Change in snow albedo
      REAL(KIND(1D0)) :: dens_change !Change in snow density
      REAL(KIND(1D0)), PARAMETER :: tau_1 = 24*60*60 !Number of seconds in a day
 
      !Initialize
      alb_change = 0
      dens_change = 0
 
      !==========================================================
      !Calculation of snow albedo by Lemonsu et al. 2010
      !(org: Verseghy (1991)&Baker et al.(1990))
      ! IF (SUM(SnowPack_prev) > 0) THEN !Check if snow on any of the surfaces
      !    IF (Temp_C < 0) THEN
      !       !alb_change = tau_a*(60*60)/tau_1
      !       alb_change = tau_a*(tstep)/tau_1
      !       SnowAlb_next = SnowAlb_prev - alb_change
      !    ELSE
      !       !alb_change = exp(-tau_f*(60*60)/tau_1)
      !       alb_change = EXP(-tau_f*(tstep)/tau_1)
      !       SnowAlb_next = (SnowAlb_prev - SnowAlbMin)*alb_change + SnowAlbMin
      !    ENDIF
      !    IF (SnowAlb_next < SnowAlbMin) SnowAlb_next = SnowAlbMin !Albedo cannot be smaller than the min albedo
      !    IF (SnowAlb_next > SnowAlbMax) SnowAlb_next = SnowAlbMax !Albedo cannot be larger than the max albedo
      !    if (SnowAlb_next < 0) print *, 'SnowAlbMin/max in SnowUpdate', SnowAlbMin, SnowAlbMax, SnowAlb_next
      ! ELSE
      !    SnowAlb_next = 0
      ! ENDIF
      ! if (SnowAlb_next < 0) print *, 'SnowAlb in SnowUpdate', SnowAlb_next
      snowalb_next = update_snow_albedo( &
                     tstep, snowpack_prev, snowalb_prev, temp_c, &
                     tau_a, tau_f, snowalbmax, snowalbmin)
 
      !Update snow density: There is a mistake in Järvi et al. (2014): tau_h should be tau_1
      ! DO is = 1, nsurf
      !    !If SnowPack existing
      !    IF (SnowPack_prev(is) > 0) THEN
      !       dens_change = EXP(-tau_r*(tstep)/tau_1)
      !       IF (SnowPack_prev(is) > 0) SnowDens_next(is) = (SnowDens_prev(is) - SnowDensMax)*dens_change + SnowDensMax
      !       IF (SnowDens_next(is) > SnowDensMax) SnowDens_next(is) = SnowDensMax
      !    ELSE
      !       SnowDens_next(is) = SnowDensMin
      !    ENDIF
      ! ENDDO
 
      snowdens_next = update_snow_dens( &
                      tstep, snowpack_prev, snowdens_prev, &
                      tau_r, snowdensmax, snowdensmin)
 
   END SUBROUTINE snowupdate
 
   FUNCTION update_snow_albedo( &
      tstep, SnowPack_prev, SnowAlb_prev, Temp_C, &
      tau_a, tau_f, SnowAlbMax, SnowAlbMin) &
      result(snowalb_next)
      IMPLICIT NONE
      INTEGER, INTENT(in) :: tstep
 
      REAL(kind(1d0)), DIMENSION(7), INTENT(in) :: snowpack_prev
      REAL(kind(1d0)), INTENT(in) :: temp_c
      REAL(kind(1d0)), INTENT(in) :: snowalb_prev
      REAL(kind(1d0)), INTENT(in) :: tau_a
      REAL(kind(1d0)), INTENT(in) :: tau_f
      REAL(kind(1d0)), INTENT(in) :: snowalbmax
      REAL(kind(1d0)), INTENT(in) :: snowalbmin
 
      REAL(kind(1d0)) :: snowalb_next
 
      REAL(kind(1d0)) :: alb_change !Change in snow albedo
      REAL(kind(1d0)), PARAMETER :: tau_1 = 24*60*60 !Number of seconds in a day
 
      !==========================================================
      !Calculation of snow albedo by Lemonsu et al. 2010
      !(org: Verseghy (1991)&Baker et al.(1990))
      IF (sum(snowpack_prev) > 0) THEN !Check if snow on any of the surfaces
         IF (temp_c < 0) THEN
            !alb_change = tau_a*(60*60)/tau_1
            alb_change = tau_a*(tstep)/tau_1
            snowalb_next = snowalb_prev - alb_change
         ELSE
            !alb_change = exp(-tau_f*(60*60)/tau_1)
            alb_change = exp(-tau_f*(tstep)/tau_1)
            snowalb_next = (snowalb_prev - snowalbmin)*alb_change + snowalbmin
         END IF
         IF (snowalb_next < snowalbmin) snowalb_next = snowalbmin !Albedo cannot be smaller than the min albedo
         IF (snowalb_next > snowalbmax) snowalb_next = snowalbmax !Albedo cannot be larger than the max albedo
         IF (snowalb_next < 0) print *, 'SnowAlbMin/max in SnowUpdate', snowalbmin, snowalbmax, snowalb_next
      ELSE
         snowalb_next = 0
      END IF
      IF (snowalb_next < 0) print *, 'SnowAlb in SnowUpdate', snowalb_next
 
   END FUNCTION update_snow_albedo
 
   FUNCTION update_snow_dens( &
      tstep, SnowPack_prev, SnowDens_prev, &
      tau_r, SnowDensMax, SnowDensMin) &
      result(snowdens_next)
      IMPLICIT NONE
      REAL(kind(1d0)), DIMENSION(nsurf), INTENT(in) :: snowpack_prev
      REAL(kind(1d0)), DIMENSION(nsurf), INTENT(in) :: snowdens_prev
      REAL(kind(1d0)), INTENT(in) :: snowdensmax
      REAL(kind(1d0)), INTENT(in) :: snowdensmin
      REAL(kind(1d0)), INTENT(in) :: tau_r
      INTEGER, INTENT(in) :: tstep
 
      ! REAL(KIND(1d0)), DIMENSION(nsurf), INTENT(out)::SnowDens_next
 
      REAL(kind(1d0)), DIMENSION(nsurf) :: snowdens_next
      REAL(kind(1d0)) :: dens_change
 
      INTEGER :: is
 
      REAL(kind(1d0)), PARAMETER :: tau_1 = 24*60*60
      ! INTEGER,parameter::nsurf=7
 
      !Update snow density: There is a mistake in Järvi et al. (2014): tau_h should be tau_1
      DO is = 1, nsurf
 
         !If SnowPack existing
         IF (snowpack_prev(is) > 0) THEN
            dens_change = exp(-tau_r*(tstep)/tau_1)
            IF (snowpack_prev(is) > 0) snowdens_next(is) = (snowdens_prev(is) - snowdensmax)*dens_change + snowdensmax
            IF (snowdens_next(is) > snowdensmax) snowdens_next(is) = snowdensmax
         ELSE
            snowdens_next(is) = snowdensmin
         END IF
      END DO
 
   END FUNCTION update_snow_dens
 
END MODULE snow_module