jb/add socrates test cases #1325

jbphyswx · 2023-02-10T16:14:43Z

Purpose

Add SOCRATES as a test case option. There are 6 flights/test cases [01, 09,10,11,12,13] forced either entirely with ERA5 data or partially with flight observations.

The forcing files are as described in https://atmos.uw.edu/~ratlas/SOCRATES-LES-cases.html following the description in https://doi.org/10.1029/2020MS002205.

The forcings are transformed from Atlas's formulation to our formulation in https://github.com/jbphyswx/SOCRATES_Single_Column_Forcings

The provided forcings are:

T Liquid Temperature $T_L$
divT Large-scale Horizontal Temperature advection
q water vapor mass-mixing ratio
divq Large-scale Horizontal water-vapor advection
lev pressure or p
omega Vertical Pressure Velocity
Ps Surface Pressure
PTend Surface pressure tendency
u zonal wind
ug geostrophic zonal wind
v meridional wind
vg geostrophic meridional wind

from which we can calculate, after transforming from pressure to altitude:

H_nudge from $T_L, q, p \rightarrow \theta_{li}$
dTdt_hadv from divT
qt_nudge from q
dqtdt_hadv from divq
subsidence from omega, p from the calculation in Appendix C of Atlas et al 2020
u_nudge from u
ug_nudge from ug
v_nudge from v
vg_nudge from vg

Other than that, Atlas et al (2020) assumes RRTMG which I haven't used here, and the surface flux scheme I have left as default Monin-Obukhov.

We create a SOCRATES subtype for each flight number/forcing type in driver/Cases.jl. These pull the forcings (; dTdt_hadv, H_nudge, dqtdt_hadv, qt_nudge, subsidence, u_nudge, v_nudge, ug_nudge, vg_nudge) from SOCRATES_Single_Column_Forcings.jl/src/process_case.jl : process_case(). process_case() interpolates the forcings onto your requested z grid, and then creates interpolation functions in time that TurbulenceConvection.jl can call at each timestep to obtain the correct profiles to relax towards.

Otherwise, the implementation in driver/Cases.jl, driver/Foricng.jl, driver/dycore.jl should be very similar to the other cases, most particularly LESData forcings, with the difference that we havent' a specific z grid predetermined; currently I'm running on the same vertical grids Atlas et al (2020) ran on.

To-do

