Merge pull request #112 from milankl/whatsnext

Whatsnext

Author: Milan K

.gitignore

@@ -2,3 +2,5 @@
 *.jld
 *.nc
 nohup.out
+progress.txt
+parameter.txt

src/Constants.jl

@@ -17,6 +17,7 @@ struct Constants{T<:AbstractFloat}
     rSST::T                 # tracer restoring timescale
     jSST::T                 # tracer consumption timescale
     SSTmin::T               # tracer minimum
+    ωyr::Float64            # # frequency [1/s] of Kelvin pumping (including 2π)
 end
 
 """Generator function for the mutable struct Constants."""
@@ -53,5 +54,8 @@ function Constants{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
     jSST = T(G.dtadvint/(P.jSST*3600*24))    # tracer consumption [1]
     SSTmin = T(P.SSTmin)
 
-    return Constants{T}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin)
+    # SURFACE FORCING
+    ωyr = -2π*P.ωyr/24/365/3600
+
+    return Constants{T}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin,ωyr)
 end

src/Continuity.jl

@@ -1,7 +1,16 @@
 """Continuity equation's right-hand side with surface relaxation
 -∂x(uh) - ∂y(vh) + γ*(η_ref - η)."""
-function continuity_surf_forc!(dη,dUdx,dVdy,η,η_ref,Fη,t)
-    m,n = size(dη) .- (2*haloη,2*haloη)
+function continuity_surf_relax!(η::AbstractMatrix,
+                                Diag::DiagnosticVars,
+                                S::ModelSetup,
+                                t::Int) where {T<:AbstractFloat}
+
+    @unpack dη = Diag.Tendencies
+    @unpack dUdx,dVdy = Diag.VolumeFluxes
+    @unpack η_ref = S.forcing
+    @unpack γ = S.constants
+
+    m,n = size(dη) .- (2,2)
     @boundscheck (m,n+2) == size(dUdx) || throw(BoundsError())
     @boundscheck (m+2,n) == size(dVdy) || throw(BoundsError())
     @boundscheck (m+2,n+2) == size(η) || throw(BoundsError())
@@ -15,14 +24,24 @@ function continuity_surf_forc!(dη,dUdx,dVdy,η,η_ref,Fη,t)
 end
 
 """Continuity equation's right-hand side with time&space dependent forcing."""
-function continuity_forcing!(dη,dUdx,dVdy,η,η_ref,Fη,t)
-    m,n = size(dη) .- (2*haloη,2*haloη)
+function continuity_forcing!(   ::Type{T},
+                                η::AbstractMatrix,
+                                Diag::DiagnosticVars,
+                                S::ModelSetup,
+                                t::Int) where {T<:AbstractFloat}
+
+    @unpack dη = Diag.Tendencies
+    @unpack dUdx,dVdy = Diag.VolumeFluxes
+    @unpack Fη = S.forcing
+
+    m,n = size(dη) .- (2,2)         # cut off halo
     @boundscheck (m,n+2) == size(dUdx) || throw(BoundsError())
     @boundscheck (m+2,n) == size(dVdy) || throw(BoundsError())
     @boundscheck (m,n) == size(Fη) || throw(BoundsError())
 
     # avoid recomputation
-    Fηtt = Fηt(t)
+    @unpack ωyr = S.constants
+    Fηtt = Fηt(T,t,ωyr)
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
@@ -32,7 +51,11 @@ function continuity_forcing!(dη,dUdx,dVdy,η,η_ref,Fη,t)
 end
 
 """Continuity equation's right-hand side -∂x(uh) - ∂y(vh) without forcing."""
-function continuity!(Diag::DiagnosticVars,S::ModelSetup)
+function continuity_itself!(::Type{T},
+                            η::AbstractMatrix,
+                            Diag::DiagnosticVars,
+                            S::ModelSetup,
+                            t::Int) where {T<:AbstractFloat}
 
     @unpack dη = Diag.Tendencies
     @unpack dUdx,dVdy = Diag.VolumeFluxes
@@ -48,10 +71,18 @@ function continuity!(Diag::DiagnosticVars,S::ModelSetup)
     end
 end
 
