relaxation, consumption, scale_sst

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/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/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)

src/tracer_advection.jl

@@ -6,6 +6,7 @@ function tracer!(   i::Integer,
                     S::ModelSetup)
 
     @unpack tracer_advection, tracer_consumption = S.parameters
+    @unpack tracer_relaxation = S.parameters
     @unpack nadvstep_half,nadvstep = S.grid
 
     # mid point (in time) velocity for the advective time step
@@ -19,17 +20,17 @@ function tracer!(   i::Integer,
     if tracer_advection && (i % nadvstep) == 0
 
         @unpack sst = Prog
-        @unpack ssti = Diag.SemiLagrange
+        @unpack ssti,sst_ref = Diag.SemiLagrange
 
         # convert to type T for mixed precision
         sstrhs = convert(Diag.PrognosticVarsRHS.sst,sst)
 
         departure!(u,v,Diag,S)
         adv_sst!(sstrhs,Diag,S)
 
-        # if tracer_relaxation
-        #     tracer_relax!(ssti,sst_ref,SSTγ)
-        # end
+        if tracer_relaxation
+            tracer_relax!(ssti,sst_ref,S)
+        end
 
         if tracer_consumption
             tracer_consumption!(ssti,S)
@@ -245,44 +246,35 @@ function bilin(f00::T,f10::T,f01::T,f11::T,x::T,y::T) where {T<:AbstractFloat}
     return f00*(oone-x)*(oone-y) + f10*x*(oone-y) + f01*(oone-x)*y + f11*x*y
 end
 
-# """Tracer relaxation."""
-# function tracer_relax!(sst::AbstractMatrix,sst_ref::AbstractMatrix,SSTγ::AbstractMatrix)
-#     m,n = size(sst)
-#     @boundscheck (m-2*halosstx,n-2*halossty) == size(sst_ref) || throw(BoundsError())
-#     @boundscheck (m-2*halosstx,n-2*halossty) == size(SSTγ) || throw(BoundsError())
-#
-#     @inbounds for j ∈ 1+halossty:n-halossty
-#         for i ∈ 1+halosstx:m-halosstx
-#             sst[i,j] += SSTγ[i-halosstx,j-halossty]*(sst_ref[i-halosstx,j-halossty] - sst[i,j])
-#         end
-#     end
-# end
-
-"""Tracer consumption via relaxation back to ."""
+"""Tracer relaxation."""
+function tracer_relax!( sst::AbstractMatrix,
+                        sst_ref::AbstractMatrix,
+                        S::ModelSetup)
+    @unpack τSST = S.constants
+    @unpack halosstx,halossty = S.grid
+
+    m,n = size(sst)
+    @boundscheck size(sst) == size(sst_ref) || throw(BoundsError())
+    
+    for j ∈ 1+halossty:n-halossty
+        for i ∈ 1+halosstx:m-halosstx
+            sst[i,j] += τSST*(sst_ref[i,j] - sst[i,j])
+        end
+    end
+end
+
+"""Tracer consumption via relaxation back to 0."""
 function tracer_consumption!(   sst::Array{T,2},
                                 S::ModelSetup) where {T<:AbstractFloat}
 
-    @unpack jSST,SSTmin = S.constants
+    @unpack jSST = S.constants
     @unpack halosstx,halossty = S.grid
 
     m,n = size(sst)
 
-    @inbounds for j ∈ 1+halossty:n-halossty
+    for j ∈ 1+halossty:n-halossty
         for i ∈ 1+halosstx:m-halosstx
-            sst[i,j] += jSST*(SSTmin - sst[i,j])
+            sst[i,j] -= jSST*sst[i,j]
         end
     end
-end
-
-# """Spatially dependent relaxation time scale."""
-# function sst_γ(x::AbstractVector,y::AbstractVector)
-#     xx,yy = meshgrid(x,y)
-#
-#     # convert from days to one over 1/s and include adv time step
-#     γ0 = dtadvint/(SST_γ0*3600*24)
-#
-#     x10E = 10*m_per_lat()   # assume Equator: lat/lon equivalence
-#     γ = γ0/2 .* (1 .- tanh.((xx.-SST_λ0)./SST_λs))
-#     γ[xx .> x10E] .= 0.0
-#     return Numtype.(γ)
-# end
+end