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

Jump to bottom

split explicit corrections #16

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

split explicit corrections #16

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a350b64
Add four split explicit changes. See PRs:
mark-petersen Apr 18, 2022
35de421
turn on gm package (temporary)
mark-petersen Sep 19, 2022
5fa93e9
add config_btr_new5_NoALE
mark-petersen Sep 19, 2022
830652e
add config_btr_new6_EdgeThickness
mark-petersen Sep 20, 2022
1a151a5
change edge thickness from refBD to zTop
mark-petersen Sep 20, 2022
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
36 changes: 36 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 @@ -1263,6 +1263,30 @@
description="number of iterations of the velocity corrector step in stage 2"
possible_values="any positive integer, but typically 1, 2, or 3"
/>
<nml_option name="config_btr_new1_resplit" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_btr_new2_ubtr_foH" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_btr_new3_sub_btrAvgbclV" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_btr_new4_VelCor" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_btr_new5_NoALE" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_btr_new6_EdgeThickness" type="logical" default_value=".false."
description="experimental flag"
possible_values=".true. or .false."
/>
<nml_option name="config_vel_correction" type="logical" default_value=".true."
description="If true, the velocity correction term is included in the horizontal advection of thickness and tracers"
possible_values=".true. or .false."
Expand Down Expand Up @@ -2722,6 +2746,18 @@
description="horizontal velocity, tangential to an edge"
packages="forwardMode;analysisMode"
/>
<var name="BtrAvgOfBaroclinicVelocity" type="real" dimensions="nEdges Time" units="m s^{-1}"
description="barotropic average of baroclinic velocity"
packages="splitTimeIntegrator"
/>
<var name="BtrAvgOfGMBolusVelocity" type="real" dimensions="nEdges Time" units="m s^{-1}"
description="barotropic average of GM Bolus velocity"
packages="splitTimeIntegrator"
/>
<var name="normalVelocityCorrection" type="real" dimensions="nEdges Time" units="m s^{-1}"
description="read it"
packages="splitTimeIntegrator"
/>
<var name="layerThicknessEdgeMean" type="real" dimensions="nVertLevels nEdges Time" units="m"
description="layer thickness averaged from cell center to edges"
packages="forwardMode;analysisMode"
Expand Down
2 changes: 2 additions & 0 deletions components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -356,6 +356,8 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
if (config_use_GM.or.config_use_Redi) then
gmActive = .true.
end if
! mrp TEMPORARY, always turn on gm package
gmActive = .true.

!
! test for use of gm
Expand Down
194 changes: 164 additions & 30 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 @@ -179,7 +179,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
flux, &! temp for computing flux for barotropic
sshEdge, &! sea surface height at edge
CoriolisTerm, &! temp for computing coriolis term (fuperp)
normalVelocityCorrection, &! velocity correction
temp, &! temp for holding vars at new time
temp_h, &! temporary for phi
temp_mask, &! temporary mask (edgeMask)
Expand Down Expand Up @@ -407,6 +406,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! Initialize * variables that are used to compute baroclinic
! tendencies below.

if (.not.config_btr_new1_resplit) then
! The baroclinic velocity needs be recomputed at the beginning
! of a timestep because the implicit vertical mixing is
! conducted on the total u. We keep normalBarotropicVelocity
Expand Down Expand Up @@ -445,11 +445,40 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
end do
#ifndef MPAS_OPENACC
!$omp end do
!$omp end parallel
#endif

else ! config_btr_resplit

!$omp parallel
!$omp do schedule(runtime) private(k, thicknessSum)
do iEdge = 1,nEdgesAll

! Compute barotropic component ubar at the beginning of every time step
thicknessSum = layerThickEdgeMean(1,iEdge)
normalBarotropicVelocityCur(iEdge) = normalVelocityCur(1,iEdge) * layerThickEdgeMean(1,iEdge)
do k = 2, maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + layerThickEdgeMean(k,iEdge)
normalBarotropicVelocityCur(iEdge) = normalBarotropicVelocityCur(iEdge) + normalVelocityCur(k,iEdge) * layerThickEdgeMean(k,iEdge)
enddo
normalBarotropicVelocityCur(iEdge) = normalBarotropicVelocityCur(iEdge) / thicknessSum

! Compute baroclinic component u' = u - ubar at the beginning of every time step
do k = 1, nVertLevels !mrp, could tighten these bounds later: maxLevelEdgeTop % array(iEdge)
normalBaroclinicVelocityCur(k,iEdge) = normalVelocityCur(k,iEdge) - normalBarotropicVelocityCur(iEdge)
normalVelocityNew(k,iEdge) = normalVelocityCur(k,iEdge)
normalBaroclinicVelocityNew(k,iEdge) = normalBaroclinicVelocityCur(k,iEdge)
end do
end do
!$omp end do
!$omp end parallel

