From d647f606e36a2c8ccc3106c8787b7948510625be Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Mon, 19 Aug 2024 16:59:32 -0600 Subject: [PATCH 1/3] Time average effectivePressure over subcycle This commit replaces 'effectivePressure' with 'effectivePressureSGH', the version of effective pressure used by the hydrology model. 'effectivePressure' is now defined as effecivePressureSGH time-averaged over the hydrology subcycle, and is what is fed into the dynamics coupler. --- .../src/Registry_subglacial_hydro.xml | 6 +- .../mode_forward/mpas_li_subglacial_hydro.F | 59 ++++++++++++++----- 2 files changed, 47 insertions(+), 18 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml index 8b3f07d4750f..ec00a63d24a3 100644 --- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml +++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml @@ -178,8 +178,10 @@ - + + 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 + + ! ============= ! ============= ! ============= @@ -682,7 +697,19 @@ 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) / (deltatSGHaccumulated + deltatSGH) + deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH + + block => block % next + end do ! ============= ! ============= ! ============= @@ -1536,7 +1563,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 @@ -1595,7 +1622,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) @@ -1625,8 +1652,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 - & @@ -1729,7 +1756,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 @@ -1744,7 +1771,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) @@ -1753,12 +1780,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 @@ -1830,7 +1857,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 @@ -1874,7 +1901,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) @@ -1938,7 +1965,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 @@ -2229,7 +2256,7 @@ subroutine ocean_connection_N(domain) type (mpas_pool_type), pointer :: hydroPool type (mpas_pool_type), pointer :: geometryPool real (kind=RKIND), dimension(:), pointer :: bedTopography - real (kind=RKIND), dimension(:), pointer :: effectivePressure + real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), pointer :: rhoi, rhoo @@ -2243,11 +2270,11 @@ subroutine ocean_connection_N(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure) + call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH) call mpas_pool_get_array(hydroPool, 'thickness', thickness) - effectivePressure = rhoi * gravity * thickness - rhoi * gravity * max(0.0_RKIND, -1.0_RKIND * rhoo/rhoi * bedTopography) - effectivePressure = max(effectivePressure, 0.0_RKIND) ! This is just to zero out N in the open ocean to avoid confusion + effectivePressureSGH = rhoi * gravity * thickness - rhoi * gravity * max(0.0_RKIND, -1.0_RKIND * rhoo/rhoi * bedTopography) + effectivePressureSGH = max(effectivePressureSGH, 0.0_RKIND) ! This is just to zero out N in the open ocean to avoid confusion block => block % next end do From 5385761638e44b68b9e590787eb5dbe9fac53490 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Thu, 29 Aug 2024 10:30:36 -0600 Subject: [PATCH 2/3] Apply suggestions/debugging from review Addresses minor bugs and edits that came up in code review --- .../src/Registry_subglacial_hydro.xml | 4 ++-- .../src/mode_forward/mpas_li_subglacial_hydro.F | 14 +++++++++----- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml index ec00a63d24a3..ac756b6d33f9 100644 --- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml +++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml @@ -180,8 +180,8 @@ - + 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 + effectivePressure = 0.0_RKIND deltatSGHaccumulated = 0.0_RKIND block => block % next @@ -700,16 +701,19 @@ subroutine li_SGH_solve(domain, err) !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) / (deltatSGHaccumulated + deltatSGH) - deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH - + effectivePressure = (effectivePressure * deltatSGHaccumulated + effectivePressureSGH * deltatSGH) / (deltatSGHaccumulated + deltatSGH) + block => block % next end do + + deltatSGHaccumulated = deltatSGHaccumulated + deltatSGH + ! ============= ! ============= ! ============= From 81c0d02530bd11e01ab3fff6ad27605edc35ecc5 Mon Sep 17 00:00:00 2001 From: Alexander Hager Date: Tue, 3 Sep 2024 12:12:52 -0600 Subject: [PATCH 3/3] Remove effectivePressureSGH from ocean connection Reverts effectivePressureSGH back to effectivePressure in ocean_connection_N --- .../src/mode_forward/mpas_li_subglacial_hydro.F | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F index cc210b6bc893..06fb2239107d 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_subglacial_hydro.F @@ -2260,7 +2260,7 @@ subroutine ocean_connection_N(domain) type (mpas_pool_type), pointer :: hydroPool type (mpas_pool_type), pointer :: geometryPool real (kind=RKIND), dimension(:), pointer :: bedTopography - real (kind=RKIND), dimension(:), pointer :: effectivePressureSGH + real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), pointer :: rhoi, rhoo @@ -2274,11 +2274,11 @@ subroutine ocean_connection_N(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'hydro', hydroPool) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(hydroPool, 'effectivePressureSGH', effectivePressureSGH) + call mpas_pool_get_array(hydroPool, 'effectivePressure', effectivePressure) call mpas_pool_get_array(hydroPool, 'thickness', thickness) - effectivePressureSGH = rhoi * gravity * thickness - rhoi * gravity * max(0.0_RKIND, -1.0_RKIND * rhoo/rhoi * bedTopography) - effectivePressureSGH = max(effectivePressureSGH, 0.0_RKIND) ! This is just to zero out N in the open ocean to avoid confusion + effectivePressure = rhoi * gravity * thickness - rhoi * gravity * max(0.0_RKIND, -1.0_RKIND * rhoo/rhoi * bedTopography) + effectivePressure = max(effectivePressure, 0.0_RKIND) ! This is just to zero out N in the open ocean to avoid confusion block => block % next end do