Skip to content

Commit

Permalink
Fix barotropic thickness on edge in split explicit
Browse files Browse the repository at this point in the history 
PR5195
E3SM-Project/E3SM#5195
  • Loading branch information
@mark-petersen
mark-petersen committed Oct 26, 2022
1 parent d4eea06 commit 240deca
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 50 deletions.
86 changes: 37 additions & 49 deletions src/core_ocean/mode_forward/mpas_ocn_time_integration_split.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,6 +164,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 @@ -182,7 +184,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 @@ -354,6 +357,10 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
call mpas_pool_get_array(tendPool, 'lowFreqDivergence', &
lowFreqDivergenceTend)

call mpas_pool_get_dimension(meshPool,'nEdgesArray', nEdgesArray)

allocate(bottomDepthEdge(nEdgesAll+1))

! Initialize * variables that are used to compute baroclinic
! tendencies below.

Expand Down Expand Up @@ -605,6 +612,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 = 1, maxLevelEdgeTop(iEdge)
! These two steps are together here:
Expand All @@ -623,6 +632,23 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!$omp end do
!$omp end parallel

!$omp parallel
!$omp do schedule(runtime) &
!$omp private(cell1, cell2, k, thicknessSum)
do iEdge = nEdgesOwned+1, nEdgesArray(4)
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
thicknessSum = layerThickEdge(1,iEdge)
do k = 2, maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + &
layerThickEdge(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 @@ -905,27 +931,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 @@ -965,16 +981,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 @@ -1180,6 +1194,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 @@ -1197,21 +1212,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 @@ -1251,21 +1252,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 @@ -2020,6 +2007,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
2 changes: 1 addition & 1 deletion src/core_ocean/shared/mpas_ocn_tracer_advection_std.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ subroutine ocn_tracer_advection_std_tend(tracers, adv_coefs, adv_coefs_3rd, nAdv

allocate(high_order_horiz_flux(nVertLevels, nEdges), &
tracer_cur(nVertLevels, nCells), &
high_order_vert_flux(nVertLevels, nCells))
high_order_vert_flux(nVertLevels+1, nCells))

! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality
do iTracer = 1, num_tracers
Expand Down

0 comments on commit 240deca

Please sign in to comment.