diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml
index dfa66a6da949..93701f1ab005 100644
--- a/components/elm/bld/namelist_files/namelist_defaults.xml
+++ b/components/elm/bld/namelist_files/namelist_defaults.xml
@@ -113,6 +113,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
.false.
+
+.false.
+
ERMM
diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml
index e5c5e9a2bcf7..5d476e3ce718 100644
--- a/components/elm/bld/namelist_files/namelist_definition.xml
+++ b/components/elm/bld/namelist_files/namelist_definition.xml
@@ -2054,6 +2054,11 @@ Runtime flag to turn on/off PETSc based thermal model.
Runtime flag to turn on/off variable soil thickness.
+
+Runtime flag to turn on/off h3D.
+
+
diff --git a/components/elm/docs/figures/h3d_schematic.jpg b/components/elm/docs/figures/h3d_schematic.jpg
new file mode 100755
index 000000000000..e8171d688e69
Binary files /dev/null and b/components/elm/docs/figures/h3d_schematic.jpg differ
diff --git a/components/elm/docs/tech-guide/h3d.md b/components/elm/docs/tech-guide/h3d.md
new file mode 100644
index 000000000000..ee6b371bf4aa
--- /dev/null
+++ b/components/elm/docs/tech-guide/h3d.md
@@ -0,0 +1,391 @@
+# Hybrid-3D hillslope hydrological model
+
+## Overview
+
+H3D represents subsurface lateral flow and groundwater dynamics within an idealized hillslope unit inside each land
+grid cell. Rather than treating each soil column as an isolated vertical profile,
+h3D connects multiple columns along a topographic gradient to explicitly simulate downslope drainage,
+water table redistribution, and interaction with the channel.
+
+This approach is conceptually one-dimensional along the hillslope,
+but retains three-dimensional realism by including slope, width, and drainage area variation.
+The model solves the Dupuit–Boussinesq groundwater flow equation for saturated thickness using an implicit finite-difference method.
+
+In order to use this model, the following variables need to be added in the surface dataset:
+
+hs_x(begg:endg, 1:nh3dc_per_lunit+1): x‐coordinates (m) of the node/edge positions along the hillslope for each grid cell g. There are N+1 nodes for N columns.
+
+hs_w(begg:endg, 1:nh3dc_per_lunit+1): width function values (m) at those nodes (planform width vs. distance).
+
+hs_area(1:nh3dc_per_lunit): Planform area between nodes k and k+1 is calculated using a trapezoid, area each column are then normalized to the total area and used later for area-weighted averages/sums across the set of hillslope columns
+within a land unit (e.g., to compute landunit-level means of fluxes or states).
+
+## Conceptual Representation
+
+Each land unit is subdivided into several h3D columns positioned along an idealized hillslope:
+
+- Lower boundary corresponds to the stream outlet.
+- Upper boundary represents the local topographic divide.
+- Intermediate nodes represent soil columns along the slope.
+
+
+
+Figure 1. The hybrid-3D hillslope hydrological model represents subsurface and surface flow along an idealized hillslope. Left panel (plan view):
+A hillslope of total length $L$ is divided into $N$ vertical soil columns, each with an equal horizontal length $\Delta x$, underlain by bedrock with a topographic slope $\alpha$.
+The width of the hillslope at a given position is denoted by $W(x)$, measured from the lower boundary of the lowest column (column 0).
+The column index $j$ refers to the center of each hillslope column.
+Land-surface properties such as bedrock depth, vegetation cover, and atmospheric forcing are assumed to be identical for all columns.
+The variable $d$ represents the depth of overland flow above the surface. Right panel (vertical section):
+Each hillslope column consists of multiple soil layers of variable thickness $\Delta z_i$.
+Vertical flow through the unsaturated and saturated zones occurs across these layers with corresponding soil water potentials $\psi_i$.
+The flow domain extends down to an impermeable bedrock boundary, which imposes a zero-flux condition at the bottom.
+Vertical flow is solved using the $\alpha$-based form of the Richards equation (Zeng & Decker, 2009).
+The layer index $i$ denotes vertical soil layers, $h$ is the height of the saturation zone, and $\nabla$ marks the position of the water table.
+
+For each node, the state variable is the saturated thickness, $h_{sat}(x,t)$, measured from the bedrock to the water table.
+The model tracks how $h_{sat}(x,t)$ evolves in time due to:
+
+- Local recharge (vertical infiltration)
+- Downslope lateral flow driven by topographic gradients
+- Variable transmissivity and drainable porosity along the slope
+
+
+Each land unit represents a single hillslope, which has consistent topographic and geometric properties:
+- **Overall slope angle:** $\alpha$ — mean hillslope angle (rad)
+- **Width function:** $w(x)$ — lateral width distribution along the hillslope (m)
+- **Distance function:** $x(i)$ — distance from the stream outlet to node $i$ (m)
+- **Total hillslope area:** $A_{hs} = \displaystyle \int_0^L w(x)\ dx$ (m²)
+
+Each column is a point along hillslope, representing a cross-section at distance x along the hillslope. Each has local properties, connected by lateral flow to adjacent columns:
+- Soil hydraulic properties (K_sat, porosity, etc.)
+- Local area contribution (hs_dA)
+- Bedrock depth
+
+Land unit captures the hillsclope-scale connectivity, while columns capture spatial variation along the flow path. Each column contributes proportionally to the area within the land unit, conserves water mass at hillslope scale, having realistic flow convergence/divergence. Column 1 = stream/outlet, column N = hillslope divide.
+
+
+## Govergning Equation
+
+The fundamental PDE is a Dupuit-style Boussinesq groundwater flow equation for saturated flow along the slope:
+
+$$
+f\frac{\partial h}{\partial t} = \frac{1}{w}\frac{\partial}{\partial x}\left(w\ k_{l}(h)\ h\left(\sin(\alpha) + \frac{\partial h}{\partial x}\cos(\alpha)\right)\right) + \cos(\alpha)\ R
+$$
+
+where $h(x,t)$ is the saturated thickness [m],
+$f_{\text{drain}}$ is the drainable porosity [–],
+$\alpha$ is the hillslope angle [rad],
+w (m) is the width of the hillslope at a given distance x (m) from the outflow point, kl(h) (m/s) is the lateral saturated hydraulic conductivity at height h. $R$ is the recharge rate between the unsaturated and saturated zones (m/s). Note - recharge rate was not in the code.
+
+## Boundary Conditions
+
+| Boundary | Condition |
+| :------------- | :------------------------------------ |
+| Lower (stream) | $\partial h/\partial x = 0$ |
+| Upper (divide) | zero lateral flow |
+
+## Constitutive Relationships (from the code)
+
+### Transmissivity
+
+The local transmissivity depends on saturated thickness and anisotropy:
+
+$$
+T(x,h)= \frac{K_{\text{aniso}}K_{\text{sat}}(z_{wt}) \ h \ w(x)}{1000}
+$$
+
+where $K_{\text{aniso}}=100$ is the horizontal/vertical anisotropy factor,
+$K_{\text{sat}}(z_{wt})$ is the saturated hydraulic conductivity at the
+water-table depth [mm s⁻¹], and $w(x)$ is the hillslope width [m].
+Division by 1000 converts from mm s⁻¹ to m s⁻¹.
+
+### Variable Drainable Porosity
+
+The specific yield varies with depth following a Brooks–Corey relation:
+
+$$
+f_{\text{drain}}
+= \alpha_{\text{sat}}
+\left[
+1 -
+\left(
+1 +
+\frac{1000\ \max(0,z_{\text{bed}}-h)}{\psi_{\text{sat}}}
+\right)^{-1/b}
+\right],
+\qquad
+f_{\text{drain}}\ge0.02
+$$
+
+where $\alpha_{\text{sat}}$ is porosity [–],
+$z_{\text{bed}}$ is bedrock depth [m],
+$\psi_{\text{sat}}$ is air-entry suction [mm],
+and $b$ is the Brooks–Corey pore-size index [–].
+
+This function allows for a smooth transition between unsaturated and fully saturated conditions and ensures stability under variable soil thickness.
+
+
+## Numerical Implementation
+
+### Spatial Discretization
+
+The PDE is solved implicitly in space and time using a tridiagonal
+system for $h_i^{t}$ at each node $i$ at time t. Node is ordered from 1 to N:
+
+$$
+a_i h_{i-1}^{t,s+1} + b_i h_i^{t,s+1} + c_i h_{i+1}^{t,s+1} = r_i^{t,s}
+$$
+
+where t and t-1 are current and previous time step, s+1 and s are current and previous iteration.
+
+#### Derivation
+
+- Lower boundary ($i=1$, stream)
+
+$$
+\begin{aligned}
+f\left(h_{1}^{t,s+1}-h_{1}^{t-1}\right)
+&= \frac{\Delta t\ \sin(\alpha)}{w_{1}\ \Delta x_{1}}
+\left(w_{\frac{3}{2}}\ k_{l_{\frac{3}{2}}}^{t,s}\ h_{\frac{3}{2}}^{t,s}\right) \\
+&\quad + \frac{\Delta t\ \cos(\alpha)}{w_{1}\ \Delta x_{1}}
+\left(
+\frac{w_{\frac{3}{2}}\ k_{l_{\frac{3}{2}}}^{t,s}\ h_{\frac{3}{2}}^{t,s}}{\Delta x_{U_{1}}}
+\left(h_{2}^{t,s+1}-h_{1}^{t,s+1}\right)
+\right) \\
+&\quad + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},1}^{t}\ .
+\end{aligned}
+$$
+
+$$
+a_1 h_0^{t,s+1} + b_1 h_1^{t,s+1} + c_1 h_2^{t,s+1} = r_1
+$$
+
+$$
+\begin{aligned}
+&\text{where:} \\
+&a_1 = 0 \\
+&b_1 = f + \frac{\Delta t\ \cos(\alpha)}{w_1\ \Delta x_1} \cdot \frac{w_{\frac{3}{2}}\ k_{l_{\frac{3}{2}}}^{t,s}\ h_{\frac{3}{2}}^{t,s}}{\Delta x_{U_1}} \\
+&c_1 = -\frac{\Delta t\ \cos(\alpha)}{w_1\ \Delta x_1} \cdot \frac{w_{\frac{3}{2}}\ k_{l_{\frac{3}{2}}}^{t,s}\ h_{\frac{3}{2}}^{t,s}}{\Delta x_{U_1}} \\
+&r_1 = f h_1^{t-1} + \frac{\Delta t\ \sin(\alpha)}{w_1\ \Delta x_1}\left(w_{\frac{3}{2}}\ k_{l_{\frac{3}{2}}}^{t,s}\ h_{\frac{3}{2}}^{t,s}\right) + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},1}^{t}
+\end{aligned}
+$$
+
+
+- Interior nodes ($i=2,\dots,N-1$)
+
+$$
+\begin{aligned}
+f\left(h_{i}^{t,s+1}-h_{i}^{t-1}\right)
+&= \frac{\Delta t\ \sin(\alpha)}{w_{i}\ \Delta x_{i}}
+\left(w_{i+\frac{1}{2}}\ k_{l_{i+\frac{1}{2}}}^{t,s}\ h_{i+\frac{1}{2}}^{t,s}
+ - w_{i-\frac{1}{2}}\ k_{l_{i-\frac{1}{2}}}^{t,s}\ h_{i-\frac{1}{2}}^{t,s}\right) \\
+&\quad + \frac{\Delta t\ \cos(\alpha)}{w_{i}\ \Delta x_{i}}
+\left(
+\frac{w_{i+\frac{1}{2}}\ k_{l_{i+\frac{1}{2}}}^{t,s}\ h_{i+\frac{1}{2}}^{t,s}}{\Delta x_{U_{i}}}
+\left(h_{i+1}^{t,s+1}-h_{i}^{t,s+1}\right) \right. \\
+&\qquad \left. - \frac{w_{i-\frac{1}{2}}\ k_{l_{i-\frac{1}{2}}}^{t,s}\ h_{i-\frac{1}{2}}^{t,s}}{\Delta x_{L_{i}}}
+\left(h_{i}^{t,s+1}-h_{i-1}^{t,s+1}\right)
+\right) \\
+&\quad + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},i}^{t}\ .
+\end{aligned}
+$$
+
+
+$$
+a_i h_{i-1}^{t,s+1} + b_i h_i^{t,s+1} + c_i h_{i+1}^{t,s+1} = r_i
+$$
+
+$$
+\begin{aligned}
+&\text{where:} \\
+&a_i = \frac{\Delta t\ \cos(\alpha)}{w_i\ \Delta x_i} \cdot \frac{w_{i-\frac{1}{2}}\ k_{l_{i-\frac{1}{2}}}^{t,s}\ h_{i-\frac{1}{2}}^{t,s}}{\Delta x_{L_i}} \\
+&b_i = f + \frac{\Delta t\ \cos(\alpha)}{w_i\ \Delta x_i} \left( \frac{w_{i+\frac{1}{2}}\ k_{l_{i+\frac{1}{2}}}^{t,s}\ h_{i+\frac{1}{2}}^{t,s}}{\Delta x_{U_i}} + \frac{w_{i-\frac{1}{2}}\ k_{l_{i-\frac{1}{2}}}^{t,s}\ h_{i-\frac{1}{2}}^{t,s}}{\Delta x_{L_i}} \right) \\
+&c_i = -\frac{\Delta t\ \cos(\alpha)}{w_i\ \Delta x_i} \cdot \frac{w_{i+\frac{1}{2}}\ k_{l_{i+\frac{1}{2}}}^{t,s}\ h_{i+\frac{1}{2}}^{t,s}}{\Delta x_{U_i}} \\
+&r_i = f h_i^{t-1} + \frac{\Delta t\ \sin(\alpha)}{w_i\ \Delta x_i} \left(w_{i+\frac{1}{2}}\ k_{l_{i+\frac{1}{2}}}^{t,s}\ h_{i+\frac{1}{2}}^{t,s} - w_{i-\frac{1}{2}}\ k_{l_{i-\frac{1}{2}}}^{t,s}\ h_{i-\frac{1}{2}}^{t,s}\right) \\
+&\quad + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},i}^{t}
+\end{aligned}
+$$
+
+- Upper boundary ($i=N$, divide)
+
+$$
+\begin{aligned}
+f\left(h_{N}^{t,s+1}-h_{N}^{t-1}\right)
+&= -\frac{\Delta t\ \sin(\alpha)}{w_{N}\ \Delta x_{N}}
+\left(w_{N-\frac{1}{2}}\ k_{l_{N-\frac{1}{2}}}^{t,s}\ h_{N-\frac{1}{2}}^{t,s}\right) \\
+&\quad -\frac{\Delta t\ \cos(\alpha)}{w_{N}\ \Delta x_{N}}
+\left(
+\frac{w_{N-\frac{1}{2}}\ k_{l_{N-\frac{1}{2}}}^{t,s}\ h_{N-\frac{1}{2}}^{t,s}}{\Delta x_{L_{N}}}
+\left(h_{N}^{t,s+1}-h_{N-1}^{t,s+1}\right)
+\right) \\
+&\quad + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},N}^{t}\ .
+\end{aligned}
+$$
+
+
+$$
+a_N h_{N-1}^{t,s+1} + b_N h_N^{t,s+1} + c_N h_{N+1}^{t,s+1} = r_N
+$$
+
+$$
+\begin{aligned}
+&\text{where:} \\
+&a_N = -\frac{\Delta t\ \cos(\alpha)}{w_N\ \Delta x_N} \cdot \frac{w_{N-\frac{1}{2}}\ k_{l_{N-\frac{1}{2}}}^{t,s}\ h_{N-\frac{1}{2}}^{t,s}}{\Delta x_{L_N}} \\
+&b_N = f + \frac{\Delta t\ \cos(\alpha)}{w_N\ \Delta x_N} \cdot \frac{w_{N-\frac{1}{2}}\ k_{l_{N-\frac{1}{2}}}^{t,s}\ h_{N-\frac{1}{2}}^{t,s}}{\Delta x_{L_N}} \\
+&c_N = 0 \\
+&r_N = f h_N^{t-1} - \frac{\Delta t\ \sin(\alpha)}{w_N\ \Delta x_N}\left(w_{N-\frac{1}{2}}\ k_{l_{N-\frac{1}{2}}}^{t,s}\ h_{N-\frac{1}{2}}^{t,s}\right) + \Delta t\ \cos(\alpha)\ R_{\mathrm{sat},N}^{t}
+\end{aligned}
+$$
+
+
+Where $\Delta t$ (seconds) is the h3d time step, i is the lateral node number. $\Delta x_{U_i}$ and ${\Delta x_{L_i}}$ are the distance (m) relative to the center of upper i + 1
+and lower i − 1 node. ${w_i}$ is the width on the center of node i. $i − {\frac{1}{2}}$ and $i + {\frac{1}{2}}$ represent the lower and upper bounds of node i.
+
+
+## IMPLEMENTATION FROM THE CODE
+
+Lower boundary ($i=1$, stream) (NOT SURE HOW THE LAST TERM IN $r_1$ WAS DERIVED)
+
+$$
+\begin{aligned}
+a_1 &= 0, \\
+c_1 &= -\frac{T_{3/2}^{s-1} \cos\alpha \ \Delta t}
+ {\Delta x_{U_i}\ \Delta x_1\ w_1}, \\
+b_1 &= f_{\text{drain},1} - c_1, \\
+r_1 &= f_{\text{drain},1} h_1^{s-1}
+ + \frac{\Delta t}{w_1\Delta x_1}
+ \left[
+ \sin\alpha\ T_{3/2}^{s-1}
+ - \frac{\cos\alpha}{\Delta x_1}
+ w_1 K_{\text{aniso}}
+ \frac{K_{\text{sat},1}}{1000}(h_1^{s-1})^2
+ \right]
+\end{aligned}
+$$
+
+Interior nodes ($i=2,\dots,N-1$)
+
+$$
+\begin{aligned}
+a_i &= -\frac{T_{i-\frac12}^{s-1} \cos\alpha \ \Delta t}
+ {\Delta x_{L_i}\ \Delta x_i\ w_i}, \\
+c_i &= -\frac{T_{i+\frac12}^{s-1} \cos\alpha \ \Delta t}
+ {\Delta x_{U_i}\ \Delta x_i\ w_i}, \\
+b_i &= f_{\text{drain},i} - (a_i + c_i), \\
+r_i &= f_{\text{drain},i} h_i^{s-1}
+ + \frac{\Delta t\sin\alpha}{w_i\Delta x_i}
+ (T_{i+\frac12}^{s-1} - T_{i-\frac12}^{s-1})
+\end{aligned}
+$$
+
+Upper boundary ($i=N$, divide)
+
+$$
+\begin{aligned}
+a_N &= -\frac{T_{N-\frac12}^{s-1} \cos\alpha \ \Delta t}
+ {\Delta x_{L_N}\ \Delta x_N\ w_N}, \\
+c_N &= 0, \\
+b_N &= f_{\text{drain},N} - a_N, \\
+r_N &= f_{\text{drain},N} h_N^{s-1}
+ - \frac{\Delta t\sin\alpha}{w_N\Delta x_N} T_{N-\frac12}^{s-1}
+\end{aligned}
+$$
+
+### Temporal Discretization
+
+A backward-Euler time step is used for stability. Nonlinear terms in transmissivity and porosity are treated by Picard iteration, updating T and $f_{drain}$ until:
+
+$$
+\max_i |h_i^{k+1} - h_i^{k}| < 10^{-4}\ \mathrm{m}
+$$
+
+If convergence fails, the time step is halved adaptively.
+
+$$
+\Delta t_{h3d}^{new} = 0.5\ \Delta t_{h3d}^{old}
+$$
+
+Sub-steps are accumulated until the total integration time equals the parent ELM time step:
+
+$$
+\sum \Delta t_{h3d} = \Delta t_{ELM}
+$$
+
+### Subsurface Runoff and Storage Change
+
+$$
+\Delta S_{\text{sat},i}
+= f_{\text{drain},i}(h_i^{t}-h_i^{t-1}),
+\qquad
+R_{\text{sub},i} = -\Delta S_{\text{sat},i},
+\qquad
+Q_{\text{sub},i} = \frac{R_{\text{sub},i}}{\Delta t}\times1000
+$$
+
+Water-Table Depth
+
+$$
+z_{wt,i} = z_{\text{bed},i} - h_i
+$$
+
+## Subroutines and workflow
+
+### DrainageH3D
+
+Top-level routine for h3D hydrology.
+Responsible for preparing column-level variables (layer thickness, water table index, slope, conductivity),
+computing area-weighted parameters, and calling the h3D solver (H3D_DRI).
+
+### H3D_DRI
+
+Performs the iterative time stepping of the hillslope system:
+
+- Initializes the saturated thickness $h_{sat}$ for each h3D column.
+
+- Computes area-weighted inputs (mean slope, width, transmissivity, decay factor).
+
+- Advances $h_{sat}$ over sub-steps by calling LateralResponse.
+
+- Converts changes in saturated storage to drainage flux:
+
+$$
+\Delta S_{sat} = f_{drain}\ (h^{t} - h^{t-1}),
+\qquad
+Q_{sub} = -\frac{\Delta S_{sat}}{\Delta t}
+$$
+
+Outputs updated water-table depth and drainage rates (qflx_drain_h3d).
+
+### LateralResponse
+
+Solves the implicit Dupuit–Boussinesq system for all h3D nodes in a landunit.
+Constructs a tridiagonal matrix from the finite-difference discretization:
+
+$$
+a_i h_{i-1}^{t,s+1} + b_i h_i^{t,s+1} + c_i h_{i+1}^{t,s+1} = r_i^{t,s}
+$$
+
+Steps:
+- Computes node-specific yield $f_{drain}(c)$ and transmissivity $wK!H(k)$.
+
+- Applies slope-dependent flux terms and boundary conditions.
+
+- Solves using the equation (Tridiagonal_h3D).
+
+- Iterates until the solution converges.
+
+## Outputs
+
+After the h3D solve, the model provides:
+
+- qflx_drain_h3d — subsurface (baseflow) drainage [mm s⁻¹]
+
+- qflx_rsub_sat_h3d — saturation-excess runoff [mm s⁻¹]
+
+- zwt_h3d — updated water-table depth [m]
+
+- f_drain — variable specific yield [–]
+
+- ΔS_sat — change in saturated storage [m]
+
+These outputs replace or augment SIMTOP drainage for h3D-active columns and are fed into the land surface river-routing components of ELM.
diff --git a/components/elm/docs/tech-guide/index.md b/components/elm/docs/tech-guide/index.md
index 6ce636548204..f57325cafb02 100644
--- a/components/elm/docs/tech-guide/index.md
+++ b/components/elm/docs/tech-guide/index.md
@@ -7,3 +7,4 @@ Shortwave radiation model
- [TOP Parameterization](top_solar_parameterization.md):
Parameterization of sub-grid topographical effects on solar radiation.
- [Longwave Radiation](longwave_radiation.md): Longwave radiation model
+- [Hybrid-3D hillslope hydrology](h3d.md): Hybrid-3D hillslope hydrological model
diff --git a/components/elm/src/biogeophys/BalanceCheckMod.F90 b/components/elm/src/biogeophys/BalanceCheckMod.F90
index c39487e5e91d..f908d16e47cd 100755
--- a/components/elm/src/biogeophys/BalanceCheckMod.F90
+++ b/components/elm/src/biogeophys/BalanceCheckMod.F90
@@ -439,6 +439,8 @@ subroutine ColWaterBalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, &
write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
write(iulog,*)'qflx_h2orof_drain = ',qflx_h2orof_drain(indexc)
write(iulog,*)'qflx_ice_runoff_xs = ',qflx_ice_runoff_xs(indexc)
+ l = col_pp%landunit(indexc)
+ write(iulog,*)'h3D soilcolumn = ',indexc,l,lun_pp%coli(l),lun_pp%colf(i)
write(iulog,*)'elm model is stopping'
call endrun(decomp_index=indexc, elmlevel=namec, msg=errmsg(__FILE__, __LINE__))
end if
diff --git a/components/elm/src/biogeophys/HydrologyDrainageMod.F90 b/components/elm/src/biogeophys/HydrologyDrainageMod.F90
index e14c229160ce..0f82b23c3716 100755
--- a/components/elm/src/biogeophys/HydrologyDrainageMod.F90
+++ b/components/elm/src/biogeophys/HydrologyDrainageMod.F90
@@ -38,6 +38,7 @@ module HydrologyDrainageMod
subroutine HydrologyDrainage(bounds, &
num_nolakec, filter_nolakec, &
num_hydrologyc, filter_hydrologyc, &
+ num_h3dc, filter_h3dc, &
num_urbanc, filter_urbanc, &
num_do_smb_c, filter_do_smb_c, &
atm2lnd_vars, glc2lnd_vars, &
@@ -57,7 +58,7 @@ subroutine HydrologyDrainage(bounds, &
use TopounitDataType , only : top_ws
use atm2lndType , only : atm2lnd_type
use elm_varpar , only : nlevgrnd, nlevurb, nlevsoi
- use SoilHydrologyMod , only : ELMVICMap, Drainage
+ use SoilHydrologyMod , only : ELMVICMap, Drainage, DrainageH3D
use elm_varctl , only : use_vsfm, use_IM2_hillslope_hydrology
!
! !ARGUMENTS:
@@ -66,6 +67,8 @@ subroutine HydrologyDrainage(bounds, &
integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points
integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points
+ integer , intent(in) :: num_h3dc ! number of column h3d points in column filter
+ integer , intent(in) :: filter_h3dc(:) ! column filter for h3d points
integer , intent(in) :: num_urbanc ! number of column urban points in column filter
integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
integer , intent(in) :: num_do_smb_c ! number of bareland columns in which SMB is calculated, in column filter
@@ -145,9 +148,15 @@ subroutine HydrologyDrainage(bounds, &
#endif
if (.not. use_vsfm) then
- call Drainage(bounds, num_hydrologyc, filter_hydrologyc, &
- num_urbanc, filter_urbanc,&
- soilhydrology_vars, soilstate_vars, dtime)
+ if (use_h3d) then
+ call DrainageH3D(bounds, num_h3dc, filter_h3dc, num_hydrologyc, &
+ filter_hydrologyc, num_urbanc, filter_urbanc, &
+ soilhydrology_vars, soilstate_vars, dtime)
+ else
+ call Drainage(bounds, num_hydrologyc, filter_hydrologyc, &
+ num_urbanc, filter_urbanc,&
+ soilhydrology_vars, soilstate_vars, dtime)
+ endif
endif
#ifndef _OPENACC
diff --git a/components/elm/src/biogeophys/SoilHydrologyMod.F90 b/components/elm/src/biogeophys/SoilHydrologyMod.F90
index e226adf96ebe..05cafd862a9e 100644
--- a/components/elm/src/biogeophys/SoilHydrologyMod.F90
+++ b/components/elm/src/biogeophys/SoilHydrologyMod.F90
@@ -36,6 +36,10 @@ module SoilHydrologyMod
public :: Drainage ! Calculate subsurface drainage
public :: DrainageVSFM ! Calculate subsurface drainage for VSFM
public :: ELMVICMap
+ public :: DrainageH3D ! Calculate subsurface drainage using h3D
+ private :: H3D_DRI ! Calculate subsurface drainage using h3D
+ private :: LateralResponse ! Calculate subsurface drainage using h3D
+ private :: Tridiagonal_h3D
!-----------------------------------------------------------------------
contains
@@ -82,6 +86,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
real(r8) :: top_moist(bounds%begc:bounds%endc) !temporary, soil moisture in top VIC layers
real(r8) :: top_max_moist(bounds%begc:bounds%endc) !temporary, maximum soil moisture in top VIC layers
real(r8) :: top_ice(bounds%begc:bounds%endc) !temporary, ice len in top VIC layers
+ integer :: jwt(bounds%begc:bounds%endc) !index of the soil layer right above the water table
character(len=32) :: subname = 'SurfaceRunoff' !subroutine name
!-----------------------------------------------------------------------
@@ -89,6 +94,7 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
snl => col_pp%snl , & ! Input: [integer (:) ] minus number of snow layers
dz => col_pp%dz , & ! Input: [real(r8) (:,:) ] layer depth (m)
nlev2bed => col_pp%nlevbed , & ! Input: [integer (:) ] number of layers to bedrock
+ zi => col_pp%zi , & ! Input: [real(r8) (:,:) ] interface level at the bottom of each layer (m)
sucsat => soilstate_vars%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
watsat => soilstate_vars%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
@@ -147,6 +153,26 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
end do
end do
+ ! The layer index of the first unsaturated layer, i.e., the layer right
+ ! above the water table
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ jwt(c) = nlevbed
+ ! allow jwt to equal zero when zwt is in top layer
+ do j = 1,nlevbed
+ if (zwt(c) <= zi(c,j)) then
+ if (zengdecker_2009_with_var_soil_thick .and. zwt(c) == zi(c,nlevbed)) then
+ exit
+ else
+ jwt(c) = j-1
+ exit
+ end if
+ end if
+ end do
+ end do
+
! Saturated fraction
do fc = 1, num_hydrologyc
@@ -156,6 +182,9 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
if (zengdecker_2009_with_var_soil_thick) then
nlevbed = nlev2bed(c)
fff(c) = 0.5_r8 * col_pp%zi(c,nlevsoi) / min(col_pp%zi(c,nlevbed), col_pp%zi(c,nlevsoi))
+ if (use_h3d) then
+ fff(c) = 2.0_r8 * col_pp%zi(c,nlevsoi) / min(col_pp%zi(c,nlevbed),col_pp%zi(c,nlevsoi))
+ end if
end if
if (use_vichydro) then
top_moist(c) = 0._r8
@@ -512,6 +541,9 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
qinmax=minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
else
qinmax=(1._r8 - fsat(c)) * minval(10._r8**(-e_ice*(icefrac(c,1:3)))*hksat(c,1:3))
+ if (use_h3d) then
+ qinmax=(1._r8 - fsat(c)) * minval(hksat(c,1:3))
+ end if
end if
end if
@@ -2181,4 +2213,1389 @@ subroutine ELMVICMap(bounds, numf, filter, &
end subroutine ELMVICMap
+ !-----------------------------------------------------------------------
+
+
+ subroutine DrainageH3D(bounds, num_h3dc,filter_h3dc,num_hydrologyc,
+filter_hydrologyc, num_urbanc, filter_urbanc, &
+ soilhydrology_vars, soilstate_vars, dtime)
+ !
+ ! !DESCRIPTION:
+ ! Calculate subsurface drainage, h3d approach
+ ! !yhfang
+ !
+ ! !USES:
+ use clm_time_manager , only : get_step_size
+ use clm_varpar , only : nlevsoi, nlevgrnd, nlayer, nlayert
+ use clm_varpar , only : nh3dc_per_lunit
+ use clm_varcon , only : pondmx, tfrz, watmin,rpi, secspday, nlvic
+ use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv
+ use abortutils , only : endrun
+ use clm_varctl , only : use_vsfm, use_var_soil_thick, use_h3d
+ use SoilWaterMovementMod, only : zengdecker_2009_with_var_soil_thick
+ use pftvarcon , only : rsub_top_globalmax
+ !
+ ! !ARGUMENTS:
+ type(bounds_type) , intent(in) :: bounds
+ integer , intent(in) :: num_h3dc ! number of column h3d points in column filter
+ integer , intent(in) :: filter_h3dc(:) ! columnfilter for h3d points
+ integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter
+ integer , intent(in) :: num_urbanc ! number of column urban points in column filter
+ integer , intent(in) :: filter_urbanc(:) ! column filter for urban points
+ integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points
+ type(temperature_type) , intent(in) :: temperature_vars
+ type(soilstate_type) , intent(in) :: soilstate_vars
+ type(soilhydrology_type) , intent(inout) :: soilhydrology_vars
+ real(r8) , intent(in) :: dtime
+ !
+ ! !LOCAL VARIABLES:
+ character(len=32) :: subname = 'DrainageH3d' ! subroutine name
+ integer :: c,j,fc,i ! indices
+ integer :: c0,l ! indices yhfang
+ integer :: nlevbed ! # layers to bedrock
+ real(r8) :: xs(bounds%begc:bounds%endc) ! water needed to bring soil moisture to watmin (mm)
+ real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevgrnd) ! layer thickness (mm)
+ integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-)
+ real(r8) :: rsub_bot(bounds%begc:bounds%endc) ! subsurface runoff - bottom drainage (mm/s)
+ real(r8) :: rsub_top(bounds%begc:bounds%endc) ! subsurface runoff - topographic control (mm/s)
+ real(r8) :: fff(bounds%begc:bounds%endc) ! decay factor (m-1)
+ real(r8) :: xsi(bounds%begc:bounds%endc) ! excess soil water above saturation at layer i (mm)
+ real(r8) :: xsia(bounds%begc:bounds%endc) ! available pore space at layer i (mm)
+ real(r8) :: xs1(bounds%begc:bounds%endc) ! excess soil water above saturation at layer 1 (mm)
+ real(r8) :: smpfz(1:nlevsoi) ! matric potential of layer right above water table (mm)
+ real(r8) :: wtsub ! summation of hk*dzmm for layers below water table (mm**2/s)
+ real(r8) :: rous ! aquifer yield (-)
+ real(r8) :: wh ! smpfz(jwt)-z(jwt) (mm)
+ real(r8) :: wh_zwt ! water head at the water table depth (mm)
+ real(r8) :: ws ! summation of pore space of layers below water table (mm)
+ real(r8) :: s_node ! soil wetness (-)
+ real(r8) :: dzsum ! summation of dzmm of layers below water table (mm)
+ real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-)
+ real(r8) :: fracice_rsub(bounds%begc:bounds%endc) ! fractional impermeability of soil layers (-)
+ real(r8) :: ka ! hydraulic conductivity of the aquifer (mm/s)
+ real(r8) :: dza ! fff*(zwt-z(jwt)) (-)
+ real(r8) :: available_h2osoi_liq ! available soil liquid water in a layer
+ real(r8) :: rsub_top_max
+ real(r8) :: h2osoi_vol
+ real(r8) :: imped
+ real(r8) :: rsub_top_tot
+ real(r8) :: rsub_top_layer
+ real(r8) :: qcharge_tot
+ real(r8) :: qcharge_layer
+ real(r8) :: theta_unsat
+ real(r8) :: f_unsat
+ real(r8) :: s_y
+ integer :: k,k_frz,k_perch
+ real(r8) :: sat_lev
+ real(r8) :: s1
+ real(r8) :: s2
+ real(r8) :: m
+ real(r8) :: b
+ real(r8) :: q_perch
+ real(r8) :: q_perch_max
+ real(r8) :: vol_ice
+ real(r8) :: dsmax_tmp(bounds%begc:bounds%endc) ! temporary variable for ARNO subsurface runoff calculation
+ real(r8) :: rsub_tmp ! temporary variable for ARNO subsurface runoff calculation
+ real(r8) :: frac ! temporary variable for ARNO subsurface runoff calculation
+ real(r8) :: rel_moist ! relative moisture, temporary variable
+ real(r8) :: wtsub_vic ! summation of hk*dzmm for layers in the third VIC layer
+
+ real(r8) :: zwt_h3d(bounds%begc:bounds%endc) ! !m
+ real(r8) :: tmp_rsub_top(bounds%begc:bounds%endc) ! !m/s
+ real(r8) :: h3d_rsub_top_max(bounds%begc:bounds%endc) ! !m/s
+
+ real(r8) :: tmp_sum_drain1 !for test only
+ real(r8) :: tmp_sum_drain2 !for test only
+ !-----------------------------------------------------------------------
+
+ associate( &
+ z => col_pp%z , &
+! Input: [real(r8) (:,:) ] layer depth (m)
+ zi => col_pp%zi , &
+! Input: [real(r8) (:,:) ] interface level below a "z" level (m)
+ dz => col_pp%dz , &
+! Input: [real(r8) (:,:) ] layer depth (m)
+ snl => col_pp%snl , &
+! Input: [integer (:) ] number of snow layers
+ nlev2bed => col_pp%nlevbed , &
+! Input: [integer (:) ] number of layers to bedrock
+
+ t_soisno => col_es%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin)
+
+ h2osfc => col_ws%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm)
+
+ bsw => soilstate_vars%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
+ hksat => soilstate_vars%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s)
+ sucsat => soilstate_vars%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
+ watsat => soilstate_vars%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
+ eff_porosity => soilstate_vars%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice
+ hk_l => soilstate_vars%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s)
+ smp_l => soilstate_vars%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matric potential (mm)
+
+ depth => soilhydrology_vars%depth_col , & ! Input: [real(r8) (:,:) ] VIC soil depth
+ c_param => soilhydrology_vars%c_param_col , & ! Input: [real(r8) (:) ] baseflow exponent (Qb)
+ Dsmax => soilhydrology_vars%dsmax_col , & ! Input: [real(r8) (:) ] max. velocity of baseflow (mm/day)
+ max_moist => soilhydrology_vars%max_moist_col , & ! Input: [real(r8) (:,:) ] maximum soil moisture (ice + liq)
+ moist => soilhydrology_vars%moist_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm)
+ Ds => soilhydrology_vars%ds_col , & ! Input: [real(r8) (:) ] fracton of Dsmax where non-linear baseflow begins
+ Wsvic => soilhydrology_vars%Wsvic_col , & ! Input: [real(r8) (:) ] fraction of maximum soil moisutre where non-liear base flow occurs
+ icefrac => soilhydrology_vars%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer
+ hkdepth => soilhydrology_vars%hkdepth_col , & ! Input: [real(r8) (:) ] decay factor (m)
+ frost_table => soilhydrology_vars%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m)
+ zwt => soilhydrology_vars%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m)
+ zwt_perched => soilhydrology_vars%zwt_perched_col , & ! Input: [real(r8) (:) ] perched water table depth (m)
+ wa => soilhydrology_vars%wa_col , & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm)
+ ice => soilhydrology_vars%ice_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm)
+ qcharge => soilhydrology_vars%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s)
+ origflag => soilhydrology_vars%origflag , & ! Input: logical
+ h2osfcflag => soilhydrology_vars%h2osfcflag , & ! Input: logical
+
+ qflx_snwcp_liq => col_wf%qflx_snwcp_liq_col , & ! Output: [real(r8) (:) ] excess rainfall due to snow capping (mm H2O /s) [+]
+ qflx_snwcp_ice => col_wf%qflx_snwcp_ice_col , & ! Output: [real(r8) (:) ] excess snowfall due to snow capping (mm H2O /s) [+]
+ !qflx_dew_grnd => waterflux_vars%qflx_dew_grnd_col , & !
+ !Output: [real(r8) (:) ] ground surface dew formation (mm H2O /s)
+ ![+]
+ !qflx_dew_snow => waterflux_vars%qflx_dew_snow_col , & !
+ !Output: [real(r8) (:) ] surface dew added to snow pack (mm H2O /s)
+ ![+]
+ !qflx_sub_snow => waterflux_vars%qflx_sub_snow_col , & !
+ !Output: [real(r8) (:) ] sublimation rate from snow pack (mm H2O /s)
+ ![+]
+ qflx_drain => col_wf%qflx_drain_col , & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s)
+ qflx_qrgwl => col_wf%qflx_qrgwl_col , & ! Output: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes (mm H2O /s)
+ qflx_rsub_sat => col_wf%qflx_rsub_sat_col , & ! Output: [real(r8) (:) ] soil saturation excess [mm h2o/s]
+ qflx_drain_perched => col_wf%qflx_drain_perched_col , & ! Output: [real(r8) (:) ] perched wt sub-surface runoff (mm H2O /s)
+
+ qflx_rsub_sat_h3d => col_wf%qflx_rsub_sat_h3dcol , & ! Output: [real(r8) (:) ] soil saturation excess from h3d [mm h2o/s]
+ qflx_drain_h3d => col_wf%qflx_drain_h3dcol , & ! Output: [real(r8) (:) ] sub-surface runoff from h3d(mm H2O /s)
+ h3d_zwt_lun => soilhydrology_vars%h3d_zwt_lun , & ! Output: [real(r8) (:,:) ] water table depth at h3d column !yhfang
+ h3d_imped => soilhydrology_vars%imped_h3d_col , & ! Inout : [real(r8) (:) ] scaling facfor due to frozen soil
+
+ h2osoi_liq => col_ws%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2)
+ h2osoi_ice => col_ws%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2)
+ )
+
+ ! Convert layer thicknesses from m to mm
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ do j = 1,nlevbed
+ dzmm(c,j) = dz(c,j)*1.e3_r8
+
+ vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice))
+ icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j))
+ end do
+ end do
+
+ ! Initial set
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ qflx_drain(c) = 0._r8
+ rsub_bot(c) = 0._r8
+ qflx_rsub_sat(c) = 0._r8
+ rsub_top(c) = 0._r8
+ fracice_rsub(c) = 0._r8
+ qflx_qrgwl(c) = 0._r8
+
+ !yhfang
+ qflx_rsub_sat_h3d(c) = 0._r8
+ qflx_drain_h3d(c) = 0._r8
+ end do
+
+ ! The layer index of the first unsaturated layer, i.e., the layer right
+ ! above
+ ! the water table
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ jwt(c) = nlevbed
+ ! allow jwt to equal zero when zwt is in top layer
+ do j = 1,nlevbed
+ if(zwt(c) <= zi(c,j)) then
+ if (zengdecker_2009_with_var_soil_thick .and. zwt(c) == zi(c,nlevbed)) then
+ exit
+ else
+ jwt(c) = j-1
+ exit
+ end if
+ end if
+ enddo
+ end do
+
+ rous = 0.2_r8
+
+ !== BASEFLOW ==================================================
+ ! perched water table code
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+
+ ! specify maximum drainage rate
+ q_perch_max = 1.e-5_r8 * sin(col_pp%topo_slope(c) * (rpi/180._r8))
+
+ ! if layer containing water table is frozen, compute the following:
+ ! frost table, perched water table, and drainage from perched
+ ! saturated layer
+
+ ! define frost table as first frozen layer with unfrozen layer above
+ ! it
+ if(t_soisno(c,1) > tfrz) then
+ k_frz=nlevbed
+ else
+ k_frz=1
+ endif
+
+ do k=2, nlevbed
+ if (t_soisno(c,k-1) > tfrz .and. t_soisno(c,k) <= tfrz) then
+ k_frz=k
+ exit
+ endif
+ enddo
+
+ frost_table(c)=z(c,k_frz)
+
+ ! initialize perched water table to frost table, and
+ ! qflx_drain_perched(c) to zero
+ zwt_perched(c)=frost_table(c)
+ qflx_drain_perched(c) = 0._r8
+
+ !=================== water table above frost table
+ !=============================
+ ! if water table is above frost table, do not use topmodel baseflow
+ ! formulation
+
+ if (zwt(c) < frost_table(c) .and. t_soisno(c,k_frz) <= tfrz &
+ .and. origflag == 0) then
+ ! compute drainage from perched saturated region
+ wtsub = 0._r8
+ q_perch = 0._r8
+ do k = jwt(c)+1, k_frz
+ imped=10._r8**(-e_ice*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi,k+1)))))
+ q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k)
+ wtsub = wtsub + dzmm(c,k)
+ end do
+ if (wtsub > 0._r8) q_perch = q_perch/wtsub
+
+ qflx_drain_perched(c) = q_perch_max * q_perch &
+ *(frost_table(c) - zwt(c))
+
+ ! remove drainage from perched saturated layers
+ rsub_top_tot = - qflx_drain_perched(c) * dtime
+ do k = jwt(c)+1, k_frz
+ rsub_top_layer=max(rsub_top_tot,-(h2osoi_liq(c,k)-watmin))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ if (use_vsfm) then
+ rsub_top_layer = 0._r8
+ endif
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+
+ h2osoi_liq(c,k) = h2osoi_liq(c,k) + rsub_top_layer
+
+ if (rsub_top_tot >= 0.) then
+ zwt(c) = zwt(c) - rsub_top_layer/eff_porosity(c,k)/1000._r8
+ exit
+ else
+ zwt(c) = zi(c,k)
+ endif
+ enddo
+
+ ! if rsub_top_tot is greater than available water (above frost
+ ! table),
+ ! then decrease qflx_drain_perched by residual amount for water
+ ! balance
+ qflx_drain_perched(c) = qflx_drain_perched(c) + rsub_top_tot/dtime
+
+ !-- recompute jwt
+ !---------------------------------------------------------
+ ! allow jwt to equal zero when zwt is in top layer
+ jwt(c) = nlevbed
+ do j = 1,nlevbed
+ if(zwt(c) <= zi(c,j)) then
+ if (zengdecker_2009_with_var_soil_thick .and. zwt(c) == zi(c,nlevbed)) then
+ exit
+ else
+ jwt(c) = j-1
+ exit
+ end if
+ end if
+ enddo
+ else
+ !=================== water table below frost table
+ !=============================
+ !-- compute possible perched water table *and* groundwater table
+ !afterwards
+ ! locate perched water table from bottom up starting at frost table
+ ! sat_lev is an arbitrary saturation level used to determine
+ ! perched water table
+ sat_lev=0.9
+
+ k_perch=1
+ do k=k_frz,1,-1
+ h2osoi_vol = h2osoi_liq(c,k)/(dz(c,k)*denh2o) &
+ + h2osoi_ice(c,k)/(dz(c,k)*denice)
+
+ if (h2osoi_vol/watsat(c,k) <= sat_lev) then
+ k_perch=k
+ exit
+ endif
+ enddo
+
+ ! if frost_table = nlevsoi, only compute perched water table if
+ ! frozen
+ if (t_soisno(c,k_frz) > tfrz) k_perch=k_frz
+
+ ! if perched water table exists
+ if (k_frz > k_perch) then
+ ! interpolate between k_perch and k_perch+1 to find perched
+ ! water table height
+ s1 = (h2osoi_liq(c,k_perch)/(dz(c,k_perch)*denh2o) &
+ + h2osoi_ice(c,k_perch)/(dz(c,k_perch)*denice))/watsat(c,k_perch)
+ s2 = (h2osoi_liq(c,k_perch+1)/(dz(c,k_perch+1)*denh2o) &
+ + h2osoi_ice(c,k_perch+1)/(dz(c,k_perch+1)*denice))/watsat(c,k_perch+1)
+
+ m=(z(c,k_perch+1)-z(c,k_perch))/(s2-s1)
+ b=z(c,k_perch+1)-m*s2
+ zwt_perched(c)=max(0._r8,m*sat_lev+b)
+
+ ! compute drainage from perched saturated region
+ wtsub = 0._r8
+ q_perch = 0._r8
+ do k = k_perch, k_frz
+ imped=10._r8**(-e_ice*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi,k+1)))))
+ q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k)
+ wtsub = wtsub + dzmm(c,k)
+ end do
+ if (wtsub > 0._r8) q_perch = q_perch/wtsub
+
+ qflx_drain_perched(c) = q_perch_max * q_perch &
+ *(frost_table(c) - zwt_perched(c))
+
+ ! no perched water table drainage if using original formulation
+ if(origflag == 1) qflx_drain_perched(c) = 0._r8
+
+ ! remove drainage from perched saturated layers
+ rsub_top_tot = - qflx_drain_perched(c) * dtime
+ do k = k_perch+1, k_frz
+ rsub_top_layer=max(rsub_top_tot,-(h2osoi_liq(c,k)-watmin))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ if (use_vsfm) rsub_top_layer = 0._r8
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+
+ h2osoi_liq(c,k) = h2osoi_liq(c,k) + rsub_top_layer
+
+ if (rsub_top_tot >= 0.) then
+ zwt_perched(c) = zwt_perched(c) - rsub_top_layer/eff_porosity(c,k)/1000._r8
+ exit
+ else
+ zwt_perched(c) = zi(c,k)
+ endif
+
+ enddo
+
+ ! if rsub_top_tot is greater than available water (above frost
+ ! table),
+ ! then decrease qflx_drain_perched by residual amount for
+ ! water balance
+ qflx_drain_perched(c) = qflx_drain_perched(c) + rsub_top_tot/dtime
+
+ else
+ qflx_drain_perched(c) = 0._r8
+ endif !k_frz > k_perch
+
+ !-- Topographic runoff
+ !----------------------------------------------------------------------
+ fff(c) = 1._r8/ hkdepth(c)
+ dzsum = 0._r8
+ icefracsum = 0._r8
+ do j = max(jwt(c),1), nlevbed
+ dzsum = dzsum + dzmm(c,j)
+ icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j)
+ end do
+ ! add ice impedance factor to baseflow
+ if(origflag == 1) then
+ if (use_vichydro) then
+ call endrun(msg="VICHYDRO is not available for origflag=1"//errmsg(__FILE__, __LINE__))
+ else
+ fracice_rsub(c) = max(0._r8,exp(-3._r8*(1._r8-(icefracsum/dzsum))) &
+ - exp(-3._r8))/(1.0_r8-exp(-3._r8))
+ imped=(1._r8 - fracice_rsub(c))
+ rsub_top_max = 5.5e-3_r8
+ end if
+ else
+ if (use_vichydro) then
+ imped=10._r8**(-e_ice*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer)))
+ dsmax_tmp(c) = Dsmax(c) * dtime/ secspday !mm/day->mm/dtime
+ rsub_top_max = dsmax_tmp(c)
+ else
+ imped=10._r8**(-e_ice*(icefracsum/dzsum))
+ rsub_top_max = min(10._r8 * sin((rpi/180.) * col_pp%topo_slope(c)), rsub_top_globalmax)
+ end if
+ endif
+ if (use_vichydro) then
+ ! ARNO model for the bottom soil layer (based on bottom soil
+ ! layer
+ ! moisture from previous time step
+ ! use watmin instead for resid_moist to be consistent with
+ ! default hydrology
+ rel_moist = (moist(c,nlayer) - watmin)/(max_moist(c,nlayer)-watmin)
+ frac = (Ds(c) * rsub_top_max )/Wsvic(c)
+ rsub_tmp = (frac * rel_moist)/dtime
+ if(rel_moist > Wsvic(c))then
+ frac = (rel_moist - Wsvic(c))/(1.0_r8 - Wsvic(c))
+ rsub_tmp = rsub_tmp + (rsub_top_max * (1.0_r8 - Ds(c)/Wsvic(c)) *frac**c_param(c))/dtime
+ end if
+ rsub_top(c) = imped * rsub_tmp
+ ! make sure baseflow isn't negative
+ rsub_top(c) = max(0._r8, rsub_top(c))
+ else
+ if (jwt(c) == nlevbed .and. zengdecker_2009_with_var_soil_thick) then
+ rsub_top(c) = 0._r8
+ else
+ rsub_top(c) = imped * rsub_top_max* exp(-fff(c)*zwt(c))
+ end if
+ end if
+
+ if (use_vsfm) rsub_top(c) = 0._r8
+
+ !no need to compute drainage for h3d soil columns
+ if (use_h3d .and. col_pp%h3d_active(c)) then
+ tmp_rsub_top(c) = rsub_top(c) / 1000._r8 !mm/s -> m/s
+ rsub_top(c) = 0._r8
+ h3d_imped(c) = imped
+ h3d_rsub_top_max(c) = rsub_top_max
+ cycle
+ end if
+
+ ! use analytical expression for aquifer specific yield
+ rous = watsat(c,nlevbed) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevbed))**(-1./bsw(c,nlevbed)))
+ rous=max(rous,0.02_r8)
+
+
+ !-- water table is below the soil column
+ !--------------------------------------
+ if(jwt(c) == nlevbed) then
+ if (zengdecker_2009_with_var_soil_thick) then
+ if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
+ zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
+ end if
+ rsub_top(c) = imped * rsub_top_max * exp(-fff(c) * zwt(c))
+ rsub_top_tot = - rsub_top(c) * dtime
+ s_y = watsat(c,nlevbed) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevbed))**(-1./bsw(c,nlevbed)))
+ s_y=max(s_y,0.02_r8)
+ rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,nlevbed) - zwt(c))*1.e3))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ h2osoi_liq(c,nlevbed) = h2osoi_liq(c,nlevbed) + rsub_top_layer
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+ if (rsub_top_tot >= 0.) then
+ zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8
+ else
+ zwt(c) = zi(c,nlevbed)
+ end if
+ if (rsub_top_tot < 0.) then
+ rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
+ rsub_top_tot = 0.
+ end if
+ else
+ wa(c) = wa(c) - rsub_top(c) * dtime
+ zwt(c) = zwt(c) + (rsub_top(c) * dtime)/1000._r8/rous
+ h2osoi_liq(c,nlevsoi) = h2osoi_liq(c,nlevsoi) + max(0._r8,(wa(c)-5000._r8))
+ wa(c) = min(wa(c), 5000._r8)
+ end if
+ else
+ !-- water table within soil layers 1-9
+ !-------------------------------------
+ !============================== RSUB_TOP
+ !=========================================
+ !-- Now remove water via rsub_top
+ rsub_top_tot = - rsub_top(c) * dtime
+ !should never be positive... but include for completeness
+ if(rsub_top_tot > 0.) then !rising water table
+
+ call endrun(msg="RSUB_TOP IS POSITIVE in Drainage!"//errmsg(__FILE__, __LINE__))
+
+ else ! deepening water table
+ if (use_vichydro) then
+ wtsub_vic = 0._r8
+ do j = (nlvic(1)+nlvic(2)+1), nlevbed
+ wtsub_vic = wtsub_vic + hk_l(c,j)*dzmm(c,j)
+ end do
+
+ do j = (nlvic(1)+nlvic(2)+1), nlevbed
+ rsub_top_layer=max(rsub_top_tot, rsub_top_tot*hk_l(c,j)*dzmm(c,j)/wtsub_vic)
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ if (use_vsfm) rsub_top_layer = 0._r8
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+ end do
+ else
+ do j = jwt(c)+1, nlevbed
+ ! use analytical expression for specific yield
+ s_y = watsat(c,j) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j)))
+ s_y=max(s_y,0.02_r8)
+
+ rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ if (use_vsfm) rsub_top_layer = 0._r8
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer
+
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+
+ if (rsub_top_tot >= 0.) then
+ zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8
+
+ exit
+ else
+ zwt(c) = zi(c,j)
+ endif
+ enddo
+ end if
+
+ !-- remove residual rsub_top
+ !---------------------------------------------
+ if (zengdecker_2009_with_var_soil_thick) then
+ if (rsub_top_tot < 0.) then
+ rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
+ rsub_top_tot = 0._r8
+ end if
+ else
+ zwt(c) = zwt(c) - rsub_top_tot/1000._r8/rous
+ wa(c) = wa(c) + rsub_top_tot
+ end if
+ endif
+
+ !-- recompute jwt
+ !---------------------------------------------------------
+ ! allow jwt to equal zero when zwt is in top layer
+ jwt(c) = nlevbed
+ do j = 1,nlevbed
+ if(zwt(c) <= zi(c,j)) then
+ if (zengdecker_2009_with_var_soil_thick .and. zwt(c) == zi(c,nlevbed)) then
+ exit
+ else
+ jwt(c) = j-1
+ exit
+ end if
+ end if
+ enddo
+ end if! end of jwt if construct
+
+ zwt(c) = max(0.0_r8,zwt(c))
+ zwt(c) = min(80._r8,zwt(c))
+
+ endif
+
+ end do
+
+!====================================================================================================================
+ !h3d code
+ !for filter_h3dc only
+ if (use_h3d) then
+
+ do fc = 1, num_h3dc
+ c = filter_h3dc(fc)
+ nlevbed = nlev2bed(c)
+ if(jwt(c) == nlevbed) then
+ if (zengdecker_2009_with_var_soil_thick) then
+ if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
+ zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
+ end if
+ end if
+ end if
+ end do
+
+ CALL H3D_DRI(bounds, num_h3dc, filter_h3dc, jwt,zwt_h3d, tmp_rsub_top,h3d_rsub_top_max,fff,&
+ soilhydrology_vars, soilstate_vars, dtime)
+
+
+ !need to double check if we need to re-compute zwt, or use zwt from h3d?
+ do fc = 1, num_h3dc
+ c = filter_h3dc(fc)
+ nlevbed = nlev2bed(c)
+ !-- water table is below the soil column
+ !--------------------------------------
+ if(jwt(c) == nlevbed) then
+ rsub_top(c) = h3d_imped(c) * qflx_drain_h3d(c) !from h3d
+ rsub_top_tot = - rsub_top(c) * dtime
+
+ if (zengdecker_2009_with_var_soil_thick) then
+ if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
+ zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
+ end if
+ !rsub_top(c) = imped * rsub_top_max * exp(-fff(c) * zwt(c))
+
+ if(rsub_top_tot > 0.) then !rising water table
+ !h3d allows water table rise
+ do j = jwt(c), 1,-1
+ !scs: use analytical expression for specific yield
+ s_y = watsat(c,j) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j)))
+ s_y=max(s_y,0.02_r8)
+
+ rsub_top_layer=min(rsub_top_tot,(s_y*(zwt(c) - zi(c,j-1))*1.e3))
+ rsub_top_layer=max(rsub_top_layer,0._r8)
+
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer
+ if(s_y > 0._r8) zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8
+
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+ if (rsub_top_tot <= 0.) then
+ if (rsub_top_tot < 0.) write(iulog,*) 'negtive rsub_top_tot',c,rsub_top_tot
+ exit
+ end if
+ enddo
+
+ if (rsub_top_tot > 0._r8) then
+ qflx_rsub_sat(c) = qflx_rsub_sat(c) + rsub_top_tot / dtime
+ end if
+
+ else
+ s_y = watsat(c,nlevbed) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevbed))**(-1./bsw(c,nlevbed)))
+ s_y=max(s_y,0.02_r8)
+
+ rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,nlevbed) - zwt(c))*1.e3))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ h2osoi_liq(c,nlevbed) = h2osoi_liq(c,nlevbed) + rsub_top_layer
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+ zwt(c) = zi(c,nlevbed)
+
+ if (rsub_top_tot < 0.) then
+ rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
+ rsub_top_tot = 0.
+ end if
+ end if
+ else
+ !rsub_top(c) = qflx_drain_h3d(c) !from h3d
+ rsub_top(c) = h3d_imped(c) * qflx_drain_h3d(c) !from h3d
+ wa(c) = wa(c) - rsub_top(c) * dtime
+ zwt(c) = zwt(c) + (rsub_top(c) * dtime)/1000._r8/rous
+ h2osoi_liq(c,nlevsoi) = h2osoi_liq(c,nlevsoi) + max(0._r8,(wa(c)-5000._r8))
+ wa(c) = min(wa(c), 5000._r8)
+ end if
+ else
+ !-- water table within soil layers 1-9
+ !-------------------------------------
+ !============================== RSUB_TOP
+ !=========================================
+ !-- Now remove water via rsub_top
+ !rsub_top(c) = qflx_drain_h3d(c) !from h3d
+ rsub_top(c) = h3d_imped(c) * qflx_drain_h3d(c) !from h3d
+ rsub_top_tot = - rsub_top(c) * dtime
+
+ if(rsub_top_tot > 0._r8) then !rising water table
+ !h3d allows water table rise
+ do j = jwt(c)+1, 1,-1
+ !scs: use analytical expression for specific yield
+ s_y = watsat(c,j) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j)))
+ s_y=max(s_y,0.02_r8)
+
+ rsub_top_layer=min(rsub_top_tot,(s_y*(zwt(c) - zi(c,j-1))*1.e3))
+ rsub_top_layer=max(rsub_top_layer,0._r8)
+
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer
+ if(s_y > 0._r8) zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8
+
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+ if (rsub_top_tot <= 0.) then
+ if (rsub_top_tot < 0.) write(iulog,*) 'negtive rsub_top_tot',c,rsub_top_tot
+ exit
+ end if
+ enddo
+
+ if (rsub_top_tot > 0._r8) then
+ qflx_rsub_sat_h3d(c) = qflx_rsub_sat_h3d(c) + rsub_top_tot/dtime
+ !qflx_rsub_sat(c) = qflx_rsub_sat(c) + rsub_top_tot / dtime
+ end if
+
+ else ! deepening water table
+ if (use_vichydro) then
+ call endrun(msg="H3D cannot work with VICHYDRO!"//errmsg(__FILE__, __LINE__))
+ else
+ do j = jwt(c)+1, nlevbed
+ ! use analytical expression for specific yield
+ s_y = watsat(c,j) &
+ * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j)))
+ s_y=max(s_y,0.02_r8)
+
+ rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3))
+ rsub_top_layer=min(rsub_top_layer,0._r8)
+ if (use_vsfm) rsub_top_layer = 0._r8
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer
+
+ rsub_top_tot = rsub_top_tot - rsub_top_layer
+
+ if (rsub_top_tot >= 0.) then
+ zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8
+
+ exit
+ else
+ zwt(c) = zi(c,j)
+ endif
+ enddo
+ end if
+
+ !-- remove residual rsub_top
+ !---------------------------------------------
+ if (zengdecker_2009_with_var_soil_thick) then
+ if (rsub_top_tot < 0.) then
+ write(iulog,*) 'adjust0 rsub_top',c,rsub_top(c)
+ rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
+ rsub_top_tot = 0._r8
+ write(iulog,*) 'adjust1 rsub_top',c,rsub_top(c)
+ end if
+ else
+ zwt(c) = zwt(c) - rsub_top_tot/1000._r8/rous
+ wa(c) = wa(c) + rsub_top_tot
+ end if
+ endif
+
+
+ !-- recompute jwt
+ !---------------------------------------------------------
+ ! allow jwt to equal zero when zwt is in top layer
+ jwt(c) = nlevbed
+ do j = 1,nlevbed
+ if(zwt(c) <= zi(c,j)) then
+ if (zengdecker_2009_with_var_soil_thick .and. zwt(c) == zi(c,nlevbed)) then
+ exit
+ else
+ jwt(c) = j-1
+ exit
+ end if
+ end if
+ enddo
+ end if! end of jwt if construct
+
+ zwt(c) = max(0.0_r8,zwt(c))
+ zwt(c) = min(80._r8,zwt(c))
+
+ end do
+
+ end if
+ !====================================================================================================================
+
+
+ ! excessive water above saturation added to the above unsaturated layer
+ ! like a bucket
+ ! if column fully saturated, excess water goes to runoff
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ do j = nlevbed,2,-1
+ xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8)
+ if (use_vsfm) then
+ xsi(c) = 0._r8
+ else
+ h2osoi_liq(c,j) = min(eff_porosity(c,j)*dzmm(c,j),h2osoi_liq(c,j))
+ h2osoi_liq(c,j-1) = h2osoi_liq(c,j-1) + xsi(c)
+ endif
+ end do
+ end do
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+
+ !scs: watmin addition to fix water balance errors
+ xs1(c) = max(max(h2osoi_liq(c,1)-watmin,0._r8)- &
+ max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1)-watmin)),0._r8)
+ if (use_vsfm) xs1(c) = 0._r8
+ h2osoi_liq(c,1) = h2osoi_liq(c,1) - xs1(c)
+
+ if (lun_pp%urbpoi(col_pp%landunit(c))) then
+ qflx_rsub_sat(c) = xs1(c) / dtime
+ else
+ if(h2osfcflag == 1) then
+ ! send this water up to h2osfc rather than sending to drainage
+ h2osfc(c) = h2osfc(c) + xs1(c)
+ qflx_rsub_sat(c) = 0._r8
+ else
+ ! use original code to send water to drainage (non-h2osfc case)
+ qflx_rsub_sat(c) = xs1(c) / dtime
+ endif
+ endif
+
+ if (use_vsfm) qflx_rsub_sat(c) = 0._r8
+
+ ! add in ice check
+ xs1(c) = max(max(h2osoi_ice(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1))),0._r8)
+ h2osoi_ice(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1)), h2osoi_ice(c,1))
+ qflx_snwcp_ice(c) = qflx_snwcp_ice(c) + xs1(c) / dtime
+ end do
+
+ ! Limit h2osoi_liq to be greater than or equal to watmin.
+ ! Get water needed to bring h2osoi_liq equal watmin from lower layer.
+ ! If insufficient water in soil layers, get from aquifer water
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ do j = 1, nlevbed-1
+ if (h2osoi_liq(c,j) < watmin) then
+ xs(c) = watmin - h2osoi_liq(c,j)
+ ! deepen water table if water is passed from below zwt layer
+ if(j == jwt(c)) then
+ zwt(c) = zwt(c) + xs(c)/eff_porosity(c,j)/1000._r8
+ endif
+ else
+ xs(c) = 0._r8
+ end if
+ h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c)
+ h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c)
+ end do
+ end do
+
+ ! Get water for bottom layer from layers above if possible
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+ nlevbed = nlev2bed(c)
+ j = nlevbed
+ if (h2osoi_liq(c,j) < watmin) then
+ xs(c) = watmin-h2osoi_liq(c,j)
+ searchforwater: do i = nlevbed-1, 1, -1
+ available_h2osoi_liq = max(h2osoi_liq(c,i)-watmin-xs(c),0._r8)
+ if (available_h2osoi_liq >= xs(c)) then
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
+ h2osoi_liq(c,i) = h2osoi_liq(c,i) - xs(c)
+ xs(c) = 0._r8
+ exit searchforwater
+ else
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + available_h2osoi_liq
+ h2osoi_liq(c,i) = h2osoi_liq(c,i) - available_h2osoi_liq
+ xs(c) = xs(c) - available_h2osoi_liq
+ end if
+ end do searchforwater
+ else
+ xs(c) = 0._r8
+ end if
+ ! Needed in case there is no water to be found
+ h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c)
+ ! Instead of removing water from aquifer where it eventually
+ ! shows up as excess drainage to the ocean, take it back out of
+ ! drainage
+ rsub_top(c) = rsub_top(c) - xs(c)/dtime
+
+ end do
+
+ ! add saturation-excess runoff from h3d to qflx_rsub_sat
+ if (use_h3d) then
+ do fc = 1, num_h3dc
+ c = filter_h3dc(fc)
+ qflx_rsub_sat(c) = qflx_rsub_sat(c) + qflx_rsub_sat_h3d(c)
+ end do
+ end if
+
+ do fc = 1, num_hydrologyc
+ c = filter_hydrologyc(fc)
+
+ ! Sub-surface runoff and drainage
+
+ qflx_drain(c) = qflx_rsub_sat(c) + rsub_top(c)
+
+ ! Set imbalance for snow capping
+
+ qflx_qrgwl(c) = qflx_snwcp_liq(c)
+ end do
+
+ ! No drainage for urban columns (except for pervious road as computed
+ ! above)
+
+ do fc = 1, num_urbanc
+ c = filter_urbanc(fc)
+ if (col_pp%itype(c) /= icol_road_perv) then
+ qflx_drain(c) = 0._r8
+ ! This must be done for roofs and impervious road (walls will be
+ ! zero)
+ qflx_qrgwl(c) = qflx_snwcp_liq(c)
+ end if
+ end do
+
+
+ !map zwt from col_pp to lun_pp for output
+ do fc = 1,num_h3dc,nh3dc_per_lunit !loop for all soil columns that belong to same land unit
+ c0 = filter_h3dc(fc)
+ l = col_pp%landunit(c0)
+
+ do k = 1,nh3dc_per_lunit
+ h3d_zwt_lun(l,k) = zwt(c0+k-1)
+ end do
+ end do
+
+ !double check if negtive drainage exists
+ if (use_h3d) then
+ do fc = 1,num_h3dc,nh3dc_per_lunit !loop for all soil columns that belong to same land unit
+ c0 = filter_h3dc(fc)
+ l = col_pp%landunit(c0)
+
+ tmp_sum_drain1 = 0._r8
+ tmp_sum_drain2 = 0._r8
+ do k = 1,nh3dc_per_lunit
+ c = c0+k-1
+ tmp_sum_drain1 = tmp_sum_drain1 + qflx_drain(c)*col_pp%wtgcell(c)
+ tmp_sum_drain2 = tmp_sum_drain2 + qflx_drain(c)*col_pp%wtgcell(c)
+
+ if (isnan(qflx_drain(c))) then
+ write(iulog,*) 'nan qflx_drain',c,col_pp%itype(c)
+ end if
+
+ end do
+
+ end do
+ end if
+
+ end associate
+
+ end subroutine DrainageH3D
+
+
+
+ !-----------------------------------------------------------------------
+ subroutine H3D_DRI(bounds, num_h3dc, filter_h3dc, jwt,zwt_h3d, tmp_rsub_top,h3d_rsub_top_max,fff,&
+ soilhydrology_vars, soilstate_vars, dtime)
+ !
+ ! !DESCRIPTION:
+ ! Calculate lateral subsurface drainage and water table variation using h3d
+ !
+ ! Note:
+ ! hillslope width function (w~x) is defined at lateral conputational node,
+ ! i.e., center of each soil column
+ !
+ ! !USES:
+ use clm_time_manager , only : get_step_size
+ use clm_varpar , only : nlevsoi, nlevgrnd, nlayer, nlayert
+ use clm_varpar , only : nh3dc_per_lunit !yhfang
+ use clm_varcon , only : pondmx, tfrz, watmin,rpi, secspday, nlvic
+ use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv
+ use abortutils , only : endrun
+ use clm_varctl , only : use_vsfm, use_var_soil_thick
+ use SoilWaterMovementMod, only : zengdecker_2009_with_var_soil_thick
+ use pftvarcon , only : rsub_top_globalmax
+ use spmdMod , only : masterproc, iam, npes, mpicom, comp_id
+ !
+ ! !ARGUMENTS:
+ type(bounds_type) , intent(in) :: bounds
+ integer , intent(in) :: num_h3dc ! number of column h3d points in column filter
+ integer , intent(in) :: filter_h3dc(:) ! columnfilter for h3d points
+ integer , intent(in) :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-)
+ real(r8) , intent(in) :: tmp_rsub_top(bounds%begc:bounds%endc)
+ real(r8) , intent(in) :: h3d_rsub_top_max(bounds%begc:bounds%endc)
+ real(r8) , intent(in) :: fff(bounds%begc:bounds%endc) !decay factor (m-1)
+ type(soilstate_type) , intent(in) :: soilstate_vars
+ type(soilhydrology_type) , intent(inout) :: soilhydrology_vars
+ real(r8) , intent(inout) :: zwt_h3d(bounds%begc:bounds%endc) ! (bounds%begc:bounds%endc)
+ real(r8) , intent(in) :: dtime
+ !
+ ! !LOCAL VARIABLES:
+ character(len=32) :: subname = 'H3D_DRI' ! subroutine name
+ integer :: c,i,j,k,l,fc ! indices
+ integer :: c0
+ integer :: h3d_begc,h3d_endc
+ integer :: h3dc
+ logical :: iter_conv
+
+ real(r8) :: dt_h3d !h3d midel time step (sec)
+ real(r8) :: rsub_top_default !default rsub_top using SIMTOP scheme [m/s]
+
+ real(r8) :: h_sat_prev(bounds%begc:bounds%endc) !Saturated zone level at previous h3d iteration step [m]
+ real(r8) :: h_sat_old (bounds%begc:bounds%endc) !Saturated zone level at previous h3d timestep [m]
+ real(r8) :: h_sat_begin (bounds%begc:bounds%endc) !Saturated zone level at previous elm timestep [m]
+ real(r8) :: h_sat (bounds%begc:bounds%endc) !Saturated zone level at current iteration step [m]
+
+ real(r8) :: hs_R_sub (bounds%begc:bounds%endc) !Subsurface runoff of h3d soil column [m]
+ real(r8) :: hs_Q_sub (bounds%begc:bounds%endc) !Subsurface runoff of h3d soil column [m/s]
+ real(r8) :: hs_R_of (bounds%begc:bounds%endc) !Soil saturation-excess of h3d soil column [m]
+ real(r8) :: hs_Q_of (bounds%begc:bounds%endc) !Soil saturation-excess of h3d soil column [m/s]
+ real(r8) :: hs_dS_sat (bounds%begc:bounds%endc) !Saturated storage change of h3d soil column [m]
+
+ real(r8) :: hs_dS_sat_tot !Total Ssturated storage change of h3d soil column in same land unit [m]
+ real(r8) :: dt_h3d_tot(bounds%begc:bounds%endc) !accumulated h3d iteration timestep (s) [m]
+
+
+ real(r8) :: zwt_h3d_avg,h3d_rsub_top_max_avg,fff_avg
+ real(r8) :: tmp_sum_drain
+
+ !-----------------------------------------------------------------------
+ associate( &
+ hs_dA => lun_pp%hs_dA , & ! Input: [real(r8) (:,:) ] surface area of h3d soil column (m2)
+ hs_area => lun_pp%hs_area , & ! Input: [real(r8) (:) ] surface area of h3d hillslope (m2) !
+ zibed => col_pp%zibed , & ! Input: [real(r8) (:) ] bedrock depth in model (interface level at nlevbed)
+ f_drain => col_pp%f_drain , & ! Inout: [real(r8) (:) ] drainable prosity, now = specific prosity s_y
+ dt_h3d => col_pp%dt_h3d , & ! Inout: [real(r8) (:) ] h3d iteration timestep (s)
+ zwt => soilhydrology_vars%zwt_col , & ! Inout: [real(r8) (:) ] water table depth (m)
+ zwtbed => soilhydrology_vars%zwtbed_h3d_col , & ! Inout: [real(r8) (:) ] col max. zwt allowed" zibed if var_soil_thichness; 25
+! otherwise
+ qflx_rsub_sat_h3d => col_wf%qflx_rsub_sat_h3dcol , & ! Output: [real(r8) (:) ] soil saturation excess from h3d [mm h2o/s]
+ qflx_drain_h3d => col_wf%qflx_drain_h3dcol , & ! Output: [real(r8) (:) ] sub-surface runoff from h3d(mm H2O /s)
+ qflx_drain => col_wf%qflx_drain_col & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s)
+ )
+
+
+ !initial set
+ do fc = 1,num_h3dc
+ c = filter_h3dc(fc)
+ dt_h3d(c) = dtime
+ hs_R_sub(c) = 0._r8
+ hs_Q_sub(c) = 0._r8
+ hs_R_of(c) = 0._r8
+ hs_Q_of(c) = 0._r8
+ hs_dS_sat(c) = 0._r8
+ dt_h3d_tot(c) = 0._r8
+
+ zwt_h3d(c) = zwt(c)
+
+ if (use_var_soil_thick) then
+ zwtbed(c) = zibed(c)
+ else
+ zwtbed(c) = 25._r8 + zibed(c)
+ end if
+
+ h_sat_begin(c) = zwtbed(c) - zwt_h3d(c)
+
+ end do
+!-----------------------------------------------------------------------------------------------------
+!(1)update h_sat (at h3d timestep)
+!-----------------------------------------------------------------------------------------------------
+ do fc = 1,num_h3dc,nh3dc_per_lunit !loop for all soil columns that belong to same land unit
+ c0 = filter_h3dc(fc)
+ h3d_begc = c0
+ h3d_endc = c0+nh3dc_per_lunit-1
+ l = col_pp%landunit(c0)
+
+ if (lun_pp%hs_area(l) == 0._r8) cycle
+
+ rsub_top_default = 0._r8
+ zwt_h3d_avg = 0._r8
+ h3d_rsub_top_max_avg = 0._r8
+ fff_avg = 0._r8
+ do c = h3d_begc,h3d_endc
+ rsub_top_default = rsub_top_default + tmp_rsub_top(c)*hs_dA(l,c-c0+1) !m3/s
+
+ zwt_h3d_avg = zwt_h3d_avg + zwt_h3d(c)*hs_dA(l,c-c0+1)
+ h3d_rsub_top_max_avg = h3d_rsub_top_max_avg + h3d_rsub_top_max(c)*hs_dA(l,c-c0+1)
+ fff_avg = fff_avg + fff(c)*hs_dA(l,c-c0+1)
+ end do
+ zwt_h3d_avg = zwt_h3d_avg / lun_pp%hs_area(l)
+ h3d_rsub_top_max_avg = h3d_rsub_top_max_avg / lun_pp%hs_area(l)
+ fff_avg = fff_avg / lun_pp%hs_area(l)
+
+ rsub_top_default = h3d_rsub_top_max_avg* exp(-fff_avg*zwt_h3d_avg) * lun_pp%hs_area(l)
+
+
+ do while(dt_h3d_tot(c0) < dtime)
+ dt_h3d_tot(c0) = dt_h3d_tot(c0) + dt_h3d(c0)
+
+ do c = h3d_begc,h3d_endc
+ h_sat_old(c) = zwtbed(c) - zwt_h3d(c)
+ end do
+
+ if (any(h_sat_old(h3d_begc:h3d_endc) == 0.)) then
+ h_sat(h3d_begc:h3d_endc) = h_sat_old(h3d_begc:h3d_endc)
+ f_drain(h3d_begc:h3d_endc) = 0.2_r8
+ cycle
+ end if
+
+
+ call LateralResponse(soilstate_vars,soilhydrology_vars,l,c0,dt_h3d(c0),&
+ rsub_top_default,jwt(h3d_begc:h3d_endc),h_sat_old(h3d_begc:h3d_endc),h_sat(h3d_begc:h3d_endc),iter_conv)
+
+ if (iter_conv) then
+ do c = h3d_begc,h3d_endc
+ h_sat(c) = max(0._r8,h_sat(c))
+ zwt_h3d(c) = zwtbed(c) - h_sat(c)
+ end do
+ else
+ dt_h3d_tot(c0) = dt_h3d_tot(c0) - dt_h3d(c0)
+ dt_h3d(c0) = 0.5_r8 * dt_h3d(c0)
+
+ if (dt_h3d(c0) < 10._r8) write(iulog,*) 'h3d time step!!!'
+ end if
+ end do
+ end do
+
+!-----------------------------------------------------------------------------------------------------
+!(2)update q_sub (at ELM timestep)
+!-----------------------------------------------------------------------------------------------------
+ do fc = 1,num_h3dc,nh3dc_per_lunit !loop for all soil columns that belong to same land unit
+ c0 = filter_h3dc(fc)
+ h3d_begc = c0
+ h3d_endc = c0+nh3dc_per_lunit-1
+ l = col_pp%landunit(c0)
+
+ if (lun_pp%hs_area(l) == 0._r8) cycle
+
+ hs_dS_sat(h3d_begc:h3d_endc) = 0._r8
+ hs_dS_sat_tot = 0._r8
+ tmp_sum_drain = 0._r8
+ do c = h3d_endc,h3d_begc,-1
+ !relative index (1...nh3dc_per_lunit)
+ h3dc = c-c0+1
+
+ hs_dS_sat(c) = f_drain(c) * (h_sat(c) - h_sat_begin(c)) !m
+ hs_R_sub(c) = -1._r8 * hs_dS_sat(c) !m
+ hs_Q_sub(c) = hs_R_sub(c) / dtime !m/s
+
+ hs_dS_sat_tot = hs_dS_sat_tot + hs_dS_sat(c)*hs_dA(l,h3dc) !m3
+
+ qflx_drain(c) = hs_Q_sub(c) * 1000._r8 !mm/s
+ qflx_drain_h3d(c) = hs_Q_sub(c) * 1000._r8 !mm/s
+
+ tmp_sum_drain = tmp_sum_drain + qflx_drain(c) * hs_dA(l,h3dc)
+ end do
+
+ if (hs_dS_sat_tot > 0._r8) then
+ write(iulog,*) '----------negtive discharge--------',c0,l
+ write(iulog,*) 'hs_dS_sat_tot',hs_dS_sat_tot / lun_pp%hs_area(l)
+ write(iulog,*) hs_dS_sat_tot / dtime
+ write(iulog,*) 'h_sat',h_sat(h3d_begc:h3d_endc)
+ write(iulog,*) 'h_sat_bgin',h_sat_begin(h3d_begc:h3d_endc)
+ write(iulog,*) 'qflx_drain_h3d',qflx_drain_h3d(h3d_begc:h3d_endc)
+ end if
+
+ if (any(isnan(qflx_drain(h3d_begc:h3d_endc)))) then
+ write(iulog,*) '----------NaN discharge--------',c0,l,dtime
+ write(iulog,*) 'hs_dS_sat_tot',hs_dS_sat_tot / lun_pp%hs_area(l)
+ write(iulog,*) hs_dS_sat_tot / dtime
+ write(iulog,*) 'h_sat',h_sat(h3d_begc:h3d_endc)
+ write(iulog,*) 'h_sat_bgin',h_sat_begin(h3d_begc:h3d_endc)
+ write(iulog,*) 'qflx_drain_h3d',qflx_drain_h3d(h3d_begc:h3d_endc)
+ write(iulog,*) 'f_drain',f_drain(h3d_begc:h3d_endc)
+ end if
+
+ end do
+
+ end associate
+
+ end subroutine H3D_DRI
+
+ !-----------------------------------------------------------------------
+ subroutine LateralResponse(soilstate_vars,soilhydrology_vars,l,c0,dt_h3d,rsub_top_default,jwt,h_sat_old,h_sat,iter_conv)
+ !
+ ! !DESCRIPTION:
+ ! Calculate lateral subsurface drainage and water table variation using h3d
+ !
+ !
+ ! !USES:
+ use clm_time_manager , only : get_step_size
+ use clm_varcon , only : rpi
+ use clm_varpar , only : nh3dc_per_lunit
+ !
+ ! !ARGUMENTS:
+ type(soilstate_type) , intent(in) :: soilstate_vars
+ type(soilhydrology_type) , intent(in) :: soilhydrology_vars
+ integer , intent(in) :: l !landunit index
+ integer , intent(in) :: c0 !soil column index pointing to the first colum in landunit l
+ real(r8) , intent(in) :: dt_h3d !h3d time step (sec)
+ integer , intent(in) :: jwt(:) ! index of the soil layer right above the water table (-)
+ real(r8) , intent(in) :: h_sat_old(:) !Saturated zone level at previous timestep [m]
+ real(r8) , intent(out) :: h_sat (:) !Saturated zone level at current iteration [m]
+ logical , intent(out) :: iter_conv !if h3d iteration converges
+ !real(r8) , intent(in) :: hs_slope (:)
+ real(r8) , intent(in) :: rsub_top_default
+ !
+ ! !LOCAL VARIABLES:
+ character(len=32) :: subname = 'LateralResponse' ! subroutine name
+ integer :: c,i,j,k ! indices
+ real(r8) :: h_sat_prev (1:nh3dc_per_lunit) !Saturated zone level at previous iteration [m]
+ real(r8) :: amx(1:nh3dc_per_lunit) ! "a" left off diagonal of tridiagonal matrix
+ real(r8) :: bmx(1:nh3dc_per_lunit) ! "b" diagonal column for tridiagonal matrix
+ real(r8) :: cmx(1:nh3dc_per_lunit) ! "c" right off diagonal tridiagonal matrix
+ real(r8) :: rmx(1:nh3dc_per_lunit) ! "r" forcing term of tridiagonal matrix
+ real(r8) :: w_kl_h(1:nh3dc_per_lunit) ! Product of the hillslope width function, saturated conductivity and saturated zone level [m^3 s^-1]
+ real(r8) :: h_sat_thres = 1.e-4_r8 ! threshold of h_sat for iteration converge [m]
+ integer :: niter_max = 20
+ integer :: niter
+ logical :: ierror
+ integer :: idx
+ real(r8) :: f_aniso ! factor to scale hk
+ real(r8) :: idx1 ! water table layer at the bottom of the column
+ !-----------------------------------------------------------------------
+ !-----------------------------------------------------------
+ !| * | * | * |
+ ! 1 2 3
+ ! hs_w_nod(1)
+ ! hs_w_itf(2)
+ ! |---hs_dx_node(2)----|
+ !|---hs_dx(1)--|
+ !----------------------------------------------------------
+ !above are 3 h3d soil columns in vegatated land unit
+ !in this case nh3dc_per_lunit = 3
+ !node 1 is lower boundary that connects to stream/wet land
+ !node 3 is upper boundary that represents hillslope divide
+ !hs_w_nod,hs_x_nod are hillslop width function defiened at node
+ !hs_w_nod,hs_x_nod are hillslop width function defiened at interface between
+ !2 columns
+ !the 2nd index is 1...nh3dc_per_lunit
+ !w_kl_h(i) => w_kl_h at interface i-1/2
+
+
+
+ associate( &
+ hs_dx => lun_pp%hs_dx , &
+! Input [real(r8) (:,:) ]
+ hs_dx_nod => lun_pp%hs_dx_node , &
+! Input [real(r8) (:,:) ]
+ hs_w_itf => lun_pp%hs_w_itf , &
+! Input [real(r8) (:,:) ]
+ hs_x_itf => lun_pp%hs_x_itf , &
+! Input [real(r8) (:,:) ]
+ hs_w_nod => lun_pp%hs_w_nod , &
+! Input [real(r8) (:,: ]
+ hs_x_nod => lun_pp%hs_x_nod , &
+! Input [real(r8) (:,:) ]
+ hs_slope => col_pp%h3d_slope , &
+! Input: [real(r8) (:) ] gridcell topographic slope (degree)
+ zibed => col_pp%zibed , &
+! Input: [real(r8) (:) ] bedrock depth in model (interface level at nlevbed)
+ f_drain => col_pp%f_drain , &
+! Inout: [real(r8) (:) ] drainable prosity, now = specific prosity s_y
+ nlev2bed => col_pp%nlevbed , &
+! Input: [integer (:) ] number of layers to bedrock
+ h3d_imped => soilhydrology_vars%imped_h3d_col , &
+! Input : [real(r8) (:) ] scaling facfor due to frozen soil
+ zwtbed => soilhydrology_vars%zwtbed_h3d_col , &
+! Input: [real(r8) (:) ] col max. zwt allowed" zibed if var_soil_thichness;25
+! otherwisw
+ hksat => soilstate_vars%hksat_col , &
+! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s)
+ watsat => soilstate_vars%watsat_col , &
+! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
+ bsw => soilstate_vars%bsw_col , &
+! Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
+ sucsat => soilstate_vars%sucsat_col &
+! Input: [real(r8) (:,:) ] minimum soil suction (mm)
+ )
+
+ f_aniso = 100._r8
+ niter = 0
+ iter_conv = .false.
+
+ h_sat_prev(:) = h_sat_old(:)
+
+ do while ((.not. iter_conv) .and. (niter < niter_max))
+
+ niter = niter + 1
+
+ do k=1,nh3dc_per_lunit
+ c = c0+k-1
+ idx = min(jwt(k)+1,nlev2bed(c))
+ f_drain(c) = watsat(c,idx) &
+ * ( 1. - (1.+1.e3*max(0._r8,(zwtbed(c) - h_sat_prev(k))) /sucsat(c,idx))**(-1./bsw(c,idx)))
+ f_drain(c)=max(f_drain(c) ,0.02_r8)
+
+
+ if (isnan(f_drain(c))) then
+ write(iulog,*) "nan_f_drain",idx,jwt(k),nlev2bed(c)
+ write(iulog,*) "nan_f_drain",zwtbed(c),h_sat_prev(k),sucsat(c,idx),bsw(c,idx)
+ end if
+
+ end do
+
+ do k=1,nh3dc_per_lunit
+ c = c0+k-1
+ idx = min(jwt(k)+1,nlev2bed(c))
+ w_kl_h(k) = f_aniso*hs_w_nod(l,k)*hksat(c,idx)*h_sat_prev(k) / 1000._r8
+ end do
+
+ k = 1
+ c = c0+k-1
+ amx(k) = 0._r8
+ cmx(k) = -1._r8 * w_kl_h(k+1) * cos(hs_slope(c)/180._r8*rpi) * dt_h3d / (hs_dx_nod(l,k+1) * hs_dx(l,k) * hs_w_nod(l,k))
+ bmx(k) = f_drain(c) - (amx(k) + cmx(k))
+ idx1 = min(jwt(k)+1,nlev2bed(c))
+ rmx(k) = f_drain(c) * h_sat_old(k) + dt_h3d / (hs_w_nod(l,k)*hs_dx(l,k)) * &
+ (sin(hs_slope(c)/180._r8*rpi) * w_kl_h(k+1) - cos(hs_slope(c)/180._r8*rpi) / hs_dx(l,k) * hs_w_nod(l,k) * &
+ f_aniso * hksat(c,idx1) / 1000._r8 * (h_sat_prev(k))**2)
+
+
+ do k=2,nh3dc_per_lunit - 1
+ c = c0+k-1
+ amx(k) = -1._r8 * w_kl_h(k) * cos(hs_slope(c)/180._r8*rpi) *dt_h3d / (hs_dx_nod(l,k) * hs_dx(l,k) * hs_w_nod(l,k))
+ cmx(k) = -1._r8 * w_kl_h(k+1) * cos(hs_slope(c)/180._r8*rpi) *dt_h3d / (hs_dx_nod(l,k+1) * hs_dx(l,k) * hs_w_nod(l,k))
+ bmx(k) = f_drain(c) - (amx(k) + cmx(k))
+ rmx(k) = f_drain(c) * h_sat_old(k) + dt_h3d * sin(hs_slope(c)/180._r8*rpi) / &
+ (hs_w_nod(l,k)*hs_dx(l,k)) * (w_kl_h(k+1) - w_kl_h(k))
+ end do
+
+
+ k = nh3dc_per_lunit
+ c = c0+k-1
+ amx(k) = -1._r8 * w_kl_h(k) * cos(hs_slope(c)/180._r8*rpi) * dt_h3d / (hs_dx_nod(l,k) * hs_dx(l,k) * hs_w_nod(l,k))
+ cmx(k) = 0._r8
+ bmx(k) = f_drain(c) - (amx(k) + cmx(k))
+ rmx(k) = f_drain(c) * h_sat_old(k) + dt_h3d * sin(hs_slope(c)/180._r8*rpi) / &
+ (hs_w_nod(l,k)*hs_dx(l,k)) * (- w_kl_h(k))
+
+ call Tridiagonal_h3D(nh3dc_per_lunit, amx, bmx, cmx, rmx, h_sat, ierror )
+
+ if ((maxval(abs(h_sat-h_sat_prev)) < h_sat_thres) .and. .not.(ierror)) then
+ iter_conv = .true.
+ end if
+
+ if (ierror) niter = niter_max
+
+ h_sat_prev(:) = h_sat(:)
+ end do
+
+
+ end associate
+
+ end subroutine LateralResponse
+
+ !-----------------------------------------------------------------------
+
+ subroutine Tridiagonal_h3D(numf, a, b, c, r, u, ierror)
+ ! !USES:
+ use shr_kind_mod , only: r8 => shr_kind_r8
+ ! !ARGUMENTS:
+ !implicit none
+ integer , intent(in) :: numf ! dimension
+ real(r8) , intent(in) :: a(1:numf)! "a" left off diagonal of
+tridiagonal matrix [j]
+ real(r8) , intent(in) :: b(1:numf)! "b" diagonal column for
+tridiagonal matrix [j]
+ real(r8) , intent(in) :: c(1:numf)! "c" right off diagonal
+tridiagonal matrix [j]
+ real(r8) , intent(in) :: r(1:numf)! "r" forcing term of
+tridiagonal matrix [j]
+ real(r8) , intent(inout) :: u(1:numf)! solution [j]
+ logical , intent(out) :: ierror ! true if error exists
+
+ real(r8) :: gam(1:numf)! temporary
+ real(r8) :: bet ! temporary
+ integer :: j !indics
+
+ character(len=255) :: subname ='Tridiagonal_h3d'
+ !-----------------------------------------------------------------------
+ ierror = .true.
+
+ if (b(1) == 0._r8) return
+
+ bet = b(1)
+ u(1) = r(1) / bet
+
+ do j=2,numf
+ gam(j) = c(j-1) / bet
+ bet = b(j) - a(j) * gam(j)
+
+ if (bet == 0._r8) return
+
+ u(j) = (r(j) - a(j)*u(j-1)) / bet
+ end do
+
+ do j = numf-1,1,-1
+ u(j) = u(j) - gam(j+1) * u(j+1)
+ end do
+
+
+ ierror = .false.
+ end subroutine Tridiagonal_h3D
+
+
end module SoilHydrologyMod
diff --git a/components/elm/src/biogeophys/SoilHydrologyType.F90 b/components/elm/src/biogeophys/SoilHydrologyType.F90
index 203e0a7d52b1..729d17215ff9 100644
--- a/components/elm/src/biogeophys/SoilHydrologyType.F90
+++ b/components/elm/src/biogeophys/SoilHydrologyType.F90
@@ -45,6 +45,8 @@ Module SoilHydrologyType
real(r8), pointer :: fcov_col (:) => null() ! col fractional impermeable area
real(r8), pointer :: fsat_col (:) => null() ! col fractional area with water table at surface
real(r8), pointer :: h2osfc_thresh_col (:) => null() ! col level at which h2osfc "percolates" (time constant)
+ real(r8), pointer :: zwt_h3d_col (:) => null() ! col water table depth from h3d
+ real(r8), pointer :: zwtbed_h3d_col (:) => null() ! col max. zwt allowed zibed if var_soil_thickness, 25 otherwise
! VIC
real(r8), pointer :: hkdepth_col (:) => null()! col VIC decay factor (m) (time constant)
@@ -68,6 +70,10 @@ Module SoilHydrologyType
real(r8), pointer :: fover (:) => null()! decay factor for surface runoff
real(r8), pointer :: pc (:) => null()! surface water threshold probability
+ !h3D
+ real(r8), pointer :: h3d_zwt_lun (:,:) => null()! lun water table depth defined at landunit
+ real(r8), pointer :: imped_h3d_col (:) => null()! scaling factor due to frozen soil
+
contains
procedure, public :: Init
@@ -106,6 +112,7 @@ subroutine InitAllocate(this, bounds)
! !USES:
!use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
use elm_varpar , only : nlevsno, nlevgrnd
+ use elm_varpar , only : nh3dc_per_lunit
!
! !ARGUMENTS:
class(soilhydrology_type) :: this
@@ -115,11 +122,13 @@ subroutine InitAllocate(this, bounds)
integer :: begp, endp
integer :: begc, endc
integer :: begg, endg
+ integer :: begl, endl
!------------------------------------------------------------------------
begp = bounds%begp; endp= bounds%endp
begc = bounds%begc; endc= bounds%endc
begg = bounds%begg; endg= bounds%endg
+ begl = bounds%begl; endl= bounds%endl
allocate(this%frost_table_col (begc:endc)) ; this%frost_table_col (:) = spval
@@ -160,6 +169,10 @@ subroutine InitAllocate(this, bounds)
allocate(this%fover (begg:endg)) ; this%fover (:) = spval
allocate(this%pc (begg:endg)) ; this%pc (:) = spval
+ allocate(this%h3d_zwt_lun (begl:endl,nh3dc_per_lunit)) ; this%h3d_zwt_lun (:,:) = spval
+ allocate(this%imped_h3d_col (begc:endc)) ; this%imped_h3d_col (:) = spval
+ allocate(this%zwtbed_h3d_col (begc:endc)) ; this%zwtbed_h3d_col (:) = spval
+
end subroutine InitAllocate
!------------------------------------------------------------------------
@@ -168,7 +181,8 @@ subroutine InitHistory(this, bounds)
! !USES:
use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
use elm_varctl , only : create_glacier_mec_landunit, use_cn, use_lch4
- use elm_varpar , only : nlevsno, crop_prog
+ use elm_varctl , only : use_h3d
+ use elm_varpar , only : nlevsno, crop_prog, nh3dc_per_lunit
use histFileMod , only : hist_addfld1d, hist_addfld2d, no_snow_normal
!
! !ARGUMENTS:
@@ -179,11 +193,13 @@ subroutine InitHistory(this, bounds)
integer :: begp, endp
integer :: begc, endc
integer :: begg, endg
+ integer :: begl, endl
!------------------------------------------------------------------------
begp = bounds%begp; endp= bounds%endp
begc = bounds%begc; endc= bounds%endc
begg = bounds%begg; endg= bounds%endg
+ begl = bounds%begl; endl= bounds%endl
this%wa_col(begc:endc) = spval
call hist_addfld1d (fname='WA', units='mm', &
@@ -220,6 +236,13 @@ subroutine InitHistory(this, bounds)
avgflag='A', long_name='perched water table depth (vegetated landunits only)', &
ptr_col=this%zwt_perched_col, l2g_scale_type='veg')
+ if (use_h3d) then
+ this%h3d_zwt_lun(begl:endl,1:nh3dc_per_lunit) = spval
+ call hist_addfld2d (fname='ZWT_h3D', units='m', type2d='h3dc', &
+ avgflag='A', long_name='water table depth (vegetated landunits only)', &
+ ptr_lunit=this%h3d_zwt_lun, l2g_scale_type='veg')
+ end if
+
end subroutine InitHistory
!-----------------------------------------------------------------------
diff --git a/components/elm/src/biogeophys/SoilWaterMovementMod.F90 b/components/elm/src/biogeophys/SoilWaterMovementMod.F90
index 6d57e7fc0c5d..1ff0fd05ca54 100644
--- a/components/elm/src/biogeophys/SoilWaterMovementMod.F90
+++ b/components/elm/src/biogeophys/SoilWaterMovementMod.F90
@@ -760,6 +760,7 @@ subroutine soilwater_zengdecker2009(bounds, num_hydrologyc, filter_hydrologyc, &
s_node = max(h2osoi_vol(c,jwt(c)+1)/watsat(c,jwt(c)+1), 0.01_r8)
s1 = min(1._r8, s_node)
+ if(isnan(imped(c,jwt(c)+1))) imped(c,jwt(c)+1)=1._r8
!scs: this is the expression for unsaturated hk
ka = imped(c,jwt(c)+1)*hksat(c,jwt(c)+1) &
*s1**(2._r8*bsw(c,jwt(c)+1)+3._r8)
diff --git a/components/elm/src/biogeophys/SurfaceResistanceMod.F90 b/components/elm/src/biogeophys/SurfaceResistanceMod.F90
index b4217890fca3..7dc2eb7a92bd 100644
--- a/components/elm/src/biogeophys/SurfaceResistanceMod.F90
+++ b/components/elm/src/biogeophys/SurfaceResistanceMod.F90
@@ -71,7 +71,9 @@ subroutine calc_soilevap_stress(bounds, num_nolakec, filter_nolakec, &
!character(len=32) :: subname = 'calc_soilevap_stress' ! subroutine name
associate( &
- soilbeta => soilstate_vars%soilbeta_col & ! Output: [real(r8) (:)] factor that reduces ground evaporation
+ soilbeta => soilstate_vars%soilbeta_col, & ! Output: [real(r8) (:)] factor that reduces ground evaporation
+ dsl => soilstate_vars%dsl_col , & ! Output: [real(r8) (:)] soil dry surface layer thickness
+ soilresis=> soilstate_vars%soilresis_col & ! Output: [real(r8) (:)] soil evaporation resistance
)
!select the right method and do the calculation
@@ -81,6 +83,16 @@ subroutine calc_soilevap_stress(bounds, num_nolakec, filter_nolakec, &
call calc_beta_leepielke1992(bounds, num_nolakec, filter_nolakec, &
soilstate_vars, soilbeta(bounds%begc:bounds%endc))
+ case (sl_14)
+ call cal_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, &
+ soilstate_vars, waterstate_vars, temperature_vars, &
+ dsl(bounds%begc:bounds%endc), soilresis(bounds%begc,bounds%endc))
+
+ case (sz_09)
+ call calc_soil_resistance_sz09(bounds, num_nolakec, filter_nolakec, &
+ soilstate_vars, waterstate_vars, &
+ dsl(bounds%begc:bounds%endc), soilresis(bounds%begc:bounds%endc))
+
case default
#ifndef _OPENACC
call endrun('calc_soilevap_stress' //':: a soilevap stress function must be specified!')
@@ -249,4 +261,268 @@ function getlblcef(rho,temp)result(cc)
return
end function getlblcef
+
+!------------------------------------------------------------------------------
+ subroutine calc_soil_resistance_sl14(bounds, num_nolakec, filter_nolakec, &
+ soilstate_inst, waterstate_inst, temperature_inst, dsl, soilresis)
+ !
+ ! DESCRIPTION
+ ! compute the lee-pielke beta factor to scal actual soil evaporation from
+ ! potential evaporation
+ ! xueyanz moved sl_14 (Swenson and Lawrence, 2014) from CLM5
+ !
+ ! USES
+ use shr_kind_mod , only : r8 => shr_kind_r8
+ use shr_const_mod , only : SHR_CONST_PI
+ use shr_log_mod , only : errMsg => shr_log_errMsg
+ use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+ use decompMod , only : bounds_type
+ use clm_varcon , only : denh2o, denice
+ use landunit_varcon , only : istice_mec, istwet, istsoil, istcrop
+ use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall
+ use column_varcon , only : icol_road_imperv, icol_road_perv
+ use ColumnType , only : col_pp
+ use LandunitType , only : lun_pp
+ !
+ implicit none
+ type(bounds_type) , intent(in) :: bounds ! bounds
+ integer , intent(in) :: num_nolakec
+ integer , intent(in) :: filter_nolakec(:)
+ type(soilstate_type) , intent(in) :: soilstate_inst
+ type(waterstate_type) , intent(in) :: waterstate_inst
+ type(temperature_type), intent(in) :: temperature_inst
+ real(r8) , intent(inout) :: dsl(bounds%begc:bounds%endc)
+ real(r8) , intent(inout) :: soilresis(bounds%begc:bounds%endc)
+
+ !local variables
+ real(r8) :: aird, eps, dg, d0, vwc_liq
+ real(r8) :: eff_por_top
+ integer :: c, l, fc !indices
+
+ SHR_ASSERT_ALL((ubound(dsl) == (/bounds%endc/)), errMsg(sourcefile,
+__LINE__))
+ SHR_ASSERT_ALL((ubound(soilresis) == (/bounds%endc/)),
+errMsg(sourcefile, __LINE__))
+
+ associate( &
+ dz => col_pp%dz , & !
+Input: [real(r8) (:,:) ] layer thickness (m)
+ watsat => soilstate_inst%watsat_col , & ! Input:
+[real(r8) (:,:)] volumetric soil water at saturation (porosity)
+ bsw => soilstate_inst%bsw_col , & !
+Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
+ sucsat => soilstate_inst%sucsat_col , & !
+Input: [real(r8) (:,:) ] minimum soil suction (mm)
+! eff_porosity => soilstate_inst%eff_porosity_col , & !
+! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice
+ t_soisno => temperature_inst%t_soisno_col , & !
+Input: [real(r8) (:,:) ] soil temperature (Kelvin)
+
+ h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input:
+[real(r8) (:,:)] ice lens (kg/m2)
+ h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input:
+[real(r8) (:,:)] liquid water (kg/m2)
+ )
+
+!xueyanz
+!open(12,file='t.txt')
+!open(13,file='eff_por_top.txt')
+!open(14,file='vwc_liq.txt')
+!open(15,file='aird.txt')
+!open(16,file='sucsat.txt')
+!open(17,file='bsw.txt')
+!open(18,file='watsat.txt')
+ do fc = 1,num_nolakec
+ c = filter_nolakec(fc)
+ l = col_pp%landunit(c)
+ if (lun_pp%itype(l)/=istwet .AND. lun_pp%itype(l)/=istice_mec) then
+ if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then
+ vwc_liq = max(h2osoi_liq(c,1),1.0e-6_r8)/(dz(c,1)*denh2o)
+! eff_porosity not calculated til SoilHydrology
+ eff_por_top = max(0.01_r8,watsat(c,1)-min(watsat(c,1),
+h2osoi_ice(c,1)/(dz(c,1)*denice)))
+
+! calculate diffusivity and air free pore space
+ aird = watsat(c,1)*(sucsat(c,1)/1.e7_r8)**(1./bsw(c,1))
+ d0 = 2.12e-5*(t_soisno(c,1)/273.15)**1.75 ![Bitelli et al., JH, 08]
+ eps = watsat(c,1) - aird
+ dg = eps*d0*(eps/watsat(c,1))**(3._r8/max(3._r8,bsw(c,1)))
+
+! dsl(c) = dzmm(c,1)*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) &
+! try arbitrary scaling (not top layer thickness)
+! dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) &
+ dsl(c) = 15._r8*max(0.001_r8,(0.8*eff_por_top - vwc_liq)) &
+ ! /max(0.001_r8,(watsat(c,1)- aird))
+ /max(0.001_r8,(0.8*watsat(c,1)- aird))
+
+ dsl(c)=max(dsl(c),0._r8)
+ dsl(c)=min(dsl(c),200._r8)
+
+ soilresis(c) = dsl(c)/(dg*eps*1.e3) + 20._r8
+ soilresis(c) = min(1.e6_r8,soilresis(c))
+! write(12,*) dg*eps/d0 !xueyanz
+! write(13,*) eff_por_top
+! write(14,*) vwc_liq
+! write(15,*) aird
+! write(16,*) sucsat(c,1)
+! write(17,*) bsw(c,1)
+! write(18,*) watsat(c,1)
+ else if (col_pp%itype(c) == icol_road_perv) then
+ soilresis(c) = 1.e6_r8
+ else if (col_pp%itype(c) == icol_sunwall .or. col_pp%itype(c) ==
+icol_shadewall) then
+ soilresis(c) = 1.e6_r8
+ else if (col_pp%itype(c) == icol_roof .or. col_pp%itype(c) ==
+icol_road_imperv) then
+ soilresis(c) = 1.e6_r8
+ endif
+ else
+ soilresis(c) = 0._r8
+ endif
+ enddo
+ end associate
+ end subroutine calc_soil_resistance_sl14
+
+!------------------------------------------------------------------------------
+ function do_soil_resistance_sl14()result(lres)
+ !
+ !DESCRIPTION
+ ! return true if the soil evaporative resistance is computed using a DSL
+ ! otherwise false
+ implicit none
+ logical :: lres
+
+ if(soil_stress_method==sl_14)then
+ lres=.true.
+ else
+ lres=.false.
+ endif
+ return
+
+ end function do_soil_resistance_sl14
+
+!------------------------------------------------------------------------------
+ subroutine calc_soil_resistance_sz09(bounds, num_nolakec, filter_nolakec, &
+ soilstate_inst, waterstate_inst,dsl, soilresis)
+ !
+ ! DESCRIPTION
+ ! compute the lee-pielke beta factor to scal actual soil evaporation from
+ ! potential evaporation
+ ! xueyanz implemented sz_09 (Sakaguchi and Zeng, 2009) DSL scheme here
+ !
+ ! USES
+ use shr_kind_mod , only : r8 => shr_kind_r8
+ use shr_const_mod , only : SHR_CONST_PI
+ use shr_log_mod , only : errMsg => shr_log_errMsg
+ use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+ use decompMod , only : bounds_type
+ use clm_varcon , only : denh2o, denice
+ use landunit_varcon , only : istice_mec, istwet, istsoil, istcrop
+ use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall
+ use column_varcon , only : icol_road_imperv, icol_road_perv
+ use ColumnType , only : col_pp
+ use LandunitType , only : lun_pp
+ !
+ implicit none
+ type(bounds_type) , intent(in) :: bounds ! bounds
+ integer , intent(in) :: num_nolakec
+ integer , intent(in) :: filter_nolakec(:)
+ type(soilstate_type) , intent(in) :: soilstate_inst
+ type(waterstate_type) , intent(in) :: waterstate_inst
+! type(temperature_type), intent(in) :: temperature_inst
+ real(r8) , intent(inout) :: dsl(bounds%begc:bounds%endc)
+ real(r8) , intent(inout) :: soilresis(bounds%begc:bounds%endc)
+
+ !local variables
+ real(r8) :: dg, vwc_liq
+ integer :: c, l, fc !indices
+
+ SHR_ASSERT_ALL((ubound(dsl) == (/bounds%endc/)), errMsg(sourcefile,
+__LINE__))
+ SHR_ASSERT_ALL((ubound(soilresis) == (/bounds%endc/)),
+errMsg(sourcefile, __LINE__))
+
+ associate( &
+ dz => col_pp%dz , & !
+Input: [real(r8) (:,:) ] layer thickness (m)
+ watsat => soilstate_inst%watsat_col , & ! Input:
+[real(r8) (:,:)] volumetric soil water at saturation (porosity)
+ watmin => soilstate_inst%watmin_col , & ! col
+minimum volumetric soil water (nlevsoi)
+ bsw => soilstate_inst%bsw_col , & !
+Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
+ h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input:
+[real(r8) (:,:)] ice lens (kg/m2)
+ h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input:
+[real(r8) (:,:)] liquid water (kg/m2)
+ )
+
+!xueyanz
+!open(14,file='vwc_liq.txt')
+!open(17,file='bsw.txt')
+!open(18,file='watsat.txt')
+!open(19,file='watmin.txt')
+ do fc = 1,num_nolakec
+ c = filter_nolakec(fc)
+ l = col_pp%landunit(c)
+ if (lun_pp%itype(l)/=istwet .AND. lun_pp%itype(l)/=istice_mec) then
+ if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then
+ vwc_liq = max(h2osoi_liq(c,1),1.0e-6_r8)/(dz(c,1)*denh2o)
+
+! calculate diffusivity and air free pore space
+
+ dg =
+2.2e-5*((watsat(c,1))**2.)*((1.-watmin(c,1)/watsat(c,1))**(2.+3*bsw(c,1)))
+
+! dsl(c) = dzmm(c,1)*max(0.001_r8,(0.8*eff_porosity(c,1) - vwc_liq)) &
+! try arbitrary scaling (not top layer thickness)
+! dsl(c) = 15._r8* (exp((1-min(1.,vwc_liq/watsat(c,1)))**5._r8)-1.) &
+! /(2.71828-1.)
+ dsl(c) = 50._r8* (exp((1-min(1.,vwc_liq/watsat(c,1)))**2.5_r8)-1.)
+&
+ /(2.71828-1.) !xueyanz
+ dsl(c)=max(dsl(c),0._r8)
+ dsl(c)=min(dsl(c),200._r8)
+
+ soilresis(c) = dsl(c)/(dg*1.e3)
+ soilresis(c) = min(1.e6_r8,soilresis(c))
+! write(14,*) vwc_liq
+! write(17,*) bsw(c,1)
+! write(18,*) watsat(c,1)
+! write(19,*) watmin(c,1)
+ else if (col_pp%itype(c) == icol_road_perv) then
+ soilresis(c) = 1.e6_r8
+ else if (col_pp%itype(c) == icol_sunwall .or. col_pp%itype(c) ==
+icol_shadewall) then
+ soilresis(c) = 1.e6_r8
+ else if (col_pp%itype(c) == icol_roof .or. col_pp%itype(c) ==
+icol_road_imperv) then
+ soilresis(c) = 1.e6_r8
+ endif
+ else
+ soilresis(c) = 0._r8
+ endif
+ enddo
+ end associate
+ end subroutine calc_soil_resistance_sz09
+!------------------------------------------------------------------------------
+
+ function do_soil_resistance_sz09()result(lres)
+ !
+ !DESCRIPTION
+ ! return true if the soil evaporative resistance is computed using a DSL
+ ! otherwise false
+ implicit none
+ logical :: lres
+
+ if(soil_stress_method==sz_09)then
+ lres=.true.
+ else
+ lres=.false.
+ endif
+ return
+
+ end function do_soil_resistance_sz09
+!------------------------------------------------------------------------------
+
end module SurfaceResistanceMod
diff --git a/components/elm/src/biogeophys/WaterfluxType.F90 b/components/elm/src/biogeophys/WaterfluxType.F90
index 8dfb1fb00c80..0a4c3c9cd2d5 100644
--- a/components/elm/src/biogeophys/WaterfluxType.F90
+++ b/components/elm/src/biogeophys/WaterfluxType.F90
@@ -119,6 +119,13 @@ module WaterfluxType
real(r8), pointer :: qflx_over_supply_col (:) ! col over supplied irrigation
integer , pointer :: n_irrig_steps_left_patch (:) ! number of time steps for which we still need to irrigate today (if 0, ignore)
+ ! For h3D
+ real(r8), pointer :: qflx_rsub_sat_h3dcol (:) ! h3dcol soil saturation excess [mm/s]
+ real(r8), pointer :: qflx_drain_h3dcol (:) ! h3dcol sub-surface runoff (mm H2O/s) [mm/s]
+ real(r8), pointer :: qflx_lateral_col (:) ! col lateral subsurface flux (mm H2O/s)
+ real(r8), pointer :: snow_sources_col (:) ! col snow sources (mm H2O/s)
+ real(r8), pointer :: snow_sinks_col (:) ! col snow sinks (mm H2O/s)
+
! For VSFM
real(r8), pointer :: mflx_infl_col_1d (:) ! infiltration source in top soil control volume (kg H2O /s)
real(r8), pointer :: mflx_dew_col_1d (:) ! liquid+snow dew source in top soil control volume (kg H2O /s)
diff --git a/components/elm/src/data_types/ColumnDataType.F90 b/components/elm/src/data_types/ColumnDataType.F90
index 8a027782a232..0de1d5827b6a 100644
--- a/components/elm/src/data_types/ColumnDataType.F90
+++ b/components/elm/src/data_types/ColumnDataType.F90
@@ -1746,6 +1746,12 @@ subroutine col_ws_init(this, begc, endc, h2osno_input, snow_depth_input, watsat_
this%h2osoi_vol(c,j) = 0.70_r8*watsat_input(c,j) !0.15_r8 to avoid very dry conditions that cause errors in FATES
else if (use_arctic_init) then
this%h2osoi_vol(c,j) = watsat_input(c,j) ! start saturated for arctic
+ else if (use_h3d) then
+ if (j < 5) then
+ this%h2osoi_vol(c,j) = watsat_col(c,j) * 0.6
+ else
+ this%h2osoi_vol(c,j) = watsat_col(c,j)
+ endif
else
this%h2osoi_vol(c,j) = 0.15_r8
endif
diff --git a/components/elm/src/data_types/ColumnType.F90 b/components/elm/src/data_types/ColumnType.F90
index 21e9d24bd8f4..2d0d60790ab8 100644
--- a/components/elm/src/data_types/ColumnType.F90
+++ b/components/elm/src/data_types/ColumnType.F90
@@ -58,6 +58,7 @@ module ColumnType
integer, pointer :: nlevbed (:) => null() ! number of layers to bedrock
real(r8), pointer :: zibed (:) => null() ! bedrock depth in model (interface level at nlevbed)
real(r8), pointer :: meangradz (:) => null() ! mean topographic gradient at the column level
+ real(r8), pointer :: h3d_slope (:) => null() ! gridcell topographic slope
! vertical levels
integer , pointer :: snl (:) => null() ! number of snow layers
@@ -69,6 +70,9 @@ module ColumnType
real(r8), pointer :: z_lake (:,:) => null() ! layer depth for lake (m)
real(r8), pointer :: lakedepth (:) => null() ! variable lake depth (m)
+ real(r8), pointer :: f_drain (:) => null() ! drainable porosity, now = specific porosity s_y
+ real(r8), pointer :: dt_h3d (:) => null() ! h3d iteration timestep (s)
+
! other column characteristics
logical , pointer :: hydrologically_active(:) => null() ! true if this column is a hydrologically active type
@@ -76,6 +80,10 @@ module ColumnType
logical, pointer :: is_fates(:) => null() ! True if this column is associated with a FATES active column
! False if otherwise. If fates is turned off, this array is
! all false
+
+ ! Is this an h3D column?
+ logical, pointer :: h3d_active(:) => null() ! true if this column is an h3D soil column
+
contains
procedure, public :: Init => col_pp_init
@@ -133,8 +141,12 @@ subroutine col_pp_init(this, begc, endc)
allocate(this%nlevbed (begc:endc)) ; this%nlevbed (:) = ispval
allocate(this%zibed (begc:endc)) ; this%zibed (:) = spval
allocate(this%meangradz (begc:endc)) ; this%meangradz (:) = spval
+ allocate(this%h3d_slope (begc:endc)) ; this%h3d_slope (:) = spval
+ allocate(this%f_drain (begc:endc)) ; this%f_drain (:) = spval
+ allocate(this%dt_h3d (begc:endc)) ; this%dt_h3d (:) = spval
allocate(this%hydrologically_active(begc:endc)) ; this%hydrologically_active(:) = .false.
+ allocate(this%h3d_active (begc:endc)) ; this%h3d_active (:) = .false.
! Assume that columns are not fates columns until fates initialization begins
allocate(this%is_fates(begc:endc)); this%is_fates(:) = .false.
@@ -176,8 +188,12 @@ subroutine col_pp_clean(this)
deallocate(this%nlevbed )
deallocate(this%zibed )
deallocate(this%meangradz )
+ deallocate(this%h3d_slope )
+ deallocate(this%f_drain )
+ deallocate(this%dt_h3d )
deallocate(this%hydrologically_active)
deallocate(this%is_fates)
+ deallocate(this%h3d_active )
end subroutine col_pp_clean
diff --git a/components/elm/src/data_types/LandunitType.F90 b/components/elm/src/data_types/LandunitType.F90
index 1631b4c3d5ba..1b5949c442f3 100644
--- a/components/elm/src/data_types/LandunitType.F90
+++ b/components/elm/src/data_types/LandunitType.F90
@@ -18,6 +18,7 @@ module LandunitType
!
use shr_kind_mod , only : r8 => shr_kind_r8
use elm_varcon , only : ispval, spval
+ use elm_varpar , only : nh3dc_per_lunit
!
! !PUBLIC TYPES:
implicit none
@@ -58,6 +59,16 @@ module LandunitType
real(r8), pointer :: z_0_town (:) => null() ! urban landunit momentum roughness length (m)
real(r8), pointer :: z_d_town (:) => null() ! urban landunit displacement height (m)
+ real(r8), pointer :: hs_w_itf (:,:) => null() ! hillslope width function defined at h3D soil column interface (N+1 values for N soil columns) (m)
+ real(r8), pointer :: hs_x_itf (:,:) => null() ! hillslope width function defined at h3D soil column interface (N+1 values for N soil columns) (m)
+ real(r8), pointer :: hs_w_nod (:,:) => null() ! hillslope width function defined at h3D soil column node (N+1 values for N soil columns) (m)
+ real(r8), pointer :: hs_x_nod (:,:) => null() ! hillslope width function defined at h3D soil column node (N+1 values for N soil columns) (m)
+
+ real(r8), pointer :: hs_dx (:,:) => null() ! dx of h3D soil column (m)
+ real(r8), pointer :: hs_dx_node (:,:) => null() ! dx between 2 adjunct h3D nodes (m)
+ real(r8), pointer :: hs_dA (:,:) => null() ! area of h3D soil column (m)
+ real(r8), pointer :: hs_area (:) => null() ! area of h3D hillslope (m)
+
contains
procedure, public :: Init => lun_pp_init
@@ -109,6 +120,16 @@ subroutine lun_pp_Init(this, begl, endl)
allocate(this%z_0_town (begl:endl)); this%z_0_town (:) = spval
allocate(this%z_d_town (begl:endl)); this%z_d_town (:) = spval
+ ! The following is set in initVerticalMod
+ allocate(this%hs_x_itf (begl:endl,1:nh3dc_per_lunit+1)); this%hs_x_itf(:,:) = spval
+ allocate(this%hs_w_itf (begl:endl,1:nh3dc_per_lunit+1)); this%hs_w_itf(:,:) = spval
+ allocate(this%hs_x_nod (begl:endl,1:nh3dc_per_lunit )); this%hs_x_nod(:,:) = spval
+ allocate(this%hs_w_nod (begl:endl,1:nh3dc_per_lunit )); this%hs_w_nod(:,:) = spval
+ allocate(this%hs_dx (begl:endl,1:nh3dc_per_lunit )); this%hs_dx (:,:) = spval
+ allocate(this%hs_dx_node (begl:endl,1:nh3dc_per_lunit )); this%hs_dx_node(:,:) = spval
+ allocate(this%hs_dA (begl:endl,1:nh3dc_per_lunit )); this%hs_dA (:,:) = spval
+ allocate(this%hs_area (begl:endl )); this%hs_area (:) = spval
+
end subroutine lun_pp_init
!------------------------------------------------------------------------
@@ -142,6 +163,13 @@ subroutine lun_pp_clean(this)
deallocate(this%wtlunit_roof )
deallocate(this%z_0_town )
deallocate(this%z_d_town )
+ deallocate(this%hs_x_itf )
+ deallocate(this%hs_w_itf )
+ deallocate(this%hs_x_nod )
+ deallocate(this%hs_dx )
+ deallocate(this%hs_dx_node )
+ deallocate(this%hs_dA )
+ deallocate(this%hs_area )
end subroutine lun_pp_clean
diff --git a/components/elm/src/data_types/h3DType.F90 b/components/elm/src/data_types/h3DType.F90
new file mode 100644
index 000000000000..61a6e2592d81
--- /dev/null
+++ b/components/elm/src/data_types/h3DType.F90
@@ -0,0 +1,121 @@
+module h3DType
+
+ !
+ use shr_kind_mod , only : r8 => shr_kind_r8
+ use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=)
+ use clm_varpar , only : nlevsno, nlevgrnd, nlevlak ,nh3dc_per_lunit
+ use clm_varcon , only : spval, ispval
+ use column_varcon , only : is_hydrologically_active
+ !
+ ! !PUBLIC TYPES:
+ implicit none
+ save
+ private
+ !
+ type, public :: h3d_type
+
+
+ real(r8), pointer :: qflx_evap_tot_lun(:,:) !
+ real(r8), pointer :: qflx_tran_veg_lun(:,:) !
+ real(r8), pointer :: qflx_rsub_sat_lun(:,:) !soil saturation excess (mm)
+ real(r8), pointer :: qflx_drain_lun(:,:) !sub-surface drainage (mm/s)
+ real(r8), pointer :: qflx_surf_lun(:,:) !surface runoff (mm/s)
+ real(r8), pointer :: qflx_charge_lun(:,:) !aquifer recharge rate (mm/s)
+ real(r8), pointer :: h2osfc_lun (:,:) !surface water (mm)
+ real(r8), pointer :: h2osoi_liq_lun (:,:) !soil liquid water (kg/m2)
+ real(r8), pointer :: eflx_lh_tot_lun (:,:) !total latent heat flux (W/m**2) [+ to atm]
+ real(r8), pointer :: eflx_sh_tot_lun (:,:) !total sensible heat flux (W/m**2) [+ to atm]
+
+ contains
+
+ procedure, public :: Init => h3d_vars_init
+ procedure, public :: Clean => h3d_vars_clean
+
+ end type h3d_type
+
+ type(h3d_type), public, target :: h3d_vars !h3d data structure
+ !------------------------------------------------------------------------
+
+contains
+
+ !------------------------------------------------------------------------
+ subroutine h3d_vars_init(this, begl, endl, begc, endc)
+ use histFileMod , only : hist_addfld1d,hist_addfld2d
+ !
+ ! !ARGUMENTS:
+ class(h3d_type) :: this
+ integer, intent(in) :: begl,endl,begc,endc
+ !------------------------------------------------------------------------
+
+ allocate(this%qflx_evap_tot_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_evap_tot_lun (:,:) = spval
+ allocate(this%qflx_tran_veg_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_tran_veg_lun (:,:) = spval
+ allocate(this%qflx_rsub_sat_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_rsub_sat_lun (:,:) = spval
+ allocate(this%qflx_drain_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_drain_lun (:,:) = spval
+ allocate(this%qflx_surf_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_surf_lun (:,:) = spval
+ allocate(this%qflx_charge_lun (begl:endl,1:nh3dc_per_lunit)) ; this%qflx_charge_lun (:,:) = spval
+ allocate(this%h2osfc_lun (begl:endl,1:nh3dc_per_lunit)) ; this%h2osfc_lun (:,:) = spval
+ allocate(this%h2osoi_liq_lun (begl:endl,1:nh3dc_per_lunit)) ; this%h2osoi_liq_lun (:,:) = spval
+ allocate(this%eflx_lh_tot_lun (begl:endl,1:nh3dc_per_lunit)) ; this%eflx_lh_tot_lun (:,:) = spval
+ allocate(this%eflx_sh_tot_lun (begl:endl,1:nh3dc_per_lunit)) ; this%eflx_sh_tot_lun (:,:) = spval
+
+ call hist_addfld2d (fname='h3d_evap_tot', units='mm SH2O/s', type2d='h3dc', &
+ avgflag='A', long_name='h3d total evapotranspiration (vegetated landunits only)', &
+ ptr_lunit=this%qflx_evap_tot_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_tran_veg', units='mm SH2O/s', type2d='h3dc', &
+ avgflag='A', long_name='h3d vegetation transpiration (vegetated landunits only)', &
+ ptr_lunit=this%qflx_tran_veg_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_rsub_sat', units='mm SH2O/s', type2d='h3dc', &
+ avgflag='A', long_name='soil saturation excess (vegetated landunits only)', &
+ ptr_lunit=this%qflx_rsub_sat_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_h2osfc', units='mm SH2O', type2d='h3dc', &
+ avgflag='A', long_name='soil saturation excess (vegetated landunits only)', &
+ ptr_lunit=this%h2osfc_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_h2osoi_liq', units='mm SH2O', type2d='h3dc', &
+ avgflag='A', long_name='soil liquid water (vegetated landunits only)', &
+ ptr_lunit=this%h2osoi_liq_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_lh_tot', units='W/m2', type2d='h3dc', &
+ avgflag='A', long_name='total latent heat flux (vegetated landunits only)', &
+ ptr_lunit=this%eflx_lh_tot_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_sh_tot', units='W/m2', type2d='h3dc', &
+ avgflag='A', long_name='total sensible heat flux (vegetated landunits only)', &
+ ptr_lunit=this%eflx_sh_tot_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_qdrain', units='mm/s', type2d='h3dc', &
+ avgflag='A', long_name='sub-surface drainage', &
+ ptr_lunit=this%qflx_drain_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_qsurf', units='mm/s', type2d='h3dc', &
+ avgflag='A', long_name='surface runoff', &
+ ptr_lunit=this%qflx_surf_lun, l2g_scale_type='veg')
+
+ call hist_addfld2d (fname='h3d_qcharge', units='mm/s', type2d='h3dc', &
+ avgflag='A', long_name='aquifer recharge rate (vegetated landunits only)', &
+ ptr_lunit=this%qflx_charge_lun, l2g_scale_type='veg')
+
+ end subroutine h3d_vars_init
+
+ !------------------------------------------------------------------------
+ subroutine h3d_vars_clean(this)
+ !
+ ! !ARGUMENTS:
+ class(h3d_type) :: this
+ !------------------------------------------------------------------------
+
+ deallocate(this%qflx_evap_tot_lun )
+ deallocate(this%qflx_tran_veg_lun )
+ deallocate(this%qflx_rsub_sat_lun )
+ deallocate(this%h2osfc_lun )
+ deallocate(this%h2osoi_liq_lun )
+ deallocate(this%eflx_lh_tot_lun )
+ deallocate(this%eflx_sh_tot_lun )
+
+ end subroutine h3d_vars_clean
+
+
+end module h3DType
diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90
index 55eddd299be1..8034e9ab458d 100755
--- a/components/elm/src/main/controlMod.F90
+++ b/components/elm/src/main/controlMod.F90
@@ -44,7 +44,7 @@ module controlMod
use FanMod , only: nh4_ads_coef
use AllocationMod , only: nu_com_phosphatase,nu_com_nfix
use elm_varctl , only: nu_com, use_var_soil_thick
- use elm_varctl , only: use_lake_wat_storage
+ use elm_varctl , only: use_lake_wat_storage, use_h3d
use seq_drydep_mod , only: drydep_method, DD_XLND, n_drydep
use elm_varctl , only: forest_fert_exp
use elm_varctl , only: ECA_Pconst_RGspin
@@ -573,6 +573,12 @@ subroutine control_init( )
endif
endif
+ ! h3D only works with use_var_soil_thick on
+ if (use_h3d .and. use_vichydro) then
+ call endrun(msg=' ERROR: use_h3d = .true. requires use_vichydro = .false.'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+
endif ! end of if-masterproc if-block
! ----------------------------------------------------------------------
@@ -750,6 +756,7 @@ subroutine control_spmd()
call mpi_bcast (use_vancouver, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_mexicocity, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_noio, 1, MPI_LOGICAL, 0, mpicom, ier)
+ call mpi_bcast (use_h3d, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_fan, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (fan_mode, len(fan_mode), MPI_CHARACTER, 0, mpicom, ier)
call mpi_bcast (fan_to_bgc_veg, 1, MPI_LOGICAL, 0, mpicom, ier)
diff --git a/components/elm/src/main/elm_driver.F90 b/components/elm/src/main/elm_driver.F90
index 18ecae88e37f..2365b52e3a92 100644
--- a/components/elm/src/main/elm_driver.F90
+++ b/components/elm/src/main/elm_driver.F90
@@ -189,6 +189,7 @@ module elm_driver
private :: elm_drv_patch2col
private :: elm_drv_init ! Initialization of variables needed from previous timestep
private :: write_diagnostic ! Write diagnostic information to log file
+ private :: prepare_h3d_vars ! Prepares h3D variables for output
!-----------------------------------------------------------------------
contains
@@ -1218,6 +1219,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
call HydrologyDrainage(bounds_clump, &
filter(nc)%num_nolakec, filter(nc)%nolakec, &
filter(nc)%num_hydrononsoic, filter(nc)%hydrononsoic, &
+ filter(nc)%num_h3dc, filter(nc)%h3dc, &
filter(nc)%num_urbanc, filter(nc)%urbanc, &
filter(nc)%num_do_smb_c, filter(nc)%do_smb_c, &
atm2lnd_vars, glc2lnd_vars, &
@@ -1228,6 +1230,7 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
call HydrologyDrainage(bounds_clump, &
filter(nc)%num_nolakec, filter(nc)%nolakec, &
filter(nc)%num_hydrologyc, filter(nc)%hydrologyc, &
+ filter(nc)%num_h3dc, filter(nc)%h3dc, &
filter(nc)%num_urbanc, filter(nc)%urbanc, &
filter(nc)%num_do_smb_c, filter(nc)%do_smb_c, &
atm2lnd_vars, glc2lnd_vars, &
@@ -1305,6 +1308,16 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
canopystate_vars, crop_vars)
end if
+ ! ============================================================================
+ ! Prepare h3D variables for output
+ ! ============================================================================
+ call t_startf('h3d_out')
+
+ call prepare_h3d_var(bounds_clump, filter(nc)%num_h3dc, filter(nc)%h3dc,&
+ waterflux_vars, waterstate_vars, energyflux_vars, h3d_vars)
+
+ call t_stopf('h3d_out')
+
! ============================================================================
! Check the energy and water balance, also carbon and nitrogen balance
@@ -1863,7 +1876,16 @@ subroutine elm_drv_patch2col (bounds, num_nolakec, filter_nolakec, &
call p2c (bounds, num_nolakec, filter_nolakec, &
qflx_evap_veg_patch(bounds%begp:bounds%endp), &
- qflx_evap_veg_col (bounds%begc:bounds%endc))
+ qflx_evap_veg_col (bounds%begc:bounds%endc) )
+
+ ! Averaging for patch energy flux variables
+ call p2c (bounds, num_nolakec, filter_nolakec, &
+ energyflux_vars%eflx_sh_tot_patch(bounds%begp:bounds%endp), &
+ energyflux_vars%eflx_sh_tot_col(bounds%begc:bounds%endc) )
+
+ call p2c (bounds, num_nolakec, filter_nolakec, &
+ energyflux_vars%eflx_lh_tot_patch(bounds%begp:bounds%endp), &
+ energyflux_vars%eflx_lh_tot_col(bounds%begc:bounds%endc) )
end associate
@@ -1939,4 +1961,122 @@ subroutine write_diagnostic (bounds, wrtdia, nstep, lnd2atm_vars)
end subroutine write_diagnostic
+ !------------------------------------------------------------------------
+ subroutine prepare_h3d_var(bounds_clump,num_h3dc,filter_h3dc,waterflux_vars,waterstate_vars,energyflux_vars,h3d_vars)
+ !
+ ! !DESCRIPTION:
+ ! arrange ELM column values to h3d column
+ !
+ ! !USES:
+ use abortutils , only : endrun
+ use shr_log_mod , only : errMsg => shr_log_errMsg
+ use WaterStateType , only : waterstate_type
+ use WaterFluxType , only : waterflux_type
+ use EnergyFluxType , only : energyflux_type
+ use H3dType , only : h3d_type
+ use clm_varpar , only : nlevgrnd, nh3dc_per_lunit
+ !
+ ! !ARGUMENTS:
+ type(bounds_type) , intent(in) :: bounds_clump
+ integer , intent(in) :: num_h3dc ! number of column
+h3d points in column filter
+ integer , intent(in) :: filter_h3dc(:) ! columnfilter for
+h3d points
+ type(waterstate_type) , intent(in) :: waterstate_vars
+ type(waterflux_type) , intent(in) :: waterflux_vars
+ type(energyflux_type) , intent(in) :: energyflux_vars
+ type(h3d_type ) , intent(inout) :: h3d_vars
+
+ integer :: c0,c,l,k,fc,j
+
+ associate( &
+ qflx_evap_tot_col => waterflux_vars%qflx_evap_tot_col ,&!
+Input(:)
+ qflx_tran_veg_col => waterflux_vars%qflx_tran_veg_col ,&!
+Input(:)
+ qflx_rsub_sat_col => waterflux_vars%qflx_rsub_sat_col ,&!
+Input: [real(r8) (:) ] soil saturation excess [mm h2o/s]
+ qflx_drain_col => waterflux_vars%qflx_drain_col ,&!
+Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s)
+ qflx_surf_col => waterflux_vars%qflx_surf_col ,&!
+Input: [real(r8) (:) ] surface runoff (mm H2O /s)
+ h2osfc_col => waterstate_vars%h2osfc_col ,&!
+Input: [real(r8) (:) ] surface water (mm)
+ h2osoi_liq_col => waterstate_vars%h2osoi_liq_col ,&!
+Input: [real(r8) (:,:)] liquid water (kg/m2)
+ eflx_lh_tot_col => energyflux_vars%eflx_lh_tot_col ,&!
+Input: [real(r8) (:) ] total latent heat flux (W/m**2) [+ to atm]
+ eflx_sh_tot_col => energyflux_vars%eflx_sh_tot_col ,&!
+Input: [real(r8) (:) ] total sensible heat flux (W/m**2) [+ to atm]
+ qcharge_col => soilhydrology_vars%qcharge_col ,&!
+Input: [real(r8) (:) ] aquifer recharge rate (mm h2o/s)
+
+ qflx_evap_tot_lun => h3d_vars%qflx_evap_tot_lun ,&! output
+(:,:)
+ qflx_tran_veg_lun => h3d_vars%qflx_tran_veg_lun ,&! output
+(:,:)
+ qflx_rsub_sat_lun => h3d_vars%qflx_rsub_sat_lun ,&! output:
+[real(r8) (:,:) ] soil saturation excess [mm h2o/s]
+ qflx_drain_lun => h3d_vars%qflx_drain_lun ,&! output:
+[real(r8) (:,:) ] sub-surface runoff [mm h2o/s]
+ qflx_surf_lun => h3d_vars%qflx_surf_lun ,&! output:
+[real(r8) (:,:) ] surface runoff [mm h2o/s]
+ h2osfc_lun => h3d_vars%h2osfc_lun ,&! Output:
+[real(r8) (:,:) ] surface water (mm)
+ h2osoi_liq_lun => h3d_vars%h2osoi_liq_lun ,&! Output:
+[real(r8) (:,:) ] liquid water (kg/m2)
+ eflx_lh_tot_lun => h3d_vars%eflx_lh_tot_lun ,&! Output:
+[real(r8) (:,: ] total latent heat flux (W/m**2) [+ to atm]
+ eflx_sh_tot_lun => h3d_vars%eflx_sh_tot_lun ,&! Output:
+[real(r8) (:,: ] total sensible heat flux (W/m**2) [+ to atm]
+ qflx_charge_lun => h3d_vars%qflx_charge_lun &! output:
+[real(r8) (:,:) ] aquifer recharge rate (mm h2o/s)
+ )
+
+
+ do fc=1,num_h3dc,nh3dc_per_lunit
+
+
+ c0 = filter_h3dc(fc)
+ l = col_pp%landunit(c0)
+
+ !write(iulog,*) '=====elm column => h3d column======='
+ !write(iulog,*) fc,num_h3dc,nh3dc_per_lunit
+ !write(iulog,*) c0,l
+
+ if (lun_pp%hs_area(l) == 0._r8) cycle
+
+ do k=1,nh3dc_per_lunit
+ c = c0+k-1
+ !water flux
+ qflx_evap_tot_lun(l,k) = qflx_evap_tot_col(c)
+ qflx_tran_veg_lun(l,k) = qflx_tran_veg_col(c)
+ qflx_rsub_sat_lun(l,k) = qflx_rsub_sat_col(c)
+ qflx_drain_lun (l,k) = qflx_drain_col(c)
+ qflx_surf_lun (l,k) = qflx_surf_col(c)
+ qflx_charge_lun (l,k) = qcharge_col(c)
+
+ !water state
+ h2osfc_lun(l,k) = h2osfc_col(c)
+ h2osoi_liq_lun(l,k) = 0._r8
+ do j = 1, nlevgrnd
+ h2osoi_liq_lun(l,k) = h2osoi_liq_lun(l,k) + h2osoi_liq_col(c,j)
+ end do
+
+ !energy state
+ eflx_lh_tot_lun(l,k) = eflx_lh_tot_col(c)
+ eflx_sh_tot_lun(l,k) = eflx_sh_tot_col(c)
+
+ !write(iulog,*) eflx_lh_tot_lun(l,k),eflx_lh_tot_col(c)
+ end do
+
+ !write(iulog,*) eflx_lh_tot_lun(l,:)
+
+ end do
+ end associate
+
+
+ end subroutine prepare_h3d_var
+
+
end module elm_driver
diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90
index 46477b58e4a4..8ea9357bed3b 100755
--- a/components/elm/src/main/elm_initializeMod.F90
+++ b/components/elm/src/main/elm_initializeMod.F90
@@ -35,6 +35,7 @@ module elm_initializeMod
use ColumnDataType , only : col_es
use VegetationType , only : veg_pp
use VegetationDataType , only : veg_es
+ use h3dType , only : h3d_vars
use elm_instMod
use WaterBudgetMod , only : WaterBudget_Reset
@@ -58,6 +59,7 @@ subroutine initialize1( )
!
! !USES:
use elm_varpar , only: elm_varpar_init, natpft_lb, natpft_ub
+ use elm_varpar , only: nh3dc_per_lunit
use elm_varpar , only: cft_lb, cft_ub, maxpatch_glcmec
use elm_varpar , only: mxpft, numveg, mxpft_nc, numpft
use elm_varpar , only: update_pft_array_bounds
@@ -66,6 +68,7 @@ subroutine initialize1( )
use landunit_varcon , only: landunit_varcon_init, max_lunit, istice_mec, max_polygon, max_non_poly_lunit
use column_varcon , only: col_itype_to_icemec_class
use elm_varctl , only: fsurdat, fatmlndfrc, flndtopo, fglcmask, noland, version
+ use elm_varctl , only: use_h3d
use pftvarcon , only: pftconrd
use soilorder_varcon , only: soilorder_conrd
use decompInitMod , only: decompInit_lnd, decompInit_clumps, decompInit_gtlcp
@@ -284,6 +287,7 @@ subroutine initialize1( )
allocate (wt_glc_mec (1,1,1))
allocate (topo_glc_mec(1,1,1))
endif
+ allocate (wt_h3dc (begg:endg,1:nh3dc_per_lunit ))
allocate (wt_polygon (begg:endg,1:max_topounits, max_polygon))
allocate (wt_tunit (begg:endg,1:max_topounits ))
allocate (elv_tunit (begg:endg,1:max_topounits ))
@@ -375,6 +379,10 @@ subroutine initialize1( )
! Initialize the vegetation (PFT) data types
call veg_pp%Init (bounds_proc%begp_all, bounds_proc%endp_all)
+ ! Initialize h3d_vars
+ call h3d_vars%Init (bounds_proc%begl_all, bounds_proc%endl_all, &
+ bounds_proc%begc_all, bounds_proc%endc_all)
+
! Initialize the cohort data types (nothing here yet)
! ...to be added later...
@@ -436,6 +444,8 @@ subroutine initialize1( )
deallocate (wt_cft, wt_glc_mec) !wt_lunit not deallocated because it is being used in CanopyHydrologyMod.F90
deallocate (wt_tunit, elv_tunit, slp_tunit, asp_tunit,num_tunit_per_grd)
deallocate (wt_polygon) ! RF - might be used elsewhere, not sure if we want to deallocate here.
+ !deallocate wt_h3dc
+ deallocate (wt_h3dc)
call t_stopf('elm_init1')
! initialize glc_topo
diff --git a/components/elm/src/main/elm_varctl.F90 b/components/elm/src/main/elm_varctl.F90
index 84e5c4734bfb..fd4299bb1481 100644
--- a/components/elm/src/main/elm_varctl.F90
+++ b/components/elm/src/main/elm_varctl.F90
@@ -410,6 +410,11 @@ module elm_varctl
character(len=32), public :: vsfm_satfunc_type = 'smooth_brooks_corey_bz3'
character(len=32), public :: vsfm_lateral_model_type = 'none'
+ !----------------------------------------------------------
+ ! h3D subsurface lateral flow
+ !----------------------------------------------------------
+ logical , public :: use_h3d = .false.
+
!----------------------------------------------------------
! PETSc-based thermal model switches
!----------------------------------------------------------
diff --git a/components/elm/src/main/elm_varpar.F90 b/components/elm/src/main/elm_varpar.F90
index 1f82d65ec7eb..2d3ec587147e 100644
--- a/components/elm/src/main/elm_varpar.F90
+++ b/components/elm/src/main/elm_varpar.F90
@@ -56,6 +56,10 @@ module elm_varpar
integer, parameter :: nlayer = 3 ! number of VIC soil layer --Added by AWang
integer :: nlayert ! number of VIC soil layer + 3 lower thermal layers
+ integer :: nh3dc_per_lunit=1 ! number of soil columns in landunit for h3d
+ integer :: nnode_per_nh3dc=1 ! number of lateral nodes in h3dc for solving hsb
+ integer :: h3d_hs_length=0.5_r8 ! h3d hillslope length
+
integer :: numpft = 50 ! actual # of patches (without bare), a total that spans LUs
integer :: numcft = 36 ! actual # of crops
logical :: crop_prog = .true. ! If prognostic crops is turned on
diff --git a/components/elm/src/main/elm_varsur.F90 b/components/elm/src/main/elm_varsur.F90
index 392336b7c82d..b0eeec60bcf4 100755
--- a/components/elm/src/main/elm_varsur.F90
+++ b/components/elm/src/main/elm_varsur.F90
@@ -46,7 +46,11 @@ module elm_varsur
real(r8), pointer :: wt_glc_mec(:,:,:)
! subgrid glacier_mec sfc elevation
- real(r8), pointer :: topo_glc_mec(:,:,:)
+ real(r8), pointer :: topo_glc_mec(:,:,:)
+
+ ! for h3D soil column, weight of each soil column on the hillslope (vegetated
+ ! land unit)
+ real(r8), pointer :: wt_h3dc(:,:)
! Topounit related poiters
real(r8), pointer :: num_tunit_per_grd(:) ! Topounit area fraction
diff --git a/components/elm/src/main/filterMod.F90 b/components/elm/src/main/filterMod.F90
index e1dde9445598..a25f79542dd8 100644
--- a/components/elm/src/main/filterMod.F90
+++ b/components/elm/src/main/filterMod.F90
@@ -66,6 +66,9 @@ module filterMod
integer, pointer :: hydrononsoic(:) ! non-soil hydrology filter (columns)
integer :: num_hydrononsoic ! number of columns in non-soil hydrology filter
+ integer, pointer :: h3dc(:) ! h3d filter (columns)
+ integer :: num_h3dc ! number of columns in h3d filter
+
integer, pointer :: urbanl(:) ! urban filter (landunits)
integer :: num_urbanl ! number of landunits in urban filter
integer, pointer :: nourbanl(:) ! non-urban filter (landunits)
@@ -203,6 +206,8 @@ subroutine allocFiltersOneGroup(this_filter)
allocate(this_filter(nc)%hydrologyc(bounds%endc-bounds%begc+1))
allocate(this_filter(nc)%hydrononsoic(bounds%endc-bounds%begc+1))
+ allocate(this_filter(nc)%h3dc(bounds%endc-bounds%begc+1))
+
allocate(this_filter(nc)%urbanp(bounds%endp-bounds%begp+1))
allocate(this_filter(nc)%nourbanp(bounds%endp-bounds%begp+1))
@@ -410,6 +415,21 @@ subroutine setFiltersOneGroup(bounds, this_filter, include_inactive, icemask_grc
this_filter(nc)%num_hydrologyc = f
this_filter(nc)%num_hydrononsoic = fn
+ ! Create column-level h3d filter (only soil cols)
+
+ f = 0
+ do c = bounds%begc,bounds%endc
+ if (col_pp%active(c) .or. include_inactive) then
+ l = col_pp%landunit(c)
+ if (lun_pp%itype(l) == istsoil) then
+ f = f + 1
+ this_filter(nc)%h3dc(f) = c
+ col_pp%h3d_active(c) = .true.
+ end if
+ end if
+ end do
+ this_filter(nc)%num_h3dc = f
+
! Create prognostic crop and soil w/o prog. crop filters at pft-level
! according to where the crop model should be used
diff --git a/components/elm/src/main/histFileMod.F90 b/components/elm/src/main/histFileMod.F90
index b3558d20cd76..5f66163b9a04 100644
--- a/components/elm/src/main/histFileMod.F90
+++ b/components/elm/src/main/histFileMod.F90
@@ -13,8 +13,10 @@ module histFileMod
use spmdMod , only : masterproc
use abortutils , only : endrun
use elm_varctl , only : iulog, use_vertsoilc, use_fates, use_extrasnowlayers
+ use elm_varctl , only : use_h3d
use elm_varcon , only : spval, ispval, dzsoi_decomp
use elm_varcon , only : grlnd, nameg, namet, namel, namec, namep
+ use elm_varpar , only : nh3dc_per_lunit
use decompMod , only : get_proc_bounds, get_proc_global, bounds_type
use GridcellType , only : grc_pp
use LandunitType , only : lun_pp
@@ -1960,6 +1962,10 @@ subroutine htape_create (t, histrest)
call ncd_defdim(lnfid, 'fates_levlulu', n_landuse_cats * n_landuse_cats, dimid)
end if
+ if (use_h3d) then
+ call ncd_defdim(lnfid, 'h3dc', nh3dc_per_lunit, dimid)
+ end if
+
if ( .not. lhistrest )then
call ncd_defdim(lnfid, 'hist_interval', 2, hist_interval_dimid)
call ncd_defdim(lnfid, 'time', ncd_unlimited, time_dimid)
@@ -2544,6 +2550,11 @@ subroutine htape_timeconst(t, mode)
end if
+ if (use_h3d) then
+ call ncd_defvar(varname='h3dc', xtype=ncd_int, dim1name='h3dc', &
+ long_name='soil column index of h3d hillslope', ncid=nfid(t))
+ end if
+
elseif (mode == 'write') then
if ( masterproc ) write(iulog, *) ' zsoi:',zsoi
call ncd_io(varname='levgrnd', data=zsoi, ncid=nfid(t), flag='write')
@@ -2599,6 +2610,10 @@ subroutine htape_timeconst(t, mode)
call ncd_io(varname='fates_pftmap_levcdpf',data=fates_hdim_pftmap_levcdpf, ncid=nfid(t), flag='write')
end if
+ if (use_h3d) then
+ call ncd_io(varname='h3dc',data=nh3dc_per_lunit, ncid=nfid(t), flag='write')
+ end if
+
endif
endif
@@ -4814,6 +4829,8 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out,
num2d = nelements_fates*nlevage_fates
case ('fates_levagefuel')
num2d = nlevage_fates*nfc_fates
+ case ('h3dc')
+ num2d = nh3dc_per_lunit
case('cft')
if (cft_size > 0) then
num2d = cft_size
diff --git a/components/elm/src/main/initVerticalMod.F90 b/components/elm/src/main/initVerticalMod.F90
index b24aad2d4622..633f7617477f 100755
--- a/components/elm/src/main/initVerticalMod.F90
+++ b/components/elm/src/main/initVerticalMod.F90
@@ -14,10 +14,11 @@ module initVerticalMod
use spmdMod , only : masterproc
use elm_varpar , only : more_vertlayers, nlevsno, nlevgrnd, nlevlak
use elm_varpar , only : toplev_equalspace, nlev_equalspace
- use elm_varpar , only : nlevsoi, nlevsoifl, nlevurb, nlevslp
+ use elm_varpar , only : nlevsoi, nlevsoifl, nlevurb, nlevslp
+ use elm_varpar , only : nh3dc_per_lunit, h3d_hs_length
use elm_varctl , only : fsurdat, iulog, use_var_soil_thick
use elm_varctl , only : use_vancouver, use_mexicocity, use_vertsoilc, use_extralakelayers, use_extrasnowlayers
- use elm_varctl , only : use_erosion, use_polygonal_tundra
+ use elm_varctl , only : use_erosion, use_polygonal_tundra, use_h3d
use elm_varcon , only : zlak, dzlak, zsoi, dzsoi, zisoi, dzsoi_decomp, spval, grlnd
use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
use landunit_varcon, only : istdlak, istice_mec
@@ -28,6 +29,8 @@ module initVerticalMod
use SnowHydrologyMod, only : InitSnowLayers
use ncdio_pio
use topounit_varcon , only : max_topounits
+ use landunit_varcon , only : istsoil
+ use domainMod , only : ldomain
use GridcellType , only : grc_pp
!
! !PUBLIC TYPES:
@@ -51,7 +54,7 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
real(r8) , intent(in) :: thick_roof(bounds%begl:)
!
! LOCAL VARAIBLES:
- integer :: c,l,t,ti,topi,g,i,j,lev ! indices
+ integer :: c,l,t,ti,topi,g,i,j,k,lev ! indices
type(file_desc_t) :: ncid ! netcdf id
logical :: readvar
integer :: dimid ! dimension id
@@ -60,6 +63,8 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
real(r8) ,pointer :: tslope (:) ! read in - topo_slope
real(r8) ,pointer :: gradz(:) ! read in - gradz (polygonal tundra only)
real(r8) ,pointer :: hslp_p10 (:,:,:) ! read in - hillslope slope percentiles
+ real(r8) ,pointer :: hs_w (:,:) ! read in - hillslope width function defined at soil column interface (N+1 values for N soil columns) (m)
+ real(r8) ,pointer :: hs_x (:,:) ! read in - hillslope width function defined at soil column interface (N+1 values for N soil columns) (m)
real(r8) ,pointer :: dtb (:,:) ! read in - DTB
real(r8) :: beddep ! temporary
integer :: nlevbed ! temporary
@@ -78,6 +83,8 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
real(r8), allocatable :: ziurb_wall(:,:) ! wall (layer interface)
real(r8), allocatable :: ziurb_roof(:,:) ! roof (layer interface)
real(r8) :: depthratio ! ratio of lake depth to standard deep lake depth
+ real(r8), pointer :: tmp_hs_x(:) ! local 1D array
+ real(r8), pointer :: tmp_hs_w(:) ! local 1D array
integer :: begc, endc
integer :: begl, endl
!------------------------------------------------------------------------
@@ -557,6 +564,102 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
end if
+ !--------------------------------------------------------
+ ! Read in hillslope width function for vegetated landunit
+ !--------------------------------------------------------
+
+ if (use_h3d) then
+ allocate(hs_x(bounds%begg:bounds%endg,1:nh3dc_per_lunit+1))
+ allocate(hs_w(bounds%begg:bounds%endg,1:nh3dc_per_lunit+1))
+
+ allocate(tmp_hs_x (1:nh3dc_per_lunit+1))
+ allocate(tmp_hs_w (1:nh3dc_per_lunit+1))
+
+ sum_tmp_hs_w = 0._r8
+ do i=1,nh3dc_per_lunit+1
+ tmp_hs_x(i) = h3d_hs_length / float(nh3dc_per_lunit) * float(i-1)
+ tmp_hs_w(i) = exp( tmp_hs_x(i)) !convergent
+ !tmp_hs_w(i) = exp(-tmp_hs_x(i)) !divergent
+ !tmp_hs_w(i) = 1._r8 !uniform
+ sum_tmp_hs_w = sum_tmp_hs_w + tmp_hs_w(i)
+ end do
+
+ do j=1,nh3dc_per_lunit+1
+ hs_x(bounds%begg:bounds%endg,j) = tmp_hs_x(j) * 1.e3_r8
+ hs_w(bounds%begg:bounds%endg,j) = tmp_hs_w(j) / sum_tmp_hs_w
+ end do
+
+ hs_x(bounds%begg:bounds%endg,2) = 0.5_r8 *
+hs_x(bounds%begg:bounds%endg,2)
+
+ call ncd_io(ncid=ncid, varname='hs_w', flag='read', data=hs_w,
+dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ call shr_sys_abort(' ERROR: HILLSLOPE WIDTH FUNCTION NOT on surfdata
+file'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+
+ call ncd_io(ncid=ncid, varname='hs_x', flag='read', data=hs_x,
+dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ call shr_sys_abort(' ERROR: HILLSLOPE WIDTH FUNCTION NOT on surfdata
+file'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+
+ !need test here
+ do l = begl,endl
+ if (lun_pp%itype(l) == istsoil) then !in current implementation
+hillslope only exist in vegetated land unit
+ g = lun_pp%gridcell(l)
+ lun_pp%hs_w_itf(l,:) = hs_w(g,:)
+ lun_pp%hs_x_itf(l,:) = hs_x(g,:)
+
+ lun_pp%hs_area(l) = 0._r8
+ do k=1,nh3dc_per_lunit
+ lun_pp%hs_dx(l,k) = hs_x(g,k+1) - hs_x(g,k)
+ lun_pp%hs_x_nod(l,k) = 0.5_r8*(hs_x(g,k+1) + hs_x(g,k))
+ lun_pp%hs_w_nod(l,k) = 0.5_r8*(hs_w(g,k+1) + hs_w(g,k))
+ lun_pp%hs_dA(l,k) = lun_pp%hs_w_nod(l,k) * lun_pp%hs_dx(l,k)
+ lun_pp%hs_area(l) = lun_pp%hs_area(l) + lun_pp%hs_dA(l,k)
+ end do
+
+ hs_w_scale = ldomain%area(g) * lun_pp%wtgcell(l) * 1.e6_r8 /
+lun_pp%hs_area(l)
+ lun_pp%hs_w_itf(l,:) = hs_w_scale * lun_pp%hs_w_itf(l,:)
+ lun_pp%hs_w_nod(l,:) = hs_w_scale * lun_pp%hs_w_nod(l,:)
+
+
+ write(*,*) 'read hs width function...'
+ write(*,*) l,ldomain%area(g),lun_pp%wtgcell(l),lun_pp%hs_area(l)
+
+
+ lun_pp%hs_area(l) = 0._r8
+ do k=1,nh3dc_per_lunit
+ lun_pp%hs_dA(l,k) = lun_pp%hs_w_nod(l,k) * lun_pp%hs_dx(l,k)
+ lun_pp%hs_area(l) = lun_pp%hs_area(l) + lun_pp%hs_dA(l,k)
+ end do
+
+ lun_pp%hs_dx_node(l,1) = 0.5*lun_pp%hs_dx(l,1)
+ do k=2,nh3dc_per_lunit
+ lun_pp%hs_dx_node(l,k) = lun_pp%hs_x_nod(l,k) -
+lun_pp%hs_x_nod(l,k-1)
+ end do
+
+ if (abs(lun_pp%hs_area(l) - 1.e6_r8*ldomain%area(g) *
+lun_pp%wtgcell(l))>1.e-1) then
+ call shr_sys_abort(' ERROR: INCONSISTANCE AREA OF H3D
+HILLSLOPE'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+ end if
+ enddo
+
+ deallocate(hs_x,hs_w)
+ deallocate(tmp_hs_x,tmp_hs_w)
+ endif
+
!-----------------------------------------------
! Read in topographic index and slope
!-----------------------------------------------
@@ -623,6 +726,41 @@ subroutine initVertical(bounds, snow_depth, thick_wall, thick_roof)
end do
end if
+ !----------------------------------------------------------
+ ! Read h3D slope map if available; use default slope if not
+ !----------------------------------------------------------
+
+ allocate(tslope(bounds%begg:bounds%endg))
+ call ncd_io(ncid=ncid, varname='H3D_SLOPE', flag='read', data=tslope, dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ do c = begc,endc
+ col_pp%h3d_slope(c) = col_pp%topo_slope(c)
+ end do
+ else
+ write(iulog,*) '-----------use h3d slope---------------'
+ do c = begc,endc
+ g = col_pp%gridcell(c)
+ ! check for near zero slopes, set minimum value
+ col_pp%h3d_slope(c) = max(tslope(g), 0.2_r8)
+ end do
+ endif
+ deallocate(tslope)
+
+ allocate(std(bounds%begg:bounds%endg))
+ call ncd_io(ncid=ncid, varname='STD_ELEV', flag='read', data=std,
+dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ call shr_sys_abort(' ERROR: TOPOGRAPHIC STDdev (STD_ELEV) NOT on
+surfdata file'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+ do c = begc,endc
+ g = col_pp%gridcell(c)
+ ! Topographic variables
+ col_pp%topo_std(c) = std(g)
+ end do
+ deallocate(std)
+
!-----------------------------------------------
! Read in depth to bedrock
!-----------------------------------------------
diff --git a/components/elm/src/main/surfrdMod.F90 b/components/elm/src/main/surfrdMod.F90
index ab30ae06a896..292bba6627d7 100755
--- a/components/elm/src/main/surfrdMod.F90
+++ b/components/elm/src/main/surfrdMod.F90
@@ -1069,9 +1069,12 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
! !USES:
use elm_varctl , only : create_crop_landunit, use_fates, use_polygonal_tundra
use elm_varctl , only : irrigate
+ use elm_varctl , only : use_h3d
use elm_varpar , only : surfpft_lb, surfpft_ub, surfpft_size, cft_lb, cft_ub, cft_size
use elm_varpar , only : crop_prog
+ use elm_varpar , only : nh3dc_per_lunit, h3d_hs_length
use elm_varsur , only : wt_lunit, wt_nat_patch, wt_cft, fert_cft, fert_p_cft, wt_polygon
+ use elm_varsur , only : wt_h3dc
use landunit_varcon , only : istsoil, istcrop
use landunit_varcon , only : istlowcenpoly, ilowcenpoly, istflatcenpoly, iflatcenpoly, isthighcenpoly, ihighcenpoly
use pftvarcon , only : nc3crop, nc3irrig, npcropmin
@@ -1090,6 +1093,7 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
integer ,intent(in) :: ntpu(:)
!
! !LOCAL VARIABLES:
+ integer :: g, k, i ! indices
integer :: nl, t ! index
integer :: dimid,varid ! netCDF id's
integer :: ier ! error status
@@ -1100,6 +1104,13 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
real(r8),pointer :: array2D(:,:,:) ! local 2D array
real(r8),pointer :: arrayNF(:,:,:)
real(r8),pointer :: arrayPF(:,:,:)
+ real(r8) :: hs_area_tot
+ real(r8),pointer :: hs_x(:,:) ! local 2D array
+ real(r8),pointer :: hs_w(:,:) ! local 2D array
+ real(r8),pointer :: hs_area(:) ! local 1D array
+ real(r8),pointer :: tmp_hs_x(:) ! local 1D array
+ real(r8),pointer :: tmp_hs_w(:) ! local 1D array
+ real(r8) :: sum_tmp_hs_w
character(len=32) :: subname = 'surfrd_veg_all' ! subroutine name
!-----------------------------------------------------------------------
@@ -1142,6 +1153,59 @@ subroutine surfrd_veg_all(begg, endg, ncid, ns,ntpu)
deallocate(arrayl)
+ !determin area weight of h3d soil column
+ if (.not. use_h3d) then
+ wt_h3dc(begg:endg,1:nh3dc_per_lunit) = 1._r8/nh3dc_per_lunit
+ else
+ allocate(hs_x (begg:endg,1:nh3dc_per_lunit+1))
+ allocate(hs_w (begg:endg,1:nh3dc_per_lunit+1))
+ allocate(hs_area( 1:nh3dc_per_lunit+1))
+
+ allocate(tmp_hs_x (1:nh3dc_per_lunit+1))
+ allocate(tmp_hs_w (1:nh3dc_per_lunit+1))
+
+ call ncd_io(ncid=ncid, varname='hs_w', flag='read',
+data=hs_w,dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ call endrun(' ERROR: HILLSLOPE WIDTH FUNCTION NOT on surfdata
+file'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+
+ call ncd_io(ncid=ncid, varname='hs_x', flag='read',
+data=hs_x,dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) then
+ call endrun(' ERROR: HILLSLOPE WIDTH FUNCTION NOT on surfdatafile'//&
+ errMsg(__FILE__, __LINE__))
+ end if
+
+ sum_tmp_hs_w = 0._r8
+ do i=1,nh3dc_per_lunit+1
+ tmp_hs_x(i) = h3d_hs_length / float(nh3dc_per_lunit) * float(i-1)
+ tmp_hs_w(i) = exp( tmp_hs_x(i)) !convergent
+ !tmp_hs_w(i) = exp(-tmp_hs_x(i)) !divergent
+ !tmp_hs_w(i) = 1._r8 !uniform
+ sum_tmp_hs_w = sum_tmp_hs_w + tmp_hs_w(i)
+ end do
+
+ do g=begg,endg
+ hs_area_tot = 0._r8
+ do k=1,nh3dc_per_lunit
+ hs_area(k) = (hs_x(g,k+1) - hs_x(g,k)) * &
+ (0.5_r8*(hs_w(g,k+1) + hs_w(g,k)))
+ hs_area_tot = hs_area_tot + hs_area(k)
+ end do
+
+ do k=1,nh3dc_per_lunit
+ wt_h3dc(g,k) = hs_area(k) / hs_area_tot
+ end do
+ end do
+
+ call check_sums_equal_1(wt_h3dc, begg, 'wt_h3dc', subname)
+ deallocate(hs_x,hs_w,hs_area)
+ deallocate(tmp_hs_x,tmp_hs_w)
+ end if
+
! Check the file format for CFT's and handle accordingly
call ncd_inqdid(ncid, 'cft', dimid, cft_dim_exists)
if ( cft_dim_exists .and. create_crop_landunit ) then