Skip to content

Commit

Permalink
Move surface and basal mass balance to separate subroutine
Browse files Browse the repository at this point in the history 
Move surface and basal mass balance out of advection subroutine
and into a separate subroutine, li_surface_basal_mass_bal, and call
this after the RK loop.
  • Loading branch information
@trhille
trhille committed Sep 21, 2023
1 parent 4f69cfe commit 709919e
Show file tree
Hide file tree
Showing 2 changed files with 243 additions and 79 deletions.
276 changes: 206 additions & 70 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 @@ -38,13 +38,31 @@ module li_advection
! Public parameters
!
!--------------------------------------------------------------------
real (kind=RKIND), dimension(:,:,:), pointer, public :: &
advectedTracers, & ! values of advected tracers
advectedTracersOld ! old values of advected tracers

real (kind=RKIND), dimension(:,:), pointer, public :: &
surfaceTracers, & ! tracer values for new ice at upper surface
basalTracers ! tracer values for new ice at lower surface

type (field3DReal), pointer, public :: &
advectedTracersField, & ! scratch field containing values of advected tracers
advectedTracersOldField ! scratch field containing old values of advected tracers

type (field2DReal), pointer, public :: &
layerThicknessOldField ! scratch field containing old values of layer thickness

type (field2DReal), pointer, public :: &
surfaceTracersField, & ! scratch field containing values of surface tracers
basalTracersField ! scratch field containing values of basal tracers
!--------------------------------------------------------------------
!
! Public member functions
!
!--------------------------------------------------------------------
public :: li_advection_thickness_tracers, &
li_surface_basal_mass_bal, &
li_layer_normal_velocity, &
li_update_geometry

Expand All @@ -65,12 +83,8 @@ module li_advection
!> \author William Lipscomb, Trevor Hillebrand
!> \date November 2015, updated Sept 2023 to include FCT
!> \details
!> This routine (1) computes new values of ice thickness and tracers under
!> horizontal advection, (2) applies surface and basal mass balance
!> terms, and (3) remaps tracers back onto the standard vertical grid.
!> Based on subroutine li_tendency_thickness by Matthew Hoffman
!> Note: The thickness and tracer fields must have had halo updates
!> before calling this subroutine.
!> This routinecomputes new values of ice thickness and tracers under
!> horizontal advection.
!-----------------------------------------------------------------------

