Merge pull request #3405 from CliMA/zs/prescribe_clouds

add option to prescribe clouds in radiation

.buildkite/pipeline.yml

@@ -885,8 +885,13 @@ steps:
       - label: "GPU: GPU dry baroclinic wave - 4 gpus"
         key: "target_gpu_implicit_baroclinic_wave_4process"
         command:
-            nsys profile --delay 100 --trace=nvtx,cuda,mpi --output=target_gpu_implicit_baroclinic_wave_4process/output_active/report-%q{PMI_RANK}
           - mkdir -p target_gpu_implicit_baroclinic_wave_4process
+          # - >
+          #   srun --cpu-bind=threads --cpus-per-task=4
+          #   nsys profile --delay 100 --trace=nvtx,cuda,mpi --output=target_gpu_implicit_baroclinic_wave_4process/output_active/report-%q{PMI_RANK}
+          #   julia --threads=3 --color=yes --project=examples examples/hybrid/driver.jl
+          #   --config_file ${GPU_CONFIG_PATH}/target_gpu_implicit_baroclinic_wave.yml
+          #   --job_id target_gpu_implicit_baroclinic_wave_4process
           - >
             srun --cpu-bind=threads --cpus-per-task=4
             julia --threads=3 --color=yes --project=examples examples/hybrid/driver.jl

Artifacts.toml

@@ -36,3 +36,10 @@ lazy = true
     [[earth_orography_60arcseconds.download]]
     sha256 = "eca66c0701d1c2b9e271742314915ffbf4a0fae92709df611c323f38e019966e"
     url = "https://caltech.box.com/shared/static/4asrxcgl6xsgenfcug9p0wkkyhtqilgk.gz"
+
+[era5_cloud]
+git-tree-sha1 = "10742e0a2e343d13bb04df379e300a83402d4106"
+
+    [[era5_cloud.download]]
+    sha256 = "bb51e2f2d315b487e05a8d38944d4ad937ee4a40c43b68541220c5d54425e24a"
+    url = "https://caltech.box.com/shared/static/b6ur4ap4vo04j09vdulem96z9fxqlgyn.gz"

NEWS.md

@@ -6,6 +6,12 @@ Main
 
 ### Features
 
+### New option for prescribing clouds in radiation
+
+When `prescribe_clouds_in_radiation` is set to true, clouds in radiation 
+is prescribed from a file (monthly cloud properties in 2010 from ERA5).
+PR [3405](https://github.com/CliMA/ClimaAtmos.jl/pull/3405)
+
 ### ETOPO2022 60arc-second topography dataset. 
 
 - Update artifacts to use 60arc-second ETOPO2022 ice-surface topography

config/default_configs/default_config.yml

@@ -99,6 +99,9 @@ dt_rad:
 idealized_clouds:
   help: "Use idealized clouds in radiation model [`false` (default), `true`]"
   value: false
+prescribe_clouds_in_radiation:
+  help: "Use prescribed clouds in radiation model. Clouds are read from ERA5 data and updated every time radiation is called. The year 2010 is used and continuously repeated. This mode only affect radiation and is only relevant for the RRTGMP mode. [`false` (default), `true`]"
+  value: false
 insolation:
   help: "Insolation used in radiation model [`idealized` (default), `timevarying`, `rcemipii`]"
   value: "idealized"

config/model_configs/gpu_aquaplanet_dyamond.yml

@@ -20,4 +20,6 @@ t_end: "8hours"
 prescribe_ozone: true
 aerosol_radiation: true
 prescribed_aerosols: ["CB1", "CB2", "DST01", "DST02", "DST03", "DST04", "OC1", "OC2", "SO4", "SOA", "SSLT01", "SSLT02", "SSLT03", "SSLT04"]
+prescribe_clouds_in_radiation: true
+radiation_reset_rng_seed: true
 toml: [toml/longrun_aquaplanet.toml]

config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml

@@ -12,6 +12,8 @@ moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
 aerosol_radiation: true
+prescribe_clouds_in_radiation: true
+radiation_reset_rng_seed: true
 insolation: "timevarying"
 non_orographic_gravity_wave: true
 orographic_gravity_wave: "gfdl_restart"

src/cache/cache.jl

@@ -182,7 +182,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
 
     radiation_args =
         atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode ?
-        (params, atmos.ozone, aerosol_names, atmos.insolation) : ()
+        (start_date, params, atmos.ozone, aerosol_names, atmos.insolation) : ()
 
     hyperdiff = hyperdiffusion_cache(Y, atmos)
     rayleigh_sponge = rayleigh_sponge_cache(Y, atmos)

src/callbacks/callbacks.jl

@@ -77,6 +77,12 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
             evaluate!(field, tv, t)
         end
     end
+    if :prescribed_clouds_field in propertynames(p.radiation)
+        for (key, tv) in pairs(p.radiation.prescribed_cloud_timevaryinginputs)
+            field = getproperty(p.radiation.prescribed_clouds_field, key)
+            evaluate!(field, tv, t)
+        end
+    end
 
     FT = Spaces.undertype(axes(Y.c))
     thermo_params = CAP.thermodynamics_params(params)
@@ -157,13 +163,25 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
             )
             # RRTMGP needs lwp and iwp in g/m^2
             kg_to_g_factor = 1000
