From 093cecfe277b09122ee87b3c77a5e05a524a5092 Mon Sep 17 00:00:00 2001
From: milankl <milankloewer@gmx.de>
Date: Tue, 5 Feb 2019 12:42:39 +0000
Subject: [PATCH] Wind forcing Fx,Fy of both u,v

---
 parameters.jl           | 34 ++++++++++++++++++----------------
 src/forcing.jl          | 38 ++++++++++++++++++++++----------------
 src/output.jl           |  4 +++-
 src/rhs.jl              | 18 +++++++++---------
 src/time_integration.jl |  5 +++--
 5 files changed, 55 insertions(+), 44 deletions(-)

diff --git a/parameters.jl b/parameters.jl
index 6984025..6da6105 100644
--- a/parameters.jl
+++ b/parameters.jl
@@ -6,25 +6,27 @@ const Numtype = Float32
 #setprecision(7)
 
 # DOMAIN RESOLUTION AND RATIO
-const nx = 400                      # number of grid cells in x-direction
-const Lx = 10000e3                   # length of the domain in x-direction
-const L_ratio = 2                   # Domain aspect ratio of Lx/Ly
+const nx = 200                      # number of grid cells in x-direction
+const Lx = 6000e3                   # length of the domain in x-direction
+const L_ratio = 1                   # Domain aspect ratio of Lx/Ly
 
 # PHYSICAL CONSTANTS
 const gravity = 0.032                # gravitational acceleration
 const water_depth = 125.            # layer thickness at rest
 const ρ = 1e3                       # density
-const ϕ = 10.                       # central latitue of the domain (for coriolis)
+const ϕ = 0.                       # central latitue of the domain (for coriolis)
 const ω = 2π/(24*3600)              # Earth's angular frequency [s^-1]
 const R = 6.371e6                   # Earth's radius [m]
 
 # WIND FORCING OPTIONS
-const wind_forcing = "constant"      # "channel", "double_gyre", "shear" or "none"
-const Fx0 = 0.01                    # wind stress strength [Pa], default 0.12
+const wind_forcing_x = "constant"      # "channel", "double_gyre", "shear","constant" or "none"
+const wind_forcing_y = "constant"      # "channel", "double_gyre", "shear","constant" or "none"
+const Fx0 = 0.0                          # wind stress strength [Pa], default 0.12
+const Fy0 = 0.01
 
 # BOTTOM TOPOGRAPHY OPTIONS
 const topography_feature = "flat"  # "ridge", "seamount", "flat", "ridges", "bathtub"
-const topofeat_height = 50.         # height of seamount
+const topofeat_height = 10.         # height of seamount
 const topofeat_width = 300e3        # horizontal scale [m] of the seamount
 
 # NEWTONIAN COOLING OPTIONS
@@ -35,14 +37,14 @@ const η_refw = 100e3                # width [m] of the tangent used for the int
 
 # TIME STEPPING OPTIONS
 const RKo = 4                       # Order of the RK time stepping scheme (3 or 4)
-const cfl = 0.9                     # CFL number (1.0 recommended for RK4, 0.6 for RK3)
-const Ndays = 600                   # number of days to integrate for
+const cfl = 1.0                     # CFL number (1.0 recommended for RK4, 0.6 for RK3)
+const Ndays = 1600                   # number of days to integrate for
 const nstep_diff = 1                # diffusive part every nstep_diff time steps.
 const nstep_advcor = 1              # advection and coriolis every nstep_advcor time steps.
 
 # BOUNDARY CONDITION OPTIONS
 const bc_x = "nonperiodic"             # "periodic" or anything else for nonperiodic
-const lbc = 1.                      # lateral boundary condition parameter
+const lbc = 0.                      # lateral boundary condition parameter
                                     # 0 free-slip, 0<lbc<2 partial-slip, 2 no-slip
 
 # MOMENTUM ADVECTION OPTIONS
@@ -56,7 +58,7 @@ const τdrag = 300.                  # bottom drag coefficient [days] for linear
 
 # DIFFUSION OPTIONS
 const diffusion = "Constant"        # "Smagorinsky" or "Constant", biharmonic in both cases
-const ν_const = 1000.0               # [m^2/s] scaling constant for Constant biharmonic diffusion
+const ν_const = 500.0               # [m^2/s] scaling constant for Constant biharmonic diffusion
 const c_smag = 0.15                 # Smagorinsky coefficient [dimensionless]
 
 # TRACER ADVECTION
@@ -75,15 +77,15 @@ const SSTϕ = 0.5                    # latitude/longitude ∈ [0,1] of sst edge
 const output = true                 # for nc output
 const output_tend = false            # ouput for tendencies as well?
 const output_diagn = false           # output diagnostic variables as well?
