From 8bccf32a29727d895f52343ea5f346d9bf8c804c Mon Sep 17 00:00:00 2001
From: Kyle <kyle.c.klenk@gmail.com>
Date: Thu, 29 Sep 2022 23:16:35 +0000
Subject: [PATCH] laugh tests produce the same result as summa-sundials
---
build/source/dshare/globalData.f90 | 2 +-
build/source/engine/coupled_em.f90 | 25 +-
build/source/engine/mDecisions.f90 | 13 +-
build/source/engine/nrtype.mod | Bin 769 -> 0 bytes
build/source/engine/opSplittin.f90 | 2149 ++++++++++++++--------------
build/source/engine/snow_utils.f90 | 120 +-
build/source/engine/soil_utils.f90 | 1218 ++++++++--------
build/source/engine/updatState.f90 | 242 ++--
8 files changed, 1879 insertions(+), 1890 deletions(-)
delete mode 100644 build/source/engine/nrtype.mod
diff --git a/build/source/dshare/globalData.f90 b/build/source/dshare/globalData.f90
index 1e1cbf7..205788e 100755
--- a/build/source/dshare/globalData.f90
+++ b/build/source/dshare/globalData.f90
@@ -315,7 +315,7 @@ MODULE globalData
! define fixed dimensions
integer(i4b),parameter,public :: nBand=2 ! number of spectral bands
- integer(i4b),parameter,public :: nTimeDelay=3000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4)
+ integer(i4b),parameter,public :: nTimeDelay=2000 ! number of hours in the time delay histogram (default: ~1 season = 24*365/4)
character(len=1024),public,save :: fname ! temporary filename
diff --git a/build/source/engine/coupled_em.f90 b/build/source/engine/coupled_em.f90
index f11a0a8..62bb6d2 100755
--- a/build/source/engine/coupled_em.f90
+++ b/build/source/engine/coupled_em.f90
@@ -65,9 +65,8 @@ USE globalData,only:globalPrintFlag ! the global print flag
! look-up values for the numerical method
USE mDecisions_module,only: &
- iterative, & ! iterative
- nonIterative, & ! non-iterative
- iterSurfEnergyBal ! iterate only on the surface energy balance
+ bEuler, & ! home-grown backward Euler solution with long time steps
+ sundials ! SUNDIALS/IDA solution
! look-up values for the maximum interception capacity
USE mDecisions_module,only: &
@@ -143,14 +142,6 @@ subroutine coupled_em(&
USE var_derive_module,only:calcHeight ! module to calculate height at layer interfaces and layer mid-point
USE computSnowDepth_module,only:computSnowDepth
! look-up values for the numerical method
- USE mDecisions_module,only: &
- iterative, & ! iterative
- nonIterative, & ! non-iterative
- iterSurfEnergyBal ! iterate only on the surface energy balance
- ! look-up values for the maximum interception capacity
- USE mDecisions_module,only: &
- stickySnow, & ! maximum interception capacity an increasing function of temerature
- lightSnow ! maximum interception capacity an inverse function of new snow density
implicit none
! model control
integer(4),intent(in) :: indxHRU ! hruId
@@ -253,14 +244,6 @@ subroutine coupled_em(&
! initialize error control
err=0; message="coupled_em/"
- ! print*, "scalarCanopyWat"
- ! print*, "mLayerVolFracWat"
- ! print*, "mLayerMatricHead"
- ! print*, ""
- ! print*, ""
- ! print*, ""
- ! print*, ""
-
! check that the decision is supported
if(model_decisions(iLookDECISIONS%groundwatr)%iDecision==bigBucket .and. &
model_decisions(iLookDECISIONS%spatial_gw)%iDecision/=localColumn)then
@@ -1177,7 +1160,9 @@ subroutine coupled_em(&
! adjust length of the sub-step (make sure that we don't exceed the step)
! TODO: This is different in Ashley's code
- dt_sub = min(data_step - dt_solv, dt_sub)
+ ! dt_sub = min(data_step - dt_solv, dt_sub)
+ dt_sub = data_step - dt_solv !min(data_step - dt_solv, dt_sub)
+
!print*, 'dt_sub = ', dt_sub
end do substeps ! (sub-step loop)
diff --git a/build/source/engine/mDecisions.f90 b/build/source/engine/mDecisions.f90
index 72258c0..4ead60c 100755
--- a/build/source/engine/mDecisions.f90
+++ b/build/source/engine/mDecisions.f90
@@ -156,7 +156,7 @@ contains
! ************************************************************************************************
! public subroutine mDecisions: save model decisions as named integers
! ************************************************************************************************
-subroutine mDecisions(err,message)
+subroutine mDecisions(num_steps,err,message)
! model time structures
USE multiconst,only:secprday ! number of seconds in a day
USE var_lookup,only:iLookTIME ! named variables that identify indices in the time structures
@@ -181,10 +181,11 @@ subroutine mDecisions(err,message)
USE time_utils_module,only:extractTime ! extract time info from units string
USE time_utils_module,only:compjulday ! compute the julian day
USE time_utils_module,only:fracDay ! compute fractional day
- USE summaFileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module
+ USE summaActors_FileManager,only: SIM_START_TM, SIM_END_TM ! time info from control file module
implicit none
! define output
+ integer(i4b),intent(out) :: num_steps
integer(i4b),intent(out) :: err ! error code
character(*),intent(out) :: message ! error message
! define local variables
@@ -287,8 +288,8 @@ subroutine mDecisions(err,message)
oldTime%var(:) = startTime%var(:)
! compute the number of time steps
- numtim = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1
- write(*,'(a,1x,i10)') 'number of time steps = ', numtim
+ num_steps = nint( (dJulianFinsh - dJulianStart)*secprday/data_step ) + 1
+ numTim = num_steps
! set Noah-MP options
@@ -675,8 +676,8 @@ subroutine readoption(err,message)
USE ascii_util_module,only:file_open ! open file
USE ascii_util_module,only:linewidth ! max character number for one line
USE ascii_util_module,only:get_vlines ! get a vector of non-comment lines
- USE summaFileManager,only:SETTINGS_PATH ! path for metadata files
- USE summaFileManager,only:M_DECISIONS ! definition of modeling options
+ USE summaActors_FileManager,only:SETTINGS_PATH ! path for metadata files
+ USE summaActors_FileManager,only:M_DECISIONS ! definition of modeling options
USE get_ixname_module,only:get_ixdecisions ! identify index of named variable
USE globalData,only:model_decisions ! model decision structure
implicit none
diff --git a/build/source/engine/nrtype.mod b/build/source/engine/nrtype.mod
deleted file mode 100644
index cb753b75bed2441f7d23ae5b89340b7015376d91..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 769
zcmV+c1OEIUiwFP!000006Wy3wZ<{a_fZzQq&YR&Om2ED8x1xkrXe~sPZ0b{#22!D=
zNek(w{rerRO~Ykvw@P6lQO*Ga{%oI{FE2bl@GTb}cH#3Tf*(<m#=9+0u2EPe(Q+S!
zuugUluub;A9;36h$%rJ$@Grik5@dp6%TE9eA1Tn{-6&AtAW1XJx32BL_F?XRalLQu
z$Z-SP4Qh4fdy`5b&h_ldnK;4d#`<al16I5=yX6uep=PJ>sC4X2u4io~o_iBm_+W|j
zH>U>8B}p#tEJY?0C1Y3(96|*e;qpj63B$LotGZLU(hg9yBS-BhJ2RD)svGQ_buwr`
zV}3V?iQt_OIzotHK!{F=F+!-rfKWRj44~iRZHS$h;_W_(w`shRUr?~8_H^#s$0+#c
zJMN7$$ze!&wQ7}j;bQJi0>^W~@_p-e;mvVG7E^m>yHnep+!nn!yaA8;L5G+!7)PBp
zfo|^i10Q+Lvi>h{)T8?-d5F_A-rm2RiL^1n`<Y0S=?_RF6F)6jq`=}i*3<Jc{9LBo
zaT@L3mK=NAmUeW+&cv9lS^YV0OctppOXO6TNo@WnD&+v;i00ax&aWDF-)cw)M%-%E
z&~DWYgTV0o7@uxzR)ngo_0@TISmRnoun0esp<`_{7>pi<rC}J_Qek~4dWNzSr#l_o
z6nuTB_c4qMh7V5&`3qxx6(l;XwZhO;RTUvO`e-mzKMlQ%;264-GT|69hrp|8w24+&
zSY{!7@qr|YyiTHOFZu2yioDI=iUtIR>U0<@symGmB7+!4W!gi5kvM%y_NVM`^(8Z-
zK7T7PfbW0qPC16_sxgkjJ~u1065^zzM!}(R&hAvbFOFyXPx<Z>15s|_Fi`4!5QgJa
z>?bhL%A%ZsURJgY#$~0&z$9(e3WKtYd6n<uS>DW-@vI>Ah@z$x`5sDYl{cZRRfz}6
z`bz*uZCnIElnhcu5md~=iy<H{T>*-swkK4SdR?NT9?UPp>6XcFLzIljJ`?}|GAngq
diff --git a/build/source/engine/opSplittin.f90 b/build/source/engine/opSplittin.f90
index ba0cbca..5fcc364 100755
--- a/build/source/engine/opSplittin.f90
+++ b/build/source/engine/opSplittin.f90
@@ -20,1094 +20,1093 @@
module opSplittin_module
- ! data types
- USE nrtype
-
- ! access the global print flag
- USE globalData,only:globalPrintFlag
-
- ! access missing values
- USE globalData,only:integerMissing ! missing integer
- USE globalData,only:realMissing ! missing double precision number
- USE globalData,only:quadMissing ! missing quadruple precision number
-
- ! access matrix information
- USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix
- USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix
- USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix
- USE globalData,only: iJac1 ! first layer of the Jacobian to print
- USE globalData,only: iJac2 ! last layer of the Jacobian to print
-
- ! domain types
- USE globalData,only:iname_cas ! named variables for the canopy air space
- USE globalData,only:iname_veg ! named variables for vegetation
- USE globalData,only:iname_snow ! named variables for snow
- USE globalData,only:iname_soil ! named variables for soil
-
- ! 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
-
- ! global metadata
- USE globalData,only:flux_meta ! metadata on the model fluxes
- USE globalData,only:diag_meta ! metadata on the model diagnostic variables
- USE globalData,only:prog_meta ! metadata on the model prognostic variables
- USE globalData,only:deriv_meta ! metadata on the model derivatives
- USE globalData,only:flux2state_orig ! metadata on flux-to-state mapping (original state variables)
- USE globalData,only:flux2state_liq ! metadata on flux-to-state mapping (liquid water state variables)
-
- ! constants
- USE multiconst,only:&
- gravity, & ! acceleration of gravity (m s-2)
- Tfreeze, & ! temperature at freezing (K)
- LH_fus, & ! latent heat of fusion (J kg-1)
- LH_vap, & ! latent heat of vaporization (J kg-1)
- LH_sub, & ! latent heat of sublimation (J kg-1)
- Cp_air, & ! specific heat of air (J kg-1 K-1)
- iden_air, & ! intrinsic density of air (kg m-3)
- iden_ice, & ! intrinsic density of ice (kg m-3)
- iden_water ! intrinsic density of liquid water (kg m-3)
-
- ! provide access to indices that define elements of the data structures
- USE var_lookup,only:iLookATTR ! named variables for structure elements
- USE var_lookup,only:iLookTYPE ! named variables for structure elements
- USE var_lookup,only:iLookPROG ! named variables for structure elements
- USE var_lookup,only:iLookDIAG ! named variables for structure elements
- USE var_lookup,only:iLookFLUX ! named variables for structure elements
- USE var_lookup,only:iLookFORCE ! named variables for structure elements
- USE var_lookup,only:iLookPARAM ! named variables for structure elements
- USE var_lookup,only:iLookINDEX ! named variables for structure elements
- USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
-
- ! look up structure for variable types
- USE var_lookup,only:iLookVarType
-
- ! provide access to the number of flux variables
- USE var_lookup,only:nFlux=>maxvarFlux ! number of model flux variables
-
- ! provide access to the derived types to define the data structures
- USE data_types,only:&
- var_i, & ! data vector (i4b)
- var_d, & ! data vector (rkind)
- var_flagVec, & ! data vector with variable length dimension (i4b)
- var_ilength, & ! data vector with variable length dimension (i4b)
- var_dlength, & ! data vector with variable length dimension (rkind)
- zLookup, & ! data vector with variable length dimension (rkind)
- model_options ! defines the model decisions
-
- ! look-up values for the choice of groundwater representation (local-column, or single-basin)
- USE mDecisions_module,only: &
- localColumn, & ! separate groundwater representation in each local soil column
- singleBasin ! single groundwater store over the entire basin
-
- ! look-up values for the choice of groundwater parameterization
- USE mDecisions_module,only: &
- qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
- bigBucket, & ! a big bucket (lumped aquifer model)
- noExplicit, & ! no explicit groundwater parameterization
- sundials, & ! SUNDIALS/IDA solution
- bEuler ! home-grown backward Euler solution with long time step
-
- ! safety: set private unless specified otherwise
+! data types
+USE nrtype
+
+! access the global print flag
+USE globalData,only:globalPrintFlag
+
+! access missing values
+USE globalData,only:integerMissing ! missing integer
+USE globalData,only:realMissing ! missing double precision number
+USE globalData,only:quadMissing ! missing quadruple precision number
+
+! access matrix information
+USE globalData,only: nBands ! length of the leading dimension of the band diagonal matrix
+USE globalData,only: ixFullMatrix ! named variable for the full Jacobian matrix
+USE globalData,only: ixBandMatrix ! named variable for the band diagonal matrix
+USE globalData,only: iJac1 ! first layer of the Jacobian to print
+USE globalData,only: iJac2 ! last layer of the Jacobian to print
+
+! domain types
+USE globalData,only:iname_cas ! named variables for the canopy air space
+USE globalData,only:iname_veg ! named variables for vegetation
+USE globalData,only:iname_snow ! named variables for snow
+USE globalData,only:iname_soil ! named variables for soil
+
+! 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
+
+! global metadata
+USE globalData,only:flux_meta ! metadata on the model fluxes
+USE globalData,only:diag_meta ! metadata on the model diagnostic variables
+USE globalData,only:prog_meta ! metadata on the model prognostic variables
+USE globalData,only:deriv_meta ! metadata on the model derivatives
+USE globalData,only:flux2state_orig ! metadata on flux-to-state mapping (original state variables)
+USE globalData,only:flux2state_liq ! metadata on flux-to-state mapping (liquid water state variables)
+
+! constants
+USE multiconst,only:&
+ gravity, & ! acceleration of gravity (m s-2)
+ Tfreeze, & ! temperature at freezing (K)
+ LH_fus, & ! latent heat of fusion (J kg-1)
+ LH_vap, & ! latent heat of vaporization (J kg-1)
+ LH_sub, & ! latent heat of sublimation (J kg-1)
+ Cp_air, & ! specific heat of air (J kg-1 K-1)
+ iden_air, & ! intrinsic density of air (kg m-3)
+ iden_ice, & ! intrinsic density of ice (kg m-3)
+ iden_water ! intrinsic density of liquid water (kg m-3)
+
+! provide access to indices that define elements of the data structures
+USE var_lookup,only:iLookATTR ! named variables for structure elements
+USE var_lookup,only:iLookTYPE ! named variables for structure elements
+USE var_lookup,only:iLookPROG ! named variables for structure elements
+USE var_lookup,only:iLookDIAG ! named variables for structure elements
+USE var_lookup,only:iLookFLUX ! named variables for structure elements
+USE var_lookup,only:iLookFORCE ! named variables for structure elements
+USE var_lookup,only:iLookPARAM ! named variables for structure elements
+USE var_lookup,only:iLookINDEX ! named variables for structure elements
+USE var_lookup,only:iLookDECISIONS ! named variables for elements of the decision structure
+
+! look up structure for variable types
+USE var_lookup,only:iLookVarType
+
+! provide access to the number of flux variables
+USE var_lookup,only:nFlux=>maxvarFlux ! number of model flux variables
+
+! provide access to the derived types to define the data structures
+USE data_types,only:&
+ var_i, & ! data vector (i4b)
+ var_d, & ! data vector (rkind)
+ var_flagVec, & ! data vector with variable length dimension (i4b)
+ var_ilength, & ! data vector with variable length dimension (i4b)
+ var_dlength, & ! data vector with variable length dimension (rkind)
+ zLookup, & ! data vector with variable length dimension (rkind)
+ model_options ! defines the model decisions
+
+! look-up values for the choice of groundwater representation (local-column, or single-basin)
+USE mDecisions_module,only: &
+ localColumn, & ! separate groundwater representation in each local soil column
+ singleBasin ! single groundwater store over the entire basin
+
+! look-up values for the choice of groundwater parameterization
+USE mDecisions_module,only: &
+ qbaseTopmodel, & ! TOPMODEL-ish baseflow parameterization
+ bigBucket, & ! a big bucket (lumped aquifer model)
+ noExplicit, & ! no explicit groundwater parameterization
+ sundials, & ! SUNDIALS/IDA solution
+ bEuler ! home-grown backward Euler solution with long time step
+
+! safety: set private unless specified otherwise
+implicit none
+private
+public::opSplittin
+
+! named variables for the coupling method
+integer(i4b),parameter :: fullyCoupled=1 ! 1st try: fully coupled solution
+integer(i4b),parameter :: stateTypeSplit=2 ! 2nd try: separate solutions for each state type
+
+! named variables for the state variable split
+integer(i4b),parameter :: nrgSplit=1 ! order in sequence for the energy operation
+integer(i4b),parameter :: massSplit=2 ! order in sequence for the mass operation
+
+! named variables for the domain type split
+integer(i4b),parameter :: vegSplit=1 ! order in sequence for the vegetation split
+integer(i4b),parameter :: snowSplit=2 ! order in sequence for the snow split
+integer(i4b),parameter :: soilSplit=3 ! order in sequence for the soil split
+integer(i4b),parameter :: aquiferSplit=4 ! order in sequence for the aquifer split
+
+! named variables for the solution method
+integer(i4b),parameter :: vector=1 ! vector solution method
+integer(i4b),parameter :: scalar=2 ! scalar solution method
+integer(i4b),parameter :: nSolutions=2 ! number of solution methods
+
+! named variables for the switch between states and domains
+integer(i4b),parameter :: fullDomain=1 ! full domain (veg+snow+soil)
+integer(i4b),parameter :: subDomain=2 ! sub domain (veg, snow, and soil separately)
+
+! maximum number of possible splits
+integer(i4b),parameter :: nStateTypes=2 ! number of state types (energy, water)
+integer(i4b),parameter :: nDomains=4 ! number of domains (vegetation, snow, soil, and aquifer)
+
+! control parameters
+real(rkind),parameter :: valueMissing=-9999._rkind ! missing value
+real(rkind),parameter :: verySmall=1.e-12_rkind ! a very small number (used to check consistency)
+real(rkind),parameter :: veryBig=1.e+20_rkind ! a very big number
+real(rkind),parameter :: dx = 1.e-8_rkind ! finite difference increment
+
+contains
+
+
+! **********************************************************************************************************
+! public subroutine opSplittin: run the coupled energy-mass model for one timestep
+!
+! The logic of the solver is as follows:
+! (1) Attempt different solutions in the following order: (a) fully coupled; (b) split by state type and by
+! domain type for a given energy and mass split (vegetation, snow, and soil); and (c) scalar solution
+! for a given state type and domain subset.
+! (2) For a given split, compute a variable number of substeps (in varSubstepSundials).
+! **********************************************************************************************************
+subroutine opSplittin(&
+ ! input: model control
+ nSnow, & ! intent(in): number of snow layers
+ nSoil, & ! intent(in): number of soil layers
+ nLayers, & ! intent(in): total number of layers
+ nState, & ! intent(in): total number of state variables
+ dt, & ! intent(inout): time step (s)
+ firstSubStep, & ! intent(in): flag to denote first sub-step
+ computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
+ ! input/output: data structures
+ type_data, & ! intent(in): type of vegetation and soil
+ attr_data, & ! intent(in): spatial attributes
+ forc_data, & ! intent(in): model forcing data
+ mpar_data, & ! intent(in): model parameters
+ indx_data, & ! intent(inout): index data
+ prog_data, & ! intent(inout): model prognostic variables for a local HRU
+ diag_data, & ! intent(inout): model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout): model fluxes for a local HRU
+ bvar_data, & ! intent(in): model variables for the local basin
+ lookup_data, & ! intent(in): lookup tables
+ model_decisions,& ! intent(in): model decisions
+ ! output: model control
+ dtMultiplier, & ! intent(out): substep multiplier (-)
+ tooMuchMelt, & ! intent(out): flag to denote that ice is insufficient to support melt
+ stepFailure, & ! intent(out): flag to denote step failure
+ ixCoupling, & ! intent(out): coupling method used in this iteration
+ err,message) ! intent(out): error code and error message
+ ! ---------------------------------------------------------------------------------------
+ ! structure allocations
+ USE allocspace_module,only:allocLocal ! allocate local data structures
+ ! simulation of fluxes and residuals given a trial state vector
+ USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content
+ USE soil_utils_module,only:liquidHead ! compute the liquid water matric potential
+ ! population/extraction of state vectors
+ USE indexState_module,only:indexSplit ! get state indices
+ USE varSubstep_module,only:varSubstep ! complete substeps for a given split
+ USE varSubstepSundials_module,only:varSubstepSundials ! complete substeps for a given split
+ ! identify name of variable type (for error message)
+ USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages
implicit none
- private
- public::opSplittin
-
- ! named variables for the coupling method
- integer(i4b),parameter :: fullyCoupled=1 ! 1st try: fully coupled solution
- integer(i4b),parameter :: stateTypeSplit=2 ! 2nd try: separate solutions for each state type
-
- ! named variables for the state variable split
- integer(i4b),parameter :: nrgSplit=1 ! order in sequence for the energy operation
- integer(i4b),parameter :: massSplit=2 ! order in sequence for the mass operation
-
- ! named variables for the domain type split
- integer(i4b),parameter :: vegSplit=1 ! order in sequence for the vegetation split
- integer(i4b),parameter :: snowSplit=2 ! order in sequence for the snow split
- integer(i4b),parameter :: soilSplit=3 ! order in sequence for the soil split
- integer(i4b),parameter :: aquiferSplit=4 ! order in sequence for the aquifer split
-
- ! named variables for the solution method
- integer(i4b),parameter :: vector=1 ! vector solution method
- integer(i4b),parameter :: scalar=2 ! scalar solution method
- integer(i4b),parameter :: nSolutions=2 ! number of solution methods
-
- ! named variables for the switch between states and domains
- integer(i4b),parameter :: fullDomain=1 ! full domain (veg+snow+soil)
- integer(i4b),parameter :: subDomain=2 ! sub domain (veg, snow, and soil separately)
-
- ! maximum number of possible splits
- integer(i4b),parameter :: nStateTypes=2 ! number of state types (energy, water)
- integer(i4b),parameter :: nDomains=4 ! number of domains (vegetation, snow, soil, and aquifer)
-
- ! control parameters
- real(rkind),parameter :: valueMissing=-9999._rkind ! missing value
- real(rkind),parameter :: verySmall=1.e-12_rkind ! a very small number (used to check consistency)
- real(rkind),parameter :: veryBig=1.e+20_rkind ! a very big number
- real(rkind),parameter :: dx = 1.e-8_rkind ! finite difference increment
-
- contains
-
-
- ! **********************************************************************************************************
- ! public subroutine opSplittin: run the coupled energy-mass model for one timestep
- !
- ! The logic of the solver is as follows:
- ! (1) Attempt different solutions in the following order: (a) fully coupled; (b) split by state type and by
- ! domain type for a given energy and mass split (vegetation, snow, and soil); and (c) scalar solution
- ! for a given state type and domain subset.
- ! (2) For a given split, compute a variable number of substeps (in varSubstepSundials).
- ! **********************************************************************************************************
- subroutine opSplittin(&
- ! input: model control
- nSnow, & ! intent(in): number of snow layers
- nSoil, & ! intent(in): number of soil layers
- nLayers, & ! intent(in): total number of layers
- nState, & ! intent(in): total number of state variables
- dt, & ! intent(inout): time step (s)
- firstSubStep, & ! intent(in): flag to denote first sub-step
- computeVegFlux, & ! intent(in): flag to denote if computing energy flux over vegetation
- ! input/output: data structures
- type_data, & ! intent(in): type of vegetation and soil
- attr_data, & ! intent(in): spatial attributes
- forc_data, & ! intent(in): model forcing data
- mpar_data, & ! intent(in): model parameters
- indx_data, & ! intent(inout): index data
- prog_data, & ! intent(inout): model prognostic variables for a local HRU
- diag_data, & ! intent(inout): model diagnostic variables for a local HRU
- flux_data, & ! intent(inout): model fluxes for a local HRU
- bvar_data, & ! intent(in): model variables for the local basin
- lookup_data, & ! intent(in): lookup tables
- model_decisions,& ! intent(in): model decisions
- ! output: model control
- dtMultiplier, & ! intent(out): substep multiplier (-)
- tooMuchMelt, & ! intent(out): flag to denote that ice is insufficient to support melt
- stepFailure, & ! intent(out): flag to denote step failure
- ixCoupling, & ! intent(out): coupling method used in this iteration
- err,message) ! intent(out): error code and error message
- ! ---------------------------------------------------------------------------------------
- ! structure allocations
- USE allocspace_module,only:allocLocal ! allocate local data structures
- ! simulation of fluxes and residuals given a trial state vector
- USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric water content
- USE soil_utils_module,only:liquidHead ! compute the liquid water matric potential
- ! population/extraction of state vectors
- USE indexState_module,only:indexSplit ! get state indices
- USE varSubstep_module,only:varSubstep ! complete substeps for a given split
- USE varSubstepSundials_module,only:varSubstepSundials ! complete substeps for a given split
- ! identify name of variable type (for error message)
- USE get_ixName_module,only:get_varTypeName ! to access type strings for error messages
- implicit none
- ! ---------------------------------------------------------------------------------------
- ! * dummy variables
- ! ---------------------------------------------------------------------------------------
- ! input: model control
- integer(i4b),intent(in) :: nSnow ! number of snow layers
- integer(i4b),intent(in) :: nSoil ! number of soil layers
- integer(i4b),intent(in) :: nLayers ! total number of layers
- integer(i4b),intent(in) :: nState ! total number of state variables
- real(rkind),intent(inout) :: dt ! time step (seconds)
- logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
- 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(var_i),intent(in) :: type_data ! type of vegetation and soil
- type(var_d),intent(in) :: attr_data ! spatial attributes
- type(var_d),intent(in) :: forc_data ! model forcing data
- type(var_dlength),intent(in) :: mpar_data ! model parameters
- type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
- type(var_dlength),intent(inout) :: prog_data ! prognostic variables for a local HRU
- type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
- type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
- type(var_dlength),intent(in) :: bvar_data ! model variables for the local basin
- type(zLookup), intent(in) :: lookup_data ! lookup tables
- type(model_options),intent(in) :: model_decisions(:) ! model decisions
- ! output: model control
- real(rkind),intent(out) :: dtMultiplier ! substep multiplier (-)
- logical(lgt),intent(out) :: tooMuchMelt ! flag to denote that ice is insufficient to support melt
- logical(lgt),intent(out) :: stepFailure ! flag to denote step failure
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! ---------------------------------------------------------------------------------------
- ! * general local variables
- ! ---------------------------------------------------------------------------------------
- character(LEN=256) :: cmessage ! error message of downwind routine
- integer(i4b) :: minLayer ! the minimum layer used in assigning flags for flux aggregations
- integer(i4b) :: iOffset ! offset to account for different indices in the soil domain
- integer(i4b) :: iMin(1),iMax(1) ! bounds of a given vector
- integer(i4b) :: iLayer,jLayer ! index of model layer
- integer(i4b) :: iSoil ! index of soil layer
- integer(i4b) :: iVar ! index of variables in data structures
- logical(lgt) :: firstSuccess ! flag to define the first success
- logical(lgt) :: firstFluxCall ! flag to define the first flux call
- logical(lgt) :: reduceCoupledStep ! flag to define the need to reduce the length of the coupled step
- type(var_dlength) :: prog_temp ! temporary model prognostic variables
- type(var_dlength) :: diag_temp ! temporary model diagnostic variables
- type(var_dlength) :: flux_temp ! temporary model fluxes
- type(var_dlength) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
- real(rkind),dimension(nLayers) :: mLayerVolFracIceInit ! initial vector for volumetric fraction of ice (-)
- ! ------------------------------------------------------------------------------------------------------
- ! * operator splitting
- ! ------------------------------------------------------------------------------------------------------
- ! minimum timestep
- real(rkind),parameter :: dtmin_coupled=1800._rkind ! minimum time step for the fully coupled solution (seconds)
- real(rkind),parameter :: dtmin_split=60._rkind ! minimum time step for the fully split solution (seconds)
- real(rkind),parameter :: dtmin_scalar=10._rkind ! minimum time step for the scalar solution (seconds)
- real(rkind) :: dt_min ! minimum time step (seconds)
- real(rkind) :: dtInit ! initial time step (seconds)
- ! explicit error tolerance (depends on state type split, so defined here)
- real(rkind),parameter :: errorTolLiqFlux=0.01_rkind ! error tolerance in the explicit solution (liquid flux)
- real(rkind),parameter :: errorTolNrgFlux=10._rkind ! error tolerance in the explicit solution (energy flux)
- ! number of substeps taken for a given split
- integer(i4b) :: nSubsteps ! number of substeps taken for a given split
- ! named variables defining the coupling and solution method
- integer(i4b) :: ixCoupling ! index of coupling method (1,2)
- integer(i4b) :: ixSolution ! index of solution method (1,2)
- integer(i4b) :: ixStateThenDomain ! switch between the state and domain (1,2)
- integer(i4b) :: tryDomainSplit ! (0,1) - flag to try the domain split
- ! actual number of splits
- integer(i4b) :: nStateTypeSplit ! number of splits for the state type
- integer(i4b) :: nDomainSplit ! number of splits for the domain
- integer(i4b) :: nStateSplit ! number of splits for the states within a given domain
- ! indices for the state type and the domain split
- integer(i4b) :: iStateTypeSplit ! index of the state type split
- integer(i4b) :: iDomainSplit ! index of the domain split
- integer(i4b) :: iStateSplit ! index of the state split
- ! flux masks
- logical(lgt) :: neededFlux(nFlux) ! .true. if flux is needed at all
- logical(lgt) :: desiredFlux ! .true. if flux is desired for a given split
- type(var_ilength) :: fluxCount ! number of times each flux is updated (should equal nSubsteps)
- type(var_flagVec) :: fluxMask ! mask defining model fluxes
- ! state masks
- integer(i4b),dimension(nState) :: stateCheck ! number of times each state variable is updated (should equal 1)
- logical(lgt),dimension(nState) :: stateMask ! mask defining desired state variables
- integer(i4b) :: nSubset ! number of selected state variables for a given split
- ! flags
- logical(lgt) :: failure ! flag to denote failure of substepping
- logical(lgt) :: doAdjustTemp ! flag to adjust temperature after the mass split
- logical(lgt) :: failedMinimumStep ! flag to denote failure of substepping for a given split
- integer(i4b) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
- integer(i4b) :: nCoupling
- real(qp) :: dt_out !
- ! ---------------------------------------------------------------------------------------
- ! point to variables in the data structures
+ ! ---------------------------------------------------------------------------------------
+ ! * dummy variables
+ ! ---------------------------------------------------------------------------------------
+ ! input: model control
+ integer(i4b),intent(in) :: nSnow ! number of snow layers
+ integer(i4b),intent(in) :: nSoil ! number of soil layers
+ integer(i4b),intent(in) :: nLayers ! total number of layers
+ integer(i4b),intent(in) :: nState ! total number of state variables
+ real(rkind),intent(inout) :: dt ! time step (seconds)
+ logical(lgt),intent(in) :: firstSubStep ! flag to indicate if we are processing the first sub-step
+ 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(var_i),intent(in) :: type_data ! type of vegetation and soil
+ type(var_d),intent(in) :: attr_data ! spatial attributes
+ type(var_d),intent(in) :: forc_data ! model forcing data
+ type(var_dlength),intent(in) :: mpar_data ! model parameters
+ type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
+ type(var_dlength),intent(inout) :: prog_data ! prognostic variables for a local HRU
+ type(var_dlength),intent(inout) :: diag_data ! diagnostic variables for a local HRU
+ type(var_dlength),intent(inout) :: flux_data ! model fluxes for a local HRU
+ type(var_dlength),intent(in) :: bvar_data ! model variables for the local basin
+ type(zLookup), intent(in) :: lookup_data ! lookup tables
+ type(model_options),intent(in) :: model_decisions(:) ! model decisions
+ ! output: model control
+ real(rkind),intent(out) :: dtMultiplier ! substep multiplier (-)
+ logical(lgt),intent(out) :: tooMuchMelt ! flag to denote that ice is insufficient to support melt
+ logical(lgt),intent(out) :: stepFailure ! flag to denote step failure
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! ---------------------------------------------------------------------------------------
+ ! * general local variables
+ ! ---------------------------------------------------------------------------------------
+ character(LEN=256) :: cmessage ! error message of downwind routine
+ integer(i4b) :: minLayer ! the minimum layer used in assigning flags for flux aggregations
+ integer(i4b) :: iOffset ! offset to account for different indices in the soil domain
+ integer(i4b) :: iMin(1),iMax(1) ! bounds of a given vector
+ integer(i4b) :: iLayer,jLayer ! index of model layer
+ integer(i4b) :: iSoil ! index of soil layer
+ integer(i4b) :: iVar ! index of variables in data structures
+ logical(lgt) :: firstSuccess ! flag to define the first success
+ logical(lgt) :: firstFluxCall ! flag to define the first flux call
+ logical(lgt) :: reduceCoupledStep ! flag to define the need to reduce the length of the coupled step
+ type(var_dlength) :: prog_temp ! temporary model prognostic variables
+ type(var_dlength) :: diag_temp ! temporary model diagnostic variables
+ type(var_dlength) :: flux_temp ! temporary model fluxes
+ type(var_dlength) :: deriv_data ! derivatives in model fluxes w.r.t. relevant state variables
+ real(rkind),dimension(nLayers) :: mLayerVolFracIceInit ! initial vector for volumetric fraction of ice (-)
+ ! ------------------------------------------------------------------------------------------------------
+ ! * operator splitting
+ ! ------------------------------------------------------------------------------------------------------
+ ! minimum timestep
+ real(rkind),parameter :: dtmin_coupled=1800._rkind ! minimum time step for the fully coupled solution (seconds)
+ real(rkind),parameter :: dtmin_split=60._rkind ! minimum time step for the fully split solution (seconds)
+ real(rkind),parameter :: dtmin_scalar=10._rkind ! minimum time step for the scalar solution (seconds)
+ real(rkind) :: dt_min ! minimum time step (seconds)
+ real(rkind) :: dtInit ! initial time step (seconds)
+ ! explicit error tolerance (depends on state type split, so defined here)
+ real(rkind),parameter :: errorTolLiqFlux=0.01_rkind ! error tolerance in the explicit solution (liquid flux)
+ real(rkind),parameter :: errorTolNrgFlux=10._rkind ! error tolerance in the explicit solution (energy flux)
+ ! number of substeps taken for a given split
+ integer(i4b) :: nSubsteps ! number of substeps taken for a given split
+ ! named variables defining the coupling and solution method
+ integer(i4b) :: ixCoupling ! index of coupling method (1,2)
+ integer(i4b) :: ixSolution ! index of solution method (1,2)
+ integer(i4b) :: ixStateThenDomain ! switch between the state and domain (1,2)
+ integer(i4b) :: tryDomainSplit ! (0,1) - flag to try the domain split
+ ! actual number of splits
+ integer(i4b) :: nStateTypeSplit ! number of splits for the state type
+ integer(i4b) :: nDomainSplit ! number of splits for the domain
+ integer(i4b) :: nStateSplit ! number of splits for the states within a given domain
+ ! indices for the state type and the domain split
+ integer(i4b) :: iStateTypeSplit ! index of the state type split
+ integer(i4b) :: iDomainSplit ! index of the domain split
+ integer(i4b) :: iStateSplit ! index of the state split
+ ! flux masks
+ logical(lgt) :: neededFlux(nFlux) ! .true. if flux is needed at all
+ logical(lgt) :: desiredFlux ! .true. if flux is desired for a given split
+ type(var_ilength) :: fluxCount ! number of times each flux is updated (should equal nSubsteps)
+ type(var_flagVec) :: fluxMask ! mask defining model fluxes
+ ! state masks
+ integer(i4b),dimension(nState) :: stateCheck ! number of times each state variable is updated (should equal 1)
+ logical(lgt),dimension(nState) :: stateMask ! mask defining desired state variables
+ integer(i4b) :: nSubset ! number of selected state variables for a given split
+ ! flags
+ logical(lgt) :: failure ! flag to denote failure of substepping
+ logical(lgt) :: doAdjustTemp ! flag to adjust temperature after the mass split
+ logical(lgt) :: failedMinimumStep ! flag to denote failure of substepping for a given split
+ integer(i4b) :: ixSaturation ! index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ integer(i4b) :: nCoupling
+ real(qp) :: dt_out !
+ ! ---------------------------------------------------------------------------------------
+ ! point to variables in the data structures
+ ! ---------------------------------------------------------------------------------------
+
+ globalVars: associate(&
+ ! model decisions
+ ixGroundwater => model_decisions(iLookDECISIONS%groundwatr)%iDecision ,& ! intent(in): [i4b] groundwater parameterization
+ ixSpatialGroundwater => model_decisions(iLookDECISIONS%spatial_gw)%iDecision ,& ! intent(in): [i4b] spatial representation of groundwater (local-column or single-basin)
+ ixNumericalMethod => model_decisions(iLookDECISIONS%num_method)%iDecision ,& ! intent(in): [i4b] choice of numerical method, backward Euler or SUNDIALS/IDA
+ ! domain boundary conditions
+ airtemp => forc_data%var(iLookFORCE%airtemp) ,& ! intent(in): [dp] temperature of the upper boundary of the snow and soil domains (K)
+ ! 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
+ ! indices of model state variables
+ ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset (missing for values not in the subset)
+ ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
+ ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
+ ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
+ ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
+ ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
+ ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
+ ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
+ ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
+ ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
+ ! numerix tracking
+ numberStateSplit => indx_data%var(iLookINDEX%numberStateSplit )%dat(1) ,& ! intent(inout): [i4b] number of state splitting solutions (-)
+ numberDomainSplitNrg => indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for energy (-)
+ numberDomainSplitMass => indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for mass (-)
+ numberScalarSolutions => indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) ,& ! intent(inout): [i4b] number of scalar solutions (-)
+ ! domain configuration
+ canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
+ mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
+ ! snow parameters
+ snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
+ ! depth-varying soil parameters
+ vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat ,& ! intent(in): [dp(:)] van Genutchen "m" parameter (-)
+ vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat ,& ! intent(in): [dp(:)] van Genutchen "n" parameter (-)
+ vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat ,& ! intent(in): [dp(:)] van Genutchen "alpha" parameter (m-1)
+ theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-)
+ theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] soil residual volumetric water content (-)
+ ! soil parameters
+ specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
+ ! model diagnostic variables (fraction of liquid water)
+ scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(out): [dp] fraction of liquid water on vegetation (-)
+ mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(out): [dp(:)] fraction of liquid water in each snow layer (-)
+ mLayerMeltFreeze => diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat ,& ! intent(out): [dp(:)] melt-freeze in each snow and soil layer (kg m-3)
+ ! model state variables (vegetation canopy)
+ scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(out): [dp] temperature of the canopy air space (K)
+ scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(out): [dp] temperature of the vegetation canopy (K)
+ scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(out): [dp] mass of ice on the vegetation canopy (kg m-2)
+ scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(out): [dp] mass of liquid water on the vegetation canopy (kg m-2)
+ scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(out): [dp] mass of total water on the vegetation canopy (kg m-2)
+ ! model state variables (snow and soil domains)
+ mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(out): [dp(:)] temperature of each snow/soil layer (K)
+ mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(out): [dp(:)] volumetric fraction of ice (-)
+ mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(out): [dp(:)] volumetric fraction of liquid water (-)
+ mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(out): [dp(:)] volumetric fraction of total water (-)
+ mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(out): [dp(:)] matric head (m)
+ mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat & ! intent(out): [dp(:)] matric potential of liquid water (m)
+ )
! ---------------------------------------------------------------------------------------
-
- globalVars: associate(&
- ! model decisions
- ixGroundwater => model_decisions(iLookDECISIONS%groundwatr)%iDecision ,& ! intent(in): [i4b] groundwater parameterization
- ixSpatialGroundwater => model_decisions(iLookDECISIONS%spatial_gw)%iDecision ,& ! intent(in): [i4b] spatial representation of groundwater (local-column or single-basin)
- ixNumericalMethod => model_decisions(iLookDECISIONS%num_method)%iDecision ,& ! intent(in): [i4b] choice of numerical method, backward Euler or SUNDIALS/IDA
- ! domain boundary conditions
- airtemp => forc_data%var(iLookFORCE%airtemp) ,& ! intent(in): [dp] temperature of the upper boundary of the snow and soil domains (K)
- ! 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
- ! indices of model state variables
- ixMapSubset2Full => indx_data%var(iLookINDEX%ixMapSubset2Full)%dat ,& ! intent(in): [i4b(:)] list of indices in the state subset (missing for values not in the subset)
- ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
- ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
- ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
- ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
- ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
- ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
- ixCasNrg => indx_data%var(iLookINDEX%ixCasNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy air space energy state variable
- ixVegNrg => indx_data%var(iLookINDEX%ixVegNrg)%dat(1) ,& ! intent(in): [i4b] index of canopy energy state variable
- ixVegHyd => indx_data%var(iLookINDEX%ixVegHyd)%dat(1) ,& ! intent(in): [i4b] index of canopy hydrology state variable (mass)
- ! numerix tracking
- numberStateSplit => indx_data%var(iLookINDEX%numberStateSplit )%dat(1) ,& ! intent(inout): [i4b] number of state splitting solutions (-)
- numberDomainSplitNrg => indx_data%var(iLookINDEX%numberDomainSplitNrg )%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for energy (-)
- numberDomainSplitMass => indx_data%var(iLookINDEX%numberDomainSplitMass)%dat(1) ,& ! intent(inout): [i4b] number of domain splitting solutions for mass (-)
- numberScalarSolutions => indx_data%var(iLookINDEX%numberScalarSolutions)%dat(1) ,& ! intent(inout): [i4b] number of scalar solutions (-)
- ! domain configuration
- canopyDepth => diag_data%var(iLookDIAG%scalarCanopyDepth)%dat(1) ,& ! intent(in): [dp] canopy depth (m)
- mLayerDepth => prog_data%var(iLookPROG%mLayerDepth)%dat ,& ! intent(in): [dp(:)] depth of each layer in the snow-soil sub-domain (m)
- ! snow parameters
- snowfrz_scale => mpar_data%var(iLookPARAM%snowfrz_scale)%dat(1) ,& ! intent(in): [dp] scaling parameter for the snow freezing curve (K-1)
- ! depth-varying soil parameters
- vGn_m => diag_data%var(iLookDIAG%scalarVGn_m)%dat ,& ! intent(in): [dp(:)] van Genutchen "m" parameter (-)
- vGn_n => mpar_data%var(iLookPARAM%vGn_n)%dat ,& ! intent(in): [dp(:)] van Genutchen "n" parameter (-)
- vGn_alpha => mpar_data%var(iLookPARAM%vGn_alpha)%dat ,& ! intent(in): [dp(:)] van Genutchen "alpha" parameter (m-1)
- theta_sat => mpar_data%var(iLookPARAM%theta_sat)%dat ,& ! intent(in): [dp(:)] soil porosity (-)
- theta_res => mpar_data%var(iLookPARAM%theta_res)%dat ,& ! intent(in): [dp(:)] soil residual volumetric water content (-)
- ! soil parameters
- specificStorage => mpar_data%var(iLookPARAM%specificStorage)%dat(1) ,& ! intent(in): [dp] specific storage coefficient (m-1)
- ! model diagnostic variables (fraction of liquid water)
- scalarFracLiqVeg => diag_data%var(iLookDIAG%scalarFracLiqVeg)%dat(1) ,& ! intent(out): [dp] fraction of liquid water on vegetation (-)
- mLayerFracLiqSnow => diag_data%var(iLookDIAG%mLayerFracLiqSnow)%dat ,& ! intent(out): [dp(:)] fraction of liquid water in each snow layer (-)
- mLayerMeltFreeze => diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat ,& ! intent(out): [dp(:)] melt-freeze in each snow and soil layer (kg m-3)
- ! model state variables (vegetation canopy)
- scalarCanairTemp => prog_data%var(iLookPROG%scalarCanairTemp)%dat(1) ,& ! intent(out): [dp] temperature of the canopy air space (K)
- scalarCanopyTemp => prog_data%var(iLookPROG%scalarCanopyTemp)%dat(1) ,& ! intent(out): [dp] temperature of the vegetation canopy (K)
- scalarCanopyIce => prog_data%var(iLookPROG%scalarCanopyIce)%dat(1) ,& ! intent(out): [dp] mass of ice on the vegetation canopy (kg m-2)
- scalarCanopyLiq => prog_data%var(iLookPROG%scalarCanopyLiq)%dat(1) ,& ! intent(out): [dp] mass of liquid water on the vegetation canopy (kg m-2)
- scalarCanopyWat => prog_data%var(iLookPROG%scalarCanopyWat)%dat(1) ,& ! intent(out): [dp] mass of total water on the vegetation canopy (kg m-2)
- ! model state variables (snow and soil domains)
- mLayerTemp => prog_data%var(iLookPROG%mLayerTemp)%dat ,& ! intent(out): [dp(:)] temperature of each snow/soil layer (K)
- mLayerVolFracIce => prog_data%var(iLookPROG%mLayerVolFracIce)%dat ,& ! intent(out): [dp(:)] volumetric fraction of ice (-)
- mLayerVolFracLiq => prog_data%var(iLookPROG%mLayerVolFracLiq)%dat ,& ! intent(out): [dp(:)] volumetric fraction of liquid water (-)
- mLayerVolFracWat => prog_data%var(iLookPROG%mLayerVolFracWat)%dat ,& ! intent(out): [dp(:)] volumetric fraction of total water (-)
- mLayerMatricHead => prog_data%var(iLookPROG%mLayerMatricHead)%dat ,& ! intent(out): [dp(:)] matric head (m)
- mLayerMatricHeadLiq => diag_data%var(iLookDIAG%mLayerMatricHeadLiq)%dat & ! intent(out): [dp(:)] matric potential of liquid water (m)
- )
- ! ---------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message="opSplittin/"
-
- ! we just solve the fully coupled problem by with Sundials for now
- select case(ixNumericalMethod)
- case(sundials); nCoupling = 1
- case(bEuler); nCoupling = 2
- case default; err=20; message=trim(message)//'expect num_method to be sundials or bEuler (or itertive, which is bEuler)'; return
- end select
-
- ! -----
- ! * initialize...
- ! ---------------
-
- ! set the global print flag
- globalPrintFlag=.false.
-
- if(globalPrintFlag)&
- print*, trim(message), dt
-
- ! initialize the first success call
- firstSuccess=.false.
-
- ! initialize the flags
- tooMuchMelt=.false. ! too much melt (merge snow layers)
- stepFailure=.false. ! step failure
-
- ! initialize flag for the success of the substepping
- failure=.false.
-
- ! initialize the flux check
- neededFlux(:) = .false.
-
- ! initialize the state check
- stateCheck(:) = 0
-
- ! compute the total water content in the vegetation canopy
- scalarCanopyWat = scalarCanopyLiq + scalarCanopyIce ! kg m-2
-
- ! save volumetric ice content at the start of the step
- ! NOTE: used for volumetric loss due to melt-freeze
- mLayerVolFracIceInit(:) = mLayerVolFracIce(:)
-
- ! compute the total water content in snow and soil
- ! NOTE: no ice expansion allowed for soil
- if(nSnow>0)&
- mLayerVolFracWat( 1:nSnow ) = mLayerVolFracLiq( 1:nSnow ) + mLayerVolFracIce( 1:nSnow )*(iden_ice/iden_water)
- mLayerVolFracWat(nSnow+1:nLayers) = mLayerVolFracLiq(nSnow+1:nLayers) + mLayerVolFracIce(nSnow+1:nLayers)
-
- ! compute the liquid water matric potential (m)
- ! NOTE: include ice content as part of the solid porosity - major effect of ice is to reduce the pore size; ensure that effSat=1 at saturation
- ! (from Zhao et al., J. Hydrol., 1997: Numerical analysis of simultaneous heat and mass transfer...)
- do iSoil=1,nSoil
- call liquidHead(mLayerMatricHead(iSoil),mLayerVolFracLiq(nSnow+iSoil),mLayerVolFracIce(nSnow+iSoil), & ! input: state variables
- vGn_alpha(iSoil),vGn_n(iSoil),theta_sat(iSoil),theta_res(iSoil),vGn_m(iSoil), & ! input: parameters
- matricHeadLiq=mLayerMatricHeadLiq(iSoil), & ! output: liquid water matric potential (m)
- err=err,message=cmessage) ! output: error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
- end do ! looping through soil layers (computing liquid water matric potential)
-
- ! allocate space for the flux mask (used to define when fluxes are updated)
- call allocLocal(flux_meta(:),fluxMask,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! allocate space for the flux count (used to check that fluxes are only updated once)
- call allocLocal(flux_meta(:),fluxCount,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! allocate space for the temporary prognostic variable structure
- call allocLocal(prog_meta(:),prog_temp,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! allocate space for the temporary diagnostic variable structure
- call allocLocal(diag_meta(:),diag_temp,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! allocate space for the temporary flux variable structure
- call allocLocal(flux_meta(:),flux_temp,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
-
- ! allocate space for the derivative structure
- call allocLocal(deriv_meta(:),deriv_data,nSnow,nSoil,err,cmessage)
- if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
-
- ! intialize the flux conter
- do iVar=1,size(flux_meta) ! loop through fluxes
- fluxCount%var(iVar)%dat(:) = 0
- end do
-
- ! initialize the model fluxes
- do iVar=1,size(flux_meta) ! loop through fluxes
- if(flux2state_orig(iVar)%state1==integerMissing .and. flux2state_orig(iVar)%state2==integerMissing) cycle ! flux does not depend on state (e.g., input)
- if(flux2state_orig(iVar)%state1==iname_watCanopy .and. .not.computeVegFlux) cycle ! use input fluxes in cases where there is no canopy
- flux_data%var(iVar)%dat(:) = 0._rkind
- end do
-
- ! initialize derivatives
- do iVar=1,size(deriv_meta)
- deriv_data%var(iVar)%dat(:) = 0._rkind
- end do
-
- ! ==========================================================================================================================================
- ! ==========================================================================================================================================
- ! ==========================================================================================================================================
- ! ==========================================================================================================================================
-
- ! loop through different coupling strategies
- coupling: do ixCoupling=1,nCoupling
-
- ! initialize the time step
- dtInit = min( merge(dt, dtmin_coupled, ixCoupling==fullyCoupled), dt) ! initial time step
- dt_min = min( merge(dtmin_coupled, dtmin_split, ixCoupling==fullyCoupled), dt) ! minimum time step
-
- ! keep track of the number of state splits
- if(ixCoupling/=fullyCoupled) numberStateSplit = numberStateSplit + 1
-
- ! define the number of operator splits for the state type
- select case(ixCoupling)
- case(fullyCoupled); nStateTypeSplit=1
- case(stateTypeSplit); nStateTypeSplit=nStateTypes
- case default; err=20; message=trim(message)//'coupling case not found'; return
- end select ! operator splitting option
-
- ! define if we wish to try the domain split
- select case(ixCoupling)
- case(fullyCoupled); tryDomainSplit=0
- case(stateTypeSplit); tryDomainSplit=1
- case default; err=20; message=trim(message)//'coupling case not found'; return
- end select ! operator splitting option
-
- ! state splitting loop
- stateTypeSplit: do iStateTypeSplit=1,nStateTypeSplit
-
- ! -----
- ! * identify state-specific variables for a given state split...
- ! --------------------------------------------------------------
-
- ! flag to adjust the temperature
- doAdjustTemp = (ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)
-
- ! modify the state type names associated with the state vector
- if(ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)then
- if(computeVegFlux)then
- where(ixStateType(ixHydCanopy)==iname_watCanopy) ixStateType(ixHydCanopy)=iname_liqCanopy
- endif
- where(ixStateType(ixHydLayer) ==iname_watLayer) ixStateType(ixHydLayer) =iname_liqLayer
- where(ixStateType(ixHydLayer) ==iname_matLayer) ixStateType(ixHydLayer) =iname_lmpLayer
- endif ! if modifying state variables for the mass split
-
- ! first try the state type split, then try the domain split within a given state type
- stateThenDomain: do ixStateThenDomain=1,1+tryDomainSplit ! 1=state type split; 2=domain split within a given state type
-
- ! keep track of the number of domain splits
- if(iStateTypeSplit==nrgSplit .and. ixStateThenDomain==subDomain) numberDomainSplitNrg = numberDomainSplitNrg + 1
- if(iStateTypeSplit==massSplit .and. ixStateThenDomain==subDomain) numberDomainSplitMass = numberDomainSplitMass + 1
-
- ! define the number of domain splits for the state type
- select case(ixStateThenDomain)
- case(fullDomain); nDomainSplit=1
- case(subDomain); nDomainSplit=nDomains
- case default; err=20; message=trim(message)//'coupling case not found'; return
- end select
-
- ! check that we haven't split the domain when we are fully coupled
- if(ixCoupling==fullyCoupled .and. nDomainSplit==nDomains)then
- message=trim(message)//'cannot split domains when fully coupled'
- err=20; return
- endif
-
- ! domain splitting loop
- domainSplit: do iDomainSplit=1,nDomainSplit
-
- ! trial with the vector then scalar solution
- solution: do ixSolution=1,nSolutions
-
- ! initialize error control
- err=0; message="opSplittin/"
-
- ! refine the time step
- if(ixSolution==scalar)then
- dtInit = min(dtmin_split, dt) ! initial time step
- dt_min = min(dtmin_scalar, dt) ! minimum time step
- endif
-
- ! initialize the first flux call
- firstFluxCall=.true.
-
- ! get the number of split layers
- select case(ixSolution)
- case(vector); nStateSplit=1
- case(scalar); nStateSplit=count(stateMask)
- case default; err=20; message=trim(message)//'unknown solution method'; return
- end select
-
-
- ! loop through layers (NOTE: nStateSplit=1 for the vector solution, hence no looping)
- stateSplit: do iStateSplit=1,nStateSplit
-
- ! -----
- ! * define state subsets for a given split...
- ! -------------------------------------------
-
- ! get the mask for the state subset
- call stateFilter(ixCoupling,ixSolution,ixStateThenDomain,iStateTypeSplit,iDomainSplit,iStateSplit,&
- indx_data,stateMask,nSubset,err,cmessage)
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif ! (check for errors)
-
- ! check that state variables exist
- if(nSubset==0) cycle domainSplit
-
- ! avoid redundant case where vector solution is of length 1
- if(ixSolution==vector .and. count(stateMask)==1) cycle solution
-
-
-
- ! -----
- ! * assemble vectors for a given split...
- ! ---------------------------------------
- ! get indices for a given split
- call indexSplit(stateMask, & ! intent(in) : logical vector (.true. if state is in the subset)
- nSnow,nSoil,nLayers,nSubset, & ! intent(in) : number of snow and soil layers, and total number of layers
- indx_data, & ! intent(inout) : index data structure
- err,cmessage) ! intent(out) : error control
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
- ! -----
- ! * define the mask of the fluxes used...
- ! ---------------------------------------
-
- ! identify the type of state for the states in the subset
- stateSubset: associate(ixStateType_subset => indx_data%var(iLookINDEX%ixStateType_subset)%dat ,& ! intent(in): [i4b(:)] indices of state types
- ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset
- ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil)
- ixLayerActive => indx_data%var(iLookINDEX%ixLayerActive)%dat ,& ! intent(in): [i4b(:)] list of indices for all active layers (inactive=integerM
- ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ) ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil)
-
- ! loop through flux variables
- do iVar=1,size(flux_meta)
-
- ! * identify flux mask for the fully coupled solution
- if(ixCoupling==fullyCoupled)then
- desiredFlux = any(ixStateType_subset==flux2state_orig(iVar)%state1) .or. any(ixStateType_subset==flux2state_orig(iVar)%state2)
+ ! initialize error control
+ err=0; message="opSplittin/"
+
+ ! we just solve the fully coupled problem by with Sundials for now
+ select case(ixNumericalMethod)
+ case(sundials); nCoupling = 1
+ case(bEuler); nCoupling = 2
+ case default; err=20; message=trim(message)//'expect num_method to be sundials or bEuler (or itertive, which is bEuler)'; return
+ end select
+
+ ! -----
+ ! * initialize...
+ ! ---------------
+
+ ! set the global print flag
+ globalPrintFlag=.false.
+
+ ! if(globalPrintFlag)&
+ print*, trim(message), dt
+
+ ! initialize the first success call
+ firstSuccess=.false.
+
+ ! initialize the flags
+ tooMuchMelt=.false. ! too much melt (merge snow layers)
+ stepFailure=.false. ! step failure
+
+ ! initialize flag for the success of the substepping
+ failure=.false.
+
+ ! initialize the flux check
+ neededFlux(:) = .false.
+
+ ! initialize the state check
+ stateCheck(:) = 0
+
+ ! compute the total water content in the vegetation canopy
+ scalarCanopyWat = scalarCanopyLiq + scalarCanopyIce ! kg m-2
+
+ ! save volumetric ice content at the start of the step
+ ! NOTE: used for volumetric loss due to melt-freeze
+ mLayerVolFracIceInit(:) = mLayerVolFracIce(:)
+
+ ! compute the total water content in snow and soil
+ ! NOTE: no ice expansion allowed for soil
+ if(nSnow>0)&
+ mLayerVolFracWat( 1:nSnow ) = mLayerVolFracLiq( 1:nSnow ) + mLayerVolFracIce( 1:nSnow )*(iden_ice/iden_water)
+ mLayerVolFracWat(nSnow+1:nLayers) = mLayerVolFracLiq(nSnow+1:nLayers) + mLayerVolFracIce(nSnow+1:nLayers)
+
+ ! compute the liquid water matric potential (m)
+ ! NOTE: include ice content as part of the solid porosity - major effect of ice is to reduce the pore size; ensure that effSat=1 at saturation
+ ! (from Zhao et al., J. Hydrol., 1997: Numerical analysis of simultaneous heat and mass transfer...)
+ do iSoil=1,nSoil
+ call liquidHead(mLayerMatricHead(iSoil),mLayerVolFracLiq(nSnow+iSoil),mLayerVolFracIce(nSnow+iSoil), & ! input: state variables
+ vGn_alpha(iSoil),vGn_n(iSoil),theta_sat(iSoil),theta_res(iSoil),vGn_m(iSoil), & ! input: parameters
+ matricHeadLiq=mLayerMatricHeadLiq(iSoil), & ! output: liquid water matric potential (m)
+ err=err,message=cmessage) ! output: error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+ end do ! looping through soil layers (computing liquid water matric potential)
+
+ ! allocate space for the flux mask (used to define when fluxes are updated)
+ call allocLocal(flux_meta(:),fluxMask,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+
+ ! allocate space for the flux count (used to check that fluxes are only updated once)
+ call allocLocal(flux_meta(:),fluxCount,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+
+ ! allocate space for the temporary prognostic variable structure
+ call allocLocal(prog_meta(:),prog_temp,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+
+ ! allocate space for the temporary diagnostic variable structure
+ call allocLocal(diag_meta(:),diag_temp,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+
+ ! allocate space for the temporary flux variable structure
+ call allocLocal(flux_meta(:),flux_temp,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; endif
+
+ ! allocate space for the derivative structure
+ call allocLocal(deriv_meta(:),deriv_data,nSnow,nSoil,err,cmessage)
+ if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
+
+ ! intialize the flux conter
+ do iVar=1,size(flux_meta) ! loop through fluxes
+ fluxCount%var(iVar)%dat(:) = 0
+ end do
+
+ ! initialize the model fluxes
+ do iVar=1,size(flux_meta) ! loop through fluxes
+ if(flux2state_orig(iVar)%state1==integerMissing .and. flux2state_orig(iVar)%state2==integerMissing) cycle ! flux does not depend on state (e.g., input)
+ if(flux2state_orig(iVar)%state1==iname_watCanopy .and. .not.computeVegFlux) cycle ! use input fluxes in cases where there is no canopy
+ flux_data%var(iVar)%dat(:) = 0._rkind
+ end do
+
+ ! initialize derivatives
+ do iVar=1,size(deriv_meta)
+ deriv_data%var(iVar)%dat(:) = 0._rkind
+ end do
+
+ ! ==========================================================================================================================================
+ ! ==========================================================================================================================================
+ ! ==========================================================================================================================================
+ ! ==========================================================================================================================================
+
+ ! loop through different coupling strategies
+ coupling: do ixCoupling=1,nCoupling
+
+ ! initialize the time step
+ dtInit = min( merge(dt, dtmin_coupled, ixCoupling==fullyCoupled), dt) ! initial time step
+ dt_min = min( merge(dtmin_coupled, dtmin_split, ixCoupling==fullyCoupled), dt) ! minimum time step
+
+ ! keep track of the number of state splits
+ if(ixCoupling/=fullyCoupled) numberStateSplit = numberStateSplit + 1
+
+ ! define the number of operator splits for the state type
+ select case(ixCoupling)
+ case(fullyCoupled); nStateTypeSplit=1
+ case(stateTypeSplit); nStateTypeSplit=nStateTypes
+ case default; err=20; message=trim(message)//'coupling case not found'; return
+ end select ! operator splitting option
+
+ ! define if we wish to try the domain split
+ select case(ixCoupling)
+ case(fullyCoupled); tryDomainSplit=0
+ case(stateTypeSplit); tryDomainSplit=1
+ case default; err=20; message=trim(message)//'coupling case not found'; return
+ end select ! operator splitting option
+
+ ! state splitting loop
+ stateTypeSplit: do iStateTypeSplit=1,nStateTypeSplit
+
+ ! -----
+ ! * identify state-specific variables for a given state split...
+ ! --------------------------------------------------------------
+
+ ! flag to adjust the temperature
+ doAdjustTemp = (ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)
+
+ ! modify the state type names associated with the state vector
+ if(ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)then
+ if(computeVegFlux)then
+ where(ixStateType(ixHydCanopy)==iname_watCanopy) ixStateType(ixHydCanopy)=iname_liqCanopy
+ endif
+ where(ixStateType(ixHydLayer) ==iname_watLayer) ixStateType(ixHydLayer) =iname_liqLayer
+ where(ixStateType(ixHydLayer) ==iname_matLayer) ixStateType(ixHydLayer) =iname_lmpLayer
+ endif ! if modifying state variables for the mass split
+
+ ! first try the state type split, then try the domain split within a given state type
+ stateThenDomain: do ixStateThenDomain=1,1+tryDomainSplit ! 1=state type split; 2=domain split within a given state type
+
+ ! keep track of the number of domain splits
+ if(iStateTypeSplit==nrgSplit .and. ixStateThenDomain==subDomain) numberDomainSplitNrg = numberDomainSplitNrg + 1
+ if(iStateTypeSplit==massSplit .and. ixStateThenDomain==subDomain) numberDomainSplitMass = numberDomainSplitMass + 1
+
+ ! define the number of domain splits for the state type
+ select case(ixStateThenDomain)
+ case(fullDomain); nDomainSplit=1
+ case(subDomain); nDomainSplit=nDomains
+ case default; err=20; message=trim(message)//'coupling case not found'; return
+ end select
+
+ ! check that we haven't split the domain when we are fully coupled
+ if(ixCoupling==fullyCoupled .and. nDomainSplit==nDomains)then
+ message=trim(message)//'cannot split domains when fully coupled'
+ err=20; return
+ endif
+
+ ! domain splitting loop
+ domainSplit: do iDomainSplit=1,nDomainSplit
+
+ ! trial with the vector then scalar solution
+ solution: do ixSolution=1,nSolutions
+
+ ! initialize error control
+ err=0; message="opSplittin/"
+
+ ! refine the time step
+ if(ixSolution==scalar)then
+ dtInit = min(dtmin_split, dt) ! initial time step
+ dt_min = min(dtmin_scalar, dt) ! minimum time step
+ endif
+
+ ! initialize the first flux call
+ firstFluxCall=.true.
+
+ ! get the number of split layers
+ select case(ixSolution)
+ case(vector); nStateSplit=1
+ case(scalar); nStateSplit=count(stateMask)
+ case default; err=20; message=trim(message)//'unknown solution method'; return
+ end select
+
+
+ ! loop through layers (NOTE: nStateSplit=1 for the vector solution, hence no looping)
+ stateSplit: do iStateSplit=1,nStateSplit
+
+ ! -----
+ ! * define state subsets for a given split...
+ ! -------------------------------------------
+
+ ! get the mask for the state subset
+ call stateFilter(ixCoupling,ixSolution,ixStateThenDomain,iStateTypeSplit,iDomainSplit,iStateSplit,&
+ indx_data,stateMask,nSubset,err,cmessage)
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif ! (check for errors)
+
+ ! check that state variables exist
+ if(nSubset==0) cycle domainSplit
+
+ ! avoid redundant case where vector solution is of length 1
+ if(ixSolution==vector .and. count(stateMask)==1) cycle solution
+
+
+
+ ! -----
+ ! * assemble vectors for a given split...
+ ! ---------------------------------------
+ ! get indices for a given split
+ call indexSplit(stateMask, & ! intent(in) : logical vector (.true. if state is in the subset)
+ nSnow,nSoil,nLayers,nSubset, & ! intent(in) : number of snow and soil layers, and total number of layers
+ indx_data, & ! intent(inout) : index data structure
+ err,cmessage) ! intent(out) : error control
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+ ! -----
+ ! * define the mask of the fluxes used...
+ ! ---------------------------------------
+
+ ! identify the type of state for the states in the subset
+ stateSubset: associate(ixStateType_subset => indx_data%var(iLookINDEX%ixStateType_subset)%dat ,& ! intent(in): [i4b(:)] indices of state types
+ ixMapFull2Subset => indx_data%var(iLookINDEX%ixMapFull2Subset)%dat ,& ! intent(in): [i4b(:)] mapping of full state vector to the state subset
+ ixControlVolume => indx_data%var(iLookINDEX%ixControlVolume)%dat ,& ! intent(in): [i4b(:)] index of control volume for different domains (veg, snow, soil)
+ ixLayerActive => indx_data%var(iLookINDEX%ixLayerActive)%dat ,& ! intent(in): [i4b(:)] list of indices for all active layers (inactive=integerM
+ ixDomainType => indx_data%var(iLookINDEX%ixDomainType)%dat ) ! intent(in): [i4b(:)] indices defining the type of the domain (iname_veg, iname_snow, iname_soil)
+
+ ! loop through flux variables
+ do iVar=1,size(flux_meta)
+
+ ! * identify flux mask for the fully coupled solution
+ if(ixCoupling==fullyCoupled)then
+ desiredFlux = any(ixStateType_subset==flux2state_orig(iVar)%state1) .or. any(ixStateType_subset==flux2state_orig(iVar)%state2)
+ fluxMask%var(iVar)%dat = desiredFlux
+
+ ! * identify flux mask for the split solution
+ else
+
+ ! identify the flux mask for a given state split
+ select case(iStateTypeSplit)
+ case(nrgSplit); desiredFlux = any(ixStateType_subset==flux2state_orig(iVar)%state1) .or. any(ixStateType_subset==flux2state_orig(iVar)%state2)
+ case(massSplit); desiredFlux = any(ixStateType_subset==flux2state_liq(iVar)%state1) .or. any(ixStateType_subset==flux2state_liq(iVar)%state2)
+ case default; err=20; message=trim(message)//'unable to identify split based on state type'; return
+ end select
+
+ ! no domain splitting
+ if(nDomains==1)then
fluxMask%var(iVar)%dat = desiredFlux
-
- ! * identify flux mask for the split solution
+
+ ! domain splitting
else
-
- ! identify the flux mask for a given state split
- select case(iStateTypeSplit)
- case(nrgSplit); desiredFlux = any(ixStateType_subset==flux2state_orig(iVar)%state1) .or. any(ixStateType_subset==flux2state_orig(iVar)%state2)
- case(massSplit); desiredFlux = any(ixStateType_subset==flux2state_liq(iVar)%state1) .or. any(ixStateType_subset==flux2state_liq(iVar)%state2)
- case default; err=20; message=trim(message)//'unable to identify split based on state type'; return
- end select
-
- ! no domain splitting
- if(nDomains==1)then
- fluxMask%var(iVar)%dat = desiredFlux
-
- ! domain splitting
- else
- print*,"split domain"
-
- ! initialize to .false.
- fluxMask%var(iVar)%dat = .false.
-
- ! only need to proceed if the flux is desired
- if(desiredFlux)then
-
- ! different domain splitting operations
- select case(iDomainSplit)
-
- ! canopy fluxes -- (:1) gets the upper boundary(0) if it exists
- case(vegSplit)
-
- ! vector solution (should only be present for energy)
- if(ixSolution==vector)then
- fluxMask%var(iVar)%dat(:1) = desiredFlux
- if(ixStateThenDomain>1 .and. iStateTypeSplit/=nrgSplit)then
- message=trim(message)//'only expect a vector solution for the vegetation domain for energy'
- err=20; return
- endif
-
- ! scalar solution
- else
- fluxMask%var(iVar)%dat(:1) = desiredFlux
+ print*,"split domain"
+
+ ! initialize to .false.
+ fluxMask%var(iVar)%dat = .false.
+
+ ! only need to proceed if the flux is desired
+ if(desiredFlux)then
+
+ ! different domain splitting operations
+ select case(iDomainSplit)
+
+ ! canopy fluxes -- (:1) gets the upper boundary(0) if it exists
+ case(vegSplit)
+
+ ! vector solution (should only be present for energy)
+ if(ixSolution==vector)then
+ fluxMask%var(iVar)%dat(:1) = desiredFlux
+ if(ixStateThenDomain>1 .and. iStateTypeSplit/=nrgSplit)then
+ message=trim(message)//'only expect a vector solution for the vegetation domain for energy'
+ err=20; return
endif
-
- ! fluxes through snow and soil
- case(snowSplit,soilSplit)
-
- ! loop through layers
- do iLayer=1,nLayers
- if(ixlayerActive(iLayer)/=integerMissing)then
-
- ! get the offset (ixLayerActive=1,2,3,...nLayers, and soil vectors nSnow+1, nSnow+2, ..., nLayers)
- iOffset = merge(nSnow, 0, flux_meta(iVar)%vartype==iLookVarType%midSoil .or. flux_meta(iVar)%vartype==iLookVarType%ifcSoil)
- jLayer = iLayer-iOffset
-
- ! identify the minimum layer
- select case(flux_meta(iVar)%vartype)
- case(iLookVarType%ifcToto, iLookVarType%ifcSnow, iLookVarType%ifcSoil); minLayer=merge(jLayer-1, jLayer, jLayer==1)
- case(iLookVarType%midToto, iLookVarType%midSnow, iLookVarType%midSoil); minLayer=jLayer
- case default; minLayer=integerMissing
- end select
-
- ! set desired layers
- select case(flux_meta(iVar)%vartype)
- case(iLookVarType%midToto,iLookVarType%ifcToto); fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
- case(iLookVarType%midSnow,iLookVarType%ifcSnow); if(iLayer<=nSnow) fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
- case(iLookVarType%midSoil,iLookVarType%ifcSoil); if(iLayer> nSnow) fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
+
+ ! scalar solution
+ else
+ fluxMask%var(iVar)%dat(:1) = desiredFlux
+ endif
+
+ ! fluxes through snow and soil
+ case(snowSplit,soilSplit)
+
+ ! loop through layers
+ do iLayer=1,nLayers
+ if(ixlayerActive(iLayer)/=integerMissing)then
+
+ ! get the offset (ixLayerActive=1,2,3,...nLayers, and soil vectors nSnow+1, nSnow+2, ..., nLayers)
+ iOffset = merge(nSnow, 0, flux_meta(iVar)%vartype==iLookVarType%midSoil .or. flux_meta(iVar)%vartype==iLookVarType%ifcSoil)
+ jLayer = iLayer-iOffset
+
+ ! identify the minimum layer
+ select case(flux_meta(iVar)%vartype)
+ case(iLookVarType%ifcToto, iLookVarType%ifcSnow, iLookVarType%ifcSoil); minLayer=merge(jLayer-1, jLayer, jLayer==1)
+ case(iLookVarType%midToto, iLookVarType%midSnow, iLookVarType%midSoil); minLayer=jLayer
+ case default; minLayer=integerMissing
+ end select
+
+ ! set desired layers
+ select case(flux_meta(iVar)%vartype)
+ case(iLookVarType%midToto,iLookVarType%ifcToto); fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
+ case(iLookVarType%midSnow,iLookVarType%ifcSnow); if(iLayer<=nSnow) fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
+ case(iLookVarType%midSoil,iLookVarType%ifcSoil); if(iLayer> nSnow) fluxMask%var(iVar)%dat(minLayer:jLayer) = desiredFlux
+ end select
+
+ ! add hydrology states for scalar variables
+ if(iStateTypeSplit==massSplit .and. flux_meta(iVar)%vartype==iLookVarType%scalarv)then
+ select case(iDomainSplit)
+ case(snowSplit); if(iLayer==nSnow) fluxMask%var(iVar)%dat = desiredFlux
+ case(soilSplit); if(iLayer==nSnow+1) fluxMask%var(iVar)%dat = desiredFlux
end select
-
- ! add hydrology states for scalar variables
- if(iStateTypeSplit==massSplit .and. flux_meta(iVar)%vartype==iLookVarType%scalarv)then
- select case(iDomainSplit)
- case(snowSplit); if(iLayer==nSnow) fluxMask%var(iVar)%dat = desiredFlux
- case(soilSplit); if(iLayer==nSnow+1) fluxMask%var(iVar)%dat = desiredFlux
- end select
- endif ! if hydrology split and scalar
-
- endif ! if the layer is active
- end do ! looping through layers
-
- ! check
- case default; err=20; message=trim(message)//'unable to identify split based on domain type'; return
- end select ! domain split
-
- endif ! if flux is desired
-
- endif ! domain splitting
- endif ! not fully coupled
-
- ! define if the flux is desired
- if(desiredFlux) neededFlux(iVar)=.true.
- !if(desiredFlux) print*, flux_meta(iVar)%varname, fluxMask%var(iVar)%dat
-
- ! * check
- if( globalPrintFlag .and. count(fluxMask%var(iVar)%dat)>0 )&
- print*, trim(flux_meta(iVar)%varname)
-
- end do ! (loop through fluxes)
-
- end associate stateSubset
-
- ! *******************************************************************************************************************************
- ! *******************************************************************************************************************************
- ! *******************************************************************************************************************************
- ! ***** trial with a given solution method...
-
- ! check that we do not attempt the scalar solution for the fully coupled case
- if(ixCoupling==fullyCoupled .and. ixSolution==scalar)then
- message=trim(message)//'only apply the scalar solution to the fully split coupling strategy'
- err=20; return
- endif
-
- ! reset the flag for the first flux call
- if(.not.firstSuccess) firstFluxCall=.true.
-
- ! save/recover copies of prognostic variables
- do iVar=1,size(prog_data%var)
- select case(failure)
- case(.false.); prog_temp%var(iVar)%dat(:) = prog_data%var(iVar)%dat(:)
- case(.true.); prog_data%var(iVar)%dat(:) = prog_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! save/recover copies of diagnostic variables
- do iVar=1,size(diag_data%var)
- select case(failure)
- case(.false.); diag_temp%var(iVar)%dat(:) = diag_data%var(iVar)%dat(:)
- case(.true.); diag_data%var(iVar)%dat(:) = diag_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! save/recover copies of model fluxes
- do iVar=1,size(flux_data%var)
- select case(failure)
- case(.false.); flux_temp%var(iVar)%dat(:) = flux_data%var(iVar)%dat(:)
- case(.true.); flux_data%var(iVar)%dat(:) = flux_temp%var(iVar)%dat(:)
- end select
- end do ! looping through variables
-
- ! -----
- ! * solve variable subset for one time step...
- ! --------------------------------------------
-
- ! keep track of the number of scalar solutions
- if(ixSolution==scalar) numberScalarSolutions = numberScalarSolutions + 1
-
- ! solve variable subset for one full time step
- select case(ixNumericalMethod)
- case(sundials)
- call varSubstepSundials(&
- ! input: model control
- dt, & ! intent(in) : time step (s)
- dtInit, & ! intent(in) : initial time step (seconds)
- dt_min, & ! intent(in) : minimum time step (seconds)
- nSubset, & ! intent(in) : total number of variables in the state subset
- doAdjustTemp, & ! intent(in) : flag to indicate if we adjust the temperature
- firstSubStep, & ! intent(in) : flag to denote first sub-step
- firstFluxCall, & ! intent(inout) : flag to indicate if we are processing the first flux call
- computeVegFlux, & ! intent(in) : flag to denote if computing energy flux over vegetation
- (ixSolution==scalar), & ! intent(in) : flag to denote computing the scalar solution
- iStateSplit, & ! intent(in) : index of the layer in the splitting operation
- fluxMask, & ! intent(in) : mask for the fluxes used in this given state subset
- fluxCount, & ! intent(inout) : number of times fluxes are updated (should equal nsubstep)
- ! input/output: data structures
- model_decisions, & ! intent(in) : model decisions
- lookup_data, & ! intent(in) : lookup tables
- type_data, & ! intent(in) : type of vegetation and soil
- attr_data, & ! intent(in) : spatial attributes
- forc_data, & ! intent(in) : model forcing data
- mpar_data, & ! intent(in) : model parameters
- indx_data, & ! intent(inout) : index data
- prog_data, & ! intent(inout) : model prognostic variables for a local HRU
- diag_data, & ! intent(inout) : model diagnostic variables for a local HRU
- flux_data, & ! intent(inout) : model fluxes for a local HRU
- deriv_data, & ! intent(inout) : derivatives in model fluxes w.r.t. relevant state variables
- bvar_data, & ! intent(in) : model variables for the local basin
- ! output: control
- ixSaturation, & ! intent(inout) : index of the lowest saturated layer (NOTE: only computed on the first iteration)
- dtMultiplier, & ! intent(out) : substep multiplier (-)
- nSubsteps, & ! intent(out) : number of substeps taken for a given split
- failedMinimumStep, & ! intent(out) : flag for failed substeps
- reduceCoupledStep, & ! intent(out) : flag to reduce the length of the coupled step
- tooMuchMelt, & ! intent(out) : flag to denote that ice is insufficient to support melt
- dt_out, & ! intent(out)
- err,cmessage) ! intent(out) : error code and error message
-
- dt = dt_out
-
- case(bEuler)
- call varSubstep(&
- ! input: model control
- dt, & ! intent(in) : time step (s)
- dtInit, & ! intent(in) : initial time step (seconds)
- dt_min, & ! intent(in) : minimum time step (seconds)
- nSubset, & ! intent(in) : total number of variables in the state subset
- doAdjustTemp, & ! intent(in) : flag to indicate if we adjust the temperature
- firstSubStep, & ! intent(in) : flag to denote first sub-step
- firstFluxCall, & ! intent(inout) : flag to indicate if we are processing the first flux call
- computeVegFlux, & ! intent(in) : flag to denote if computing energy flux over vegetation
- (ixSolution==scalar), & ! intent(in) : flag to denote computing the scalar solution
- iStateSplit, & ! intent(in) : index of the layer in the splitting operation
- fluxMask, & ! intent(in) : mask for the fluxes used in this given state subset
- fluxCount, & ! intent(inout) : number of times fluxes are updated (should equal nsubstep)
- ! input/output: data structures
- model_decisions, & ! intent(in) : model decisions
- lookup_data, & ! intent(in) : lookup tables
- type_data, & ! intent(in) : type of vegetation and soil
- attr_data, & ! intent(in) : spatial attributes
- forc_data, & ! intent(in) : model forcing data
- mpar_data, & ! intent(in) : model parameters
- indx_data, & ! intent(inout) : index data
- prog_data, & ! intent(inout) : model prognostic variables for a local HRU
- diag_data, & ! intent(inout) : model diagnostic variables for a local HRU
- flux_data, & ! intent(inout) : model fluxes for a local HRU
- deriv_data, & ! intent(inout) : derivatives in model fluxes w.r.t. relevant state variables
- bvar_data, & ! intent(in) : model variables for the local basin
- ! output: control
- ixSaturation, & ! intent(inout) : index of the lowest saturated layer (NOTE: only computed on the first iteration)
- dtMultiplier, & ! intent(out) : substep multiplier (-)
- nSubsteps, & ! intent(out) : number of substeps taken for a given split
- failedMinimumStep, & ! intent(out) : flag for failed substeps
- reduceCoupledStep, & ! intent(out) : flag to reduce the length of the coupled step
- tooMuchMelt, & ! intent(out) : flag to denote that ice is insufficient to support melt
- err,cmessage) ! intent(out) : error code and error message
- ! check
- case default; err=20; message=trim(message)//'expect num_method to be sundials or bEuler (or itertive, which is bEuler)'; return
+ endif ! if hydrology split and scalar
+
+ endif ! if the layer is active
+ end do ! looping through layers
+
+ ! check
+ case default; err=20; message=trim(message)//'unable to identify split based on domain type'; return
+ end select ! domain split
+
+ endif ! if flux is desired
+
+ endif ! domain splitting
+ endif ! not fully coupled
+
+ ! define if the flux is desired
+ if(desiredFlux) neededFlux(iVar)=.true.
+ !if(desiredFlux) print*, flux_meta(iVar)%varname, fluxMask%var(iVar)%dat
+
+ ! * check
+ if( globalPrintFlag .and. count(fluxMask%var(iVar)%dat)>0 )&
+ print*, trim(flux_meta(iVar)%varname)
+
+ end do ! (loop through fluxes)
+
+ end associate stateSubset
+
+ ! *******************************************************************************************************************************
+ ! *******************************************************************************************************************************
+ ! *******************************************************************************************************************************
+ ! ***** trial with a given solution method...
+
+ ! check that we do not attempt the scalar solution for the fully coupled case
+ if(ixCoupling==fullyCoupled .and. ixSolution==scalar)then
+ message=trim(message)//'only apply the scalar solution to the fully split coupling strategy'
+ err=20; return
+ endif
+
+ ! reset the flag for the first flux call
+ if(.not.firstSuccess) firstFluxCall=.true.
+
+ ! save/recover copies of prognostic variables
+ do iVar=1,size(prog_data%var)
+ select case(failure)
+ case(.false.); prog_temp%var(iVar)%dat(:) = prog_data%var(iVar)%dat(:)
+ case(.true.); prog_data%var(iVar)%dat(:) = prog_temp%var(iVar)%dat(:)
end select
-
- if(err/=0)then
- message=trim(message)//trim(cmessage)
- if(err>0) return
- endif ! (check for errors)
-
- ! reduce coupled step if failed the minimum step for the scalar solution
- if(failedMinimumStep .and. ixSolution==scalar) reduceCoupledStep=.true.
-
- ! if too much melt (or some other need to reduce the coupled step) then return
- ! NOTE: need to go all the way back to coupled_em and merge snow layers, as all splitting operations need to occur with the same layer geometry
- if(tooMuchMelt .or. reduceCoupledStep)then
- stepFailure=.true.
- err=0 ! recovering
- return
- endif
-
- ! define failure
- failure = (failedMinimumStep .or. err<0)
- if(.not.failure) firstSuccess=.true.
-
- ! if failed, need to reset the flux counter
- if(failure)then
- do iVar=1,size(flux_meta)
- iMin=lbound(flux_data%var(iVar)%dat)
- iMax=ubound(flux_data%var(iVar)%dat)
- do iLayer=iMin(1),iMax(1)
- if(fluxMask%var(iVar)%dat(iLayer)) fluxCount%var(iVar)%dat(iLayer) = fluxCount%var(iVar)%dat(iLayer) - nSubsteps
- end do
+ end do ! looping through variables
+
+ ! save/recover copies of diagnostic variables
+ do iVar=1,size(diag_data%var)
+ select case(failure)
+ case(.false.); diag_temp%var(iVar)%dat(:) = diag_data%var(iVar)%dat(:)
+ case(.true.); diag_data%var(iVar)%dat(:) = diag_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ ! save/recover copies of model fluxes
+ do iVar=1,size(flux_data%var)
+ select case(failure)
+ case(.false.); flux_temp%var(iVar)%dat(:) = flux_data%var(iVar)%dat(:)
+ case(.true.); flux_data%var(iVar)%dat(:) = flux_temp%var(iVar)%dat(:)
+ end select
+ end do ! looping through variables
+
+ ! -----
+ ! * solve variable subset for one time step...
+ ! --------------------------------------------
+
+ ! keep track of the number of scalar solutions
+ if(ixSolution==scalar) numberScalarSolutions = numberScalarSolutions + 1
+
+ ! solve variable subset for one full time step
+ select case(ixNumericalMethod)
+ case(sundials)
+ call varSubstepSundials(&
+ ! input: model control
+ dt, & ! intent(in) : time step (s)
+ dtInit, & ! intent(in) : initial time step (seconds)
+ dt_min, & ! intent(in) : minimum time step (seconds)
+ nSubset, & ! intent(in) : total number of variables in the state subset
+ doAdjustTemp, & ! intent(in) : flag to indicate if we adjust the temperature
+ firstSubStep, & ! intent(in) : flag to denote first sub-step
+ firstFluxCall, & ! intent(inout) : flag to indicate if we are processing the first flux call
+ computeVegFlux, & ! intent(in) : flag to denote if computing energy flux over vegetation
+ (ixSolution==scalar), & ! intent(in) : flag to denote computing the scalar solution
+ iStateSplit, & ! intent(in) : index of the layer in the splitting operation
+ fluxMask, & ! intent(in) : mask for the fluxes used in this given state subset
+ fluxCount, & ! intent(inout) : number of times fluxes are updated (should equal nsubstep)
+ ! input/output: data structures
+ model_decisions, & ! intent(in) : model decisions
+ lookup_data, & ! intent(in) : lookup tables
+ type_data, & ! intent(in) : type of vegetation and soil
+ attr_data, & ! intent(in) : spatial attributes
+ forc_data, & ! intent(in) : model forcing data
+ mpar_data, & ! intent(in) : model parameters
+ indx_data, & ! intent(inout) : index data
+ prog_data, & ! intent(inout) : model prognostic variables for a local HRU
+ diag_data, & ! intent(inout) : model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout) : model fluxes for a local HRU
+ deriv_data, & ! intent(inout) : derivatives in model fluxes w.r.t. relevant state variables
+ bvar_data, & ! intent(in) : model variables for the local basin
+ ! output: control
+ ixSaturation, & ! intent(inout) : index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ dtMultiplier, & ! intent(out) : substep multiplier (-)
+ nSubsteps, & ! intent(out) : number of substeps taken for a given split
+ failedMinimumStep, & ! intent(out) : flag for failed substeps
+ reduceCoupledStep, & ! intent(out) : flag to reduce the length of the coupled step
+ tooMuchMelt, & ! intent(out) : flag to denote that ice is insufficient to support melt
+ dt_out, & ! intent(out)
+ err,cmessage) ! intent(out) : error code and error message
+
+ dt = dt_out
+
+ case(bEuler)
+ call varSubstep(&
+ ! input: model control
+ dt, & ! intent(in) : time step (s)
+ dtInit, & ! intent(in) : initial time step (seconds)
+ dt_min, & ! intent(in) : minimum time step (seconds)
+ nSubset, & ! intent(in) : total number of variables in the state subset
+ doAdjustTemp, & ! intent(in) : flag to indicate if we adjust the temperature
+ firstSubStep, & ! intent(in) : flag to denote first sub-step
+ firstFluxCall, & ! intent(inout) : flag to indicate if we are processing the first flux call
+ computeVegFlux, & ! intent(in) : flag to denote if computing energy flux over vegetation
+ (ixSolution==scalar), & ! intent(in) : flag to denote computing the scalar solution
+ iStateSplit, & ! intent(in) : index of the layer in the splitting operation
+ fluxMask, & ! intent(in) : mask for the fluxes used in this given state subset
+ fluxCount, & ! intent(inout) : number of times fluxes are updated (should equal nsubstep)
+ ! input/output: data structures
+ model_decisions, & ! intent(in) : model decisions
+ lookup_data, & ! intent(in) : lookup tables
+ type_data, & ! intent(in) : type of vegetation and soil
+ attr_data, & ! intent(in) : spatial attributes
+ forc_data, & ! intent(in) : model forcing data
+ mpar_data, & ! intent(in) : model parameters
+ indx_data, & ! intent(inout) : index data
+ prog_data, & ! intent(inout) : model prognostic variables for a local HRU
+ diag_data, & ! intent(inout) : model diagnostic variables for a local HRU
+ flux_data, & ! intent(inout) : model fluxes for a local HRU
+ deriv_data, & ! intent(inout) : derivatives in model fluxes w.r.t. relevant state variables
+ bvar_data, & ! intent(in) : model variables for the local basin
+ ! output: control
+ ixSaturation, & ! intent(inout) : index of the lowest saturated layer (NOTE: only computed on the first iteration)
+ dtMultiplier, & ! intent(out) : substep multiplier (-)
+ nSubsteps, & ! intent(out) : number of substeps taken for a given split
+ failedMinimumStep, & ! intent(out) : flag for failed substeps
+ reduceCoupledStep, & ! intent(out) : flag to reduce the length of the coupled step
+ tooMuchMelt, & ! intent(out) : flag to denote that ice is insufficient to support melt
+ err,cmessage) ! intent(out) : error code and error message
+ ! check
+ case default; err=20; message=trim(message)//'expect num_method to be sundials or bEuler (or itertive, which is bEuler)'; return
+ end select
+
+ if(err/=0)then
+ message=trim(message)//trim(cmessage)
+ if(err>0) return
+ endif ! (check for errors)
+
+ ! reduce coupled step if failed the minimum step for the scalar solution
+ if(failedMinimumStep .and. ixSolution==scalar) reduceCoupledStep=.true.
+
+ ! if too much melt (or some other need to reduce the coupled step) then return
+ ! NOTE: need to go all the way back to coupled_em and merge snow layers, as all splitting operations need to occur with the same layer geometry
+ if(tooMuchMelt .or. reduceCoupledStep)then
+ stepFailure=.true.
+ err=0 ! recovering
+ return
+ endif
+
+ ! define failure
+ failure = (failedMinimumStep .or. err<0)
+ if(.not.failure) firstSuccess=.true.
+
+ ! if failed, need to reset the flux counter
+ if(failure)then
+ do iVar=1,size(flux_meta)
+ iMin=lbound(flux_data%var(iVar)%dat)
+ iMax=ubound(flux_data%var(iVar)%dat)
+ do iLayer=iMin(1),iMax(1)
+ if(fluxMask%var(iVar)%dat(iLayer)) fluxCount%var(iVar)%dat(iLayer) = fluxCount%var(iVar)%dat(iLayer) - nSubsteps
end do
- endif
-
- ! try the fully split solution if failed to converge with a minimum time step in the coupled solution
- if(ixCoupling==fullyCoupled .and. failure) cycle coupling
-
- ! try the scalar solution if failed to converge with a minimum time step in the split solution
- if(ixCoupling/=fullyCoupled)then
- select case(ixStateThenDomain)
- case(fullDomain); if(failure) cycle stateThenDomain
- case(subDomain); if(failure) cycle solution
- case default; err=20; message=trim(message)//'unknown ixStateThenDomain case'
- end select
- endif
-
- ! check that state variables updated
- where(stateMask) stateCheck = stateCheck+1
- if(any(stateCheck>1))then
- message=trim(message)//'state variable updated more than once!'
+ end do
+ endif
+
+ ! try the fully split solution if failed to converge with a minimum time step in the coupled solution
+ if(ixCoupling==fullyCoupled .and. failure) cycle coupling
+
+ ! try the scalar solution if failed to converge with a minimum time step in the split solution
+ if(ixCoupling/=fullyCoupled)then
+ select case(ixStateThenDomain)
+ case(fullDomain); if(failure) cycle stateThenDomain
+ case(subDomain); if(failure) cycle solution
+ case default; err=20; message=trim(message)//'unknown ixStateThenDomain case'
+ end select
+ endif
+
+ ! check that state variables updated
+ where(stateMask) stateCheck = stateCheck+1
+ if(any(stateCheck>1))then
+ message=trim(message)//'state variable updated more than once!'
+ err=20; return
+ endif
+
+ ! success = exit solution
+ if(.not.failure)then
+ select case(ixStateThenDomain)
+ case(fullDomain); if(iStateSplit==nStateSplit) exit stateThenDomain
+ case(subDomain); if(iStateSplit==nStateSplit) exit solution
+ case default; err=20; message=trim(message)//'unknown ixStateThenDomain case'
+ end select
+ else
+
+ ! check that we did not fail for the scalar solution (last resort)
+ if(ixSolution==scalar)then
+ message=trim(message)//'failed the minimum step for the scalar solution'
err=20; return
- endif
-
- ! success = exit solution
- if(.not.failure)then
- select case(ixStateThenDomain)
- case(fullDomain); if(iStateSplit==nStateSplit) exit stateThenDomain
- case(subDomain); if(iStateSplit==nStateSplit) exit solution
- case default; err=20; message=trim(message)//'unknown ixStateThenDomain case'
- end select
+
+ ! check for an unexpected failure
else
-
- ! check that we did not fail for the scalar solution (last resort)
- if(ixSolution==scalar)then
- message=trim(message)//'failed the minimum step for the scalar solution'
- err=20; return
-
- ! check for an unexpected failure
- else
- message=trim(message)//'unexpected failure'
- err=20; return
- endif
-
- endif ! success check
-
- end do stateSplit ! solution with split layers
-
- end do solution ! trial with the full layer solution then the split layer solution
-
- ! ***** trial with a given solution method...
- ! *******************************************************************************************************************************
- ! *******************************************************************************************************************************
- ! *******************************************************************************************************************************
-
- end do domainSplit ! domain type splitting loop
-
-
- end do stateThenDomain ! switch between the state and the domain
-
-
- ! -----
- ! * reset state variables for the mass split...
- ! ---------------------------------------------
- ! modify the state type names associated with the state vector
- if(ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)then
- if(computeVegFlux)then
- where(ixStateType(ixHydCanopy)==iname_liqCanopy) ixStateType(ixHydCanopy)=iname_watCanopy
- endif
- where(ixStateType(ixHydLayer) ==iname_liqLayer) ixStateType(ixHydLayer) =iname_watLayer
- where(ixStateType(ixHydLayer) ==iname_lmpLayer) ixStateType(ixHydLayer) =iname_matLayer
- endif ! if modifying state variables for the mass split
-
- end do stateTypeSplit ! state type splitting loop
-
- ! success = exit the coupling loop
- if(ixCoupling==fullyCoupled .and. .not.failure) exit coupling
-
- end do coupling ! coupling method
-
- ! check that all state variables were updated
- if(any(stateCheck==0))then
- message=trim(message)//'some state variables were not updated!'
+ message=trim(message)//'unexpected failure'
+ err=20; return
+ endif
+
+ endif ! success check
+
+ end do stateSplit ! solution with split layers
+
+ end do solution ! trial with the full layer solution then the split layer solution
+
+ ! ***** trial with a given solution method...
+ ! *******************************************************************************************************************************
+ ! *******************************************************************************************************************************
+ ! *******************************************************************************************************************************
+
+ end do domainSplit ! domain type splitting loop
+
+
+ end do stateThenDomain ! switch between the state and the domain
+
+
+ ! -----
+ ! * reset state variables for the mass split...
+ ! ---------------------------------------------
+ ! modify the state type names associated with the state vector
+ if(ixCoupling/=fullyCoupled .and. iStateTypeSplit==massSplit)then
+ if(computeVegFlux)then
+ where(ixStateType(ixHydCanopy)==iname_liqCanopy) ixStateType(ixHydCanopy)=iname_watCanopy
+ endif
+ where(ixStateType(ixHydLayer) ==iname_liqLayer) ixStateType(ixHydLayer) =iname_watLayer
+ where(ixStateType(ixHydLayer) ==iname_lmpLayer) ixStateType(ixHydLayer) =iname_matLayer
+ endif ! if modifying state variables for the mass split
+
+ end do stateTypeSplit ! state type splitting loop
+
+ ! success = exit the coupling loop
+ if(ixCoupling==fullyCoupled .and. .not.failure) exit coupling
+
+ end do coupling ! coupling method
+
+ ! check that all state variables were updated
+ if(any(stateCheck==0))then
+ message=trim(message)//'some state variables were not updated!'
+ err=20; return
+ endif
+
+ ! check that the desired fluxes were computed
+ do iVar=1,size(flux_meta)
+ if(neededFlux(iVar) .and. any(fluxCount%var(iVar)%dat==0))then
+ print*, 'fluxCount%var(iVar)%dat = ', fluxCount%var(iVar)%dat
+ message=trim(message)//'flux '//trim(flux_meta(iVar)%varname)//' was not computed'
err=20; return
endif
-
- ! check that the desired fluxes were computed
- do iVar=1,size(flux_meta)
- if(neededFlux(iVar) .and. any(fluxCount%var(iVar)%dat==0))then
- print*, 'fluxCount%var(iVar)%dat = ', fluxCount%var(iVar)%dat
- message=trim(message)//'flux '//trim(flux_meta(iVar)%varname)//' was not computed'
- err=20; return
- endif
- end do
-
- ! use step halving if unable to complete the fully coupled solution in one substep
- if(ixCoupling/=fullyCoupled .or. nSubsteps>1) dtMultiplier=0.5_rkind
-
- ! compute the melt in each snow and soil layer
- if(nSnow>0)then
- diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(1:nSnow) = -( mLayerVolFracIce(1:nSnow) - mLayerVolFracIceInit(1:nSnow) ) * iden_ice
- diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(nSnow+1:nLayers) = -(mLayerVolFracIce(nSnow+1:nLayers) - mLayerVolFracIceInit(nSnow+1:nLayers))*iden_water
- endif
-
- ! end associate statements
- end associate globalVars
-
- end subroutine opSplittin
-
-
- ! **********************************************************************************************************
- ! private subroutine stateFilter: get a mask for the desired state variables
- ! **********************************************************************************************************
- subroutine stateFilter(ixCoupling,ixSolution,ixStateThenDomain,iStateTypeSplit,iDomainSplit,iStateSplit,&
- indx_data,stateMask,nSubset,err,message)
-
- USE indexState_module,only:indxSubset ! get state indices
- implicit none
- ! input
- integer(i4b),intent(in) :: ixCoupling ! index of coupling method (1,2)
- integer(i4b),intent(in) :: ixSolution ! index of solution method (1,2)
- integer(i4b),intent(in) :: ixStateThenDomain ! switch between full domain and sub domains
- integer(i4b),intent(in) :: iStateTypeSplit ! index of the state type split
- integer(i4b),intent(in) :: iDomainSplit ! index of the domain split
- integer(i4b),intent(in) :: iStateSplit ! index of the layer split
- type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
- ! output
- logical(lgt),intent(out) :: stateMask(:) ! mask defining desired state variables
- integer(i4b),intent(out) :: nSubset ! number of selected state variables for a given split
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! local
- integer(i4b),allocatable :: ixSubset(:) ! list of indices in the state subset
- character(len=256) :: cmessage ! error message
+ end do
+
+ ! use step halving if unable to complete the fully coupled solution in one substep
+ if(ixCoupling/=fullyCoupled .or. nSubsteps>1) dtMultiplier=0.5_rkind
+
+ ! compute the melt in each snow and soil layer
+ if(nSnow>0)then
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(1:nSnow) = -( mLayerVolFracIce(1:nSnow) - mLayerVolFracIceInit(1:nSnow) ) * iden_ice
+ diag_data%var(iLookDIAG%mLayerMeltFreeze)%dat(nSnow+1:nLayers) = -(mLayerVolFracIce(nSnow+1:nLayers) - mLayerVolFracIceInit(nSnow+1:nLayers))*iden_water
+ endif
+
+ ! end associate statements
+ end associate globalVars
+
+end subroutine opSplittin
+
+
+! **********************************************************************************************************
+! private subroutine stateFilter: get a mask for the desired state variables
+! **********************************************************************************************************
+subroutine stateFilter(ixCoupling,ixSolution,ixStateThenDomain,iStateTypeSplit,iDomainSplit,iStateSplit,&
+ indx_data,stateMask,nSubset,err,message)
+
+ USE indexState_module,only:indxSubset ! get state indices
+ implicit none
+ ! input
+ integer(i4b),intent(in) :: ixCoupling ! index of coupling method (1,2)
+ integer(i4b),intent(in) :: ixSolution ! index of solution method (1,2)
+ integer(i4b),intent(in) :: ixStateThenDomain ! switch between full domain and sub domains
+ integer(i4b),intent(in) :: iStateTypeSplit ! index of the state type split
+ integer(i4b),intent(in) :: iDomainSplit ! index of the domain split
+ integer(i4b),intent(in) :: iStateSplit ! index of the layer split
+ type(var_ilength),intent(inout) :: indx_data ! indices for a local HRU
+ ! output
+ logical(lgt),intent(out) :: stateMask(:) ! mask defining desired state variables
+ integer(i4b),intent(out) :: nSubset ! number of selected state variables for a given split
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! local
+ integer(i4b),allocatable :: ixSubset(:) ! list of indices in the state subset
+ character(len=256) :: cmessage ! error message
+ ! --------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ ! data structures
+ associate(&
+ ! indices of model state variables
+ ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
+ ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
+ ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
+ ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
+ ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
+ ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
+ ixWatAquifer => indx_data%var(iLookINDEX%ixWatAquifer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for water storage in the aquifer
+ ixAllState => indx_data%var(iLookINDEX%ixAllState)%dat ,& ! intent(in): [i4b(:)] list of indices for all model state variables (1,2,3,...nState)
+ ! 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
+ ) ! data structures
! --------------------------------------------------------------------------------------------------------------------------------------------------------------------------
- ! data structures
- associate(&
- ! indices of model state variables
- ixStateType => indx_data%var(iLookINDEX%ixStateType)%dat ,& ! intent(in): [i4b(:)] indices defining the type of the state (ixNrgState...)
- ixNrgCanair => indx_data%var(iLookINDEX%ixNrgCanair)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in canopy air space domain
- ixNrgCanopy => indx_data%var(iLookINDEX%ixNrgCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the canopy domain
- ixHydCanopy => indx_data%var(iLookINDEX%ixHydCanopy)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the canopy domain
- ixNrgLayer => indx_data%var(iLookINDEX%ixNrgLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for energy states in the snow+soil domain
- ixHydLayer => indx_data%var(iLookINDEX%ixHydLayer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for hydrology states in the snow+soil domain
- ixWatAquifer => indx_data%var(iLookINDEX%ixWatAquifer)%dat ,& ! intent(in): [i4b(:)] indices IN THE FULL VECTOR for water storage in the aquifer
- ixAllState => indx_data%var(iLookINDEX%ixAllState)%dat ,& ! intent(in): [i4b(:)] list of indices for all model state variables (1,2,3,...nState)
- ! 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
- ) ! data structures
- ! --------------------------------------------------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message='stateFilter/'
-
- ! identify splitting option
- select case(ixCoupling)
-
- ! -----
- ! - fully coupled...
- ! ------------------
-
- ! use all state variables
- case(fullyCoupled); stateMask(:) = .true.
-
- ! -----
- ! - splitting by state type...
- ! ----------------------------
-
- ! initial split by state type
- case(stateTypeSplit)
-
- ! switch between full domain and sub domains
- select case(ixStateThenDomain)
-
- ! split into energy and mass
- case(fullDomain)
- select case(iStateTypeSplit)
- case(nrgSplit); stateMask = (ixStateType==iname_nrgCanair .or. ixStateType==iname_nrgCanopy .or. ixStateType==iname_nrgLayer)
- case(massSplit); stateMask = (ixStateType==iname_liqCanopy .or. ixStateType==iname_liqLayer .or. ixStateType==iname_lmpLayer .or. ixStateType==iname_watAquifer)
- case default; err=20; message=trim(message)//'unable to identify split based on state type'; return
- end select
-
- ! split into vegetation, snow, and soil
- case(subDomain)
-
- ! define state mask
- stateMask=.false. ! (initialize state mask)
- select case(iStateTypeSplit)
-
- ! define mask for energy
- case(nrgSplit)
- select case(iDomainSplit)
- case(vegSplit)
- if(ixNrgCanair(1)/=integerMissing) stateMask(ixNrgCanair) = .true. ! energy of the canopy air space
- if(ixNrgCanopy(1)/=integerMissing) stateMask(ixNrgCanopy) = .true. ! energy of the vegetation canopy
- stateMask(ixNrgLayer(1)) = .true. ! energy of the upper-most layer in the snow+soil domain
- case(snowSplit); if(nSnow>1) stateMask(ixNrgLayer(2:nSnow)) = .true. ! NOTE: (2:) because the top layer in the snow+soil domain included in vegSplit
- case(soilSplit); stateMask(ixNrgLayer(max(2,nSnow+1):nLayers)) = .true. ! NOTE: max(2,nSnow+1) gives second layer unless more than 2 snow layers
- case(aquiferSplit) ! do nothing: no energy state variable for the aquifer domain
- case default; err=20; message=trim(message)//'unable to identify model sub-domain'; return
- end select
-
- ! define mask for water
- case(massSplit)
- select case(iDomainSplit)
- case(vegSplit); if(ixHydCanopy(1)/=integerMissing) stateMask(ixHydCanopy) = .true. ! hydrology of the vegetation canopy
- case(snowSplit); stateMask(ixHydLayer(1:nSnow)) = .true. ! snow hydrology
- case(soilSplit); stateMask(ixHydLayer(nSnow+1:nLayers)) = .true. ! soil hydrology
- case(aquiferSplit); if(ixWatAquifer(1)/=integerMissing) stateMask(ixWatAquifer) = .true. ! aquifer storage
- case default; err=20; message=trim(message)//'unable to identify model sub-domain'; return
- end select
-
- ! check
- case default; err=20; message=trim(message)//'unable to identify the state type'; return
- end select ! (split based on state type)
-
- ! check
- case default; err=20; message=trim(message)//'unable to identify the switch between full domains and sub domains'; return
- end select ! (switch between full domains and sub domains)
-
+ ! initialize error control
+ err=0; message='stateFilter/'
+
+ ! identify splitting option
+ select case(ixCoupling)
+
+ ! -----
+ ! - fully coupled...
+ ! ------------------
+
+ ! use all state variables
+ case(fullyCoupled); stateMask(:) = .true.
+
+ ! -----
+ ! - splitting by state type...
+ ! ----------------------------
+
+ ! initial split by state type
+ case(stateTypeSplit)
+
+ ! switch between full domain and sub domains
+ select case(ixStateThenDomain)
+
+ ! split into energy and mass
+ case(fullDomain)
+ select case(iStateTypeSplit)
+ case(nrgSplit); stateMask = (ixStateType==iname_nrgCanair .or. ixStateType==iname_nrgCanopy .or. ixStateType==iname_nrgLayer)
+ case(massSplit); stateMask = (ixStateType==iname_liqCanopy .or. ixStateType==iname_liqLayer .or. ixStateType==iname_lmpLayer .or. ixStateType==iname_watAquifer)
+ case default; err=20; message=trim(message)//'unable to identify split based on state type'; return
+ end select
+
+ ! split into vegetation, snow, and soil
+ case(subDomain)
+
+ ! define state mask
+ stateMask=.false. ! (initialize state mask)
+ select case(iStateTypeSplit)
+
+ ! define mask for energy
+ case(nrgSplit)
+ select case(iDomainSplit)
+ case(vegSplit)
+ if(ixNrgCanair(1)/=integerMissing) stateMask(ixNrgCanair) = .true. ! energy of the canopy air space
+ if(ixNrgCanopy(1)/=integerMissing) stateMask(ixNrgCanopy) = .true. ! energy of the vegetation canopy
+ stateMask(ixNrgLayer(1)) = .true. ! energy of the upper-most layer in the snow+soil domain
+ case(snowSplit); if(nSnow>1) stateMask(ixNrgLayer(2:nSnow)) = .true. ! NOTE: (2:) because the top layer in the snow+soil domain included in vegSplit
+ case(soilSplit); stateMask(ixNrgLayer(max(2,nSnow+1):nLayers)) = .true. ! NOTE: max(2,nSnow+1) gives second layer unless more than 2 snow layers
+ case(aquiferSplit) ! do nothing: no energy state variable for the aquifer domain
+ case default; err=20; message=trim(message)//'unable to identify model sub-domain'; return
+ end select
+
+ ! define mask for water
+ case(massSplit)
+ select case(iDomainSplit)
+ case(vegSplit); if(ixHydCanopy(1)/=integerMissing) stateMask(ixHydCanopy) = .true. ! hydrology of the vegetation canopy
+ case(snowSplit); stateMask(ixHydLayer(1:nSnow)) = .true. ! snow hydrology
+ case(soilSplit); stateMask(ixHydLayer(nSnow+1:nLayers)) = .true. ! soil hydrology
+ case(aquiferSplit); if(ixWatAquifer(1)/=integerMissing) stateMask(ixWatAquifer) = .true. ! aquifer storage
+ case default; err=20; message=trim(message)//'unable to identify model sub-domain'; return
+ end select
+
+ ! check
+ case default; err=20; message=trim(message)//'unable to identify the state type'; return
+ end select ! (split based on state type)
+
! check
- case default; err=20; message=trim(message)//'unable to identify coupling method'; return
- end select ! (selecting solution method)
-
- ! identify scalar solutions
- if(ixSolution==scalar)then
-
- ! get the subset of indices
- call indxSubset(ixSubset, ixAllState, stateMask, err, cmessage)
- if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
-
- ! get the mask
- stateMask(:) = .false.
- stateMask( ixSubset(iStateSplit) ) = .true.
-
+ case default; err=20; message=trim(message)//'unable to identify the switch between full domains and sub domains'; return
+ end select ! (switch between full domains and sub domains)
+
! check
- if(count(stateMask)/=1)then
- message=trim(message)//'expect size=1 (scalar)'
- err=20; return
- endif
-
+ case default; err=20; message=trim(message)//'unable to identify coupling method'; return
+ end select ! (selecting solution method)
+
+ ! identify scalar solutions
+ if(ixSolution==scalar)then
+
+ ! get the subset of indices
+ call indxSubset(ixSubset, ixAllState, stateMask, err, cmessage)
+ if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+
+ ! get the mask
+ stateMask(:) = .false.
+ stateMask( ixSubset(iStateSplit) ) = .true.
+
+ ! check
+ if(count(stateMask)/=1)then
+ message=trim(message)//'expect size=1 (scalar)'
+ err=20; return
endif
-
- ! get the number of selected state variables
- nSubset = count(stateMask)
-
- ! end associations
- end associate
-
- end subroutine stateFilter
-
- end module opSplittin_module
-
\ No newline at end of file
+
+ endif
+
+ ! get the number of selected state variables
+ nSubset = count(stateMask)
+
+ ! end associations
+ end associate
+
+end subroutine stateFilter
+
+end module opSplittin_module
diff --git a/build/source/engine/snow_utils.f90 b/build/source/engine/snow_utils.f90
index 496ecfc..2008a9a 100755
--- a/build/source/engine/snow_utils.f90
+++ b/build/source/engine/snow_utils.f90
@@ -42,73 +42,73 @@ public::tcond_snow
contains
- ! ***********************************************************************************************************
- ! public function fracliquid: compute fraction of liquid water
- ! ***********************************************************************************************************
- function fracliquid(Tk,fc_param)
- implicit none
- real(dp),intent(in) :: Tk ! temperature (K)
- real(dp),intent(in) :: fc_param ! freezing curve parameter (K-1)
- real(dp) :: fracliquid ! fraction of liquid water (-)
- ! compute fraction of liquid water (-)
- fracliquid = 1._dp / ( 1._dp + (fc_param*( Tfreeze - min(Tk,Tfreeze) ))**2._dp )
- end function fracliquid
+! ***********************************************************************************************************
+! public function fracliquid: compute fraction of liquid water
+! ***********************************************************************************************************
+function fracliquid(Tk,fc_param)
+ implicit none
+ real(rkind),intent(in) :: Tk ! temperature (K)
+ real(rkind),intent(in) :: fc_param ! freezing curve parameter (K-1)
+ real(rkind) :: fracliquid ! fraction of liquid water (-)
+ ! compute fraction of liquid water (-)
+ fracliquid = 1._rkind / ( 1._rkind + (fc_param*( Tfreeze - min(Tk,Tfreeze) ))**2._rkind )
+end function fracliquid
- ! ***********************************************************************************************************
- ! public function templiquid: invert the fraction of liquid water function
- ! ***********************************************************************************************************
- function templiquid(fracliquid,fc_param)
- implicit none
- real(dp),intent(in) :: fracliquid ! fraction of liquid water (-)
- real(dp),intent(in) :: fc_param ! freezing curve parameter (K-1)
- real(dp) :: templiquid ! temperature (K)
- ! compute temperature based on the fraction of liquid water (K)
- templiquid = Tfreeze - ((1._dp/fracliquid - 1._dp)/fc_param**2._dp)**(0.5_dp)
- end function templiquid
+! ***********************************************************************************************************
+! public function templiquid: invert the fraction of liquid water function
+! ***********************************************************************************************************
+function templiquid(fracliquid,fc_param)
+ implicit none
+ real(rkind),intent(in) :: fracliquid ! fraction of liquid water (-)
+ real(rkind),intent(in) :: fc_param ! freezing curve parameter (K-1)
+ real(rkind) :: templiquid ! temperature (K)
+ ! compute temperature based on the fraction of liquid water (K)
+ templiquid = Tfreeze - ((1._rkind/fracliquid - 1._rkind)/fc_param**2._rkind)**(0.5_rkind)
+end function templiquid
- ! ***********************************************************************************************************
- ! public function dFracLiq_dTk: differentiate the freezing curve
- ! ***********************************************************************************************************
- function dFracLiq_dTk(Tk,fc_param)
- implicit none
- ! dummies
- real(dp),intent(in) :: Tk ! temperature (K)
- real(dp),intent(in) :: fc_param ! freezing curve parameter (K-1)
- real(dp) :: dFracLiq_dTk ! differentiate the freezing curve (K-1)
- ! locals
- real(dp) :: Tdep ! temperature depression (K)
- real(dp) :: Tdim ! dimensionless temperature (-)
- ! compute local variables (just to make things more efficient)
- Tdep = Tfreeze - min(Tk,Tfreeze)
- Tdim = fc_param*Tdep
- ! differentiate the freezing curve w.r.t temperature
- dFracLiq_dTk = (fc_param*2._dp*Tdim) / ( ( 1._dp + Tdim**2._dp)**2._dp )
- end function dFracLiq_dTk
+! ***********************************************************************************************************
+! public function dFracLiq_dTk: differentiate the freezing curve
+! ***********************************************************************************************************
+function dFracLiq_dTk(Tk,fc_param)
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: Tk ! temperature (K)
+ real(rkind),intent(in) :: fc_param ! freezing curve parameter (K-1)
+ real(rkind) :: dFracLiq_dTk ! differentiate the freezing curve (K-1)
+ ! locals
+ real(rkind) :: Tdep ! temperature depression (K)
+ real(rkind) :: Tdim ! dimensionless temperature (-)
+ ! compute local variables (just to make things more efficient)
+ Tdep = Tfreeze - min(Tk,Tfreeze)
+ Tdim = fc_param*Tdep
+ ! differentiate the freezing curve w.r.t temperature
+ dFracLiq_dTk = (fc_param*2._rkind*Tdim) / ( ( 1._rkind + Tdim**2._rkind)**2._rkind )
+end function dFracLiq_dTk
- ! ***********************************************************************************************************
- ! public subroutine tcond_snow: compute thermal conductivity of snow
- ! ***********************************************************************************************************
- subroutine tcond_snow(BulkDenIce,thermlcond,err,message)
- implicit none
- real(dp),intent(in) :: BulkDenIce ! bulk density of ice (kg m-3)
- real(dp),intent(out) :: thermlcond ! thermal conductivity of snow (W m-1 K-1)
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! initialize error control
- err=0; message="tcond_snow/"
- ! compute thermal conductivity of snow
- select case(model_decisions(iLookDECISIONS%thCondSnow)%iDecision)
- case(Yen1965); thermlcond = 3.217d-6 * BulkDenIce**2._dp ! Yen (1965)
- case(Mellor1977); thermlcond = 2.576d-6 * BulkDenIce**2._dp + 7.4d-2 ! Mellor (1977)
- case(Jordan1991); thermlcond = lambda_air + (7.75d-5*BulkDenIce + 1.105d-6*(BulkDenIce**2._dp)) &
- * (lambda_ice-lambda_air) ! Jordan (1991)
+! ***********************************************************************************************************
+! public subroutine tcond_snow: compute thermal conductivity of snow
+! ***********************************************************************************************************
+subroutine tcond_snow(BulkDenIce,thermlcond,err,message)
+ implicit none
+ real(rkind),intent(in) :: BulkDenIce ! bulk density of ice (kg m-3)
+ real(rkind),intent(out) :: thermlcond ! thermal conductivity of snow (W m-1 K-1)
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! initialize error control
+ err=0; message="tcond_snow/"
+ ! compute thermal conductivity of snow
+ select case(model_decisions(iLookDECISIONS%thCondSnow)%iDecision)
+ case(Yen1965); thermlcond = 3.217d-6 * BulkDenIce**2._rkind ! Yen (1965)
+ case(Mellor1977); thermlcond = 2.576d-6 * BulkDenIce**2._rkind + 7.4d-2 ! Mellor (1977)
+ case(Jordan1991); thermlcond = lambda_air + (7.75d-5*BulkDenIce + 1.105d-6*(BulkDenIce**2._rkind)) &
+ * (lambda_ice-lambda_air) ! Jordan (1991)
case default
- err=10; message=trim(message)//"unknownOption"; return
- end select
- end subroutine tcond_snow
+ err=10; message=trim(message)//"unknownOption"; return
+ end select
+end subroutine tcond_snow
end module snow_utils_module
diff --git a/build/source/engine/soil_utils.f90 b/build/source/engine/soil_utils.f90
index 747618e..baab0fb 100755
--- a/build/source/engine/soil_utils.f90
+++ b/build/source/engine/soil_utils.f90
@@ -24,9 +24,9 @@ module soil_utils_module
USE nrtype
USE multiconst,only: gravity, & ! acceleration of gravity (m s-2)
- Tfreeze, & ! temperature at freezing (K)
- LH_fus, & ! latent heat of fusion (J kg-1, or m2 s-2)
- R_wv ! gas constant for water vapor (J kg-1 K-1; [J = Pa m3])
+ Tfreeze, & ! temperature at freezing (K)
+ LH_fus, & ! latent heat of fusion (J kg-1, or m2 s-2)
+ R_wv ! gas constant for water vapor (J kg-1 K-1; [J = Pa m3])
! privacy
implicit none
@@ -52,666 +52,664 @@ public::liquidHead
public::gammp
! constant parameters
-real(dp),parameter :: valueMissing=-9999._dp ! missing value parameter
-real(dp),parameter :: verySmall=epsilon(1.0_dp) ! a very small number (used to avoid divide by zero)
-real(dp),parameter :: dx=-1.e-12_dp ! finite difference increment
+real(rkind),parameter :: valueMissing=-9999._rkind ! missing value parameter
+real(rkind),parameter :: verySmall=epsilon(1.0_rkind) ! a very small number (used to avoid divide by zero)
+real(rkind),parameter :: dx=-1.e-12_rkind ! finite difference increment
contains
- ! ******************************************************************************************************************************
- ! public subroutine iceImpede: compute the ice impedence factor
- ! ******************************************************************************************************************************
- subroutine iceImpede(volFracIce,f_impede, & ! input
- iceImpedeFactor,dIceImpede_dLiq) ! output
- ! computes the ice impedence factor (separate function, as used multiple times)
- implicit none
- ! input variables
- real(dp),intent(in) :: volFracIce ! volumetric fraction of ice (-)
- real(dp),intent(in) :: f_impede ! ice impedence parameter (-)
- ! output variables
- real(dp) :: iceImpedeFactor ! ice impedence factor (-)
- real(dp) :: dIceImpede_dLiq ! derivative in ice impedence factor w.r.t. volumetric liquid water content (-)
- ! compute ice impedance factor as a function of volumetric ice content
- iceImpedeFactor = 10._dp**(-f_impede*volFracIce)
- dIceImpede_dLiq = 0._dp
-
- end subroutine iceImpede
-
-
- ! ******************************************************************************************************************************
- ! public subroutine dIceImpede_dTemp: compute the derivative in the ice impedence factor w.r.t. temperature
- ! ******************************************************************************************************************************
- subroutine dIceImpede_dTemp(volFracIce,dTheta_dT,f_impede,dIceImpede_dT)
- ! computes the derivative in the ice impedance factor w.r.t. temperature
- implicit none
- ! input variables
- real(dp),intent(in) :: volFracIce ! volumetric fraction of ice (-)
- real(dp),intent(in) :: dTheta_dT ! derivative in volumetric liquid water content w.r.t temperature (K-1)
- real(dp),intent(in) :: f_impede ! ice impedence parameter (-)
- ! output variables
- real(dp) :: dIceImpede_dT ! derivative in the ice impedance factor w.r.t. temperature (K-1)
- ! --
- dIceImpede_dT = log(10._dp)*f_impede*(10._dp**(-f_impede*volFracIce))*dTheta_dT
- end subroutine dIceImpede_dTemp
-
-
- ! ******************************************************************************************************************************
- ! public subroutine: compute the liquid water matric potential (and the derivatives w.r.t. total matric potential and temperature)
- ! ******************************************************************************************************************************
- subroutine liquidHead(&
- ! input
- matricHeadTotal ,& ! intent(in) : total water matric potential (m)
- volFracLiq ,& ! intent(in) : volumetric fraction of liquid water (-)
- volFracIce ,& ! intent(in) : volumetric fraction of ice (-)
- vGn_alpha,vGn_n,theta_sat,theta_res,vGn_m,& ! intent(in) : soil parameters
- dVolTot_dPsi0 ,& ! intent(in) : derivative in the soil water characteristic (m-1)
- dTheta_dT ,& ! intent(in) : derivative in volumetric total water w.r.t. temperature (K-1)
- ! output
- matricHeadLiq ,& ! intent(out) : liquid water matric potential (m)
- dPsiLiq_dPsi0 ,& ! intent(out) : derivative in the liquid water matric potential w.r.t. the total water matric potential (-)
- dPsiLiq_dTemp ,& ! intent(out) : derivative in the liquid water matric potential w.r.t. temperature (m K-1)
- err,message) ! intent(out) : error control
- ! computes the liquid water matric potential (and the derivatives w.r.t. total matric potential and temperature)
- implicit none
- ! input
- real(dp),intent(in) :: matricHeadTotal ! total water matric potential (m)
- real(dp),intent(in) :: volFracLiq ! volumetric fraction of liquid water (-)
- real(dp),intent(in) :: volFracIce ! volumetric fraction of ice (-)
- real(dp),intent(in) :: vGn_alpha,vGn_n,theta_sat,theta_res,vGn_m ! soil parameters
- real(dp),intent(in) ,optional :: dVolTot_dPsi0 ! derivative in the soil water characteristic (m-1)
- real(dp),intent(in) ,optional :: dTheta_dT ! derivative in volumetric total water w.r.t. temperature (K-1)
- ! output
- real(dp),intent(out) :: matricHeadLiq ! liquid water matric potential (m)
- real(dp),intent(out) ,optional :: dPsiLiq_dPsi0 ! derivative in the liquid water matric potential w.r.t. the total water matric potential (-)
- real(dp),intent(out) ,optional :: dPsiLiq_dTemp ! derivative in the liquid water matric potential w.r.t. temperature (m K-1)
- ! output: error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! local
- real(dp) :: xNum,xDen ! temporary variables (numeratir, denominator)
- real(dp) :: effSat ! effective saturation (-)
- real(dp) :: dPsiLiq_dEffSat ! derivative in liquid water matric potential w.r.t. effective saturation (m)
- real(dp) :: dEffSat_dTemp ! derivative in effective saturation w.r.t. temperature (K-1)
- ! ------------------------------------------------------------------------------------------------------------------------------
- ! initialize error control
- err=0; message='liquidHead/'
-
- ! ** partially frozen soil
- if(volFracIce > verySmall .and. matricHeadTotal < 0._dp)then ! check that ice exists and that the soil is unsaturated
-
- ! -----
- ! - compute liquid water matric potential...
- ! ------------------------------------------
-
- ! - compute effective saturation
- ! NOTE: include ice content as part of the solid porosity - major effect of ice is to reduce the pore size; ensure that effSat=1 at saturation
- ! (from Zhao et al., J. Hydrol., 1997: Numerical analysis of simultaneous heat and mass transfer...)
- xNum = volFracLiq - theta_res
- xDen = theta_sat - volFracIce - theta_res
- effSat = xNum/xDen ! effective saturation
-
- ! - matric head associated with liquid water
- matricHeadLiq = matricHead(effSat,vGn_alpha,0._dp,1._dp,vGn_n,vGn_m) ! argument is effective saturation, so theta_res=0 and theta_sat=1
-
- ! compute derivative in liquid water matric potential w.r.t. effective saturation (m)
- if(present(dPsiLiq_dPsi0).or.present(dPsiLiq_dTemp))then
- dPsiLiq_dEffSat = dPsi_dTheta(effSat,vGn_alpha,0._dp,1._dp,vGn_n,vGn_m)
- endif
-
- ! -----
- ! - compute derivative in the liquid water matric potential w.r.t. the total water matric potential...
- ! ----------------------------------------------------------------------------------------------------
-
- ! check if the derivative is desired
- if(present(dPsiLiq_dTemp))then
-
- ! (check required input derivative is present)
- if(.not.present(dVolTot_dPsi0))then
- message=trim(message)//'dVolTot_dPsi0 argument is missing'
- err=20; return
- endif
-
- ! (compute derivative in the liquid water matric potential w.r.t. the total water matric potential)
- dPsiLiq_dPsi0 = dVolTot_dPsi0*dPsiLiq_dEffSat*xNum/(xDen**2._dp)
-
- endif ! if dPsiLiq_dTemp is desired
-
- ! -----
- ! - compute the derivative in the liquid water matric potential w.r.t. temperature...
- ! -----------------------------------------------------------------------------------
-
- ! check if the derivative is desired
- if(present(dPsiLiq_dTemp))then
-
- ! (check required input derivative is present)
- if(.not.present(dTheta_dT))then
- message=trim(message)//'dTheta_dT argument is missing'
- err=20; return
- endif
-
- ! (compute the derivative in the liquid water matric potential w.r.t. temperature)
- dEffSat_dTemp = -dTheta_dT*xNum/(xDen**2._dp) + dTheta_dT/xDen
- dPsiLiq_dTemp = dPsiLiq_dEffSat*dEffSat_dTemp
-
- endif ! if dPsiLiq_dTemp is desired
-
- ! ** unfrozen soil
- else ! (no ice)
- matricHeadLiq = matricHeadTotal
- if(present(dPsiLiq_dTemp)) dPsiLiq_dPsi0 = 1._dp ! derivative=1 because values are identical
- if(present(dPsiLiq_dTemp)) dPsiLiq_dTemp = 0._dp ! derivative=0 because no impact of temperature for unfrozen conditions
- end if ! (if ice exists)
-
- end subroutine liquidHead
-
- ! ******************************************************************************************************************************
- ! public function hydCondMP_liq: compute the hydraulic conductivity of macropores as a function of liquid water content (m s-1)
- ! ******************************************************************************************************************************
- function hydCondMP_liq(volFracLiq,theta_sat,theta_mp,mpExp,satHydCond_ma,satHydCond_mi)
- ! computes hydraulic conductivity given volFracLiq and soil hydraulic parameters
- ! theta_sat, theta_mp, mpExp, satHydCond_ma, and satHydCond_mi
- implicit none
- ! dummies
- real(dp),intent(in) :: volFracLiq ! volumetric liquid water content (-)
- real(dp),intent(in) :: theta_sat ! soil porosity (-)
- real(dp),intent(in) :: theta_mp ! minimum volumetric liquid water content for macropore flow (-)
- real(dp),intent(in) :: mpExp ! empirical exponent in macropore flow equation (-)
- real(dp),intent(in) :: satHydCond_ma ! saturated hydraulic conductivity for macropores (m s-1)
- real(dp),intent(in) :: satHydCond_mi ! saturated hydraulic conductivity for micropores (m s-1)
- real(dp) :: hydCondMP_liq ! hydraulic conductivity (m s-1)
- ! locals
- real(dp) :: theta_e ! effective soil moisture
- if(volFracLiq > theta_mp)then
+! ******************************************************************************************************************************
+! public subroutine iceImpede: compute the ice impedence factor
+! ******************************************************************************************************************************
+subroutine iceImpede(volFracIce,f_impede, & ! input
+ iceImpedeFactor,dIceImpede_dLiq) ! output
+ ! computes the ice impedence factor (separate function, as used multiple times)
+ implicit none
+ ! input variables
+ real(rkind),intent(in) :: volFracIce ! volumetric fraction of ice (-)
+ real(rkind),intent(in) :: f_impede ! ice impedence parameter (-)
+ ! output variables
+ real(rkind) :: iceImpedeFactor ! ice impedence factor (-)
+ real(rkind) :: dIceImpede_dLiq ! derivative in ice impedence factor w.r.t. volumetric liquid water content (-)
+ ! compute ice impedance factor as a function of volumetric ice content
+ iceImpedeFactor = 10._rkind**(-f_impede*volFracIce)
+ dIceImpede_dLiq = 0._rkind
+
+end subroutine iceImpede
+
+
+! ******************************************************************************************************************************
+! public subroutine dIceImpede_dTemp: compute the derivative in the ice impedence factor w.r.t. temperature
+! ******************************************************************************************************************************
+subroutine dIceImpede_dTemp(volFracIce,dTheta_dT,f_impede,dIceImpede_dT)
+ ! computes the derivative in the ice impedance factor w.r.t. temperature
+ implicit none
+ ! input variables
+ real(rkind),intent(in) :: volFracIce ! volumetric fraction of ice (-)
+ real(rkind),intent(in) :: dTheta_dT ! derivative in volumetric liquid water content w.r.t temperature (K-1)
+ real(rkind),intent(in) :: f_impede ! ice impedence parameter (-)
+ ! output variables
+ real(rkind) :: dIceImpede_dT ! derivative in the ice impedance factor w.r.t. temperature (K-1)
+ ! --
+ dIceImpede_dT = log(10._rkind)*f_impede*(10._rkind**(-f_impede*volFracIce))*dTheta_dT
+end subroutine dIceImpede_dTemp
+
+
+! ******************************************************************************************************************************
+! public subroutine: compute the liquid water matric potential (and the derivatives w.r.t. total matric potential and temperature)
+! ******************************************************************************************************************************
+subroutine liquidHead(&
+ ! input
+ matricHeadTotal ,& ! intent(in) : total water matric potential (m)
+ volFracLiq ,& ! intent(in) : volumetric fraction of liquid water (-)
+ volFracIce ,& ! intent(in) : volumetric fraction of ice (-)
+ vGn_alpha,vGn_n,theta_sat,theta_res,vGn_m,& ! intent(in) : soil parameters
+ dVolTot_dPsi0 ,& ! intent(in) : derivative in the soil water characteristic (m-1)
+ dTheta_dT ,& ! intent(in) : derivative in volumetric total water w.r.t. temperature (K-1)
+ ! output
+ matricHeadLiq ,& ! intent(out) : liquid water matric potential (m)
+ dPsiLiq_dPsi0 ,& ! intent(out) : derivative in the liquid water matric potential w.r.t. the total water matric potential (-)
+ dPsiLiq_dTemp ,& ! intent(out) : derivative in the liquid water matric potential w.r.t. temperature (m K-1)
+ err,message) ! intent(out) : error control
+ ! computes the liquid water matric potential (and the derivatives w.r.t. total matric potential and temperature)
+ implicit none
+ ! input
+ real(rkind),intent(in) :: matricHeadTotal ! total water matric potential (m)
+ real(rkind),intent(in) :: volFracLiq ! volumetric fraction of liquid water (-)
+ real(rkind),intent(in) :: volFracIce ! volumetric fraction of ice (-)
+ real(rkind),intent(in) :: vGn_alpha,vGn_n,theta_sat,theta_res,vGn_m ! soil parameters
+ real(rkind),intent(in) ,optional :: dVolTot_dPsi0 ! derivative in the soil water characteristic (m-1)
+ real(rkind),intent(in) ,optional :: dTheta_dT ! derivative in volumetric total water w.r.t. temperature (K-1)
+ ! output
+ real(rkind),intent(out) :: matricHeadLiq ! liquid water matric potential (m)
+ real(rkind),intent(out) ,optional :: dPsiLiq_dPsi0 ! derivative in the liquid water matric potential w.r.t. the total water matric potential (-)
+ real(rkind),intent(out) ,optional :: dPsiLiq_dTemp ! derivative in the liquid water matric potential w.r.t. temperature (m K-1)
+ ! output: error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! local
+ real(rkind) :: xNum,xDen ! temporary variables (numeratir, denominator)
+ real(rkind) :: effSat ! effective saturation (-)
+ real(rkind) :: dPsiLiq_dEffSat ! derivative in liquid water matric potential w.r.t. effective saturation (m)
+ real(rkind) :: dEffSat_dTemp ! derivative in effective saturation w.r.t. temperature (K-1)
+ ! ------------------------------------------------------------------------------------------------------------------------------
+ ! initialize error control
+ err=0; message='liquidHead/'
+
+ ! ** partially frozen soil
+ if(volFracIce > verySmall .and. matricHeadTotal < 0._rkind)then ! check that ice exists and that the soil is unsaturated
+
+ ! -----
+ ! - compute liquid water matric potential...
+ ! ------------------------------------------
+
+ ! - compute effective saturation
+ ! NOTE: include ice content as part of the solid porosity - major effect of ice is to reduce the pore size; ensure that effSat=1 at saturation
+ ! (from Zhao et al., J. Hydrol., 1997: Numerical analysis of simultaneous heat and mass transfer...)
+ xNum = volFracLiq - theta_res
+ xDen = theta_sat - volFracIce - theta_res
+ effSat = xNum/xDen ! effective saturation
+
+ ! - matric head associated with liquid water
+ matricHeadLiq = matricHead(effSat,vGn_alpha,0._rkind,1._rkind,vGn_n,vGn_m) ! argument is effective saturation, so theta_res=0 and theta_sat=1
+
+ ! compute derivative in liquid water matric potential w.r.t. effective saturation (m)
+ if(present(dPsiLiq_dPsi0).or.present(dPsiLiq_dTemp))then
+ dPsiLiq_dEffSat = dPsi_dTheta(effSat,vGn_alpha,0._rkind,1._rkind,vGn_n,vGn_m)
+ endif
+
+ ! -----
+ ! - compute derivative in the liquid water matric potential w.r.t. the total water matric potential...
+ ! ----------------------------------------------------------------------------------------------------
+
+ ! check if the derivative is desired
+ if(present(dPsiLiq_dTemp))then
+
+ ! (check required input derivative is present)
+ if(.not.present(dVolTot_dPsi0))then
+ message=trim(message)//'dVolTot_dPsi0 argument is missing'
+ err=20; return
+ endif
+
+ ! (compute derivative in the liquid water matric potential w.r.t. the total water matric potential)
+ dPsiLiq_dPsi0 = dVolTot_dPsi0*dPsiLiq_dEffSat*xNum/(xDen**2._rkind)
+
+ endif ! if dPsiLiq_dTemp is desired
+
+ ! -----
+ ! - compute the derivative in the liquid water matric potential w.r.t. temperature...
+ ! -----------------------------------------------------------------------------------
+
+ ! check if the derivative is desired
+ if(present(dPsiLiq_dTemp))then
+
+ ! (check required input derivative is present)
+ if(.not.present(dTheta_dT))then
+ message=trim(message)//'dTheta_dT argument is missing'
+ err=20; return
+ endif
+
+ ! (compute the derivative in the liquid water matric potential w.r.t. temperature)
+ dEffSat_dTemp = -dTheta_dT*xNum/(xDen**2._rkind) + dTheta_dT/xDen
+ dPsiLiq_dTemp = dPsiLiq_dEffSat*dEffSat_dTemp
+
+ endif ! if dPsiLiq_dTemp is desired
+
+ ! ** unfrozen soil
+ else ! (no ice)
+ matricHeadLiq = matricHeadTotal
+ if(present(dPsiLiq_dTemp)) dPsiLiq_dPsi0 = 1._rkind ! derivative=1 because values are identical
+ if(present(dPsiLiq_dTemp)) dPsiLiq_dTemp = 0._rkind ! derivative=0 because no impact of temperature for unfrozen conditions
+ end if ! (if ice exists)
+
+end subroutine liquidHead
+
+! ******************************************************************************************************************************
+! public function hydCondMP_liq: compute the hydraulic conductivity of macropores as a function of liquid water content (m s-1)
+! ******************************************************************************************************************************
+function hydCondMP_liq(volFracLiq,theta_sat,theta_mp,mpExp,satHydCond_ma,satHydCond_mi)
+ ! computes hydraulic conductivity given volFracLiq and soil hydraulic parameters
+ ! theta_sat, theta_mp, mpExp, satHydCond_ma, and satHydCond_mi
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: volFracLiq ! volumetric liquid water content (-)
+ real(rkind),intent(in) :: theta_sat ! soil porosity (-)
+ real(rkind),intent(in) :: theta_mp ! minimum volumetric liquid water content for macropore flow (-)
+ real(rkind),intent(in) :: mpExp ! empirical exponent in macropore flow equation (-)
+ real(rkind),intent(in) :: satHydCond_ma ! saturated hydraulic conductivity for macropores (m s-1)
+ real(rkind),intent(in) :: satHydCond_mi ! saturated hydraulic conductivity for micropores (m s-1)
+ real(rkind) :: hydCondMP_liq ! hydraulic conductivity (m s-1)
+ ! locals
+ real(rkind) :: theta_e ! effective soil moisture
+ if(volFracLiq > theta_mp)then
theta_e = (volFracLiq - theta_mp) / (theta_sat - theta_mp)
hydCondMP_liq = (satHydCond_ma - satHydCond_mi) * (theta_e**mpExp)
- else
- hydCondMP_liq = 0._dp
- end if
- !write(*,'(a,4(f9.3,1x),2(e20.10))') 'in soil_utils: theta_mp, theta_sat, volFracLiq, hydCondMP_liq, satHydCond_ma, satHydCond_mi = ', &
- ! theta_mp, theta_sat, volFracLiq, hydCondMP_liq, satHydCond_ma, satHydCond_mi
- end function hydCondMP_liq
-
-
- ! ******************************************************************************************************************************
- ! public function hydCond_psi: compute the hydraulic conductivity as a function of matric head (m s-1)
- ! ******************************************************************************************************************************
- function hydCond_psi(psi,k_sat,alpha,n,m)
- ! computes hydraulic conductivity given psi and soil hydraulic parameters k_sat, alpha, n, and m
- implicit none
- ! dummies
- real(dp),intent(in) :: psi ! soil water suction (m)
- real(dp),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: hydCond_psi ! hydraulic conductivity (m s-1)
- if(psi<0._dp)then
+ else
+ hydCondMP_liq = 0._rkind
+ end if
+end function hydCondMP_liq
+
+
+! ******************************************************************************************************************************
+! public function hydCond_psi: compute the hydraulic conductivity as a function of matric head (m s-1)
+! ******************************************************************************************************************************
+function hydCond_psi(psi,k_sat,alpha,n,m)
+ ! computes hydraulic conductivity given psi and soil hydraulic parameters k_sat, alpha, n, and m
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: psi ! soil water suction (m)
+ real(rkind),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: hydCond_psi ! hydraulic conductivity (m s-1)
+ if(psi<0._rkind)then
hydCond_psi = k_sat * &
- ( ( (1._dp - (psi*alpha)**(n-1._dp) * (1._dp + (psi*alpha)**n)**(-m))**2._dp ) &
- / ( (1._dp + (psi*alpha)**n)**(m/2._dp) ) )
- else
+ ( ( (1._rkind - (psi*alpha)**(n-1._rkind) * (1._rkind + (psi*alpha)**n)**(-m))**2._rkind ) &
+ / ( (1._rkind + (psi*alpha)**n)**(m/2._rkind) ) )
+ else
hydCond_psi = k_sat
- end if
- end function hydCond_psi
-
-
- ! ******************************************************************************************************************************
- ! public function hydCond_liq: compute the hydraulic conductivity as a function of volumetric liquid water content (m s-1)
- ! ******************************************************************************************************************************
- function hydCond_liq(volFracLiq,k_sat,theta_res,theta_sat,m)
- ! computes hydraulic conductivity given volFracLiq and soil hydraulic parameters k_sat, theta_sat, theta_res, and m
- implicit none
- ! dummies
- real(dp),intent(in) :: volFracLiq ! volumetric liquid water content (-)
- real(dp),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
- real(dp),intent(in) :: theta_res ! residual volumetric liquid water content (-)
- real(dp),intent(in) :: theta_sat ! soil porosity (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: hydCond_liq ! hydraulic conductivity (m s-1)
- ! locals
- real(dp) :: theta_e ! effective soil moisture
- if(volFracLiq < theta_sat)then
+ end if
+end function hydCond_psi
+
+
+! ******************************************************************************************************************************
+! public function hydCond_liq: compute the hydraulic conductivity as a function of volumetric liquid water content (m s-1)
+! ******************************************************************************************************************************
+function hydCond_liq(volFracLiq,k_sat,theta_res,theta_sat,m)
+ ! computes hydraulic conductivity given volFracLiq and soil hydraulic parameters k_sat, theta_sat, theta_res, and m
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: volFracLiq ! volumetric liquid water content (-)
+ real(rkind),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric liquid water content (-)
+ real(rkind),intent(in) :: theta_sat ! soil porosity (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: hydCond_liq ! hydraulic conductivity (m s-1)
+ ! locals
+ real(rkind) :: theta_e ! effective soil moisture
+ if(volFracLiq < theta_sat)then
theta_e = (volFracLiq - theta_res) / (theta_sat - theta_res)
- hydCond_liq = k_sat*theta_e**(1._dp/2._dp) * (1._dp - (1._dp - theta_e**(1._dp/m) )**m)**2._dp
- else
+ hydCond_liq = k_sat*theta_e**(1._rkind/2._rkind) * (1._rkind - (1._rkind - theta_e**(1._rkind/m) )**m)**2._rkind
+ else
hydCond_liq = k_sat
- end if
- end function hydCond_liq
-
-
- ! ******************************************************************************************************************************
- ! public function volFracLiq: compute the volumetric liquid water content (-)
- ! ******************************************************************************************************************************
- function volFracLiq(psi,alpha,theta_res,theta_sat,n,m)
- ! computes the volumetric liquid water content given psi and soil hydraulic parameters theta_res, theta_sat, alpha, n, and m
- implicit none
- real(dp),intent(in) :: psi ! soil water suction (m)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: theta_res ! residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: volFracLiq ! volumetric liquid water content (-)
- if(psi<0._dp)then
- volFracLiq = theta_res + (theta_sat - theta_res)*(1._dp + (alpha*psi)**n)**(-m)
- else
+ end if
+end function hydCond_liq
+
+
+! ******************************************************************************************************************************
+! public function volFracLiq: compute the volumetric liquid water content (-)
+! ******************************************************************************************************************************
+function volFracLiq(psi,alpha,theta_res,theta_sat,n,m)
+ ! computes the volumetric liquid water content given psi and soil hydraulic parameters theta_res, theta_sat, alpha, n, and m
+ implicit none
+ real(rkind),intent(in) :: psi ! soil water suction (m)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: volFracLiq ! volumetric liquid water content (-)
+ if(psi<0._rkind)then
+ volFracLiq = theta_res + (theta_sat - theta_res)*(1._rkind + (alpha*psi)**n)**(-m)
+ else
volFracLiq = theta_sat
- end if
- end function volFracLiq
-
-
- ! ******************************************************************************************************************************
- ! public function matricHead: compute the matric head (m) based on the volumetric liquid water content
- ! ******************************************************************************************************************************
- function matricHead(theta,alpha,theta_res,theta_sat,n,m)
- ! computes the volumetric liquid water content given psi and soil hydraulic parameters theta_res, theta_sat, alpha, n, and m
- implicit none
- ! dummy variables
- real(dp),intent(in) :: theta ! volumetric liquid water content (-)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: theta_res ! residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: matricHead ! matric head (m)
- ! local variables
- real(dp) :: effSat ! effective saturation (-)
- real(dp),parameter :: verySmall=epsilon(1._dp) ! a very small number (avoid effective saturation of zero)
- ! compute effective saturation
- effSat = max(verySmall, (theta - theta_res) / (theta_sat - theta_res))
- ! compute matric head
- if (effSat < 1._dp .and. effSat > 0._dp)then
- matricHead = (1._dp/alpha)*( effSat**(-1._dp/m) - 1._dp)**(1._dp/n)
- else
- matricHead = 0._dp
- end if
- end function matricHead
-
-
- ! ******************************************************************************************************************************
- ! public function dTheta_dPsi: compute the derivative of the soil water characteristic (m-1)
- ! ******************************************************************************************************************************
- function dTheta_dPsi(psi,alpha,theta_res,theta_sat,n,m)
- implicit none
- real(dp),intent(in) :: psi ! soil water suction (m)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: theta_res ! residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: dTheta_dPsi ! derivative of the soil water characteristic (m-1)
- if(psi<=0._dp)then
+ end if
+end function volFracLiq
+
+
+! ******************************************************************************************************************************
+! public function matricHead: compute the matric head (m) based on the volumetric liquid water content
+! ******************************************************************************************************************************
+function matricHead(theta,alpha,theta_res,theta_sat,n,m)
+ ! computes the volumetric liquid water content given psi and soil hydraulic parameters theta_res, theta_sat, alpha, n, and m
+ implicit none
+ ! dummy variables
+ real(rkind),intent(in) :: theta ! volumetric liquid water content (-)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: matricHead ! matric head (m)
+ ! local variables
+ real(rkind) :: effSat ! effective saturation (-)
+ real(rkind),parameter :: verySmall=epsilon(1._rkind) ! a very small number (avoid effective saturation of zero)
+ ! compute effective saturation
+ effSat = max(verySmall, (theta - theta_res) / (theta_sat - theta_res))
+ ! compute matric head
+ if (effSat < 1._rkind .and. effSat > 0._rkind)then
+ matricHead = (1._rkind/alpha)*( effSat**(-1._rkind/m) - 1._rkind)**(1._rkind/n)
+ else
+ matricHead = 0._rkind
+ end if
+end function matricHead
+
+
+! ******************************************************************************************************************************
+! public function dTheta_dPsi: compute the derivative of the soil water characteristic (m-1)
+! ******************************************************************************************************************************
+function dTheta_dPsi(psi,alpha,theta_res,theta_sat,n,m)
+ implicit none
+ real(rkind),intent(in) :: psi ! soil water suction (m)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: dTheta_dPsi ! derivative of the soil water characteristic (m-1)
+ if(psi<=0._rkind)then
dTheta_dPsi = (theta_sat-theta_res) * &
- (-m*(1._dp + (psi*alpha)**n)**(-m-1._dp)) * n*(psi*alpha)**(n-1._dp) * alpha
+ (-m*(1._rkind + (psi*alpha)**n)**(-m-1._rkind)) * n*(psi*alpha)**(n-1._rkind) * alpha
if(abs(dTheta_dPsi) < epsilon(psi)) dTheta_dPsi = epsilon(psi)
- else
+ else
dTheta_dPsi = epsilon(psi)
- end if
- end function dTheta_dPsi
-
-
- ! ******************************************************************************************************************************
- ! public function dPsi_dTheta: compute the derivative of the soil water characteristic (m-1)
- ! ******************************************************************************************************************************
- function dPsi_dTheta(volFracLiq,alpha,theta_res,theta_sat,n,m)
- implicit none
- ! dummies
- real(dp),intent(in) :: volFracLiq ! volumetric liquid water content (-)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: theta_res ! residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: dPsi_dTheta ! derivative of the soil water characteristic (m)
- ! locals
- real(dp) :: y1,d1 ! 1st function and derivative
- real(dp) :: y2,d2 ! 2nd function and derivative
- real(dp) :: theta_e ! effective soil moisture
- ! check if less than saturation
- if(volFracLiq < theta_sat)then
+ end if
+end function dTheta_dPsi
+
+
+! ******************************************************************************************************************************
+! public function dPsi_dTheta: compute the derivative of the soil water characteristic (m-1)
+! ******************************************************************************************************************************
+function dPsi_dTheta(volFracLiq,alpha,theta_res,theta_sat,n,m)
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: volFracLiq ! volumetric liquid water content (-)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: dPsi_dTheta ! derivative of the soil water characteristic (m)
+ ! locals
+ real(rkind) :: y1,d1 ! 1st function and derivative
+ real(rkind) :: y2,d2 ! 2nd function and derivative
+ real(rkind) :: theta_e ! effective soil moisture
+ ! check if less than saturation
+ if(volFracLiq < theta_sat)then
! compute effective water content
theta_e = max(0.001,(volFracLiq - theta_res) / (theta_sat - theta_res))
! compute the 1st function and derivative
- y1 = theta_e**(-1._dp/m) - 1._dp
- d1 = (-1._dp/m)*theta_e**(-1._dp/m - 1._dp) / (theta_sat - theta_res)
+ y1 = theta_e**(-1._rkind/m) - 1._rkind
+ d1 = (-1._rkind/m)*theta_e**(-1._rkind/m - 1._rkind) / (theta_sat - theta_res)
! compute the 2nd function and derivative
- y2 = y1**(1._dp/n)
- d2 = (1._dp/n)*y1**(1._dp/n - 1._dp)
+ y2 = y1**(1._rkind/n)
+ d2 = (1._rkind/n)*y1**(1._rkind/n - 1._rkind)
! compute the final function value
dPsi_dTheta = d1*d2/alpha
- else
- dPsi_dTheta = 0._dp
- end if
- end function dPsi_dTheta
-
-
- ! ******************************************************************************************************************************
- ! public function dPsi_dTheta2: compute the derivative of dPsi_dTheta (m-1)
- ! ******************************************************************************************************************************
- function dPsi_dTheta2(volFracLiq,alpha,theta_res,theta_sat,n,m,lTangent)
- implicit none
- ! dummies
- real(dp),intent(in) :: volFracLiq ! volumetric liquid water content (-)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: theta_res ! residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
- real(dp) :: dPsi_dTheta2 ! derivative of the soil water characteristic (m)
- ! locals for analytical derivatives
- real(dp) :: xx ! temporary variable
- real(dp) :: y1,d1 ! 1st function and derivative
- real(dp) :: y2,d2 ! 2nd function and derivative
- real(dp) :: theta_e ! effective soil moisture
- ! locals for numerical derivative
- real(dp) :: func0,func1 ! function evaluations
- ! check if less than saturation
- if(volFracLiq < theta_sat)then
+ else
+ dPsi_dTheta = 0._rkind
+ end if
+end function dPsi_dTheta
+
+
+! ******************************************************************************************************************************
+! public function dPsi_dTheta2: compute the derivative of dPsi_dTheta (m-1)
+! ******************************************************************************************************************************
+function dPsi_dTheta2(volFracLiq,alpha,theta_res,theta_sat,n,m,lTangent)
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: volFracLiq ! volumetric liquid water content (-)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: theta_res ! residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
+ real(rkind) :: dPsi_dTheta2 ! derivative of the soil water characteristic (m)
+ ! locals for analytical derivatives
+ real(rkind) :: xx ! temporary variable
+ real(rkind) :: y1,d1 ! 1st function and derivative
+ real(rkind) :: y2,d2 ! 2nd function and derivative
+ real(rkind) :: theta_e ! effective soil moisture
+ ! locals for numerical derivative
+ real(rkind) :: func0,func1 ! function evaluations
+ ! check if less than saturation
+ if(volFracLiq < theta_sat)then
! ***** compute analytical derivatives
if(lTangent)then
- ! compute the effective saturation
- theta_e = (volFracLiq - theta_res) / (theta_sat - theta_res)
- ! get the first function and derivative
- y1 = (-1._dp/m)*theta_e**(-1._dp/m - 1._dp) / (theta_sat - theta_res)
- d1 = ( (m + 1._dp) / (m**2._dp * (theta_sat - theta_res)**2._dp) ) * theta_e**(-1._dp/m - 2._dp)
- ! get the second function and derivative
- xx = theta_e**(-1._dp/m) - 1._dp
- y2 = (1._dp/n)*xx**(1._dp/n - 1._dp)
- d2 = ( -(1._dp - n)/((theta_sat - theta_res)*m*n**2._dp) ) * xx**(1._dp/n - 2._dp) * theta_e**(-1._dp/m - 1._dp)
- ! return the derivative
- dPsi_dTheta2 = (d1*y2 + y1*d2)/alpha
+ ! compute the effective saturation
+ theta_e = (volFracLiq - theta_res) / (theta_sat - theta_res)
+ ! get the first function and derivative
+ y1 = (-1._rkind/m)*theta_e**(-1._rkind/m - 1._rkind) / (theta_sat - theta_res)
+ d1 = ( (m + 1._rkind) / (m**2._rkind * (theta_sat - theta_res)**2._rkind) ) * theta_e**(-1._rkind/m - 2._rkind)
+ ! get the second function and derivative
+ xx = theta_e**(-1._rkind/m) - 1._rkind
+ y2 = (1._rkind/n)*xx**(1._rkind/n - 1._rkind)
+ d2 = ( -(1._rkind - n)/((theta_sat - theta_res)*m*n**2._rkind) ) * xx**(1._rkind/n - 2._rkind) * theta_e**(-1._rkind/m - 1._rkind)
+ ! return the derivative
+ dPsi_dTheta2 = (d1*y2 + y1*d2)/alpha
! ***** compute numerical derivatives
else
- func0 = dPsi_dTheta(volFracLiq, alpha,theta_res,theta_sat,n,m)
- func1 = dPsi_dTheta(volFracLiq+dx,alpha,theta_res,theta_sat,n,m)
- dPsi_dTheta2 = (func1 - func0)/dx
+ func0 = dPsi_dTheta(volFracLiq, alpha,theta_res,theta_sat,n,m)
+ func1 = dPsi_dTheta(volFracLiq+dx,alpha,theta_res,theta_sat,n,m)
+ dPsi_dTheta2 = (func1 - func0)/dx
end if
- ! (case where volumetric liquid water content exceeds porosity)
- else
- dPsi_dTheta2 = 0._dp
- end if
- end function dPsi_dTheta2
-
-
- ! ******************************************************************************************************************************
- ! public function dHydCond_dPsi: compute the derivative in hydraulic conductivity w.r.t. matric head (s-1)
- ! ******************************************************************************************************************************
- function dHydCond_dPsi(psi,k_sat,alpha,n,m,lTangent)
- ! computes the derivative in hydraulic conductivity w.r.t matric head,
- ! given psi and soil hydraulic parameters k_sat, alpha, n, and m
- implicit none
- ! dummies
- real(dp),intent(in) :: psi ! soil water suction (m)
- real(dp),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
- real(dp),intent(in) :: alpha ! scaling parameter (m-1)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
- real(dp) :: dHydCond_dPsi ! derivative in hydraulic conductivity w.r.t. matric head (s-1)
- ! locals for analytical derivatives
- real(dp) :: f_x1 ! f(x) for part of the numerator
- real(dp) :: f_x2 ! f(x) for part of the numerator
- real(dp) :: f_nm ! f(x) for the numerator
- real(dp) :: f_dm ! f(x) for the denominator
- real(dp) :: d_x1 ! df(x)/dpsi for part of the numerator
- real(dp) :: d_x2 ! df(x)/dpsi for part of the numerator
- real(dp) :: d_nm ! df(x)/dpsi for the numerator
- real(dp) :: d_dm ! df(x)/dpsi for the denominator
- ! locals for numerical derivatives
- real(dp) :: hydCond0 ! hydraulic condictivity value for base case
- real(dp) :: hydCond1 ! hydraulic condictivity value for perturbed case
- ! derivative is zero if saturated
- if(psi<0._dp)then
+ ! (case where volumetric liquid water content exceeds porosity)
+ else
+ dPsi_dTheta2 = 0._rkind
+ end if
+end function dPsi_dTheta2
+
+
+! ******************************************************************************************************************************
+! public function dHydCond_dPsi: compute the derivative in hydraulic conductivity w.r.t. matric head (s-1)
+! ******************************************************************************************************************************
+function dHydCond_dPsi(psi,k_sat,alpha,n,m,lTangent)
+ ! computes the derivative in hydraulic conductivity w.r.t matric head,
+ ! given psi and soil hydraulic parameters k_sat, alpha, n, and m
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: psi ! soil water suction (m)
+ real(rkind),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
+ real(rkind),intent(in) :: alpha ! scaling parameter (m-1)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
+ real(rkind) :: dHydCond_dPsi ! derivative in hydraulic conductivity w.r.t. matric head (s-1)
+ ! locals for analytical derivatives
+ real(rkind) :: f_x1 ! f(x) for part of the numerator
+ real(rkind) :: f_x2 ! f(x) for part of the numerator
+ real(rkind) :: f_nm ! f(x) for the numerator
+ real(rkind) :: f_dm ! f(x) for the denominator
+ real(rkind) :: d_x1 ! df(x)/dpsi for part of the numerator
+ real(rkind) :: d_x2 ! df(x)/dpsi for part of the numerator
+ real(rkind) :: d_nm ! df(x)/dpsi for the numerator
+ real(rkind) :: d_dm ! df(x)/dpsi for the denominator
+ ! locals for numerical derivatives
+ real(rkind) :: hydCond0 ! hydraulic condictivity value for base case
+ real(rkind) :: hydCond1 ! hydraulic condictivity value for perturbed case
+ ! derivative is zero if saturated
+ if(psi<0._rkind)then
! ***** compute analytical derivatives
if(lTangent)then
- ! compute the derivative for the numerator
- f_x1 = (psi*alpha)**(n - 1._dp)
- f_x2 = (1._dp + (psi*alpha)**n)**(-m)
- d_x1 = alpha * (n - 1._dp)*(psi*alpha)**(n - 2._dp)
- d_x2 = alpha * n*(psi*alpha)**(n - 1._dp) * (-m)*(1._dp + (psi*alpha)**n)**(-m - 1._dp)
- f_nm = (1._dp - f_x1*f_x2)**2._dp
- d_nm = (-d_x1*f_x2 - f_x1*d_x2) * 2._dp*(1._dp - f_x1*f_x2)
- ! compute the derivative for the denominator
- f_dm = (1._dp + (psi*alpha)**n)**(m/2._dp)
- d_dm = alpha * n*(psi*alpha)**(n - 1._dp) * (m/2._dp)*(1._dp + (psi*alpha)**n)**(m/2._dp - 1._dp)
- ! and combine
- dHydCond_dPsi = k_sat*(d_nm*f_dm - d_dm*f_nm) / (f_dm**2._dp)
+ ! compute the derivative for the numerator
+ f_x1 = (psi*alpha)**(n - 1._rkind)
+ f_x2 = (1._rkind + (psi*alpha)**n)**(-m)
+ d_x1 = alpha * (n - 1._rkind)*(psi*alpha)**(n - 2._rkind)
+ d_x2 = alpha * n*(psi*alpha)**(n - 1._rkind) * (-m)*(1._rkind + (psi*alpha)**n)**(-m - 1._rkind)
+ f_nm = (1._rkind - f_x1*f_x2)**2._rkind
+ d_nm = (-d_x1*f_x2 - f_x1*d_x2) * 2._rkind*(1._rkind - f_x1*f_x2)
+ ! compute the derivative for the denominator
+ f_dm = (1._rkind + (psi*alpha)**n)**(m/2._rkind)
+ d_dm = alpha * n*(psi*alpha)**(n - 1._rkind) * (m/2._rkind)*(1._rkind + (psi*alpha)**n)**(m/2._rkind - 1._rkind)
+ ! and combine
+ dHydCond_dPsi = k_sat*(d_nm*f_dm - d_dm*f_nm) / (f_dm**2._rkind)
+ else
+ ! ***** compute numerical derivatives
+ hydcond0 = hydCond_psi(psi, k_sat,alpha,n,m)
+ hydcond1 = hydCond_psi(psi+dx,k_sat,alpha,n,m)
+ dHydCond_dPsi = (hydcond1 - hydcond0)/dx
+ end if
else
- ! ***** compute numerical derivatives
- hydcond0 = hydCond_psi(psi, k_sat,alpha,n,m)
- hydcond1 = hydCond_psi(psi+dx,k_sat,alpha,n,m)
- dHydCond_dPsi = (hydcond1 - hydcond0)/dx
+ dHydCond_dPsi = 0._rkind
end if
- else
- dHydCond_dPsi = 0._dp
- end if
- end function dHydCond_dPsi
-
-
- ! ******************************************************************************************************************************
- ! public function dHydCond_dLiq: compute the derivative in hydraulic conductivity w.r.t. volumetric liquid water content (m s-1)
- ! ******************************************************************************************************************************
- ! computes the derivative in hydraulic conductivity w.r.t the volumetric fraction of liquid water,
- ! given volFracLiq and soil hydraulic parameters k_sat, theta_sat, theta_res, and m
- ! ******************************************************************************************************************************
- function dHydCond_dLiq(volFracLiq,k_sat,theta_res,theta_sat,m,lTangent)
- implicit none
- ! dummies
- real(dp),intent(in) :: volFracLiq ! volumetric fraction of liquid water (-)
- real(dp),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
- real(dp),intent(in) :: theta_res ! soil residual volumetric water content (-)
- real(dp),intent(in) :: theta_sat ! soil porosity (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
- real(dp) :: dHydCond_dLiq ! derivative in hydraulic conductivity w.r.t. matric head (s-1)
- ! locals for analytical derivatives
- real(dp) :: theta_e ! effective soil moisture
- real(dp) :: f1 ! f(x) for the first function
- real(dp) :: d1 ! df(x)/dLiq for the first function
- real(dp) :: x1,x2 ! f(x) for different parts of the second function
- real(dp) :: p1,p2,p3 ! df(x)/dLiq for different parts of the second function
- real(dp) :: f2 ! f(x) for the second function
- real(dp) :: d2 ! df(x)/dLiq for the second function
- ! locals for numerical derivatives
- real(dp) :: hydCond0 ! hydraulic condictivity value for base case
- real(dp) :: hydCond1 ! hydraulic condictivity value for perturbed case
- ! derivative is zero if super-saturated
- if(volFracLiq < theta_sat)then
+end function dHydCond_dPsi
+
+
+! ******************************************************************************************************************************
+! public function dHydCond_dLiq: compute the derivative in hydraulic conductivity w.r.t. volumetric liquid water content (m s-1)
+! ******************************************************************************************************************************
+! computes the derivative in hydraulic conductivity w.r.t the volumetric fraction of liquid water,
+! given volFracLiq and soil hydraulic parameters k_sat, theta_sat, theta_res, and m
+! ******************************************************************************************************************************
+function dHydCond_dLiq(volFracLiq,k_sat,theta_res,theta_sat,m,lTangent)
+ implicit none
+ ! dummies
+ real(rkind),intent(in) :: volFracLiq ! volumetric fraction of liquid water (-)
+ real(rkind),intent(in) :: k_sat ! saturated hydraulic conductivity (m s-1)
+ real(rkind),intent(in) :: theta_res ! soil residual volumetric water content (-)
+ real(rkind),intent(in) :: theta_sat ! soil porosity (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ logical(lgt),intent(in) :: lTangent ! method used to compute derivative (.true. = analytical)
+ real(rkind) :: dHydCond_dLiq ! derivative in hydraulic conductivity w.r.t. matric head (s-1)
+ ! locals for analytical derivatives
+ real(rkind) :: theta_e ! effective soil moisture
+ real(rkind) :: f1 ! f(x) for the first function
+ real(rkind) :: d1 ! df(x)/dLiq for the first function
+ real(rkind) :: x1,x2 ! f(x) for different parts of the second function
+ real(rkind) :: p1,p2,p3 ! df(x)/dLiq for different parts of the second function
+ real(rkind) :: f2 ! f(x) for the second function
+ real(rkind) :: d2 ! df(x)/dLiq for the second function
+ ! locals for numerical derivatives
+ real(rkind) :: hydCond0 ! hydraulic condictivity value for base case
+ real(rkind) :: hydCond1 ! hydraulic condictivity value for perturbed case
+ ! derivative is zero if super-saturated
+ if(volFracLiq < theta_sat)then
! ***** compute analytical derivatives
if(lTangent)then
- ! compute the effective saturation
- theta_e = (volFracLiq - theta_res) / (theta_sat - theta_res)
- ! compute the function and derivative of the first fuction
- f1 = k_sat*theta_e**0.5_dp
- d1 = k_sat*0.5_dp*theta_e**(-0.5_dp) / (theta_sat - theta_res)
- ! compute the function and derivative of the second function
- ! (first part)
- x1 = 1._dp - theta_e**(1._dp/m)
- p1 = (-1._dp/m)*theta_e**(1._dp/m - 1._dp) / (theta_sat - theta_res) ! differentiate (1.d - theta_e**(1.d/m)
- ! (second part)
- x2 = x1**m
- p2 = m*x1**(m - 1._dp)
- ! (final)
- f2 = (1._dp - x2)**2._dp
- p3 = -2._dp*(1._dp - x2)
- ! (combine)
- d2 = p1*p2*p3
- ! pull it all together
- dHydCond_dLiq = (d1*f2 + d2*f1)
+ ! compute the effective saturation
+ theta_e = (volFracLiq - theta_res) / (theta_sat - theta_res)
+ ! compute the function and derivative of the first fuction
+ f1 = k_sat*theta_e**0.5_rkind
+ d1 = k_sat*0.5_rkind*theta_e**(-0.5_rkind) / (theta_sat - theta_res)
+ ! compute the function and derivative of the second function
+ ! (first part)
+ x1 = 1._rkind - theta_e**(1._rkind/m)
+ p1 = (-1._rkind/m)*theta_e**(1._rkind/m - 1._rkind) / (theta_sat - theta_res) ! differentiate (1.d - theta_e**(1.d/m)
+ ! (second part)
+ x2 = x1**m
+ p2 = m*x1**(m - 1._rkind)
+ ! (final)
+ f2 = (1._rkind - x2)**2._rkind
+ p3 = -2._rkind*(1._rkind - x2)
+ ! (combine)
+ d2 = p1*p2*p3
+ ! pull it all together
+ dHydCond_dLiq = (d1*f2 + d2*f1)
else
- ! ***** compute numerical derivatives
- hydcond0 = hydCond_liq(volFracLiq, k_sat,theta_res,theta_sat,m)
- hydcond1 = hydCond_liq(volFracLiq+dx,k_sat,theta_res,theta_sat,m)
- dHydCond_dLiq = (hydcond1 - hydcond0)/dx
+ ! ***** compute numerical derivatives
+ hydcond0 = hydCond_liq(volFracLiq, k_sat,theta_res,theta_sat,m)
+ hydcond1 = hydCond_liq(volFracLiq+dx,k_sat,theta_res,theta_sat,m)
+ dHydCond_dLiq = (hydcond1 - hydcond0)/dx
end if
- else
- dHydCond_dLiq = 0._dp
- end if
- end function dHydCond_dLiq
-
-
- ! ******************************************************************************************************************************
- ! public function RH_soilair: compute relative humidity of air in soil pore space
- ! ******************************************************************************************************************************
- function RH_soilair(matpot,Tk)
- implicit none
- real(dp),intent(in) :: matpot ! soil water suction -- matric potential (m)
- real(dp),intent(in) :: Tk ! temperature (K)
- real(dp) :: RH_soilair ! relative humidity of air in soil pore space
- ! compute relative humidity (UNITS NOTE: Pa = kg m-1 s-2, so R_wv units = m2 s-2 K-1)
- RH_soilair = exp( (gravity*matpot) / (R_wv*Tk) )
- end function RH_soilair
-
-
- ! ******************************************************************************************************************************
- ! public function crit_soilT: compute the critical temperature above which all water is unfrozen
- ! ******************************************************************************************************************************
- function crit_soilT(psi)
- implicit none
- real(dp),intent(in) :: psi ! matric head (m)
- real(dp) :: crit_soilT ! critical soil temperature (K)
- crit_soilT = Tfreeze + min(psi,0._dp)*gravity*Tfreeze/LH_fus
- end function crit_soilT
-
-
- ! ******************************************************************************************************************************
- ! public function dTheta_dTk: differentiate the freezing curve w.r.t. temperature
- ! ******************************************************************************************************************************
- function dTheta_dTk(Tk,theta_res,theta_sat,alpha,n,m)
- implicit none
- real(dp),intent(in) :: Tk ! temperature (K)
- real(dp),intent(in) :: theta_res ! residual liquid water content (-)
- real(dp),intent(in) :: theta_sat ! porosity (-)
- real(dp),intent(in) :: alpha ! vGn scaling parameter (m-1)
- real(dp),intent(in) :: n ! vGn "n" parameter (-)
- real(dp),intent(in) :: m ! vGn "m" parameter (-)
- real(dp) :: dTheta_dTk ! derivative of the freezing curve w.r.t. temperature (K-1)
- ! local variables
- real(dp) :: kappa ! constant (m K-1)
- real(dp) :: xtemp ! alpha*kappa*(Tk-Tfreeze) -- dimensionless variable (used more than once)
- ! compute kappa (m K-1)
- kappa = LH_fus/(gravity*Tfreeze) ! NOTE: J = kg m2 s-2
- ! define a tempory variable that is used more than once (-)
- xtemp = alpha*kappa*(Tk-Tfreeze)
- ! differentiate the freezing curve w.r.t. temperature -- making use of the chain rule
- dTheta_dTk = (alpha*kappa) * n*xtemp**(n - 1._dp) * (-m)*(1._dp + xtemp**n)**(-m - 1._dp) * (theta_sat - theta_res)
- end function dTheta_dTk
-
-
- ! ******************************************************************************************************************************
- ! public function gammp: compute cumulative probability using the Gamma distribution
- ! ******************************************************************************************************************************
- FUNCTION gammp(a,x)
- IMPLICIT NONE
- REAL(DP), INTENT(IN) :: a,x
- REAL(DP) :: gammp
- if (x<a+1.0_dp) then
+ else
+ dHydCond_dLiq = 0._rkind
+ end if
+end function dHydCond_dLiq
+
+
+! ******************************************************************************************************************************
+! public function RH_soilair: compute relative humidity of air in soil pore space
+! ******************************************************************************************************************************
+function RH_soilair(matpot,Tk)
+ implicit none
+ real(rkind),intent(in) :: matpot ! soil water suction -- matric potential (m)
+ real(rkind),intent(in) :: Tk ! temperature (K)
+ real(rkind) :: RH_soilair ! relative humidity of air in soil pore space
+ ! compute relative humidity (UNITS NOTE: Pa = kg m-1 s-2, so R_wv units = m2 s-2 K-1)
+ RH_soilair = exp( (gravity*matpot) / (R_wv*Tk) )
+end function RH_soilair
+
+
+! ******************************************************************************************************************************
+! public function crit_soilT: compute the critical temperature above which all water is unfrozen
+! ******************************************************************************************************************************
+function crit_soilT(psi)
+ implicit none
+ real(rkind),intent(in) :: psi ! matric head (m)
+ real(rkind) :: crit_soilT ! critical soil temperature (K)
+ crit_soilT = Tfreeze + min(psi,0._rkind)*gravity*Tfreeze/LH_fus
+end function crit_soilT
+
+
+! ******************************************************************************************************************************
+! public function dTheta_dTk: differentiate the freezing curve w.r.t. temperature
+! ******************************************************************************************************************************
+function dTheta_dTk(Tk,theta_res,theta_sat,alpha,n,m)
+ implicit none
+ real(rkind),intent(in) :: Tk ! temperature (K)
+ real(rkind),intent(in) :: theta_res ! residual liquid water content (-)
+ real(rkind),intent(in) :: theta_sat ! porosity (-)
+ real(rkind),intent(in) :: alpha ! vGn scaling parameter (m-1)
+ real(rkind),intent(in) :: n ! vGn "n" parameter (-)
+ real(rkind),intent(in) :: m ! vGn "m" parameter (-)
+ real(rkind) :: dTheta_dTk ! derivative of the freezing curve w.r.t. temperature (K-1)
+ ! local variables
+ real(rkind) :: kappa ! constant (m K-1)
+ real(rkind) :: xtemp ! alpha*kappa*(Tk-Tfreeze) -- dimensionless variable (used more than once)
+ ! compute kappa (m K-1)
+ kappa = LH_fus/(gravity*Tfreeze) ! NOTE: J = kg m2 s-2
+ ! define a tempory variable that is used more than once (-)
+ xtemp = alpha*kappa*(Tk-Tfreeze)
+ ! differentiate the freezing curve w.r.t. temperature -- making use of the chain rule
+ dTheta_dTk = (alpha*kappa) * n*xtemp**(n - 1._rkind) * (-m)*(1._rkind + xtemp**n)**(-m - 1._rkind) * (theta_sat - theta_res)
+end function dTheta_dTk
+
+
+! ******************************************************************************************************************************
+! public function gammp: compute cumulative probability using the Gamma distribution
+! ******************************************************************************************************************************
+FUNCTION gammp(a,x)
+ IMPLICIT NONE
+ real(rkind), INTENT(IN) :: a,x
+ real(rkind) :: gammp
+ if (x<a+1.0_rkind) then
gammp=gser(a,x)
- else
- gammp=1.0_dp-gcf(a,x)
- end if
- END FUNCTION gammp
-
-
- ! ******************************************************************************************************************************
- ! private function gcf: continued fraction development of the incomplete Gamma function
- ! ******************************************************************************************************************************
- FUNCTION gcf(a,x,gln)
- IMPLICIT NONE
- REAL(DP), INTENT(IN) :: a,x
- REAL(DP), OPTIONAL, INTENT(OUT) :: gln
- REAL(DP) :: gcf
- INTEGER(I4B), PARAMETER :: ITMAX=100
- REAL(DP), PARAMETER :: EPS=epsilon(x),FPMIN=tiny(x)/EPS
- INTEGER(I4B) :: i
- REAL(DP) :: an,b,c,d,del,h
- if (x == 0.0) then
+ else
+ gammp=1.0_rkind-gcf(a,x)
+ end if
+END FUNCTION gammp
+
+
+! ******************************************************************************************************************************
+! private function gcf: continued fraction development of the incomplete Gamma function
+! ******************************************************************************************************************************
+FUNCTION gcf(a,x,gln)
+ IMPLICIT NONE
+ real(rkind), INTENT(IN) :: a,x
+ real(rkind), OPTIONAL, INTENT(OUT) :: gln
+ real(rkind) :: gcf
+ INTEGER(I4B), PARAMETER :: ITMAX=100
+ real(rkind), PARAMETER :: EPS=epsilon(x),FPMIN=tiny(x)/EPS
+ INTEGER(I4B) :: i
+ real(rkind) :: an,b,c,d,del,h
+ if (x == 0.0) then
gcf=1.0
RETURN
- end if
- b=x+1.0_dp-a
- c=1.0_dp/FPMIN
- d=1.0_dp/b
- h=d
- do i=1,ITMAX
+ end if
+ b=x+1.0_rkind-a
+ c=1.0_rkind/FPMIN
+ d=1.0_rkind/b
+ h=d
+ do i=1,ITMAX
an=-i*(i-a)
- b=b+2.0_dp
+ b=b+2.0_rkind
d=an*d+b
if (abs(d) < FPMIN) d=FPMIN
c=b+an/c
if (abs(c) < FPMIN) c=FPMIN
- d=1.0_dp/d
+ d=1.0_rkind/d
del=d*c
h=h*del
- if (abs(del-1.0_dp) <= EPS) exit
- end do
- if (i > ITMAX) stop 'a too large, ITMAX too small in gcf'
- if (present(gln)) then
+ if (abs(del-1.0_rkind) <= EPS) exit
+ end do
+ if (i > ITMAX) stop 'a too large, ITMAX too small in gcf'
+ if (present(gln)) then
gln=gammln(a)
gcf=exp(-x+a*log(x)-gln)*h
- else
+ else
gcf=exp(-x+a*log(x)-gammln(a))*h
- end if
- END FUNCTION gcf
-
-
- ! ******************************************************************************************************************************
- ! private function gser: series development of the incomplete Gamma function
- ! ******************************************************************************************************************************
- FUNCTION gser(a,x,gln)
- IMPLICIT NONE
- REAL(DP), INTENT(IN) :: a,x
- REAL(DP), OPTIONAL, INTENT(OUT) :: gln
- REAL(DP) :: gser
- INTEGER(I4B), PARAMETER :: ITMAX=100
- REAL(DP), PARAMETER :: EPS=epsilon(x)
- INTEGER(I4B) :: n
- REAL(DP) :: ap,del,summ
- if (x == 0.0) then
+ end if
+END FUNCTION gcf
+
+
+! ******************************************************************************************************************************
+! private function gser: series development of the incomplete Gamma function
+! ******************************************************************************************************************************
+FUNCTION gser(a,x,gln)
+ IMPLICIT NONE
+ real(rkind), INTENT(IN) :: a,x
+ real(rkind), OPTIONAL, INTENT(OUT) :: gln
+ real(rkind) :: gser
+ INTEGER(I4B), PARAMETER :: ITMAX=100
+ real(rkind), PARAMETER :: EPS=epsilon(x)
+ INTEGER(I4B) :: n
+ real(rkind) :: ap,del,summ
+ if (x == 0.0) then
gser=0.0
RETURN
- end if
- ap=a
- summ=1.0_dp/a
- del=summ
- do n=1,ITMAX
- ap=ap+1.0_dp
+ end if
+ ap=a
+ summ=1.0_rkind/a
+ del=summ
+ do n=1,ITMAX
+ ap=ap+1.0_rkind
del=del*x/ap
summ=summ+del
if (abs(del) < abs(summ)*EPS) exit
- end do
- if (n > ITMAX) stop 'a too large, ITMAX too small in gser'
- if (present(gln)) then
+ end do
+ if (n > ITMAX) stop 'a too large, ITMAX too small in gser'
+ if (present(gln)) then
gln=gammln(a)
gser=summ*exp(-x+a*log(x)-gln)
- else
+ else
gser=summ*exp(-x+a*log(x)-gammln(a))
- end if
- END FUNCTION gser
-
-
- ! ******************************************************************************************************************************
- ! private function gammln: gamma function
- ! ******************************************************************************************************************************
- FUNCTION gammln(xx)
- USE nr_utility_module,only:arth ! use to build vectors with regular increments
- IMPLICIT NONE
- REAL(DP), INTENT(IN) :: xx
- REAL(DP) :: gammln
- REAL(DP) :: tmp,x
- REAL(DP) :: stp = 2.5066282746310005_dp
- REAL(DP), DIMENSION(6) :: coef = (/76.18009172947146_dp,&
- -86.50532032941677_dp,24.01409824083091_dp,&
- -1.231739572450155_dp,0.1208650973866179e-2_dp,&
- -0.5395239384953e-5_dp/)
- if(xx <= 0._dp) stop 'xx > 0 in gammln'
- x=xx
- tmp=x+5.5_dp
- tmp=(x+0.5_dp)*log(tmp)-tmp
- gammln=tmp+log(stp*(1.000000000190015_dp+&
- sum(coef(:)/arth(x+1.0_dp,1.0_dp,size(coef))))/x)
- END FUNCTION gammln
+ end if
+END FUNCTION gser
+
+
+! ******************************************************************************************************************************
+! private function gammln: gamma function
+! ******************************************************************************************************************************
+FUNCTION gammln(xx)
+ USE nr_utility_module,only:arth ! use to build vectors with regular increments
+ IMPLICIT NONE
+ real(rkind), INTENT(IN) :: xx
+ real(rkind) :: gammln
+ real(rkind) :: tmp,x
+ real(rkind) :: stp = 2.5066282746310005_rkind
+ real(rkind), DIMENSION(6) :: coef = (/76.18009172947146_rkind,&
+ -86.50532032941677_rkind,24.01409824083091_rkind,&
+ -1.231739572450155_rkind,0.1208650973866179e-2_rkind,&
+ -0.5395239384953e-5_rkind/)
+ if(xx <= 0._rkind) stop 'xx > 0 in gammln'
+ x=xx
+ tmp=x+5.5_rkind
+ tmp=(x+0.5_rkind)*log(tmp)-tmp
+ gammln=tmp+log(stp*(1.000000000190015_rkind+&
+ sum(coef(:)/arth(x+1.0_rkind,1.0_rkind,size(coef))))/x)
+END FUNCTION gammln
end module soil_utils_module
diff --git a/build/source/engine/updatState.f90 b/build/source/engine/updatState.f90
index 698c8b1..29c82f1 100755
--- a/build/source/engine/updatState.f90
+++ b/build/source/engine/updatState.f90
@@ -35,121 +35,127 @@ public::updateSoil
contains
- ! *************************************************************************************************************
- ! public subroutine updateSnow: compute phase change impacts on volumetric liquid water and ice
- ! *************************************************************************************************************
- subroutine updateSnow(&
- ! input
- mLayerTemp ,& ! intent(in): temperature (K)
- mLayerTheta ,& ! intent(in): volume fraction of total water (-)
- snowfrz_scale ,& ! intent(in): scaling parameter for the snow freezing curve (K-1)
- ! output
- mLayerVolFracLiq ,& ! intent(out): volumetric fraction of liquid water (-)
- mLayerVolFracIce ,& ! intent(out): volumetric fraction of ice (-)
- fLiq ,& ! intent(out): fraction of liquid water (-)
- err,message) ! intent(out): error control
- ! utility routines
- USE snow_utils_module,only:fracliquid ! compute volumetric fraction of liquid water
- implicit none
- ! input variables
- real(dp),intent(in) :: mLayerTemp ! temperature (K)
- real(dp),intent(in) :: mLayerTheta ! volume fraction of total water (-)
- real(dp),intent(in) :: snowfrz_scale ! scaling parameter for the snow freezing curve (K-1)
- ! output variables
- real(dp),intent(out) :: mLayerVolFracLiq ! volumetric fraction of liquid water (-)
- real(dp),intent(out) :: mLayerVolFracIce ! volumetric fraction of ice (-)
- real(dp),intent(out) :: fLiq ! fraction of liquid water (-)
- ! error control
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! initialize error control
- err=0; message="updateSnow/"
-
- ! compute the volumetric fraction of liquid water and ice (-)
- fLiq = fracliquid(mLayerTemp,snowfrz_scale)
- mLayerVolFracLiq = fLiq*mLayerTheta
- mLayerVolFracIce = (1._dp - fLiq)*mLayerTheta*(iden_water/iden_ice)
- !print*, 'mLayerTheta - (mLayerVolFracIce*(iden_ice/iden_water) + mLayerVolFracLiq) = ', mLayerTheta - (mLayerVolFracIce*(iden_ice/iden_water) + mLayerVolFracLiq)
- !write(*,'(a,1x,4(f20.10,1x))') 'in updateSnow: fLiq, mLayerTheta, mLayerVolFracIce = ', &
- ! fLiq, mLayerTheta, mLayerVolFracIce
- !pause
-
- end subroutine updateSnow
-
- ! *************************************************************************************************************
- ! public subroutine updateSoil: compute phase change impacts on matric head and volumetric liquid water and ice
- ! *************************************************************************************************************
- subroutine updateSoil(&
- ! input
- mLayerTemp ,& ! intent(in): temperature vector (K)
- mLayerMatricHead ,& ! intent(in): matric head (m)
- vGn_alpha ,& ! intent(in): van Genutchen "alpha" parameter
- vGn_n ,& ! intent(in): van Genutchen "n" parameter
- theta_sat ,& ! intent(in): soil porosity (-)
- theta_res ,& ! intent(in): soil residual volumetric water content (-)
- vGn_m ,& ! intent(in): van Genutchen "m" parameter (-)
- ! output
- mLayerVolFracWat ,& ! intent(out): volumetric fraction of total water (-)
- mLayerVolFracLiq ,& ! intent(out): volumetric fraction of liquid water (-)
- mLayerVolFracIce ,& ! intent(out): volumetric fraction of ice (-)
- err,message) ! intent(out): error control
- ! utility routines
- USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water based on matric head
- USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric liquid water content
- implicit none
- ! input variables
- real(dp),intent(in) :: mLayerTemp ! estimate of temperature (K)
- real(dp),intent(in) :: mLayerMatricHead ! matric head (m)
- real(dp),intent(in) :: vGn_alpha ! van Genutchen "alpha" parameter
- real(dp),intent(in) :: vGn_n ! van Genutchen "n" parameter
- real(dp),intent(in) :: theta_sat ! soil porosity (-)
- real(dp),intent(in) :: theta_res ! soil residual volumetric water content (-)
- real(dp),intent(in) :: vGn_m ! van Genutchen "m" parameter (-)
- ! output variables
- real(dp),intent(out) :: mLayerVolFracWat ! fractional volume of total water (-)
- real(dp),intent(out) :: mLayerVolFracLiq ! volumetric fraction of liquid water (-)
- real(dp),intent(out) :: mLayerVolFracIce ! volumetric fraction of ice (-)
- integer(i4b),intent(out) :: err ! error code
- character(*),intent(out) :: message ! error message
- ! define local variables
- real(dp) :: TcSoil ! critical soil temperature when all water is unfrozen (K)
- real(dp) :: xConst ! constant in the freezing curve function (m K-1)
- real(dp) :: mLayerPsiLiq ! liquid water matric potential (m)
- ! initialize error control
- err=0; message="updateSoil/"
-
- ! compute fractional **volume** of total water (liquid plus ice)
- mLayerVolFracWat = volFracLiq(mLayerMatricHead,vGn_alpha,theta_res,theta_sat,vGn_n,vGn_m)
- if(mLayerVolFracWat > theta_sat)then; err=20; message=trim(message)//'volume of liquid and ice exceeds porosity'; return; end if
-
- ! compute the critical soil temperature where all water is unfrozen (K)
- ! (eq 17 in Dall'Amico 2011)
- TcSoil = Tfreeze + min(mLayerMatricHead,0._dp)*gravity*Tfreeze/LH_fus ! (NOTE: J = kg m2 s-2, so LH_fus is in units of m2 s-2)
-
- ! *** compute volumetric fraction of liquid water and ice for partially frozen soil
- if(mLayerTemp < TcSoil)then ! (check if soil temperature is less than the critical temperature)
-
- ! - volumetric liquid water content (-)
- ! NOTE: mLayerPsiLiq is the liquid water matric potential from the Clapeyron equation, used to separate the total water into liquid water and ice
- ! mLayerPsiLiq is DIFFERENT from the liquid water matric potential used in the flux calculations
- xConst = LH_fus/(gravity*Tfreeze) ! m K-1 (NOTE: J = kg m2 s-2)
- mLayerPsiLiq = xConst*(mLayerTemp - Tfreeze) ! liquid water matric potential from the Clapeyron eqution
- mLayerVolFracLiq = volFracLiq(mLayerPsiLiq,vGn_alpha,theta_res,theta_sat,vGn_n,vGn_m)
-
- ! - volumetric ice content (-)
- mLayerVolFracIce = mLayerVolFracWat - mLayerVolFracLiq
-
- ! *** compute volumetric fraction of liquid water and ice for unfrozen soil
- else
-
- ! all water is unfrozen
- mLayerPsiLiq = mLayerMatricHead
- mLayerVolFracLiq = mLayerVolFracWat
- mLayerVolFracIce = 0._dp
-
- end if ! (check if soil is partially frozen)
-
- end subroutine updateSoil
-
-
-end module updatState_module
+! *************************************************************************************************************
+! public subroutine updateSnow: compute phase change impacts on volumetric liquid water and ice (veg or soil)
+! *************************************************************************************************************
+subroutine updateSnow(&
+ ! input
+ mLayerTemp ,& ! intent(in): temperature (K)
+ mLayerTheta ,& ! intent(in): volume fraction of total water (-)
+ snowfrz_scale ,& ! intent(in): scaling parameter for the snow freezing curve (K-1)
+ ! output
+ mLayerVolFracLiq ,& ! intent(out): volumetric fraction of liquid water (-)
+ mLayerVolFracIce ,& ! intent(out): volumetric fraction of ice (-)
+ fLiq ,& ! intent(out): fraction of liquid water (-)
+ err,message) ! intent(out): error control
+ ! utility routines
+ USE snow_utils_module,only:fracliquid ! compute volumetric fraction of liquid water
+ implicit none
+ ! input variables
+ real(rkind),intent(in) :: mLayerTemp ! temperature (K)
+ real(rkind),intent(in) :: mLayerTheta ! volume fraction of total water (-)
+ real(rkind),intent(in) :: snowfrz_scale ! scaling parameter for the snow freezing curve (K-1)
+ ! output variables
+ real(rkind),intent(out) :: mLayerVolFracLiq ! volumetric fraction of liquid water (-)
+ real(rkind),intent(out) :: mLayerVolFracIce ! volumetric fraction of ice (-)
+ real(rkind),intent(out) :: fLiq ! fraction of liquid water (-)
+ ! error control
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! initialize error control
+ err=0; message="updateSnow/"
+
+ ! compute the volumetric fraction of liquid water and ice (-)
+ fLiq = fracliquid(mLayerTemp,snowfrz_scale)
+ mLayerVolFracLiq = fLiq*mLayerTheta
+ mLayerVolFracIce = (1._rkind - fLiq)*mLayerTheta*(iden_water/iden_ice)
+end subroutine updateSnow
+
+! *************************************************************************************************************
+! public subroutine updateSoil: compute phase change impacts on matric head and volumetric liquid water and ice
+! *************************************************************************************************************
+subroutine updateSoil(&
+ ! input
+ mLayerTemp ,& ! intent(in): temperature vector (K)
+ mLayerMatricHead ,& ! intent(in): matric head (m)
+ vGn_alpha ,& ! intent(in): van Genutchen "alpha" parameter
+ vGn_n ,& ! intent(in): van Genutchen "n" parameter
+ theta_sat ,& ! intent(in): soil porosity (-)
+ theta_res ,& ! intent(in): soil residual volumetric water content (-)
+ vGn_m ,& ! intent(in): van Genutchen "m" parameter (-)
+ ! output
+ mLayerVolFracWat ,& ! intent(out): volumetric fraction of total water (-)
+ mLayerVolFracLiq ,& ! intent(out): volumetric fraction of liquid water (-)
+ mLayerVolFracIce ,& ! intent(out): volumetric fraction of ice (-)
+ err,message) ! intent(out): error control
+ ! utility routines
+ USE soil_utils_module,only:volFracLiq ! compute volumetric fraction of liquid water based on matric head
+ USE soil_utils_module,only:matricHead ! compute the matric head based on volumetric liquid water content
+ implicit none
+ ! input variables
+ real(rkind),intent(in) :: mLayerTemp ! estimate of temperature (K)
+ real(rkind),intent(in) :: mLayerMatricHead ! matric head (m)
+ real(rkind),intent(in) :: vGn_alpha ! van Genutchen "alpha" parameter
+ real(rkind),intent(in) :: vGn_n ! van Genutchen "n" parameter
+ real(rkind),intent(in) :: theta_sat ! soil porosity (-)
+ real(rkind),intent(in) :: theta_res ! soil residual volumetric water content (-)
+ real(rkind),intent(in) :: vGn_m ! van Genutchen "m" parameter (-)
+ ! output variables
+ real(rkind),intent(out) :: mLayerVolFracWat ! fractional volume of total water (-)
+ real(rkind),intent(out) :: mLayerVolFracLiq ! volumetric fraction of liquid water (-)
+ real(rkind),intent(out) :: mLayerVolFracIce ! volumetric fraction of ice (-)
+ integer(i4b),intent(out) :: err ! error code
+ character(*),intent(out) :: message ! error message
+ ! define local variables
+ real(rkind) :: TcSoil ! critical soil temperature when all water is unfrozen (K)
+ real(rkind) :: xConst ! constant in the freezing curve function (m K-1)
+ real(rkind) :: mLayerPsiLiq ! liquid water matric potential (m)
+ real(rkind),parameter :: tinyVal=epsilon(1._rkind) ! used in balance check
+ ! initialize error control
+ err=0; message="updateSoil/"
+
+ ! compute fractional **volume** of total water (liquid plus ice)
+ mLayerVolFracWat = volFracLiq(mLayerMatricHead,vGn_alpha,theta_res,theta_sat,vGn_n,vGn_m)
+ if(mLayerVolFracWat > (theta_sat + tinyVal)) then
+ err=20
+ message=trim(message)//'volume of liquid and ice (mLayerVolFracWat) exceeds porosity'
+ print*, 'mLayerVolFracWat = ', mLayerVolFracWat
+ print*, 'theta_sat (porosity) = ', theta_sat
+ print*, 'mLayerMatricHead = ', mLayerMatricHead
+ print*, 'theta_res = ', theta_res
+ print*, 'vGn_alpha = ', vGn_alpha
+ print*, 'vGn_n = ', vGn_n
+ print*, 'vGn_m = ', vGn_m
+ return
+ end if
+
+ ! compute the critical soil temperature where all water is unfrozen (K)
+ ! (eq 17 in Dall'Amico 2011)
+ TcSoil = Tfreeze + min(mLayerMatricHead,0._rkind)*gravity*Tfreeze/LH_fus ! (NOTE: J = kg m2 s-2, so LH_fus is in units of m2 s-2)
+
+ ! *** compute volumetric fraction of liquid water and ice for partially frozen soil
+ if(mLayerTemp < TcSoil)then ! (check if soil temperature is less than the critical temperature)
+
+ ! - volumetric liquid water content (-)
+ ! NOTE: mLayerPsiLiq is the liquid water matric potential from the Clapeyron equation, used to separate the total water into liquid water and ice
+ ! mLayerPsiLiq is DIFFERENT from the liquid water matric potential used in the flux calculations
+ xConst = LH_fus/(gravity*Tfreeze) ! m K-1 (NOTE: J = kg m2 s-2)
+ mLayerPsiLiq = xConst*(mLayerTemp - Tfreeze) ! liquid water matric potential from the Clapeyron eqution
+ mLayerVolFracLiq = volFracLiq(mLayerPsiLiq,vGn_alpha,theta_res,theta_sat,vGn_n,vGn_m)
+
+ ! - volumetric ice content (-)
+ mLayerVolFracIce = mLayerVolFracWat - mLayerVolFracLiq
+
+ ! *** compute volumetric fraction of liquid water and ice for unfrozen soil
+ else
+
+ ! all water is unfrozen
+ mLayerPsiLiq = mLayerMatricHead
+ mLayerVolFracLiq = mLayerVolFracWat
+ mLayerVolFracIce = 0._rkind
+
+ end if ! (check if soil is partially frozen)
+
+end subroutine updateSoil
+
+end module updatState_module
\ No newline at end of file
--
GitLab