From add9e240ef00626119ea0b722662a8786ec44594 Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Fri, 10 May 2024 17:29:05 -0500
Subject: [PATCH 1/5] Add framework routines to interp. 2d vector fields

... from cells to edges
---
 .../src/operators/mpas_vector_operations.F    | 154 ++++++++++++++++++
 1 file changed, 154 insertions(+)

diff --git a/components/mpas-framework/src/operators/mpas_vector_operations.F b/components/mpas-framework/src/operators/mpas_vector_operations.F
index 64c848944183..fd7ad78e3470 100644
--- a/components/mpas-framework/src/operators/mpas_vector_operations.F
+++ b/components/mpas-framework/src/operators/mpas_vector_operations.F
@@ -46,6 +46,8 @@ module mpas_vector_operations
              mpas_cross_product_in_r3, &
              mpas_tangential_velocity, &
              mpas_tangential_vector_1d, &
+             mpas_vector_cell_to_edge_anisotropic, &
+             mpas_vector_cell_to_edge_isotropic, &
              mpas_vector_R3Cell_to_2DEdge, &
              mpas_vector_R3Cell_to_normalVectorEdge, &
              mpas_vector_R3_to_LonLatR, &
@@ -288,6 +290,158 @@ subroutine mpas_vector_R3Cell_to_normalVectorEdge(vectorR3Cell, &
 
    end subroutine mpas_vector_R3Cell_to_normalVectorEdge!}}}
 
