Merge pull request #157 from milankl/scale

Scale sst, tracer relaxation & consumption, checkerboard initial conditions, ridges

Author: milankl

src/constants.jl

@@ -46,16 +46,17 @@ struct Constants{T<:AbstractFloat,Tprog<:AbstractFloat}
     γ::T                    # frequency of interface relaxation
     cSmag::T                # Smagorinsky constant
     νB::T                   # biharmonic diffusion coefficient
-    rSST::T                 # tracer restoring timescale
+    τSST::T                 # tracer restoring timescale
     jSST::T                 # tracer consumption timescale
-    SSTmin::T               # tracer minimum
     ωFη::Float64            # frequency [1/s] of seasonal surface forcing incl 2π
     ωFx::Float64            # frequency [1/s] of seasonal wind x incl 2π
     ωFy::Float64            # frequency [1/2] of seasonal wind y incl 2π
 
     # SCALING
     scale::T                # multiplicative constant for low-precision arithmetics
-    scale_inv::T            # and it's inverse
+    scale_inv::T            # and its inverse
+    scale_sst::T            # scale for sst
+    scale_sst_inv::T        # and its inverse
 end
 
 """Generator function for the mutable struct Constants."""
@@ -105,22 +106,24 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
     νB = T(-P.νB/30000)           # linear scaling based on 540m^s/s at Δ=30km
 
     # TRACER ADVECTION
-    rSST = T(G.dtadvint/(P.τSST*3600*24))    # tracer restoring [1]
+    τSST = T(G.dtadvint/(P.τSST*3600*24))    # tracer restoring [1]
     jSST = T(G.dtadvint/(P.jSST*3600*24))    # tracer consumption [1]
-    SSTmin = T(P.SSTmin)
 
     # TIME DEPENDENT FORCING
     ωFη = -2π*P.ωFη/24/365.25/3600
     ωFx = 2π*P.ωFx/24/365.25/3600
     ωFy = 2π*P.ωFy/24/365.25/3600
 
     # SCALING
-    scale = T(P.scale)
-    scale_inv = T(1/P.scale)
+    scale = convert(T,P.scale)
+    scale_inv = convert(T,1/P.scale)
+    scale_sst = convert(T,P.scale_sst)
+    scale_sst_inv = convert(T,1/P.scale_sst)
 
     return Constants{T,Tprog}(  RKaΔt,RKbΔt,Δt_Δs,Δt_Δ,Δt_Δ_half,
                                 SSPRK3c,one_minus_α,
-                                g,cD,rD,γ,cSmag,νB,rSST,
-                                jSST,SSTmin,ωFη,ωFx,ωFy,
-                                scale,scale_inv)
+                                g,cD,rD,γ,cSmag,νB,τSST,jSST,
+                                ωFη,ωFx,ωFy,
+                                scale,scale_inv,
+                                scale_sst,scale_sst_inv)
 end

src/default_parameters.jl

@@ -20,20 +20,22 @@
     R::Real=6.371e6                     # Earth's radius [m]
 
     # SCALE
-    scale::Real=1                       # multiplicative scale for the momentum equations u,v
+    scale::Real=2^6                     # multiplicative scale for the momentum equations u,v
+    scale_sst::Real=2^15                # multiplicative scale for sst
 
     # WIND FORCING OPTIONS
-    wind_forcing_x::String="channel"    # "channel", "double_gyre", "shear","constant" or "none"
+    wind_forcing_x::String="shear"      # "channel", "double_gyre", "shear","constant" or "none"
     wind_forcing_y::String="constant"   # "channel", "double_gyre", "shear","constant" or "none"
     Fx0::Real=0.12                      # wind stress strength [Pa] in x-direction
     Fy0::Real=0.0                       # wind stress strength [Pa] in y-direction
-    seasonal_wind_x::Bool=false         # Change the wind stress with a sine of frequency ωFx,ωFy
+    seasonal_wind_x::Bool=true          # Change the wind stress with a sine of frequency ωFx,ωFy
     seasonal_wind_y::Bool=false         # same for y-component