subroutine li_advection_thickness_tracers(&
Expand Down Expand Up @@ -168,29 +182,6 @@ subroutine li_advection_thickness_tracers(&
layerThicknessEdge, & ! layer thickness on upstream edge of cell
normalVelocity ! horizontal velocity on interfaces

real (kind=RKIND), dimension(:,:,:), pointer :: &
advectedTracers, & ! values of advected tracers
advectedTracersOld ! old values of advected tracers

real (kind=RKIND), dimension(:,:), pointer :: &
surfaceTracers, & ! tracer values for new ice at upper surface
basalTracers ! tracer values for new ice at lower surface

type (field3DReal), pointer :: &
advectedTracersField, & ! scratch field containing values of advected tracers
advectedTracersOldField ! scratch field containing old values of advected tracers

type (field2DReal), pointer :: &
layerThicknessOldField ! scratch field containing old values of layer thickness

type (field2DReal), pointer :: &
surfaceTracersField, & ! scratch field containing values of surface tracers
basalTracersField ! scratch field containing values of basal tracers

type (field1DInteger), pointer :: &
cellMaskTemporaryField, & ! scratch field containing old values of cellMask
thermalCellMaskField

! Allocatable arrays need for flux-corrected transport advection
real (kind=RKIND), dimension(:,:,:), allocatable :: tend
real (kind=RKIND), dimension(:,:,:), allocatable :: activeTracerHorizontalAdvectionEdgeFlux
Expand Down Expand Up @@ -349,13 +340,6 @@ subroutine li_advection_thickness_tracers(&
call mpas_allocate_scratch_field(basalTracersField, .true.)
basalTracers => basalTracersField % array

call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField)
call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.)

call mpas_pool_get_field(scratchPool, 'iceCellMask', thermalCellMaskField)
call mpas_allocate_scratch_field(thermalCellMaskField, .true.)
thermalCellMask => thermalCellMaskField % array

! define max number of advection cells as max number of edges*2
nAdvCellsMax = maxEdges2

Expand All @@ -367,7 +351,7 @@ subroutine li_advection_thickness_tracers(&
err = ior(err, err_tmp)

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

!-----------------------------------------------------------------
! Horizontal transport of thickness and tracers
Expand Down Expand Up @@ -528,17 +512,194 @@ subroutine li_advection_thickness_tracers(&
enddo
endif

! Deallocate arrays for fct
if ( (trim(config_thickness_advection) .eq. 'fct') .or. &
(trim(config_tracer_advection) .eq. 'fct') ) then
deallocate( nAdvCellsForEdge, &
advCellsForEdge, &
advCoefs, &
advCoefs3rd, &
advMaskHighOrder, &
advMask2ndOrder, &
tend, &
activeTracerHorizontalAdvectionEdgeFlux)
endif

! 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
elseif (trim(config_thickness_advection) .eq. 'none') then

! Do nothing

! Update the thickness and cellMask before applying the mass balance.
! The update is needed because the SMB and BMB depend on whether ice is present.
thickness = sum(layerThickness, 1)
call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
else

call mpas_log_write(trim(config_thickness_advection) // ' is not a valid option for thickness/tracer advection.', &
MPAS_LOG_ERR)
err_tmp = 1

end if ! config_thickness_advection

err = ior(err,err_tmp)

! === error check
if (err > 0) then
call mpas_log_write("An error has occurred in li_advection_thickness_tracers", MPAS_LOG_ERR)
endif

!--------------------------------------------------------------------
end subroutine li_advection_thickness_tracers

!***********************************************************************
!
! subroutine li_surface_basal_mass_bal
!
!> \brief Applies surface and basal mass balances after advection
!> \author William Lipscomb, Trevor Hillebrand
!> \date November 2015, moved to separate subroutine Sept 2023
!> \details
!> Applies surface and basal mass balance terms and remaps tracers
!> back onto the standard vertical grid.
!> Based on subroutine li_tendency_thickness by Matthew Hoffman
!> Note: The thickness and tracer fields must have had halo updates
!> before calling this subroutine.
!-----------------------------------------------------------------------

subroutine li_surface_basal_mass_bal(domain, err, advectTracersIn)

logical, intent(in), optional :: &
advectTracersIn !< Input: if true, then advect tracers as well as thickness

type (domain_type), intent(inout) :: &
domain !< Input/Output: domain object

integer, intent(out) :: &
err !< Output: error flag

type (mpas_pool_type), pointer :: &
meshPool, & !< Input: mesh information
velocityPool , & !< Input/output: velocity information
geometryPool, & !< Input/output: geometry information to be updated
thermalPool, & !< Input/output: temperature/enthalpy information to be updated
scratchPool

real (kind=RKIND), dimension(:), pointer :: &
thickness, &
bedTopography, &
sfcMassBal, &
sfcMassBalApplied, &
groundedSfcMassBalApplied, &
basalMassBal, &
basalMassBalApplied, &
groundedBasalMassBal, &
floatingBasalMassBal, &
groundedBasalMassBalApplied, &
floatingBasalMassBalApplied, &
groundedToFloatingThickness, & ! thickness changing from grounded to floating or vice versa
fluxAcrossGroundingLine, & ! magnitude of flux across GL
fluxAcrossGroundingLineOnCells ! magnitude of flux across GL, calculated on cells

real (kind=RKIND), dimension(:,:), pointer :: &
layerThickness
real (kind=RKIND), dimension(:,:), pointer :: &
layerNormalVelocity, & ! normal component of velocity on cell edges at layer midpoints
layerThicknessEdge ! layer thickness on upstream edge of cell

real (kind=RKIND), dimension(:,:), pointer :: &
temperature, & ! interior ice temperature
waterFrac, & ! interior water fraction
drainedInternalMeltRate, & ! excess internal melt drained to the bed
enthalpy ! interior ice enthalpy

real(kind=RKIND) :: GLfluxSign, thicknessFluxEdge
integer, dimension(:), pointer :: indexToCellID, indexToEdgeID
integer, pointer :: nCells, nEdges, nVertLevels
integer, pointer :: config_stats_cell_ID
integer :: iEdge, iEdgeOnCell, iCell
integer :: err_tmp
integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnCell
integer :: iCell1, iCell2, theGroundedCell, k
logical :: advectTracers

logical, pointer :: config_thermal_calculate_bmb
real (kind=RKIND), pointer :: &
config_thermal_thickness, &
config_ice_density, &
dt

type (field1DInteger), pointer :: &
cellMaskTemporaryField, & ! scratch field containing old values of cellMask
thermalCellMaskField

integer, dimension(:), pointer :: cellMask, thermalCellMask, edgeMask

call mpas_pool_get_subpool(domain % blocklist % structs, 'mesh', meshPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'thermal', thermalPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool)
call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool)

! get dimensions
call mpas_pool_get_dimension(meshPool, 'nCells', nCells)

! get arrays from the geometry pool
call mpas_pool_get_array(geometryPool, 'thickness', thickness)
call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography)
call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal)
call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied)
call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied)
call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal)
call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal)
call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal)
call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied)
call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied)
call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied)
call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness)
call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)

