first PR comment

CliMA · Sep 4, 2024 · 96c90b4 · 96c90b4
1 parent 9e72938
commit 96c90b4
Show file tree

Hide file tree

Showing 3 changed files with 69 additions and 56 deletions.
diff --git a/experiments/long_runs/land.jl b/experiments/long_runs/land.jl
@@ -67,7 +67,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         depth = depth,
         nelements = nelements,
         npolynomial = 1,
-        dz_tuple = FT.((10.0, 0.5)),# top layer should ideally be only a few cm!
+        dz_tuple = FT.((10.0, 0.05)),# top layer should ideally be only a few cm!
     )
     surface_space = domain.space.surface
     subsurface_space = domain.space.subsurface
@@ -604,7 +604,8 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         t0,
         ref_time;
         output_writer = nc_writer,
-        output_vars = :long,
+        output_vars = :short,
+        average_period = :monthly,
     )
 
     diagnostic_handler =
@@ -619,7 +620,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 7 # keep short until it runs! * 365
+    tf = 60 * 60.0 * 24 * 7
     Δt = 900.0
     nelements = (101, 15)
     if greet
@@ -648,7 +649,7 @@ setup_and_solve_problem(; greet = true);
 # read in diagnostics and make some plots!
 #### ClimaAnalysis ####
 simdir = ClimaAnalysis.SimDir(outdir)
-short_names = ["gpp", "swc", "si", "sie"]
+short_names = ["gpp", "ct", "lai", "swc", "si"]
 for short_name in short_names
     var = get(simdir; short_name)
     times = ClimaAnalysis.times(var)

diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl
@@ -425,6 +425,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         t0,
         ref_time;
         output_writer = nc_writer,
+        average_period = :monthly,
     )
 
     diagnostic_handler =
@@ -439,7 +440,7 @@ end
 function setup_and_solve_problem(; greet = false)
 
     t0 = 0.0
-    tf = 60 * 60.0 * 24 * 180 # keep short until it runs! * 365
+    tf = 60 * 60.0 * 24 * 180
     Δt = 900.0
     nelements = (101, 15)
     if greet
@@ -472,8 +473,7 @@ short_names = ["swc", "si", "sie"]
 for short_name in short_names
     var = get(simdir; short_name)
     times = ClimaAnalysis.times(var)
-    id = 1:10:length(times)
-    for t in times[id]
+    for t in times
         var = get(simdir; short_name)
         fig = CairoMakie.Figure(size = (800, 600))
         kwargs = ClimaAnalysis.has_altitude(var) ? Dict(:z => 1) : Dict()

diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl
@@ -802,10 +802,10 @@ function turbulent_fluxes(
     )
     return soil_turbulent_fluxes_at_a_point.(
         T_sfc,
-        h_sfc,
-        d_sfc,
         θ_l_sfc,
         θ_i_sfc,
+        h_sfc,
+        d_sfc,
         hydrology_cm_sfc,
         ν_sfc,
         θ_r_sfc,
@@ -824,30 +824,38 @@ function turbulent_fluxes(
 end
 
 """
-    soil_turbulent_fluxes_at_a_point(T_sfc::FT,
-                                h_sfc::FT,
-                                d_sfc::FT,
+    soil_turbulent_fluxes_at_a_point(
+                                T_sfc::FT,
                                 θ_l_sfc::FT,
                                 θ_i_sfc::FT,
+                                h_sfc::FT,
+                                d_sfc::FT,
                                 hydrology_cm_sfc::C,
-                                ν_sfc,
-                                θ_r_sfc,
-                                K_sat_sfc,
-                                ts_in,
-                                u::FT,
-                                h::FT,
+                                ν_sfc::FT,
+                                θ_r_sfc::FT,
+                                K_sat_sfc::FT,
+                                thermal_state_air::Thermodynamics.PhaseEquil{FT},
+                                u_air::FT,
+                                h_air::FT,
                                 gustiness::FT,
-                                params, P
-                               ) where {FT <: AbstractFloat, C, P}
+                                z_0m::FT,
+                                z_0b::FT,
+                                Ω::FT,
+                                γ::FT,
+                                γT_ref,::FT
+                                earth_param_set::EP
+                               ) where {FT <: AbstractFloat, C, EP}
 
 Computes turbulent surface fluxes for soil at a point on a surface given
-(1) the surface temperature (T_sfc)
-(2) Other surface properties, such as the surface hydraulic conductivity, the surface ice content, 
-    and the hydrologic parameters, as well as
-(3) the topographical height of the surface (h_sfc),  the displacement height `d_sfc`.
-(4) soil scalar parameters contained in `params`
-(5) the prescribed atmospheric state, `ts_in`, u, h the height
-    at which these measurements are made, and the gustiness parameter (m/s).
+(1) Surface state conditions (`T_sfc`, `θ_l_sfc`, `θ_i_sfc`)
+(2) Surface properties, such as the topographical height of the surface (`h_sfc`),  
+    the displacement height (`d_sfc`), hydraulic parameters (`hydrology_cm_sfc`,
+    `ν_sfc, `θ_r_sfc`, `K_sat_sfc`)
+(4) Atmospheric state conditions (`thermal_state_air`, `u_air`)
+(5) Height corresponding to where atmospheric state is measured (`h_air`)
+(6) Parameters: `gustiness`, roughness lengths `z_0m`, `z_0b`, several
+    required to compute the soil conductivity `Ω`, γ`, γT_ref`, and the
+    `earth_param_set`
 
 This returns an energy flux and a liquid water volume flux, stored in
 a tuple with self explanatory keys.
@@ -862,9 +870,9 @@ function soil_turbulent_fluxes_at_a_point(
     ν_sfc::FT,
     θ_r_sfc::FT,
     K_sat_sfc::FT,
-    ts_in::Thermodynamics.PhaseEquil{FT},
-    u::FT,
-    h::FT,
+    thermal_state_air::Thermodynamics.PhaseEquil{FT},
+    u_air::FT,
+    h_air::FT,
     gustiness::FT,
     z_0m::FT,
     z_0b::FT,
@@ -875,7 +883,7 @@ function soil_turbulent_fluxes_at_a_point(
 ) where {FT <: AbstractFloat, P}
     thermo_params = LP.thermodynamic_parameters(earth_param_set)
     # Estimate surface air density
-    ρ_sfc::FT = ClimaLand.compute_ρ_sfc(thermo_params, ts_in, T_sfc)
+    ρ_sfc::FT = ClimaLand.compute_ρ_sfc(thermo_params, thermal_state_air, T_sfc)
     # Compute saturated specific humidities at surface over ice and liquid water
     q_sat_ice::FT = Thermodynamics.q_vap_saturation_generic(
         thermo_params,
@@ -894,25 +902,25 @@ function soil_turbulent_fluxes_at_a_point(
     ts_sfc::Thermodynamics.PhaseEquil{FT} =
         Thermodynamics.PhaseEquil_ρTq(thermo_params, ρ_sfc, T_sfc, q_sat_liq)# use to get potential evaporation E0, r_ae (weak dependence on q).
 
-    # SurfaceFluxes.jl expects a relative difference between where u = 0
+    # SurfaceFluxes.jl expects a relative difference between where u_air = 0
     # and the atmosphere height. Here, we assume h and h_sfc are measured
     # relative to a common reference. Then d_sfc + h_sfc + z_0m is the apparent
     # source of momentum, and
-    # Δh ≈ h - d_sfc - h_sfc is the relative height difference between the
+    # Δh ≈ h_air - d_sfc - h_sfc is the relative height difference between the
     # apparent source of momentum and the atmosphere height.
 
     # In this we have neglected z_0m and z_0b (i.e. assumed they are small
     # compared to Δh).
     state_sfc = SurfaceFluxes.StateValues(FT(0), SVector{2, FT}(0, 0), ts_sfc)
-    state_in = SurfaceFluxes.StateValues(
-        h - d_sfc - h_sfc,
-        SVector{2, FT}(u, 0),
-        ts_in,
+    state_air = SurfaceFluxes.StateValues(
+        h_air - d_sfc - h_sfc,
+        SVector{2, FT}(u_air, 0),
+        thermal_state_air,
     )
 
     # State containers
-    sc = SurfaceFluxes.ValuesOnly(
-        state_in,
+    states = SurfaceFluxes.ValuesOnly(
+        state_air,
         state_sfc,
         z_0m,
         z_0b,
@@ -921,18 +929,22 @@ function soil_turbulent_fluxes_at_a_point(
     surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
     conditions = SurfaceFluxes.surface_conditions(
         surface_flux_params,
-        sc;
+        states;
         tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
     )
     _LH_v0::FT = LP.LH_v0(earth_param_set)
     _ρ_liq::FT = LP.ρ_cloud_liq(earth_param_set)
-    cp_m::FT = Thermodynamics.cp_m(thermo_params, ts_in)
-    T_in::FT = Thermodynamics.air_temperature(thermo_params, ts_in)
-    ΔT::FT = T_in - T_sfc
-    r_ae_liq::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
-    ρ_air::FT = Thermodynamics.air_density(thermo_params, ts_in)
-    q_air::FT = Thermodynamics.total_specific_humidity(thermo_params, ts_in)
-    E0::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch)# potential evaporation rate, mass flux
+    cp_m::FT = Thermodynamics.cp_m(thermo_params, thermal_state_air)
+    T_air::FT = Thermodynamics.air_temperature(thermo_params, thermal_state_air)
+    ρ_air::FT = Thermodynamics.air_density(thermo_params, thermal_state_air)
+    q_air::FT =
+        Thermodynamics.total_specific_humidity(thermo_params, thermal_state_air)
+
+    ΔT::FT = T_air - T_sfc
+    r_ae_liq::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states))
+
+    E0::FT =
+        SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)# potential evaporation rate, mass flux
     Ẽ0::FT = E0 / _ρ_liq
     Ẽ::FT = Ẽ0
     if q_air < q_sat_liq # adjust potential evaporation rate to account for soil resistance
@@ -952,9 +964,9 @@ function soil_turbulent_fluxes_at_a_point(
         x::FT = 4 * K_sfc * (1 + Ẽ0 / (4 * K_c))
         Ẽ *= x / (Ẽ0 + x)
     else
-        Ẽ *= 0
+        Ẽ *= 0 # condensation, set to zero
     end
-    
+
     H_liq::FT = -ρ_air * cp_m * ΔT / r_ae_liq
 
     # Sublimation:
@@ -964,28 +976,28 @@ function soil_turbulent_fluxes_at_a_point(
     if q_air < q_sat_ice # sublimation, adjust β
         β_i *= (θ_i_sfc / ν_sfc)^3
     else
-        β_i *= 0
+        β_i *= 0 # frost, set to zero
     end
 
     state_sfc = SurfaceFluxes.StateValues(FT(0), SVector{2, FT}(0, 0), ts_sfc)
-    sc = SurfaceFluxes.ValuesOnly(
-        state_in,
+    states = SurfaceFluxes.ValuesOnly(
+        state_air,
         state_sfc,
         z_0m,
         z_0b,
         beta = β_i,
         gustiness = gustiness,
     )
-    surface_flux_params = LP.surface_fluxes_parameters(earth_param_set)
     conditions = SurfaceFluxes.surface_conditions(
         surface_flux_params,
-        sc;
+        states;
         tol_neutral = SFP.cp_d(surface_flux_params) / 100000,
     )
-    S::FT = SurfaceFluxes.evaporation(surface_flux_params, sc, conditions.Ch)# sublimation rate, mass flux
+    S::FT =
+        SurfaceFluxes.evaporation(surface_flux_params, states, conditions.Ch)# sublimation rate, mass flux
     S̃::FT = S / _ρ_liq
 
-    r_ae_ice::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(sc))
+    r_ae_ice::FT = 1 / (conditions.Ch * SurfaceFluxes.windspeed(states))
     H_ice::FT = -ρ_air * cp_m * ΔT / r_ae_ice
 
     # Heat fluxes