From 709919e1b8198c8c07f4a9e68531c5d5299bf9b3 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 20 Sep 2023 20:22:21 -0600 Subject: [PATCH] Move surface and basal mass balance to separate subroutine 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. --- .../src/mode_forward/mpas_li_advection.F | 276 +++++++++++++----- .../mpas_li_time_integration_fe_rk.F | 46 ++- 2 files changed, 243 insertions(+), 79 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 index 0a0f7036daa2..bfa9c3c820f5 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F @@ -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 @@ -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(& @@ -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 @@ -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 @@ -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 @@ -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 @@ -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.) @@ -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 !*********************************************************************** diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F index 17714dfb6ffc..fcbd798cf076 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F @@ -82,6 +82,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) use li_subglacial_hydro use li_velocity use li_bedtopo + use li_mask !----------------------------------------------------------------- ! input variables @@ -104,12 +105,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) !----------------------------------------------------------------- type (block_type), pointer :: block integer :: err_tmp - type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool + type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool, velocityPool logical, pointer :: config_restore_calving_front logical, pointer :: config_calculate_damage logical, pointer :: config_finalize_damage_after_advection - character (len=StrKIND), pointer :: config_thickness_advection + logical :: advectTracers + character (len=StrKIND), pointer :: config_thickness_advection, config_tracer_advection character (len=StrKIND), pointer :: config_thermal_solver character (len=StrKIND), pointer :: config_time_integration integer, pointer :: config_rk_order, config_rk3_stages @@ -118,7 +120,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) temperature, & enthalpy, & waterFrac - real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d + real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d, dynamicThickening real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & @@ -143,6 +145,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_config(liConfigs, 'config_calculate_damage',config_calculate_damage) call mpas_pool_get_config(liConfigs, 'config_finalize_damage_after_advection', config_finalize_damage_after_advection) call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) + call mpas_pool_get_config(liConfigs, 'config_tracer_advection', config_tracer_advection) call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages) @@ -151,6 +154,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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, 'mesh', meshPool) + call mpas_pool_get_subpool(domain % blocklist % structs, 'velocity', velocityPool) call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) call mpas_pool_get_array(meshPool, 'deltat', deltat) @@ -223,6 +227,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! *** TODO: Should basal melt rate calculation, column physics, and hydrology go inside RK loop? *** call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) call mpas_pool_get_array(geometryPool, 'damage', damage) @@ -267,8 +272,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSSPweights(4) = 0.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 3) ) then + ! Standard 3-stage Strong Stability Preserving RK3 if (config_rk3_stages == 3) then - ! use Strong Stability Preserving RK3 nRKstages = 3 rkSubstepWeights(:) = 1.0_RKIND @@ -298,8 +303,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) endif ! *** Start RK loop *** do rkStage = 1, nRKstages - call mpas_log_write('beginning rk step $i', & - intArgs=(/rkStage/)) + call mpas_log_write('beginning rk stage $i of $i', & + intArgs=(/rkStage, nRKstages/)) deltat = deltatFull * rkSubstepWeights(rkStage) ! === calculate damage =========== @@ -390,8 +395,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo - endif ! if config_rk_order == 3 - + endif ! if 1 < rkStage < 4 + ! === Ensure damage is within bounds after full time integration === if ( (rkStage .eq. nRKstages) .and. (config_finalize_damage_after_advection) ) then call mpas_timer_start("finalize damage") @@ -408,8 +413,31 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! *** end RK loop *** enddo -! Reset time step to full length after RK loop + ! Reset time step to full length after RK loop deltat = deltatFull + + ! Calculate dynamicThickening + dynamicThickening = (sum(layerThickness, 1) - sum(layerThicknessPrev,1)) / deltat * scyr ! units of m/yr + + + ! 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) + +! === Surface and basal mass balance === +! TODO: Determine whether halo updates are needed prior to calling li_surface_basal_mass_bal + call mpas_timer_start("smb_bmb") + if ( trim(config_tracer_advection) .eq. 'none' ) then + advectTracers = .false. + else + advectTracers = .true. + endif + + call li_surface_basal_mass_bal(domain, err_tmp, advectTracers) + err = ior(err, err_tmp) + call mpas_timer_stop("smb_bmb") + ! === Calve ice ======================== call mpas_timer_start("calve_ice")