! get arrays from the velocity pool
call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity)
call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine)
call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells)

! get arrays from the thermal pool
call mpas_pool_get_array(thermalPool, 'temperature', temperature)
call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac)
call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy)
call mpas_pool_get_array(thermalPool, 'drainedInternalMeltRate', drainedInternalMeltRate)

! get config variables
call mpas_pool_get_config(liConfigs, 'config_thermal_calculate_bmb', config_thermal_calculate_bmb)
call mpas_pool_get_config(liConfigs, 'config_thermal_thickness', config_thermal_thickness)
call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density)

call mpas_pool_get_dimension(meshPool, 'nCells', nCells)
call mpas_pool_get_array(meshPool, 'deltat', dt)
call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges)
call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels)
call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell)
call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell)
call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID)
call mpas_pool_get_array(meshPool, 'indexToEdgeID', indexToEdgeID)
call mpas_pool_get_config(liConfigs, 'config_stats_cell_ID', config_stats_cell_ID)
call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge)

call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField)
call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.)

call mpas_pool_get_field(scratchPool, 'iceCellMask', thermalCellMaskField)
call mpas_allocate_scratch_field(thermalCellMaskField, .true.)
thermalCellMask => thermalCellMaskField % array

if (present(advectTracersIn)) then
advectTracers = advectTracersIn
else
advectTracers = .false.
endif
! Set up tracers again using post-advection fields, or make tracer variables public?
! First trying with public variables

!-----------------------------------------------------------------
! Add the surface and basal mass balance to the layer thickness
Expand Down Expand Up @@ -691,33 +852,8 @@ subroutine li_advection_thickness_tracers(&
advectedTracers)
endif
elseif (trim(config_thickness_advection) .eq. 'none') then
! Do nothing
else
call mpas_log_write(trim(config_thickness_advection) // ' is not a valid option for thickness/tracer advection.', &
MPAS_LOG_ERR)
err_tmp = 1
end if ! config_thickness_advection
err = ior(err,err_tmp)
! Deallocate arrays for fct
if ( (trim(config_thickness_advection) .eq. 'fct') .or. &
(trim(config_tracer_advection) .eq. 'fct') ) then
deallocate( nAdvCellsForEdge, &
advCellsForEdge, &
advCoefs, &
advCoefs3rd, &
advMaskHighOrder, &
advMask2ndOrder, &
tend, &
activeTracerHorizontalAdvectionEdgeFlux)
endif
! clean up
call mpas_deallocate_scratch_field(advectedTracersField, .true.)
call mpas_deallocate_scratch_field(advectedTracersOldField, .true.)
Expand All @@ -729,11 +865,11 @@ subroutine li_advection_thickness_tracers(&
! === error check
if (err > 0) then
call mpas_log_write("An error has occurred in li_advection_thickness_tracers", MPAS_LOG_ERR)
call mpas_log_write("An error has occurred in li_surface_basal_mass_bal", MPAS_LOG_ERR)
endif
!--------------------------------------------------------------------
end subroutine li_advection_thickness_tracers
end subroutine li_surface_basal_mass_bal
!***********************************************************************
Expand Down
Loading

0 comments on commit 709919e

Please sign in to comment.