-const output_progn_vars = ["u","v","eta"]
+const output_progn_vars = ["eta"]
 const output_tend_vars = ["du","dv","deta","Bu","Bv","LLu1","LLu2","LLv1","LLv2",
                             "qhv","qhu","dpdx","dpdy","dUdx","dVdy"]
 #const output_tend_vars = ["qhv","qhu","dpdx","dpdy","dUdx","dVdy"]
-const output_diagn_vars = ["q","relvort"]#,"q","p","dudx","dvdy","dudy","dvdx","Lu","Lv","xd","yd"]
+const output_diagn_vars = ["q"]#,"q","p","dudx","dvdy","dudy","dvdx","Lu","Lv","xd","yd"]
 
-const output_dt = 24                 # output time step in hours
-#const outpath = "/network/aopp/chaos/pred/kloewer/julsdata/ssthr/"
-const outpath = "/Users/milan/phd/"
+const output_dt = 72                 # output time step in hours
+const outpath = "/network/aopp/chaos/pred/kloewer/julsdata/wave/"
+#const outpath = "/Users/milan/phd/"
 
 # INITIAL CONDITIONS
 const initial_cond = "rest"       # "rest" or "ncfile"
diff --git a/src/forcing.jl b/src/forcing.jl
index 965c9d1..684b1ff 100644
--- a/src/forcing.jl
+++ b/src/forcing.jl
@@ -25,6 +25,14 @@ function constant_wind()
     return Numtype.(Fx)
 end
 
+"""Returns the constant in space forcing matrix Fy."""
+function constant_wind_y()
+    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
+    xx_v,yy_v = meshgrid(x_v,y_v)
+    Fy = (Δ*Fy0/ρ/water_depth)*ones(size(xx_v))
+    return Numtype.(Fy)
+end
+
 """Returns the constant forcing matrix Fx that varies only meriodionally
 with a superposition of sin & cos for a double gyre circulation.
 See Cooper&Zanna 2015 or Kloewer et al 2018."""
@@ -35,12 +43,6 @@ function double_gyre_wind()
     return Numtype.(Fx)
 end
 
-"""Returns a zero matrix for no wind forcing."""
-function no_wind()
-    return zeros(Numtype,nux,nuy)
-end
-
-
 """Returns a reference state for Newtonian cooling/surface relaxation shaped as a
 hyperbolic tangent to force the continuity equation."""
 function interface_relaxation()
@@ -50,16 +52,20 @@ function interface_relaxation()
     return Numtype.(η_ref)
 end
 
-if wind_forcing == "channel"
-    wind = channel_wind
-elseif wind_forcing == "shear"
-    wind = shear_wind
-elseif wind_forcing == "double_gyre"
-    wind = double_gyre_wind
-elseif wind_forcing == "constant"
-    wind = constant_wind
-elseif wind_forcing == "none"
-    wind = no_wind
+if wind_forcing_x == "channel"
+    windx = channel_wind
+elseif wind_forcing_x == "shear"
+    windx = shear_wind
+elseif wind_forcing_x == "double_gyre"
+    windx = double_gyre_wind
+elseif wind_forcing_x == "constant"
+    windx = constant_wind
+else
+    throw(error("Wind forcing not correctly specified."))
+end
+
+if wind_forcing_y == "constant"
+    windy = constant_wind_y
 else
     throw(error("Wind forcing not correctly specified."))
 end
diff --git a/src/output.jl b/src/output.jl
index 527b1da..f0525ce 100644
--- a/src/output.jl
+++ b/src/output.jl
@@ -378,8 +378,10 @@ function output_dict()
     Dictu["phi"] = ϕ
     Dictu["density"] = ρ
 
-    Dictu["wind_forcing"] = wind_forcing
+    Dictu["wind_forcing_x"] = wind_forcing_x
+    Dictu["wind_forcing_y"] = wind_forcing_y
     Dictu["Fx0"] = Fx0
+    Dictu["Fy0"] = Fy0
 
     Dictu["topography_feature"] = topography_feature
     Dictu["topofeat_height"] = topofeat_height
diff --git a/src/rhs.jl b/src/rhs.jl
index 3508f7e..0f62c2f 100644
--- a/src/rhs.jl
+++ b/src/rhs.jl
@@ -1,7 +1,7 @@
 """Tendencies du,dv,dη of
 
         ∂u/∂t = qhv - ∂(1/2*(u²+v²) + gη)/∂x + Fx
-        ∂v/∂t = -qhu - ∂(1/2*(u²+v²) + gη)/∂y
+        ∂v/∂t = -qhu - ∂(1/2*(u²+v²) + gη)/∂y + Fy
         ∂η/∂t =  -∂(uh)/∂x - ∂(vh)/∂y + γ(η_ref-η)
 
 where non-linear terms
@@ -9,7 +9,7 @@ where non-linear terms
         q, (u²+v²)
 
 of the mom eq. are only updated outside the function (slowly varying terms)."""
