Skip to content

Commit

Permalink
Merge branch 'alexolinhager/subcycleAverageEffectivePressure' into MA…
Browse files Browse the repository at this point in the history 
…LI-Dev/develop

This PR defines a new variable, effectivePressureSGH, that is identical
to the previous effectivePressure variable. effectivePressureSGH is
calculated every hydrology timestep and is only used by the hydrology
model. effectivePressure is now redefined as the average
effectivePressureSGH throughout a dynamics coupling interval, and is
only used by the dynamics model when running in a coupled configuration.
The purpose of this PR is to increase stability when coupling the
hydrology and dynamics models by making the dynamics model respond to
the effectivePressure field throughout the hydrology subcycle, instead
of an arbitrary sampling at its last timestep. This reduces the
likelihood of pressure waves building resonance between the two model
components.

* origin/alexolinhager/subcycleAverageEffectivePressure:
  Remove effectivePressureSGH from ocean connection
  Apply suggestions/debugging from review
  Time average effectivePressure over subcycle
  • Loading branch information
@matthewhoffman
matthewhoffman committed Sep 6, 2024
2 parents d6d0f1f + 81c0d02 commit e7d8d43
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 13 deletions.
4 changes: 3 additions & 1 deletion components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,8 +178,10 @@
<var name="frictionAngle" type="real" dimensions="nCells Time" units="None"
description="subglacial till friction angle" />
<!-- convenience variables -->
<var name="effectivePressureSGH" type="real" dimensions="nCells Time" units="Pa"
description="effective ice pressure in subglacial hydrology system" />
<var name="effectivePressure" type="real" dimensions="nCells Time" units="Pa" packages="higherOrderVelocity"
description="effective ice pressure in subglacial hydrology system" />
description="Effective pressure used in basal friction calculation. If subglacial hydrology model is active, this will be effectivePressureSGH averaged over the subglacial hydrology model timestepping subcycles. If subglacial hydrology model is inactive, this will come from a file or a parameterization."/>
<var name="hydropotential" type="real" dimensions="nCells Time" units="Pa"
description="hydropotential in subglacial hydrology system" />
<var name="waterFlux" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
Expand Down
55 changes: 43 additions & 12 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -316,6 +316,8 @@ subroutine li_SGH_solve(domain, err)
real (kind=RKIND), dimension(:), pointer :: waterVelocity
real (kind=RKIND), dimension(:), pointer :: waterVelocityCellX
real (kind=RKIND), dimension(:), pointer :: waterVelocityCellY
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), allocatable :: cellJunk
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell
Expand All @@ -330,6 +332,7 @@ subroutine li_SGH_solve(domain, err)
integer, pointer :: nCells
integer :: iCell, iEdge, iEdgeOnCell
real (kind=RKIND) :: timeLeft ! subcycling time remaining until master dt is reached
real (kind=RKIND) :: deltatSGHaccumulated
integer :: numSubCycles ! number of subcycles
integer :: err_tmp
Expand Down Expand Up @@ -383,6 +386,19 @@ subroutine li_SGH_solve(domain, err)
call mpas_pool_get_array(meshPool, 'deltat', masterDeltat)
timeLeft = masterDeltat ! in seconds
numSubCycles = 0
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
effectivePressure = 0.0_RKIND
deltatSGHaccumulated = 0.0_RKIND
block => block % next
end do
! =============
! =============
! =============
Expand Down Expand Up @@ -682,7 +698,22 @@ subroutine li_SGH_solve(domain, err)
call mpas_dmpar_field_halo_exch(domain, 'waterPressureSmooth')
call mpas_timer_stop("halo updates")
!Average effectivePressureSGH over coupling interval for use in dynamics model
block => domain % blocklist
do while (associated(block))
call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH)
effectivePressure = (effectivePressure * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH)
block => block % next
end do
deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH
! =============
! =============
! =============
Expand Down Expand Up @@ -1536,7 +1567,7 @@ subroutine calc_pressure(block, err)
real (kind=RKIND), dimension(:), pointer :: waterPressureOld
real (kind=RKIND), dimension(:), pointer :: waterPressureTendency
real (kind=RKIND), dimension(:), pointer :: waterThickness
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: zeroOrderSum
real (kind=RKIND), dimension(:), pointer :: closingRate
real (kind=RKIND), dimension(:), pointer :: openingRate
Expand Down Expand Up @@ -1595,7 +1626,7 @@ subroutine calc_pressure(block, err)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
call mpas_pool_get_array(hydroPool, 'waterPressureOld', waterPressureOld)
call mpas_pool_get_array(hydroPool, 'waterPressureTendency', waterPressureTendency)
Expand Down Expand Up @@ -1625,8 +1656,8 @@ subroutine calc_pressure(block, err)
openingRate = max(0.0_RKIND, openingRate)
closingRate(:) = creepCoeff * flowParamA(nVertLevels, :) * &
effectivePressure(:)**3 * waterThickness(:)
! closingRate(:) = waterThickness(:) * effectivePressure(:) / 1.0e13_RKIND
effectivePressureSGH(:)**3 * waterThickness(:)
! closingRate(:) = waterThickness(:) * effectivePressureSGH(:) / 1.0e13_RKIND
! ! Hewitt 2011 creep closure form. Denominator is ice viscosity
zeroOrderSum = closingRate - openingRate + (basalMeltInput + externalWaterInput) / rho_water - &
Expand Down Expand Up @@ -1729,7 +1760,7 @@ subroutine calc_pressure_diag_vars(block, err)
real (kind=RKIND), dimension(:), pointer :: hydropotentialBase
real (kind=RKIND), dimension(:), pointer :: waterThickness
real (kind=RKIND), dimension(:), pointer :: hydropotential
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: iceThicknessHydro
integer, dimension(:), pointer :: cellMask
real (kind=RKIND), pointer :: config_sea_level
Expand All @@ -1744,7 +1775,7 @@ subroutine calc_pressure_diag_vars(block, err)
call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level)
call mpas_pool_get_config(liConfigs, 'config_ocean_density', rhoo)

