Merge pull request #430 from SpeedyWeather/mk/parametrizations

Adjust time step with radius

src/dynamics/spectral_grid.jl

@@ -1,6 +1,9 @@
 const DEFAULT_NF = Float32
 const DEFAULT_MODEL = PrimitiveDry
 const DEFAULT_GRID = OctahedralGaussianGrid
+const DEFAULT_RADIUS = 6.371e6
+const DEFAULT_TRUNC = 31
+const DEFAULT_NLEV = 8
 
 """
 Defines the horizontal spectral resolution and corresponding grid and the
@@ -15,7 +18,7 @@ Base.@kwdef struct SpectralGrid
 
     # HORIZONTAL
     "horizontal resolution as the maximum degree of spherical harmonics"
-    trunc::Int = 31
+    trunc::Int = DEFAULT_TRUNC
 
     "horizontal grid used for calculations in grid-point space"
     Grid::Type{<:AbstractGrid} = DEFAULT_GRID
@@ -24,7 +27,7 @@ Base.@kwdef struct SpectralGrid
     dealiasing::Float64 = 2
 
     "radius of the sphere [m]"
-    radius::Float64 = 6.371e6
+    radius::Float64 = DEFAULT_RADIUS
 
     # SIZE OF GRID from trunc, Grid, dealiasing:
     "number of latitude rings on one hemisphere (Equator incl)"
@@ -35,7 +38,7 @@ Base.@kwdef struct SpectralGrid
 
     # VERTICAL
     "number of vertical levels"
-    nlev::Int = 8
+    nlev::Int = DEFAULT_NLEV
 
     "coordinates used to discretize the vertical"
     vertical_coordinates::VerticalCoordinates = SigmaCoordinates(;nlev)
@@ -67,7 +70,7 @@ function Base.show(io::IO,SG::SpectralGrid)
     lat = get_lat(Grid,nlat_half)
     res_eq_y = (lat[nlat_half] - lat[nlat_half+1])*radius/1000
 
-    s(x) = @sprintf("%.3g",x)
+    s(x) = x > 1000 ? @sprintf("%i",x) : @sprintf("%.3g",x)
 
     println(io,"$(typeof(SG)):")
     println(io,"├ Spectral:   T$trunc LowerTriangularMatrix{Complex{$NF}}, radius = $radius m")

src/dynamics/time_integration.jl

@@ -13,10 +13,10 @@ Base.@kwdef mutable struct Leapfrog{NF} <: TimeStepper{NF}
     Δt_at_T31::Second = Minute(30)
 
     "radius of sphere [m], used for scaling"
-    radius::NF = 6.371e6
+    radius::NF = DEFAULT_RADIUS
 
     "adjust Δt_at_T31 with the output_dt to reach output_dt exactly in integer time steps"
-    adjust_with_output::Bool = true 
+    adjust_with_output::Bool = true
 
     # NUMERICS
     "Robert (1966) time filter coefficeint to suppress comput. mode"
@@ -27,7 +27,7 @@ Base.@kwdef mutable struct Leapfrog{NF} <: TimeStepper{NF}
 
     # DERIVED FROM OPTIONS    
     "time step Δt [ms] at specified resolution"
-    Δt_millisec::Millisecond = get_Δt_millisec(Second(Δt_at_T31), trunc, adjust_with_output)
+    Δt_millisec::Millisecond = get_Δt_millisec(Second(Δt_at_T31), trunc, radius, adjust_with_output)
 
     "time step Δt [s] at specified resolution"
     Δt_sec::NF = Δt_millisec.value/1000
@@ -46,12 +46,18 @@ adjusted to the closest divisor of `output_dt` so that the output time axis is k
 function get_Δt_millisec(
     Δt_at_T31::Dates.TimePeriod,
     trunc,
+    radius,
     adjust_with_output::Bool,
     output_dt::Dates.TimePeriod = DEFAULT_OUTPUT_DT,
 )
     # linearly scale Δt with trunc+1 (which are often powers of two)
-    resolution_factor = (32/(trunc+1))
-    Δt_at_trunc = Second(Δt_at_T31).value*resolution_factor
+    resolution_factor = (DEFAULT_TRUNC+1)/(trunc+1)
+
+    # radius also affects grid spacing, scale proportionally
+    radius_factor = radius/DEFAULT_RADIUS
+
+    # maybe rename to _at_trunc_and_radius?
+    Δt_at_trunc = Second(Δt_at_T31).value * resolution_factor * radius_factor
 
     if adjust_with_output && (output_dt > Millisecond(0))
         k = round(Int,Second(output_dt).value / Δt_at_trunc)
@@ -103,7 +109,7 @@ function initialize!(L::Leapfrog,model::ModelSetup)
 
     if L.adjust_with_output
         # take actual output dt from model.output and recalculate timestep
-        L.Δt_millisec = get_Δt_millisec(L.Δt_at_T31, L.trunc, L.adjust_with_output, output_dt)
+        L.Δt_millisec = get_Δt_millisec(L.Δt_at_T31, L.trunc, L.radius, L.adjust_with_output, output_dt)
         L.Δt_sec = L.Δt_millisec.value/1000
         L.Δt = L.Δt_sec/L.radius
     end