+            cloud_liquid_water_content =
+                radiation_mode.cloud isa PrescribedCloudInRadiation ?
+                p.radiation.prescribed_clouds_field.clwc :
+                cloud_diagnostics_tuple.q_liq
+            cloud_ice_water_content =
+                radiation_mode.cloud isa PrescribedCloudInRadiation ?
+                p.radiation.prescribed_clouds_field.ciwc :
+                cloud_diagnostics_tuple.q_ice
+            cloud_fraction =
+                radiation_mode.cloud isa PrescribedCloudInRadiation ?
+                p.radiation.prescribed_clouds_field.cc :
+                cloud_diagnostics_tuple.cf
             @. ᶜlwp =
-                kg_to_g_factor * Y.c.ρ * cloud_diagnostics_tuple.q_liq * ᶜΔz /
-                max(cloud_diagnostics_tuple.cf, eps(FT))
+                kg_to_g_factor * Y.c.ρ * cloud_liquid_water_content * ᶜΔz /
+                max(cloud_fraction, eps(FT))
             @. ᶜiwp =
-                kg_to_g_factor * Y.c.ρ * cloud_diagnostics_tuple.q_ice * ᶜΔz /
-                max(cloud_diagnostics_tuple.cf, eps(FT))
-            @. ᶜfrac = cloud_diagnostics_tuple.cf
+                kg_to_g_factor * Y.c.ρ * cloud_ice_water_content * ᶜΔz /
+                max(cloud_fraction, eps(FT))
+            @. ᶜfrac = cloud_fraction
         end
     end
 

src/parameterized_tendencies/radiation/RRTMGPInterface.jl

@@ -1,5 +1,7 @@
 module RRTMGPInterface
 
+import ..AbstractCloudInRadiation
+
 using RRTMGP
 import RRTMGP.AtmosphericStates as AS
 using ClimaCore: DataLayouts, Spaces, Fields
@@ -21,20 +23,36 @@ struct ClearSkyRadiation <: AbstractRRTMGPMode
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool
 end
-struct AllSkyRadiation <: AbstractRRTMGPMode
+struct AllSkyRadiation{ACR <: AbstractCloudInRadiation} <: AbstractRRTMGPMode
     idealized_h2o::Bool
     idealized_clouds::Bool
