Merge pull request #55 from milankl/ensmode

Fixing wrong scaling for bottom friction and biharmonic viscosity

ens_juls.jl

@@ -7,7 +7,7 @@ using NetCDF
 using FileIO
 
 #using MPI
-#using SigmoidNumbers
+using SigmoidNumbers
 
 # Finite16nonu
 #include("/home/kloewer/julia/FiniteFloats.jl/src/FiniteFloats.jl")
@@ -36,7 +36,7 @@ include("src/feedback.jl")
 include("src/output.jl")
 
 # INITIALISE
-for ens_mem in 1:80
+for ens_mem in 1:10
 	global run_id,runpath = get_run_id_path("fill")
 	starti = load(initpath*"starti.jld2")["starti"][run_id+1]
 	println("Ensemble member $ens_mem, start from $starti")

parameters.jl

@@ -1,6 +1,6 @@
 # define constants
-const Numtype = Float64
-#const Numtype = Posit{16,1}
+const Numtype = Float32
+#const Numtype = Posit{16,0}
 #const Numtype = Main.FiniteFloats.Finite16
 
 const nx = 100                  # number of grid cells in x-direction
@@ -24,27 +24,27 @@ const η_refh = 5.               # height difference [m] of the interface relaxa
 const η_refw = 50e3             # width [m] of the tangent used for the interface relaxation
 
 const cfl = 1.0                 # CFL number
-const Ndays = 100               # number of days to integrate for
+const Ndays = 2000               # number of days to integrate for
 
 # boundary condtions
 const bc_x = "periodic"         # "periodic" or anything else for nonperiodic
 const lbc = 1.                  # lateral boundary condition parameter
                                 # 0 free-slip, 0<lbc<2 partial-slip, 2 no-slip
 
-const adv_scheme = "Sadourny"   # "Sadourny" or "ArakawaHsu"
+const adv_scheme = "ArakawaHsu"   # "Sadourny" or "ArakawaHsu"
 
 const bottom_friction = "linear" # "linear" or "quadratic"
 const drag = 1e-5               # bottom drag coefficient [dimensionless] for quadratic
-const τdrag = 10.               # bottom drag coefficient [days] for linear
+const τdrag = 300.               # bottom drag coefficient [days] for linear
 
 const diffusion = "Constant"    # "Smagorinsky" or "Constant", biharmonic in both cases
-const ν_const = 540             # [m^2/s] scaling constant for Constant biharmonic diffusion
+const ν_const = 500             # [m^2/s] scaling constant for Constant biharmonic diffusion
 const c_smag = 0.15             # Smagorinsky coefficient [dimensionless]
 
 const output = 1                # 1 for nc output 0 for none
-const output_vars = ["eta"]
+const output_vars = ["u"]
 const output_dt = 6             # output time step in hours
-const outpath = "/network/aopp/cirrus/pred/kloewer/julsdata/forecast/Float64Sad/"
+const outpath = "/network/aopp/cirrus/pred/kloewer/julsdata/stabilitytests/"
 
 const initial_cond = "rest"   # "rest" or "ncfile"
 const initpath = "/network/aopp/cirrus/pred/kloewer/julsdata/forecast/"

run_juls.jl

@@ -4,7 +4,7 @@ etc, and runs the model =#
 using Dates
 using Printf
 using NetCDF
-#using FileIO
+using FileIO
 
 #using MPI
 #using SigmoidNumbers
@@ -34,10 +34,8 @@ include("src/preallocate.jl")
 # OUTPUT AND FEEDBACK
 include("src/feedback.jl")
 include("src/output.jl")
-global run_id,runpath = get_run_id_path("fill")
+global run_id,runpath = get_run_id_path()
 
-# INITIALISE
+# INITIALISE & RUN
 u,v,η = initial_conditions()
-
-# RUN
 u,v,η = time_integration(u,v,η)

src/constants.jl

@@ -15,16 +15,16 @@ const one_quarter = Numtype(0.25)
 # will be used for the Bernoulli potential
 const g = Numtype(gravity)
 
-# for the bottom friction
-const c_D = Numtype(drag)
-const r_B = Numtype(1/(τdrag*24*3600))  # [1/s]
+# for the bottom friction (include grid spacing as gradient operators are dimensionless)
+const c_D = Numtype(Δ*drag)
+const r_B = Numtype(Δ/(τdrag*24*3600))  # [1/s]
 
 # frequency [1/s] of interface relaxation (for non-dimensional gradients, γ contains the grid spacing Δ)
 const γ = Numtype(Δ/(t_relax*3600*24))
 
 # for biharmonic diffusion
 const cSmag = Numtype(-c_smag)
-const νB = Numtype(ν_const/(30000*Δ))   # linear scaling based on 540m^s/s at Δ=30km
+const νB = Numtype(-ν_const/30000)   # linear scaling based on 540m^s/s at Δ=30km
 
 # for analysis the old operators with boundary conditions are kept. Constants for these
 const one_over_Δ = Numtype(1/Δ)