Merge pull request #126 from milankl/KelvinPump

surface forcing latitude adjustment

Author: Milan K

src/Constants.jl

@@ -55,7 +55,7 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
     SSTmin = T(P.SSTmin)
 
     # SURFACE FORCING
-    ωyr = -2π*P.ωyr/24/365/3600
+    ωyr = -2π*P.ωyr/24/365.25/3600
 
     return Constants{T,Tprog}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin,ωyr)
 end

src/DefaultParameters.jl

@@ -14,7 +14,7 @@
     g::Real=10.                         # gravitational acceleration [m/s]
     H::Real=500.                        # layer thickness at rest [m]
     ρ::Real=1e3                         # water density [kg/m^3]
-    ϕ::Real=45.                         # central latitue of the domain (for coriolis) [°]
+    ϕ::Real=45.                         # central latitude of the domain (for coriolis) [°]
     ω::Real=2π/(24*3600)                # Earth's angular frequency [s^-1]
     R::Real=6.371e6                     # Earth's radius [m]
 
@@ -35,10 +35,12 @@
     η_refh::Real=5.                     # height difference [m] of the interface relaxation profile
     η_refw::Real=50e3                   # width [m] of the tangent used for the interface relaxation
 
-    # SURFACE FORCING (Currently only Kelvin wave pumping at Eq.)
+    # SURFACE FORCING
     surface_forcing::Bool=false         # yes?
     ωyr::Real=1.0                       # (annual) frequency [1/year]
     A::Real=3e-5                        # Amplitude [m/s]
+    ϕk::Real=ϕ                          # Central latitude of Kelvin wave pumping
+    wk::Real=10e3                       # width [m] in y of Gaussian used for surface forcing
 
     # TIME STEPPING OPTIONS
     RKo::Int=4                          # Order of the RK time stepping scheme (3 or 4)

src/Feedback.jl

@@ -139,6 +139,7 @@ function feedback!(Prog::PrognosticVars,feedback::Feedback,S::ModelSetup)
             if feedback.nans_detected
                 println(" NaNs detected at time step $i")
                 if output
+                    @unpack progress_txt = feedback
                     write(progress_txt," NaNs detected at time step $i")
                     flush(progress_txt)
                 end

src/Forcing.jl

@@ -194,25 +194,19 @@ function FlatBottom(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
     return fill(T(H),(nx+2*haloη,ny+2*haloη))
 end
 
-"""Returns Kelvin wave pumping forcing of the continuity equation at the equator."""
+"""Returns Kelvin wave pumping forcing of the continuity equation."""
 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
+    @unpack R,ϕ,Δ = G
+    @unpack A,ϕk,wk = 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ϕ
+    mϕ = 2π*R/360.          # meters per degree latitude
+    y0 = Ly/2 - (ϕ-ϕk)*mϕ   # y[m] for central latitude of pumping
 
-    Fη = A*Δ*exp.(-β*(yy_T.-y_eq).^2/(2c))
-
-    Fη[yy_T .< y_15S] .= 0.0
-    Fη[yy_T .> y_15N] .= 0.0
+    Fη = A*Δ*exp.(-(yy_T.-y0).^2/(2*wk^2))
 
     return T.(Fη)
 end