diff --git a/.cirrus.yml b/.cirrus.yml
deleted file mode 100644
index 80a837b..0000000
--- a/.cirrus.yml
+++ /dev/null
@@ -1,17 +0,0 @@
-freebsd_instance:
-  image: freebsd-12-0-release-amd64
-task:
-  name: FreeBSD
-  env:
-    matrix:
-      - JULIA_VERSION: 1.0
-      - JULIA_VERSION: 1.3
-      - JULIA_VERSION: nightly
-  install_script:
-    - sh -c "$(fetch https://raw.githubusercontent.com/ararslan/CirrusCI.jl/master/bin/install.sh -o -)"
-  build_script:
-    - cirrusjl build
-  test_script:
-    - cirrusjl test
-  coverage_script:
-    - cirrusjl coverage codecov coveralls
diff --git a/src/Constants.jl b/src/Constants.jl
index 390319b..d45d52b 100644
--- a/src/Constants.jl
+++ b/src/Constants.jl
@@ -17,7 +17,9 @@ struct Constants{T<:AbstractFloat,Tprog<: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π)
+    ω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π
 end
 
 """Generator function for the mutable struct Constants."""
@@ -54,8 +56,12 @@ function Constants{T,Tprog}(P::Parameter,G::Grid) where {T<:AbstractFloat,Tprog<
     jSST = T(G.dtadvint/(P.jSST*3600*24))    # tracer consumption [1]
     SSTmin = T(P.SSTmin)
 
-    # SURFACE FORCING
-    ωyr = -2π*P.ωyr/24/365.25/3600
+    # 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
 
-    return Constants{T,Tprog}(RKaΔt,RKbΔt,one_minus_α,g,cD,rD,γ,cSmag,νB,rSST,jSST,SSTmin,ωyr)
+    return Constants{T,Tprog}(  RKaΔt,RKbΔt,one_minus_α,
+                                g,cD,rD,γ,cSmag,νB,rSST,
+                                jSST,SSTmin,ωFη,ωFx,ωFy)
 end
diff --git a/src/Continuity.jl b/src/Continuity.jl
index 24794b1..547609a 100644
--- a/src/Continuity.jl
+++ b/src/Continuity.jl
@@ -38,12 +38,12 @@ function continuity_forcing!(   Diag::DiagnosticVars{T,Tprog},
     @boundscheck (m,n) == size(Fη) || throw(BoundsError())
 
     # avoid recomputation
-    @unpack ωyr = S.constants
-    Fηtt = Fηt(T,t,ωyr)
+    @unpack ωFη = S.constants
+    Fηt = Ftime(T,t,ωFη)
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(Fηtt*Fη[i,j])
+            dη[i+1,j+1] = -(Tprog(dUdx[i,j+1]) + Tprog(dVdy[i+1,j])) + Tprog(Fηt*Fη[i,j])
         end
     end
 end
diff --git a/src/DefaultParameters.jl b/src/DefaultParameters.jl
index 7d8758b..eba28f6 100644
--- a/src/DefaultParameters.jl
+++ b/src/DefaultParameters.jl
@@ -23,6 +23,10 @@
     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_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
 
     # BOTTOM TOPOGRAPHY OPTIONS
     topography::String="ridge"          # "ridge", "seamount", "flat", "ridges", "bathtub"
@@ -37,7 +41,7 @@
 
     # SURFACE FORCING
     surface_forcing::Bool=false         # yes?
-    ωyr::Real=1.0                       # (annual) frequency [1/year]
+    ωFη::Real=1.0                       # frequency [1/year] for surfance forcing
     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
@@ -114,7 +118,6 @@
     @assert topo_width > 0.0    "topo_width has to be >0, $topo_width given."
     @assert t_relax > 0.0       "t_relax has to be >0, $t_relax given."
     @assert η_refw > 0.0        "η_refw has to be >0, $η_refw given."
-    @assert ωyr > 0.0           "ωyr has to be >0, $ωyr given."
     @assert RKo in [3,4]        "RKo has to be 3 or 4, $RKo given."
     @assert Ndays > 0.0         "Ndays has to be >0, $Ndays given."
     @assert nstep_diff > 0      "nstep_diff has to be >0, $nstep_diff given."
@@ -144,7 +147,7 @@ Creates a Parameter struct with following options and default values
     T::DataType=Float32                 # number format
 
     Tprog::DataType=T                   # number format for prognostic variables
-    Tcomm::DataType=T                   # number format for ghost-point copies
+    Tcomm::DataType=Tprog               # number format for ghost-point copies
 
     # DOMAIN RESOLUTION AND RATIO
     nx::Int=100                         # number of grid cells in x-direction
@@ -155,7 +158,7 @@ Creates a Parameter struct with following options and default values
     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]
 
@@ -164,6 +167,10 @@ Creates a Parameter struct with following options and default values
     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_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
 
     # BOTTOM TOPOGRAPHY OPTIONS
     topography::String="ridge"          # "ridge", "seamount", "flat", "ridges", "bathtub"
@@ -176,10 +183,12 @@ Creates a Parameter struct with following options and default values
     η_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]
