Skip to content

Commit

Permalink
Merge branch 'mark-petersen/ocean/barotropic-edge-thickness' (PR #5195)
Browse files Browse the repository at this point in the history 
Fix barotropic thickness on edge in split explicit subcycling

Barotropic thickness at the edge must match the sum of baroclinic
thicknesses at that edge. This PR corrects a long-standing problem of
noise in the vertVelocityTop variable, which is the Eulerian
(fixed-frame) velocity for momentum. The vertTransportVelocityTop
variable for tracer and thickness transport was smooth, both before and
after this change.

Currently in the code, the thickness at an edge used for volume (SSH)
transport in the barotropic subcycling is
   thicknessSum = sshEdge + min(bottomDepth(cell1), bottomDepth(cell2))
That is, the thickness at an edge is taken as the shallower of the two
neighboring cells. This was inconsistent with the baroclinic mode, which
takes the average thickness of the two neighboring cells.

[CC]
  • Loading branch information
@jonbob
jonbob committed Oct 5, 2022
2 parents 0075116 + ef076d4 commit 4dcc8db
Showing 1 changed file with 37 additions and 49 deletions.
86 changes: 37 additions & 49 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
edgeHaloComputeCounter, &! halo counters to reduce
cellHaloComputeCounter ! halo updates during cycling

integer, dimension(:), pointer :: nEdgesArray

real (kind=RKIND) :: &
normalThicknessFluxSum, &! sum of thickness flux in column
thicknessSum, &! sum of thicknesses in column
Expand All @@ -191,7 +193,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
uTemp

real (kind=RKIND), dimension(:), allocatable :: &
btrvel_temp
btrvel_temp, &
bottomDepthEdge

! State Array Pointers
real (kind=RKIND), dimension(:), pointer :: &
Expand Down Expand Up @@ -382,6 +385,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

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

call mpas_pool_get_dimension(meshPool,'nEdgesArray', nEdgesArray)

allocate(bottomDepthEdge(nEdgesAll+1))

if (config_transport_tests_flow_id > 0) then
! This is a transport test. Write advection velocity from prescribed
! flow field.
Expand Down Expand Up @@ -768,6 +775,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
enddo
barotropicForcing(iEdge) = splitFact* &
normalThicknessFluxSum/thicknessSum/dt
bottomDepthEdge(iEdge) = thicknessSum &
- 0.5_RKIND*(sshNew(cell1) + sshNew(cell2))
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
! These two steps are together here:
Expand All @@ -786,6 +795,23 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp end do
!$omp end parallel

!$omp parallel
!$omp do schedule(runtime) &
!$omp private(k, thicknessSum)
do iEdge = nEdgesOwned+1, nEdgesArray(4)
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
thicknessSum = layerThickEdgeFlux(minLevelEdgeBot(iEdge),iEdge)
do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + &
layerThickEdgeFlux(k,iEdge)
enddo
bottomDepthEdge(iEdge) = thicknessSum &
- 0.5_RKIND*(sshNew(cell1) + sshNew(cell2))
enddo ! iEdge
!$omp end do
!$omp end parallel

deallocate(uTemp)

call mpas_timer_start("se halo normalBaroclinicVelocity")
Expand Down Expand Up @@ -1089,27 +1115,17 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
sshTend(iCell) = 0.0_RKIND
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
if (maxLevelEdgeTop(iEdge).eq.0) cycle

cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)

sshEdge = 0.5_RKIND*(sshSubcycleCur(cell1) + &
sshSubcycleCur(cell2))

! method 0: orig, works only without pbc:
!thicknessSum = sshEdge + &
! refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

! method 1, matches method 0 without pbcs, works with pbcs.
thicknessSum = sshEdge + min(bottomDepth(cell1), &
bottomDepth(cell2))

! method 2: may be better than method 1.
! Take average of full thickness at two
! neighboring cells.
!thicknessSum = sshEdge + 0.5* &
! (bottomDepth(cell1) + bottomDepth(cell2))

! Compute barotropic thickness at the edge. Note this
! matches the sum of baroclinic thicknesses at this edge.
thicknessSum = sshEdge + bottomDepthEdge(iEdge)

flux = ((1.0-config_btr_gam1_velWt1)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1149,16 +1165,14 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp do schedule(runtime) &
!$omp private(cell1, cell2, sshEdge,thicknessSum,flux)
do iEdge = 1, nEdges
if (maxLevelEdgeTop(iEdge).eq.0) cycle
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)

sshEdge = 0.5_RKIND*(sshSubcycleCur(cell1) + &
sshSubcycleCur(cell2))

! method 1, matches method 0 without pbcs,
! works with pbcs.
thicknessSum = sshEdge + min(bottomDepth(cell1), &
bottomDepth(cell2))
thicknessSum = sshEdge + bottomDepthEdge(iEdge)

flux = ((1.0-config_btr_gam1_velWt1)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1370,6 +1384,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
sshTend(iCell) = 0.0_RKIND
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnCell(i, iCell)
if (maxLevelEdgeTop(iEdge).eq.0) cycle

cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
Expand All @@ -1387,21 +1402,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

sshEdge = 0.5_RKIND*(sshCell1 + sshCell2)

! method 0: orig, works only without pbc:
!thicknessSum = sshEdge + &
!refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

! method 1, matches method 0 without pbcs,
! but works with pbcs.
thicknessSum = sshEdge + &
min(bottomDepth(cell1), &
bottomDepth(cell2))

! method 2: may be better than method 1.
! take average of full thickness at two
! neighboring cells
!thicknessSum = sshEdge + 0.5* &
! (bottomDepth(cell1) + bottomDepth (cell2))
thicknessSum = sshEdge + bottomDepthEdge(iEdge)

flux = ((1.0-config_btr_gam3_velWt2)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1441,21 +1442,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
sshSubcycleNew(cell2)
sshEdge = 0.5_RKIND * (sshCell1 + sshCell2)

! method 0: orig, works only without pbc:
!thicknessSum = sshEdge + &
!refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

! method 1, matches method 0 without pbcs,
! but works with pbcs.
thicknessSum = sshEdge + min(bottomDepth(cell1),&
bottomDepth(cell2))

! method 2, better, I think.
! take average of full thickness at two
! neighboring cells
!thicknessSum = sshEdge + 0.5* &
! (bottomDepth(cell1) + &
! bottomDepth(cell2) )
thicknessSum = sshEdge + bottomDepthEdge(iEdge)

flux = ((1.0-config_btr_gam3_velWt2) * &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -2695,6 +2682,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)
call mpas_timer_stop('se fini')
call mpas_timer_stop("se timestep")
Expand Down

0 comments on commit 4dcc8db

Please sign in to comment.