remove CairoMakie, alter US-MOz params, allocation removal (#957)

.buildkite/Project.toml

@@ -29,7 +29,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Poppler_jll = "9c32591e-4766-534b-9725-b71a8799265b"
 ProfileCanvas = "efd6af41-a80b-495e-886c-e51b0c7d77a3"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
@@ -47,3 +46,4 @@ cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd"
 ClimaAnalysis = "0.5.12"
 ClimaTimeSteppers = "0.7"
 Statistics = "1"
+Flux = "0.15"

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

@@ -15,12 +15,12 @@ lat = FT(38.7441) # degree
 long = FT(-92.2000) # degree
 
 # Soil parameters
-soil_ν = FT(0.5) # m3/m3
-soil_K_sat = FT(4e-7) # m/s, matches Natan
-soil_S_s = FT(1e-3) # 1/m, guess
-soil_vg_n = FT(2.05) # unitless
-soil_vg_α = FT(0.04) # inverse meters
-θ_r = FT(0.067) # m3/m3, from Wang et al. 2021 https://doi.org/10.5194/gmd-14-6741-2021
+soil_ν = FT(0.55) # m3/m3
+soil_K_sat = FT(4e-7) # m/s
+soil_S_s = FT(1e-2) # 1/m, guess
+soil_vg_n = FT(2.0) # unitless
+soil_vg_α = FT(0.05) # inverse meters
+θ_r = FT(0.04) # m3/m3, 
 
 # Soil makeup
 ν_ss_quartz = FT(0.1)
@@ -44,23 +44,25 @@ G_Function = CLMGFunction(χl)
 ϵ_canopy = FT(0.97)
 
 # Energy Balance model
-ac_canopy = FT(7.5e3)
+ac_canopy = FT(5e2)
 
 # Conductance Model
 g1 = FT(141) # Wang et al: 141 sqrt(Pa) for Medlyn model; Natan used 300.
 Drel = FT(1.6)
 g0 = FT(1e-4)
 
 #Photosynthesis model
-Vcmax25 = FT(9e-5) # from Yujie's paper 4.5e-5 , Natan used 9e-5
+Vcmax25 = FT(4.5e-5) # from Yujie's paper 4.5e-5
 
 # Plant Hydraulics and general plant parameters
+pc = FT(-2.5e6)
+sc = FT(5e-6)
 SAI = FT(1.0) # m2/m2 or: estimated from Wang et al, FT(0.00242) ?
 f_root_to_shoot = FT(3.5)
-K_sat_plant = 5e-9 # m/s # seems much too small?
+K_sat_plant = 3e-9 # m/s 
 ψ63 = FT(-4 / 0.0098) # / MPa to m, Holtzman's original parameter value is -4 MPa
 Weibull_param = FT(4) # unitless, Holtzman's original c param value
-a = FT(0.05 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
+a = FT(0.1 * 0.0098) # Holtzman's original parameter for the bulk modulus of elasticity
 conductivity_model =
     PlantHydraulics.Weibull{FT}(K_sat_plant, ψ63, Weibull_param)
 retention_model = PlantHydraulics.LinearRetentionCurve{FT}(a)

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

@@ -6,7 +6,7 @@ variables for running the simulation."""
 
 # Column dimensions - separation of layers at the top and bottom of the column:
 dz_bottom = FT(1.5)
-dz_top = FT(0.025)
+dz_top = FT(0.1)
 
 # Stem and leaf compartments and their heights:
 n_stem = Int64(1)
@@ -20,4 +20,4 @@ h_leaf = FT(9.5) # m
 t0 = Float64(0)
 
 # Time step size:
-dt = Float64(450)
+dt = Float64(1800)

experiments/integrated/fluxnet/data_tools.jl

@@ -109,9 +109,9 @@ function filter_column(driver_data::Matrix, column_name::String, units::String)
     # Check that the data exists and read it in if so
     col_dat, status = column_status(driver_data, column_name)
     # Set missing data threshold above which column is treated as absent
-    missing_threshold = 0.5
+    missing_threshold = 10.0
     # Set poor quality threshold above which the column undergoes no replacement using the quality flag
-    quality_threshold = 0.5
+    quality_threshold = 10.0
     # if it does not exist, exit
     if status == absent
         @info "Warning: Data for $column_name is absent"

experiments/integrated/fluxnet/ozark_pft.jl

@@ -9,7 +9,7 @@ import SciMLBase
 import ClimaTimeSteppers as CTS
 using ClimaCore
 import ClimaParams as CP
-using Plots
+using CairoMakie
 using Statistics
 using Dates
 using Insolation
@@ -299,7 +299,7 @@ output_dir = ClimaUtilities.OutputPathGenerator.generate_output_path(outdir)
 
 d_writer = ClimaDiagnostics.Writers.DictWriter()
 
-ref_time = DateTime(2005) # random. not sure what it should be
+ref_time = DateTime(2010)
 
 diags = ClimaLand.default_diagnostics(
     land,
@@ -353,7 +353,6 @@ short_names_1D = [
     "ghf", # G
     "rn", # Rn
 ]
-
 short_names_2D = [
     "swc", # swc_sfc or swc_5 or swc_10
     "tsoil", # soil_T_sfc or soil_T_5 or soil_T_10
@@ -557,78 +556,88 @@ end
 # Water content in soil and snow
 # Soil water content
 # Current resolution has the first layer at 0.1 cm, the second at 5cm.
-plt1 = Plots.plot(size = (1500, 800))
-Plots.plot!(
-    plt1,
-    model_times ./ 3600 ./ 24,
-    swc,
-    label = "1.25cm",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    ylim = [0.05, 0.55],
-    xlabel = "Days",
-    ylabel = "SWC [m/m]",
-    color = "blue",
-    margin = 10Plots.mm,
+fig = Figure(size = (1500, 800), fontsize = 20)
+ax1 = Axis(fig[2, 1], xlabel = "Days", ylabel = "SWC [m/m]")
+limits!(
+    ax1,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    0.05,
+    0.6,
 )
-
-plot!(
-    plt1,
+lines!(ax1, model_times ./ 3600 ./ 24, swc, label = "1.25cm", color = "blue")
+lines!(
+    ax1,
     model_times ./ 3600 ./ 24,
     si,
     color = "cyan",
     label = "Ice, 1.25cm",
 )
 
 if drivers.SWC.status != absent
-    Plots.plot!(
-        plt1,
+    lines!(
+        ax1,
         data_times ./ 3600 ./ 24,
         drivers.SWC.values[data_id_post_spinup],
         label = "Data",
     )
 end
-plt2 = Plots.plot(
-    data_times ./ 3600 ./ 24,
-    (drivers.P.values .* (-1e3 * 24 * 3600))[data_id_post_spinup],
-    label = "",
+
+ax2 = Axis(
+    fig[1, 1],
+    xlabel = "",
     ylabel = "Precipitation [mm/day]",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    margin = 10Plots.mm,
-    ylim = [-200, 0],
-    size = (1500, 400),
+    xticklabelsvisible = false,
+)
+limits!(
+    ax2,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    -500,
+    0.0,
 )
 
-Plots.plot(plt2, plt1, layout = grid(2, 1, heights = [0.3, 0.7]))
-Plots.savefig(joinpath(savedir, "ground_water_content.png"))
+lines!(
+    ax2,
+    data_times ./ 3600 ./ 24,
+    (drivers.P.values .* (-1e3 * 24 * 3600) .* (1 .- snow_frac))[data_id_post_spinup],
+    label = "Rain (data)",
+)
+lines!(
+    ax2,
+    data_times ./ 3600 ./ 24,
+    (drivers.P.values .* (-1e3 * 24 * 3600) .* snow_frac)[data_id_post_spinup],
+    label = "Snow (data)",
+)
+CairoMakie.save(joinpath(savedir, "ground_water_content.png"), fig)
+
 
-plt1 = Plots.plot(size = (1500, 800))
-Plots.plot!(
-    plt1,
+
+fig2 = Figure(size = (1500, 800))
+ax12 = Axis(fig2[1, 1], xlabel = "Days", ylabel = "Temperature (K)")
+limits!(
+    ax12,
+    minimum(model_times ./ 3600 ./ 24),
+    maximum(model_times ./ 3600 ./ 24),
+    265,
+    315,
+)
+
+lines!(
+    ax12,
     model_times ./ 3600 ./ 24,
     soil_T,
     label = "1.25cm",
-    xlim = [
-        minimum(model_times ./ 3600 ./ 24),
-        maximum(model_times ./ 3600 ./ 24),
-    ],
-    xlabel = "Days",
-    ylabel = "Temperature (K)",
     color = "blue",
-    margin = 10Plots.mm,
 )
+
 if drivers.TS.status != absent
-    Plots.plot!(
-        plt1,
+    lines!(
+        ax12,
         data_times ./ 3600 ./ 24,
         drivers.TS.values[data_id_post_spinup],
         label = "Data",
     )
 end
 
-Plots.savefig(joinpath(savedir, "soil_temperature.png"))
+CairoMakie.save(joinpath(savedir, "soil_temperature.png"), fig2)

experiments/integrated/fluxnet/plot_utils.jl

@@ -1,7 +1,7 @@
 """Plotting utilities for the integrated fluxnet site experiments"""
 
 using Interpolations
-using Plots
+using CairoMakie
 using StatsBase
 
 S_PER_DAY = 86400 # Number of seconds in a day
@@ -61,17 +61,16 @@ function plot_daily_avg(
         compute_diurnal_avg(data, [0:data_dt:(num_days * S_PER_DAY);], num_days)
 
     # Plot the data diurnal cycle
-    plt = Plots.plot(size = (800, 400))
-    Plots.plot!(
-        plt,
-        0.5:0.5:24,
-        data_hh_avg,
-        label = label,
+    fig = CairoMakie.Figure(size = (800, 400))
+    ax = CairoMakie.Axis(
+        fig[1, 1],
         xlabel = "Hour of day",
         ylabel = "$var_name $(unit)",
         title = "$var_name",
     )
-    Plots.savefig(joinpath(savedir, "$(var_name)_avg.png"))
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = label)
+    axislegend(ax, position = :lt)
+    CairoMakie.save(joinpath(savedir, "$(var_name)_avg.png"), fig)
 end
 
 """This function will be used to plot the comparison of the diurnal average of a 
@@ -101,17 +100,17 @@ function plot_avg_comp(
     R² = StatsBase.cor(model_hh_avg, data_hh_avg)^2
 
     # Plot the model and data diurnal cycles
-    plt = Plots.plot(size = (800, 400))
-    Plots.plot!(
-        plt,
-        0.5:0.5:24,
-        model_hh_avg,
-        label = "Model",
+    fig = CairoMakie.Figure(size = (800, 400))
+    ax = CairoMakie.Axis(
+        fig[1, 1],
         xlabel = "Hour of day",
         ylabel = "$var_name $(units)",
         title = "$var_name: RMSD = $(round(RMSD, digits = 2)), R² = $(round(R²[1][1], digits = 2))",
-        legend = :topleft,
     )
-    Plots.plot!(plt, 0.5:0.5:24, data_hh_avg, label = "Data")
-    Plots.savefig(joinpath(savedir, "$(var_name)_avg.png"))
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), model_hh_avg[:], label = "Model")
+
+    CairoMakie.lines!(ax, Array(0.5:0.5:24), data_hh_avg[:], label = "Data")
+    axislegend(ax, position = :lt)
+
+    CairoMakie.save(joinpath(savedir, "$(var_name)_avg.png"), fig)
 end