Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Introduce isotropic surface stress to MPAS-Ocean #6917

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Introduce isotropic surface stress to MPAS-Ocean #6917

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
add9e24
Add framework routines to interp. 2d vector fields
xylar May 10, 2024
1bbdea9
Add OpenACC and OpenMP directives
xylar May 10, 2024
b792406
Add option to use isotropic interp for surf. wind stress
xylar May 11, 2024
4b7397f
Add config_bulk_wind_stress_interp_isotropic to E3SM nml
xylar May 11, 2024
4f9c3ed
Add isotropic_wind_stress stealth test
xylar May 11, 2024
2498846
Merge pull request #1 from xylar/mpas/isotropic-interp-normal-wind-st…
proteanplanet May 13, 2024
e88528d
Switching to quad surface ocean stress interpolation as the default.
proteanplanet Jan 23, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
183 changes: 183 additions & 0 deletions components/mpas-framework/src/operators/mpas_vector_operations.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -288,6 +290,187 @@ 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

#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)

zonalEdge(iEdge) = 0.5_RKIND * (zonalCell(cell1) + &
zonalCell(cell2))
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!}}}



!***********************************************************************
!
! 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

#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
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
#ifndef MPAS_OPENACC
!$omp end do
!$omp end parallel
#endif

end subroutine mpas_vector_cell_to_edge_isotropic!}}}



!***********************************************************************
!
! routine mpas_tangential_velocity
Expand Down
1 change: 1 addition & 0 deletions components/mpas-ocean/bld/build-namelist
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -692,6 +692,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');
Expand Down
1 change: 1 addition & 0 deletions components/mpas-ocean/bld/build-namelist-section
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,6 +221,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');
Expand Down
1 change: 1 addition & 0 deletions components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -332,6 +332,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>
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A question for @cbegeman and @vanroekel: Given the extensive testing in B-cases, should we immediately just switch on isotropic wind stress averaging as the default in master, thus changing, this to true.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would support making isotropic default behavior based on your testing

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would as well.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with changing the default to true

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done! Just left for @cbegeman to approve.

<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>
Expand Down
8 changes: 8 additions & 0 deletions components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1164,6 +1164,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.
Expand Down
4 changes: 4 additions & 0 deletions components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/README
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
1 change: 1 addition & 0 deletions ...s/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/isotropic_wind_stress/user_nl_mpaso
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
config_bulk_wind_stress_interp_isotropic = .true.
4 changes: 4 additions & 0 deletions components/mpas-ocean/src/Registry.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -690,6 +690,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."
Expand Down
60 changes: 37 additions & 23 deletions components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
!-----------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
Loading