|
| 1 | +# # Cloudless Aquaplanet |
| 2 | + |
| 3 | +redirect_stderr(IOContext(stderr, :stacktrace_types_limited => Ref(false))) |
| 4 | + |
| 5 | +#= |
| 6 | +## Configuration initialization |
| 7 | +=# |
| 8 | + |
| 9 | +#= |
| 10 | +### Package Import |
| 11 | +=# |
| 12 | + |
| 13 | +## standard packages |
| 14 | +import Dates |
| 15 | +import YAML |
| 16 | + |
| 17 | +# ## ClimaESM packages |
| 18 | +import ClimaAtmos as CA |
| 19 | +import ClimaComms |
| 20 | +import ClimaCore as CC |
| 21 | + |
| 22 | +# ## Coupler specific imports |
| 23 | +import ClimaCoupler |
| 24 | +import ClimaCoupler: |
| 25 | + BCReader, |
| 26 | + ConservationChecker, |
| 27 | + Checkpointer, |
| 28 | + Diagnostics, |
| 29 | + FieldExchanger, |
| 30 | + FluxCalculator, |
| 31 | + Interfacer, |
| 32 | + Regridder, |
| 33 | + TimeManager, |
| 34 | + Utilities |
| 35 | + |
| 36 | +pkg_dir = pkgdir(ClimaCoupler) |
| 37 | + |
| 38 | +#= |
| 39 | +### Helper Functions |
| 40 | +=# |
| 41 | + |
| 42 | +## helpers for component models |
| 43 | +include("components/atmosphere/climaatmos.jl") |
| 44 | +include("components/ocean/slab_ocean.jl") |
| 45 | + |
| 46 | +## helpers for user-specified IO |
| 47 | +include("user_io/user_diagnostics.jl") |
| 48 | +include("user_io/user_logging.jl") |
| 49 | + |
| 50 | +include("user_io/io_helpers.jl") |
| 51 | + |
| 52 | +#= |
| 53 | +### Setup simulation parameters |
| 54 | +Here we follow Frierson et al. (2006) with diffusion for vertical transport. |
| 55 | +=# |
| 56 | + |
| 57 | +## run names |
| 58 | +run_name = "cloudless_aquaplanet" |
| 59 | +coupler_output_dir = "$run_name" |
| 60 | +const FT = Float64 |
| 61 | +restart_dir = "unspecified" |
| 62 | +restart_t = Int(0) |
| 63 | + |
| 64 | +## coupler simulation specific configuration |
| 65 | +Δt_cpl = Float64(400) |
| 66 | +t_end = "1000days" |
| 67 | +tspan = (Float64(0.0), Float64(time_to_seconds(t_end))) |
| 68 | +start_date = "19790301" |
| 69 | +hourly_checkpoint = true |
| 70 | +dt_rad = "6hours" |
| 71 | + |
| 72 | +## namelist |
| 73 | +config_dict = Dict( |
| 74 | + # general |
| 75 | + "FLOAT_TYPE" => string(FT), |
| 76 | + # file paths |
| 77 | + "atmos_config_file" => nothing, |
| 78 | + "coupler_toml_file" => nothing, |
| 79 | + "coupler_output_dir" => coupler_output_dir, |
| 80 | + "mode_name" => "", |
| 81 | + "run_name" => run_name, |
| 82 | + "atmos_config_repo" => "ClimaAtmos", |
| 83 | + # timestepping |
| 84 | + "dt" => "$(Δt_cpl)secs", |
| 85 | + "dt_save_to_sol" => "1days", |
| 86 | + "t_end" => t_end, |
| 87 | + "start_date" => "19790301", |
| 88 | + # domain |
| 89 | + "h_elem" => 4, |
| 90 | + "z_elem" => 10, |
| 91 | + "z_max" => 30000.0, # semi-high top |
| 92 | + "dz_bottom" => 300.0, |
| 93 | + "nh_poly" => 4, |
| 94 | + # output |
| 95 | + "dt_save_to_sol" => "1days", |
| 96 | + # numerics |
| 97 | + "apply_limiter" => false, |
| 98 | + "viscous_sponge" => false, |
| 99 | + "rayleigh_sponge" => false, |
| 100 | + "vert_diff" => "FriersonDiffusion", |
| 101 | + "hyperdiff" => "ClimaHyperdiffusion", |
| 102 | + # run |
| 103 | + "job_id" => run_name, |
| 104 | + "surface_setup" => "PrescribedSurface", |
| 105 | + # diagnostic (nested with period and short_name) |
| 106 | + "output_default_diagnostics" => false, |
| 107 | + "diagnostics" => [ |
| 108 | + Dict( |
| 109 | + "short_name" => |
| 110 | + ["mse", "lr", "mass_strf", "stab", "vt", "egr", "ua", "va", "wa", "ta", "rhoa", "pfull"], |
| 111 | + "period" => "6hours", |
| 112 | + "reduction" => "inst", |
| 113 | + ), |
| 114 | + ], |
| 115 | + # cloudless aquaplanet-specific |
| 116 | + "precip_model" => "0M", |
| 117 | + "moist" => "equil", |
| 118 | + "prognostic_surface" => "PrescribedSurfaceTemperature", |
| 119 | + "turb_flux_partition" => "CombinedStateFluxesMOST", |
| 120 | + "rad" => "gray", |
| 121 | + "dt_rad" => dt_rad, |
| 122 | + "albedo_model" => "CouplerAlbedo", |
| 123 | +) |
| 124 | + |
| 125 | +## merge dictionaries of command line arguments, coupler dictionary and component model dictionaries |
| 126 | +atmos_config_dict, config_dict = get_atmos_config_dict(config_dict) |
| 127 | +atmos_config_object = CA.AtmosConfig(atmos_config_dict) |
| 128 | + |
| 129 | +#= |
| 130 | +## Setup Communication Context |
| 131 | +=# |
| 132 | + |
| 133 | +comms_ctx = Utilities.get_comms_context(Dict("device" => "auto")) |
| 134 | +ClimaComms.init(comms_ctx) |
| 135 | + |
| 136 | +#= |
| 137 | +### I/O Directory Setup |
| 138 | +=# |
| 139 | + |
| 140 | +dir_paths = setup_output_dirs(output_dir = coupler_output_dir, comms_ctx = comms_ctx) |
| 141 | +ClimaComms.iamroot(comms_ctx) ? @info(config_dict) : nothing |
| 142 | + |
| 143 | +#= |
| 144 | +## Component Model Initialization |
| 145 | +=# |
| 146 | + |
| 147 | +#= |
| 148 | +### Atmosphere |
| 149 | +This uses the `ClimaAtmos.jl` model, with parameterization options specified in the `config_dict_atmos` dictionary. |
| 150 | +=# |
| 151 | + |
| 152 | +## init atmos model component |
| 153 | +atmos_sim = atmos_init(atmos_config_object); |
| 154 | +thermo_params = get_thermo_params(atmos_sim) |
| 155 | + |
| 156 | +#= |
| 157 | +### Boundary Space |
| 158 | +=# |
| 159 | + |
| 160 | +## init a 2D boundary space at the surface |
| 161 | +boundary_space = CC.Spaces.horizontal_space(atmos_sim.domain.face_space) # TODO: specify this in the coupler and pass it to all component models #665 |
| 162 | + |
| 163 | +#= |
| 164 | +### Surface Model: Slab Ocean |
| 165 | +=# |
| 166 | + |
| 167 | +ocean_sim = ocean_init( |
| 168 | + FT; |
| 169 | + tspan = tspan, |
| 170 | + dt = Δt_cpl, |
| 171 | + space = boundary_space, |
| 172 | + saveat = Float64(time_to_seconds(config_dict["dt_save_to_sol"])), |
| 173 | + area_fraction = ones(boundary_space), |
| 174 | + thermo_params = thermo_params, |
| 175 | + evolving = true, |
| 176 | +) |
| 177 | + |
| 178 | +#= |
| 179 | +## Coupler Initialization |
| 180 | +=# |
| 181 | + |
| 182 | +## coupler exchange fields |
| 183 | +coupler_field_names = ( |
| 184 | + :T_S, |
| 185 | + :z0m_S, |
| 186 | + :z0b_S, |
| 187 | + :ρ_sfc, |
| 188 | + :q_sfc, |
| 189 | + :surface_direct_albedo, |
| 190 | + :surface_diffuse_albedo, |
| 191 | + :beta, |
| 192 | + :F_turb_energy, |
| 193 | + :F_turb_moisture, |
| 194 | + :F_turb_ρτxz, |
| 195 | + :F_turb_ρτyz, |
| 196 | + :F_radiative, |
| 197 | + :P_liq, |
| 198 | + :P_snow, |
| 199 | + :radiative_energy_flux_toa, |
| 200 | + :P_net, |
| 201 | + :temp1, |
| 202 | + :temp2, |
| 203 | +) |
| 204 | +coupler_fields = |
| 205 | + NamedTuple{coupler_field_names}(ntuple(i -> CC.Fields.zeros(boundary_space), length(coupler_field_names))) |
| 206 | +Utilities.show_memory_usage(comms_ctx) |
| 207 | + |
| 208 | +## model simulations |
| 209 | +model_sims = (atmos_sim = atmos_sim, ocean_sim = ocean_sim); |
| 210 | + |
| 211 | +## dates |
| 212 | +date0 = date = Dates.DateTime(start_date, Dates.dateformat"yyyymmdd") |
| 213 | +dates = (; date = [date], date0 = [date0], date1 = [Dates.firstdayofmonth(date0)], new_month = [false]) |
| 214 | + |
| 215 | +#= |
| 216 | +## Initialize Callbacks |
| 217 | +=# |
| 218 | + |
| 219 | +checkpoint_cb = TimeManager.HourlyCallback( |
| 220 | + dt = FT(480), |
| 221 | + func = checkpoint_sims, |
| 222 | + ref_date = [dates.date[1]], |
| 223 | + active = hourly_checkpoint, |
| 224 | +) # 20 days |
| 225 | +update_firstdayofmonth!_cb = TimeManager.MonthlyCallback( |
| 226 | + dt = FT(1), |
| 227 | + func = TimeManager.update_firstdayofmonth!, |
| 228 | + ref_date = [dates.date1[1]], |
| 229 | + active = true, |
| 230 | +) |
| 231 | +dt_water_albedo = parse(FT, filter(x -> !occursin(x, "hours"), dt_rad)) |
| 232 | +albedo_cb = TimeManager.HourlyCallback( |
| 233 | + dt = dt_water_albedo, |
| 234 | + func = FluxCalculator.water_albedo_from_atmosphere!, |
| 235 | + ref_date = [dates.date[1]], |
| 236 | + active = true, |
| 237 | +) |
| 238 | +callbacks = |
| 239 | + (; checkpoint = checkpoint_cb, update_firstdayofmonth! = update_firstdayofmonth!_cb, water_albedo = albedo_cb) |
| 240 | + |
| 241 | +#= |
| 242 | +## Initialize turbulent fluxes |
| 243 | +
|
| 244 | +Decide on the type of turbulent flux partition (see `FluxCalculator` documentation for more details). |
| 245 | +=# |
| 246 | + |
| 247 | +turbulent_fluxes = nothing |
| 248 | +if config_dict["turb_flux_partition"] == "CombinedStateFluxesMOST" |
| 249 | + turbulent_fluxes = FluxCalculator.CombinedStateFluxesMOST() |
| 250 | +else |
| 251 | + error("turb_flux_partition must be CombinedStateFluxesMOST") |
| 252 | +end |
| 253 | + |
| 254 | +#= |
| 255 | +## Initialize Coupled Simulation |
| 256 | +=# |
| 257 | + |
| 258 | +cs = Interfacer.CoupledSimulation{FT}( |
| 259 | + comms_ctx, |
| 260 | + dates, |
| 261 | + boundary_space, |
| 262 | + coupler_fields, |
| 263 | + config_dict, |
| 264 | + nothing, # conservation checks |
| 265 | + [tspan[1], tspan[2]], |
| 266 | + atmos_sim.integrator.t, |
| 267 | + Δt_cpl, |
| 268 | + (; land = zeros(boundary_space), ocean = ones(boundary_space), ice = zeros(boundary_space)), |
| 269 | + model_sims, |
| 270 | + (;), # mode_specifics |
| 271 | + (), # coupler diagnostics |
| 272 | + callbacks, |
| 273 | + dir_paths, |
| 274 | + turbulent_fluxes, |
| 275 | + thermo_params, |
| 276 | +); |
| 277 | + |
| 278 | +#= |
| 279 | +## Restart component model states if specified in the config_dict |
| 280 | +=# |
| 281 | + |
| 282 | +if restart_dir !== "unspecified" |
| 283 | + for sim in cs.model_sims |
| 284 | + if Checkpointer.get_model_prog_state(sim) !== nothing |
| 285 | + Checkpointer.restart_model_state!(sim, comms_ctx, restart_t; input_dir = restart_dir) |
| 286 | + end |
| 287 | + end |
| 288 | +end |
| 289 | + |
| 290 | +#= |
| 291 | +## Initialize Component Model Exchange |
| 292 | +=# |
| 293 | + |
| 294 | +# 1.surface density (`ρ_sfc`): calculated by the coupler by adiabatically extrapolating atmospheric thermal state to the surface. |
| 295 | +# For this, we need to import surface and atmospheric fields. The model sims are then updated with the new surface density. |
| 296 | +FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims, cs.turbulent_fluxes) |
| 297 | +FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, cs.turbulent_fluxes) |
| 298 | +FieldExchanger.update_model_sims!(cs.model_sims, cs.fields, cs.turbulent_fluxes) |
| 299 | + |
| 300 | +# 2.surface vapor specific humidity (`q_sfc`): step surface models with the new surface density to calculate their respective `q_sfc` internally |
| 301 | +Interfacer.step!(ocean_sim, Δt_cpl) |
| 302 | + |
| 303 | +# 3.turbulent fluxes |
| 304 | +## import the new surface properties into the coupler (note the atmos state was also imported in step 3.) |
| 305 | +FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims, cs.turbulent_fluxes) # i.e. T_sfc, albedo, z0, beta, q_sfc |
| 306 | +## calculate turbulent fluxes inside the atmos cache based on the combined surface state in each grid box |
| 307 | +FluxCalculator.combined_turbulent_fluxes!(cs.model_sims, cs.fields, cs.turbulent_fluxes) # this updates the atmos thermo state, sfc_ts |
| 308 | + |
| 309 | +# 4.reinitialize models + radiative flux: prognostic states and time are set to their initial conditions. |
| 310 | +FieldExchanger.reinit_model_sims!(cs.model_sims) |
| 311 | + |
| 312 | +# 5.update all fluxes: coupler re-imports updated atmos fluxes |
| 313 | +FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, cs.turbulent_fluxes) |
| 314 | +FieldExchanger.update_model_sims!(cs.model_sims, cs.fields, cs.turbulent_fluxes) |
| 315 | + |
| 316 | +#= |
| 317 | +## Coupling Loop |
| 318 | +=# |
| 319 | + |
| 320 | +function solve_coupler!(cs) |
| 321 | + (; Δt_cpl, tspan, comms_ctx) = cs |
| 322 | + |
| 323 | + ClimaComms.iamroot(comms_ctx) && @info("Starting coupling loop") |
| 324 | + ## step in time |
| 325 | + for t in ((tspan[begin] + Δt_cpl):Δt_cpl:tspan[end]) |
| 326 | + |
| 327 | + cs.dates.date[1] = TimeManager.current_date(cs, t) |
| 328 | + |
| 329 | + ## print date on the first of month |
| 330 | + if cs.dates.date[1] >= cs.dates.date1[1] |
| 331 | + ClimaComms.iamroot(comms_ctx) && @show(cs.dates.date[1]) |
| 332 | + end |
| 333 | + |
| 334 | + ClimaComms.barrier(comms_ctx) |
| 335 | + |
| 336 | + ## update water albedo from wind at dt_water_albedo (this will be extended to a radiation callback from the coupler) |
| 337 | + TimeManager.trigger_callback!(cs, cs.callbacks.water_albedo) |
| 338 | + |
| 339 | + ## run component models sequentially for one coupling timestep (Δt_cpl) |
| 340 | + FieldExchanger.update_model_sims!(cs.model_sims, cs.fields, cs.turbulent_fluxes) |
| 341 | + |
| 342 | + ## step sims |
| 343 | + FieldExchanger.step_model_sims!(cs.model_sims, t) |
| 344 | + |
| 345 | + ## exchange combined fields and (if specified) calculate fluxes using combined states |
| 346 | + FieldExchanger.import_combined_surface_fields!(cs.fields, cs.model_sims, cs.turbulent_fluxes) # i.e. T_sfc, surface_albedo, z0, beta |
| 347 | + FluxCalculator.combined_turbulent_fluxes!(cs.model_sims, cs.fields, cs.turbulent_fluxes) |
| 348 | + |
| 349 | + FieldExchanger.import_atmos_fields!(cs.fields, cs.model_sims, cs.boundary_space, cs.turbulent_fluxes) # radiative and/or turbulent |
| 350 | + |
| 351 | + ## callback to update the fist day of month if needed (for BCReader) |
| 352 | + TimeManager.trigger_callback!(cs, cs.callbacks.update_firstdayofmonth!) |
| 353 | + |
| 354 | + ## callback to checkpoint model state |
| 355 | + TimeManager.trigger_callback!(cs, cs.callbacks.checkpoint) |
| 356 | + |
| 357 | + end |
| 358 | + |
| 359 | + return nothing |
| 360 | +end |
| 361 | + |
| 362 | +solve_coupler!(cs) |
0 commit comments