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.

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

Switch wind stress to use isotropic interpolation from cells to edge normals #98

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
Original file line number Diff line number Diff line change
Expand Up @@ -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');
Expand Down
1 change: 1 addition & 0 deletions components/mpas-ocean/bld/build-namelist-section
Original file line number Diff line number Diff line change
Expand Up @@ -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');
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
Expand Down
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.
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
Original file line number Diff line number Diff line change
Expand Up @@ -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."
Expand Down
60 changes: 37 additions & 23 deletions components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F
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