Calculate grounding line flux after advection but before updating thi…

…ckness Calculate grounding line flux after advection but before updating thickness, which is required to close mass budgets.
MALI-Dev · Nov 17, 2023 · fe5b000 · fe5b000
1 parent 40723e7
commit fe5b000
Showing 1 changed file with 42 additions and 41 deletions.
diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
@@ -381,47 +381,6 @@ subroutine li_advection_thickness_tracers(&
 
       if ( (trim(config_thickness_advection) .eq. 'fo') .or. &
           (trim(config_thickness_advection) .eq. 'fct') ) then
-         ! Calculate flux across grounding line
-         ! Do this after new thickness & mask have been calculated, including SMB/BMB
-         fluxAcrossGroundingLine(:) = 0.0_RKIND
-         fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND
-         do iEdge = 1, nEdgesSolve
-            if (li_mask_is_grounding_line(edgeMask(iEdge))) then
-               iCell1 = cellsOnEdge(1,iEdge)
-               iCell2 = cellsOnEdge(2,iEdge)
-               if (li_mask_is_grounded_ice(cellMask(iCell1))) then
-                  GLfluxSign = 1.0_RKIND ! edge sign convention is positive from iCell1 to iCell2 on an edge
-                  theGroundedCell = iCell1
-               else
-                  GLfluxSign = -1.0_RKIND
-                  theGroundedCell = iCell2
-               endif
-               if (trim(config_thickness_advection) == 'fct') then
-                  do k = 1, nVertLevels
-                     fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
-                  enddo
-               else
-                  do k = 1, nVertLevels
-                     layerThicknessEdgeFlux(k,iEdge) = layerNormalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge)
-                     fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
-                  enddo
-               endif
-               ! assign to grounded cell in fluxAcrossGroundingLineOnCells
-               if (thickness(theGroundedCell) <= 0.0_RKIND) then
-                  ! This should never be the case, but checking to avoid possible divide by zero
-                  call mpas_log_write("thickness at a grounding line is unexepectedly <=0", MPAS_LOG_ERR)
-                  err = ior(err, 1)
-                  return
-               endif
-               fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + &
-                       fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density  ! adjust to correct units
-            endif
-         enddo ! edges
-
-         call mpas_timer_start("halo updates")
-         call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLine')
-         call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLineOnCells')
-         call mpas_timer_stop("halo updates")
 
          if (config_print_thickness_advection_info) then
             call mpas_log_write('Using first-order upwind for thickness advection')
@@ -582,6 +541,48 @@ subroutine li_advection_thickness_tracers(&
          call mpas_dmpar_field_halo_exch(domain, 'layerThickness')
          call mpas_timer_stop("halo updates")
 
+         ! Calculate flux across grounding line
+         ! Do this before new masks and thickness have been calculated, and before SMB/BMB
+         fluxAcrossGroundingLine(:) = 0.0_RKIND
+         fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND
+         do iEdge = 1, nEdgesSolve
+            if (li_mask_is_grounding_line(edgeMask(iEdge))) then
+               iCell1 = cellsOnEdge(1,iEdge)
+               iCell2 = cellsOnEdge(2,iEdge)
+               if (li_mask_is_grounded_ice(cellMask(iCell1))) then
+                  GLfluxSign = 1.0_RKIND ! edge sign convention is positive from iCell1 to iCell2 on an edge
+                  theGroundedCell = iCell1
+               else
+                  GLfluxSign = -1.0_RKIND
+                  theGroundedCell = iCell2
+               endif
+               if (trim(config_thickness_advection) == 'fct') then
+                  do k = 1, nVertLevels
+                     fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
+                  enddo
+               else
+                  do k = 1, nVertLevels
+                     layerThicknessEdgeFlux(k,iEdge) = layerNormalVelocity(k,iEdge) * layerThicknessEdge(k,iEdge)
+                     fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * layerThicknessEdgeFlux(k,iEdge)
+                  enddo
+               endif
+               ! assign to grounded cell in fluxAcrossGroundingLineOnCells
+               if (thickness(theGroundedCell) <= 0.0_RKIND) then
+                  ! This should never be the case, but checking to avoid possible divide by zero
+                  call mpas_log_write("thickness at a grounding line is unexepectedly <=0", MPAS_LOG_ERR)
+                  err = ior(err, 1)
+                  return
+               endif
+               fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + &
+                       fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density  ! adjust to correct units
+            endif
+         enddo ! edges
+
+         call mpas_timer_start("halo updates")
+         call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLine')
+         call mpas_dmpar_field_halo_exch(domain, 'fluxAcrossGroundingLineOnCells')
+         call mpas_timer_stop("halo updates")
+
          ! Calculate dynamicThickening (layerThickness is updated by advection at this point, while thickness is still old)
          dynamicThickening = (sum(layerThickness, 1) - thickness) / dt * scyr  ! units of m/yr