Merge pull request #732 from CliMA/kd/evaporation_issues

New evaporation scheme

.buildkite/longruns_gpu/pipeline.yml

@@ -37,7 +37,7 @@ steps:
         artifact_paths: "land_longrun_gpu/*png"
         agents:
           slurm_gpus: 1
-          slurm_time: 01:00:00
+          slurm_time: 03:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
@@ -47,6 +47,6 @@ steps:
         artifact_paths: "soil_longrun_gpu/*png"
         agents:
           slurm_gpus: 1
-          slurm_time: 01:00:00
+          slurm_time: 03:00:00
         env:
           CLIMACOMMS_DEVICE: "CUDA"

docs/src/APIs/Soil.md

@@ -53,14 +53,6 @@ ClimaLand.Soil.phase_change_source
 ClimaLand.Soil.thermal_time
 ```
 
-## Soil Surface Parameterizations
-
-```@docs
-ClimaLand.soil.soil_resistance
-ClimaLand.Soil.dry_soil_layer_thickness
-ClimaLand.Soil.soil_tortuosity
-```
-
 ## Soil Runoff Types and Methods
 
 ```@docs

docs/tutorials/standalone/Soil/evaporation.jl

@@ -223,10 +223,8 @@ sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat);
 
 # Extract the evaporation at each saved step
 evap = [
-    parent(
-        sv.saveval[k].soil.turbulent_fluxes.vapor_flux .*
-        (1 .- sv.saveval[k].soil.ice_frac),
-    )[1] for k in 1:length(sol.t)
+    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol.t)
 ]
 savepath = joinpath(pkgdir(ClimaLand), "docs/tutorials/standalone/Soil/")
 evaporation_data =

docs/tutorials/standalone/Soil/evaporation_gilat_loess.jl

@@ -37,7 +37,7 @@ thermo_params = LP.thermodynamic_parameters(earth_param_set);
 
 # Parameters
 K_sat = FT(0.01 / 3600 / 24)
-vg_n = FT(1.45)
+vg_n = FT(1.55)
 vg_α = FT(1.5)
 hcm = vanGenuchten{FT}(; α = vg_α, n = vg_n)
 ν = FT(0.4)
@@ -256,10 +256,8 @@ driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
 cb = SciMLBase.CallbackSet(driver_cb, saving_cb)
 sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat)
 evap = [
-    parent(
-        sv.saveval[k].soil.turbulent_fluxes.vapor_flux .*
-        (1 .- sv.saveval[k].soil.ice_frac),
-    )[1] for k in 1:length(sol.t)
+    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol.t)
 ];
 
 ## Repeat with no drainage (Ksat = 0, different BC), and with evaporation, in shorter domain
@@ -323,10 +321,8 @@ cb = SciMLBase.CallbackSet(driver_cb, saving_cb)
 sol_no_drainage =
     SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat)
 evap_no_drainage = [
-    parent(
-        sv.saveval[k].soil.turbulent_fluxes.vapor_flux .*
-        (1 .- sv.saveval[k].soil.ice_frac),
-    )[1] for k in 1:length(sol.t)
+    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol.t)
 ];
 
 

docs/tutorials/standalone/Soil/sublimation.jl

@@ -209,16 +209,12 @@ sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat);
 
 # Extract the evaporation at each saved step
 evap = [
-    parent(
-        sv.saveval[k].soil.turbulent_fluxes.vapor_flux .*
-        (1 .- sv.saveval[k].soil.ice_frac),
-    )[1] for k in 1:length(sol.t)
+    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for
+    k in 1:length(sol.t)
 ]
 sub = [
-    parent(
-        sv.saveval[k].soil.turbulent_fluxes.vapor_flux .*
-        sv.saveval[k].soil.ice_frac,
-    )[1] for k in 1:length(sol.t)
+    parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice)[1] for
+    k in 1:length(sol.t)
 ]
 
 savepath = joinpath(pkgdir(ClimaLand), "docs/tutorials/standalone/Soil/")

experiments/benchmarks/bucket.jl

@@ -224,7 +224,7 @@ if ClimaComms.device() isa ClimaComms.CUDADevice
 end
 
 if get(ENV, "BUILDKITE_PIPELINE_SLUG", nothing) == "climaland-benchmark"