-    ωFx::Real=1.0                       # frequency [1/year] for x component
-    ωFy::Real=1.0                       # frequency [1/year] for y component
+    ωFx::Real=2                         # frequency [1/year] for x component
+    ωFy::Real=2                         # frequency [1/year] for y component
 
     # BOTTOM TOPOGRAPHY OPTIONS
     topography::String="ridges"         # "ridge", "seamount", "flat", "ridges", "bathtub"
+    topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9]
     topo_height::Real=100.               # height of seamount [m]
     topo_width::Real=300e3              # horizontal scale [m] of the seamount
 
@@ -83,23 +85,21 @@
 
     # TRACER ADVECTION
     tracer_advection::Bool=true         # yes?
-    tracer_relaxation::Bool=false       # yes?
+    tracer_relaxation::Bool=true        # yes?
     tracer_consumption::Bool=false      # yes?
-    tracer_pumping::Bool=false          # yes?
-    injection_region::String="west"     # "west", "south", "rect" or "flat"
-    sst_initial::String="south"         # "west", "south", "rect", "flat" or "restart"
+    sst_initial::String="waves"         # "west", "south", "linear", "waves","rect", "flat" or "restart"
     sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                         # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
     Uadv::Real=0.2                      # Velocity scale [m/s] for tracer advection
-    SSTmax::Real=1.                     # tracer (sea surface temperature) max for restoring
-    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for restoring
-    τSST::Real=500.                     # tracer restoring time scale [days]
-    jSST::Real=365.                     # tracer consumption [days]
-    SST_λ0::Real=222e3                  # [m] transition position of relaxation timescale
-    SST_λs::Real=111e3                  # [m] transition width of relaxation timescale
-    SST_γ0::Real=8.35                   # [days] injection time scale
+    SSTmax::Real=1.                     # tracer (sea surface temperature) max for initial conditions
+    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for initial conditions
+    τSST::Real=100                      # tracer restoring time scale [days]
+    jSST::Real=365                      # tracer consumption [days]
     SSTw::Real=5e5                      # width [m] of the tangent used for the IC and interface relaxation
     SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
+    SSTwaves_ny::Real=4                 # wave crests/troughs in y
+    SSTwaves_nx::Real=SSTwaves_ny*L_ratio  # wave crests/troughs in x
+    SSTwaves_p::Real=1/2                # power for rectangles (p<1)/smootheness(p>=1) of waves
 
     # OUTPUT OPTIONS
     output::Bool=false                  # netcdf output?
@@ -143,7 +143,6 @@
     @assert diffusion in ["Smagorinsky", "constant"] "Diffusion '$diffusion' unsupported."
     @assert νB > 0.0     "Diffusion scaling constant νB has to be >0, $νB given."
     @assert cSmag > 0.0  "Smagorinsky coefficient cSmag has to be >0, $cSmag given."
-    @assert injection_region in ["west","south"] "Injection region '$injection_region' unsupported."
     @assert Uadv > 0.0   "Advection velocity scale Uadv has to be >0, $Uadv given."
     @assert output_dt > 0   "Output time step has to be >0, $output_dt given."
     @assert initial_cond in ["rest", "ncfile"] "Initial conditions '$initial_cond' unsupported."
@@ -175,20 +174,22 @@ Creates a Parameter struct with following options and default values
     R::Real=6.371e6                     # Earth's radius [m]
 
     # SCALE
-    scale::Real=1                       # multiplicative scale for the momentum equations u,v
+    scale::Real=2^6                     # multiplicative scale for the momentum equations u,v
+    scale_sst::Real=2^15                # multiplicative scale for sst
 
     # WIND FORCING OPTIONS
-    wind_forcing_x::String="channel"    # "channel", "double_gyre", "shear","constant" or "none"
+    wind_forcing_x::String="shear"      # "channel", "double_gyre", "shear","constant" or "none"
     wind_forcing_y::String="constant"   # "channel", "double_gyre", "shear","constant" or "none"
     Fx0::Real=0.12                      # wind stress strength [Pa] in x-direction
     Fy0::Real=0.0                       # wind stress strength [Pa] in y-direction