+
+
+!***********************************************************************
+!
+!  routine mpas_vector_cell_to_edge_anisotropic
+!
+!> \brief   Interpolate a 2d vector field from cell centers to edges
+!> \author  Xylar Asay-Davis
+!> \date    May 2024
+!> \details
+!>  This subroutine interpolates a 2D vector field on the sphere to edges
+!>  using an anisotropic approach that only uses the 2 closest cell centers
+!
+!-----------------------------------------------------------------------
+
+   subroutine mpas_vector_cell_to_edge_anisotropic(zonalCell, meridionalCell, &
+      nEdges, cellsOnEdge, zonalEdge, meridionalEdge)!{{{
+
+      !-----------------------------------------------------------------
+      !
+      ! input variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(in) :: &
+         zonalCell, &    !< Input: zonal component of the vector at cell centers
+         meridionalCell  !< Input: meridional component of the vector at cell centers
+
+      integer, intent(in) :: &
+         nEdges !< Input: The number of edges (level of halo) to include in the computation
+
+      integer, dimension(:,:), intent(in) :: &
+         cellsOnEdge  !< Input: cells adjacent to each mesh edge
+
+      !-----------------------------------------------------------------
+      !
+      ! output variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(out) :: &
+         zonalEdge, &    !< Input: zonal component of the vector at cell centers
+         meridionalEdge  !< Input: meridional component of the vector at cell centers
+
+      !-----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !-----------------------------------------------------------------
+
+      integer :: iEdge, cell1, cell2
+
+      do iEdge = 1, nEdges
+        cell1 = cellsOnEdge(1, iEdge)
+        cell2 = cellsOnEdge(2, iEdge)
+
+        zonalEdge(iEdge) = 0.5_RKIND * (zonalCell(cell1) + &
+                                        zonalCell(cell2))
+        meridionalEdge(iEdge) = 0.5_RKIND * (meridionalCell(cell1) + &
+                                             meridionalCell(cell2))
+      end do
+
+   end subroutine mpas_vector_cell_to_edge_anisotropic!}}}
+
+
+
+!***********************************************************************
+!
+!  routine mpas_vector_cell_to_edge_isotropic
+!
+!> \brief   Interpolate a 2d vector field from cell centers to edges
+!> \author  Xylar Asay-Davis
+!> \date    May 2024
+!> \details
+!>  This subroutine interpolates a 2D vector field on the sphere to edges
+!>  using an isotropic approach that uses 4 adjacent cell centers
+!
+!-----------------------------------------------------------------------
+
+   subroutine mpas_vector_cell_to_edge_isotropic(zonalCell, meridionalCell, &
+      nEdges, vertexDegree, verticesOnEdge, cellsOnVertex, kiteAreasOnVertex, &
+      zonalEdge, meridionalEdge)!{{{
+
+      !-----------------------------------------------------------------
+      !
+      ! input variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(in) :: &
+         zonalCell, &    !< Input: zonal component of the vector at cell centers
+         meridionalCell  !< Input: meridional component of the vector at cell centers
+
+      integer, intent(in) :: &
+         nEdges, &    !< Input: The number of edges (level of halo) to include in the computation
+         vertexDegree !< Input: The number of edges and cells adjacent to each vertex
+
+      integer, dimension(:,:), intent(in) :: &
+         verticesOnEdge, &  !< Input: vertices adjacent to each mesh edge
+         cellsOnVertex      !< Input: cells adjacent to each mesh vertex
+
+      real (kind=RKIND), dimension(:, :), intent(in) :: &
+         kiteAreasOnVertex  !< Input: The area of elements bounded by a cell, a vertex and
+                            !<        the centers of the two edges they share
+
+      !-----------------------------------------------------------------
+      !
+      ! output variables
+      !
+      !-----------------------------------------------------------------
+
+      real (kind=RKIND), dimension(:), intent(out) :: &
+         zonalEdge, &    !< Input: zonal component of the vector at cell centers
+         meridionalEdge  !< Input: meridional component of the vector at cell centers
+
+      !-----------------------------------------------------------------
+      !
+      ! local variables
+      !
+      !-----------------------------------------------------------------
+
+      integer :: iEdge, iVert, iCell, v, c
+
+      real (kind=RKIND) :: area, areaSum, zonal, meridional
+
+      do iEdge = 1, nEdges
+         zonalEdge(iEdge) = 0.0_RKIND
+         meridionalEdge(iEdge) = 0.0_RKIND
+         areaSum = 0.0_RKIND
+         do v = 1, 2
+            iVert = verticesOnEdge(v, iEdge)
+            do c = 1, vertexDegree
+               iCell = cellsOnVertex(c, iVert)
+               ! kite areas are zero for invalid cells on vertex
+               area = kiteAreasOnVertex(c, iVert)
+
+               zonalEdge(iEdge) = zonalEdge(iEdge) &
+                                + zonalCell(iCell) * area
+               meridionalEdge(iEdge) = meridionalEdge(iEdge) &
+                                     + meridionalCell(iCell) * area
+
+               areaSum = areaSum + area
+            end do
+         end do
+         zonalEdge(iEdge) = zonalEdge(iEdge) / areaSum
+         meridionalEdge(iEdge) = meridionalEdge(iEdge) / areaSum
+      end do
+
+   end subroutine mpas_vector_cell_to_edge_isotropic!}}}
+
+
+
 !***********************************************************************
 !
 !  routine mpas_tangential_velocity

From 1bbdea9a8dd5665e0680ab56a6f9f8f44cb2ec9a Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Fri, 10 May 2024 18:28:56 -0500
Subject: [PATCH 2/5] Add OpenACC and OpenMP directives

---
 .../src/operators/mpas_vector_operations.F    | 29 +++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/components/mpas-framework/src/operators/mpas_vector_operations.F b/components/mpas-framework/src/operators/mpas_vector_operations.F
index fd7ad78e3470..108b28786e29 100644
--- a/components/mpas-framework/src/operators/mpas_vector_operations.F
+++ b/components/mpas-framework/src/operators/mpas_vector_operations.F
@@ -342,6 +342,16 @@ subroutine mpas_vector_cell_to_edge_anisotropic(zonalCell, meridionalCell, &
 
       integer :: iEdge, cell1, cell2
 
+#ifdef MPAS_OPENACC
+      !$acc parallel loop &
+      !$acc    present(cellsOnEdge, zonalCell, meridionalCell, &
+      !$acc            zonalEdge, meridionalEdge) &
+      !$acc    private(cell1, cell2)
+#else
+      !$omp parallel
+      !$omp do schedule(runtime) &
+      !$omp    private(cell1, cell2)
+#endif
       do iEdge = 1, nEdges
         cell1 = cellsOnEdge(1, iEdge)
         cell2 = cellsOnEdge(2, iEdge)
@@ -351,6 +361,10 @@ subroutine mpas_vector_cell_to_edge_anisotropic(zonalCell, meridionalCell, &
         meridionalEdge(iEdge) = 0.5_RKIND * (meridionalCell(cell1) + &
                                              meridionalCell(cell2))
       end do
+#ifndef MPAS_OPENACC
+      !$omp end do
+      !$omp end parallel
+#endif
 
    end subroutine mpas_vector_cell_to_edge_anisotropic!}}}
 