+    cloud::ACR
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool
-    "Reset the RNG seed before calling RRTGMP to a known value (the timestep number). When modeling cloud optics, RRTGMP uses a random number generator. Resetting the seed every time RRTGMP is called to a deterministic value ensures that the simulation is fully reproducible and can be restarted in a reproducible way. Disable this option when running production runs."
+    """
+    Reset the RNG seed before calling RRTGMP to a known value (the timestep number). 
+    When modeling cloud optics, RRTGMP uses a random number generator. 
+    Resetting the seed every time RRTGMP is called to a deterministic value ensures that 
+    the simulation is fully reproducible and can be restarted in a reproducible way. 
+    Disable this option when running production runs.
+    """
     reset_rng_seed::Bool
 end
-struct AllSkyRadiationWithClearSkyDiagnostics <: AbstractRRTMGPMode
+struct AllSkyRadiationWithClearSkyDiagnostics{
+    ACR <: AbstractCloudInRadiation,
+} <: AbstractRRTMGPMode
     idealized_h2o::Bool
     idealized_clouds::Bool
+    cloud::ACR
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool
-    "Reset the RNG seed before calling RRTGMP to a known value (the timestep number). When modeling cloud optics, RRTGMP uses a random number generator. Resetting the seed every time RRTGMP is called to a deterministic value ensures that the simulation is fully reproducible and can be restarted in a reproducible way. Disable this option when running production runs."
+    """
+    Reset the RNG seed before calling RRTGMP to a known value (the timestep number). 
+    When modeling cloud optics, RRTGMP uses a random number generator. 
+    Resetting the seed every time RRTGMP is called to a deterministic value ensures that 
+    the simulation is fully reproducible and can be restarted in a reproducible way. 
+    Disable this option when running production runs.
+    """
     reset_rng_seed::Bool
 end
 

src/parameterized_tendencies/radiation/radiation.jl

@@ -10,11 +10,14 @@ import .Parameters as CAP
 import RRTMGP
 import .RRTMGPInterface as RRTMGPI
 
-import Dates: Year
+import Dates: Year, Date
 import ClimaUtilities.TimeVaryingInputs:
-    TimeVaryingInput, LinearPeriodFillingInterpolation
+    TimeVaryingInput,
+    PeriodicCalendar,
+    LinearPeriodFillingInterpolation,
+    LinearInterpolation
 
-import Interpolations
+import Interpolations as Intp
 using Statistics: mean
 
 radiation_model_cache(Y, atmos::AtmosModel, args...) =
