Merge pull request #146 from milankl/ssprk

Keep UV and PVadv on every substep

Author: milankl

src/Continuity.jl

@@ -80,8 +80,7 @@ function continuity!(   u::AbstractMatrix,
     @unpack nstep_advcor = S.grid
     @unpack time_scheme,surface_relax,surface_forcing = S.parameters
 
-    if nstep_advcor > 0 ||              # then UV wasn't calculated inside rhs!
-        time_scheme != "RK"             # then u,v changed before evaluating semi-implciti continuity
+    if time_scheme != "RK"             # then u,v changed before evaluating semi-implicit continuity
         UVfluxes!(u,v,η,Diag,S)
     end
 

src/DefaultParameters.jl

@@ -49,7 +49,7 @@
     # TIME STEPPING OPTIONS
     time_scheme::String="RK"            # Runge-Kutta ("RK") or strong-stability preserving RK ("SSPRK2","SSPRK3")
     RKo::Int=4                          # Order of the RK time stepping scheme (2, 3 or 4)
-    RKs::Int=2                          # Number of stages for SSPRK2
+    RKs::Int=3                          # Number of stages for SSPRK2
     cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
     Ndays::Real=10.0                    # number of days to integrate for
     nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
@@ -148,10 +148,9 @@ end
 """
 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=Tprog               # number format for ghost-point copies
+    T=Float32                 # number format
+    Tprog=T                   # number format for prognostic variables
+    Tcomm=Tprog               # number format for ghost-point copies
 
     # DOMAIN RESOLUTION AND RATIO
     nx::Int=100                         # number of grid cells in x-direction
@@ -195,8 +194,10 @@ Creates a Parameter struct with following options and default values
     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 (2,3 or 4)
-    cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3, 0.1 for RK2)
+    time_scheme::String="RK"            # Runge-Kutta ("RK") or strong-stability preserving RK ("SSPRK2","SSPRK3")
+    RKo::Int=4                          # Order of the RK time stepping scheme (2, 3 or 4)
+    RKs::Int=2                          # Number of stages for SSPRK2
+    cfl::Real=1.0                       # CFL number (1.0 recommended for RK4, 0.6 for RK3)
     Ndays::Real=10.0                    # number of days to integrate for
     nstep_diff::Int=1                   # diffusive part every nstep_diff time steps.
     nstep_advcor::Int=0                 # advection and coriolis update every nstep_advcor time steps.
@@ -243,7 +244,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"]  # which variables to output? q,ζ,du,dv,dη also allowed.
+    output_vars::Array{String,1}=["u","v","η","sst"]  # which variables to output? "du","dv","dη" also allowed
     output_dt::Real=6                   # output time step [hours]
     outpath::String=pwd()               # path to output folder
 

src/TimeIntegration.jl

@@ -163,7 +163,6 @@ function time_integration(  Prog::PrognosticVars{Tprog},
         if dynamics == "nonlinear" && nstep_advcor > 0 && (i % nstep_advcor) == 0
             UVfluxes!(u0rhs,v0rhs,η0rhs,Diag,S)
             advection_coriolis!(u0rhs,v0rhs,η0rhs,Diag,S)
-            PVadvection!(Diag,S)
         end
 
         # DIFFUSIVE TERMS - SEMI-IMPLICIT EULER

src/rhs.jl

@@ -32,12 +32,12 @@ function rhs_nonlinear!(u::AbstractMatrix,
     @unpack H = S.forcing
     @unpack ep = S.grid
 
-    if S.grid.nstep_advcor == 0              # then evaluate every RK substep
-        UVfluxes!(u,v,η,Diag,S)              # U,V needed for PV advection and in the continuity equation
+    UVfluxes!(u,v,η,Diag,S)              # U,V needed for PV advection and in the continuity equation
+    if S.grid.nstep_advcor == 0              # evaluate every RK substep
         advection_coriolis!(u,v,η,Diag,S)    # PV and non-linear Bernoulli terms
-        PVadvection!(Diag,S)                 # advect the PV with u,v
     end
-
+    PVadvection!(Diag,S)                 # advect the PV with U,V
+    
     # Bernoulli potential - recalculate for new η, KEu,KEv are only updated in advection_coriolis
     @unpack p,KEu,KEv,dpdx,dpdy = Diag.Bernoulli
     @unpack g = S.constants