Figure out how to store the test case data and allow it to be interpolated to an arbitrary vertical grid while still having functions to be evaluated in time, while not necessarily relying on another package or something else that ( should we pick one vertical grid and store those data in AtmosphericProfilesLibrary.jl or something?

Content

I have read and checked the items on the review checklist.

jbphyswx · 2023-02-10T16:16:38Z

This won't build because there's no way to specify a registry in the toml and I used my registry at https://github.com/jbphyswx/MyRegistry to register SOCRATES_Single_Column_Forcings

Ideas to get around that welcome... I dont want to add all the code in that package to this repo.

trontrytel · 2023-02-13T23:06:41Z

Hi @costachris and @charleskawczynski - could you take a look here?

@jbphyswx is adding test cases that are based on SOCRATES papers. The test cases need profiles not only for initial conditions but also for forcings during the simulations. I know that for LES driven SCMs we just store the relevant data online in a box. But alternatively, we also have the Atmospheric Profiles Library repository, where we store the initial conditions from different paper based cases.

Would you suggest following the LES driven SCM cases setup? Or could we try to extend the APL.jl to include the SOCRATES data? @jbphyswx so far put the data in a separate repo but we would probably like to avoid adding one more repo to the stack.

Project.toml

trontrytel · 2023-02-13T23:09:52Z

@jbphyswx - could you add in the PR description the list of forcings that are defined for the SOCRATES cases?

driver/Cases.jl

integration_tests/Manifest.toml

driver/generate_namelist.jl

trontrytel · 2023-02-13T23:34:59Z

driver/dycore.jl

@@ -335,18 +339,20 @@ function compute_gm_tendencies!(
    C123 = CCG.Covariant123Vector
    C12 = CCG.Contravariant12Vector
    lg = CC.Fields.local_geometry_field(axes(ρ_c))
-    @. tendencies_gm_uₕ -= f × (C12(C123(prog_gm_uₕ)) - C12(C123(aux_gm_uₕ_g)))
+    @. tendencies_gm_uₕ -= f × (C12(C123(prog_gm_uₕ)) - C12(C123(aux_gm_uₕ_g))) # geostrophic matters here...


Is this comment needed?

it was a refernence to how it's hard to tell where the geostrophic wind forcing matters relative to the actual wind nudging tendency and where and how each should be applied since they both operate on winds... it's not clear how they interplay in general so i made a note of that. I've removed the note now

driver/dycore.jl

charleskawczynski · 2023-02-13T23:50:32Z

It looks like the things being done in this branch are pretty complicated. Maybe a short zoom meeting would be best to discuss (@jbphyswx), if you'd like?

driver/dycore.jl

trontrytel · 2023-02-14T00:07:12Z

driver/common_spaces.jl

@@ -65,6 +65,19 @@ function construct_mesh(namelist; FT = Float64)
        end
        nz = isnothing(nz) ? Int(zmax ÷ Δz) : Int(nz)
        Δz = isnothing(Δz) ? FT(zmax ÷ nz) : FT(Δz)
+    elseif typeof(Cases.get_case(namelist)) <: Cases.SOCRATES # (we dont have a Cases.SOCRATES instance so improvise) # you can't easily create your own mesh in ClimaCore.Meshes so we'll use some representation of the one they have...


Would it be possible to keep the pattern here and only define Δz and nz here? The z_mesh is then created for all cases in lines 102-108

I think I did this because the grid is not uniform but an exact copy of Atlas's -- if you want to just switch to a different vertical grid to keep the same format, we could.

driver/dycore.jl

trontrytel · 2023-02-14T01:08:31Z

driver/Forcing.jl

@@ -33,3 +33,31 @@ function initialize(::ForcingBase{ForcingLES}, grid, state, LESDat::LESData)
        aux_gm.v_nudge[k] = nt.v_nudge[k]
    end
 end
+
+function initialize(forcing::ForcingBase{T}, grid, state, param_set) where {T<:ForcingSOCRATES}# added param_set so we can calculate stuff....
+    FT = eltype(param_set)


param_set is not used anywhere. You can get the FT type from aux_gm

driver/Forcing.jl

driver/Cases.jl

trontrytel · 2023-02-14T01:32:01Z

driver/Cases.jl

@@ -731,7 +783,7 @@ function initialize_forcing(::ARM_SGP, forcing, grid::Grid, state, param_set)
    return nothing
 end

-function update_forcing(::ARM_SGP, grid, state, t::Real, param_set)
+function update_forcing(::ARM_SGP, grid, state, t::Real, param_set) #- should these be in Forcing.jl?


I think it's ok to keep it in Cases. We want to apply minimal number of changes to the code that would get us to SOCRATES setup

trontrytel · 2023-02-14T01:40:24Z

Hi @jbphyswx thank you for opening the PR. I left some initial comments and will keep working on this tomorrow.

The main question I have right now is how to simplify the case specification, so that we don't have to create the sub-types for different research flights and forcing types? How much of that could be hidden by just providing the path to the dataset that specifies the forcing? Similar as it is done with LES driven SCM cases, where the same case is used to run different simulations by just swapping the path to data file?

trontrytel · 2023-02-24T22:28:30Z

The Socrates dependency got merged #1328 🎉

Could you merge that into your PR and see how things look @jbphyswx ?

trontrytel · 2023-02-28T20:59:27Z

Would it be possible for this PR to not change anything in the toml files? Are there any more dependencies that you need updated @jbphyswx?

costachris · 2023-02-28T23:37:50Z

driver/Cases.jl

-get_case(::Val{:GABLS}) = GABLS()
-get_case(::Val{:DryBubble}) = DryBubble()
-get_case(::Val{:LES_driven_SCM}) = LES_driven_SCM()
+get_case(namelist::Dict, param_set::APS)  = get_case(namelist["meta"]["casename"], namelist, param_set) 


I don't think you're using param_set anymore (right?), so these function definitions for all the other cases can be reverted (and a new one added for SOCRATES).

the one i have does use param_set for socrates -- the constructor uses it to fill out the fields around line 108. The alternative is to keep that call the same and use like a SOCRATESDat file -- but i think then we'd have to pass some more information to surface_ref_state, and add a dispatch for that separately for socrates. There are some other small complications with each so it would be up to yall to say which is better.

costachris · 2023-02-28T23:44:36Z

driver/Cases.jl

+        aux_gm.q_tot[k]     = IC.qt_nudge[k]
+        prog_gm_u[k]        = IC.u_nudge[k]
+        prog_gm_v[k]        = IC.v_nudge[k]
+        aux_gm.tke[k]       = 0 # what is this supposed to be? unset for interactive? maybe we dont need to set it initially at all idk... didn't check yet but these get updated later interactively so hope it's fine


You may have to experiment with this. If results don't match what you expect try initializing with a TKE profile like LES_driven_SCM

driver/dycore.jl

costachris · 2023-03-30T00:59:14Z

This is somehow changing the behavior of non-socrates cases in the CI. See ie. life_cycle_tan2018 in the buildkite tests.

…tie-ins to namelist args later

…socrates type formulation

…ntained in SOCRATES object and don't generate a SOCRATESDat array using param_set -- now get_case has it

…initialize since we don't need it anyway since our forcing funcs already used param_set in the SOCRATES object

…e call in main.jl

…ing param_set from aux_data_kwarg

…ironment I guess)

…ow that SOCRATES is a smaller container

… comparison of mse results, and use a default empty mse directory until we get a good one for each case later if ever (and also we could add truth for each socrates subtype later)

…rentiated betweenthe cases yet, and 2, we don't have a truth yet for calculating MSE anyway

…nt_repl_scripts calling parsed_args_from_command_line_flags doesn't crash on the latters assert iseven parsed_args_list flag