end if ! config_btr_resplit

#ifdef MPAS_OPENACC
!$acc loop gang
#else
!$omp parallel
!$omp do schedule(runtime) private(k)
#endif
do iCell = 1, nCellsAll
Expand Down Expand Up @@ -1083,12 +1112,13 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{

!$omp parallel
!$omp do schedule(runtime) &
!$omp private(i, iEdge, cell1, cell2, sshEdge, &
!$omp private(i, k, iEdge, cell1, cell2, sshEdge, &
!$omp thicknessSum, flux)
do iCell = 1, nCells
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 @@ -1100,16 +1130,25 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!thicknessSum = sshEdge + &
! refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

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

else
! mrp note, this works for z-level right now.
! for z-star, could use -zTop(k+1) for refBottomDepth(k)
! but that would not work for general coordinate
! 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))

k = maxLevelEdgeTop(iEdge)
thicknessSum = sshEdge + 0.5* &
( min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+ min(bottomDepth(cell2), -zTop(k+1,cell2)) )
! first attempt, for z-level only:
!( min(bottomDepth(cell1), refBottomDepth(k)) &
! + min(bottomDepth(cell2), refBottomDepth(k)) )
endif

flux = ((1.0-config_btr_gam1_velWt1)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1149,16 +1188,29 @@ 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.
if (.not.config_btr_new6_EdgeThickness) then
! method 1, matches method 0 without pbcs, works with pbcs.
thicknessSum = sshEdge + min(bottomDepth(cell1), &
bottomDepth(cell2))
else
! method 2: may be better than method 1.
! Take average of full thickness at two
! neighboring cells.
k = maxLevelEdgeTop(iEdge)
thicknessSum = sshEdge + 0.5* &
( min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+ min(bottomDepth(cell2), -zTop(k+1,cell2)) )
! first attempt, for z-level only:
!( min(bottomDepth(cell1), refBottomDepth(k)) &
! + min(bottomDepth(cell2), refBottomDepth(k)) )
endif

flux = ((1.0-config_btr_gam1_velWt1)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1370,6 +1422,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 @@ -1391,17 +1444,24 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!thicknessSum = sshEdge + &
!refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

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

else
! 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))
k = maxLevelEdgeTop(iEdge)
thicknessSum = sshEdge + 0.5* &
( min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+ min(bottomDepth(cell2), -zTop(k+1,cell2)) )
! first attempt, for z-level only:
!( min(bottomDepth(cell1), refBottomDepth(k)) &
! + min(bottomDepth(cell2), refBottomDepth(k)) )
endif

flux = ((1.0-config_btr_gam3_velWt2)* &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1445,17 +1505,22 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
!thicknessSum = sshEdge + &
!refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)

! method 1, matches method 0 without pbcs,
! but works with pbcs.
thicknessSum = sshEdge + min(bottomDepth(cell1),&
if (.not.config_btr_new6_EdgeThickness) then
! method 1, matches method 0 without pbcs, 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) )
else
! method 2: may be better than method 1.
! Take average of full thickness at two
! neighboring cells.
k = maxLevelEdgeTop(iEdge)
thicknessSum = sshEdge + 0.5* &
( min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+ min(bottomDepth(cell2), -zTop(k+1,cell2)) )
! first attempt, for z-level only:
!( min(bottomDepth(cell1), refBottomDepth(k)) &
! + min(bottomDepth(cell2), refBottomDepth(k)) )
endif

flux = ((1.0-config_btr_gam3_velWt2) * &
normalBarotropicVelocitySubcycleCur(iEdge) &
Expand Down Expand Up @@ -1556,7 +1621,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! efficiency: This next block of code is really a check for
! debugging, and can be removed later.
call mpas_timer_start('btr se ssh verif')

if (.not.config_btr_new4_VelCor) then
! Correction velocity
! normalVelocityCorrection = (Flux - Sum(h u*))/H
! or, for the full latex version:
Expand Down Expand Up @@ -1624,7 +1689,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
layerThickEdgeFlux(k,iEdge)
enddo

normalVelocityCorrection = useVelocityCorrection* &
normalVelocityCorrection(iEdge) = useVelocityCorrection* &
((barotropicThicknessFlux(iEdge) - &
normalThicknessFluxSum)/thicknessSum)

Expand All @@ -1639,15 +1704,81 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! and tracers in tendency calls in stage 3.
normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge)*( &
normalTransportVelocity(k,iEdge) + &
normalVelocityCorrection)
normalVelocityCorrection(iEdge))
enddo

