Variable Coriolis (#2201)

AMLattanzi · web-flow · commit 74070a3d0291 · 2025-03-19T08:00:27.000-07:00
* Variable Coriolis.

* Make coriolis vars not depend upon netcdf since we need their pointer in other places.

* clean up style.

* Add docs and make user choose to employ var_cor.
diff --git a/Docs/sphinx_doc/Inputs.rst b/Docs/sphinx_doc/Inputs.rst
@@ -1169,6 +1169,10 @@ List of Parameters
 | **erf.use_coriolis**                | Include Coriolis       | true / false      | false               |
 |                                     | forcing                |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+
+| **erf.variable_coriolis**           | Include Coriolis       | true / false      | false               |
+|                                     | forcing that varies    |                   |                     |
+|                                     | with latitude          |                   |                     |
++-------------------------------------+------------------------+-------------------+---------------------+
 | **erf.rotational_time_period**      | Used to calculate the  | Real              | 86400.0             |
 |                                     | Coriolis frequency     |                   |                     |
 +-------------------------------------+------------------------+-------------------+---------------------+
diff --git a/Docs/sphinx_doc/theory/Forcings.rst b/Docs/sphinx_doc/theory/Forcings.rst
@@ -54,6 +54,14 @@ period (measured in seconds), and :math:`\phi` the latitude.
 Values for ``erf.rotational_time_period``, ``erf.latitude``, and ``erf.coriolis_3d``; the first two are used
 to compute the Coriolis frequency and the last of these determines whether to include the z-component in the Coriolis forcing.
 
+When initializing from a ``wrfinput`` or ``met_em`` file, the latitude at the grid cell centers will be known. For this case, a user may specify
+
+::
+
+      variable_coriolis == true
+
+to use the grid latitude, :math:`\phi(y)`, when computing the sine and cosine coefficients above.
+
 There is no dependence on the radial distance from the center of the earth, thus the curvature of the earth is neglected.
 
 Rayleigh Damping
diff --git a/Source/BoundaryConditions/ERF_ABLMost.H b/Source/BoundaryConditions/ERF_ABLMost.H
@@ -468,7 +468,7 @@ public:
     amrex::Real
     get_zref () { return m_ma.get_zref(); }
 
-    const amrex::FArrayBox*
+    amrex::FArrayBox*
     get_z0 (const int& lev) { return &z_0[lev]; }
 
     bool have_variable_sea_roughness() { return m_var_z0; }
diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H
@@ -326,6 +326,8 @@ struct SolverChoice {
 
         // Include Coriolis forcing?
         pp.query("use_coriolis", use_coriolis);
+        pp.query("has_lat_lon", has_lat_lon);
+        if (has_lat_lon) { pp.query("variable_coriolis", variable_coriolis); }
 
         // Include Rayleigh damping (separate flags for each variable)
         pp.query("rayleigh_damp_U", rayleigh_damp_U);
@@ -369,7 +371,7 @@ struct SolverChoice {
 
         if (use_coriolis)
         {
-            build_coriolis_forcings(pp_prefix);
+            build_coriolis_forcings_const_lat(pp_prefix);
         }
 
         pp.query("add_custom_rhotheta_forcing", custom_rhotheta_forcing);
@@ -591,7 +593,7 @@ struct SolverChoice {
         }
     }
 
-    void build_coriolis_forcings(std::string pp_prefix)
+    void build_coriolis_forcings_const_lat (std::string pp_prefix)
     {
         amrex::ParmParse pp(pp_prefix);
 
@@ -765,6 +767,9 @@ struct SolverChoice {
     std::string abl_geo_wind_table;
     bool have_geo_wind_profile {false};
 
+    bool has_lat_lon{false};
+    bool variable_coriolis{false};
+
     int ave_plane {2};
     // Microphysics params
     bool use_moist_background {false};
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -678,9 +678,11 @@ private:
     amrex::Vector<amrex::Vector<amrex::FArrayBox>> bdy_data_yhi;
 
     amrex::Real bdy_time_interval;
-    amrex::Vector<std::unique_ptr<amrex::MultiFab>> lat_m, lon_m;
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> lat_m, lon_m;       // Latitude and Longitude on the grid
 #endif // ERF_USE_NETCDF
 
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> sinPhi_m, cosPhi_m; // Coriolis factors on the grid
+
     // Struct for working with the sounding data we take as an input
     InputSoundingData input_sounding_data;
 
diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -344,13 +344,20 @@ ERF::ERF_shared ()
     // Size lat long arrays if using netcdf
     lat_m.resize(nlevs_max);
     lon_m.resize(nlevs_max);
-    for (int lev = 0; lev < max_level; ++lev)
-    {
+    for (int lev = 0; lev < max_level; ++lev) {
         lat_m[lev] = nullptr;
         lon_m[lev] = nullptr;
     }
 #endif
 
+    // Variable coriolis
+    sinPhi_m.resize(nlevs_max);
+    cosPhi_m.resize(nlevs_max);
+    for (int lev = 0; lev < max_level; ++lev) {
+        sinPhi_m[lev] = nullptr;
+        cosPhi_m[lev] = nullptr;
+    }
+
     // Initialize tagging criteria for mesh refinement
     refinement_criteria_setup();
 
diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp
@@ -228,7 +228,34 @@ ERF::WriteCheckpointFile () const
             }
             VisMF::Write(z0, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "Z0"));
         }
-    }
+
+#ifdef ERF_USE_NETCDF
+         // Write lat/lon if it exists
+        if (lat_m[lev] && lon_m[lev] && solverChoice.has_lat_lon) {
+            amrex::Print() << "Writing Lat/Lon variables" << std::endl;
+            IntVect ngv = ng; ngv[2] = 0;
+            MultiFab lat(ba2d,dmap[lev],1,ngv);
+            MultiFab lon(ba2d,dmap[lev],1,ngv);
+            MultiFab::Copy(lat,*lat_m[lev],0,0,1,ngv);
+            MultiFab::Copy(lon,*lon_m[lev],0,0,1,ngv);
+            VisMF::Write(lat, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "LAT"));
+            VisMF::Write(lon, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "LON"));
+        }
+
+        // Write sinPhi and cosPhi if it exists
+        if (cosPhi_m[lev] && sinPhi_m[lev] && solverChoice.variable_coriolis) {
+            amrex::Print() << "Writing Coriolis factors" << std::endl;
+            IntVect ngv = ng; ngv[2] = 0;
+            MultiFab sphi(ba2d,dmap[lev],1,ngv);
+            MultiFab cphi(ba2d,dmap[lev],1,ngv);
+            MultiFab::Copy(sphi,*sinPhi_m[lev],0,0,1,ngv);
+            MultiFab::Copy(cphi,*cosPhi_m[lev],0,0,1,ngv);
+            VisMF::Write(sphi, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "SinPhi"));
+            VisMF::Write(cphi, MultiFabFileFullPrefix(lev, checkpointname, "Level_", "CosPhi"));
+        }
+#endif
+
+    } // for lev
 
 #ifdef ERF_USE_PARTICLES
    particleData.Checkpoint(checkpointname);