-    PREVIOUS_BEST_TIME = 3.6
+    PREVIOUS_BEST_TIME = 2.8
     if average_timing_s > PREVIOUS_BEST_TIME + std_timing_s
         @info "Possible performance regression, previous average time was $(PREVIOUS_BEST_TIME)"
     elseif average_timing_s < PREVIOUS_BEST_TIME - std_timing_s

experiments/integrated/fluxnet/ozark_pft.jl

@@ -455,8 +455,10 @@ T =
     ] .* (1e3 * 24 * 3600)
 E =
     [
-        parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for
-        k in 1:length(sol.t)
+        parent(
+            sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+
+            sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice,
+        )[1] for k in 1:length(sol.t)
     ] .* (1e3 * 24 * 3600)
 ET_model = T .+ E
 if drivers.LE.status == absent

experiments/integrated/fluxnet/run_fluxnet.jl

@@ -399,8 +399,10 @@ T =
     ] .* (1e3 * 24 * 3600)
 E =
     [
-        parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for
-        k in 1:length(sol.t)
+        parent(
+            sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+
+            sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice,
+        )[1] for k in 1:length(sol.t)
     ] .* (1e3 * 24 * 3600)
 ET_model = T .+ E
 if drivers.LE.status == absent

experiments/integrated/performance/conservation/ozark_conservation.jl

@@ -104,8 +104,10 @@ for float_type in (Float32, Float64)
 
         # Evaporation
         E = [
-            parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux)[1] for
-            k in 2:length(sol.t)
+            parent(
+                sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq .+
+                sv.saveval[k].soil.turbulent_fluxes.vapor_flux_ice,
+            )[1] for k in 2:length(sol.t)
         ]
         # Root sink term: a positive root extraction is a sink term for soil; add minus sign
         root_sink =

experiments/long_runs/land.jl

@@ -613,7 +613,7 @@ function setup_prob(t0, tf, Δt; outdir = outdir, nelements = (101, 15))
         t0,
         ref_time;
         output_writer = nc_writer,
-        output_vars = :long,
+        output_vars = :short,
     )
 
     diagnostic_handler =
@@ -628,7 +628,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 * 14
     Δt = 900.0
     nelements = (101, 15)
     if greet
@@ -641,12 +641,12 @@ function setup_and_solve_problem(; greet = false)
     prob, cb = setup_prob(t0, tf, Δt; nelements)
 
     # Define timestepper and ODE algorithm
-    stepper = CTS.ARS343()
+    stepper = CTS.ARS111()
     ode_algo = CTS.IMEXAlgorithm(
         stepper,
         CTS.NewtonsMethod(
-            max_iters = 1,
-            update_j = CTS.UpdateEvery(CTS.NewTimeStep),
+            max_iters = 3,
+            update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
         ),
     )
     SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb, adaptive = false)
@@ -657,7 +657,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)

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 * 60 # keep short until it runs! * 365
+    tf = 60 * 60.0 * 24 * 340
     Δt = 900.0
     nelements = (101, 15)
     if greet

lib/ClimaLandSimulations/src/utilities/climaland_output_dataframe.jl

@@ -56,7 +56,7 @@ function make_output_df(
         (1, :soil, :turbulent_fluxes, :lhf),
         collect(map(i -> (i, :soil, :T), 1:20)), # 20 shouldn't be hard-coded, but an arg, equal to n layers
         collect(map(i -> (i, :soil, :θ_l), 1:20)),
-        (1, :soil, :turbulent_fluxes, :vapor_flux),
+        (1, :soil, :turbulent_fluxes, :vapor_flux_liq),
         (1, :canopy, :sif, :SIF),
     )
 

lib/ClimaLandSimulations/src/utilities/makie_plots.jl

