! SUMMA - Structure for Unifying Multiple Modeling Alternatives
! Copyright (C) 2014-2020 NCAR/RAL; University of Saskatchewan; University of Washington
!
! This file is part of SUMMA
!
! For more information see: http://www.ral.ucar.edu/projects/summa
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
module canopySnow_module
! data types
USE nrtype
! derived types to define the data structures
USE data_types,only:&
var_i, & ! data vector (i4b)
var_d, & ! data vector (dp)
var_dlength, & ! data vector with variable length dimension (dp)
model_options ! defines the model decisions
! physical constants
USE multiconst,only:Tfreeze ! freezing point of pure water (K)
! named variables defining elements in the data structures
USE var_lookup,only:iLookFORCE,iLookPARAM,iLookDIAG,iLookPROG,iLookFLUX ! named variables for structure elements
USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
! model decisions
USE mDecisions_module,only: &
stickySnow, & ! maximum interception capacity an increasing function of temerature
lightSnow, & ! maximum interception capacity an inverse function of new snow density
meltDripUnload, & ! Hedstrom and Pomeroy (1998), Storck et al 2002 (snowUnloadingCoeff & ratioDrip2Unloading)
windUnload ! Roesch et al 2001, formulate unloading based on wind and temperature
! privacy
implicit none
private
public::canopySnow
contains
! ************************************************************************************************
! public subroutine canopySnow: compute change in snow stored on the vegetation canopy
! ************************************************************************************************
subroutine canopySnow(&
! input: model control
dt, & ! intent(in): time step (seconds)
exposedVAI, & ! intent(in): exposed vegetation area index (m2 m-2)
computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
! input/output: data structures
model_decisions, & ! intent(in): model decisions
forc_data, & ! intent(in): model forcing data
mpar_data, & ! intent(in): model parameters
diag_data, & ! intent(in): model diagnostic variables for a local HRU
prog_data, & ! intent(inout): model prognostic variables for a local HRU
flux_data, & ! intent(inout): model flux variables
! output: error control
err,message) ! intent(out): error control
! ------------------------------------------------------------------------------------------------
implicit none
! ------------------------------------------------------------------------------------------------
! input: model control
real(dp),intent(in) :: dt ! time step (seconds)
real(dp),intent(in) :: exposedVAI ! exposed vegetation area index -- leaf + stem -- after burial by snow (m2 m-2)
logical(lgt),intent(in) :: computeVegFlux ! flag to indicate if we are computing fluxes over vegetation (.false. means veg is buried with snow)
! input/output: data structures
type(model_options),intent(in) :: model_decisions(:) ! model decisions
type(var_d),intent(in) :: forc_data ! model forcing data
type(var_dlength),intent(in) :: mpar_data ! model parameters
type(var_dlength),intent(in) :: diag_data ! model diagnostic variables for a local HRU
type(var_dlength),intent(inout) :: prog_data ! model prognostic variables for a local HRU
type(var_dlength),intent(inout) :: flux_data ! model flux variables
! output: error control
integer(i4b),intent(out) :: err ! error code
character(*),intent(out) :: message ! error message
! local variables
real(dp),parameter :: valueMissing=-9999._dp ! missing value
integer(i4b) :: iter ! iteration index
integer(i4b),parameter :: maxiter=50 ! maximum number of iterations
real(dp) :: unloading_melt ! unloading associated with canopy drip (kg m-2 s-1)
real(dp) :: airtemp_degC ! value of air temperature in degrees Celcius
real(dp) :: leafScaleFactor ! scaling factor for interception based on temperature (-)
real(dp) :: leafInterceptCapSnow ! storage capacity for snow per unit leaf area (kg m-2)
real(dp) :: canopyIceScaleFactor ! capacity scaling factor for throughfall (kg m-2)
real(dp) :: throughfallDeriv ! derivative in throughfall flux w.r.t. canopy storage (s-1)
real(dp) :: unloadingDeriv ! derivative in unloading flux w.r.t. canopy storage (s-1)
real(dp) :: scalarCanopyIceIter ! trial value for mass of ice on the vegetation canopy (kg m-2) (kg m-2)
real(dp) :: flux ! net flux (kg m-2 s-1)
real(dp) :: delS ! change in storage (kg m-2)
real(dp) :: resMass ! residual in mass equation (kg m-2)
real(dp) :: tempUnloadingFun ! temperature unloading functions, Eq. 14 in Roesch et al. 2001
real(dp) :: windUnloadingFun ! temperature unloading functions, Eq. 15 in Roesch et al. 2001
real(dp),parameter :: convTolerMass=0.0001_dp ! convergence tolerance for mass (kg m-2)
! -------------------------------------------------------------------------------------------------------------------------------
! initialize error control
err=0; message='canopySnow/'
! ------------------------------------------------------------------------------------------------
! associate variables in the data structure
associate(&
! model decisions
ixSnowInterception => model_decisions(iLookDECISIONS%snowIncept)%iDecision, & ! intent(in): [i4b] choice of option to determine maximum snow interception capacity
ixSnowUnload => model_decisions(iLookDECISIONS%snowUnload)%iDecision, & ! intent(in): [i4b] choice of option to determing how snow unloads from canopy
! model forcing data
scalarAirtemp => forc_data%var(iLookFORCE%airtemp), & ! intent(in): [dp] air temperature (K)
! model parameters
refInterceptCapSnow => mpar_data%var(iLookPARAM%refInterceptCapSnow)%dat(1), & ! intent(in): [dp] reference canopy interception capacity for snow per unit leaf area (kg m-2)
ratioDrip2Unloading => mpar_data%var(iLookPARAM%ratioDrip2Unloading)%dat(1), & ! intent(in): [dp] ratio of canopy drip to snow unloading (-)
snowUnloadingCoeff => mpar_data%var(iLookPARAM%snowUnloadingCoeff)%dat(1), & ! intent(in): [dp] time constant for unloading of snow from the forest canopy (s-1)
minTempUnloading => mpar_data%var(iLookPARAM%minTempUnloading)%dat(1), & ! constant describing the minimum temperature for snow unloading in windySnow parameterization (K)
minWindUnloading => mpar_data%var(iLookPARAM%minWindUnloading)%dat(1), & ! constant describing the minimum temperature for snow unloading in windySnow parameterization (K)
rateTempUnloading => mpar_data%var(iLookPARAM%rateTempUnloading)%dat(1), & ! constant describing how quickly snow will unload due to temperature in windySnow parameterization (K s)
rateWindUnloading => mpar_data%var(iLookPARAM%rateWindUnloading)%dat(1), & ! constant describing how quickly snow will unload due to wind in windySnow parameterization (K s)
! model diagnostic variables
scalarNewSnowDensity => diag_data%var(iLookDIAG%scalarNewSnowDensity)%dat(1), & ! intent(in): [dp] density of new snow (kg m-3)
! model prognostic variables (input/output)
scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1), & ! intent(inout): [dp] mass of ice on the vegetation canopy (kg m-2)
! model fluxes (input)
scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1), & ! intent(in): [dp] temperature of the canopy air space (k)
scalarSnowfall => flux_data%var(iLookFLUX%scalarSnowfall)%dat(1), & ! intent(in): [dp] computed snowfall rate (kg m-2 s-1)
scalarCanopyLiqDrainage => flux_data%var(iLookFLUX%scalarCanopyLiqDrainage)%dat(1), & ! intent(in): [dp] liquid drainage from the vegetation canopy (kg m-2 s-1)
scalarWindspdCanopyTop => flux_data%var(iLookFLUX%scalarWindspdCanopyTop)%dat(1), & ! intent(in): [dp] windspeed at the top of the canopy (m s-1)
! model variables (output)
scalarThroughfallSnow => flux_data%var(iLookFLUX%scalarThroughfallSnow)%dat(1), & ! intent(out): [dp] snow that reaches the ground without ever touching the canopy (kg m-2 s-1)
scalarCanopySnowUnloading => flux_data%var(iLookFLUX%scalarCanopySnowUnloading)%dat(1) & ! intent(out): [dp] unloading of snow from the vegetion canopy (kg m-2 s-1)
) ! associate variables in the data structures
! -----------------------------------------------------------------------------------------------------------------------------------------------------
! compute unloading due to melt drip...
! *************************************
if(computeVegFlux)then
unloading_melt = min(ratioDrip2Unloading*scalarCanopyLiqDrainage, scalarCanopyIce/dt) ! kg m-2 s-1
else
unloading_melt = 0._dp
end if
scalarCanopyIce = scalarCanopyIce - unloading_melt*dt
! *****
! compute the ice balance due to snowfall and unloading...
! ********************************************************
! check for early returns
if(.not.computeVegFlux .or. (scalarSnowfall<tiny(dt) .and. scalarCanopyIce<tiny(dt)))then
scalarThroughfallSnow = scalarSnowfall ! throughfall of snow through the canopy (kg m-2 s-1)
scalarCanopySnowUnloading = unloading_melt ! unloading of snow from the canopy (kg m-2 s-1)
return
end if
! get a trial value for canopy storage
scalarCanopyIceIter = scalarCanopyIce
do iter=1,maxiter
! ** compute unloading
if (ixSnowUnload==meltDripUnload) then
scalarCanopySnowUnloading = snowUnloadingCoeff*scalarCanopyIceIter
unloadingDeriv = snowUnloadingCoeff
else if (ixSnowUnload==windUnload) then
tempUnloadingFun = max(scalarCanairTemp - minTempUnloading, 0._dp) / rateTempUnloading ! (s-1)
if (scalarWindspdCanopyTop >= minWindUnloading) then
windUnloadingFun = abs(scalarWindspdCanopyTop) / rateWindUnloading ! (s-1)
else
windUnloadingFun = 0._dp ! (s-1)
end if
! implement the "windySnow" Roesch et al. 2001 parameterization, Eq. 13 in Roesch et al. 2001
scalarCanopySnowUnloading = scalarCanopyIceIter * (tempUnloadingFun + windUnloadingFun)
unloadingDeriv = tempUnloadingFun + windUnloadingFun
end if
! no snowfall
if(scalarSnowfall<tiny(dt))then ! no snow
scalarThroughfallSnow = scalarSnowfall ! throughfall (kg m-2 s-1)
canopyIceScaleFactor = valueMissing ! not used
throughfallDeriv = 0._dp
else
! ** process different options for maximum branch snow interception
select case(ixSnowInterception)
case(lightSnow)
! (check new snow density is valid)
if(scalarNewSnowDensity < 0._dp)then; err=20; message=trim(message)//'invalid new snow density'; return; end if
! (compute storage capacity of new snow)
leafScaleFactor = 0.27_dp + 46._dp/scalarNewSnowDensity
leafInterceptCapSnow = refInterceptCapSnow*leafScaleFactor ! per unit leaf area (kg m-2)
case(stickySnow)
airtemp_degC = scalarAirtemp - Tfreeze
if (airtemp_degC > -1._dp) then
leafScaleFactor = 4.0_dp
elseif(airtemp_degC > -3._dp) then
leafScaleFactor = 1.5_dp*airtemp_degC + 5.5_dp
else
leafScaleFactor = 1.0_dp
end if
leafInterceptCapSnow = refInterceptCapSnow*leafScaleFactor
case default
message=trim(message)//'unable to identify option for maximum branch interception capacity'
err=20; return
end select
! compute maximum interception capacity for the canopy
canopyIceScaleFactor = leafInterceptCapSnow*exposedVAI
! (compute throughfall)
scalarThroughfallSnow = scalarSnowfall*(scalarCanopyIceIter/canopyIceScaleFactor)
throughfallDeriv = scalarSnowfall/canopyIceScaleFactor
end if ! (if snow is falling)
! ** compute iteration increment
flux = scalarSnowfall - scalarThroughfallSnow - scalarCanopySnowUnloading ! net flux (kg m-2 s-1)
delS = (flux*dt - (scalarCanopyIceIter - scalarCanopyIce))/(1._dp + (throughfallDeriv + unloadingDeriv)*dt)
! ** check for convergence
resMass = scalarCanopyIceIter - (scalarCanopyIce + flux*dt)
if(abs(resMass) < convTolerMass)exit
! ** check for non-convengence
if(iter==maxiter)then; err=20; message=trim(message)//'failed to converge [mass]'; return; end if
! ** update value
scalarCanopyIceIter = scalarCanopyIceIter + delS
end do ! iterating
! add the unloading associated with melt drip (kg m-2 s-1)
scalarCanopySnowUnloading = scalarCanopySnowUnloading + unloading_melt
! *****
! update mass of ice on the canopy (kg m-2)
scalarCanopyIce = scalarCanopyIceIter
! end association to variables in the data structure
end associate
end subroutine canopySnow
end module canopySnow_module