runs

experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl

@@ -20,7 +20,7 @@ h_stem = FT(14) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
+dt = Float64(450)
 
 # Number of timesteps between saving output:
-n = 45
+n = 4

experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl

@@ -20,7 +20,7 @@ h_leaf = FT(9.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(120)
+dt = Float64(900)
 
 # Number of timesteps between saving output:
-n = 15
+n = 2

experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl

@@ -20,7 +20,7 @@ h_stem = FT(7.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
+dt = Float64(450)
 
 # Number of timesteps between saving output:
-n = 45
+n = 4

experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl

@@ -20,7 +20,7 @@ h_stem = FT(0) # m
 t0 = Float64(21 * 3600 * 24)# start day 21 of the year
 
 # Time step size:
-dt = Float64(20)
+dt = Float64(900)
 
 # Number of timesteps between saving output:
-n = 45
+n = 1

experiments/integrated/fluxnet/fluxnet_simulation.jl

@@ -13,7 +13,7 @@ timestepper = CTS.ARS343();
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
-        update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
+        max_iters = 2,
+        update_j = CTS.UpdateEvery(CTS.NewTimeStep),
     ),
 );

experiments/integrated/performance/profile_allocations.jl

@@ -332,15 +332,15 @@ imp_tendency! = make_imp_tendency(land);
 tendency_jacobian! = make_tendency_jacobian(land);
 
 # Set up timestepping and simulation callbacks
-dt = Float64(150)
+dt = Float64(900)
 tf = Float64(t0 + 100dt)
 saveat = Array(t0:dt:tf)
 updateat = Array(t0:dt:tf)
 timestepper = CTS.ARS343()
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 2,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );

src/standalone/Vegetation/canopy_boundary_fluxes.jl

@@ -173,7 +173,8 @@ function canopy_boundary_fluxes!(
     shf .= canopy_tf.shf
     lhf .= canopy_tf.lhf
     r_ae .= canopy_tf.r_ae
-
+    p.canopy.energy.∂LHF∂qc .= canopy_tf.∂LHF∂qc
+    p.canopy.energy.∂SHF∂Tc .= canopy_tf.∂SHF∂Tc
     # Transpiration is per unit ground area, not leaf area (mult by LAI)
     fa.:($i_end) .= PlantHydraulics.transpiration_per_ground_area(
         canopy.hydraulics.transpiration,
@@ -334,7 +335,7 @@ function canopy_turbulent_fluxes_at_a_point(
     T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in)
     ΔT = T_in - T_sfc
     r_ae::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
-    ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in)
+    ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_sfc)
     ustar::FT = conditions.ustar
     r_b::FT = FT(1 / 0.01 * (ustar / 0.04)^(-1 / 2)) # CLM 5, tech note Equation 5.122
     leaf_r_b = r_b / LAI
@@ -344,5 +345,14 @@ function canopy_turbulent_fluxes_at_a_point(
     Ẽ = E / _ρ_liq
     H = -ρ_air * cp_m * ΔT / (canopy_r_b + r_ae) # CLM 5, tech note Equation 5.88, setting H_v = H and solving to remove T_s
     LH = _LH_v0 * E
-    return (lhf = LH, shf = H, vapor_flux = Ẽ, r_ae = r_ae)
+    ∂LHF∂qc = ρ_air * _LH_v0 / (leaf_r_b + leaf_r_stomata + r_ae)
+    ∂SHF∂Tc = ρ_air * cp_m / (canopy_r_b + r_ae)
+    return (
+        lhf = LH,
+        shf = H,
+        vapor_flux = Ẽ,
+        r_ae = r_ae,
+        ∂LHF∂qc = ∂LHF∂qc,
+        ∂SHF∂Tc = ∂SHF∂Tc,
+    )
 end

src/standalone/Vegetation/canopy_energy.jl

@@ -10,12 +10,28 @@ abstract type AbstractCanopyEnergyModel{FT} <: AbstractCanopyComponent{FT} end
 ClimaLand.name(model::AbstractCanopyEnergyModel) = :energy
 
 
-ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) =
-    (:shf, :lhf, :fa_energy_roots, :r_ae)
+ClimaLand.auxiliary_vars(model::AbstractCanopyEnergyModel) = (
+    :shf,
+    :lhf,
+    :fa_energy_roots,
+    :r_ae,
+    :∂SHF∂Tc,
+    :∂LHF∂qc,
+    :∂LW_n∂Tc,
+    :∂qc∂Tc,
+)
 ClimaLand.auxiliary_types(model::AbstractCanopyEnergyModel{FT}) where {FT} =
