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

Update AB2 (mpas_ocn_time_integration_split_ab2.F) code #6323

Closed
wants to merge 1 commit into from
Closed

Update AB2 (mpas_ocn_time_integration_split_ab2.F) code #6323

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
135 changes: 116 additions & 19 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,9 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -276,6 +277,10 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
! Forcing pool
real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta

real (kind=RKIND) :: totalDepth
real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue
real (kind=RKIND), dimension(:,:), pointer :: restingThickness

real (kind=RKIND), dimension(:), pointer :: &
seaIcePressure, atmosphericPressure

Expand Down Expand Up @@ -406,7 +411,7 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -908,20 +913,18 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -1015,6 +1018,8 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -1089,6 +1094,30 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{

end if ! tidal potential forcing

! direct application of tidal boundary condition
if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then
call mpas_pool_get_array(forcingPool, &
'sshSubcycleCurWithTides', &
sshSubcycleCurWithTides)
call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask)
call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue)
!$omp parallel
!$omp do schedule(runtime)
do iCell=1, nCellsAll
! boolean mask for now, could generalize to tappered sponge layer
! only modify layer thicknesses on tidal boundary
do k = minLevelCell(iCell), maxLevelCell(iCell)
sshSubcycleCurWithTides(iCell) = sshSubcycleCur(iCell) + &
(tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleCur(iCell)))
end do
end do
!$omp end do
!$omp end parallel
call mpas_pool_get_array(forcingPool, &
'sshSubcycleCurWithTides', &
sshSubcycleCur)
end if

if (edgeHaloComputeCounter <= 1) then
nEdges = nEdgesOwned
else
Expand Down Expand Up @@ -1187,8 +1216,16 @@ subroutine ocn_time_integrator_split_ab2(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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@hyungyukang commented in a meeting that this is the critical change in this PR.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks from #6047 that this would be BFB when the thickness is centered, so I don't think this will affect the instability we encountered in the v3 HR G-case.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cbegeman , thanks for your comment. When I tested this branch very first time, I was able to see some differences in simple test cases. So I thought this change is critical. But today, I cleaned and compiled again after the meeting, and it is very similar now. Maybe I did something wrong during early tests. However, for code uniformity between split_explicit_ab2.F and split_explicit.F, this PR is worth merging.

This PR is not BFB with very small differences (~1e-9).

Just in case, I will test this PR a bit more.


!$omp parallel
!$omp do schedule(runtime) &
Expand Down Expand Up @@ -1316,6 +1353,38 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
sshSubcycleNew)
end if ! tidal potential forcing

! direct application of tidal boundary condition
if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then
call mpas_pool_get_array(forcingPool, &
'sshSubcycleCurWithTides', &
sshSubcycleCurWithTides)
call mpas_pool_get_array(forcingPool, &
'sshSubcycleNewWithTides', &
sshSubcycleNewWithTides)
call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask)
call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue)
!$omp parallel
!$omp do schedule(runtime)
do iCell=1, nCellsAll
! boolean mask for now, could generalize to tappered sponge layer
! only modify layer thicknesses on tidal boundary
do k = minLevelCell(iCell), maxLevelCell(iCell)
sshSubcycleCurWithTides(iCell) = sshSubcycleCur(iCell) + &
(tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleCur(iCell)))
sshSubcycleNewWithTides(iCell) = sshSubcycleNew(iCell) + &
(tidalInputMask(iCell)*(tidalBCValue(iCell) - sshSubcycleNew(iCell)))
end do
end do
!$omp end do
!$omp end parallel
call mpas_pool_get_array(forcingPool, &
'sshSubcycleCurWithTides', &
sshSubcycleCur)
call mpas_pool_get_array(forcingPool, &
'sshSubcycleNewWithTides', &
sshSubcycleNew)
end if

! Need to initialize btr_vel_temp over one more halo
! than we are computing over
nEdges = nEdgesHalo(min(edgeHaloComputeCounter, &
Expand Down Expand Up @@ -1437,12 +1506,13 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -1458,8 +1528,8 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -1489,18 +1559,19 @@ subroutine ocn_time_integrator_split_ab2(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 @@ -1542,7 +1613,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{

!$omp parallel
!$omp do schedule(runtime)
do iEdge = 1, nEdgesAll
do iEdge = 1, nEdgesOwned
normalBarotropicVelocityNew(iEdge) = &
normalBarotropicVelocityNew(iEdge) + &
normalBarotropicVelocitySubcycleNew(iEdge)
Expand All @@ -1559,6 +1630,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END Barotropic subcycle loop
Expand All @@ -1570,7 +1642,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{

!$omp parallel
!$omp do schedule(runtime)
do iEdge = 1, nEdgesAll
do iEdge = 1, nEdgesOwned
barotropicThicknessFlux(iEdge) = &
barotropicThicknessFlux(iEdge) &
/(nBtrSubcycles*config_btr_subcycle_loop_factor)
Expand Down Expand Up @@ -2831,6 +2903,31 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
#endif
end if

! direct application of tidal boundary condition
if (config_use_tidal_forcing .and. trim(config_tidal_forcing_type) == 'direct') then
call mpas_pool_get_array(verticalMeshPool, 'restingThickness', restingThickness)
call mpas_pool_get_array(forcingPool, 'tidalInputMask', tidalInputMask)
call mpas_pool_get_array(forcingPool, 'tidalBCValue', tidalBCValue)
do iCell=1, nCells
! artificially assumes boolean mask for now, could generalize to tappered sponge layer
if (tidalInputMask(iCell) == 1.0_RKIND) then
! compute total depth for relative thickness contribution
totalDepth = 0.0_RKIND
do k = minLevelCell(iCell), maxLevelCell(iCell)
totalDepth = totalDepth + restingThickness(k,iCell)
end do

! only modify layer thicknesses on tidal boundary
do k = minLevelCell(iCell), maxLevelCell(iCell)
layerThicknessNew(k, iCell) = tidalInputMask(iCell)* &
(tidalBCValue(iCell) + bottomDepth(iCell))* &
(restingThickness(k,iCell)/totalDepth)
!(1.0_RKIND - tidalInputMask(iCell))*layerThicknessNew(k, iCell) ! generalized tappered assumption code
end do
end if
end do
end if

call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, &
scratchPool, tracersPool, 2)

Expand Down Expand Up @@ -2953,7 +3050,7 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
call mpas_dmpar_exch_halo_field(effectiveDensityField)
call mpas_timer_stop("se effective density halo")
end if
deallocate(bottomDepthEdge)
deallocate(baroclinicThickness)

! get current time & setup num TS iteration
nowTime = mpas_get_clock_time(domain%clock, MPAS_NOW, err)
Expand Down