-function rhs_nonlin!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
+function rhs_nonlin!(du,dv,dη,u,v,η,Fx,Fy,f_u,f_v,f_q,H,η_ref,
             dvdx,dudy,dpdx,dpdy,
             p,u²,v²,KEu,KEv,dUdx,dVdy,
             h,h_u,h_v,h_q,U,V,U_v,V_u,u_v,v_u,
@@ -43,7 +43,7 @@ function rhs_nonlin!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
 
     # adding the terms
     momentum_u!(du,qhv,dpdx,Fx)
-    momentum_v!(dv,qhu,dpdy)
+    momentum_v!(dv,qhu,dpdy,Fy)
     continuity!(dη,dUdx,dVdy,η,η_ref)
 end
 
@@ -54,7 +54,7 @@ end
         ∂η/∂t =  -∂(uH)/∂x - ∂(vH)/∂y + γ(η_ref-η),
 
 the linear shallow water equations."""
-function rhs_lin!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
+function rhs_lin!(du,dv,dη,u,v,η,Fx,Fy,f_u,f_v,f_q,H,η_ref,
             dvdx,dudy,dpdx,dpdy,
             p,u²,v²,KEu,KEv,dUdx,dVdy,
             h,h_u,h_v,h_q,U,V,U_v,V_u,u_v,v_u,
@@ -81,7 +81,7 @@ function rhs_lin!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
 
     # adding the terms
     momentum_u!(du,qhv,dpdx,Fx)
-    momentum_v!(dv,qhu,dpdy)
+    momentum_v!(dv,qhu,dpdy,Fy)
     continuity!(dη,dUdx,dVdy,η,η_ref)
 end
 
@@ -92,7 +92,7 @@ end
         ∂η/∂t =  -∂(uh)/∂x - ∂(vh)/∂y + γ(η_ref-η),
 
 the shallow water equations which are linear in momentum but non-linear in continuity."""
-function rhs_linmom!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
+function rhs_linmom!(du,dv,dη,u,v,η,Fx,Fy,f_u,f_v,f_q,H,η_ref,
             dvdx,dudy,dpdx,dpdy,
             p,KEu,KEv,dUdx,dVdy,
             h,h_u,h_v,h_q,U,V,U_v,V_u,u_v,v_u,
@@ -124,7 +124,7 @@ function rhs_linmom!(du,dv,dη,u,v,η,Fx,f_u,f_v,f_q,H,η_ref,
 
     # adding the terms
     momentum_u!(du,qhv,dpdx,Fx)
-    momentum_v!(dv,qhu,dpdy)
+    momentum_v!(dv,qhu,dpdy,Fy)
     continuity!(dη,dUdx,dVdy,η,η_ref)
 end
 
@@ -250,14 +250,14 @@ function momentum_u!(du,qhv,dpdx,Fx)
 end
 
 """Sum up the tendencies of the non-diffusive right-hand side for the v-component."""
-function momentum_v!(dv,qhu,dpdy)
+function momentum_v!(dv,qhu,dpdy,Fy)
     m,n = size(dv) .- (2*halo,2*halo) # cut off the halo
     @boundscheck (m,n) == size(qhu) || throw(BoundsError())
     @boundscheck (m+2,n+2) == size(dpdy) || throw(BoundsError())
 
     @inbounds for j ∈ 1:n
         for i ∈ 1:m
-             dv[i+2,j+2] = -qhu[i,j] - dpdy[i+1,j+1]
+             dv[i+2,j+2] = -qhu[i,j] - dpdy[i+1,j+1] + Fy[i,j]
         end
     end
 end
diff --git a/src/time_integration.jl b/src/time_integration.jl
index 493079d..80a9001 100644
--- a/src/time_integration.jl
+++ b/src/time_integration.jl
@@ -1,7 +1,8 @@
 function time_integration(u,v,η,sst)
 
     # FORCING
-    Fx = wind()
+    Fx = windx()
+    Fy = windy()
     f_u,f_v,f_q = beta_plane()
     H = topography()
     η_ref = interface_relaxation()
@@ -49,7 +50,7 @@ function time_integration(u,v,η,sst)
                 ghost_points!(u1,v1,η1)
             end
 
-            rhs!(du,dv,dη,u1,v1,η1,Fx,f_u,f_v,f_q,H,η_ref,
+            rhs!(du,dv,dη,u1,v1,η1,Fx,Fy,f_u,f_v,f_q,H,η_ref,
                 dvdx,dudy,dpdx,dpdy,
                 p,u²,v²,KEu,KEv,dUdx,dVdy,
                 h,h_u,h_v,h_q,U,V,U_v,V_u,u_v,v_u,