@@ -532,7 +559,6 @@ ERF::ReadCheckpointFile ()
         }
         BoxArray ba2d(std::move(bl2d));
 
-        {
         IntVect ng = mapfac_m[lev]->nGrowVect();
         MultiFab mf_m(ba2d,dmap[lev],1,ng);
         VisMF::Read(mf_m, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_m"));
@@ -547,8 +573,50 @@ ERF::ReadCheckpointFile ()
         MultiFab mf_v(convert(ba2d,IntVect(0,1,0)),dmap[lev],1,ng);
         VisMF::Read(mf_v, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "MapFactor_v"));
         MultiFab::Copy(*mapfac_v[lev],mf_v,0,0,1,ng);
+
+        if (m_most && m_most->have_variable_sea_roughness())  {
+            amrex::Print() << "Reading variable surface roughness" << std::endl;
+            ng = vars_new[lev][Vars::cons].nGrowVect(); ng[2]=0;
+            MultiFab z0(ba2d,dmap[lev],1,ng);
+            VisMF::Read(z0, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "Z0"));
+            for (amrex::MFIter mfi(z0); mfi.isValid(); ++mfi) {
+                const Box& bx = mfi.growntilebox();
+                FArrayBox* most_z0 = (m_most->get_z0(lev));
+                most_z0->copy<RunOn::Host>(z0[mfi], bx);
+            }
         }