@@ -415,6 +429,17 @@ subroutine mpas_vector_cell_to_edge_isotropic(zonalCell, meridionalCell, &
 
       real (kind=RKIND) :: area, areaSum, zonal, meridional
 
+#ifdef MPAS_OPENACC
+      !$acc parallel loop gang vector &
+      !$acc    present(verticesOnEdge, cellsOnVertex, kiteAreasOnVertex, &
+      !$acc            zonalCell, meridionalCell, &
+      !$acc            zonalEdge, meridionalEdge) &
+      !$acc    private(areaSum, area, v, c, iVert, iCell)
+#else
+      !$omp parallel
+      !$omp do schedule(runtime) &
+      !$omp    private(areaSum, area, v, c, iVert, iCell)
+#endif
       do iEdge = 1, nEdges
          zonalEdge(iEdge) = 0.0_RKIND
          meridionalEdge(iEdge) = 0.0_RKIND
@@ -437,6 +462,10 @@ subroutine mpas_vector_cell_to_edge_isotropic(zonalCell, meridionalCell, &
          zonalEdge(iEdge) = zonalEdge(iEdge) / areaSum
          meridionalEdge(iEdge) = meridionalEdge(iEdge) / areaSum
       end do
+#ifndef MPAS_OPENACC
+      !$omp end do
+      !$omp end parallel
+#endif
 
    end subroutine mpas_vector_cell_to_edge_isotropic!}}}
 

From b7924064c8e2bccc44654fbdd50551fe564bdf3d Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Fri, 10 May 2024 21:43:13 -0500
Subject: [PATCH 3/5] Add option to use isotropic interp for surf. wind stress

---
 components/mpas-ocean/src/Registry.xml        |  4 ++
 .../shared/mpas_ocn_surface_bulk_forcing.F    | 60 ++++++++++++-------
 2 files changed, 41 insertions(+), 23 deletions(-)

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index f8b56c083d44..d198ea033bb0 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -679,6 +679,10 @@
 					description="Controls if zonal and meridional components of windstress are used to build surface wind stress."
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_bulk_wind_stress_interp_isotropic" type="logical" default_value=".false."
+					description="Controls if windstress is interpolated to edges from 4 adjacent cells rather than just 2."
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_use_bulk_thickness_flux" type="logical" default_value=".false."
 					description="Controls if a bulk thickness flux will be computed for surface forcing."
 					possible_values=".true. or .false."
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
index e28c8f1e5819..49db5593aa38 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
@@ -25,6 +25,7 @@ module ocn_surface_bulk_forcing
    use mpas_derived_types
    use mpas_pool_routines
    use mpas_timekeeping
+   use mpas_vector_operations
    use ocn_constants
    use ocn_config
    use ocn_mesh
