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

Fix split-explicit upwinded advection #6047

Merged
merged 2 commits into from
Dec 14, 2023
Merged

Fix split-explicit upwinded advection #6047

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5353901
Use deltaSSH and baroclinicThickness not bottomDepthEdge
cbegeman Oct 23, 2023
8ef6785
Use nCellsAll for deltaSSH
cbegeman Oct 23, 2023
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
45 changes: 28 additions & 17 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,8 +199,9 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
real (kind=RKIND), dimension(:), allocatable :: &
btrvel_temp, &
sshTemp, &
deltaSSH, &
wctEdge, &! flux water column thickness at edge
bottomDepthEdge
baroclinicThickness

! State Array Pointers
real (kind=RKIND), dimension(:), pointer :: &
Expand Down Expand Up @@ -391,7 +392,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2)

allocate(bottomDepthEdge(nEdgesAll+1))
allocate(baroclinicThickness(nEdgesAll+1))

if (config_transport_tests_flow_id > 0) then
! This is a transport test. Write advection velocity from prescribed
Expand Down Expand Up @@ -798,20 +799,18 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp end do
!$omp end parallel

bottomDepthEdge(:) = 0.0_RKIND
baroclinicThickness(:) = 0.0_RKIND
!$omp parallel
!$omp do schedule(runtime) &
!$omp private(cell1, cell2, k, thicknessSum)
!$omp private(cell1, cell2, k)
do iEdge = 1, nEdgesHalo(config_num_halos+1)
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
thicknessSum = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge)
baroclinicThickness(iEdge) = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge)
do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + &
baroclinicThickness(iEdge) = baroclinicThickness(iEdge) + &
layerThickEdgeFlux(k,iEdge)
enddo
bottomDepthEdge(iEdge) = thicknessSum &
- 0.5_RKIND*(sshNew(cell1) + sshNew(cell2))
enddo ! iEdge
!$omp end do
!$omp end parallel
Expand Down Expand Up @@ -954,6 +953,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
btrvel_temp(:) = 0.0_RKIND
allocate(wctEdge(nEdgesAll+1))
wctEdge(:) = 0.0_RKIND
allocate(deltaSSH(nCellsAll))
deltaSSH(:) = 0.0_RKIND

cellHaloComputeCounter = 0
edgeHaloComputeCounter = 0
Expand Down Expand Up @@ -1126,8 +1127,15 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

! Compute barotropic thickness at the edge. Note this
! matches the sum of baroclinic thicknesses at this edge.
call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshSubcycleCur, &
bottomDepthEdge, wctEdge)
!$omp parallel
!$omp do schedule(runtime)
do iCell = 1, nCellsAll
deltaSSH(iCell) = sshSubcycleCur(iCell) - sshNew(iCell)
end do
!$omp end do
!$omp end parallel
call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, &
baroclinicThickness, wctEdge)

!$omp parallel
!$omp do schedule(runtime) &
Expand Down Expand Up @@ -1376,12 +1384,13 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

!$omp parallel
!$omp do schedule(runtime)
do iCell = 1, nCells
do iCell = 1, nCellsAll
! SSH is linear combination of SSHold and SSHnew
sshTemp(iCell) = (1-config_btr_gam2_SSHWt1)* &
sshSubcycleCur(iCell) + &
config_btr_gam2_SSHWt1 * &
sshSubcycleNew(iCell)
deltaSSH(iCell) = sshTemp(iCell) - sshNew(iCell)
end do ! cell loop for ssh
!$omp end do
!$omp end parallel
Expand All @@ -1397,8 +1406,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp end do
!$omp end parallel

call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshTemp, &
bottomDepthEdge, wctEdge)
call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, &
baroclinicThickness, wctEdge)

!$omp parallel
!$omp do schedule(runtime) &
Expand Down Expand Up @@ -1428,18 +1437,19 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

!$omp parallel
!$omp do schedule(runtime)
do iCell = 1, nCells
do iCell = 1, nCellsAll
! SSH is linear combination of SSHold and SSHnew
sshTemp(iCell) = (1-config_btr_gam2_SSHWt1)* &
sshSubcycleCur(iCell) + &
config_btr_gam2_SSHWt1 * &
sshSubcycleNew(iCell)
deltaSSH(iCell) = sshTemp(iCell) - sshNew(iCell)
end do ! cell loop for ssh
!$omp end do
!$omp end parallel

call ocn_diagnostic_solve_wctEdge(btrvel_temp, sshTemp, &
bottomDepthEdge, wctEdge)
call ocn_diagnostic_solve_wctEdge(btrvel_temp, deltaSSH, &
baroclinicThickness, wctEdge)

! compute barotropic thickness flux on edges
!$omp parallel
Expand Down Expand Up @@ -1498,6 +1508,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

