Generalize geostrophic forcing (erf-model#1718)

* Add coriolis_3d flag (default: true) * Remove hard-coded coriolis factor when using custom geo wind profiles
ewquon · Aug 1, 2024 · 7b812de · 7b812de
1 parent 71431cc
commit 7b812de
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 9 deletions.
diff --git a/Source/DataStructs/DataStruct.H b/Source/DataStructs/DataStruct.H
@@ -430,10 +430,14 @@ struct SolverChoice {
         amrex::Real latitude = 90.0;
         pp.query("latitude", latitude);
 
+        pp.query("coriolis_3d", coriolis_3d);
+
         // Convert to radians
         latitude *= (PI/180.);
         sinphi = std::sin(latitude);
-        cosphi = std::cos(latitude);
+        if (coriolis_3d) {
+            cosphi = std::cos(latitude);
+        }
 
         amrex::Print() << "Coriolis frequency, f = " << coriolis_factor * sinphi << " 1/s" << std::endl;
 
@@ -478,6 +482,7 @@ struct SolverChoice {
     // Specify what additional physics/forcing modules we use
     bool        use_gravity            = false;
     bool        use_coriolis           = false;
+    bool        coriolis_3d            = true;
 
     bool        rayleigh_damp_U        = false;
     bool        rayleigh_damp_V        = false;
@@ -503,7 +508,7 @@ struct SolverChoice {
 
     // Coriolis forcing
     amrex::Real coriolis_factor = 0.0;
-    amrex::Real cosphi          = 0.0 ;
+    amrex::Real cosphi          = 0.0;
     amrex::Real sinphi          = 0.0;
 
     // User-specified forcings in problem definition

diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp
@@ -93,7 +93,7 @@ void make_mom_sources (int /*level*/,
     // Flag for Geostrophic forcing
     // *****************************************************************************
     auto abl_geo_forcing  = solverChoice.abl_geo_forcing;
-    auto geostrophic_wind = solverChoice.custom_geostrophic_profile;
+    auto geo_wind_profile = solverChoice.custom_geostrophic_profile;
 
     // *****************************************************************************
     // Data for Rayleigh damping
@@ -286,20 +286,18 @@ void make_mom_sources (int /*level*/,
         // *****************************************************************************
         // Add height-dependent GEOSTROPHIC forcing
         // *****************************************************************************
-        if (geostrophic_wind) {
+        if (geo_wind_profile) {
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) );
-                Real rho_v_loc = 0.25 * (rho_v(i,j+1,k) + rho_v(i,j,k) + rho_v(i-1,j+1,k) + rho_v(i-1,j,k));
-                xmom_src_arr(i, j, k) += -0.376E-4 * (rho_on_u_face * dptr_v_geos[k] - rho_v_loc);
+                xmom_src_arr(i, j, k) -= coriolis_factor * rho_on_u_face * dptr_v_geos[k] * sinphi;
             });
             ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) );
-                Real rho_u_loc = 0.25 * (rho_u(i+1,j,k) + rho_u(i,j,k) + rho_u(i+1,j-1,k) + rho_u(i,j-1,k));
-                ymom_src_arr(i, j, k) += 0.376E-4 * (rho_on_v_face * dptr_u_geos[k] - rho_u_loc);
+                ymom_src_arr(i, j, k) += coriolis_factor * rho_on_v_face * dptr_u_geos[k] * sinphi;
             });
-        } // geostrophic_wind
+        } // geo_wind_profile
 
         // *****************************************************************************
         // Add custom SUBSIDENCE terms