@@ -172,16 +173,13 @@ subroutine ocn_surface_bulk_forcing_vel( &
 
       integer ::         &
          iEdge, iCell,   &! loop indices for edge, cell loops
-         nCells, nEdges, &! number of cells, edges participating
-         cell1, cell2     ! neighbor cell indices across edge
-
-      real (kind=RKIND) :: &
-         meridWSEdge,      &! meridional wind stress avg to edge
-         zonalWSEdge        ! zonal wind stress averaged to edge
+         nCells, nEdges   ! number of cells, edges participating
 
       real (kind=RKIND), dimension(:), pointer :: &
          windStressZonal, &! zonal wind stress from forcing
-         windStressMerid   ! meridional wind stress from forcing
+         windStressMerid, &! meridional wind stress from forcing
+         meridWSEdge,     &! meridional wind stress avg to edge
+         zonalWSEdge       ! zonal wind stress averaged to edge
 
       ! End preamble
       !-----------------------------------------------------------------
@@ -211,28 +209,38 @@ subroutine ocn_surface_bulk_forcing_vel( &
       !*** ocean model wind stress by averaging cell-centered wind
       !*** stress to edges and then rotating vectors
 
+      allocate(meridWSEdge(nEdges), &
+               zonalWSEdge(nEdges))
+
+#ifdef MPAS_OPENACC
+      !$acc enter data create(meridWSEdge, zonalWSEdge)
+#endif
+
+      if (config_bulk_wind_stress_interp_isotropic) then
+
+        call mpas_vector_cell_to_edge_isotropic(windStressZonal, windStressMerid, &
+          nEdges, vertexDegree, verticesOnEdge, cellsOnVertex, kiteAreasOnVertex, &
+          zonalWSEdge, meridWSEdge)
+
+      else
+
+        call mpas_vector_cell_to_edge_anisotropic(windStressZonal, windStressMerid, &
+          nEdges, cellsOnEdge, zonalWSEdge, meridWSEdge)
+
+      end if
+
 #ifdef MPAS_OPENACC
       !$acc parallel loop &
-      !$acc    present(cellsOnEdge, angleEdge, sfcStress, &
-      !$acc            windStressZonal, windStressMerid) &
-      !$acc    private(cell1, cell2, zonalWSEdge, meridWSEdge)
+      !$acc    present(angleEdge, sfcStress, &
+      !$acc            zonalWSEdge, meridWSEdge)
 #else
       !$omp parallel
-      !$omp do schedule(runtime) &
-      !$omp    private(cell1, cell2, zonalWSEdge, meridWSEdge)
+      !$omp do schedule(runtime)
 #endif
       do iEdge = 1, nEdges
-        cell1 = cellsOnEdge(1, iEdge)
-        cell2 = cellsOnEdge(2, iEdge)
-
-        zonalWSEdge = 0.5_RKIND * (windStressZonal(cell1) + &
-                                   windStressZonal(cell2))
-        meridWSEdge = 0.5_RKIND * (windStressMerid(cell1) + &
-                                   windStressMerid(cell2))
-
         sfcStress(iEdge) = sfcStress(iEdge) + &
-                           cos(angleEdge(iEdge))*zonalWSEdge + &
-                           sin(angleEdge(iEdge))*meridWSEdge
+                           cos(angleEdge(iEdge))*zonalWSEdge(iEdge) + &
+                           sin(angleEdge(iEdge))*meridWSEdge(iEdge)
       end do
 #ifndef MPAS_OPENACC
       !$omp end do
@@ -257,6 +265,12 @@ subroutine ocn_surface_bulk_forcing_vel( &
       !$omp end parallel
 #endif
 
+#ifdef MPAS_OPENACC
+      !$acc exit data delete(zonalWSEdge, meridWSEdge)
+#endif
+
+      deallocate(zonalWSEdge, meridWSEdge)
+
       !*** Transfer any needed data back to host and delete inputs
       !$acc exit data &
       !$acc    delete(windStressZonal, windStressMerid)
@@ -565,7 +579,7 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer
                       * max(tracerGroup(index_temperature_flux,minLevelCell(iCell),iCell), 0.0_RKIND) / rho_sw
 
            ! Accumulate fluxes that use the freezing point
-! mrp performance note: should call ocn_freezing_temperature just once here 
+! mrp performance note: should call ocn_freezing_temperature just once here
            seaIceTemperatureFlux(iCell) = seaIceFreshWaterFlux(iCell) * &
                ocn_freezing_temperature( tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell), pressure=0.0_RKIND, &
                                          inLandIceCavity=.false.) / rho_sw

From 4b7397f2fc20d42e7fbbae2f50b5c3b032a1361c Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Sat, 11 May 2024 11:14:07 -0500
Subject: [PATCH 4/5] Add config_bulk_wind_stress_interp_isotropic to E3SM nml

---
 components/mpas-ocean/bld/build-namelist               |  1 +
 components/mpas-ocean/bld/build-namelist-section       |  1 +
 .../bld/namelist_files/namelist_defaults_mpaso.xml     |  1 +
 .../bld/namelist_files/namelist_definition_mpaso.xml   | 10 +++++++++-
 4 files changed, 12 insertions(+), 1 deletion(-)

diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 0b8d3c782fa7..54bc9d586d87 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -690,6 +690,7 @@ add_default($nl, 'config_gotm_constant_bottom_drag_coeff');
 
 add_default($nl, 'config_use_variable_drag');
 add_default($nl, 'config_use_bulk_wind_stress');
+add_default($nl, 'config_bulk_wind_stress_interp_isotropic');
 add_default($nl, 'config_use_bulk_thickness_flux');
 add_default($nl, 'config_flux_attenuation_coefficient');
 add_default($nl, 'config_flux_attenuation_coefficient_runoff');
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index 28030392442d..f2935c7c25c7 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -219,6 +219,7 @@ add_default($nl, 'config_gotm_constant_bottom_drag_coeff');
 
 add_default($nl, 'config_use_variable_drag');
 add_default($nl, 'config_use_bulk_wind_stress');
+add_default($nl, 'config_bulk_wind_stress_interp_isotropic');
 add_default($nl, 'config_use_bulk_thickness_flux');
 add_default($nl, 'config_flux_attenuation_coefficient');
 add_default($nl, 'config_flux_attenuation_coefficient_runoff');
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
index 4069eb144cff..5596662aa469 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -284,6 +284,7 @@
 <!-- forcing -->
 <config_use_variable_drag>.false.</config_use_variable_drag>
 <config_use_bulk_wind_stress>.true.</config_use_bulk_wind_stress>
+<config_bulk_wind_stress_interp_isotropic>.false.</config_bulk_wind_stress_interp_isotropic>
 <config_use_bulk_thickness_flux>.true.</config_use_bulk_thickness_flux>
 <config_flux_attenuation_coefficient>0.001</config_flux_attenuation_coefficient>
 <config_flux_attenuation_coefficient_runoff>10.0</config_flux_attenuation_coefficient_runoff>
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
index 01456cd6583a..f4147c7e281c 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -1148,6 +1148,14 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_bulk_wind_stress_interp_isotropic" type="logical"
+	category="forcing" group="forcing">
+Controls if windstress is interpolated to edges from 4 adjacent cells rather than just 2.
+
+Valid values: .true. or .false.
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_use_bulk_thickness_flux" type="logical"
 	category="forcing" group="forcing">
 Controls if a bulk thickness flux will be computed for surface forcing.
@@ -1188,7 +1196,7 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_sw_absorption_type" type="char*1024"
 	category="shortwaveRadiation" group="shortwaveRadiation">
-Name of shortwave absorption type used in simulation. 
+Name of shortwave absorption type used in simulation.
 
 Valid values: 'jerlov' or 'ohlmann00' or 'none'
 Default: Defined in namelist_defaults.xml

From 4f9c3eda9647bc1c9fb562ff35d11ba10a116964 Mon Sep 17 00:00:00 2001
From: Xylar Asay-Davis <xylarstorm@gmail.com>
Date: Sat, 11 May 2024 11:14:39 -0500
Subject: [PATCH 5/5] Add isotropic_wind_stress stealth test

---
 .../testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README | 4 ++++
 .../testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso   | 1 +
 2 files changed, 5 insertions(+)
 create mode 100644 components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README
 create mode 100644 components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso

diff --git a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README
new file mode 100644
index 000000000000..a2889bc0afc2
--- /dev/null
+++ b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README
@@ -0,0 +1,4 @@
+This testdef is used to test a stealth feature in mpaso introduced by
+PR #XXXX. It simplies changes one mpaso namelist variable,
+    config_bulk_wind_stress_interp_isotropic
+from its default value of false to true.
diff --git a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso
new file mode 100644
index 000000000000..2780fe82a5a5
--- /dev/null
+++ b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso
@@ -0,0 +1 @@
+ config_bulk_wind_stress_interp_isotropic = .true.