deallocate(btrvel_temp)
deallocate(wctEdge)
deallocate(deltaSSH)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END Barotropic subcycle loop
Expand Down Expand Up @@ -2704,7 +2715,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
call mpas_dmpar_exch_halo_field(effectiveDensityField)
call mpas_timer_stop("se effective density halo")
end if
deallocate(bottomDepthEdge)
deallocate(baroclinicThickness)

call mpas_timer_stop('se fini')
call mpas_timer_stop("se timestep")
Expand Down
39 changes: 19 additions & 20 deletions components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -681,20 +681,19 @@ end subroutine ocn_relativeVorticity_circulation!}}}
!
!-----------------------------------------------------------------------

subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, sshCell, &
bottomDepthEdge, wctEdge)!{{{
subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, deltaSSH, &
baroclinicThickness, wctEdge)!{{{

!-----------------------------------------------------------------
! input variables
!-----------------------------------------------------------------

real (kind=RKIND), dimension(:), intent(in) :: &
btrVelocity, &!< barotropic velocity
sshCell, &!< ssh at cells, used for centered water column
!< thickness
bottomDepthEdge !< bottom depth centered at edges, used for
!< centered water column thickness
! bottomDepth is already present
btrVelocity, &!< barotropic velocity
deltaSSH, &!< ssh change from baroclinic step to barotropic
!< substep
baroclinicThickness !< water column thickness on edges at the last
!< baroclinic step (centered or upwinded)

!-----------------------------------------------------------------
! output variables
Expand Down Expand Up @@ -725,9 +724,9 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, sshCell, &

! Compute mean layer thickness at edge
#ifdef MPAS_OPENACC
!$acc enter data copyin(bottomDepthEdge, sshCell, wctEdge)
!$acc enter data copyin(baroclinicThickness, deltaSSH, wctEdge)
!$acc parallel loop &
!$acc present(bottomDepthEdge, sshCell, wctEdge, cellsOnEdge) &
!$acc present(baroclinicThickness, deltaSSH, wctEdge, cellsOnEdge) &
!$acc private(iEdge, cell1, cell2)
#else
!$omp parallel
Expand All @@ -737,25 +736,25 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, sshCell, &
if ( maxLevelEdgeTop(iEdge) == 0 ) cycle
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
wctEdge(iEdge) = 0.5_RKIND*(sshCell(cell1) + sshCell(cell2)) &
+ bottomDepthEdge(iEdge)
wctEdge(iEdge) = 0.5_RKIND*(deltaSSH(cell1) + deltaSSH(cell2)) &
+ baroclinicThickness(iEdge)
end do
#ifndef MPAS_OPENACC
!$omp end do
!$omp end parallel
#endif

#ifdef MPAS_OPENACC
!$acc exit data copyout(wctEdge) delete(bottomDepthEdge, sshCell)
!$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH)
#endif

case (thickEdgeFluxUpwind)

! Use upwind thickness as the edge flux value
#ifdef MPAS_OPENACC
!$acc enter data copyin(bottomDepthEdge, sshCell, btrVelocity, wctEdge)
!$acc enter data copyin(baroclinicThickness, deltaSSH, btrVelocity, wctEdge)
!$acc parallel loop &
!$acc present(bottomDepthEdge, btrVelocity, sshCell, cellsOnEdge, maxLevelEdgeTop, wctEdge) &
!$acc present(baroclinicThickness, btrVelocity, deltaSSH, cellsOnEdge, maxLevelEdgeTop, wctEdge) &
!$acc private(cell1, cell2, iEdge)
#else
!$omp parallel
Expand All @@ -766,12 +765,12 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, sshCell, &
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
if (btrVelocity(iEdge) > 0.0_RKIND) then
wctEdge(iEdge) = sshCell(cell1) + bottomDepthEdge(iEdge)
wctEdge(iEdge) = deltaSSH(cell1) + baroclinicThickness(iEdge)
elseif (btrVelocity(iEdge) < 0.0_RKIND) then
wctEdge(iEdge) = sshCell(cell2) + bottomDepthEdge(iEdge)
wctEdge(iEdge) = deltaSSH(cell2) + baroclinicThickness(iEdge)
else
wctEdge(iEdge) = max(sshCell(cell1), sshCell(cell2)) &
+ bottomDepthEdge(iEdge)
wctEdge(iEdge) = max(deltaSSH(cell1), deltaSSH(cell2)) &
+ baroclinicThickness(iEdge)
end if
end do
#ifndef MPAS_OPENACC
Expand All @@ -780,7 +779,7 @@ subroutine ocn_diagnostic_solve_wctEdge(btrVelocity, sshCell, &
#endif

#ifdef MPAS_OPENACC
!$acc exit data copyout(wctEdge) delete(bottomDepthEdge, sshCell, btrVelocity)
!$acc exit data copyout(wctEdge) delete(baroclinicThickness, deltaSSH, btrVelocity)
#endif

case default
Expand Down