Skip to content

Commit

Permalink
Update 'mpas_ocn_time_integration_split_ab2.F' code
Browse files Browse the repository at this point in the history 
* Some codes & routines were updated in split.F, but not in split_ab2.F
* 'ocn_diagnostic_solve_wctEdge' routine was updated, but split_ab2.F was using a call sentence of the old version of the routine.
   - 'deltaSSH' computation was missing in split_ab2.F.
   - 'baroclinicThickness' is computed instead of 'bottomDepthEdge'.
* 'Direct application of tidal boundary condition', which was only in split.F, is added to split_ab2.F.
  • Loading branch information
@hyungyukang
hyungyukang committed Apr 3, 2024
1 parent 6b9ecaa commit eb6b6be
Showing 1 changed file with 116 additions and 19 deletions.
135 changes: 116 additions & 19 deletions components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
View file
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)

!$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

0 comments on commit eb6b6be

Please sign in to comment.