k = minLevelEdgeBot(iEdge)
thicknessSum = 1.0e-15_RKIND
BtrAvgOfBaroclinicVelocity(iEdge) = 0.0_RKIND
BtrAvgOfGMBolusVelocity(iEdge) = 0.0_RKIND
! note, BtrAvgOfGMBolusVelocity could be changed from an array
! to a scalar if it is not needed to save for later inspection
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + layerThickEdgeFlux(k,iEdge)
BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) + normalBaroclinicVelocityNew(k,iEdge) * layerThickEdgeFlux(k,iEdge)
BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) + normalGMBolusVelocity(k,iEdge) * layerThickEdgeFlux(k,iEdge)
enddo
BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) / thicknessSum
BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) / thicknessSum
if (config_btr_new2_ubtr_foH) then
normalBarotropicVelocityNew(iEdge) = barotropicThicknessFlux(iEdge) / thicknessSum
end if

end do ! iEdge
!$omp end do
!$omp end parallel

deallocate(uTemp)
else ! if (config_btr_new4_VelCor) then
nEdges = nEdgesHalo(config_num_halos+1 ) !mrp check min w nEdgesPtr in PR

!$omp parallel
!$omp do schedule(runtime) private(k, thicknessSum)
do iEdge = 1, nEdges

! thicknessSum is initialized outside the loop because on land boundaries
! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a
! nonzero value to avoid a NaN.
k = minLevelEdgeBot(iEdge)
thicknessSum = 1.0e-15_RKIND
BtrAvgOfBaroclinicVelocity(iEdge) = 0.0_RKIND
BtrAvgOfGMBolusVelocity(iEdge) = 0.0_RKIND
! note, BtrAvgOfGMBolusVelocity could be changed from an array
! to a scalar if it is not needed to save for later inspection
do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
thicknessSum = thicknessSum + layerThickEdgeFlux(k,iEdge)
BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) + normalBaroclinicVelocityNew(k,iEdge) * layerThickEdgeFlux(k,iEdge)
BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) + normalGMBolusVelocity(k,iEdge) * layerThickEdgeFlux(k,iEdge)
enddo
BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) / thicknessSum
BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) / thicknessSum

! Compute barotropic velocity from thickness flux.
! mrp changed back. With this commented out, barotropic velocity
! is computed from average of subcycles above.
! Before, we had
! Set momentum velocity (normalVelocity) to:
! u = F/H + u' - btrAvg( u' )
! with this commented out, we have
! u = subcycleAvg(ubar) + u' - btrAvg( u' )
if (config_btr_new2_ubtr_foH) then
normalBarotropicVelocityNew(iEdge) = barotropicThicknessFlux(iEdge) / thicknessSum
end if

do k = 1, nVertLevels

! Compute transport velocity, used to transport thickness and tracers.
! uTr = F/H + u' - btrAvg( u' ) + uBolus - btrAvg( uBolus )
normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge) &
*( barotropicThicknessFlux(iEdge) / thicknessSum &
+ normalBaroclinicVelocityNew(k,iEdge) - BtrAvgOfBaroclinicVelocity(iEdge) &
+ normalGMBolusVelocity(k,iEdge) - BtrAvgOfGMBolusVelocity(iEdge) )
enddo

end do ! iEdge
!$omp end do
!$omp end parallel
end if
call mpas_timer_stop('btr se ssh verif')

endif ! split_explicit
Expand Down Expand Up @@ -1909,9 +2040,17 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
! u = normalBarotropicVelocity + normalBaroclinicVelocity
! here normalBaroclinicVelocity is at time n+1/2
! This is u used in next iteration or step
if (.not.config_btr_new3_sub_btrAvgbclV) then
normalVelocityNew(k,iEdge) = edgeMask(k,iEdge)* &
( normalBarotropicVelocityNew(iEdge) + &
normalBaroclinicVelocityNew(k,iEdge) )
else !if (.not.config_btr_new3_sub_btrAvgbclV) then

normalVelocityNew(k,iEdge) = edgeMask(k,iEdge)* &
( normalBarotropicVelocityNew(iEdge) + &
normalBaroclinicVelocityNew(k,iEdge) - &
BtrAvgOfBaroclinicVelocity(iEdge) )
end if

enddo
end do ! iEdge
Expand Down Expand Up @@ -2794,12 +2933,7 @@ subroutine ocn_time_integration_split_init(domain)!{{{
nBtrSubcycles = nBtrSubcycles + 1
end if

!*** Set mask for using velocity correction
if (config_vel_correction) then
useVelocityCorrection = 1.0_RKIND
else
useVelocityCorrection = 0.0_RKIND
endif
useVelocityCorrection = 1.0_RKIND ! mrp

!*** Compute some halo needs to reduce the number
!*** of halo updates during subcycling
Expand Down
Loading