Merge pull request #39 from jagoosw/jsw/add-biogeochemistry-interface

Adds biogeochemistry interface + fix example to work out of the box

Author: navidcy

Project.toml

@@ -2,7 +2,7 @@ name = "ClimaOcean"
 uuid = "0376089a-ecfe-4b0e-a64f-9c555d74d754"
 license = "MIT"
 authors = ["Climate Modeling Alliance and contributors"]
-version = "0.1.0"
+version = "0.1.2"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

examples/one_degree_near_global_simulation.jl

@@ -8,16 +8,20 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
 using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities:
     MixingLength, TurbulentKineticEnergyEquation, CATKEVerticalDiffusivity
 
-using ParameterEstimocean.Parameters: closure_with_parameters
+#using ParameterEstimocean.Parameters: closure_with_parameters
 using JLD2
 
+const years = 365days
+
+architecuture = GPU()
+
 #####
 ##### Boundary layer turbulence closure options
 #####
 
 # "Ri-based" --- uses calibrated defaults in Oceananigans
 ri_based = RiBasedVerticalDiffusivity() 
-
+#=
 # CATKE with calibrated parameters loaded from file
 neutral_default_mixing_length_parameters = (Cᵇu = Inf, Cᵇc = Inf, Cᵇe = Inf,
                                             Cˢu = Inf, Cˢc = Inf, Cˢe = Inf,
@@ -36,8 +40,9 @@ catke_parameters = load(catke_parameters_filename, "mean") # load ensemble mean
 
 catke = closure_with_parameters(neutral_catke, catke_parameters)
 
+=#
 # Choose closure
-boundary_layer_turbulence_closure = catke
+boundary_layer_turbulence_closure = ri_based#catke
 
 @show boundary_layer_turbulence_closure
 
@@ -49,7 +54,8 @@ start_time = 345days
 stop_time = start_time + 2years
 with_isopycnal_skew_symmetric_diffusivity = true
 
-simulation = one_degree_near_global_simulation(; start_time, stop_time,
+simulation = one_degree_near_global_simulation(architecure;
+    start_time, stop_time,
     with_isopycnal_skew_symmetric_diffusivity,
     boundary_layer_turbulence_closure,
     isopycnal_κ_skew = 900.0,
@@ -63,7 +69,7 @@ slices_save_interval = 2day
 fields_save_interval = 10days
 Nx, Ny, Nz = size(simulation.model.grid)
 
-dir = "/storage1/greg" 
+dir = "." 
 closure_name = typeof(boundary_layer_turbulence_closure).name.wrapper
 output_prefix = "near_global_$(Nx)_$(Ny)_$(Nz)_$closure_name"
 with_isopycnal_skew_symmetric_diffusivity || (output_prefix *= "_no_gm")

src/NearGlobalSimulations/NearGlobalSimulations.jl

@@ -144,6 +144,8 @@ function leith_viscosity(formulation; C_vort = 3.0, C_div = 3.0)
                                        parameters = (; Cζ, Cδ))
 end
 
+NoBiogeochemistry(args...; kwargs...) = nothing
+
 include("one_degree_global_simulation.jl")
 include("quarter_degree_global_simulation.jl")
 

src/NearGlobalSimulations/one_degree_global_simulation.jl

@@ -31,6 +31,9 @@ function one_degree_near_global_simulation(architecture = GPU();
     initial_conditions                           = datadep"near_global_one_degree/initial_conditions_month_01_360_150_48.jld2",
     bathymetry_path                              = datadep"near_global_one_degree/bathymetry_lat_lon_360_150.jld2",
     surface_boundary_conditions_path             = datadep"near_global_one_degree/surface_boundary_conditions_12_months_360_150.jld2",
+    biogeochemistry = NoBiogeochemistry,
+    biogeochemistry_kwargs = (),
+    progress_callback = false,
     )
 
     size == (360, 150, 48) || throw(ArgumentError("Only size = (360, 150, 48) is supported."))
@@ -104,7 +107,7 @@ function one_degree_near_global_simulation(architecture = GPU();
 
     vitd = VerticallyImplicitTimeDiscretization()
 
-    horizontal_κ = (T=0, S=0, e=horizontal_tke_diffusivity)
+    horizontal_κ = merge(NamedTuple{tracers}([0 for tracer in tracers]), (e = horizontal_tke_diffusivity, ))
     horizontal_diffusivity = HorizontalScalarDiffusivity(ν=horizontal_viscosity, κ=horizontal_κ)
     vertical_viscosity   = VerticalScalarDiffusivity(vitd, ν=νz, κ=background_vertical_diffusivity)
 
@@ -186,7 +189,7 @@ function one_degree_near_global_simulation(architecture = GPU();
 
     equation_of_state = TEOS10EquationOfState(; reference_density)
     buoyancy = SeawaterBuoyancy(; equation_of_state)
-    coriolis = HydrostaticSphericalCoriolis(scheme = ActiveCellEnstrophyConservingScheme())
+    coriolis = HydrostaticSphericalCoriolis()
     free_surface = ImplicitFreeSurface()
 
     @info "Building a model..."; start=time_ns()
@@ -195,7 +198,8 @@ function one_degree_near_global_simulation(architecture = GPU();
                                         momentum_advection = VectorInvariant(), 
                                         tracer_advection = WENO(underlying_grid),
                                         closure = closures,
-                                        boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs))
+                                        boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs),
+                                        biogeochemistry = biogeochemistry(; grid, biogeochemistry_kwargs...))
 
     @info "... built $model."
     @info "Model building time: " * prettytime(1e-9 * (time_ns() - start))
@@ -237,7 +241,7 @@ function one_degree_near_global_simulation(architecture = GPU();
         return nothing
     end
 
-    simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
+    progress_callback && (simulation.callbacks[:progress] = Callback(progress, IterationInterval(10)))
 
     return simulation
 end

src/NearGlobalSimulations/quarter_degree_global_simulation.jl

@@ -30,6 +30,8 @@ function quarter_degree_near_global_simulation(architecture = GPU();
         salt_surface_boundary_conditions_path        = datadep"near_global_quarter_degree/salt-1440x600-latitude-75.jld2",
         u_stress_surface_boundary_conditions_path    = datadep"near_global_quarter_degree/tau_x-1440x600-latitude-75.jld2",
         v_stress_surface_boundary_conditions_path    = datadep"near_global_quarter_degree/tau_y-1440x600-latitude-75.jld2",
+        biogeochemistry = NoBiogeochemistry,
+        biogeochemistry_kwargs = ()
 )
 
     bathymetry_file = jldopen(bathymetry_path)
@@ -164,7 +166,8 @@ function quarter_degree_near_global_simulation(architecture = GPU();
                                                                                vertical_scheme    = WENO()),
                                           closure = closures,
                                           boundary_conditions = (u=u_bcs, v=v_bcs, T=T_bcs, S=S_bcs),
-                                          tracer_advection = WENO(underlying_grid))
+                                          tracer_advection = WENO(underlying_grid),
+                                          biogeochemistry = biogeochemistry(; grid, biogeochemistry_kwargs...))
 
     @info "... built $model."
     @info "Model building time: " * prettytime(1e-9 * (time_ns() - start))