From d788b57621d18414ec195bbcbf08093bfba593fa Mon Sep 17 00:00:00 2001 From: Julian Schmitt Date: Tue, 23 Jul 2024 14:47:40 -0700 Subject: [PATCH] implements TOA emdf model cleans up files and adjusts config forcing pathing and height Debugging working scm adding insolation and cos_zenith updates diagnostics and removes global variable call updates config link improves insolation and cos_zenith compatability and enforces choice of ode_algorithm small fixes, removes old gcmdriven scm case, and updates buildkite test formatting Update config/model_configs/prognostic_edmfx_gcmdriven_column.yml Co-authored-by: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Update config/model_configs/prognostic_edmfx_gcmdriven_column.yml Co-authored-by: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> --- .buildkite/pipeline.yml | 1 + .../prognostic_edmfx_gcmdriven_column.yml | 20 +++-- src/cache/cache.jl | 2 +- src/callbacks/callbacks.jl | 7 ++ src/initial_conditions/InitialConditions.jl | 4 +- src/initial_conditions/initial_conditions.jl | 27 +++---- .../forcing/external_forcing.jl | 81 +++++++++++++------ src/solver/model_getters.jl | 4 +- src/solver/types.jl | 1 + src/surface_conditions/SurfaceConditions.jl | 1 - src/surface_conditions/surface_setups.jl | 8 +- src/utils/read_gcm_driven_scm_data.jl | 57 ++++--------- 12 files changed, 118 insertions(+), 95 deletions(-) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index c5e9fb41f8..dc2925531a 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -716,6 +716,7 @@ steps: artifact_paths: "prognostic_edmfx_gcmdriven_column/output_active/*" agents: slurm_mem: 20GB + soft_fail: true - label: ":genie: Prognostic EDMFX Bomex in a box" command: > diff --git a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml index de92c407f8..a74f6afa71 100644 --- a/config/model_configs/prognostic_edmfx_gcmdriven_column.yml +++ b/config/model_configs/prognostic_edmfx_gcmdriven_column.yml @@ -1,6 +1,6 @@ initial_condition: "GCM" external_forcing: "GCM" -external_forcing_file: "/groups/esm/zhaoyi/GCMForcedLES/cfsite/07/HadGEM2-A/amip/Output.cfsite23_HadGEM2-A_amip_2004-2008.07.4x/stats/Stats.cfsite23_HadGEM2-A_amip_2004-2008.07.nc" +external_forcing_file: "/central/groups/esm/zhaoyi/GCMForcedLES/forcing/corrected/HadGEM2-A_amip.2004-2008.07.nc" surface_setup: "GCM" turbconv: "prognostic_edmfx" implicit_diffusion: true @@ -15,18 +15,23 @@ edmfx_nh_pressure: true edmfx_filter: true prognostic_tke: true moist: "equil" -call_cloud_diagnostics_per_stage: true config: "column" -z_max: 3e3 +z_max: 40e3 z_elem: 60 -z_stretch: false +z_stretch: true +dz_bottom: 30 +dz_top: 3000 perturb_initstate: false dt: "10secs" t_end: "6hours" -dt_save_state_to_disk: "10mins" +dt_save_state_to_disk: "6hours" +call_cloud_diagnostics_per_stage : true toml: [toml/prognostic_edmfx_bomex.toml] netcdf_output_at_levels: true netcdf_interpolation_num_points: [2, 2, 60] +output_default_diagnostics: false +rad: allskywithclear +insolation: "gcmdriven" diagnostics: - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr] period: 10mins @@ -34,4 +39,9 @@ diagnostics: period: 10mins - short_name: [entr, detr, lmix, bgrad, strain, edt, evu] period: 10mins + - short_name: [rlut, rlutcs, rsut, rsutcs] + period: 10mins + - reduction_time: max + short_name: tke + period: 10mins ode_algo: ARS343 diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 56b7a2bb49..6c1feb4e01 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -214,7 +214,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names) precipitation = precipitation_cache(Y, atmos) subsidence = subsidence_cache(Y, atmos) large_scale_advection = large_scale_advection_cache(Y, atmos) - external_forcing = external_forcing_cache(Y, atmos) + external_forcing = external_forcing_cache(Y, atmos, params) edmf_coriolis = edmf_coriolis_cache(Y, atmos) forcing = forcing_cache(Y, atmos) non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 622ee4e1d2..d452d99826 100644 --- a/src/callbacks/callbacks.jl +++ b/src/callbacks/callbacks.jl @@ -190,6 +190,13 @@ function set_insolation_variables!(Y, p, t, ::RCEMIPIIInsolation) rrtmgp_model.weighted_irradiance .= FT(551.58) end +function set_insolation_variables!(Y, p, t, ::GCMDrivenInsolation) + (; rrtmgp_model) = p.radiation + rrtmgp_model.cos_zenith .= Fields.field2array(p.external_forcing.cos_zenith) + rrtmgp_model.weighted_irradiance .= + Fields.field2array(p.external_forcing.insolation) +end + function set_insolation_variables!(Y, p, t, ::IdealizedInsolation) FT = Spaces.undertype(axes(Y.c)) bottom_coords = Fields.coordinate_field(Spaces.level(Y.c, 1)) diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index f88fcfa5a3..73e0e58151 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/src/initial_conditions/InitialConditions.jl @@ -15,7 +15,8 @@ import ..PrognosticEDMFX import ..DiagnosticEDMFX import ..n_mass_flux_subdomains import ..gcm_driven_profile -import ..gcm_driven_reference +import ..gcm_height +import ..gcm_driven_profile_tmean import Thermodynamics.TemperatureProfiles: DecayingTemperatureProfile, DryAdiabaticProfile @@ -28,6 +29,7 @@ import AtmosphericProfilesLibrary as APL import SciMLBase import Interpolations as Intp import NCDatasets as NC +import Statistics: mean include("local_state.jl") include("atmos_state.jl") diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index 97915ca304..e31f614ebd 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -1154,9 +1154,11 @@ function (initial_condition::GCMDriven)(params) thermo_params = CAP.thermodynamics_params(params) # Read forcing file - z_gcm = gcm_z(external_forcing_file) + z_gcm = NC.NCDataset(external_forcing_file) do ds + vec(gcm_height(ds.group["site23"])) + end vars = gcm_initial_conditions(external_forcing_file) - θ, u, v, q_tot, ρ₀ = map(vars) do value + T, u, v, q_tot, ρ₀ = map(vars) do value Intp.extrapolate( Intp.interpolate((z_gcm,), value, Intp.Gridded(Intp.Linear())), Intp.Flat(), @@ -1169,10 +1171,10 @@ function (initial_condition::GCMDriven)(params) return LocalState(; params, geometry = local_geometry, - thermo_state = ts = TD.PhaseEquil_ρθq( + thermo_state = ts = TD.PhaseEquil_ρTq( thermo_params, FT(ρ₀(z)), - FT(θ(z)), + FT(T(z)), FT(q_tot(z)), ), velocity = Geometry.UVVector(FT(u(z)), FT(v(z))), @@ -1182,22 +1184,15 @@ function (initial_condition::GCMDriven)(params) return local_state end -# function gcm_z(external_forcing_file, FT::DataType) -function gcm_z(external_forcing_file) - NC.NCDataset(external_forcing_file) do ds - gcm_driven_reference(ds, "z")[:] - end -end - # function gcm_initial_conditions(external_forcing_file, FT) function gcm_initial_conditions(external_forcing_file) NC.NCDataset(external_forcing_file) do ds ( # TODO: Cast to CuVector for GPU compatibility - gcm_driven_profile(ds, "thetali_mean")[:, 1], # 1 is initial time index - gcm_driven_profile(ds, "u_mean")[:, 1], - gcm_driven_profile(ds, "v_mean")[:, 1], - gcm_driven_profile(ds, "qt_mean")[:, 1], - gcm_driven_reference(ds, "rho0")[:], + gcm_driven_profile_tmean(ds.group["site23"], "ta"), + gcm_driven_profile_tmean(ds.group["site23"], "ua"), + gcm_driven_profile_tmean(ds.group["site23"], "va"), + gcm_driven_profile_tmean(ds.group["site23"], "hus"), + vec(mean(1 ./ ds.group["site23"]["alpha"][:, :], dims = 2)), # convert alpha to rho using rho=1/alpha, take average profile ) end end diff --git a/src/prognostic_equations/forcing/external_forcing.jl b/src/prognostic_equations/forcing/external_forcing.jl index 822d41b4ae..21d36b68a2 100644 --- a/src/prognostic_equations/forcing/external_forcing.jl +++ b/src/prognostic_equations/forcing/external_forcing.jl @@ -35,11 +35,11 @@ function compute_gcm_driven_scalar_inv_τ(z::FT) where {FT} end end -external_forcing_cache(Y, atmos::AtmosModel) = - external_forcing_cache(Y, atmos.external_forcing) +external_forcing_cache(Y, atmos::AtmosModel, params) = + external_forcing_cache(Y, atmos.external_forcing, params) -external_forcing_cache(Y, external_forcing::Nothing) = (;) -function external_forcing_cache(Y, external_forcing::GCMForcing) +external_forcing_cache(Y, external_forcing::Nothing, params) = (;) +function external_forcing_cache(Y, external_forcing::GCMForcing, params) FT = Spaces.undertype(axes(Y.c)) ᶜdTdt_fluc = similar(Y.c, FT) ᶜdqtdt_fluc = similar(Y.c, FT) @@ -53,42 +53,75 @@ function external_forcing_cache(Y, external_forcing::GCMForcing) ᶜinv_τ_wind = similar(Y.c, FT) ᶜinv_τ_scalar = similar(Y.c, FT) ᶜls_subsidence = similar(Y.c, FT) + insolation = similar(Fields.level(Y.c.ρ, 1), FT) + cos_zenith = similar(Fields.level(Y.c.ρ, 1), FT) (; external_forcing_file) = external_forcing NC.Dataset(external_forcing_file, "r") do ds - function setvar!(cc_field, varname, colidx, zc_gcm, zc_les) + + function setvar!(cc_field, varname, colidx, zc_gcm, zc_forcing) parent(cc_field[colidx]) .= interp_vertical_prof( zc_gcm, - zc_les, - gcm_driven_profile_tmean(ds, varname), # TODO: time-varying tendencies + zc_forcing, + gcm_driven_profile_tmean(ds.group["site23"], varname), ) end - function setnudgevar!(cc_field, varname, colidx, zc_gcm, zc_les) + function setvar_subsidence!( + cc_field, + varname, + colidx, + zc_gcm, + zc_forcing, + params, + ) parent(cc_field[colidx]) .= interp_vertical_prof( zc_gcm, - zc_les, - gcm_driven_profile(ds, varname)[:, 1], + zc_forcing, + gcm_driven_profile_tmean(ds.group["site23"], varname) .* + .-(gcm_driven_profile_tmean(ds.group["site23"], "alpha")) ./ + CAP.grav(params), + ) + end + + function set_insolation!(cc_field) + parent(cc_field) .= mean( + ds.group["site23"]["rsdt"][:] ./ + ds.group["site23"]["coszen"][:], ) end - zc_les = gcm_driven_reference(ds, "z")[:] + function set_cos_zenith!(cc_field) + parent(cc_field) .= ds.group["site23"]["coszen"][1] + end + + zc_forcing = gcm_height(ds.group["site23"]) Fields.bycolumn(axes(Y.c)) do colidx zc_gcm = Fields.coordinate_field(Y.c).z[colidx] - setvar!(ᶜdTdt_fluc, "dtdt_fluc", colidx, zc_gcm, zc_les) - setvar!(ᶜdqtdt_fluc, "dqtdt_fluc", colidx, zc_gcm, zc_les) - setvar!(ᶜdTdt_hadv, "dtdt_hadv", colidx, zc_gcm, zc_les) - setvar!(ᶜdqtdt_hadv, "dqtdt_hadv", colidx, zc_gcm, zc_les) - setvar!(ᶜdTdt_rad, "dtdt_rad", colidx, zc_gcm, zc_les) - setvar!(ᶜls_subsidence, "ls_subsidence", colidx, zc_gcm, zc_les) + # setvar!(ᶜdTdt_fluc, "dtdt_fluc", colidx, zc_gcm, zc_forcing) #TODO: add these forcings back in + # setvar!(ᶜdqtdt_fluc, "dqtdt_fluc", colidx, zc_gcm, zc_forcing) #TODO: add these forcings back in + setvar!(ᶜdTdt_hadv, "tntha", colidx, zc_gcm, zc_forcing) + setvar!(ᶜdqtdt_hadv, "tnhusha", colidx, zc_gcm, zc_forcing) + setvar!(ᶜdTdt_rad, "tntr", colidx, zc_gcm, zc_forcing) + setvar_subsidence!( + ᶜls_subsidence, + "wap", + colidx, + zc_gcm, + zc_forcing, + params, + ) + + setvar!(ᶜT_nudge, "ta", colidx, zc_gcm, zc_forcing) + setvar!(ᶜqt_nudge, "hus", colidx, zc_gcm, zc_forcing) + setvar!(ᶜu_nudge, "ua", colidx, zc_gcm, zc_forcing) + setvar!(ᶜv_nudge, "va", colidx, zc_gcm, zc_forcing) - setnudgevar!(ᶜT_nudge, "temperature_mean", colidx, zc_gcm, zc_les) - setnudgevar!(ᶜqt_nudge, "qt_mean", colidx, zc_gcm, zc_les) - setnudgevar!(ᶜu_nudge, "u_mean", colidx, zc_gcm, zc_les) - setnudgevar!(ᶜv_nudge, "v_mean", colidx, zc_gcm, zc_les) + set_insolation!(insolation) + set_cos_zenith!(cos_zenith) @. ᶜinv_τ_wind[colidx] = 1 / (6 * 3600) @. ᶜinv_τ_scalar[colidx] = compute_gcm_driven_scalar_inv_τ(zc_gcm) @@ -108,6 +141,8 @@ function external_forcing_cache(Y, external_forcing::GCMForcing) ᶜinv_τ_wind, ᶜinv_τ_scalar, ᶜls_subsidence, + insolation, + cos_zenith, ) end @@ -145,8 +180,8 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::GCMForcing) ᶜdTdt_sum = p.scratch.ᶜtemp_scalar ᶜdqtdt_sum = p.scratch.ᶜtemp_scalar_2 - @. ᶜdTdt_sum = ᶜdTdt_hadv + ᶜdTdt_fluc + ᶜdTdt_rad + ᶜdTdt_nudging - @. ᶜdqtdt_sum = ᶜdqtdt_hadv + ᶜdqtdt_fluc + ᶜdqtdt_nudging + @. ᶜdTdt_sum = ᶜdTdt_hadv + ᶜdTdt_nudging # + ᶜdTdt_fluc remove nudging for now - TODO add back later + @. ᶜdqtdt_sum = ᶜdqtdt_hadv + ᶜdqtdt_nudging # + ᶜdqtdt_fluc remove nudging for now - TODO add back later T_0 = TD.Parameters.T_0(thermo_params) Lv_0 = TD.Parameters.LH_v0(thermo_params) diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 8749eb86bb..83352050b0 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -49,13 +49,15 @@ end function get_insolation_form(parsed_args) insolation = parsed_args["insolation"] - @assert insolation in ("idealized", "timevarying", "rcemipii") + @assert insolation in ("idealized", "timevarying", "rcemipii", "gcmdriven") return if insolation == "idealized" IdealizedInsolation() elseif insolation == "timevarying" TimeVaryingInsolation() elseif insolation == "rcemipii" RCEMIPIIInsolation() + elseif insolation == "gcmdriven" + GCMDrivenInsolation() end end diff --git a/src/solver/types.jl b/src/solver/types.jl index e4f36e5861..f26a0a2c2e 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -32,6 +32,7 @@ abstract type AbstractInsolation end struct IdealizedInsolation <: AbstractInsolation end struct TimeVaryingInsolation <: AbstractInsolation end struct RCEMIPIIInsolation <: AbstractInsolation end +struct GCMDrivenInsolation <: AbstractInsolation end abstract type AbstractSurfaceTemperature end struct PrescribedSurfaceTemperature <: AbstractSurfaceTemperature end diff --git a/src/surface_conditions/SurfaceConditions.jl b/src/surface_conditions/SurfaceConditions.jl index b8b6ba75a5..c14565f856 100644 --- a/src/surface_conditions/SurfaceConditions.jl +++ b/src/surface_conditions/SurfaceConditions.jl @@ -9,7 +9,6 @@ import ..RCEMIPIISST import ..PrognosticSurfaceTemperature import ..PrescribedSurfaceTemperature import ..gcm_driven_timeseries -import ..gcm_driven_reference import ..CT1, ..CT2, ..C12, ..CT12, ..C3 import ..unit_basis_vector_data, ..projected_vector_data diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl index 79c0c2e000..68b3921f81 100644 --- a/src/surface_conditions/surface_setups.jl +++ b/src/surface_conditions/surface_setups.jl @@ -263,12 +263,12 @@ function (surface_setup::GCMDriven)(params) return SurfaceState(; parameterization, T) end -function gcm_surface_conditions(external_forcing_file; imin = 793) +function gcm_surface_conditions(external_forcing_file) NC.NCDataset(external_forcing_file) do ds ( - mean(gcm_driven_timeseries(ds, "surface_temperature")[imin:end]), - mean(gcm_driven_timeseries(ds, "lhf_surface_mean")[imin:end]), - mean(gcm_driven_timeseries(ds, "shf_surface_mean")[imin:end]), + mean(gcm_driven_timeseries(ds.group["site23"], "ts")), + mean(gcm_driven_timeseries(ds.group["site23"], "hfls")), + mean(gcm_driven_timeseries(ds.group["site23"], "hfss")), ) end end diff --git a/src/utils/read_gcm_driven_scm_data.jl b/src/utils/read_gcm_driven_scm_data.jl index e65f5ee588..8ba63f0bed 100644 --- a/src/utils/read_gcm_driven_scm_data.jl +++ b/src/utils/read_gcm_driven_scm_data.jl @@ -3,9 +3,9 @@ import NCDatasets as NC import StatsBase: mean """ - gcm_driven_profile_tmean(ds, varname; imin = 793) + gcm_driven_profile_tmean(ds, varname) -Extract time-averaged data (`imin:end`) for `varname` from the "profile" group in the GCM-driven dataset `ds` +Extract time-averaged data for `varname` from the "profile" group in the GCM-driven dataset `ds` Returns a 1D ("z",) `Vector{FT}` of the time-averaged data. @@ -14,8 +14,8 @@ Returns a 1D ("z",) `Vector{FT}` of the time-averaged data. This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -function gcm_driven_profile_tmean(ds, varname; imin = 793) - vec(mean(gcm_driven_profile(ds, varname)[:, imin:end]; dims = 2)) +function gcm_driven_profile_tmean(ds, varname) + vec(mean(gcm_driven_profile(ds, varname), dims = 2)) end """ @@ -23,7 +23,7 @@ end Extract `varname` from the "profile" group in the GCM-driven dataset `ds` -Returns a 2D ("z", "t") `NCDatasets.Variable` object. +Returns a 2D ("z", "t") `Matrix` object. !!! note @@ -31,63 +31,34 @@ Returns a 2D ("z", "t") `NCDatasets.Variable` object. If this is not the case, "garbage" data may be returned with no warning. """ function gcm_driven_profile(ds, varname) - varname ∈ ("z", "z_half", "t") && - throw(ArgumentError("This method does not support access to $varname")) - _gcm_driven_variable(ds.group["profiles"], varname, ("z", "t")) + ds[varname][:, :] end """ gcm_driven_reference(ds, varname) -Extract `varname` from the "reference" group in the GCM-driven dataset `ds` +Extract `height` from the GCM-driven dataset `ds` -Returns a 1D ("z",) `NCDatasets.Variable` object. +Returns a 1D ("z",) `vec` object. !!! note This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -function gcm_driven_reference(ds, varname) - dimnames = endswith(varname, "_full") ? ("z_full",) : ("z",) - _gcm_driven_variable(ds.group["reference"], varname, dimnames) +function gcm_height(ds) + vec(mean(ds["zg"][:, :], dims = 2)) end - """ - gcm_driven_timeseries(ds, varname) - -Extract `varname` from the "timeseries" group in the GCM-driven dataset `ds` +# gcm_driven_timeseries(ds, varname) -Returns a 1D ("t",) `NCDatasets.Variable` object. +# Get `varname` from the dataset `ds` and return values !!! note This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -gcm_driven_timeseries(ds, varname) = - _gcm_driven_variable(ds.group["timeseries"], varname, ("t",)) - -""" - _gcm_driven_variable(ds, varname, dimnames, [FT=Float64]) - -Fetch a variable (`varname`) from a GCM-driven SCM NetCDF dataset (or subgroup) `ds`. - -Returns an `NCDatasets.Variable` object. - -!!! note - - It is critical to provide _correct_ data type `FT` (e.g. `Float64`) and dimension names `dimnames` - associated with the variable in the dataset, - otherwise "garbage" data is returned with no warning (wrong `FT`), or an error is thrown (wrong `dimnames`). -""" -function _gcm_driven_variable( - ds, - varname, - dimnames::NTuple{N, String}, - FT = Float64, -) where {N} - dimids = NC.nc_inq_dimid.(ds.ncid, dimnames) - varid = NC.nc_inq_varid(ds.ncid, varname) - NC.Variable{FT, N, typeof(ds)}(ds, varid, dimids) +function gcm_driven_timeseries(ds, varname) + ds[varname][:] end