Skip to content

Commit

Permalink
Merge branch 'trhille/higher_order_flux_across_grounding_line' into M…
Browse files Browse the repository at this point in the history 
…ALI-Dev/develop

When using higher-order thickness advection, calculate higher-order flux across grounding line.

* trhille/higher_order_flux_across_grounding_line:
  Generalize layerThicknessEdgeFlux for all thickness advection
  Calculate higher-order fluxAcrossGroundingLine
  Add higher order thickness flux on edges to Registry
  • Loading branch information
@matthewhoffman
matthewhoffman committed Oct 9, 2023
2 parents f5d2e27 + 8be7fbb commit 1343d1d
Show file tree
Hide file tree
Showing 3 changed files with 31 additions and 22 deletions.
3 changes: 3 additions & 0 deletions components/mpas-albany-landice/src/Registry.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1433,6 +1433,9 @@ is the value of that variable from the *previous* time level!
<var name="layerNormalVelocity" type="real" dimensions="nVertLevels nEdges Time" units="m s^{-1}"
description="horizonal velocity, normal component to an edge, layer midpoint"
/>
<var name="layerThicknessEdgeFlux" type="real" dimensions="nVertLevels nEdges Time" units="m^2 s^{-1}"
description="layer-normal thickness flux on edges"
/>
<var name="normalVelocityInitial" type="real" dimensions="nVertInterfaces nEdges Time" units="m s^{-1}"
description="horizonal velocity, normal component to an edge, computed at initialization"
/>
Expand Down
35 changes: 20 additions & 15 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,8 @@ subroutine li_advection_thickness_tracers(&
layerThickness, & ! thickness of each layer
layerThicknessOld, & ! old layer thickness
layerThicknessEdge, & ! layer thickness on upstream edge of cell
normalVelocity ! horizontal velocity on interfaces
normalVelocity, & ! horizontal velocity on interfaces
layerThicknessEdgeFlux ! higher order thickness flux on layers and edges

real (kind=RKIND), dimension(:,:,:), pointer :: &
advectedTracers, & ! values of advected tracers
Expand All @@ -193,7 +194,6 @@ subroutine li_advection_thickness_tracers(&

! Allocatable arrays need for flux-corrected transport advection
real (kind=RKIND), dimension(:,:,:), allocatable :: tend
real (kind=RKIND), dimension(:,:,:), allocatable :: activeTracerHorizontalAdvectionEdgeFlux

integer, dimension(:), pointer :: &
cellMask, & ! integer bitmask for cells
Expand Down Expand Up @@ -292,6 +292,7 @@ subroutine li_advection_thickness_tracers(&

! get arrays from the velocity pool
call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity)
call mpas_pool_get_array(velocityPool, 'layerThicknessEdgeFlux', layerThicknessEdgeFlux)
call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity)
call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine)
call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells)
Expand Down Expand Up @@ -368,6 +369,7 @@ subroutine li_advection_thickness_tracers(&

! save old copycellMask for determining cells changing from grounded to floating and vice versa
cellMaskTemporaryField % array(:) = cellMask(:)
layerThicknessEdgeFlux(:,:) = 0.0_RKIND

!-----------------------------------------------------------------
! Horizontal transport of thickness and tracers
Expand Down Expand Up @@ -413,8 +415,6 @@ subroutine li_advection_thickness_tracers(&
(trim(config_tracer_advection) .eq. 'fct' ) ) then
allocate(tend(nTracers,nVertLevels,nCells+1))
tend(:,:,:) = 0.0_RKIND
allocate(activeTracerHorizontalAdvectionEdgeFlux(nTracers,nVertLevels+1,nEdges+1))
activeTracerHorizontalAdvectionEdgeFlux(:,:,:) = 0.0_RKIND
endif

! Transport thickness and tracers
Expand Down Expand Up @@ -481,25 +481,25 @@ subroutine li_advection_thickness_tracers(&
call li_tracer_advection_fct_tend(&
tend, advectedTracers, layerThicknessOld, &
layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, &
nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)!{{{
nTracers, layerThicknessEdgeFlux, computeBudgets=.false.)!{{{
elseif (trim(config_thickness_advection) .eq. 'fct') then
! Call fct routine for thickness first, and use activeTracerHorizontalAdvectionEdgeFlux
! Call fct routine for thickness first, and use layerThicknessEdgeFlux
! returned by that call as normalThicknessFlux for call to tracer fct
call li_tracer_advection_fct_tend(&
tend(nTracers:,:,:), advectedTracers(nTracers:,:,:), layerThicknessOld * 0.0_RKIND + 1.0_RKIND, &
layerNormalVelocity, 0.0_RKIND * normalVelocity, dt, &
1, activeTracerHorizontalAdvectionEdgeFlux(nTracers:,:,:), computeBudgets=.false.)
1, layerThicknessEdgeFlux, computeBudgets=.false.)
! layerThickness is last tracer. However, for some reason
! this: layerThickness(:,:) = advectedTracers(nTracers,:,:) does not conserve mass!
! This does conserve mass:
layerThickness(:,:) = layerThickness(:,:) + tend(nTracers,:,:) * dt

if (trim(config_tracer_advection) .eq. 'fct') then
! Call fct for tracers, using activeTracerHorizontalAdvectionEdgeFlux
! Call fct for tracers, using layerThicknessEdgeFlux
! from fct thickness advection as normalThicknessFlux
call li_tracer_advection_fct_tend(&
tend(1:nTracers-1,:,:), advectedTracers(1:nTracers-1,:,:), layerThicknessOld, &
activeTracerHorizontalAdvectionEdgeFlux(nTracers,:,:), 0.0_RKIND * normalVelocity, dt, &
layerThicknessEdgeFlux, 0.0_RKIND * normalVelocity, dt, &
nTracers-1, computeBudgets=.false.)
elseif (trim(config_tracer_advection) .eq. 'none') then
! do nothing
Expand Down Expand Up @@ -642,10 +642,16 @@ subroutine li_advection_thickness_tracers(&
GLfluxSign = -1.0_RKIND
theGroundedCell = iCell2
endif
do k = 1, nVertLevels
thicknessFluxEdge = layerNormalVelocity(k,iEdge) * dvEdge(iEdge) * layerThicknessEdge(k,iEdge)
fluxAcrossGroundingLine(iEdge) = fluxAcrossGroundingLine(iEdge) + GLfluxSign * thicknessFluxEdge / dvEdge(iEdge)
enddo
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
Expand Down Expand Up @@ -714,8 +720,7 @@ subroutine li_advection_thickness_tracers(&
advCoefs3rd, &
advMaskHighOrder, &
advMask2ndOrder, &
tend, &
activeTracerHorizontalAdvectionEdgeFlux)
tend)
endif
! clean up
Expand Down
15 changes: 8 additions & 7 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ subroutine li_tracer_advection_fct_tend( &
tend, tracers, layerThickness, &
normalThicknessFlux, w, dt, &
nTracers, &
activeTracerHorizontalAdvectionEdgeFlux, &
layerThicknessEdgeFlux, &
computeBudgets)!{{{
use li_mesh
!-----------------------------------------------------------------
Expand All @@ -74,8 +74,8 @@ subroutine li_tracer_advection_fct_tend( &

real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
tend !< [inout] Tracer tendency to which advection added
real (kind=RKIND), dimension(:,:,:), intent(inout), optional :: &
activeTracerHorizontalAdvectionEdgeFlux !< [inout] used to compute higher order normalThicknessFlux
real (kind=RKIND), dimension(:,:), intent(inout), optional :: &
layerThicknessEdgeFlux !< [inout] used to compute higher order normalThicknessFlux
!-----------------------------------------------------------------
! Input parameters
!-----------------------------------------------------------------
Expand Down Expand Up @@ -462,18 +462,19 @@ subroutine li_tracer_advection_fct_tend( &
#endif
! Compute budget and monotonicity diagnostics if needed
! Use activeTracerHorizontalAdvectionEdgeFlux from the call to fct for thickness
! Use layerThicknessEdgeFlux from the call to fct for thickness
! advection as the higher-order normalThicknessFlux for fct tracer advection.
if (present(activeTracerHorizontalAdvectionEdgeFlux)) then
if (present(layerThicknessEdgeFlux)) then
do iEdge = 1,nEdges
do k = 1, nVertLevels
! Save u*h*T flux on edge for analysis. This variable will be
! divided by h at the end of the time step.
! average normal velocities from layer interfaces to layer midpoints
if (dvEdge(iEdge) > 0.0_RKIND) then
activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = &
layerThicknessEdgeFlux(k,iEdge) = &
(lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge)
else
activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = 0.0_RKIND
layerThicknessEdgeFlux(k,iEdge) = 0.0_RKIND
endif
enddo
enddo
Expand Down

0 comments on commit 1343d1d

Please sign in to comment.