-    seasonal_wind_x::Bool=false         # Change the wind stress with a sine of frequency ωFx,ωFy
+    seasonal_wind_x::Bool=true          # Change the wind stress with a sine of frequency ωFx,ωFy
     seasonal_wind_y::Bool=false         # same for y-component
-    ωFx::Real=1.0                       # frequency [1/year] for x component
-    ωFy::Real=1.0                       # frequency [1/year] for y component
+    ωFx::Real=2                         # frequency [1/year] for x component
+    ωFy::Real=2                         # frequency [1/year] for y component
 
     # BOTTOM TOPOGRAPHY OPTIONS
     topography::String="ridges"         # "ridge", "seamount", "flat", "ridges", "bathtub"
+    topo_ridges_positions::Vector = [0.05,0.25,0.45,0.9]
     topo_height::Real=100.               # height of seamount [m]
     topo_width::Real=300e3              # horizontal scale [m] of the seamount
 
@@ -238,29 +239,28 @@ Creates a Parameter struct with following options and default values
 
     # TRACER ADVECTION
     tracer_advection::Bool=true         # yes?
-    tracer_relaxation::Bool=false       # yes?
+    tracer_relaxation::Bool=true        # yes?
     tracer_consumption::Bool=false      # yes?
-    tracer_pumping::Bool=false          # yes?
-    injection_region::String="west"     # "west", "south", "rect" or "flat"
-    sst_initial::String="south"         # "west", "south", "rect", "flat" or "restart"
+    sst_initial::String="waves"         # "west", "south", "linear", "waves","rect", "flat" or "restart"
     sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                         # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
     Uadv::Real=0.2                      # Velocity scale [m/s] for tracer advection
-    SSTmax::Real=1.                     # tracer (sea surface temperature) max for restoring
-    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for restoring
-    τSST::Real=500.                     # tracer restoring time scale [days]
-    jSST::Real=365.                     # tracer consumption [days]
-    SST_λ0::Real=222e3                  # [m] transition position of relaxation timescale
-    SST_λs::Real=111e3                  # [m] transition width of relaxation timescale
-    SST_γ0::Real=8.35                   # [days] injection time scale
+    SSTmax::Real=1.                     # tracer (sea surface temperature) max for initial conditions
+    SSTmin::Real=-1.                    # tracer (sea surface temperature) min for initial conditions
+    τSST::Real=100                      # tracer restoring time scale [days]
+    jSST::Real=365                      # tracer consumption [days]
     SSTw::Real=5e5                      # width [m] of the tangent used for the IC and interface relaxation
     SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
+    SSTwaves_ny::Real=4                 # wave crests/troughs in y
+    SSTwaves_nx::Real=SSTwaves_ny*L_ratio  # wave crests/troughs in x
+    SSTwaves_p::Real=1/2                # power for rectangles (p<1)/smootheness(p>=1) of waves
 
     # OUTPUT OPTIONS
     output::Bool=false                  # netcdf output?
     output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? "du","dv","dη" also allowed
     output_dt::Real=24                  # output time step [hours]
     outpath::String=pwd()               # path to output folder
+    compression_level::Int=3            # compression level
 
     # INITIAL CONDITIONS
     initial_cond::String="rest"         # "rest" or "ncfile" for restart from file

src/feedback.jl

@@ -56,15 +56,25 @@ function duration_estimate(feedback::Feedback,S::ModelSetup)
     end
 end
 
+"""Returns the number of NaNs in an array."""
+function countnans(A::AbstractArray)
+    n_nan = 0
+    for a in A
+        n_nan += isnan(a) ? 1 : 0
+    end
+    return n_nan
+end
+
 """Returns a boolean whether the prognostic variables contains a NaN."""
 function nan_detection!(Prog::PrognosticVars,feedback::Feedback)
 