-# if surface_relax
-#     continuity! = continuity_surf_forc!
-# elseif surface_forcing
-#     continuity! = continuity_forcing!
-# else
-#     continuity! = continuity_itself!
-# end
+
+"""Transit function to call the specified continuity function."""
+function continuity!(   η::Array{T,2},
+                        Diag::DiagnosticVars,
+                        S::ModelSetup,
+                        t::Int) where {T<:AbstractFloat}
+
+    if S.parameters.surface_relax
+        continuity_surf_relax!(T,η,Diag,S,t)
+    elseif S.parameters.surface_forcing
+        continuity_forcing!(T,η,Diag,S,t)
+    else
+        continuity_itself!(T,η,Diag,S,t)
+    end
+end

src/DefaultParameters.jl

@@ -30,9 +30,9 @@
     surface_relax::Bool=false           # yes?
     t_relax::Real=100.                  # time scale of the relaxation [days]
     η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
-    η_refw::Real=100e3                  # width [m] of the tangent used for the interface relaxation
+    η_refw::Real=50e3                  # width [m] of the tangent used for the interface relaxation
 
-    # SURFACE FORCING
+    # SURFACE FORCING (Currently only Kelvin wave pumping at Eq.)
     surface_forcing::Bool=false         # yes?
     ωyr::Real=1.0                       # (annual) frequency [1/year]
     A::Real=3e-5                        # Amplitude [m/s]

src/Feedback.jl

@@ -97,6 +97,11 @@ function feedback_init(S::ModelSetup)
         write(txt,"Boundary conditions are $bc with lbc=$α.\n")
         write(txt,"Number format is "*string(T)*".\n")
         write(txt,"\nAll data will be stored in $runpath\n")
+
+        # Parameter.txt
+        ptxt = open(joinpath(runpath,"parameter.txt"),"w")
+        print(ptxt,S.parameters)
+        close(ptxt)
     else
         println("Starting Juls on "*Dates.format(now(),Dates.RFC1123Format)*" without output.")
         txt = nothing

src/Forcing.jl

@@ -3,6 +3,7 @@ struct Forcing{T<:AbstractFloat}
     Fy::Array{T,2}
     H::Array{T,2}
     η_ref::Array{T,2}
+    Fη::Array{T,2}
     #sst_ref::Array{T,2}
     #sst_γ::Array{T,2}
 end