-    (FT, FT, FT, FT)
-ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) =
-    (:surface, :surface, :surface, :surface)
+    (FT, FT, FT, FT, FT, FT, FT, FT)
+ClimaLand.auxiliary_domain_names(model::AbstractCanopyEnergyModel) = (
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+    :surface,
+)
 
 """
     PrescribedCanopyTempModel{FT} <: AbstractCanopyEnergyModel{FT}
@@ -161,7 +177,53 @@ function ClimaLand.make_update_jacobian(
 
         # The derivative of the residual with respect to the prognostic variable
         ∂Tres∂T = matrix[@name(canopy.energy.T), @name(canopy.energy.T)]
-        @. ∂Tres∂T = dtγ *  MatrixFields.DiagonalMatrixRow(Y.canopy.energy.T) - (I,)
+        ∂LHF∂qc = p.canopy.energy.∂LHF∂qc
+        ∂SHF∂Tc = p.canopy.energy.∂SHF∂Tc
+        ∂LW_n∂Tc = p.canopy.energy.∂LW_n∂Tc
+        ∂qc∂Tc = p.canopy.energy.∂qc∂Tc
+        ϵ_c = p.canopy.radiative_transfer.ϵ
+        area_index = p.canopy.hydraulics.area_index
+        ac_canopy = model.parameters.ac_canopy
+        earth_param_set = canopy.parameters.earth_param_set
+        _σ = FT(LP.Stefan(earth_param_set))
+        @. ∂LW_n∂Tc = -2 * 4 * _σ * ϵ_c * Y.canopy.energy.T^3 # ≈ ϵ_ground = 1
+        @. ∂qc∂Tc = partial_q_sat_partial_T_liq(p.drivers.P, Y.canopy.energy.T)# use atmos air pressure as approximation for surface air pressure
+        @. ∂Tres∂T =
+            dtγ * MatrixFields.DiagonalMatrixRow(
+                (∂LW_n∂Tc - ∂SHF∂Tc - ∂LHF∂qc * ∂qc∂Tc) /
+                (ac_canopy * max(area_index.leaf + area_index.stem, eps(FT))),
+            ) - (I,)
     end
     return update_jacobian!
 end
+
+"""
+    partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT}
+
+Computes the quantity ∂q_sat∂T at temperature T and pressure P,
+over liquid water. 
+
+Uses the polynomial approximation from Flatau et al. (1992).
+"""
+function partial_q_sat_partial_T_liq(P::FT, T::FT) where {FT}
+    esat = FT(
+        6.11213476e2 +
+        4.44007856e1 * T +
+        1.43064234 * T^2 +
+        2.64461437e-2 * T^3 +
+        3.05903558e-4 * T^4 +
+        1.96237241e-6 * T^5 +
+        8.92344772e-9 * T^6 - 3.73208410e-11 * T^7 + 2.09339997e-14 * T^8,
+    )
+    desatdT = FT(
+        4.44017302e1 +
+        2.86064092 * T +
+        7.94683137e-2 * T^2 +
+        1.21211669e-3 * T^3 +
+        1.03354611e-5 * T^4 +
+        4.04125005e-8 * T^5 - 7.88037859e-11 * T^6 - 1.14596802e-12 * T^7 +
+        3.81294516e-15 * T^8,
+    )
+
+    return FT(0.622) * P / (P - FT(0.378) * esat)^2 * desatdT
+end

src/standalone/Vegetation/radiation.jl

@@ -236,8 +236,6 @@ function canopy_radiant_energy_fluxes!(
     @. p.canopy.radiative_transfer.SW_n =
         (energy_per_photon_PAR * N_a * APAR) +
         (energy_per_photon_NIR * N_a * ANIR)
-    LAI = p.canopy.hydraulics.area_index.leaf
-    SAI = p.canopy.hydraulics.area_index.stem
     ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI
     # Long wave: use soil conditions from the PrescribedSoil driver
     T_soil::FT = s.T(t)