call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
call mpas_pool_get_array(hydroPool, 'waterPressure', waterPressure)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase)
Expand All @@ -1753,12 +1784,12 @@ subroutine calc_pressure_diag_vars(block, err)
call mpas_pool_get_array(hydroPool, 'iceThicknessHydro', iceThicknessHydro)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)

effectivePressure = rhoi * gravity * iceThicknessHydro - waterPressure
effectivePressureSGH = rhoi * gravity * iceThicknessHydro - waterPressure
! < this should evalute to 0 for floating ice if Pw set correctly there.
where (.not. (li_mask_is_grounded_ice(cellMask)))
effectivePressure = 0.0_RKIND ! zero effective pressure where no ice to avoid confusion
effectivePressureSGH = 0.0_RKIND ! zero effective pressure where no ice to avoid confusion
end where

hydropotentialBase = rho_water * gravity * bedTopography + waterPressure
where ((.not. li_mask_is_grounded_ice(cellMask)) .and. (bedTopography <= config_sea_level))
hydropotentialBase = 0.0_RKIND !zero hydropotential in ocean
Expand Down Expand Up @@ -1830,7 +1861,7 @@ subroutine update_channel(block, err)
real (kind=RKIND), dimension(:), pointer :: channelChangeRate
real (kind=RKIND), dimension(:), pointer :: flowParamAChannel
real (kind=RKIND), dimension(:), pointer :: channelEffectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressure
real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH
real (kind=RKIND), dimension(:), pointer :: channelDiffusivity
real (kind=RKIND), dimension(:), pointer :: waterThicknessEdgeUpwind
integer, dimension(:), pointer :: waterFluxMask
Expand Down Expand Up @@ -1874,7 +1905,7 @@ subroutine update_channel(block, err)
call mpas_pool_get_array(hydroPool, 'channelEffectivePressure', channelEffectivePressure)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA)
call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure)
call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH)
call mpas_pool_get_array(hydroPool, 'waterFluxMask', waterFluxMask)
call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
call mpas_pool_get_array(hydroPool, 'channelDiffusivity', channelDiffusivity)
Expand Down Expand Up @@ -1938,7 +1969,7 @@ subroutine update_channel(block, err)
! Not sure if these ought to be upwind average, but using centered
flowParamAChannel(iEdge) = 0.5_RKIND * (flowParamA(nVertLevels, cell1) + flowParamA(nVertLevels, cell2) )
channelEffectivePressure(iEdge) = 0.5_RKIND * (effectivePressure(cell1) + effectivePressure(cell2))
channelEffectivePressure(iEdge) = 0.5_RKIND * (effectivePressureSGH(cell1) + effectivePressureSGH(cell2))
end do
channelClosingRate(:) = creep_coeff * channelArea(:) * flowParamAChannel(:) * channelEffectivePressure(:)**3
Expand Down

0 comments on commit e7d8d43

Please sign in to comment.