@@ -72,7 +72,7 @@ function timeseries_fluxes_fig(
     p_ET_m = lines!(
         ax_W,
         datetime2unix.(climaland.DateTime),
-        (climaland.vapor_flux .* 1e3 .* 24 .* 3600) .+
+        (climaland.vapor_flux_liq .* 1e3 .* 24 .* 3600) .+
         (climaland.transpiration .* 1e3 .* 24 .* 3600),
         color = :blue,
     ) # not sure about units

src/diagnostics/land_compute_methods.jl

@@ -136,7 +136,7 @@ end
 @diagnostic_compute "soil_aerodynamic_resistance" SoilCanopyModel p.soil.turbulent_fluxes.r_ae
 @diagnostic_compute "soil_latent_heat_flux" SoilCanopyModel p.soil.turbulent_fluxes.lhf
 @diagnostic_compute "soil_sensible_heat_flux" SoilCanopyModel p.soil.turbulent_fluxes.shf
-@diagnostic_compute "vapor_flux" SoilCanopyModel p.soil.turbulent_fluxes.vapor_flux
+@diagnostic_compute "vapor_flux" SoilCanopyModel p.soil.turbulent_fluxes.vapor_flux_liq # should add ice here
 
 # Soil - SoilCO2
 function compute_heterotrophic_respiration!(

src/integrated/soil_canopy_model.jl

@@ -418,11 +418,8 @@ function soil_boundary_fluxes!(
     bc = soil.boundary_conditions.top
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, soil, Y, p, t)
     Soil.Runoff.update_runoff!(p, bc.runoff, Y, t, soil)
-    # Multiply the vapor flux by 1 - p.soil.ice_frac to get
-    # the approximated evaporation of liquid water
     @. p.soil.top_bc.water =
-        p.soil.infiltration +
-        p.soil.turbulent_fluxes.vapor_flux * (1 - p.soil.ice_frac)
+        p.soil.infiltration + p.soil.turbulent_fluxes.vapor_flux_liq
     @. p.soil.top_bc.heat =
         -p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf
 end

src/standalone/Soil/Soil.jl

@@ -64,7 +64,9 @@ import ClimaCore.MatrixFields: @name, ⋅
 import ..Parameters as LP
 import ClimaCore: Fields, Operators, Geometry, Spaces
 using Thermodynamics
-
+using SurfaceFluxes
+using StaticArrays
+import SurfaceFluxes.Parameters as SFP
 import ClimaLand.Domains: Column, HybridBox, SphericalShell
 import ClimaLand:
     AbstractImExModel,
@@ -90,6 +92,7 @@ import ClimaLand:
     surface_emissivity,
     surface_height,
     surface_resistance,
+    turbulent_fluxes,
     get_drivers
 export RichardsModel,
     RichardsParameters,

src/standalone/Soil/boundary_conditions.jl

@@ -662,7 +662,6 @@ boundary_vars(
     ::ClimaLand.TopBoundary,
 ) = (
     :turbulent_fluxes,
-    :ice_frac,
     :R_n,
     :top_bc,
     :sfc_scratch,
@@ -692,7 +691,6 @@ boundary_var_domain_names(
     :surface,
     :surface,
     :surface,
-    :surface,
     Runoff.runoff_var_domain_names(bc.runoff)...,
 )
 """
@@ -717,8 +715,10 @@ boundary_var_types(
     },
     ::ClimaLand.TopBoundary,
 ) where {FT} = (
-    NamedTuple{(:lhf, :shf, :vapor_flux, :r_ae), Tuple{FT, FT, FT, FT}},
-    FT,
+    NamedTuple{
+        (:lhf, :shf, :vapor_flux_liq, :r_ae, :vapor_flux_ice),
+        Tuple{FT, FT, FT, FT, FT},
+    },
     FT,
     NamedTuple{(:water, :heat), Tuple{FT, FT}},
     FT,
@@ -764,12 +764,9 @@ function soil_boundary_fluxes!(
     p.soil.turbulent_fluxes .= turbulent_fluxes(bc.atmos, model, Y, p, t)
     p.soil.R_n .= net_radiation(bc.radiation, model, Y, p, t)
     update_runoff!(p, bc.runoff, Y, t, model)
-    # We do not model the energy flux from infiltration. We multiply
-    # the vapor flux by the ice fraction in order to get the contribution
-    # from liquid water
+    # We do not model the energy flux from infiltration.
     @. p.soil.top_bc.water =
-        p.soil.infiltration +
-        p.soil.turbulent_fluxes.vapor_flux * (1 - p.soil.ice_frac)
+        p.soil.infiltration + p.soil.turbulent_fluxes.vapor_flux_liq
     @. p.soil.top_bc.heat =
         p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf
 end