-    }
+
+#ifdef ERF_USE_NETCDF
+        // Read lat/lon if it exists
+        if (solverChoice.has_lat_lon) {
+            amrex::Print() << "Reading Lat/Lon variables" << std::endl;
+            IntVect ngv = ng; ngv[2] = 0;
+            MultiFab lat(ba2d,dmap[lev],1,ngv);
+            MultiFab lon(ba2d,dmap[lev],1,ngv);
+            VisMF::Read(lat, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "LAT"));
+            VisMF::Read(lon, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "LON"));
+            lat_m[lev] = std::make_unique<MultiFab>(ba2d,dmap[lev],1,ngv);
+            lon_m[lev] = std::make_unique<MultiFab>(ba2d,dmap[lev],1,ngv);
+            MultiFab::Copy(*lat_m[lev],lat,0,0,1,ngv);
+            MultiFab::Copy(*lon_m[lev],lon,0,0,1,ngv);
+        }
+
+        // Read sinPhi and cosPhi if it exists
+        if (solverChoice.variable_coriolis) {
+            amrex::Print() << "Reading Coriolis factors" << std::endl;
+            IntVect ngv = ng; ngv[2] = 0;
+            MultiFab sphi(ba2d,dmap[lev],1,ngv);
+            MultiFab cphi(ba2d,dmap[lev],1,ngv);
+            VisMF::Read(sphi, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "SinPhi"));
+            VisMF::Read(cphi, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "CosPhi"));
+            sinPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dmap[lev],1,ngv);
+            cosPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dmap[lev],1,ngv);
+            MultiFab::Copy(*sinPhi_m[lev],sphi,0,0,1,ngv);
+            MultiFab::Copy(*cosPhi_m[lev],cphi,0,0,1,ngv);
+        }
+#endif
+
+    } // for lev
 
 #ifdef ERF_USE_PARTICLES
     restartTracers((ParGDBBase*)GetParGDB(),restart_chkfile);
diff --git a/Source/Initialization/ERF_InitFromMetgrid.cpp b/Source/Initialization/ERF_InitFromMetgrid.cpp
@@ -205,18 +205,27 @@ ERF::init_from_metgrid (int lev)
         }
     }
 
-    lat_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+    solverChoice.has_lat_lon       = true;
+    lat_m[lev]    = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+    sinPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+    cosPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
     for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
         Box gtbx = mfi.growntilebox();
         FArrayBox& dst = (*(lat_m[lev]))[mfi];
         FArrayBox& src = NC_LAT_fab[0];
+        const Array4<      Real>& sin_arr = (sinPhi_m[lev])->array(mfi);
+        const Array4<      Real>& cos_arr = (cosPhi_m[lev])->array(mfi);
         const Array4<      Real>& dst_arr = dst.array();
         const Array4<const Real>& src_arr = src.const_array();
         ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
         {
             int li = min(max(i, i_lo), i_hi);
             int lj = min(max(j, j_lo), j_hi);
             dst_arr(i,j,0) = src_arr(li,lj,0);
+
+            Real lat_rad = dst_arr(i,j,0) * (PI/180.);
+            sin_arr(i,j,0) = std::sin(lat_rad);
+            cos_arr(i,j,0) = std::cos(lat_rad);
         });
     }
 
diff --git a/Source/Initialization/ERF_InitFromWRFInput.cpp b/Source/Initialization/ERF_InitFromWRFInput.cpp
@@ -375,65 +375,74 @@ ERF::init_from_wrfinput (int lev)
           int i_lo = geom[lev].Domain().smallEnd(0); int i_hi = geom[lev].Domain().bigEnd(0);
           int j_lo = geom[lev].Domain().smallEnd(1); int j_hi = geom[lev].Domain().bigEnd(1);
 
