From e5b7f0ca619f97fe4ad752bc277fb5fe5f52da60 Mon Sep 17 00:00:00 2001
From: Alexander Hager <ahager@lanl.gov>
Date: Mon, 12 Jun 2023 11:01:31 -0600
Subject: [PATCH] Create hydroTerrestrialMarginMask

Creates hydroTerrestrialMarginMask variable, which designates the
terrestrial margins of the active subglacial hydrology model
---
 .../src/Registry_subglacial_hydro.xml            |  4 +++-
 .../src/mode_forward/mpas_li_subglacial_hydro.F  | 16 ++++++++++++++++
 2 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
index 602b7c11f6d5..66a4a306bbf3 100644
--- a/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
+++ b/components/mpas-albany-landice/src/Registry_subglacial_hydro.xml
@@ -174,7 +174,9 @@
                      description="mask indicating how to handle fluxes on each edge: 0=calculate based on hydropotential gradient; 1=allow outflow based on hydropotential gradient, but no inflow (NOT YET IMPLEMENTED); 2=zero flux" />
                 <var name="hydroMarineMarginMask" type="integer" dimensions="nEdges Time" units="none"
                      description="mask indicating the marine boundary of the active subglacial hydrology domain" />
-                <var name="waterFluxAdvec" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
+		<var name="hydroTerrestrialMarginMask" type="integer" dimensions="nEdges Time" units="none"
+	 	     description="mask indicating the terrestrial boundary of the active subglacial hydrology domain" />
+		<var name="waterFluxAdvec" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
                      description="advective water flux in subglacial hydrology system" />
                 <var name="waterFluxDiffu" type="real" dimensions="nEdges Time" units="m{^2} s^{-1}"
                      description="diffusive water flux in subglacial hydrology system" />
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 f7118b3fe850..f73177ab3907 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
@@ -2183,6 +2183,7 @@ subroutine calc_hydro_mask(domain)
       type (mpas_pool_type), pointer :: meshPool
       real (kind=RKIND), dimension(:), pointer :: bedTopography
       integer, dimension(:), pointer :: hydroMarineMarginMask
+      integer, dimension(:), pointer :: hydroTerrestrialMarginMask
       integer, dimension(:,:), pointer :: cellsOnEdge
       integer, dimension(:), pointer :: cellMask
       integer, pointer :: nEdgesSolve
@@ -2199,6 +2200,7 @@ subroutine calc_hydro_mask(domain)
          call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
          call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
          call mpas_pool_get_array(hydroPool, 'hydroMarineMarginMask', hydroMarineMarginMask)
+         call mpas_pool_get_array(hydroPool, 'hydroTerrestrialMarginMask', hydroTerrestrialMarginMask)
          call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
          call mpas_pool_get_dimension(meshPool, 'nEdgesSolve', nEdgesSolve)
 
@@ -2218,6 +2220,20 @@ subroutine calc_hydro_mask(domain)
             endif
          enddo
 
+         hydroTerrestrialMarginMask(:) = 0
+         do iEdge = 1, nEdgesSolve
+            cell1 = cellsOnEdge(1, iEdge)
+            cell2 = cellsOnEdge(2, iEdge)
+            !Look for edges with 1 cell on grounding ice and the other cell on land without ice
+            if ((li_mask_is_grounded_ice(cellMask(cell1))) .and. ( .not. li_mask_is_ice(cellMask(cell2))) &
+               .and. (bedTopography(cell2) > config_sea_level)) then
+               hydroTerrestrialMarginMask(iEdge) = 1
+            elseif ((li_mask_is_grounded_ice(cellMask(cell2))) .and. ( .not. li_mask_is_ice(cellMask(cell1))) &
+               .and. (bedTopography(cell1) > config_sea_level)) then
+               hydroTerrestrialMarginMask(iEdge) = 1
+            endif
+         enddo
+
          block => block % next
       end do