Merge pull request #116 from milankl/forecasts

Forecasts

Author: Milan K

experiments/forecast_ens.jl

@@ -0,0 +1,49 @@
+include("/home/kloewer/git/Juls.jl/src/Juls.jl")
+using .Juls
+using FileIO
+
+"""Finds the first gap in a list of integers."""
+function gap(a::Array{Int,1})
+    try
+        return minimum([i for i in minimum(a):maximum(a) if ~(i in a)])
+    catch
+        return maximum(a)+1
+    end
+end
+
+function get_run_id(path::String,order::String="continue")
+    runlist = filter(x->startswith(x,"run"),readdir(path))
+    existing_runs = [parse(Int,id[4:end]) for id in runlist]
+    if length(existing_runs) == 0           # if no runfolder exists yet
+        run_id = 0
+    else                                    # create next folder
+        if order == "fill"  # find the smallest gap in runfolders
+            run_id = gap(existing_runs)
+        elseif order == "continue" # find largest folder and count one up
+            run_id = maximum(existing_runs)+1
+        end
+    end
+    return run_id
+end
+
+path = "/network/aopp/chaos/pred/kloewer/julsdata/forecast/"
+startis = load(joinpath(path,"starti.jld2"))["starti"]
+outpath = joinpath(path,"Float64")
+
+for i in 1:10
+    run_id = get_run_id(outpath)
+
+    starti = startis[run_id+1]
+
+    RunJuls(Float64,
+            output=true,
+            Ndays=100.0,
+            outpath=outpath,
+            initial_cond="ncfile",
+            output_vars=["u","v","η","sst"],
+            initpath=path,
+            init_run_id=1,
+            init_starti=starti
+            )
+
+end

experiments/forecast_master.jl

@@ -6,7 +6,7 @@ RunJuls(Float32,
         outpath="/network/aopp/chaos/pred/kloewer/julsdata/forecast",
         output_vars=["u","v","η"],
         α=1.,
-        diffusion="Constant")
+        diffusion="constant")
 
 RunJuls(Float32,
         output=true,
@@ -16,4 +16,4 @@ RunJuls(Float32,
         initial_cond="ncfile",
         init_run_id=0,
         α=1.,
-        diffusion="Constant")
+        diffusion="constant")

src/DefaultParameters.jl

@@ -76,15 +76,15 @@
     sst_initial::String="south"         # "west", "south", "rect", "flat" or "restart"
     sst_rect_coords::Array{Float64,1}=[0.,0.15,0.,1.0]
                                         # (x0,x1,y0,y1) are the size of the rectangle in [0,1]
-    Uadv::Real=0.25                     # Velocity scale [m/s] for tracer advection
+    Uadv::Real=0.2                      # Velocity scale [m/s] for tracer advection
     SSTmax::Real=1.                     # tracer (sea surface temperature) max for restoring
     SSTmin::Real=0.                     # tracer (sea surface temperature) min for restoring
     τSST::Real=500.                     # tracer restoring time scale [days]
     jSST::Real=365.                     # tracer consumption [days]
     SST_λ0::Real=222e3                  # [m] transition position of relaxation timescale
     SST_λs::Real=111e3                  # [m] transition width of relaxation timescale
     SST_γ0::Real=8.35                   # [days] injection time scale
-    SSTw::Real=50e3                     # width [m] of the tangent used for the IC and interface relaxation
+    SSTw::Real=5e3                      # width [m] of the tangent used for the IC and interface relaxation
     SSTϕ::Real=0.5                      # latitude/longitude fraction ∈ [0,1] of sst edge
 
     # OUTPUT OPTIONS
@@ -111,7 +111,7 @@
     @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, $Ndyas 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."
     @assert nstep_advcor >= 0   "nstep_advcor has to be >=0, $nstep_advcor given."
     @assert bc in ["periodic","nonperiodic"]    "boundary condition '$bc' unsupported."

src/TracerAdvection.jl

@@ -148,7 +148,7 @@ function interp_uv!(uvi::AbstractMatrix,
     elseif (m+4,n+2) == (muv,nuv)    # v case
         # halo sizes 1: left/bottom, 2: right/top
         xh1 = 2
-        xh2 = 1         # don't fully understand why this is 1 ...
+        xh2 = 2
         yh1 = 1
         yh2 = 1
     else