@@ -19,6 +20,8 @@ function Forcing{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
         Fx,_ = DoubleGyreWind(T,P,G)
     elseif wind_forcing_x == "constant"
         Fx,_ = ConstantWind(T,P,G)
+    elseif wind_forcing_x == "none"
+        Fx,_ = NoWind(T,P,G)
     end
 
     if wind_forcing_y == "channel"
@@ -29,6 +32,8 @@ function Forcing{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
         _,Fy = DoubleGyreWind(T,P,G)
     elseif wind_forcing_y == "constant"
         _,Fy = ConstantWind(T,P,G)
+    elseif wind_forcing_x == "none"
+        _,Fy = NoWind(T,P,G)
     end
 
     if topography == "ridge"
@@ -42,8 +47,9 @@ function Forcing{T}(P::Parameter,G::Grid) where {T<:AbstractFloat}
     end
 
     η_ref = InterfaceRelaxation(T,P,G)
+    Fη = KelvinPump(T,P,G)
 
-    return Forcing{T}(Fx,Fy,H,η_ref)
+    return Forcing{T}(Fx,Fy,H,η_ref,Fη)
 end
 
 """Returns the constant forcing matrices Fx,Fy that vary only meriodionally/zonally
@@ -104,6 +110,17 @@ function ConstantWind(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
     return Fx,Fy
 end
 
+"""Returns constant in in space forcing matrices Fx,Fy."""
+function NoWind(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
+
+    @unpack nux,nuy,nvx,nvy = G
+
+    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
+    Fx = zeros(T,nux,nuy)
+    Fy = zeros(T,nvx,nvy)
+
+    return Fx,Fy
+end
 
 """Returns a reference state for Newtonian cooling/surface relaxation shaped as a
 hyperbolic tangent to force the continuity equation."""
@@ -176,3 +193,31 @@ function FlatBottom(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
     @unpack H = P
     return fill(T(H),(nx+2*haloη,ny+2*haloη))
 end
+
+"""Returns Kelvin wave pumping forcing of the continuity equation at the equator."""
+function KelvinPump(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
+
+    @unpack x_T,y_T,Lx,Ly = G
+    @unpack c,β,R,ϕ,Δ = G
+    @unpack A = P
+
+    xx_T,yy_T = meshgrid(x_T,y_T)
+
+    # y-coordinate of the Equator
+    mϕ = 2π*R/360.      # meters per degree latitude
+    y_eq = Ly/2 - ϕ*mϕ
+    y_15S = Ly/2 - (ϕ+15)*mϕ
+    y_15N = Ly/2 - (ϕ-15)*mϕ
+
+    Fη = A*Δ*exp.(-β*(yy_T.-y_eq).^2/(2c))
+
+    Fη[yy_T .< y_15S] .= 0.0
+    Fη[yy_T .> y_15N] .= 0.0
+
+    return T.(Fη)
+end
+
+"""Time evolution of surface forcing."""
+function Fηt(::Type{T},t::Int,ωyr::AbstractFloat) where {T<:AbstractFloat}
+    return T(sin(ωyr*t))
+end

src/TimeIntegration.jl

@@ -47,7 +47,7 @@ function time_integration(  ::Type{T},
                 ghost_points!(u1,v1,η1,S)
             end
 
-            rhs!(u1,v1,η1,Diag,S)
+            rhs!(u1,v1,η1,Diag,S,t)
 
             if rki < RKo
                 caxb!(u1,u,RKbΔt[rki],du)   #u1 .= u .+ RKb[rki]*Δt*du

src/rhs.jl

@@ -3,14 +3,15 @@ function rhs!(  u::AbstractMatrix,
                 v::AbstractMatrix,
                 η::AbstractMatrix,
                 Diag::DiagnosticVars,
-                S::ModelSetup)
+                S::ModelSetup,
+                t::Int)
 
     @unpack dynamics = S.parameters
 
     if dynamics == "linear"
-        rhs_linear!(u,v,η,Diag,S)
+        rhs_linear!(u,v,η,Diag,S,t)
     else
-        rhs_nonlinear!(u,v,η,Diag,S)
+        rhs_nonlinear!(u,v,η,Diag,S,t)
     end
 end
 
@@ -25,7 +26,8 @@ function rhs_nonlinear!(u::AbstractMatrix,
                         v::AbstractMatrix,
                         η::AbstractMatrix,
                         Diag::DiagnosticVars,
-                        S::ModelSetup)
+                        S::ModelSetup,
+                        t::Int)
 
     @unpack h,h_u,h_v,U,V,dUdx,dVdy = Diag.VolumeFluxes
     @unpack H = S.forcing
@@ -55,14 +57,13 @@ function rhs_nonlinear!(u::AbstractMatrix,
     ∂x!(dpdx,p)
     ∂y!(dpdy,p)
 
-    #TODO check whether the order of PV,PVadvection etc. is correct
     # Potential vorticity and advection thereof
     PVadvection!(Diag,S)
 
     # adding the terms
     momentum_u!(Diag,S)
     momentum_v!(Diag,S)
-    continuity!(Diag,S)
+    continuity!(η,Diag,S,t)
 end
 
 """Tendencies du,dv,dη of
@@ -76,7 +77,8 @@ function rhs_linear!(   u::AbstractMatrix,
                         v::AbstractMatrix,
                         η::AbstractMatrix,
                         Diag::DiagnosticVars,
-                        S::ModelSetup)
+                        S::ModelSetup,
+                        t::Int)
 
     @unpack h,h_u,h_v,U,V,dUdx,dVdy = Diag.VolumeFluxes
     @unpack g = S.constants
@@ -106,7 +108,7 @@ function rhs_linear!(   u::AbstractMatrix,
     # adding the terms
     momentum_u!(Diag,S)
     momentum_v!(Diag,S)
-    continuity!(Diag,S)
+    continuity!(η,Diag,S,t)
 end
 
 """ Update advective and Coriolis tendencies."""