-    #TODO include a check for Posits
-    #TODO include check for sst
-
     @unpack u,v,η,sst = Prog
 
-    n_nan = sum(isnan.(u)) + sum(isnan.(v)) + sum(isnan.(η)) + sum(isnan.(sst))
+    n_nan = countnans(u)
+    n_nan += countnans(v)
+    n_nan += countnans(η)
+    n_nan += countnans(sst)
+
     if n_nan > 0
         feedback.nans_detected = true
     end

src/forcing.jl

@@ -167,24 +167,24 @@ function Ridge(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
     return T.(Hx),T.(Hy)
 end
 
-"""Same as Ridge() but for 3 ridges at 1/4,1/2,3/4 of the domain."""
+"""Same as Ridge() but for n ridges at various x positions."""
 function Ridges(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
 
     @unpack x_T_halo,y_T_halo,Lx,Ly = G
     @unpack topo_width,topo_height,H = P
+    @unpack topo_ridges_positions = P
+    n_ridges = length(topo_ridges_positions)
 
     xx_T,yy_T = meshgrid(x_T_halo,y_T_halo)
 
-    # bumps in x direction
-    # shift slightly left/right to avoid a symmetric solution
-    b0x = exp.(-(xx_T.^2)/(2*topo_width^2))
-    b1x = exp.(-((xx_T .- 0.99*Lx/4).^2)/(2*topo_width^2))
-    b2x = exp.(-((xx_T .- 1.01*Lx/2).^2)/(2*topo_width^2))
-    b3x = exp.(-((xx_T .- 0.99*3*Lx/4).^2)/(2*topo_width^2))
-    b4x = exp.(-((xx_T .- Lx).^2)/(2*topo_width^2))
+    # loop over bumps in x direction
+    R = zero(xx_T) .+ H
 
-    th = topo_height    # for convenience
-    return T.(H .- th*b0x .- th*b1x .- th*b2x .- th*b3x .- th*b4x)
+    for i in 1:n_ridges
+        R .-= topo_height*exp.(-((xx_T .- topo_ridges_positions[i]*Lx).^2)/(2*topo_width^2))
+    end
+
+    return T.(R)
 end
 
 """Returns a matrix of constant water depth H."""
@@ -213,5 +213,5 @@ end
 
 """Time evolution of forcing."""
 function Ftime(::Type{T},t::Int,ω::Real) where {T<:AbstractFloat}
-    return T(sin(ω*t))
-end
+    return convert(T,sin(ω*t))
+end

src/ghost_points.jl

@@ -22,10 +22,10 @@ function add_halo(  u::Array{T,2},
     sst = cat(zeros(T,nx+2*halosstx,halossty),sst,zeros(T,nx+2*halosstx,halossty),dims=2)
 
     # SCALING
-    @unpack scale = S.constants
+    @unpack scale,scale_sst = S.constants
     u *= scale
     v *= scale
-    sst *= scale
+    sst *= scale_sst
 
     ghost_points!(u,v,η,S)
     ghost_points_sst!(sst,S)
@@ -40,13 +40,13 @@ function remove_halo(   u::Array{T,2},
                         S::ModelSetup) where {T<:AbstractFloat}
 
     @unpack halo,haloη,halosstx,halossty = S.grid
-    @unpack scale_inv = S.constants
+    @unpack scale_inv,scale_sst_inv = S.constants
 
     # undo scaling as well
     ucut = scale_inv*u[halo+1:end-halo,halo+1:end-halo]
     vcut = scale_inv*v[halo+1:end-halo,halo+1:end-halo]
     ηcut = η[haloη+1:end-haloη,haloη+1:end-haloη]
-    sstcut = scale_inv*sst[halosstx+1:end-halosstx,halossty+1:end-halossty]
+    sstcut = scale_sst_inv*sst[halosstx+1:end-halosstx,halossty+1:end-halossty]
 
     return ucut,vcut,ηcut,sstcut
 end

src/grid.jl

@@ -148,8 +148,8 @@ function meshgrid(x::AbstractVector,y::AbstractVector)
     m,n = length(x),length(y)
 
     # preallocate preserving the data type of x,y
-    xx = zeros(typeof(x[1]),m,n)
-    yy = zeros(typeof(y[1]),m,n)
+    xx = zeros(eltype(x),m,n)
+    yy = zeros(eltype(y),m,n)
 
     for i in 1:m
         xx[i,:] .= x[i]

src/initial_conditions.jl

@@ -119,6 +119,7 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
     ## SST
 
     @unpack SSTmin, SSTmax, SSTw, SSTϕ = S.parameters
+    @unpack SSTwaves_nx,SSTwaves_ny,SSTwaves_p = S.parameters
     @unpack sst_initial,scale = S.parameters
     @unpack x_T,y_T,Lx,Ly = S.grid
 
@@ -128,6 +129,10 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
         sst = (SSTmin+SSTmax)/2 .+ tanh.(2π*(Ly/(4*SSTw))*(yy_T/Ly .- SSTϕ))*(SSTmin-SSTmax)/2
     elseif sst_initial == "west"
         sst = (SSTmin+SSTmax)/2 .+ tanh.(2π*(Lx/(4*SSTw))*(xx_T/Lx .- SSTϕ))*(SSTmin-SSTmax)/2
+    elseif sst_initial == "linear"
+        sst = SSTmin .+ yy_T/Ly*(SSTmax-SSTmin)
+    elseif sst_initial == "waves"
+        sst = waves(xx_T/Lx,yy_T/Ly,SSTwaves_nx,SSTwaves_ny,SSTwaves_p)
     elseif sst_initial == "flat"
         sst = fill(SSTmin,size(xx_T))
     elseif sst_initial == "rect"
@@ -158,3 +163,13 @@ function initial_conditions(::Type{T},S::ModelSetup) where {T<:AbstractFloat}
 
     return PrognosticVars{T}(u,v,η,sst)
 end
+
+"""Create a wave-checkerboard pattern over xx,yy like a nx x ny checkerboard.
+p is the power to which the waves are raised. Choose p<1 for rectangles, and
+p > 1 for more smootheness."""
+function waves(xx::AbstractMatrix,yy::AbstractMatrix,nx::Real,ny::Real,p::Real)
+    @boundscheck size(xx) == size(yy) || throw(BoundsError())
+    w = sin.(nx*π*xx) .* sin.(ny*π*yy)
+    s = sign.(w)
+    return s.*abs.(w).^p
+end

src/output.jl

@@ -90,7 +90,7 @@ function output_nc!(i::Int,
 
         @unpack halo,haloη,halosstx,halossty = S.grid
         @unpack f_q,ep,dtint = S.grid
-        @unpack scale,scale_inv = S.constants
+        @unpack scale,scale_inv,scale_sst_inv = S.constants
 
         # CUT OFF HALOS
         # As output is before copyto!(u,u0), take u0,v0,η0
@@ -109,7 +109,7 @@ function output_nc!(i::Int,
             NetCDF.putvar(ncs.η,"eta",η,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.sst != nothing
-            @views sst = Float32.(scale_inv*Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
+            @views sst = Float32.(scale_sst_inv*Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty])
             NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
         end
         if ncs.q != nothing

src/preallocate.jl

@@ -421,6 +421,7 @@ end
     vinterp::Array{T,2} = zeros(T,nx,ny)                    # v interpolated on mid-point xd,yd
 
     ssti::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst interpolated on departure points
+    sst_ref::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst initial conditions for relaxation
 end
 
 """Generator function for SemiLagrange VarCollection."""

src/time_integration.jl

@@ -35,6 +35,9 @@ function time_integration(  Prog::PrognosticVars{Tprog},
     copyto!(v0,v)
     copyto!(η0,η)
 
+    # store initial conditions of sst for relaxation
+    copyto!(Diag.SemiLagrange.sst_ref,sst)
+
     # feedback, output initialisation and storing initial conditions
     feedback = feedback_init(S)
     netCDFfiles = NcFiles(feedback,S)