From 90be91bb198538caa2c6e351aca79047b44b7794 Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Tue, 27 Sep 2022 18:30:55 +0000
Subject: [PATCH] copied from non-actors summa
---
.../source/engine/sundials/computHeatCap.f90 | 961 +++++++++---------
1 file changed, 480 insertions(+), 481 deletions(-)
diff --git a/build/source/engine/sundials/computHeatCap.f90 b/build/source/engine/sundials/computHeatCap.f90
index 53997c8..b1f242c 100644
--- a/build/source/engine/sundials/computHeatCap.f90
+++ b/build/source/engine/sundials/computHeatCap.f90
@@ -20,492 +20,491 @@
module computHeatCap_module
- ! data types
- USE nrtype
-
- ! derived types to define the data structures
- USE data_types,only:&
- var_d, & ! data vector (rkind)
- var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength ! data vector with variable length dimension (rkind)
-
- ! named variables defining elements in the data structures
- USE var_lookup,only:iLookPARAM,iLookDIAG,iLookINDEX ! named variables for structure elements
-
- ! physical constants
- USE multiconst,only:&
- Tfreeze, & ! freezing point of water (K)
- iden_air, & ! intrinsic density of air (kg m-3)
- iden_ice, & ! intrinsic density of ice (kg m-3)
- iden_water, & ! intrinsic density of water (kg m-3)
- ! specific heat
- Cp_air, & ! specific heat of air (J kg-1 K-1)
- Cp_ice, & ! specific heat of ice (J kg-1 K-1)
- Cp_soil, & ! specific heat of soil (J kg-1 K-1)
- Cp_water ! specific heat of liquid water (J kg-1 K-1)
- ! named variables to describe the state variable type
- USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space
- USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy
- USE globalData,only:iname_watCanopy ! named variable defining the mass of total water on the vegetation canopy
- USE globalData,only:iname_liqCanopy ! named variable defining the mass of liquid water on the vegetation canopy
- USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers
- USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers
- USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers
- USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers
- USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers
- USE globalData,only:iname_watAquifer ! named variable defining the water storage in the aquifer
-
- ! missing values
- USE globalData,only:integerMissing ! missing integer
-
- ! named variables that define the layer type
- USE globalData,only:iname_snow ! snow
- USE globalData,only:iname_soil ! soil
-
-
- ! privacy
- implicit none
- private
- public::computHeatCap
- public::computStatMult
- public::computHeatCapAnalytic
- public::computCm
-
- contains
-
-
- ! **********************************************************************************************************
- ! public subroutine computHeatCap: compute diagnostic energy variables (heat capacity)
- ! **********************************************************************************************************
- subroutine computHeatCap(&
- ! input: control variables
- nLayers, & ! intent(in): number of layers (soil+snow)
- computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
- canopyDepth, & ! intent(in): canopy depth (m)
- ! input data structures
- mpar_data, & ! intent(in): model parameters
- indx_data, & ! intent(in): model layer indices
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- ! input: state variables
- scalarCanopyIce, & ! intent(in)
- scalarCanopyLiquid, & ! intent(in)
- scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K)
- scalarCanopyTempPrev, & ! intent(in): previous value of canopy temperature (K)
- scalarCanopyEnthalpyTrial, & ! intent(in): trial enthalpy of the vegetation canopy (J m-3)
- scalarCanopyEnthalpyPrev, & ! intent(in): previous enthalpy of the vegetation canopy (J m-3)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-)
- mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
- mLayerTempTrial, & ! intent(in): trial temperature
- mLayerTempPrev, & ! intent(in): previous temperature
- mLayerEnthalpyTrial, & ! intent(in): trial enthalpy for snow and soil
- mLayerEnthalpyPrev, & ! intent(in): previous enthalpy for snow and soil
- ! output
- heatCapVeg, &
- mLayerHeatCap, & ! intent(out): heat capacity for snow and soil
- ! output: error control
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------------
- ! input: control variables
- logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
- real(rkind),intent(in) :: canopyDepth ! depth of the vegetation canopy (m)
- ! input/output: data structures
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_ilength),intent(in) :: indx_data ! model layer indices
- ! input:
- integer(i4b),intent(in) :: nLayers
- type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
- real(rkind),intent(in) :: scalarCanopyIce ! trial value of canopy ice content (kg m-2)
- real(rkind),intent(in) :: scalarCanopyLiquid
- real(rkind),intent(in) :: scalarCanopyTempTrial ! trial value of canopy temperature
- real(rkind),intent(in) :: scalarCanopyEnthalpyTrial ! trial enthalpy of the vegetation canopy (J m-3)
- real(rkind),intent(in) :: scalarCanopyEnthalpyPrev ! intent(in): previous enthalpy of the vegetation canopy (J m-3)
- real(rkind),intent(in) :: scalarCanopyTempPrev ! Previous value of canopy temperature
- real(rkind),intent(in) :: mLayerVolFracLiq(:) ! trial vector of volumetric liquid water content (-)
- real(rkind),intent(in) :: mLayerVolFracIce(:) ! trial vector of volumetric ice water content (-)
- real(rkind),intent(in) :: mLayerTempTrial(:)
- real(rkind),intent(in) :: mLayerTempPrev(:)
- real(rkind),intent(in) :: mLayerEnthalpyTrial(:)
- real(rkind),intent(in) :: mLayerEnthalpyPrev(:)
- ! output:
- real(qp),intent(out) :: heatCapVeg
- real(qp),intent(out) :: mLayerHeatCap(:)
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- integer(i4b) :: iLayer ! index of model layer
- real(rkind) :: delT
- real(rkind) :: delEnt
- integer(i4b) :: iSoil ! index of soil layer
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! associate variables in data structure
- associate(&
- ! input: coordinate variables
- nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
- layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
- ! input: heat capacity and thermal conductivity
- specificHeatVeg => mpar_data%var(iLookPARAM%specificHeatVeg)%dat(1), & ! intent(in): specific heat of vegetation (J kg-1 K-1)
- maxMassVegetation => mpar_data%var(iLookPARAM%maxMassVegetation)%dat(1), & ! intent(in): maximum mass of vegetation (kg m-2)
- ! input: depth varying soil parameters
- iden_soil => mpar_data%var(iLookPARAM%soil_dens_intr)%dat, & ! intent(in): intrinsic density of soil (kg m-3)
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat & ! intent(in): soil porosity (-)
- ) ! end associate statemen
- ! initialize error control
- err=0; message="computHeatCap/"
-
- ! initialize the soil layer
- iSoil=integerMissing
-
- ! compute the bulk volumetric heat capacity of vegetation (J m-3 K-1)
- if(computeVegFlux)then
- delT = scalarCanopyTempTrial - scalarCanopyTempPrev
- if(abs(delT) <= 1e-14_rkind)then
- heatCapVeg = specificHeatVeg*maxMassVegetation/canopyDepth + & ! vegetation component
+! data types
+USE nrtype
+
+! derived types to define the data structures
+USE data_types,only:&
+ var_d, & ! data vector (rkind)
+ var_ilength, & ! data vector with variable length dimension (i4b)
+ var_dlength ! data vector with variable length dimension (rkind)
+
+! named variables defining elements in the data structures
+USE var_lookup,only:iLookPARAM,iLookDIAG,iLookINDEX ! named variables for structure elements
+
+! physical constants
+USE multiconst,only:&
+ Tfreeze, & ! freezing point of water (K)
+ iden_air, & ! intrinsic density of air (kg m-3)
+ iden_ice, & ! intrinsic density of ice (kg m-3)
+ iden_water, & ! intrinsic density of water (kg m-3)
+ ! specific heat
+ Cp_air, & ! specific heat of air (J kg-1 K-1)
+ Cp_ice, & ! specific heat of ice (J kg-1 K-1)
+ Cp_soil, & ! specific heat of soil (J kg-1 K-1)
+ Cp_water ! specific heat of liquid water (J kg-1 K-1)
+! named variables to describe the state variable type
+USE globalData,only:iname_nrgCanair ! named variable defining the energy of the canopy air space
+USE globalData,only:iname_nrgCanopy ! named variable defining the energy of the vegetation canopy
+USE globalData,only:iname_watCanopy ! named variable defining the mass of total water on the vegetation canopy
+USE globalData,only:iname_liqCanopy ! named variable defining the mass of liquid water on the vegetation canopy
+USE globalData,only:iname_nrgLayer ! named variable defining the energy state variable for snow+soil layers
+USE globalData,only:iname_watLayer ! named variable defining the total water state variable for snow+soil layers
+USE globalData,only:iname_liqLayer ! named variable defining the liquid water state variable for snow+soil layers
+USE globalData,only:iname_matLayer ! named variable defining the matric head state variable for soil layers
+USE globalData,only:iname_lmpLayer ! named variable defining the liquid matric potential state variable for soil layers
+USE globalData,only:iname_watAquifer ! named variable defining the water storage in the aquifer
+
+! missing values
+USE globalData,only:integerMissing ! missing integer
+
+! named variables that define the layer type
+USE globalData,only:iname_snow ! snow
+USE globalData,only:iname_soil ! soil
+
+
+! privacy
+implicit none
+private
+public::computHeatCap
+public::computStatMult
+public::computHeatCapAnalytic
+public::computCm
+
+contains
+
+
+! **********************************************************************************************************
+! public subroutine computHeatCap: compute diagnostic energy variables (heat capacity)
+! **********************************************************************************************************
+subroutine computHeatCap(&
+ ! input: control variables
+ nLayers, & ! intent(in): number of layers (soil+snow)
+ computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
+ canopyDepth, & ! intent(in): canopy depth (m)
+ ! input data structures
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(in): model layer indices
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ ! input: state variables
+ scalarCanopyIce, & ! intent(in)
+ scalarCanopyLiquid, & ! intent(in)
+ scalarCanopyTempTrial, & ! intent(in): trial value of canopy temperature (K)
+ scalarCanopyTempPrev, & ! intent(in): previous value of canopy temperature (K)
+ scalarCanopyEnthalpyTrial, & ! intent(in): trial enthalpy of the vegetation canopy (J m-3)
+ scalarCanopyEnthalpyPrev, & ! intent(in): previous enthalpy of the vegetation canopy (J m-3)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-)
+ mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
+ mLayerTempTrial, & ! intent(in): trial temperature
+ mLayerTempPrev, & ! intent(in): previous temperature
+ mLayerEnthalpyTrial, & ! intent(in): trial enthalpy for snow and soil
+ mLayerEnthalpyPrev, & ! intent(in): previous enthalpy for snow and soil
+ ! output
+ heatCapVeg, &
+ mLayerHeatCap, & ! intent(out): heat capacity for snow and soil
+ ! output: error control
+ err,message) ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------------------------------------
+ ! input: control variables
+ logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
+ real(rkind),intent(in) :: canopyDepth ! depth of the vegetation canopy (m)
+ ! input/output: data structures
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(in) :: indx_data ! model layer indices
+ ! input:
+ integer(i4b),intent(in) :: nLayers
+ type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
+ real(rkind),intent(in) :: scalarCanopyIce ! trial value of canopy ice content (kg m-2)
+ real(rkind),intent(in) :: scalarCanopyLiquid
+ real(rkind),intent(in) :: scalarCanopyTempTrial ! trial value of canopy temperature
+ real(rkind),intent(in) :: scalarCanopyEnthalpyTrial ! trial enthalpy of the vegetation canopy (J m-3)
+ real(rkind),intent(in) :: scalarCanopyEnthalpyPrev ! intent(in): previous enthalpy of the vegetation canopy (J m-3)
+ real(rkind),intent(in) :: scalarCanopyTempPrev ! Previous value of canopy temperature
+ real(rkind),intent(in) :: mLayerVolFracLiq(:) ! trial vector of volumetric liquid water content (-)
+ real(rkind),intent(in) :: mLayerVolFracIce(:) ! trial vector of volumetric ice water content (-)
+ real(rkind),intent(in) :: mLayerTempTrial(:)
+ real(rkind),intent(in) :: mLayerTempPrev(:)
+ real(rkind),intent(in) :: mLayerEnthalpyTrial(:)
+ real(rkind),intent(in) :: mLayerEnthalpyPrev(:)
+ ! output:
+ real(qp),intent(out) :: heatCapVeg
+ real(qp),intent(out) :: mLayerHeatCap(:)
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ integer(i4b) :: iLayer ! index of model layer
+ real(rkind) :: delT
+ real(rkind) :: delEnt
+ integer(i4b) :: iSoil ! index of soil layer
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! associate variables in data structure
+ associate(&
+ ! input: coordinate variables
+ nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
+ layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
+ ! input: heat capacity and thermal conductivity
+ specificHeatVeg => mpar_data%var(iLookPARAM%specificHeatVeg)%dat(1), & ! intent(in): specific heat of vegetation (J kg-1 K-1)
+ maxMassVegetation => mpar_data%var(iLookPARAM%maxMassVegetation)%dat(1), & ! intent(in): maximum mass of vegetation (kg m-2)
+ ! input: depth varying soil parameters
+ iden_soil => mpar_data%var(iLookPARAM%soil_dens_intr)%dat, & ! intent(in): intrinsic density of soil (kg m-3)
+ theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat & ! intent(in): soil porosity (-)
+ ) ! end associate statemen
+ ! initialize error control
+ err=0; message="computHeatCap/"
+
+ ! initialize the soil layer
+ iSoil=integerMissing
+
+ ! compute the bulk volumetric heat capacity of vegetation (J m-3 K-1)
+ if(computeVegFlux)then
+ delT = scalarCanopyTempTrial - scalarCanopyTempPrev
+ if(abs(delT) <= 1e-14_rkind)then
+ heatCapVeg = specificHeatVeg*maxMassVegetation/canopyDepth + & ! vegetation component
Cp_water*scalarCanopyLiquid/canopyDepth + & ! liquid water component
Cp_ice*scalarCanopyIce/canopyDepth ! ice component
- else
- delEnt = scalarCanopyEnthalpyTrial - scalarCanopyEnthalpyPrev
- heatCapVeg = delEnt / delT
- end if
- end if
-
- ! loop through layers
- do iLayer=1,nLayers
- delT = mLayerTempTrial(iLayer) - mLayerTempPrev(iLayer)
- if(abs(delT) <= 1e-14_rkind)then
+ else
+ delEnt = scalarCanopyEnthalpyTrial - scalarCanopyEnthalpyPrev
+ heatCapVeg = delEnt / delT
+ end if
+ end if
+
+ ! loop through layers
+ do iLayer=1,nLayers
+ delT = mLayerTempTrial(iLayer) - mLayerTempPrev(iLayer)
+ if(abs(delT) <= 1e-14_rkind)then
! get the soil layer
if(iLayer>nSnow) iSoil = iLayer-nSnow
- select case(layerType(iLayer))
- ! * soil
- case(iname_soil)
- mLayerHeatCap(iLayer) = iden_soil(iSoil) * Cp_soil * ( 1._rkind - theta_sat(iSoil) ) + & ! soil component
- iden_ice * Cp_Ice * mLayerVolFracIce(iLayer) + & ! ice component
- iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
- iden_air * Cp_air * ( theta_sat(iSoil) - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) )! air component
- case(iname_snow)
- mLayerHeatCap(iLayer) = iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
- iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
- iden_air * Cp_air * ( 1._rkind - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) ) ! air component
- case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute olumetric heat capacity'; return
- end select
- else
+ select case(layerType(iLayer))
+ ! * soil
+ case(iname_soil)
+ mLayerHeatCap(iLayer) = iden_soil(iSoil) * Cp_soil * ( 1._rkind - theta_sat(iSoil) ) + & ! soil component
+ iden_ice * Cp_Ice * mLayerVolFracIce(iLayer) + & ! ice component
+ iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
+ iden_air * Cp_air * ( theta_sat(iSoil) - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) )! air component
+ case(iname_snow)
+ mLayerHeatCap(iLayer) = iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
+ iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
+ iden_air * Cp_air * ( 1._rkind - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) ) ! air component
+ case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute olumetric heat capacity'; return
+ end select
+ else
delEnt = mLayerEnthalpyTrial(iLayer) - mLayerEnthalpyPrev(iLayer)
mLayerHeatCap(iLayer) = delEnt / delT
- endif
- end do ! looping through layers
-
- end associate
-
- end subroutine computHeatCap
-
- ! **********************************************************************************************************
- ! public subroutine computStatMult: get scale factors
- ! **********************************************************************************************************
- subroutine computStatMult(&
- heatCapVeg, &
- mLayerHeatCap, &
- ! input: data structures
- diag_data, & ! intent(in): model diagnostic variables for a local HRU
- indx_data, & ! intent(in): indices defining model states and layers
- ! output
- sMul, & ! intent(out): multiplier for state vector (used in the residual calculations)
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------
- USE nr_utility_module,only:arth ! get a sequence of numbers arth(start, incr, count)
- USE f2008funcs_module,only:findIndex ! finds the index of the first value within a vector
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! input: data structures
- real(qp),intent(out) :: heatCapVeg
- real(qp),intent(out) :: mLayerHeatCap(:)
- type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
- type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers
- ! output: state vectors
- real(qp),intent(out) :: sMul(:) ! NOTE: qp ! multiplier for state vector (used in the residual calculations)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! state subsets
- integer(i4b) :: iLayer ! index of layer within the snow+soil domain
- integer(i4b) :: ixStateSubset ! index within the state subset
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! make association with variables in the data structures
- associate(&
- ! model diagnostic variables
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
- volHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1),& ! intent(in) : [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1)
- ! indices defining specific model states
- ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy air space energy state variable
- ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy energy state variable
- ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy hydrology state variable (mass)
- ! vector of energy and hydrology indices for the snow and soil domains
- ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in) : [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
- ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in) : [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
- nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in) : [i4b] number of energy state variables in the snow+soil domain
- nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in) : [i4b] number of hydrology state variables in the snow+soil domain
- ! type of model state variabless
- ixStateType_subset => indx_data%var(iLookINDEX%ixStateType_subset)%dat ,& ! intent(in) : [i4b(:)] [state subset] type of desired model state variables
- ! number of layers
- nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1) ,& ! intent(in) : [i4b] number of snow layers
- nSoil => indx_data%var(iLookINDEX%nSoil)%dat(1) ,& ! intent(in) : [i4b] number of soil layers
- nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1) & ! intent(in) : [i4b] total number of layers
- ) ! end association with variables in the data structures
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message='computStatMult/'
-
- ! -----
- ! * define components of derivative matrices that are constant over a time step (substep)...
- ! ------------------------------------------------------------------------------------------
-
- ! define the multiplier for the state vector for residual calculations (vegetation canopy)
- ! NOTE: Use the "where" statement to generalize to multiple canopy layers (currently one canopy layer)
-
- where(ixStateType_subset==iname_nrgCanair) sMul = Cp_air*iden_air ! volumetric heat capacity of air (J m-3 K-1)
- where(ixStateType_subset==iname_nrgCanopy) sMul = heatCapVeg ! volumetric heat capacity of the vegetation (J m-3 K-1)
- where(ixStateType_subset==iname_watCanopy) sMul = 1._rkind ! nothing else on the left hand side
- where(ixStateType_subset==iname_liqCanopy) sMul = 1._rkind ! nothing else on the left hand side
-
-
- ! define the energy multiplier for the state vector for residual calculations (snow-soil domain)
- if(nSnowSoilNrg>0)then
- do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain)
- ixStateSubset = ixSnowSoilNrg(iLayer) ! index within the state vector
- sMul(ixStateSubset) = mLayerHeatCap(iLayer) ! transfer volumetric heat capacity to the state multiplier
- end do ! looping through non-missing energy state variables in the snow+soil domain
- endif
-
- ! define the hydrology multiplier and diagonal elements for the state vector for residual calculations (snow-soil domain)
- if(nSnowSoilHyd>0)then
- do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain)
- ixStateSubset = ixSnowSoilHyd(iLayer) ! index within the state vector
- sMul(ixStateSubset) = 1._rkind ! state multiplier = 1 (nothing else on the left-hand-side)
- end do ! looping through non-missing energy state variables in the snow+soil domain
- endif
-
- ! define the scaling factor and diagonal elements for the aquifer
- where(ixStateType_subset==iname_watAquifer) sMul = 1._rkind
-
- ! ------------------------------------------------------------------------------------------
- ! ------------------------------------------------------------------------------------------
-
- end associate
- ! end association to variables in the data structure where vector length does not change
- end subroutine computStatMult
-
- ! **********************************************************************************************************
- ! public subroutine computHeatCapAnalytic: compute diagnostic energy variables (heat capacity)
- ! **********************************************************************************************************
- subroutine computHeatCapAnalytic(&
- ! input: control variables
- computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
- canopyDepth, & ! intent(in): canopy depth (m)
- ! input: state variables
- scalarCanopyIce, & ! intent(in)
- scalarCanopyLiquid, & ! intent(in)
- mLayerVolFracIce, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-)
- mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
- ! input data structures
- mpar_data, & ! intent(in): model parameters
- indx_data, & ! intent(in): model layer indices
- ! output
- heatCapVeg, &
- mLayerHeatCap, &
- ! output: error control
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------------
- ! --------------------------------------------------------------------------------------------------------------------------------------
- ! input: model control
- logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
- real(rkind),intent(in) :: canopyDepth ! depth of the vegetation canopy (m)
- real(rkind),intent(in) :: scalarCanopyIce ! trial value of canopy ice content (kg m-2)
- real(rkind),intent(in) :: scalarCanopyLiquid
- real(rkind),intent(in) :: mLayerVolFracLiq(:) ! trial vector of volumetric liquid water content (-)
- real(rkind),intent(in) :: mLayerVolFracIce(:) ! trial vector of volumetric ice water content (-)
- ! input/output: data structures
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_ilength),intent(in) :: indx_data ! model layer indices
- ! output: error control
- real(qp),intent(out) :: heatCapVeg
- real(qp),intent(out) :: mLayerHeatCap(:)
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- character(LEN=256) :: cmessage ! error message of downwind routine
- integer(i4b) :: iLayer ! index of model layer
- integer(i4b) :: iSoil ! index of soil layer
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! associate variables in data structure
- associate(&
- ! input: coordinate variables
- nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
- nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1), & ! intent(in): total number of layers
- layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
- ! input: heat capacity and thermal conductivity
- specificHeatVeg => mpar_data%var(iLookPARAM%specificHeatVeg)%dat(1), & ! intent(in): specific heat of vegetation (J kg-1 K-1)
- maxMassVegetation => mpar_data%var(iLookPARAM%maxMassVegetation)%dat(1), & ! intent(in): maximum mass of vegetation (kg m-2)
- ! input: depth varying soil parameters
- iden_soil => mpar_data%var(iLookPARAM%soil_dens_intr)%dat, & ! intent(in): intrinsic density of soil (kg m-3)
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat & ! intent(in): soil porosity (-)
- ) ! end associate statement
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="computHeatCapAnalytic/"
-
- ! initialize the soil layer
- iSoil=integerMissing
-
- ! compute the bulk volumetric heat capacity of vegetation (J m-3 K-1)
- if(computeVegFlux)then
- heatCapVeg = specificHeatVeg*maxMassVegetation/canopyDepth + & ! vegetation component
- Cp_water*scalarCanopyLiquid/canopyDepth + & ! liquid water component
- Cp_ice*scalarCanopyIce/canopyDepth ! ice component
- end if
-
- ! loop through layers
- do iLayer=1,nLayers
-
- ! get the soil layer
- if(iLayer>nSnow) iSoil = iLayer-nSnow
-
- ! *****
- ! * compute the volumetric heat capacity of each layer (J m-3 K-1)...
- ! *******************************************************************
- select case(layerType(iLayer))
- ! * soil
- case(iname_soil)
- mLayerHeatCap(iLayer) = iden_soil(iSoil) * Cp_soil * ( 1._rkind - theta_sat(iSoil) ) + & ! soil component
- iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
- iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
- iden_air * Cp_air * ( theta_sat(iSoil) - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) )! air component
- case(iname_snow)
- mLayerHeatCap(iLayer) = iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
- iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
- iden_air * Cp_air * ( 1._rkind - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) ) ! air component
- case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute olumetric heat capacity'; return
- end select
-
- end do ! looping through layers
- !pause
-
- ! end association to variables in the data structure
- end associate
-
- end subroutine computHeatCapAnalytic
-
- ! **********************************************************************************************************
- ! public subroutine computCm
- ! **********************************************************************************************************
- subroutine computCm(&
- ! input: control variables
- computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
- ! input: state variables
- scalarCanopyTemp, & ! intent(in)
- mLayerTemp, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
- mLayerMatricHead, & ! intent(in)
- ! input data structures
- mpar_data, & ! intent(in): model parameters
- indx_data, & ! intent(in): model layer indices
- ! output
- scalarCanopyCm, & ! intent(out): Cm for vegetation
- mLayerCm, & ! intent(out): Cm for soil and snow
- ! output: error control
- err,message) ! intent(out): error control
- ! --------------------------------------------------------------------------------------------------------------------------------------
- ! provide access to external subroutines
- USE soil_utils_module,only:crit_soilT ! compute critical temperature below which ice exists
- ! --------------------------------------------------------------------------------------------------------------------------------------
- ! input: model control
- logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
- real(rkind),intent(in) :: scalarCanopyTemp ! value of canopy ice content (kg m-2)
- real(rkind),intent(in) :: mLayerTemp(:) ! vector of volumetric liquid water content (-)
- real(rkind),intent(in) :: mLayerMatricHead(:) ! vector of total water matric potential (m)
- ! input/output: data structures
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_ilength),intent(in) :: indx_data ! model layer indices
- ! output: error control
- real(qp),intent(out) :: scalarCanopyCm
- real(qp),intent(out) :: mLayerCm(:)
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! local variables
- character(LEN=256) :: cmessage ! error message of downwind routine
- integer(i4b) :: iLayer ! index of model layer
- integer(i4b) :: iSoil ! index of soil layer
- real(rkind) :: g1
- real(rkind) :: g2
- real(rkind) :: Tcrit ! temperature where all water is unfrozen (K)
-
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! associate variables in data structure
- associate(&
- ! input: coordinate variables
- nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
- nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1), & ! intent(in): total number of layers
- layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
- snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) & ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
- ) ! end associate statement
- ! --------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="computCm/"
-
- ! initialize the soil layer
- iSoil=integerMissing
-
- ! compute Cm of vegetation
- ! Note that scalarCanopyCm/iden_water is computed
- if(computeVegFlux)then
- g2 = scalarCanopyTemp - Tfreeze
- g1 = (1._rkind/snowfrz_scale) * atan(snowfrz_scale * g2)
- if(scalarCanopyTemp < Tfreeze)then
- scalarCanopyCm = Cp_water * g1 + Cp_ice * (g2 - g1)
- else
- scalarCanopyCm = Cp_water * g2
- end if
- end if
-
- ! loop through layers
- do iLayer=1,nLayers
-
- ! get the soil layer
- if(iLayer>nSnow) iSoil = iLayer-nSnow
-
- ! *****
- ! * compute Cm of of each layer
- ! *******************************************************************
- select case(layerType(iLayer))
- ! * soil
- case(iname_soil)
- g2 = mLayerTemp(iLayer) - Tfreeze
- Tcrit = crit_soilT( mLayerMatricHead(iSoil) )
- if( mLayerTemp(iLayer) < Tcrit)then
- mLayerCm(iLayer) = (iden_ice * Cp_ice - iden_air * Cp_air) * g2
- else
- mLayerCm(iLayer) = (iden_water * Cp_water - iden_air * Cp_air) * g2
- end if
-
- case(iname_snow)
- g2 = mLayerTemp(iLayer) - Tfreeze
+ endif
+ end do ! looping through layers
+
+ end associate
+
+end subroutine computHeatCap
+
+! **********************************************************************************************************
+! public subroutine computStatMult: get scale factors
+! **********************************************************************************************************
+subroutine computStatMult(&
+ heatCapVeg, &
+ mLayerHeatCap, &
+ ! input: data structures
+ diag_data, & ! intent(in): model diagnostic variables for a local HRU
+ indx_data, & ! intent(in): indices defining model states and layers
+ ! output
+ sMul, & ! intent(out): multiplier for state vector (used in the residual calculations)
+ err,message) ! intent(out): error control
+! --------------------------------------------------------------------------------------------------------------------------------
+USE nr_utility_module,only:arth ! get a sequence of numbers arth(start, incr, count)
+USE f2008funcs_module,only:findIndex ! finds the index of the first value within a vector
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! input: data structures
+ real(qp),intent(out) :: heatCapVeg
+ real(qp),intent(out) :: mLayerHeatCap(:)
+ type(var_dlength),intent(in) :: diag_data ! diagnostic variables for a local HRU
+ type(var_ilength),intent(in) :: indx_data ! indices defining model states and layers
+ ! output: state vectors
+ real(qp),intent(out) :: sMul(:) ! NOTE: qp ! multiplier for state vector (used in the residual calculations)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! state subsets
+ integer(i4b) :: iLayer ! index of layer within the snow+soil domain
+ integer(i4b) :: ixStateSubset ! index within the state subset
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! make association with variables in the data structures
+ associate(&
+ ! model diagnostic variables
+ canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
+ volHeatCapVeg => diag_data%var(iLookDIAG%scalarBulkVolHeatCapVeg)%dat(1),& ! intent(in) : [dp] bulk volumetric heat capacity of vegetation (J m-3 K-1)
+ ! indices defining specific model states
+ ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy air space energy state variable
+ ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy energy state variable
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat ,& ! intent(in) : [i4b(:)] [length=1] index of canopy hydrology state variable (mass)
+ ! vector of energy and hydrology indices for the snow and soil domains
+ ixSnowSoilNrg => indx_data%var(iLookINDEX%ixSnowSoilNrg)%dat ,& ! intent(in) : [i4b(:)] index in the state subset for energy state variables in the snow+soil domain
+ ixSnowSoilHyd => indx_data%var(iLookINDEX%ixSnowSoilHyd)%dat ,& ! intent(in) : [i4b(:)] index in the state subset for hydrology state variables in the snow+soil domain
+ nSnowSoilNrg => indx_data%var(iLookINDEX%nSnowSoilNrg )%dat(1) ,& ! intent(in) : [i4b] number of energy state variables in the snow+soil domain
+ nSnowSoilHyd => indx_data%var(iLookINDEX%nSnowSoilHyd )%dat(1) ,& ! intent(in) : [i4b] number of hydrology state variables in the snow+soil domain
+ ! type of model state variabless
+ ixStateType_subset => indx_data%var(iLookINDEX%ixStateType_subset)%dat ,& ! intent(in) : [i4b(:)] [state subset] type of desired model state variables
+ ! number of layers
+ nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1) ,& ! intent(in) : [i4b] number of snow layers
+ nSoil => indx_data%var(iLookINDEX%nSoil)%dat(1) ,& ! intent(in) : [i4b] number of soil layers
+ nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1) & ! intent(in) : [i4b] total number of layers
+ ) ! end association with variables in the data structures
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message='computStatMult/'
+
+ ! -----
+ ! * define components of derivative matrices that are constant over a time step (substep)...
+ ! ------------------------------------------------------------------------------------------
+
+ ! define the multiplier for the state vector for residual calculations (vegetation canopy)
+ ! NOTE: Use the "where" statement to generalize to multiple canopy layers (currently one canopy layer)
+
+ where(ixStateType_subset==iname_nrgCanair) sMul = Cp_air*iden_air ! volumetric heat capacity of air (J m-3 K-1)
+ where(ixStateType_subset==iname_nrgCanopy) sMul = heatCapVeg ! volumetric heat capacity of the vegetation (J m-3 K-1)
+ where(ixStateType_subset==iname_watCanopy) sMul = 1._rkind ! nothing else on the left hand side
+ where(ixStateType_subset==iname_liqCanopy) sMul = 1._rkind ! nothing else on the left hand side
+
+
+ ! define the energy multiplier for the state vector for residual calculations (snow-soil domain)
+ if(nSnowSoilNrg>0)then
+ do concurrent (iLayer=1:nLayers,ixSnowSoilNrg(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain)
+ ixStateSubset = ixSnowSoilNrg(iLayer) ! index within the state vector
+ sMul(ixStateSubset) = mLayerHeatCap(iLayer) ! transfer volumetric heat capacity to the state multiplier
+ end do ! looping through non-missing energy state variables in the snow+soil domain
+ endif
+
+ ! define the hydrology multiplier and diagonal elements for the state vector for residual calculations (snow-soil domain)
+ if(nSnowSoilHyd>0)then
+ do concurrent (iLayer=1:nLayers,ixSnowSoilHyd(iLayer)/=integerMissing) ! (loop through non-missing energy state variables in the snow+soil domain)
+ ixStateSubset = ixSnowSoilHyd(iLayer) ! index within the state vector
+ sMul(ixStateSubset) = 1._rkind ! state multiplier = 1 (nothing else on the left-hand-side)
+ end do ! looping through non-missing energy state variables in the snow+soil domain
+ endif
+
+ ! define the scaling factor and diagonal elements for the aquifer
+ where(ixStateType_subset==iname_watAquifer) sMul = 1._rkind
+
+ ! ------------------------------------------------------------------------------------------
+ ! ------------------------------------------------------------------------------------------
+
+ end associate
+! end association to variables in the data structure where vector length does not change
+end subroutine computStatMult
+
+! **********************************************************************************************************
+! public subroutine computHeatCapAnalytic: compute diagnostic energy variables (heat capacity)
+! **********************************************************************************************************
+subroutine computHeatCapAnalytic(&
+ ! input: control variables
+ computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
+ canopyDepth, & ! intent(in): canopy depth (m)
+ ! input: state variables
+ scalarCanopyIce, & ! intent(in)
+ scalarCanopyLiquid, & ! intent(in)
+ mLayerVolFracIce, & ! intent(in): volumetric fraction of ice at the start of the sub-step (-)
+ mLayerVolFracLiq, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
+ ! input data structures
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(in): model layer indices
+ ! output
+ heatCapVeg, &
+ mLayerHeatCap, &
+ ! output: error control
+ err,message) ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------------------------------------
+ ! --------------------------------------------------------------------------------------------------------------------------------------
+ ! input: model control
+ logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
+ real(rkind),intent(in) :: canopyDepth ! depth of the vegetation canopy (m)
+ real(rkind),intent(in) :: scalarCanopyIce ! trial value of canopy ice content (kg m-2)
+ real(rkind),intent(in) :: scalarCanopyLiquid
+ real(rkind),intent(in) :: mLayerVolFracLiq(:) ! trial vector of volumetric liquid water content (-)
+ real(rkind),intent(in) :: mLayerVolFracIce(:) ! trial vector of volumetric ice water content (-)
+ ! input/output: data structures
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(in) :: indx_data ! model layer indices
+ ! output: error control
+ real(qp),intent(out) :: heatCapVeg
+ real(qp),intent(out) :: mLayerHeatCap(:)
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ integer(i4b) :: iLayer ! index of model layer
+ integer(i4b) :: iSoil ! index of soil layer
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! associate variables in data structure
+ associate(&
+ ! input: coordinate variables
+ nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
+ nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1), & ! intent(in): total number of layers
+ layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
+ ! input: heat capacity and thermal conductivity
+ specificHeatVeg => mpar_data%var(iLookPARAM%specificHeatVeg)%dat(1), & ! intent(in): specific heat of vegetation (J kg-1 K-1)
+ maxMassVegetation => mpar_data%var(iLookPARAM%maxMassVegetation)%dat(1), & ! intent(in): maximum mass of vegetation (kg m-2)
+ ! input: depth varying soil parameters
+ iden_soil => mpar_data%var(iLookPARAM%soil_dens_intr)%dat, & ! intent(in): intrinsic density of soil (kg m-3)
+ theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat & ! intent(in): soil porosity (-)
+ ) ! end associate statement
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message="computHeatCapAnalytic/"
+
+ ! initialize the soil layer
+ iSoil=integerMissing
+
+ ! compute the bulk volumetric heat capacity of vegetation (J m-3 K-1)
+ if(computeVegFlux)then
+ heatCapVeg = specificHeatVeg*maxMassVegetation/canopyDepth + & ! vegetation component
+ Cp_water*scalarCanopyLiquid/canopyDepth + & ! liquid water component
+ Cp_ice*scalarCanopyIce/canopyDepth ! ice component
+ end if
+
+ ! loop through layers
+ do iLayer=1,nLayers
+
+ ! get the soil layer
+ if(iLayer>nSnow) iSoil = iLayer-nSnow
+
+ ! *****
+ ! * compute the volumetric heat capacity of each layer (J m-3 K-1)...
+ ! *******************************************************************
+ select case(layerType(iLayer))
+ ! * soil
+ case(iname_soil)
+ mLayerHeatCap(iLayer) = iden_soil(iSoil) * Cp_soil * ( 1._rkind - theta_sat(iSoil) ) + & ! soil component
+ iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
+ iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
+ iden_air * Cp_air * ( theta_sat(iSoil) - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) )! air component
+ case(iname_snow)
+ mLayerHeatCap(iLayer) = iden_ice * Cp_ice * mLayerVolFracIce(iLayer) + & ! ice component
+ iden_water * Cp_water * mLayerVolFracLiq(iLayer) + & ! liquid water component
+ iden_air * Cp_air * ( 1._rkind - (mLayerVolFracIce(iLayer) + mLayerVolFracLiq(iLayer)) ) ! air component
+ case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute olumetric heat capacity'; return
+ end select
+
+ end do ! looping through layers
+ !pause
+
+ ! end association to variables in the data structure
+ end associate
+end subroutine computHeatCapAnalytic
+
+! **********************************************************************************************************
+! public subroutine computCm
+! **********************************************************************************************************
+subroutine computCm(&
+ ! input: control variables
+ computeVegFlux, & ! intent(in): flag to denote if computing the vegetation flux
+ ! input: state variables
+ scalarCanopyTemp, & ! intent(in)
+ mLayerTemp, & ! intent(in): volumetric fraction of liquid water at the start of the sub-step (-)
+ mLayerMatricHead, & ! intent(in)
+ ! input data structures
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(in): model layer indices
+ ! output
+ scalarCanopyCm, & ! intent(out): Cm for vegetation
+ mLayerCm, & ! intent(out): Cm for soil and snow
+ ! output: error control
+ err,message) ! intent(out): error control
+ ! --------------------------------------------------------------------------------------------------------------------------------------
+ ! provide access to external subroutines
+ USE soil_utils_module,only:crit_soilT ! compute critical temperature below which ice exists
+ ! --------------------------------------------------------------------------------------------------------------------------------------
+ ! input: model control
+ logical(lgt),intent(in) :: computeVegFlux ! logical flag to denote if computing the vegetation flux
+ real(rkind),intent(in) :: scalarCanopyTemp ! value of canopy ice content (kg m-2)
+ real(rkind),intent(in) :: mLayerTemp(:) ! vector of volumetric liquid water content (-)
+ real(rkind),intent(in) :: mLayerMatricHead(:) ! vector of total water matric potential (m)
+ ! input/output: data structures
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(in) :: indx_data ! model layer indices
+ ! output: error control
+ real(qp),intent(out) :: scalarCanopyCm
+ real(qp),intent(out) :: mLayerCm(:)
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! local variables
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ integer(i4b) :: iLayer ! index of model layer
+ integer(i4b) :: iSoil ! index of soil layer
+ real(rkind) :: g1
+ real(rkind) :: g2
+ real(rkind) :: Tcrit ! temperature where all water is unfrozen (K)
+
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! associate variables in data structure
+ associate(&
+ ! input: coordinate variables
+ nSnow => indx_data%var(iLookINDEX%nSnow)%dat(1), & ! intent(in): number of snow layers
+ nLayers => indx_data%var(iLookINDEX%nLayers)%dat(1), & ! intent(in): total number of layers
+ layerType => indx_data%var(iLookINDEX%layerType)%dat, & ! intent(in): layer type (iname_soil or iname_snow)
+ snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) & ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
+ ) ! end associate statement
+ ! --------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message="computCm/"
+
+ ! initialize the soil layer
+ iSoil=integerMissing
+
+ ! compute Cm of vegetation
+ ! Note that scalarCanopyCm/iden_water is computed
+ if(computeVegFlux)then
+ g2 = scalarCanopyTemp - Tfreeze
g1 = (1._rkind/snowfrz_scale) * atan(snowfrz_scale * g2)
- mLayerCm(iLayer) = (iden_ice * Cp_ice - iden_air * Cp_air * iden_water/iden_ice) * ( g2 - g1 ) &
- + (iden_water * Cp_water - iden_air * Cp_air) * g1
-
- case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute Cm'; return
- end select
-
- end do ! looping through layers
- !pause
-
- ! end association to variables in the data structure
- end associate
-
- end subroutine computCm
-
-
- end module computHeatCap_module
+ if(scalarCanopyTemp < Tfreeze)then
+ scalarCanopyCm = Cp_water * g1 + Cp_ice * (g2 - g1)
+ else
+ scalarCanopyCm = Cp_water * g2
+ end if
+ end if
+
+ ! loop through layers
+ do iLayer=1,nLayers
+
+ ! get the soil layer
+ if(iLayer>nSnow) iSoil = iLayer-nSnow
+
+ ! *****
+ ! * compute Cm of of each layer
+ ! *******************************************************************
+ select case(layerType(iLayer))
+ ! * soil
+ case(iname_soil)
+ g2 = mLayerTemp(iLayer) - Tfreeze
+ Tcrit = crit_soilT( mLayerMatricHead(iSoil) )
+ if( mLayerTemp(iLayer) < Tcrit)then
+ mLayerCm(iLayer) = (iden_ice * Cp_ice - iden_air * Cp_air) * g2
+ else
+ mLayerCm(iLayer) = (iden_water * Cp_water - iden_air * Cp_air) * g2
+ end if
+
+ case(iname_snow)
+ g2 = mLayerTemp(iLayer) - Tfreeze
+ g1 = (1._rkind/snowfrz_scale) * atan(snowfrz_scale * g2)
+ mLayerCm(iLayer) = (iden_ice * Cp_ice - iden_air * Cp_air * iden_water/iden_ice) * ( g2 - g1 ) &
+ + (iden_water * Cp_water - iden_air * Cp_air) * g1
+
+ case default; err=20; message=trim(message)//'unable to identify type of layer (snow or soil) to compute Cm'; return
+ end select
+
+ end do ! looping through layers
+ !pause
+
+ ! end association to variables in the data structure
+ end associate
+
+end subroutine computCm
+
+
+end module computHeatCap_module
\ No newline at end of file
--
GitLab