-          // Initialize Latitude
+          // Initialize Latitude & Coriolis factors
           if ( var_name == "XLAT_V" ) {
-            lat_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
-            for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-              Box gtbx = mfi.growntilebox();
-              const Array4<      Real>& dst_arr = (lat_m[lev])->array(mfi);
-              const Array4<const Real>& src_arr = var_fab.const_array();
-              ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
-              {
-                int li = amrex::min(amrex::max(i, i_lo), i_hi);
-                int lj = amrex::min(amrex::max(j, j_lo), j_hi);
-                dst_arr(i,j,0) = src_arr(li,lj,0);
-              });
-            }
+              solverChoice.has_lat_lon       = true;
+              lat_m[lev]    = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+              sinPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+              cosPhi_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+              for ( MFIter mfi(*(lat_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+                  Box gtbx = mfi.growntilebox();
+                  const Array4<      Real>& sin_arr = (sinPhi_m[lev])->array(mfi);
+                  const Array4<      Real>& cos_arr = (cosPhi_m[lev])->array(mfi);
+                  const Array4<      Real>& dst_arr = (lat_m[lev])->array(mfi);
+                  const Array4<const Real>& src_arr = var_fab.const_array();
+                  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
+                  {
+                      int li = amrex::min(amrex::max(i, i_lo), i_hi);
+                      int lj = amrex::min(amrex::max(j, j_lo), j_hi);
+                      dst_arr(i,j,0) = src_arr(li,lj,0);
+
+                      Real lat_rad = dst_arr(i,j,0) * (PI/180.);
+                      sin_arr(i,j,0) = std::sin(lat_rad);
+                      cos_arr(i,j,0) = std::cos(lat_rad);
+                  });
+              }
           }
 
           // Initialize Longitude
           if ( var_name == "XLONG_U" ) {
-            lon_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
-            for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-              Box gtbx = mfi.growntilebox();
-              const Array4<      Real>& dst_arr = (lon_m[lev])->array(mfi);
-              const Array4<const Real>& src_arr = var_fab.const_array();
-              ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
-              {
-                int li = amrex::min(amrex::max(i, i_lo), i_hi);
-                int lj = amrex::min(amrex::max(j, j_lo), j_hi);
-                dst_arr(i,j,0) = src_arr(li,lj,0);
-              });
-            }
+              lon_m[lev] = std::make_unique<MultiFab>(ba2d,dm,1,ngv);
+              for ( MFIter mfi(*(lon_m[lev]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+                  Box gtbx = mfi.growntilebox();
+                  const Array4<      Real>& dst_arr = (lon_m[lev])->array(mfi);
+                  const Array4<const Real>& src_arr = var_fab.const_array();
+                  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
+                  {
+                      int li = amrex::min(amrex::max(i, i_lo), i_hi);
+                      int lj = amrex::min(amrex::max(j, j_lo), j_hi);
+                      dst_arr(i,j,0) = src_arr(li,lj,0);
+                  });
+              }
           }
 
           // Initialize Landmask
           if ( var_name == "LANDMASK" ) {
-            for ( MFIter mfi(*(lmask_lev[lev][0]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
-              Box gtbx = mfi.growntilebox();
-              const Array4<       int>& dst_arr = lmask_lev[lev][0]->array(mfi);
-              const Array4<const Real>& src_arr = var_fab.const_array();
-              ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
-              {
-                int li = amrex::min(amrex::max(i, i_lo), i_hi);
-                int lj = amrex::min(amrex::max(j, j_lo), j_hi);
-                dst_arr(i,j,0) = static_cast<int>(src_arr(li,lj,0));
-              });
-            }
-            (lmask_lev[lev])[0]->FillBoundary(geom[lev].periodicity());
+              for ( MFIter mfi(*(lmask_lev[lev][0]), TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+                  Box gtbx = mfi.growntilebox();
+                  const Array4<       int>& dst_arr = lmask_lev[lev][0]->array(mfi);
+                  const Array4<const Real>& src_arr = var_fab.const_array();
+                  ParallelFor(gtbx, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
+                  {
+                      int li = amrex::min(amrex::max(i, i_lo), i_hi);
+                      int lj = amrex::min(amrex::max(j, j_lo), j_hi);
+                      dst_arr(i,j,0) = static_cast<int>(src_arr(li,lj,0));
+                  });
+              }
+              (lmask_lev[lev])[0]->FillBoundary(geom[lev].periodicity());
           }
 
           // Initialize MapFac U
           if ( var_name == "MAPFAC_UY" ) {
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
-            for ( MFIter mfi(*mapfac_u[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
-            {
-              // Define fabs for holding the initial data
-              FArrayBox &msf_fab = (*mapfac_u[lev])[mfi];
-              msf_fab.template copy<RunOn::Device>(var_fab);
-            }
+              for ( MFIter mfi(*mapfac_u[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
+              {
+                  // Define fabs for holding the initial data
+                  FArrayBox &msf_fab = (*mapfac_u[lev])[mfi];
+                  msf_fab.template copy<RunOn::Device>(var_fab);
+              }
           }
 
           // Initialize MapFac V
@@ -454,20 +463,20 @@ ERF::init_from_wrfinput (int lev)
 #ifdef _OPENMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
-            for ( MFIter mfi(*mapfac_m[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
-            {
-              // Define fabs for holding the initial data
-              FArrayBox &msf_fab = (*mapfac_m[lev])[mfi];
-              msf_fab.template copy<RunOn::Device>(var_fab);
-            }
+              for ( MFIter mfi(*mapfac_m[lev], TilingIfNotGPU()); mfi.isValid(); ++mfi )
+              {
+                  // Define fabs for holding the initial data
+                  FArrayBox &msf_fab = (*mapfac_m[lev])[mfi];
+                  msf_fab.template copy<RunOn::Device>(var_fab);
+              }
           }
 
           if (success) {
               var_fab.clear();
               Print() << " DONE\n";
           }
         } // ivar
-      Print() << "\n";
+        Print() << "\n";
     } // idx
 
     // Convert the velocities using the map factors
diff --git a/Source/SourceTerms/ERF_MakeMomSources.cpp b/Source/SourceTerms/ERF_MakeMomSources.cpp
diff --git a/Source/SourceTerms/ERF_SrcHeaders.H b/Source/SourceTerms/ERF_SrcHeaders.H
diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H

Original file line number	Diff line number	Diff line change
`@@ -326,6 +326,8 @@ struct SolverChoice {`
`326`	`326`
`327`	`327`	`// Include Coriolis forcing?`
`328`	`328`	`pp.query("use_coriolis", use_coriolis);`
	`329`	`+ pp.query("has_lat_lon", has_lat_lon);`
	`330`	`+ if (has_lat_lon) { pp.query("variable_coriolis", variable_coriolis); }`
`329`	`331`
`330`	`332`	`// Include Rayleigh damping (separate flags for each variable)`
`331`	`333`	`pp.query("rayleigh_damp_U", rayleigh_damp_U);`
`@@ -369,7 +371,7 @@ struct SolverChoice {`
`369`	`371`
`370`	`372`	`if (use_coriolis)`
`371`	`373`	`{`
`372`		`- build_coriolis_forcings(pp_prefix);`
	`374`	`+ build_coriolis_forcings_const_lat(pp_prefix);`
`373`	`375`	`}`
`374`	`376`
`375`	`377`	`pp.query("add_custom_rhotheta_forcing", custom_rhotheta_forcing);`
`@@ -591,7 +593,7 @@ struct SolverChoice {`
`591`	`593`	`}`
`592`	`594`	`}`
`593`	`595`
`594`		`- void build_coriolis_forcings(std::string pp_prefix)`
	`596`	`+ void build_coriolis_forcings_const_lat (std::string pp_prefix)`
`595`	`597`	`{`
`596`	`598`	`amrex::ParmParse pp(pp_prefix);`
`597`	`599`
`@@ -765,6 +767,9 @@ struct SolverChoice {`
`765`	`767`	`std::string abl_geo_wind_table;`
`766`	`768`	`bool have_geo_wind_profile {false};`
`767`	`769`
	`770`	`+ bool has_lat_lon{false};`
	`771`	`+ bool variable_coriolis{false};`
	`772`	`+`
`768`	`773`	`int ave_plane {2};`
`769`	`774`	`// Microphysics params`
`770`	`775`	`bool use_moist_background {false};`