@@ -86,6 +89,7 @@ end
 function radiation_model_cache(
     Y,
     radiation_mode::RRTMGPI.AbstractRRTMGPMode,
+    start_date,
     params,
     ozone,
     aerosol_names,
@@ -257,10 +261,49 @@ function radiation_model_cache(
             kwargs...,
         )
     end
+    cloud_cache = (;)
+    if (radiation_mode isa RRTMGPI.AllSkyRadiation) ||
+       (radiation_mode isa RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics)
+        cloud_cache = get_cloud_cache(radiation_mode.cloud, Y, start_date)
+    end
     return merge(
         (; rrtmgp_model, ᶠradiation_flux = similar(Y.f, Geometry.WVector{FT})),
         insolation_cache(insolation_mode, Y),
+        cloud_cache,
+    )
+end
+
+get_cloud_cache(_, _, _) = (;)
+function get_cloud_cache(::PrescribedCloudInRadiation, Y, start_date)
+    target_space = axes(Y.c)
+    prescribed_cloud_names = ("cc", "clwc", "ciwc")
+    prescribed_cloud_names_as_symbols = Symbol.(prescribed_cloud_names)
+    extrapolation_bc = (Intp.Periodic(), Intp.Flat(), Intp.Flat())
+    timevaryinginputs = [
+        TimeVaryingInput(
+            joinpath(
+                @clima_artifact("era5_cloud", ClimaComms.context(Y.c)),
+                "era5_cloud.nc",
+            ),
+            name,
+            target_space;
+            reference_date = start_date,
+            regridder_type = :InterpolationsRegridder,
+            regridder_kwargs = (; extrapolation_bc),
+            method = LinearInterpolation(PeriodicCalendar(Year(1), Date(2010))),
+        ) for name in prescribed_cloud_names
+    ]
+
+    prescribed_clouds_field = similar(
+        Y.c,
+        NamedTuple{
+            prescribed_cloud_names_as_symbols,
+            NTuple{length(prescribed_cloud_names_as_symbols), eltype(Y.c.ρ)},
+        },
     )
+    prescribed_cloud_timevaryinginputs =
+        (; zip(prescribed_cloud_names_as_symbols, timevaryinginputs)...)
+    return (; prescribed_clouds_field, prescribed_cloud_timevaryinginputs)
 end
 
 insolation_cache(_, _) = (;)

src/solver/model_getters.jl

@@ -221,6 +221,12 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
     @assert idealized_h2o in (true, false)
     idealized_clouds = parsed_args["idealized_clouds"]
     @assert idealized_clouds in (true, false)
+    cloud = get_cloud_in_radiation(parsed_args)
+    if idealized_clouds && (cloud isa PrescribedCloudInRadiation)
+        error(
+            "idealized_clouds and prescribe_clouds_in_radiation cannot be true at the same time",
+        )
+    end
     add_isothermal_boundary_layer = parsed_args["add_isothermal_boundary_layer"]
     @assert add_isothermal_boundary_layer in (true, false)
     aerosol_radiation = parsed_args["aerosol_radiation"]
@@ -242,6 +248,10 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
     if !(radiation_name in ("allsky", "allskywithclear")) && reset_rng_seed
         @warn "reset_rng_seed does not have any effect with $radiation_name radiation option"
     end
+    if !(radiation_name in ("allsky", "allskywithclear")) &&
+       (cloud isa PrescribedCloudInRadiation)
+        @warn "prescribe_clouds_in_radiation does not have any effect with $radiation_name radiation option"
+    end
     return if radiation_name == "gray"
         RRTMGPI.GrayRadiation(add_isothermal_boundary_layer)
     elseif radiation_name == "clearsky"
@@ -254,6 +264,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
         RRTMGPI.AllSkyRadiation(
             idealized_h2o,
             idealized_clouds,
+            cloud,
             add_isothermal_boundary_layer,
             aerosol_radiation,
             reset_rng_seed,
@@ -262,6 +273,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
         RRTMGPI.AllSkyRadiationWithClearSkyDiagnostics(
             idealized_h2o,
             idealized_clouds,
+            cloud,
             add_isothermal_boundary_layer,
             aerosol_radiation,
             reset_rng_seed,
@@ -308,6 +320,12 @@ function get_ozone(parsed_args)
     return parsed_args["prescribe_ozone"] ? PrescribedOzone() : IdealizedOzone()
 end
 
+function get_cloud_in_radiation(parsed_args)
+    isnothing(parsed_args["prescribe_clouds_in_radiation"]) && return nothing
+    return parsed_args["prescribe_clouds_in_radiation"] ?
+           PrescribedCloudInRadiation() : InteractiveCloudInRadiation()
+end
+
 function get_forcing_type(parsed_args)
     forcing = parsed_args["forcing"]
     @assert forcing in (nothing, "held_suarez")

src/solver/types.jl

@@ -65,6 +65,28 @@ Refer to ClimaArtifacts for more information on how to obtain the artifact.
 """
 struct PrescribedOzone <: AbstractOzone end
 
+"""
+    AbstractCloudInRadiation
+
+Describe how cloud properties should be set in radiation.
+
+This is only relevant for RRTGMP.
+"""
+abstract type AbstractCloudInRadiation end
+
+"""
+    InteractiveCloudInRadiation
+
+Use cloud properties computed in the model
+"""
+struct InteractiveCloudInRadiation <: AbstractCloudInRadiation end
+
+"""
+    PrescribedCloudInRadiation
+
+Use monthly-average cloud properties from ERA5.
+"""
+struct PrescribedCloudInRadiation <: AbstractCloudInRadiation end
 
 abstract type AbstractSurfaceTemperature end
 struct PrescribedSurfaceTemperature <: AbstractSurfaceTemperature end