+    ωFη::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)
@@ -230,7 +239,7 @@ Creates a Parameter struct with following options and default values
 
     # OUTPUT OPTIONS
     output::Bool=false                  # netcdf output?
-    output_vars::Array{String,1}=["u","v","η","sst","q","ζ"]  # which variables to output?
+    output_vars::Array{String,1}=["u","v","η","sst","q","ζ"]  # which variables to output? "du","dv","dη" also allowed
     output_dt::Real=6                   # output time step [hours]
     outpath::String=pwd()               # path to output folder
 
@@ -240,5 +249,8 @@ Creates a Parameter struct with following options and default values
     init_run_id::Int=0                  # run id for restart from run number
     init_starti::Int=-1                 # timestep to start from (-1 meaning last)
     get_id_mode::String="continue"      # How to determine the run id: "continue" or "fill"
+    run_id::Int=-1                      # Output with a specific run id
+    init_interpolation::Bool=true       # Interpolate the initial conditions in case grids don't match?
+
 """
 Parameter
diff --git a/src/Forcing.jl b/src/Forcing.jl
index 114a8c7..d20f5e8 100644
--- a/src/Forcing.jl
+++ b/src/Forcing.jl
@@ -211,7 +211,7 @@ function KelvinPump(::Type{T},P::Parameter,G::Grid) where {T<:AbstractFloat}
     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))
+"""Time evolution of forcing."""
+function Ftime(::Type{T},t::Int,ω::Real) where {T<:AbstractFloat}
+    return T(sin(ω*t))
 end
diff --git a/src/rhs.jl b/src/rhs.jl
index 21c8787..2cdae39 100644
--- a/src/rhs.jl
+++ b/src/rhs.jl
@@ -61,8 +61,8 @@ function rhs_nonlinear!(u::AbstractMatrix,
     PVadvection!(Diag,S)
 
     # adding the terms
-    momentum_u!(Diag,S)
-    momentum_v!(Diag,S)
+    momentum_u!(Diag,S,t)
+    momentum_v!(Diag,S,t)
     continuity!(η,Diag,S,t)
 end
 
@@ -106,8 +106,8 @@ function rhs_linear!(   u::AbstractMatrix,
     fu!(qhu,f_v,u_v)
 
     # adding the terms
-    momentum_u!(Diag,S)
-    momentum_v!(Diag,S)
+    momentum_u!(Diag,S,t)
+    momentum_v!(Diag,S,t)
     continuity!(η,Diag,S,t)
 end
 
@@ -281,7 +281,9 @@ function fu!(   qhu::AbstractMatrix,
 end
 
 """Sum up the tendencies of the non-diffusive right-hand side for the u-component."""
-function momentum_u!(Diag::DiagnosticVars{T,Tprog},S::ModelSetup) where {T,Tprog}
+function momentum_u!(   Diag::DiagnosticVars{T,Tprog},
+                        S::ModelSetup,
+                        t::Int) where {T,Tprog}
 
     @unpack du = Diag.Tendencies
     @unpack qhv = Diag.Vorticity
@@ -294,15 +296,24 @@ function momentum_u!(Diag::DiagnosticVars{T,Tprog},S::ModelSetup) where {T,Tprog
     @boundscheck (m+2-ep,n+2) == size(dpdx) || throw(BoundsError())
     @boundscheck (m,n) == size(Fx) || throw(BoundsError())
 
+    if S.parameters.seasonal_wind_x
+        @unpack ωFx = S.constants
+        Fxt = Ftime(T,t,ωFx)
+    else
+        Fxt = one(T)
+    end
+
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-            du[i+2,j+2] = Tprog(qhv[i,j]) - Tprog(dpdx[i+1-ep,j+1]) + Tprog(Fx[i,j])
+            du[i+2,j+2] = Tprog(qhv[i,j]) - Tprog(dpdx[i+1-ep,j+1]) + Tprog(Fxt*Fx[i,j])
         end
     end
 end
 
 """Sum up the tendencies of the non-diffusive right-hand side for the v-component."""
-function momentum_v!(Diag::DiagnosticVars{T,Tprog},S::ModelSetup) where {T,Tprog}
+function momentum_v!(   Diag::DiagnosticVars{T,Tprog},
+                        S::ModelSetup,
+                        t::Int) where {T,Tprog}
 
     @unpack dv = Diag.Tendencies
     @unpack qhu = Diag.Vorticity
@@ -314,9 +325,16 @@ function momentum_v!(Diag::DiagnosticVars{T,Tprog},S::ModelSetup) where {T,Tprog
     @boundscheck (m,n) == size(qhu) || throw(BoundsError())
     @boundscheck (m+2,n+2) == size(dpdy) || throw(BoundsError())
 
+    if S.parameters.seasonal_wind_y
+        @unpack ωFy = S.constants
+        Fyt = Ftime(T,t,ωFy)
+    else
+        Fyt = one(T)
+    end
+
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-             dv[i+2,j+2] = -Tprog(qhu[i,j]) - Tprog(dpdy[i+1,j+1]) + Tprog(Fy[i,j])
+             dv[i+2,j+2] = -Tprog(qhu[i,j]) - Tprog(dpdy[i+1,j+1]) + Tprog(Fyt*Fy[i,j])
         end
     end
 end