From 4ec4101c73b8659d8a040f1a5e2f09849fa865e5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 30 Aug 2023 19:16:29 -0700 Subject: [PATCH 01/36] Add RK2 loop in time integration module Add RK2 loop in time integration module. This is the bare bones of RK2 time stepping. Compiles but likely will not run. --- .../mpas-albany-landice/src/Registry.xml | 3 + .../src/mode_forward/mpas_li_calving.F | 13 +- .../mpas_li_time_integration_fe.F | 179 +++++++++++++++--- 3 files changed, 161 insertions(+), 34 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 22a863096ad4..671e260ca9ea 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1183,6 +1183,9 @@ is the value of that variable from the *previous* time level! + diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 36302a1a78c3..3fbc14c04a65 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -3314,7 +3314,7 @@ end subroutine damage_calving !----------------------------------------------------------------------- - subroutine li_calculate_damage(domain, err) + subroutine li_calculate_damage(domain, rkTend, err) !----------------------------------------------------------------- ! input variables @@ -3330,7 +3330,7 @@ subroutine li_calculate_damage(domain, err) ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - + real (kind=RKIND), dimension(:,:), intent(out) :: rkTend !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -3368,8 +3368,8 @@ subroutine li_calculate_damage(domain, err) real (kind=RKIND), pointer :: deltat !< time step (s) integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell - integer, pointer :: nCells - integer :: iCell, jCell, iNeighbor, n_damage_downstream + integer, pointer :: nCells, nVertLevels + integer :: iCell, jCell, iNeighbor, n_damage_downstream, k real(kind=RKIND) :: damage_downstream real(kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY real(kind=RKIND), dimension(6) :: localMinInfo, localMaxInfo, globalMinInfo, globalMaxInfo @@ -3396,6 +3396,7 @@ subroutine li_calculate_damage(domain, err) ! get fields call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_array(meshPool, 'deltat', deltat) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) @@ -3455,6 +3456,10 @@ subroutine li_calculate_damage(domain, err) ! damageSource(:) * (1+exp(-uReconstructX(1,:)*seconds/1000.0))) ddamagedt(:) = damageSource(:) * damage(:) + ! Update damage tendency for RK time integration + do k = 1, nVertLevels + rkTend(k,:) = rkTend(k,:) + ddamagedt(:) + enddo damage(:) = damage(:) + ddamagedt(:) * deltat diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index 0b3cd8427fd7..ca16d208a7ee 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -30,11 +30,12 @@ module li_time_integration_fe use li_advection use li_calving, only: li_calve_ice, li_restore_calving_front, li_calculate_damage, li_finalize_damage_after_advection - use li_thermal, only: li_thermal_solver + use li_thermal, only: li_thermal_solver, li_enthalpy_to_temperature_kelvin use li_iceshelf_melt use li_diagnostic_vars use li_setup use li_constants + use li_mesh implicit none private @@ -103,17 +104,57 @@ subroutine li_time_integrator_forwardeuler(domain, err) !----------------------------------------------------------------- type (block_type), pointer :: block integer :: err_tmp - + type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool + 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 + character (len=StrKIND), pointer :: config_thermal_solver + integer :: rkStep, iCell, iTracer, k + real (kind=RKIND), dimension(:,:,:), pointer :: rkTend + real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & + temperature, & + enthalpy, & + waterFrac + real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d + real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions + real (kind=RKIND), dimension(:), allocatable :: damageProv, passiveTracer2dProv + real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessProv, & + temperatureProv, & + enthalpyProv + integer, pointer :: nVertLevels + integer, pointer :: nCells + real (kind=RKIND), pointer :: deltat err = 0 err_tmp = 0 call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) 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_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_thermal_solver', config_thermal_solver) + + 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_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + call mpas_pool_get_array(meshPool, 'layerCenterSigma', layerCenterSigma) + + allocate(temperatureProv(nVertLevels, nCells+1)) + allocate(enthalpyProv(nVertLevels, nCells+1)) + allocate(layerThicknessProv(nVertLevels, nCells+1)) + allocate(damageProv(nCells+1)) + allocate(passiveTracer2dProv(nCells+1)) + + temperatureProv(:,:) = 0.0_RKIND + enthalpyProv(:,:) = 0.0_RKIND + layerThicknessProv(:,:) = 0.0_RKIND + damageProv(:) = 0.0_RKIND + passiveTracer2dProv(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -135,6 +176,8 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("face melting for grounded ice") +! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? *** + ! === Basal melting for floating ice =========== call mpas_timer_start("basal melting for floating ice") call li_basal_melt_floating_ice(domain, err_tmp) @@ -147,38 +190,113 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("vertical therm") -! === calculate damage =========== - if (config_calculate_damage) then - call mpas_timer_start("damage") - call li_calculate_damage(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("damage") - endif - -! === Compute new state for prognostic variables ================================== - call mpas_timer_start("advect thickness and tracers") - call advection_solver(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("advect thickness and tracers") - -! === finalize damage after advection =========== - if (config_finalize_damage_after_advection) then - call mpas_timer_start("finalize damage") - call li_finalize_damage_after_advection(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("finalize damage") - endif - ! === Update subglacial hydrology =========== ! It's not clear where the best place to call this should be. ! Seems sensible to put it after thermal evolution is complete to get updated basal melting source term. ! Also seems (might be?) better to put it after geometry evolution. ! We want it before the velocity solve since the hydro model can control the velo basal b.c. +! This was moved to be before geometry evolution to simplify RK time integration. call mpas_timer_start("subglacial hydro") call li_SGH_solve(domain, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("subglacial hydro") +! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? *** + call mpas_pool_get_array(geometryPool, 'rkTend', rkTend) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) + call mpas_pool_get_array(geometryPool, 'damage', damage) + + call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) + call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + ! Save relevant fields before RK loop, to be used in update at the end + temperatureProv(:,:) = temperature(:,:) + enthalpyProv(:,:) = enthalpy(:,:) + layerThicknessProv(:,:) = layerThickness(:,:) + damageProv(:) = damage(:) + passiveTracer2dProv(:) = passiveTracer2d(:) + rkTend(:,:,:) = 0.0_RKIND +! *** Start RK loop *** + do rkStep = 1, 2 + ! === calculate damage =========== + if (config_calculate_damage) then + call mpas_timer_start("damage") + call li_calculate_damage(domain, rkTend(2,:,:), err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("damage") + endif + + ! === Compute new state for prognostic variables === + call mpas_timer_start("advect thickness and tracers") + call advection_solver(domain, rkTend, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advect thickness and tracers") + + ! === finalize damage after advection =========== + if (config_finalize_damage_after_advection) then + call mpas_timer_start("finalize damage") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") + endif + + ! Update thickness, damage, temperature or enthalpy, passiveTracer2d + ! before final velocity solve + if (rkStep .eq. 2) then ! or whichever is last rk step + ! y_{n+1} = y_n + (k1 + k2) / 2 + if (trim(config_thermal_solver) == 'enthalpy') then + enthalpy(:,:) = enthalpyProv(:,:) + deltat * rkTend(1,:,:) / 2.0_RKIND + ! given the enthalpy, compute the temperature and water fraction + do iCell = 1, nCells + + call li_enthalpy_to_temperature_kelvin(& + layerCenterSigma, & + thickness(iCell), & + enthalpy(:,iCell), & + temperature(:,iCell), & + waterFrac(:,iCell)) + + enddo + elseif (trim(config_thermal_solver) == 'temperature') then + temperature(:,:) = temperatureProv(:,:) + deltat * rkTend(1,:,:) / 2.0_RKIND + endif + + if (config_calculate_damage) then + do iCell = 1, nCells + damage(iCell) = damageProv(iCell) + deltat * sum(rkTend(2,:,iCell) & + * layerThicknessFractions) / 2.0_RKIND + enddo + endif + + do iCell = 1, nCells + passiveTracer2d(iCell) = passiveTracer2dProv(iCell) + deltat * & + sum(rkTend(3,:,iCell) * layerThicknessFractions) / 2.0_RKIND + enddo + + if (trim(config_thickness_advection) == 'fct') then + layerThickness(:,:) = layerThicknessProv(:,:) + deltat * rkTend(4,:,:) / 2.0_RKIND + thickness = sum(layerThickness, 1) + endif + + ! === Ensure damage is within bounds after full time integration === + if (config_finalize_damage_after_advection) then + call mpas_timer_start("finalize damage") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") + endif + endif + + ! Update velocity for each RK step + ! === Solve Velocity ===================== + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + +! *** end RK loop *** + enddo + ! === Calve ice ======================== call mpas_timer_start("calve_ice") @@ -207,10 +325,11 @@ subroutine li_time_integrator_forwardeuler(domain, err) call li_bedtopo_solve(domain, err=err_tmp) err = ior(err, err_tmp) +! TODO: Determine whether we need to update velocities again ! === Solve Velocity ===================== - ! During time-stepping, we always solveVelo - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) - err = ior(err, err_tmp) +! ! During time-stepping, we always solveVelo +! call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) +! err = ior(err, err_tmp) ! === Calculate diagnostic variables for new state ===================== @@ -591,7 +710,7 @@ end subroutine prepare_advection !> in calculate tendencies are now in prepare_advection. !----------------------------------------------------------------------- - subroutine advection_solver(domain, err) + subroutine advection_solver(domain, rkTend, err) use mpas_timekeeping use li_mask @@ -604,7 +723,7 @@ subroutine advection_solver(domain, err) ! input/output variables !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - + real (kind=RKIND), dimension(:,:,:), intent(inout) :: rkTend !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- From bf947b6ccb4305c6f7a11978cabf243219dd743e Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 30 Aug 2023 20:17:54 -0700 Subject: [PATCH 02/36] Pass rkTend to advection routine Pass rkTend to advection routine to accumulate tendencies --- .../src/mode_forward/mpas_li_advection.F | 38 +++++++++++++++---- .../mpas_li_time_integration_fe.F | 2 + 2 files changed, 33 insertions(+), 7 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 7d93f114b797..32d58d4af6d7 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 @@ -81,6 +81,7 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & scratchPool, & + rkTend, & err, & advectTracersIn) @@ -110,6 +111,9 @@ subroutine li_advection_thickness_tracers(& !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + rkTend !< Input/output: accumulated tendencies for RK time integration + type (mpas_pool_type), intent(inout) :: & velocityPool !< Input/output: velocity information ! (needs to be inout for li_calculate_mask call @@ -207,12 +211,14 @@ subroutine li_advection_thickness_tracers(& real (kind=RKIND), dimension(:), pointer :: dvEdge character (len=StrKIND), pointer :: & - config_thickness_advection ! method for advecting thickness and tracers + config_thickness_advection, & ! method for advecting thickness and tracers + config_thermal_solver logical, pointer :: & config_print_thickness_advection_info, & !TODO - change to config_print_advection_info? config_print_thermal_info, & - config_thermal_calculate_bmb + config_thermal_calculate_bmb, & + config_calculate_damage real (kind=RKIND), pointer :: & config_ice_density, & ! rhoi @@ -225,7 +231,7 @@ subroutine li_advection_thickness_tracers(& logical :: advectTracers ! if true, then advect tracers as well as thickness - integer :: iCell1, iCell2, theGroundedCell + integer :: iCell1, iCell2, theGroundedCell, iTracer integer, parameter :: horizAdvOrder = 4 integer, parameter :: vertAdvOrder = 1 @@ -315,6 +321,8 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) call mpas_pool_get_config(liConfigs, 'config_num_halos', config_num_halos) + call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) + call mpas_pool_get_config(liConfigs, 'config_calculate_damage', config_calculate_damage) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) @@ -514,11 +522,27 @@ subroutine li_advection_thickness_tracers(& call mpas_log_write(trim(config_tracer_advection) // & ' tracer advection is not currently supported with fct thickness advection.', MPAS_LOG_ERR) endif + + ! Update tendencies for RK integration + iTracer = 0 + if ( (trim(config_thermal_solver) == 'enthalpy') .or. (trim(config_thermal_solver) == 'temperature') ) then + iTracer = iTracer + 1 + rkTend(1,:,:) = rkTend(1,:,:) + tend(iTracer,:,:) + endif + if (config_calculate_damage) then + iTracer = iTracer + 1 + rkTend(2,:,:) = rkTend(2,:,:) + tend(iTracer,:,:) + endif + iTracer = iTracer + 1 + rkTend(3,:,:) = rkTend(3,:,:) + tend(iTracer,:,:) ! passiveTracer2d + iTracer = iTracer + 1 + rkTend(4,:,:) = rkTend(4,:,:) + tend(iTracer,:,:) ! layerThickness + else - err_tmp = 1 - call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & - ", config_tracer_advection = " // trim(config_tracer_advection) // & - " is not a supported combination.", MPAS_LOG_ERR) + err_tmp = 1 + call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & + ", config_tracer_advection = " // trim(config_tracer_advection) // & + " is not a supported combination.", MPAS_LOG_ERR) endif if (config_print_thickness_advection_info) then diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index ca16d208a7ee..c36c49186c50 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -841,6 +841,7 @@ subroutine advection_solver(domain, rkTend, err) geometryPool, & thermalPool, & scratchPool, & + rkTend, & err_tmp, & advectTracersIn = .true.) @@ -860,6 +861,7 @@ subroutine advection_solver(domain, rkTend, err) geometryPool, & thermalPool, & scratchPool, & + rkTend, & err_tmp, & advectTracersIn = .false.) From c703988d4a83d55d468a0a03578ec8c6f29f8755 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 31 Aug 2023 07:59:15 -0700 Subject: [PATCH 03/36] Finalize thickness and tracers at end of RK loop --- .../src/mode_forward/mpas_li_advection.F | 1 + .../src/mode_forward/mpas_li_calving.F | 4 +- .../mpas_li_time_integration_fe.F | 79 ++++++++++++++----- 3 files changed, 64 insertions(+), 20 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 32d58d4af6d7..d1098bc52f75 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 @@ -524,6 +524,7 @@ subroutine li_advection_thickness_tracers(& endif ! Update tendencies for RK integration + ! How do we deal with these? should they be divided by updated layerThickness now? iTracer = 0 if ( (trim(config_thermal_solver) == 'enthalpy') .or. (trim(config_thermal_solver) == 'temperature') ) then iTracer = iTracer + 1 diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 3fbc14c04a65..99ee7a01e39e 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -3348,6 +3348,7 @@ subroutine li_calculate_damage(domain, rkTend, err) real(kind=RKIND), pointer :: config_flowLawExponent logical, pointer :: config_print_calving_info + real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: eMax real (kind=RKIND), dimension(:), pointer :: tauMax, tauMin @@ -3402,6 +3403,7 @@ subroutine li_calculate_damage(domain, rkTend, err) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'damage', damage) call mpas_pool_get_array(geometryPool, 'ddamagedt', ddamagedt) call mpas_pool_get_array(geometryPool, 's0', s0) @@ -3458,7 +3460,7 @@ subroutine li_calculate_damage(domain, rkTend, err) ddamagedt(:) = damageSource(:) * damage(:) ! Update damage tendency for RK time integration do k = 1, nVertLevels - rkTend(k,:) = rkTend(k,:) + ddamagedt(:) + rkTend(k,:) = rkTend(k,:) + ddamagedt(:) * layerThickness(k,:) enddo damage(:) = damage(:) + ddamagedt(:) * deltat diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index c36c49186c50..8e6c2d233fc5 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -119,10 +119,13 @@ subroutine li_time_integrator_forwardeuler(domain, err) waterFrac real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions - real (kind=RKIND), dimension(:), allocatable :: damageProv, passiveTracer2dProv real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessProv, & temperatureProv, & - enthalpyProv + enthalpyProv, & + passiveTracer3d, & + passiveTracer3dProv, & + damage3d, & + damage3dProv integer, pointer :: nVertLevels integer, pointer :: nCells real (kind=RKIND), pointer :: deltat @@ -147,14 +150,18 @@ subroutine li_time_integrator_forwardeuler(domain, err) allocate(temperatureProv(nVertLevels, nCells+1)) allocate(enthalpyProv(nVertLevels, nCells+1)) allocate(layerThicknessProv(nVertLevels, nCells+1)) - allocate(damageProv(nCells+1)) - allocate(passiveTracer2dProv(nCells+1)) + allocate(damage3dProv(nVertLevels, nCells+1)) + allocate(damage3d(nVertLevels, nCells+1)) + allocate(passiveTracer3dProv(nVertLevels, nCells+1)) + allocate(passiveTracer3d(nVertLevels, nCells+1)) temperatureProv(:,:) = 0.0_RKIND enthalpyProv(:,:) = 0.0_RKIND layerThicknessProv(:,:) = 0.0_RKIND - damageProv(:) = 0.0_RKIND - passiveTracer2dProv(:) = 0.0_RKIND + damage3dProv(:,:) = 0.0_RKIND + damage3d(:,:) = 0.0_RKIND + passiveTracer3dProv(:,:) = 0.0_RKIND + passiveTracer3d(:,:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -215,8 +222,10 @@ subroutine li_time_integrator_forwardeuler(domain, err) temperatureProv(:,:) = temperature(:,:) enthalpyProv(:,:) = enthalpy(:,:) layerThicknessProv(:,:) = layerThickness(:,:) - damageProv(:) = damage(:) - passiveTracer2dProv(:) = passiveTracer2d(:) + do k = 1, nVertLevels + damage3dProv(k,:) = damage(:) + passiveTracer3dProv(k,:) = passiveTracer2d(:) + enddo rkTend(:,:,:) = 0.0_RKIND ! *** Start RK loop *** do rkStep = 1, 2 @@ -246,8 +255,19 @@ subroutine li_time_integrator_forwardeuler(domain, err) ! before final velocity solve if (rkStep .eq. 2) then ! or whichever is last rk step ! y_{n+1} = y_n + (k1 + k2) / 2 + + ! Need layerThickness first in order to update other tracers + if (trim(config_thickness_advection) == 'fct') then + layerThickness(:,:) = layerThicknessProv(:,:) + deltat * rkTend(4,:,:) / 2.0_RKIND + thickness = sum(layerThickness, 1) + endif + if (trim(config_thermal_solver) == 'enthalpy') then - enthalpy(:,:) = enthalpyProv(:,:) + deltat * rkTend(1,:,:) / 2.0_RKIND + where (layerThickness(:,:) > 0.0_RKIND) + enthalpy(:,:) = ( enthalpyProv(:,:) * layerThicknessProv(:,:) + & + deltat * rkTend(1,:,:) / 2.0_RKIND ) / layerThickness(:,:) + ! elsewhere, what? + end where ! given the enthalpy, compute the temperature and water fraction do iCell = 1, nCells @@ -259,27 +279,41 @@ subroutine li_time_integrator_forwardeuler(domain, err) waterFrac(:,iCell)) enddo + elseif (trim(config_thermal_solver) == 'temperature') then - temperature(:,:) = temperatureProv(:,:) + deltat * rkTend(1,:,:) / 2.0_RKIND + where (layerThickness(:,:) > 0.0_RKIND) + temperature(:,:) = ( temperatureProv(:,:) * layerThicknessProv(:,:) + & + deltat * rkTend(1,:,:) / 2.0_RKIND ) / layerThickness(:,:) + ! elsewhere, what? + end where endif if (config_calculate_damage) then do iCell = 1, nCells - damage(iCell) = damageProv(iCell) + deltat * sum(rkTend(2,:,iCell) & - * layerThicknessFractions) / 2.0_RKIND + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + damage3d(k,iCell) = ( damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + deltat * rkTend(2,k,iCell) / 2.0_RKIND ) / layerThickness(k,iCell) + else + damage3d(k,iCell) = 0.0_RKIND + endif + enddo + damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) enddo endif do iCell = 1, nCells - passiveTracer2d(iCell) = passiveTracer2dProv(iCell) + deltat * & - sum(rkTend(3,:,iCell) * layerThicknessFractions) / 2.0_RKIND + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + passiveTracer3d(k,iCell) = ( passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + deltat * rkTend(3,k,iCell) / 2.0_RKIND ) / layerThickness(k,iCell) + else + passiveTracer3d(k,iCell) = 0.0_RKIND + endif + enddo + passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo - if (trim(config_thickness_advection) == 'fct') then - layerThickness(:,:) = layerThicknessProv(:,:) + deltat * rkTend(4,:,:) / 2.0_RKIND - thickness = sum(layerThickness, 1) - endif - ! === Ensure damage is within bounds after full time integration === if (config_finalize_damage_after_advection) then call mpas_timer_start("finalize damage") @@ -344,6 +378,13 @@ subroutine li_time_integrator_forwardeuler(domain, err) call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler.", MPAS_LOG_ERR) endif + deallocate(temperatureProv) + deallocate(enthalpyProv) + deallocate(layerThicknessProv) + deallocate(damage3dProv) + deallocate(damage3d) + deallocate(passiveTracer3dProv) + deallocate(passiveTracer3d) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler From bf88eb89e9df6d6bbbba71641ad9fb8f4678e9ed Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 31 Aug 2023 09:52:35 -0700 Subject: [PATCH 04/36] Correct dimensions of rkTend Correct dimensions of rkTend, which was inadvertently 2D instead of 3D. --- components/mpas-albany-landice/src/Registry.xml | 2 +- .../src/mode_forward/mpas_li_time_integration_fe.F | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 671e260ca9ea..0a9880107647 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1183,7 +1183,7 @@ is the value of that variable from the *previous* time level! - Date: Thu, 31 Aug 2023 10:11:38 -0700 Subject: [PATCH 05/36] Update halos for rkTend in each step of RK loop Update halos for rkTend in each step of RK loop --- .../src/mode_forward/mpas_li_time_integration_fe.F | 2 ++ 1 file changed, 2 insertions(+) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index bc0578728e0a..90fc8ba9829a 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -245,6 +245,8 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advect thickness and tracers") + call mpas_dmpar_field_halo_exch(domain, 'rkTend') + ! === finalize damage after advection =========== if (config_finalize_damage_after_advection) then call mpas_timer_start("finalize damage") From 9ccdde68a97e16a83a1005c6bb61f82b5aebebf5 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 31 Aug 2023 10:44:00 -0700 Subject: [PATCH 06/36] Ensure damage >= damageNye after RK integration Ensure damage >= damageNye after RK integration. Previously, li_finalize_damage_after_advection did not impose the damageNye minumum. --- .../src/mode_forward/mpas_li_calving.F | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 99ee7a01e39e..500bdd87db7a 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -3683,14 +3683,18 @@ subroutine li_finalize_damage_after_advection(domain, err) endif ! put the damageMax value back to preserve the damage value (no heal) - where (damage < 0.0_RKIND) - damage = 0.0_RKIND - end where - - where (damage > 1.0_RKIND) - damage = 1.0_RKIND - end where - + ! limit damage to [damageNye, 1.0] + do iCell = 1, nCells + if (damage(iCell) < 0.0_RKIND) then + damage(iCell) = 0.0_RKIND + end if + if (damage(iCell) < damageNye(iCell)) then + damage(iCell) = damageNye(iCell) + end if + if (damage(iCell) > 1.0_RKIND) then + damage(iCell) = 1.0_RKIND + end if + end do do iCell = 1, nCells if (li_mask_is_grounded_ice(cellMask(iCell)) .or. .not. li_mask_is_ice(cellMask(iCell))) then From 5d1e5736a5e71db7f64bbdd64d38c20bdd9e7976 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Sun, 3 Sep 2023 11:04:01 -0700 Subject: [PATCH 07/36] Generalize RK time integration Generalize RK time integration to use order defined in namelist. Currently supports config_rk_order = 1, 2, 4. Third order to be added soon. --- .../mpas-albany-landice/src/Registry.xml | 6 +- .../src/mode_forward/mpas_li_advection.F | 9 ++- .../src/mode_forward/mpas_li_calving.F | 6 +- .../mode_forward/mpas_li_time_integration.F | 5 +- .../mpas_li_time_integration_fe.F | 68 ++++++++++++++++--- 5 files changed, 73 insertions(+), 21 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 0a9880107647..9555a7a2b592 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -486,8 +486,12 @@ /> + 0.0_RKIND) enthalpy(:,:) = ( enthalpyProv(:,:) * layerThicknessProv(:,:) + & - deltat * rkTend(1,:,:) / 2.0_RKIND ) / layerThickness(:,:) + deltatFull * rkTend(1,:,:) ) / layerThickness(:,:) ! elsewhere, what? end where ! given the enthalpy, compute the temperature and water fraction @@ -287,7 +329,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) elseif (trim(config_thermal_solver) == 'temperature') then where (layerThickness(:,:) > 0.0_RKIND) temperature(:,:) = ( temperatureProv(:,:) * layerThicknessProv(:,:) + & - deltat * rkTend(1,:,:) / 2.0_RKIND ) / layerThickness(:,:) + deltatFull * rkTend(1,:,:) ) / layerThickness(:,:) ! elsewhere, what? end where endif @@ -297,7 +339,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then damage3d(k,iCell) = ( damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltat * rkTend(2,k,iCell) / 2.0_RKIND ) / layerThickness(k,iCell) + deltatFull * rkTend(2,k,iCell) ) / layerThickness(k,iCell) else damage3d(k,iCell) = 0.0_RKIND endif @@ -310,7 +352,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then passiveTracer3d(k,iCell) = ( passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltat * rkTend(3,k,iCell) / 2.0_RKIND ) / layerThickness(k,iCell) + deltatFull * rkTend(3,k,iCell) ) / layerThickness(k,iCell) else passiveTracer3d(k,iCell) = 0.0_RKIND endif @@ -335,6 +377,8 @@ subroutine li_time_integrator_forwardeuler(domain, err) ! *** end RK loop *** enddo +! Reset time step to full length after RK loop + deltat = deltatFull ! === Calve ice ======================== call mpas_timer_start("calve_ice") @@ -755,7 +799,7 @@ end subroutine prepare_advection !> in calculate tendencies are now in prepare_advection. !----------------------------------------------------------------------- - subroutine advection_solver(domain, rkTend, err) + subroutine advection_solver(domain, rkTend, rkWeight, err) use mpas_timekeeping use li_mask @@ -763,7 +807,7 @@ subroutine advection_solver(domain, rkTend, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- - + real (kind=RKIND), intent(in) :: rkWeight !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- @@ -887,6 +931,7 @@ subroutine advection_solver(domain, rkTend, err) thermalPool, & scratchPool, & rkTend, & + rkWeight, & err_tmp, & advectTracersIn = .true.) @@ -907,6 +952,7 @@ subroutine advection_solver(domain, rkTend, err) thermalPool, & scratchPool, & rkTend, & + rkWeight, & err_tmp, & advectTracersIn = .false.) From 278a5d7836e086f9cbb0b2a80e2de92b4a75970d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 4 Sep 2023 11:08:53 -0700 Subject: [PATCH 08/36] Fix forward Euler time integration The previous Runge-Kutta commits left forward Euler time integration non-functioning. This makes FE time integration function again. --- .../src/mode_forward/mpas_li_advection.F | 1 - .../src/mode_forward/mpas_li_time_integration_fe.F | 11 ++++++----- 2 files changed, 6 insertions(+), 6 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 24eac8effbf9..88dcc0af9af5 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 @@ -527,7 +527,6 @@ subroutine li_advection_thickness_tracers(& endif ! Update tendencies for RK integration - ! How do we deal with these? should they be divided by updated layerThickness now? iTracer = 0 if ( (trim(config_thermal_solver) == 'enthalpy') .or. (trim(config_thermal_solver) == 'temperature') ) then iTracer = iTracer + 1 diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index 8c94d601364a..0e7ba4d51ddf 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -216,7 +216,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("subglacial hydro") -! *** TODO: Should basal melt rate calculation and column physics go inside RK loop? *** +! *** TODO: Should basal melt rate calculation, column physics, and hydrology go inside RK loop? *** call mpas_pool_get_array(geometryPool, 'rkTend', rkTend) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) @@ -240,6 +240,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then + config_rk_order = 1 rkWeights(:) = 1.0_RKIND rkSubstepWeights(:) = 1.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & @@ -297,10 +298,10 @@ subroutine li_time_integrator_forwardeuler(domain, err) call mpas_timer_stop("finalize damage") endif - ! Update thickness, damage, temperature or enthalpy, passiveTracer2d - ! before final velocity solve - if (rkStep .eq. config_rk_order) then ! or whichever is last rk step - ! y_{n+1} = y_n + (k1 + k2) / 2 + ! If using Runge-Kutta time integration, update thickness, damage, + ! temperature or enthalpy, passiveTracer2d before final velocity solve. + ! If using forward Euler, these have already been updated above. + if ( (rkStep .eq. config_rk_order) .and. (config_rk_order > 1) ) then ! Need layerThickness first in order to update other tracers if (trim(config_thickness_advection) == 'fct') then From af2170588a51c0b7dd10af04a785d3c2ede4466c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 4 Sep 2023 11:28:56 -0700 Subject: [PATCH 09/36] Rename files, modules, and routines to include Runge-Kutta Rename files, modules, and routines to include Runge-Kutta --- .../mpas-albany-landice/src/landice.cmake | 2 +- .../src/mode_forward/Makefile | 6 ++-- .../mode_forward/mpas_li_time_integration.F | 6 ++-- ..._fe.F => mpas_li_time_integration_fe_rk.F} | 28 +++++++++---------- 4 files changed, 21 insertions(+), 21 deletions(-) rename components/mpas-albany-landice/src/mode_forward/{mpas_li_time_integration_fe.F => mpas_li_time_integration_fe_rk.F} (98%) diff --git a/components/mpas-albany-landice/src/landice.cmake b/components/mpas-albany-landice/src/landice.cmake index 3b550fd34067..e4f0c7d4c7d1 100644 --- a/components/mpas-albany-landice/src/landice.cmake +++ b/components/mpas-albany-landice/src/landice.cmake @@ -57,7 +57,7 @@ list(APPEND RAW_SOURCES core_landice/mode_forward/mpas_li_core.F core_landice/mode_forward/mpas_li_core_interface.F core_landice/mode_forward/mpas_li_time_integration.F - core_landice/mode_forward/mpas_li_time_integration_fe.F + core_landice/mode_forward/mpas_li_time_integration_fe_rk.F core_landice/mode_forward/mpas_li_diagnostic_vars.F core_landice/mode_forward/mpas_li_advection.F core_landice/mode_forward/mpas_li_calving.F diff --git a/components/mpas-albany-landice/src/mode_forward/Makefile b/components/mpas-albany-landice/src/mode_forward/Makefile index c2c9d880ca6b..380be8cb0202 100644 --- a/components/mpas-albany-landice/src/mode_forward/Makefile +++ b/components/mpas-albany-landice/src/mode_forward/Makefile @@ -4,7 +4,7 @@ OBJS = mpas_li_core.o \ mpas_li_core_interface.o \ mpas_li_time_integration.o \ - mpas_li_time_integration_fe.o \ + mpas_li_time_integration_fe_rk.o \ mpas_li_diagnostic_vars.o \ mpas_li_advection.o \ mpas_li_advection_fct_shared.o \ @@ -35,9 +35,9 @@ mpas_li_core.o: mpas_li_time_integration.o \ mpas_li_statistics.o \ mpas_li_calving.o -mpas_li_time_integration.o: mpas_li_time_integration_fe.o +mpas_li_time_integration.o: mpas_li_time_integration_fe_rk.o -mpas_li_time_integration_fe.o: mpas_li_advection.o \ +mpas_li_time_integration_fe_rk.o: mpas_li_advection.o \ mpas_li_calving.o \ mpas_li_thermal.o \ mpas_li_iceshelf_melt.o \ diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F index f85393029bce..d8bd4ae94a67 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F @@ -26,7 +26,7 @@ module li_time_integration use mpas_timekeeping use mpas_log - use li_time_integration_fe + use li_time_integration_fe_rk use li_setup use li_constants @@ -166,9 +166,9 @@ subroutine li_timestep(domain, err) !call mpas_log_write('Using ' // trim(config_time_integration) // ' time integration.') select case (config_time_integration) case ('forward_euler') - call li_time_integrator_forwardeuler(domain, err_tmp) + call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp) case ('runge_kutta') - call li_time_integrator_forwardeuler(domain, err_tmp) + call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp) case default call mpas_log_write(trim(config_time_integration) // ' is not a valid land ice time integration option.', MPAS_LOG_ERR) err_tmp = 1 diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F similarity index 98% rename from components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F rename to components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F index 0e7ba4d51ddf..c6f1ae2c352f 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F @@ -9,17 +9,17 @@ !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! li_time_integration_fe +! li_time_integration_fe_rk ! -!> \brief MPAS land ice Forward Euler time integration scheme -!> \author Matt Hoffman -!> \date 17 April 2011 +!> \brief MPAS land ice Forward Euler and Runge-Kutta time integration schemes +!> \author Matt Hoffman, Trevor Hillebrand +!> \date 17 April 2011, Runge-Kutta added Sept 2023 !> \details -!> This module contains the Forward Euler time integration scheme +!> This module contains the Forward Euler and Runge-Kutta time integration schemes ! !----------------------------------------------------------------------- -module li_time_integration_fe +module li_time_integration_fe_rk use mpas_derived_types use mpas_pool_routines @@ -52,7 +52,7 @@ module li_time_integration_fe ! !-------------------------------------------------------------------- - public :: li_time_integrator_forwardeuler + public :: li_time_integrator_forwardeuler_rungekutta !-------------------------------------------------------------------- ! @@ -68,16 +68,16 @@ module li_time_integration_fe !*********************************************************************** ! -! routine li_time_integrator_forwardeuler +! routine li_time_integrator_forwardeuler_rungekutta ! -!> \brief Forward Euler time integration scheme +!> \brief Forward Euler and Runge-Kutta time integration schemes !> \author Matthew Hoffman !> \date 10 January 2012 !> \details -!> This routine performs Forward Euler time integration. +!> This routine performs Forward Euler and Runge-Kutta time integration. ! !----------------------------------------------------------------------- - subroutine li_time_integrator_forwardeuler(domain, err) + subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) use li_subglacial_hydro use li_velocity @@ -424,7 +424,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) ! === error check if (err == 1) then - call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler.", MPAS_LOG_ERR) + call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler_rungekutta.", MPAS_LOG_ERR) endif deallocate(temperatureProv) @@ -435,7 +435,7 @@ subroutine li_time_integrator_forwardeuler(domain, err) deallocate(passiveTracer3dProv) deallocate(passiveTracer3d) !-------------------------------------------------------------------- - end subroutine li_time_integrator_forwardeuler + end subroutine li_time_integrator_forwardeuler_rungekutta @@ -1414,5 +1414,5 @@ subroutine advance_clock(domain, err) end subroutine advance_clock -end module li_time_integration_fe +end module li_time_integration_fe_rk From 63a169be3de0a5dc7d5f1fe3842a6d7d42a79ad8 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 4 Sep 2023 16:20:40 -0700 Subject: [PATCH 10/36] Add SSP-RK3 time integration Add support for config_rk_order = 3, using the Strong Stability- Preserving RK3 method. This follows a slightly different algorithm from the classic RK2 and RK4 methods, and so requires significantl increasing the complexity of the RK loop. See eq 2.48 (pg 56) in Durran (2010): Numerical Methods for Fluid Dynamics for the description of this method. --- .../mpas_li_time_integration_fe_rk.F | 125 ++++++++++++++++-- 1 file changed, 111 insertions(+), 14 deletions(-) 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 c6f1ae2c352f..4c99a138c6fd 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 @@ -113,7 +113,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) character (len=StrKIND), pointer :: config_thermal_solver character (len=StrKIND), pointer :: config_time_integration integer, pointer :: config_rk_order - integer :: rkStep, iCell, iTracer, k + integer :: rkStage, iCell, iTracer, k real (kind=RKIND), dimension(:,:,:), pointer :: rkTend real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & temperature, & @@ -122,6 +122,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessProv, & + layerThicknessTmp, & temperatureProv, & enthalpyProv, & passiveTracer3d, & @@ -133,6 +134,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND), pointer :: deltat real (kind=RKIND) :: deltatFull real (kind=RKIND), dimension(4) :: rkWeights, rkSubstepWeights + real (kind=RKIND), dimension(3) :: rkSSPweights err = 0 err_tmp = 0 @@ -158,6 +160,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(temperatureProv(nVertLevels, nCells+1)) allocate(enthalpyProv(nVertLevels, nCells+1)) allocate(layerThicknessProv(nVertLevels, nCells+1)) + allocate(layerThicknessTmp(nVertLevels, nCells+1)) allocate(damage3dProv(nVertLevels, nCells+1)) allocate(damage3d(nVertLevels, nCells+1)) allocate(passiveTracer3dProv(nVertLevels, nCells+1)) @@ -166,6 +169,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) temperatureProv(:,:) = 0.0_RKIND enthalpyProv(:,:) = 0.0_RKIND layerThicknessProv(:,:) = 0.0_RKIND + layerThicknessTmp(:,:) = 0.0_RKIND damage3dProv(:,:) = 0.0_RKIND damage3d(:,:) = 0.0_RKIND passiveTracer3dProv(:,:) = 0.0_RKIND @@ -237,6 +241,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkTend(:,:,:) = 0.0_RKIND deltatFull = deltat ! Save deltat in order to reset it at end of RK loop + ! Set RK weights based on desired time integration method. Note + ! that rkSubstepWeights are used to update at each sub-step, and + ! are thus offset from the typical writing of the coefficients + ! by one index. e.g., the coefficients for RK4 are usually written + ! (0, 1/2, 1/2, 1), while we use (1/2, 1/2, 1). The last entry of 1.0 + ! is simply for ease of implementation. + rkSSPweights(:) = 1.0_RKIND if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then @@ -249,9 +260,15 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! add config option to use midpoint rkWeights(:) = 0.5_RKIND rkSubstepWeights(:) = 1.0_RKIND - !elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & - ! (config_rk_order == 3) ) then - ! use Strong Stability Preserving RK3? + elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & + (config_rk_order == 3) ) then + ! use Strong Stability Preserving RK3 + rkWeights(:) = 1.0_RKIND + rkSubstepWeights(:) = 1.0_RKIND + + rkSSPweights(1) = 1.0_RKIND + rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND + rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 4) ) then rkWeights(1) = 1.0_RKIND / 6.0_RKIND @@ -262,7 +279,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSubstepWeights(1) = 0.5_RKIND rkSubstepWeights(2) = 0.5_RKIND rkSubstepWeights(3) = 1.0_RKIND - rkSubstepWeights(4) = 1.0_RKIND + rkSubstepWeights(4) = 1.0_RKIND ! not used else err = 1 call mpas_log_write('config_time_integration = ' // trim(config_time_integration) & @@ -270,21 +287,22 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) intArgs=(/config_rk_order/), messageType=MPAS_LOG_ERR) endif ! *** Start RK loop *** - do rkStep = 1, config_rk_order + do rkStage = 1, config_rk_order call mpas_log_write('beginning rk step $i; rkSubstepWeight = $r', & - intArgs=(/rkStep/), realArgs=(/rkSubstepWeights(rkStep)/)) - deltat = deltatFull * rkSubstepWeights(rkStep) + intArgs=(/rkStage/), realArgs=(/rkSubstepWeights(rkStage)/)) + deltat = deltatFull * rkSubstepWeights(rkStage) + ! === calculate damage =========== if (config_calculate_damage) then call mpas_timer_start("damage") - call li_calculate_damage(domain, rkTend(2,:,:), rkWeights(rkStep), err_tmp) + call li_calculate_damage(domain, rkTend(2,:,:), rkWeights(rkStage), err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("damage") endif ! === Compute new state for prognostic variables === call mpas_timer_start("advect thickness and tracers") - call advection_solver(domain, rkTend, rkWeights(rkStep), err_tmp) + call advection_solver(domain, rkTend, rkWeights(rkStage), err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("advect thickness and tracers") @@ -298,10 +316,88 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("finalize damage") endif - ! If using Runge-Kutta time integration, update thickness, damage, + ! If using SSP RK3, then update thickness and tracers incrementally. + ! For first RK stage, thickness and tracer updates above are sufficient + if ( (config_rk_order == 3) .and. (rkStage > 1 ) ) then + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) + call mpas_pool_get_array(geometryPool, 'damage', damage) + + call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) + call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + + layerThicknessTmp(:,:) = layerThickness(:,:) + layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessProv(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) + thickness = sum(layerThickness, 1) + + if (trim(config_thermal_solver) == 'enthalpy') then + where (layerThickness(:,:) > 0.0_RKIND) + enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyProv(:,:) * layerThicknessProv(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * layerThicknessTmp(:,:) ) / layerThickness(:,:) + ! elsewhere, what? + end where + ! given the enthalpy, compute the temperature and water fraction + do iCell = 1, nCells + + call li_enthalpy_to_temperature_kelvin(& + layerCenterSigma, & + thickness(iCell), & + enthalpy(:,iCell), & + temperature(:,iCell), & + waterFrac(:,iCell)) + + enddo + + elseif (trim(config_thermal_solver) == 'temperature') then + where (layerThickness(:,:) > 0.0_RKIND) + temperature(:,:) = ( rkSSPweights(rkStage) * temperatureProv(:,:) * layerThicknessProv(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * layerThicknessTmp(:,:) ) / layerThickness(:,:) + ! elsewhere, what? + end where + endif + + if (config_calculate_damage) then + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + else + damage3d(k,iCell) = 0.0_RKIND + endif + enddo + damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) + enddo + endif + + do iCell = 1, nCells + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + else + passiveTracer3d(k,iCell) = 0.0_RKIND + endif + enddo + passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) + enddo + + ! === Ensure damage is within bounds after full time integration === + if ( (rkStage .eq. config_rk_order) .and. (config_finalize_damage_after_advection) ) then + call mpas_timer_start("finalize damage") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") + endif + endif ! if config_rk_order == 3 + + ! If using RK2 or RK4 time integration, update thickness, damage, ! temperature or enthalpy, passiveTracer2d before final velocity solve. - ! If using forward Euler, these have already been updated above. - if ( (rkStep .eq. config_rk_order) .and. (config_rk_order > 1) ) then + ! If using forward Euler or SSPRK3, these have already been updated above. + if ( (rkStage .eq. config_rk_order) .and. ( (config_rk_order == 2) .or. (config_rk_order == 4) ) ) then ! Need layerThickness first in order to update other tracers if (trim(config_thickness_advection) == 'fct') then @@ -313,7 +409,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) where (layerThickness(:,:) > 0.0_RKIND) enthalpy(:,:) = ( enthalpyProv(:,:) * layerThicknessProv(:,:) + & deltatFull * rkTend(1,:,:) ) / layerThickness(:,:) - ! elsewhere, what? + ! elsewhere, what? end where ! given the enthalpy, compute the temperature and water fraction do iCell = 1, nCells @@ -430,6 +526,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(temperatureProv) deallocate(enthalpyProv) deallocate(layerThicknessProv) + deallocate(layerThicknessTmp) deallocate(damage3dProv) deallocate(damage3d) deallocate(passiveTracer3dProv) From b920b334e394319aab83abfe0c60b7c99cdd286d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 6 Sep 2023 14:12:43 -0700 Subject: [PATCH 11/36] Update RK2 and RK4 substeps using 0th step solution Previous commits inadvertently updated sub-step solutions incrementally (i.e., using y(k-1) to solve for y(k), where k is the substep), when they should have been using y0. This commit fixes that, although tests show that this treatment is not complete and solutions are unstable. --- .../mpas-albany-landice/src/Registry.xml | 3 + .../src/mode_forward/mpas_li_advection.F | 15 +-- .../src/mode_forward/mpas_li_calving.F | 7 +- .../mpas_li_time_integration_fe_rk.F | 91 ++++++++++--------- 4 files changed, 64 insertions(+), 52 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9555a7a2b592..ff92c44e6b36 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -1188,6 +1188,9 @@ is the value of that variable from the *previous* time level! description="diagnostic field of dynamic thickening rate (calculated as negative of flux divergence)" /> + 0.0_RKIND) enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyProv(:,:) * layerThicknessProv(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * layerThicknessTmp(:,:) ) / layerThickness(:,:) + (1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * & + layerThicknessTmp(:,:) ) / layerThickness(:,:) ! elsewhere, what? end where ! given the enthalpy, compute the temperature and water fraction @@ -354,7 +356,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) elseif (trim(config_thermal_solver) == 'temperature') then where (layerThickness(:,:) > 0.0_RKIND) temperature(:,:) = ( rkSSPweights(rkStage) * temperatureProv(:,:) * layerThicknessProv(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * layerThicknessTmp(:,:) ) / layerThickness(:,:) + (1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * & + layerThicknessTmp(:,:) ) / layerThickness(:,:) ! elsewhere, what? end where endif @@ -364,7 +367,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) else damage3d(k,iCell) = 0.0_RKIND endif @@ -377,7 +381,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) else passiveTracer3d(k,iCell) = 0.0_RKIND endif @@ -385,30 +390,27 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo - ! === Ensure damage is within bounds after full time integration === - if ( (rkStage .eq. config_rk_order) .and. (config_finalize_damage_after_advection) ) then - call mpas_timer_start("finalize damage") - call li_finalize_damage_after_advection(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("finalize damage") - endif endif ! if config_rk_order == 3 ! If using RK2 or RK4 time integration, update thickness, damage, ! temperature or enthalpy, passiveTracer2d before final velocity solve. ! If using forward Euler or SSPRK3, these have already been updated above. - if ( (rkStage .eq. config_rk_order) .and. ( (config_rk_order == 2) .or. (config_rk_order == 4) ) ) then + if ( (config_rk_order == 2) .or. (config_rk_order == 4) ) then + if (rkStage == nRKstages) then + deltat = deltatFull + rkTend(:,:,:) = rkTendAccum(:,:,:) + endif ! Need layerThickness first in order to update other tracers if (trim(config_thickness_advection) == 'fct') then - layerThickness(:,:) = layerThicknessProv(:,:) + deltatFull * rkTend(4,:,:) + layerThickness(:,:) = layerThicknessProv(:,:) + deltat * rkTend(4,:,:) thickness = sum(layerThickness, 1) endif if (trim(config_thermal_solver) == 'enthalpy') then where (layerThickness(:,:) > 0.0_RKIND) enthalpy(:,:) = ( enthalpyProv(:,:) * layerThicknessProv(:,:) + & - deltatFull * rkTend(1,:,:) ) / layerThickness(:,:) + deltat * rkTend(1,:,:) ) / layerThickness(:,:) ! elsewhere, what? end where ! given the enthalpy, compute the temperature and water fraction @@ -426,17 +428,18 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) elseif (trim(config_thermal_solver) == 'temperature') then where (layerThickness(:,:) > 0.0_RKIND) temperature(:,:) = ( temperatureProv(:,:) * layerThicknessProv(:,:) + & - deltatFull * rkTend(1,:,:) ) / layerThickness(:,:) + deltat * rkTend(1,:,:) ) / layerThickness(:,:) ! elsewhere, what? end where endif + damage3d(:,:) = 0.0_RKIND if (config_calculate_damage) then do iCell = 1, nCells do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then damage3d(k,iCell) = ( damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltatFull * rkTend(2,k,iCell) ) / layerThickness(k,iCell) + deltat * rkTend(2,k,iCell) ) / layerThickness(k,iCell) else damage3d(k,iCell) = 0.0_RKIND endif @@ -445,11 +448,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo endif + passiveTracer3d(:,:) = 0.0_RKIND do iCell = 1, nCells do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then passiveTracer3d(k,iCell) = ( passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltatFull * rkTend(3,k,iCell) ) / layerThickness(k,iCell) + deltat * rkTend(3,k,iCell) ) / layerThickness(k,iCell) else passiveTracer3d(k,iCell) = 0.0_RKIND endif @@ -457,15 +461,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo - ! === Ensure damage is within bounds after full time integration === - if (config_finalize_damage_after_advection) then - call mpas_timer_start("finalize damage") - call li_finalize_damage_after_advection(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("finalize damage") - endif + endif ! config_rk_order == 2 or 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") + call li_finalize_damage_after_advection(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("finalize damage") endif - + ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) @@ -897,7 +902,7 @@ end subroutine prepare_advection !> in calculate tendencies are now in prepare_advection. !----------------------------------------------------------------------- - subroutine advection_solver(domain, rkTend, rkWeight, err) + subroutine advection_solver(domain, rkTendAccum, rkTend, rkWeight, err) use mpas_timekeeping use li_mask @@ -910,7 +915,7 @@ subroutine advection_solver(domain, rkTend, rkWeight, err) ! input/output variables !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - real (kind=RKIND), dimension(:,:,:), intent(inout) :: rkTend + real (kind=RKIND), dimension(:,:,:), intent(inout) :: rkTendAccum, rkTend !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- @@ -1028,6 +1033,7 @@ subroutine advection_solver(domain, rkTend, rkWeight, err) geometryPool, & thermalPool, & scratchPool, & + rkTendAccum, & rkTend, & rkWeight, & err_tmp, & @@ -1049,6 +1055,7 @@ subroutine advection_solver(domain, rkTend, rkWeight, err) geometryPool, & thermalPool, & scratchPool, & + rkTendAccum, & rkTend, & rkWeight, & err_tmp, & From dc33a678da295453c9a7aa160c0fa6c753a7d39b Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 7 Sep 2023 14:57:15 -0700 Subject: [PATCH 12/36] Only use SSP RK2 and RK3 methods Only use SSP RK2 and RK3 methods, and remove support for classical RK2 and RK4. --- .../mpas-albany-landice/src/Registry.xml | 10 +- .../src/mode_forward/mpas_li_advection.F | 42 +----- .../src/mode_forward/mpas_li_calving.F | 38 ++---- .../mpas_li_time_integration_fe_rk.F | 127 +++--------------- 4 files changed, 37 insertions(+), 180 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index ff92c44e6b36..688761cb8f94 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -489,8 +489,8 @@ possible_values="'forward_euler' or 'runge_kutta'" /> - - 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 95f5ead00ba1..7d93f114b797 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 @@ -81,9 +81,6 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & scratchPool, & - rkTendAccum, & - tend, & - rkWeight, & err, & advectTracersIn) @@ -99,8 +96,6 @@ subroutine li_advection_thickness_tracers(& real (kind=RKIND), intent(in) :: & dt !< Input: time step (s) - real (kind=RKIND), intent(in) :: rkWeight - type (mpas_pool_type), intent(in) :: & meshPool !< Input: mesh information @@ -115,9 +110,6 @@ subroutine li_advection_thickness_tracers(& !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - real (kind=RKIND), dimension(:,:,:), intent(inout) :: & - rkTendAccum !< Input/output: accumulated tendencies for RK time integration - type (mpas_pool_type), intent(inout) :: & velocityPool !< Input/output: velocity information ! (needs to be inout for li_calculate_mask call @@ -137,7 +129,6 @@ subroutine li_advection_thickness_tracers(& ! !----------------------------------------------------------------- - real (kind=RKIND), dimension(:,:,:), intent(out) :: tend integer, intent(out) :: & err !< Output: error flag @@ -216,14 +207,12 @@ subroutine li_advection_thickness_tracers(& real (kind=RKIND), dimension(:), pointer :: dvEdge character (len=StrKIND), pointer :: & - config_thickness_advection, & ! method for advecting thickness and tracers - config_thermal_solver + config_thickness_advection ! method for advecting thickness and tracers logical, pointer :: & config_print_thickness_advection_info, & !TODO - change to config_print_advection_info? config_print_thermal_info, & - config_thermal_calculate_bmb, & - config_calculate_damage + config_thermal_calculate_bmb real (kind=RKIND), pointer :: & config_ice_density, & ! rhoi @@ -236,7 +225,7 @@ subroutine li_advection_thickness_tracers(& logical :: advectTracers ! if true, then advect tracers as well as thickness - integer :: iCell1, iCell2, theGroundedCell, iTracer + integer :: iCell1, iCell2, theGroundedCell integer, parameter :: horizAdvOrder = 4 integer, parameter :: vertAdvOrder = 1 @@ -326,8 +315,6 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) call mpas_pool_get_config(liConfigs, 'config_num_halos', config_num_halos) - call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) - call mpas_pool_get_config(liConfigs, 'config_calculate_damage', config_calculate_damage) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) @@ -527,26 +514,11 @@ subroutine li_advection_thickness_tracers(& call mpas_log_write(trim(config_tracer_advection) // & ' tracer advection is not currently supported with fct thickness advection.', MPAS_LOG_ERR) endif - - ! Update tendencies for RK integration - iTracer = 0 - if ( (trim(config_thermal_solver) == 'enthalpy') .or. (trim(config_thermal_solver) == 'temperature') ) then - iTracer = iTracer + 1 - rkTendAccum(1,:,:) = rkTendAccum(1,:,:) + tend(iTracer,:,:) * rkWeight - endif - if (config_calculate_damage) then - iTracer = iTracer + 1 - rkTendAccum(2,:,:) = rkTendAccum(2,:,:) + tend(iTracer,:,:) * rkWeight - endif - iTracer = iTracer + 1 - rkTendAccum(3,:,:) = rkTendAccum(3,:,:) + tend(iTracer,:,:) * rkWeight ! passiveTracer2d - iTracer = iTracer + 1 - rkTendAccum(4,:,:) = rkTendAccum(4,:,:) + tend(iTracer,:,:) * rkWeight ! layerThickness else - err_tmp = 1 - call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & - ", config_tracer_advection = " // trim(config_tracer_advection) // & - " is not a supported combination.", MPAS_LOG_ERR) + err_tmp = 1 + call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & + ", config_tracer_advection = " // trim(config_tracer_advection) // & + " is not a supported combination.", MPAS_LOG_ERR) endif if (config_print_thickness_advection_info) then diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 51ebea915302..36302a1a78c3 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -3314,12 +3314,12 @@ end subroutine damage_calving !----------------------------------------------------------------------- - subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) + subroutine li_calculate_damage(domain, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- - real (kind=RKIND), intent(in) :: rkWeight + !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- @@ -3330,7 +3330,7 @@ subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - real (kind=RKIND), dimension(:,:), intent(out) :: rkTend, rkTendAccum + !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -3348,7 +3348,6 @@ subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) real(kind=RKIND), pointer :: config_flowLawExponent logical, pointer :: config_print_calving_info - real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:), pointer :: thickness real (kind=RKIND), dimension(:), pointer :: eMax real (kind=RKIND), dimension(:), pointer :: tauMax, tauMin @@ -3369,8 +3368,8 @@ subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) real (kind=RKIND), pointer :: deltat !< time step (s) integer, dimension(:), pointer :: nEdgesOnCell ! number of cells that border each cell integer, dimension(:,:), pointer :: cellsOnCell ! list of cells that neighbor each cell - integer, pointer :: nCells, nVertLevels - integer :: iCell, jCell, iNeighbor, n_damage_downstream, k + integer, pointer :: nCells + integer :: iCell, jCell, iNeighbor, n_damage_downstream real(kind=RKIND) :: damage_downstream real(kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY real(kind=RKIND), dimension(6) :: localMinInfo, localMaxInfo, globalMinInfo, globalMaxInfo @@ -3397,13 +3396,11 @@ subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) ! get fields call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_array(meshPool, 'deltat', deltat) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'damage', damage) call mpas_pool_get_array(geometryPool, 'ddamagedt', ddamagedt) call mpas_pool_get_array(geometryPool, 's0', s0) @@ -3458,11 +3455,6 @@ subroutine li_calculate_damage(domain, rkTendAccum, rkTend, rkWeight, err) ! damageSource(:) * (1+exp(-uReconstructX(1,:)*seconds/1000.0))) ddamagedt(:) = damageSource(:) * damage(:) - ! Update damage tendency for RK time integration - do k = 1, nVertLevels - rkTend(k,:) = ddamagedt(:) * layerThickness(k,:) - rkTendAccum(k,:) = rkTendAccum(k,:) + rkTend(k,:) * rkWeight - enddo damage(:) = damage(:) + ddamagedt(:) * deltat @@ -3684,18 +3676,14 @@ subroutine li_finalize_damage_after_advection(domain, err) endif ! put the damageMax value back to preserve the damage value (no heal) - ! limit damage to [damageNye, 1.0] - do iCell = 1, nCells - if (damage(iCell) < 0.0_RKIND) then - damage(iCell) = 0.0_RKIND - end if - if (damage(iCell) < damageNye(iCell)) then - damage(iCell) = damageNye(iCell) - end if - if (damage(iCell) > 1.0_RKIND) then - damage(iCell) = 1.0_RKIND - end if - end do + where (damage < 0.0_RKIND) + damage = 0.0_RKIND + end where + + where (damage > 1.0_RKIND) + damage = 1.0_RKIND + end where + do iCell = 1, nCells if (li_mask_is_grounded_ice(cellMask(iCell)) .or. .not. li_mask_is_ice(cellMask(iCell))) then 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 acc8a897c07c..5375348298eb 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 @@ -71,8 +71,8 @@ module li_time_integration_fe_rk ! routine li_time_integrator_forwardeuler_rungekutta ! !> \brief Forward Euler and Runge-Kutta time integration schemes -!> \author Matthew Hoffman -!> \date 10 January 2012 +!> \author Matthew Hoffman, Trevor Hillebrand +!> \date 10 January 2012, updated Sept 2023 for Runge-Kutta !> \details !> This routine performs Forward Euler and Runge-Kutta time integration. ! @@ -114,7 +114,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) character (len=StrKIND), pointer :: config_time_integration integer, pointer :: config_rk_order integer :: rkStage, iCell, iTracer, k - real (kind=RKIND), dimension(:,:,:), pointer :: rkTendAccum, rkTend real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & temperature, & enthalpy, & @@ -133,7 +132,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) integer, pointer :: nCells real (kind=RKIND), pointer :: deltat real (kind=RKIND) :: deltatFull - real (kind=RKIND), dimension(4) :: rkWeights, rkSubstepWeights + real (kind=RKIND), dimension(4) :: rkSubstepWeights real (kind=RKIND), dimension(3) :: rkSSPweights integer :: nRKstages @@ -222,8 +221,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("subglacial hydro") ! *** TODO: Should basal melt rate calculation, column physics, and hydrology go inside RK loop? *** - call mpas_pool_get_array(geometryPool, 'rkTendAccum', rkTendAccum) - call mpas_pool_get_array(geometryPool, 'rkTend', rkTend) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) @@ -240,8 +237,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3dProv(k,:) = damage(:) passiveTracer3dProv(k,:) = passiveTracer2d(:) enddo - rkTendAccum(:,:,:) = 0.0_RKIND - rkTend(:,:,:) = 0.0_RKIND deltatFull = deltat ! Save deltat in order to reset it at end of RK loop ! Set RK weights based on desired time integration method. Note @@ -256,38 +251,27 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) (config_rk_order == 1)) ) then config_rk_order = 1 nRKstages = 1 - rkWeights(:) = 1.0_RKIND rkSubstepWeights(:) = 1.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 2) ) then ! use classical (i.e., endpoint, or Heun's method) RK2. Could also - ! add config option to use midpoint - nRKstages = 1 - rkWeights(:) = 0.5_RKIND + ! add config option to use standard endpoint or midpoint methods, but + ! these are algorithmically more complex and less suitable for our domains + nRKstages = 2 rkSubstepWeights(:) = 1.0_RKIND + + rkSSPweights(1) = 1.0_RKIND + rkSSPweights(2) = 0.5_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 3) ) then ! use Strong Stability Preserving RK3 nRKstages = 3 - rkWeights(:) = 1.0_RKIND rkSubstepWeights(:) = 1.0_RKIND rkSSPweights(1) = 1.0_RKIND rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND ! Add option for 4-stage SSP-RK3? Allows for CFL=2. - elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & - (config_rk_order == 4) ) then - nRKstages = 4 - rkWeights(1) = 1.0_RKIND / 6.0_RKIND - rkWeights(2) = 1.0_RKIND / 3.0_RKIND - rkWeights(3) = 1.0_RKIND / 3.0_RKIND - rkWeights(4) = 1.0_RKIND / 6.0_RKIND - - rkSubstepWeights(1) = 0.5_RKIND - rkSubstepWeights(2) = 0.5_RKIND - rkSubstepWeights(3) = 1.0_RKIND - rkSubstepWeights(4) = 1.0_RKIND ! not used else err = 1 call mpas_log_write('config_time_integration = ' // trim(config_time_integration) & @@ -299,27 +283,24 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_log_write('beginning rk step $i; rkSubstepWeight = $r', & intArgs=(/rkStage/), realArgs=(/rkSubstepWeights(rkStage)/)) deltat = deltatFull * rkSubstepWeights(rkStage) - rkTend(:,:,:) = 0.0_RKIND ! === calculate damage =========== if (config_calculate_damage) then call mpas_timer_start("damage") - call li_calculate_damage(domain, rkTendAccum(2,:,:), rkTend(2,:,:), rkWeights(rkStage), err_tmp) + call li_calculate_damage(domain, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("damage") endif ! === Compute new state for prognostic variables === call mpas_timer_start("advect thickness and tracers") - call advection_solver(domain, rkTendAccum, rkTend, rkWeights(rkStage), err_tmp) + call advection_solver(domain, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("advect thickness and tracers") - call mpas_dmpar_field_halo_exch(domain, 'rkTendAccum') - - ! If using SSP RK3, then update thickness and tracers incrementally. + ! If using SSP RK, then update thickness and tracers incrementally. ! For first RK stage, thickness and tracer updates above are sufficient - if ( (config_rk_order == 3) .and. (rkStage > 1 ) ) then + if ( rkStage > 1 ) then call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) @@ -391,77 +372,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo endif ! if config_rk_order == 3 - - ! If using RK2 or RK4 time integration, update thickness, damage, - ! temperature or enthalpy, passiveTracer2d before final velocity solve. - ! If using forward Euler or SSPRK3, these have already been updated above. - if ( (config_rk_order == 2) .or. (config_rk_order == 4) ) then - if (rkStage == nRKstages) then - deltat = deltatFull - rkTend(:,:,:) = rkTendAccum(:,:,:) - endif - - ! Need layerThickness first in order to update other tracers - if (trim(config_thickness_advection) == 'fct') then - layerThickness(:,:) = layerThicknessProv(:,:) + deltat * rkTend(4,:,:) - thickness = sum(layerThickness, 1) - endif - - if (trim(config_thermal_solver) == 'enthalpy') then - where (layerThickness(:,:) > 0.0_RKIND) - enthalpy(:,:) = ( enthalpyProv(:,:) * layerThicknessProv(:,:) + & - deltat * rkTend(1,:,:) ) / layerThickness(:,:) - ! elsewhere, what? - end where - ! given the enthalpy, compute the temperature and water fraction - do iCell = 1, nCells - - call li_enthalpy_to_temperature_kelvin(& - layerCenterSigma, & - thickness(iCell), & - enthalpy(:,iCell), & - temperature(:,iCell), & - waterFrac(:,iCell)) - - enddo - - elseif (trim(config_thermal_solver) == 'temperature') then - where (layerThickness(:,:) > 0.0_RKIND) - temperature(:,:) = ( temperatureProv(:,:) * layerThicknessProv(:,:) + & - deltat * rkTend(1,:,:) ) / layerThickness(:,:) - ! elsewhere, what? - end where - endif - - damage3d(:,:) = 0.0_RKIND - if (config_calculate_damage) then - do iCell = 1, nCells - do k = 1, nVertLevels - if (layerThickness(k,iCell) > 0.0_RKIND) then - damage3d(k,iCell) = ( damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltat * rkTend(2,k,iCell) ) / layerThickness(k,iCell) - else - damage3d(k,iCell) = 0.0_RKIND - endif - enddo - damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) - enddo - endif - - passiveTracer3d(:,:) = 0.0_RKIND - do iCell = 1, nCells - do k = 1, nVertLevels - if (layerThickness(k,iCell) > 0.0_RKIND) then - passiveTracer3d(k,iCell) = ( passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & - deltat * rkTend(3,k,iCell) ) / layerThickness(k,iCell) - else - passiveTracer3d(k,iCell) = 0.0_RKIND - endif - enddo - passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) - enddo - - endif ! config_rk_order == 2 or 4 ! === Ensure damage is within bounds after full time integration === if ( (rkStage .eq. nRKstages) .and. (config_finalize_damage_after_advection) ) then @@ -902,7 +812,7 @@ end subroutine prepare_advection !> in calculate tendencies are now in prepare_advection. !----------------------------------------------------------------------- - subroutine advection_solver(domain, rkTendAccum, rkTend, rkWeight, err) + subroutine advection_solver(domain, err) use mpas_timekeeping use li_mask @@ -910,12 +820,11 @@ subroutine advection_solver(domain, rkTendAccum, rkTend, rkWeight, err) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- - real (kind=RKIND), intent(in) :: rkWeight + !----------------------------------------------------------------- ! input/output variables !----------------------------------------------------------------- type (domain_type), intent(inout) :: domain !< Input/Output: domain object - real (kind=RKIND), dimension(:,:,:), intent(inout) :: rkTendAccum, rkTend !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- @@ -1033,9 +942,6 @@ subroutine advection_solver(domain, rkTendAccum, rkTend, rkWeight, err) geometryPool, & thermalPool, & scratchPool, & - rkTendAccum, & - rkTend, & - rkWeight, & err_tmp, & advectTracersIn = .true.) @@ -1055,9 +961,6 @@ subroutine advection_solver(domain, rkTendAccum, rkTend, rkWeight, err) geometryPool, & thermalPool, & scratchPool, & - rkTendAccum, & - rkTend, & - rkWeight, & err_tmp, & advectTracersIn = .false.) From 74c014f9daf19263d7d51cd0ebae796226ea6b0c Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 7 Sep 2023 15:42:32 -0700 Subject: [PATCH 13/36] Some cleanup Rename some local variables, update comments, update log write --- .../mpas_li_time_integration_fe_rk.F | 66 +++++++++---------- 1 file changed, 33 insertions(+), 33 deletions(-) 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 5375348298eb..0854534a110b 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 @@ -120,14 +120,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) waterFrac real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions - real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessProv, & + real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & - temperatureProv, & - enthalpyProv, & + temperaturePrev, & + enthalpyPrev, & passiveTracer3d, & - passiveTracer3dProv, & + passiveTracer3dPrev, & damage3d, & - damage3dProv + damage3dPrev integer, pointer :: nVertLevels integer, pointer :: nCells real (kind=RKIND), pointer :: deltat @@ -157,22 +157,22 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) - allocate(temperatureProv(nVertLevels, nCells+1)) - allocate(enthalpyProv(nVertLevels, nCells+1)) - allocate(layerThicknessProv(nVertLevels, nCells+1)) + allocate(temperaturePrev(nVertLevels, nCells+1)) + allocate(enthalpyPrev(nVertLevels, nCells+1)) + allocate(layerThicknessPrev(nVertLevels, nCells+1)) allocate(layerThicknessTmp(nVertLevels, nCells+1)) - allocate(damage3dProv(nVertLevels, nCells+1)) + allocate(damage3dPrev(nVertLevels, nCells+1)) allocate(damage3d(nVertLevels, nCells+1)) - allocate(passiveTracer3dProv(nVertLevels, nCells+1)) + allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) - temperatureProv(:,:) = 0.0_RKIND - enthalpyProv(:,:) = 0.0_RKIND - layerThicknessProv(:,:) = 0.0_RKIND + temperaturePrev(:,:) = 0.0_RKIND + enthalpyPrev(:,:) = 0.0_RKIND + layerThicknessPrev(:,:) = 0.0_RKIND layerThicknessTmp(:,:) = 0.0_RKIND - damage3dProv(:,:) = 0.0_RKIND + damage3dPrev(:,:) = 0.0_RKIND damage3d(:,:) = 0.0_RKIND - passiveTracer3dProv(:,:) = 0.0_RKIND + passiveTracer3dPrev(:,:) = 0.0_RKIND passiveTracer3d(:,:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== @@ -230,12 +230,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) ! Save relevant fields before RK loop, to be used in update at the end - temperatureProv(:,:) = temperature(:,:) - enthalpyProv(:,:) = enthalpy(:,:) - layerThicknessProv(:,:) = layerThickness(:,:) + temperaturePrev(:,:) = temperature(:,:) + enthalpyPrev(:,:) = enthalpy(:,:) + layerThicknessPrev(:,:) = layerThickness(:,:) do k = 1, nVertLevels - damage3dProv(k,:) = damage(:) - passiveTracer3dProv(k,:) = passiveTracer2d(:) + damage3dPrev(k,:) = damage(:) + passiveTracer3dPrev(k,:) = passiveTracer2d(:) enddo deltatFull = deltat ! Save deltat in order to reset it at end of RK loop @@ -254,7 +254,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSubstepWeights(:) = 1.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 2) ) then - ! use classical (i.e., endpoint, or Heun's method) RK2. Could also + ! use Strong Stability Preserving RK2. Could also ! add config option to use standard endpoint or midpoint methods, but ! these are algorithmically more complex and less suitable for our domains nRKstages = 2 @@ -280,8 +280,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; rkSubstepWeight = $r', & - intArgs=(/rkStage/), realArgs=(/rkSubstepWeights(rkStage)/)) + call mpas_log_write('beginning rk step $i', & + intArgs=(/rkStage/)) deltat = deltatFull * rkSubstepWeights(rkStage) ! === calculate damage =========== @@ -311,13 +311,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) layerThicknessTmp(:,:) = layerThickness(:,:) - layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessProv(:,:) + & + layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) thickness = sum(layerThickness, 1) if (trim(config_thermal_solver) == 'enthalpy') then where (layerThickness(:,:) > 0.0_RKIND) - enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyProv(:,:) * layerThicknessProv(:,:) + & + enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyPrev(:,:) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * & layerThicknessTmp(:,:) ) / layerThickness(:,:) ! elsewhere, what? @@ -336,7 +336,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) elseif (trim(config_thermal_solver) == 'temperature') then where (layerThickness(:,:) > 0.0_RKIND) - temperature(:,:) = ( rkSSPweights(rkStage) * temperatureProv(:,:) * layerThicknessProv(:,:) + & + temperature(:,:) = ( rkSSPweights(rkStage) * temperaturePrev(:,:) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * & layerThicknessTmp(:,:) ) / layerThickness(:,:) ! elsewhere, what? @@ -347,7 +347,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) do iCell = 1, nCells do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then - damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) else @@ -361,7 +361,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) do iCell = 1, nCells do k = 1, nVertLevels if (layerThickness(k,iCell) > 0.0_RKIND) then - passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dProv(k,iCell) * layerThicknessProv(k,iCell) + & + passiveTracer3d(k,iCell) = ( rkSSPweights(rkStage) * passiveTracer3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & (1.0_RKIND - rkSSPweights(rkStage)) * passiveTracer2d(iCell) * & layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) else @@ -438,13 +438,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_log_write("An error has occurred in li_time_integrator_forwardeuler_rungekutta.", MPAS_LOG_ERR) endif - deallocate(temperatureProv) - deallocate(enthalpyProv) - deallocate(layerThicknessProv) + deallocate(temperaturePrev) + deallocate(enthalpyPrev) + deallocate(layerThicknessPrev) deallocate(layerThicknessTmp) - deallocate(damage3dProv) + deallocate(damage3dPrev) deallocate(damage3d) - deallocate(passiveTracer3dProv) + deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From e938e9251797bfe7c12dc46ec3a3b1906cb93a33 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 19 Sep 2023 22:05:35 -0600 Subject: [PATCH 14/36] Add option for 4-stage SSPRK3 Add option for 4-stage SSPRK3, which involves four velocities solves instead of three, but also allows for a CFL fraction of 2. To use the 4-stage formulation, set 'config_rk_order = 3' and 'config_rk3_stages = 4'. --- .../mpas-albany-landice/src/Registry.xml | 4 ++ .../mpas_li_time_integration_fe_rk.F | 58 ++++++++++++++----- 2 files changed, 47 insertions(+), 15 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 688761cb8f94..f2559440c5ef 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -492,6 +492,10 @@ description="Order of Runge-Kutta time integration to use. A value of 1 is equivalent to forward euler. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods." possible_values="1, 2, 3" /> + block % next end do + ! If using 4-stage SSPRK3, CFL number of 2 is theoretically allowed + if ( (trim(config_time_integration) == 'runge_kutta') .and. & + (config_rk_order == 3) .and. (config_rk3_stages == 4) ) then + allowableAdvecDtOnProc = allowableAdvecDtOnProc * 2.0_RKIND + allowableDiffDtOnProc = allowableDiffDtOnProc * 2.0_RKIND + endif ! Local advective CFL info call mpas_set_timeInterval(allowableAdvecDtOnProcInterval, dt=allowableAdvecDtOnProc, ierr=err_tmp) From f0c9a1fd718fde40f09583adb9d15da514635e28 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 21 Sep 2023 17:39:55 -0600 Subject: [PATCH 15/36] Include SMB and BMB within RK loop Include SMB and BMB within RK loop by calculating a weighted average based on each specific RK formulation. --- .../mpas_li_time_integration_fe_rk.F | 71 +++++++++++++++++-- 1 file changed, 66 insertions(+), 5 deletions(-) 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 2191e2b6ebab..02ffc94e23b6 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 @@ -119,6 +119,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enthalpy, & waterFrac real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d + real (kind=RKIND), dimension(:), pointer :: sfcMassBalApplied, & + basalMassBalApplied, & + groundedSfcMassBalApplied, & + groundedBasalMassBalApplied, & + floatingBasalMassBalApplied + real (kind=RKIND), dimension(:), allocatable :: sfcMassBalAccum, & + basalMassBalAccum, & + groundedSfcMassBalAccum, & + groundedBasalMassBalAccum, & + floatingBasalMassBalAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & @@ -134,6 +144,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND) :: deltatFull real (kind=RKIND), dimension(4) :: rkSubstepWeights real (kind=RKIND), dimension(4) :: rkSSPweights + real (kind=RKIND), dimension(4) :: rkTendWeights ! Weights used for calculating budget terms integer :: nRKstages err = 0 @@ -167,6 +178,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) + allocate(sfcMassBalAccum(nCells+1)) + allocate(basalMassBalAccum(nCells+1)) + allocate(groundedSfcMassBalAccum(nCells+1)) + allocate(groundedBasalMassBalAccum(nCells+1)) + allocate(floatingBasalMassBalAccum(nCells+1)) + temperaturePrev(:,:) = 0.0_RKIND enthalpyPrev(:,:) = 0.0_RKIND layerThicknessPrev(:,:) = 0.0_RKIND @@ -176,6 +193,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) passiveTracer3dPrev(:,:) = 0.0_RKIND passiveTracer3d(:,:) = 0.0_RKIND + sfcMassBalAccum(:) = 0.0_RKIND + basalMassBalAccum(:) = 0.0_RKIND + groundedSfcMassBalAccum(:) = 0.0_RKIND + groundedBasalMassBalAccum(:) = 0.0_RKIND + floatingBasalMassBalAccum(:) = 0.0_RKIND + ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! call mpas_timer_start("advection prep") @@ -246,7 +269,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! by one index. e.g., the coefficients for RK4 are usually written ! (0, 1/2, 1/2, 1), while we use (1/2, 1/2, 1). The last entry of 1.0 ! is simply for ease of implementation. - rkSSPweights(:) = 1.0_RKIND + rkSSPweights(:) = 1.0_RKIND ! updated for each case below + rkTendWeights(:) = 0.0_RKIND ! updated for each case below if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then @@ -265,6 +289,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSSPweights(2) = 0.5_RKIND rkSSPweights(3) = 0.0_RKIND rkSSPweights(4) = 0.0_RKIND + + rkTendWeights(1) = 0.5_RKIND + rkTendWeights(2) = 0.5_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 3) ) then if (config_rk3_stages == 3) then @@ -276,6 +303,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND rkSSPweights(4) = 0.0_RKIND + + rkTendWeights(1) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(2) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(3) = 2.0_RKIND / 3.0_RKIND ! 4-stage SSP-RK3? Allows for CFL=2. elseif (config_rk3_stages == 4) then nRKstages = 4 @@ -285,6 +316,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) rkSSPweights(2) = 0.0_RKIND rkSSPweights(3) = 2.0_RKIND / 3.0_RKIND rkSSPweights(4) = 0.0_RKIND + + rkTendWeights(1) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(2) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(3) = 1.0_RKIND / 6.0_RKIND + rkTendWeights(4) = 1.0_RKIND / 2.0_RKIND else err = 1 call mpas_log_write('config_rk3_stages must 3 or 4', & @@ -390,16 +426,28 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo - endif ! if config_rk_order == 3 + endif - ! === Ensure damage is within bounds after full time integration === - if ( (rkStage .eq. nRKstages) .and. (config_finalize_damage_after_advection) ) then + ! === Ensure damage is within bounds before velocity solve === + if ( config_finalize_damage_after_advection ) then call mpas_timer_start("finalize damage") call li_finalize_damage_after_advection(domain, err_tmp) err = ior(err, err_tmp) call mpas_timer_stop("finalize damage") endif - + + call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) + call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) + ! update budgets + sfcMassBalAccum = sfcMassBalAccum + rkTendWeights(rkStage) * sfcMassBalApplied + groundedSfcMassBalAccum = groundedSfcMassBalAccum + rkTendWeights(rkStage) * groundedSfcMassBalApplied + basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied + groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied + floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied + ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) @@ -408,8 +456,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! *** end RK loop *** enddo +! Finalize budget updates + sfcMassBalApplied(:) = sfcMassBalAccum(:) + groundedSfcMassBalApplied(:) = sfcMassBalAccum(:) + basalMassBalApplied(:) = basalMassBalAccum(:) + groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) + floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) + ! Reset time step to full length after RK loop deltat = deltatFull + ! === Calve ice ======================== call mpas_timer_start("calve_ice") @@ -465,6 +521,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(damage3d) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) + deallocate(sfcMassBalAccum) + deallocate(basalMassBalAccum) + deallocate(groundedSfcMassBalAccum) + deallocate(groundedBasalMassBalAccum) + deallocate(floatingBasalMassBalAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From 79b4635376aa93cb642a8e4fbecb26b75f15d632 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 21 Sep 2023 18:02:50 -0600 Subject: [PATCH 16/36] Fix SMB and BMB budget for forward euler Fix SMB and BMB budget for forward euler, as previous commit caused these fields to be written out as all zeroes. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 02ffc94e23b6..5d1e14952ac4 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 @@ -270,7 +270,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! (0, 1/2, 1/2, 1), while we use (1/2, 1/2, 1). The last entry of 1.0 ! is simply for ease of implementation. rkSSPweights(:) = 1.0_RKIND ! updated for each case below - rkTendWeights(:) = 0.0_RKIND ! updated for each case below + rkTendWeights(:) = 1.0_RKIND ! updated for each case below if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then From 413f304acdca0e12c72e7d30ff26e199f8ada9da Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 3 Oct 2023 16:27:27 -0700 Subject: [PATCH 17/36] Some more cleanup Correct a description in Registry, update comments and log writes. --- components/mpas-albany-landice/src/Registry.xml | 2 +- .../mode_forward/mpas_li_time_integration_fe_rk.F | 14 ++++---------- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index f2559440c5ef..ac119be2e7f3 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -485,7 +485,7 @@ possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS', but limited by CFL condition." /> Date: Wed, 4 Oct 2023 10:02:43 -0700 Subject: [PATCH 18/36] Fix bug in groundedSfcMassBalApplied after Runge-Kutta loop Fix bug in groundedSfcMassBalApplied after Runge-Kutta loop. It was inadvertently being assigned the enitre applied SMB field instead of just the grounded portion. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 cabe37cbdc74..e8e8855c765c 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 @@ -458,7 +458,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Finalize budget updates sfcMassBalApplied(:) = sfcMassBalAccum(:) - groundedSfcMassBalApplied(:) = sfcMassBalAccum(:) + groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) From 606cbb0e0d19f30978635012932217b8ffc91a23 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 4 Oct 2023 09:02:51 -0700 Subject: [PATCH 19/36] Calculate groundedToFloatingThickness and grounding line flux after RK Calculate the total groundedToFloatingThickness and fluxAcrossGroundingLine at the end of the RK loop. Budgets do not yet close, but are reasonably close. --- .../src/mode_forward/mpas_li_advection.F | 9 +-- .../mpas_li_time_integration_fe_rk.F | 70 +++++++++++++++++-- 2 files changed, 68 insertions(+), 11 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 7d93f114b797..5e0b90846735 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 @@ -46,7 +46,8 @@ module li_advection !-------------------------------------------------------------------- public :: li_advection_thickness_tracers, & li_layer_normal_velocity, & - li_update_geometry + li_update_geometry, & + li_grounded_to_floating !-------------------------------------------------------------------- ! @@ -631,7 +632,7 @@ subroutine li_advection_thickness_tracers(& ! Calculate volume converted from grounded to floating ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state - call grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells) + call li_grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells) ! Calculate flux across grounding line ! Do this after new thickness & mask have been calculated, including SMB/BMB @@ -2107,7 +2108,7 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers end subroutine vertical_remap - subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) + subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells) !----------------------------------------------------------------- ! input variables !----------------------------------------------------------------- @@ -2151,7 +2152,7 @@ subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, grounde endif enddo - end subroutine grounded_to_floating + end subroutine li_grounded_to_floating !*********************************************************************** 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 e8e8855c765c..513e7c599be1 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 @@ -36,6 +36,7 @@ module li_time_integration_fe_rk use li_setup use li_constants use li_mesh + use li_mask implicit none private @@ -82,6 +83,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) use li_subglacial_hydro use li_velocity use li_bedtopo + use li_mask + use li_advection, only: li_grounded_to_floating !----------------------------------------------------------------- ! input variables @@ -104,7 +107,7 @@ 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 @@ -123,13 +126,18 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) basalMassBalApplied, & groundedSfcMassBalApplied, & groundedBasalMassBalApplied, & - floatingBasalMassBalApplied + floatingBasalMassBalApplied, & + fluxAcrossGroundingLine, & + fluxAcrossGroundingLineOnCells, & + groundedToFloatingThickness real (kind=RKIND), dimension(:), allocatable :: sfcMassBalAccum, & basalMassBalAccum, & groundedSfcMassBalAccum, & groundedBasalMassBalAccum, & - floatingBasalMassBalAccum + floatingBasalMassBalAccum, & + fluxAcrossGroundingLineAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions + integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & temperaturePrev, & @@ -139,8 +147,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3d, & damage3dPrev integer, pointer :: nVertLevels - integer, pointer :: nCells - real (kind=RKIND), pointer :: deltat + integer, pointer :: nCells, nEdges + integer :: iCell1, iCell2, iEdge, theGroundedCell + integer, dimension(:), pointer :: edgeMask, cellMask + real (kind=RKIND), pointer :: deltat, config_ice_density real (kind=RKIND) :: deltatFull real (kind=RKIND), dimension(4) :: rkSubstepWeights real (kind=RKIND), dimension(4) :: rkSSPweights @@ -158,16 +168,22 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_config(liConfigs, 'config_rk_order', config_rk_order) call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages) call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) + call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) 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) call mpas_pool_get_array(meshPool, 'layerCenterSigma', layerCenterSigma) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) + call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) + + call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) + call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) allocate(temperaturePrev(nVertLevels, nCells+1)) allocate(enthalpyPrev(nVertLevels, nCells+1)) @@ -177,12 +193,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(damage3d(nVertLevels, nCells+1)) allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) + allocate(cellMaskPrev(nCells+1)) allocate(sfcMassBalAccum(nCells+1)) allocate(basalMassBalAccum(nCells+1)) allocate(groundedSfcMassBalAccum(nCells+1)) allocate(groundedBasalMassBalAccum(nCells+1)) allocate(floatingBasalMassBalAccum(nCells+1)) + allocate(fluxAcrossGroundingLineAccum(nEdges+1)) temperaturePrev(:,:) = 0.0_RKIND enthalpyPrev(:,:) = 0.0_RKIND @@ -191,6 +209,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3dPrev(:,:) = 0.0_RKIND damage3d(:,:) = 0.0_RKIND passiveTracer3dPrev(:,:) = 0.0_RKIND + cellMaskPrev(:) = 0 passiveTracer3d(:,:) = 0.0_RKIND sfcMassBalAccum(:) = 0.0_RKIND @@ -198,6 +217,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum(:) = 0.0_RKIND groundedBasalMassBalAccum(:) = 0.0_RKIND floatingBasalMassBalAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -249,6 +269,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(thermalPool, 'temperature', temperature) call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) @@ -257,6 +278,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) temperaturePrev(:,:) = temperature(:,:) enthalpyPrev(:,:) = enthalpy(:,:) layerThicknessPrev(:,:) = layerThickness(:,:) + cellMaskPrev(:) = cellMask(:) do k = 1, nVertLevels damage3dPrev(k,:) = damage(:) passiveTracer3dPrev(k,:) = passiveTracer2d(:) @@ -427,7 +449,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo endif - + + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + err = ior(err, err_tmp) + + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) + ! === Ensure damage is within bounds before velocity solve === if ( config_finalize_damage_after_advection ) then call mpas_timer_start("finalize damage") @@ -441,13 +469,17 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) + ! update budgets sfcMassBalAccum = sfcMassBalAccum + rkTendWeights(rkStage) * sfcMassBalApplied groundedSfcMassBalAccum = groundedSfcMassBalAccum + rkTendWeights(rkStage) * groundedSfcMassBalApplied basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied - + fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine + ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) @@ -457,12 +489,34 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo ! Finalize budget updates + ! Calculate volume converted from grounded to floating + ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + call li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) sfcMassBalApplied(:) = sfcMassBalAccum(:) groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) + fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) + fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND + do iEdge = 1, nEdges + 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 + theGroundedCell = iCell1 + else + theGroundedCell = iCell2 + endif + if ( thickness(theGroundedCell) > 0.0_RKIND ) then + fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + & + fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density ! adjust to correct units + else + fluxAcrossGroundingLineOnCells(theGroundedCell) = 0.0_RKIND + endif + endif + enddo ! edges ! Reset time step to full length after RK loop deltat = deltatFull @@ -515,11 +569,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(damage3d) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) + deallocate(cellMaskPrev) deallocate(sfcMassBalAccum) deallocate(basalMassBalAccum) deallocate(groundedSfcMassBalAccum) deallocate(groundedBasalMassBalAccum) deallocate(floatingBasalMassBalAccum) + deallocate(fluxAcrossGroundingLineAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From 35fe8fa5fde6b754afa61ef253aa40ce27086fb9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 4 Oct 2023 21:09:13 -0600 Subject: [PATCH 20/36] Reset layerThickness when using config_restore_thickness_after_advection Reset layerThickness when using config_restore_thickness_after_advection. Restoring thickness alone is not sufficient when using Runge-Kutta time integration because the advected layerThickness gets propagated through all stages of the Runge-Kutta loop, which leads to advection and surface/basal mass balances causing thickness change. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 4 ++++ 1 file changed, 4 insertions(+) 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 513e7c599be1..6fd1eae52c59 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 @@ -985,6 +985,7 @@ subroutine advection_solver(domain, err) real (kind=RKIND), dimension(:), pointer :: thicknessOld real (kind=RKIND), dimension(:), pointer :: thicknessNew real (kind=RKIND), dimension(:), pointer :: thickness + real (kind=RKIND), dimension(:,:), pointer :: layerThickness real (kind=RKIND), dimension(:,:), pointer :: temperature real (kind=RKIND), dimension(:,:), pointer :: waterFrac @@ -1173,7 +1174,10 @@ subroutine advection_solver(domain, err) if ( config_restore_thickness_after_advection ) then call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'thicknessOld', thicknessOld) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) thickness(:) = thicknessOld(:) + + call li_calculate_layerThickness(meshPool, thickness, layerThickness) endif call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) From 07fdc99acad95f26694085295b3ebd9f390f7986 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 6 Oct 2023 16:04:44 -0600 Subject: [PATCH 21/36] Update stresses and strain rates after calving Update stresses and strain rates after calving. This allows restart test with damage to pass validation, but may not be the most rigorous fix. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 5 +++++ 1 file changed, 5 insertions(+) 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 6fd1eae52c59..162749c1d81a 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 @@ -541,6 +541,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") + ! Update stresses etc + ! === Solve Velocity ===================== + call li_velocity_solve(domain, solveVelo=.false., err=err_tmp) + err = ior(err, err_tmp) + ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice From f875198f61a86c566f724c7709293f6bbf4b22d5 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 6 Oct 2023 17:45:32 -0600 Subject: [PATCH 22/36] Update thickness and tracer halos before velocity solve Update thickness and tracer halos before velocity solve. This allows eismint2 decomposition test to pass validation. --- .../mode_forward/mpas_li_time_integration_fe_rk.F | 13 +++++++++++++ 1 file changed, 13 insertions(+) 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 162749c1d81a..ace741af863f 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 @@ -480,6 +480,19 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine + ! Halo updates + call mpas_timer_start("halo updates") + + call mpas_dmpar_field_halo_exch(domain, 'thickness') + call mpas_dmpar_field_halo_exch(domain, 'temperature') + call mpas_dmpar_field_halo_exch(domain, 'waterFrac') + call mpas_dmpar_field_halo_exch(domain, 'enthalpy') + call mpas_dmpar_field_halo_exch(domain, 'damage') + call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + + call mpas_timer_stop("halo updates") + + ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) From 20fdec7651e439895d889a9dd1a0521242e1c2cd Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 9 Oct 2023 17:21:37 -0600 Subject: [PATCH 23/36] Update temperature and waterFrac in RK loop, but not enthalpy Update temperature and waterFrac in RK loop, but not enthalpy. Updating enthalpy resulted in negative temperatures, which caused an error in li_thermal. Updating temperature and waterFrac instead after enthalpy is advected allows enthalpy tests in full_integration suite to pass execution and validation. --- .../mpas_li_time_integration_fe_rk.F | 48 ++++++------------- 1 file changed, 15 insertions(+), 33 deletions(-) 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 ace741af863f..171eeb1199fc 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 @@ -119,7 +119,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) integer :: rkStage, iCell, iTracer, k real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & temperature, & - enthalpy, & waterFrac real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d real (kind=RKIND), dimension(:), pointer :: sfcMassBalApplied, & @@ -141,7 +140,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & temperaturePrev, & - enthalpyPrev, & + waterFracPrev, & passiveTracer3d, & passiveTracer3dPrev, & damage3d, & @@ -186,7 +185,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) allocate(temperaturePrev(nVertLevels, nCells+1)) - allocate(enthalpyPrev(nVertLevels, nCells+1)) + allocate(waterFracPrev(nVertLevels, nCells+1)) allocate(layerThicknessPrev(nVertLevels, nCells+1)) allocate(layerThicknessTmp(nVertLevels, nCells+1)) allocate(damage3dPrev(nVertLevels, nCells+1)) @@ -203,7 +202,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) temperaturePrev(:,:) = 0.0_RKIND - enthalpyPrev(:,:) = 0.0_RKIND layerThicknessPrev(:,:) = 0.0_RKIND layerThicknessTmp(:,:) = 0.0_RKIND damage3dPrev(:,:) = 0.0_RKIND @@ -272,11 +270,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(thermalPool, 'temperature', temperature) - call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) ! Save relevant fields before RK loop, to be used in update at the end temperaturePrev(:,:) = temperature(:,:) - enthalpyPrev(:,:) = enthalpy(:,:) + waterFracPrev(:,:) = waterFrac(:,:) layerThicknessPrev(:,:) = layerThickness(:,:) cellMaskPrev(:) = cellMask(:) do k = 1, nVertLevels @@ -384,7 +381,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'damage', damage) call mpas_pool_get_array(thermalPool, 'temperature', temperature) - call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) layerThicknessTmp(:,:) = layerThickness(:,:) @@ -392,32 +388,19 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) thickness = sum(layerThickness, 1) - if (trim(config_thermal_solver) == 'enthalpy') then - where (layerThickness(:,:) > 0.0_RKIND) - enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyPrev(:,:) * layerThicknessPrev(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * & - layerThicknessTmp(:,:) ) / layerThickness(:,:) - ! elsewhere, what? - end where - ! given the enthalpy, compute the temperature and water fraction + if (trim(config_thermal_solver) .ne. 'none') then do iCell = 1, nCells - - call li_enthalpy_to_temperature_kelvin(& - layerCenterSigma, & - thickness(iCell), & - enthalpy(:,iCell), & - temperature(:,iCell), & - waterFrac(:,iCell)) - + do k = 1, nVertLevels + if (layerThickness(k,iCell) > 0.0_RKIND) then + temperature(k,iCell) = ( rkSSPweights(rkStage) * temperaturePrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * temperature(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + waterFrac(k,iCell) = ( rkSSPweights(rkStage) * waterFracPrev(k,iCell) * layerThicknessPrev(k,iCell) + & + (1.0_RKIND - rkSSPweights(rkStage)) * waterFrac(k,iCell) * & + layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell) + endif + enddo enddo - - elseif (trim(config_thermal_solver) == 'temperature') then - where (layerThickness(:,:) > 0.0_RKIND) - temperature(:,:) = ( rkSSPweights(rkStage) * temperaturePrev(:,:) * layerThicknessPrev(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * & - layerThicknessTmp(:,:) ) / layerThickness(:,:) - ! elsewhere, what? - end where endif if (config_calculate_damage) then @@ -486,7 +469,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'thickness') call mpas_dmpar_field_halo_exch(domain, 'temperature') call mpas_dmpar_field_halo_exch(domain, 'waterFrac') - call mpas_dmpar_field_halo_exch(domain, 'enthalpy') call mpas_dmpar_field_halo_exch(domain, 'damage') call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') @@ -580,7 +562,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) endif deallocate(temperaturePrev) - deallocate(enthalpyPrev) + deallocate(waterFracPrev) deallocate(layerThicknessPrev) deallocate(layerThicknessTmp) deallocate(damage3dPrev) From fbb6d6def295e9c95114e2d97758c6283c90640e Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 27 Oct 2023 11:06:54 -0600 Subject: [PATCH 24/36] Simple clean-up after code review Simple clean-up after code review, including moving call to mask calculation, initializing waterFracPrev to zeroes, adding return statements after error messages. Also add a halo update on cellMask and edgeMask before calculating groundedToFloatingThickness and fluxAcrossGroundingLineOnCells. --- .../mpas_li_time_integration_fe_rk.F | 21 ++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) 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 171eeb1199fc..7e82021c2eb0 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 @@ -202,6 +202,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) temperaturePrev(:,:) = 0.0_RKIND + waterFracPrev(:,:) = 0.0_RKIND layerThicknessPrev(:,:) = 0.0_RKIND layerThicknessTmp(:,:) = 0.0_RKIND damage3dPrev(:,:) = 0.0_RKIND @@ -293,7 +294,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then - config_rk_order = 1 nRKstages = 1 rkSubstepWeights(:) = 1.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & @@ -344,12 +344,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = 1 call mpas_log_write('config_rk3_stages must 3 or 4', & messageType=MPAS_LOG_ERR) + return endif else err = 1 call mpas_log_write('config_time_integration = ' // trim(config_time_integration) & // ' is not supported with config_rk_order = $i', & intArgs=(/config_rk_order/), messageType=MPAS_LOG_ERR) + return endif ! *** Start RK loop *** do rkStage = 1, nRKstages @@ -387,6 +389,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) thickness = sum(layerThickness, 1) + ! Calculate masks after updating thickness + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + err = ior(err, err_tmp) if (trim(config_thermal_solver) .ne. 'none') then do iCell = 1, nCells @@ -433,12 +438,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) endif - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) - err = ior(err, err_tmp) - - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) - call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) - ! === Ensure damage is within bounds before velocity solve === if ( config_finalize_damage_after_advection ) then call mpas_timer_start("finalize damage") @@ -486,6 +485,14 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Finalize budget updates ! Calculate volume converted from grounded to floating ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'cellMask') + call mpas_dmpar_field_halo_exch(domain, 'edgeMask') + call mpas_timer_stop("halo updates") + call li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) sfcMassBalApplied(:) = sfcMassBalAccum(:) groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) From d7eac5ca41285b1962574191c98e9182cc443e49 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 30 Oct 2023 09:51:13 -0700 Subject: [PATCH 25/36] Initialize RK weights to large negative values Initialize the various RK weights (rkSSPweights, rkTendWeights, rkSubstepWeights) to large negative values which are overwritten for each case. This will make it obvious if a weight is being used when it should not be. --- .../mpas_li_time_integration_fe_rk.F | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) 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 7e82021c2eb0..bac6ae785b51 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 @@ -289,25 +289,28 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! by one index. e.g., the coefficients for RK4 are usually written ! (0, 1/2, 1/2, 1), while we use (1/2, 1/2, 1). The last entry of 1.0 ! is simply for ease of implementation. - rkSSPweights(:) = 1.0_RKIND ! updated for each case below - rkTendWeights(:) = 1.0_RKIND ! updated for each case below + rkSSPweights(:) = -9999.0_RKIND ! Initialized to this value to make it obvious if + ! a weight is used that should not be. Appropriate + ! weights are updated for each case below + rkTendWeights(:) = -9999.0_RKIND + rkSubstepWeights(:) = -9999.0_RKIND if ( (trim(config_time_integration) == 'forward_euler') .or. & ((trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 1)) ) then nRKstages = 1 - rkSubstepWeights(:) = 1.0_RKIND + rkSSPweights(1) = 1.0_RKIND + rkTendWeights(1) = 1.0_RKIND + rkSubstepWeights(1) = 1.0_RKIND elseif ( (trim(config_time_integration) == 'runge_kutta') .and. & (config_rk_order == 2) ) then ! use Strong Stability Preserving RK2. Could also ! add config option to use standard endpoint or midpoint methods, but ! these are algorithmically more complex and less suitable for our domains nRKstages = 2 - rkSubstepWeights(:) = 1.0_RKIND + rkSubstepWeights(1:2) = 1.0_RKIND rkSSPweights(1) = 1.0_RKIND rkSSPweights(2) = 0.5_RKIND - rkSSPweights(3) = 0.0_RKIND - rkSSPweights(4) = 0.0_RKIND rkTendWeights(1) = 0.5_RKIND rkTendWeights(2) = 0.5_RKIND @@ -316,12 +319,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) if (config_rk3_stages == 3) then ! use three-stage Strong Stability Preserving RK3 nRKstages = 3 - rkSubstepWeights(:) = 1.0_RKIND + rkSubstepWeights(1:3) = 1.0_RKIND rkSSPweights(1) = 1.0_RKIND rkSSPweights(2) = 3.0_RKIND / 4.0_RKIND rkSSPweights(3) = 1.0_RKIND / 3.0_RKIND - rkSSPweights(4) = 0.0_RKIND rkTendWeights(1) = 1.0_RKIND / 6.0_RKIND rkTendWeights(2) = 1.0_RKIND / 6.0_RKIND From 3bcd1d9793ed763d335731244b1ac9e64bb5bb44 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 30 Oct 2023 12:31:28 -0600 Subject: [PATCH 26/36] Restore calving front within RK loop Restore calving front within RK loop to prevent velocity solutions on extended ice geometry. --- .../mode_forward/mpas_li_time_integration_fe_rk.F | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) 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 bac6ae785b51..60da47d13c71 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 @@ -109,7 +109,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) integer :: err_tmp type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool, velocityPool - logical, pointer :: config_restore_calving_front + logical, pointer :: config_restore_calving_front, & + config_restore_calving_front_prevent_retreat logical, pointer :: config_calculate_damage logical, pointer :: config_finalize_damage_after_advection character (len=StrKIND), pointer :: config_thickness_advection @@ -160,6 +161,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err_tmp = 0 call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) + call mpas_pool_get_config(liConfigs, 'config_restore_calving_front_prevent_retreat', config_restore_calving_front_prevent_retreat) 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) @@ -475,6 +477,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("halo updates") + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif ! Update velocity for each RK step ! === Solve Velocity ===================== @@ -531,8 +538,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call li_calve_ice(domain, err_tmp) err = ior(err, err_tmp) - if (config_restore_calving_front) then - ! restore the calving front to its initial position; calving options are ignored + if (config_restore_calving_front .and. config_restore_calving_front_prevent_retreat) then + ! restore the calving front to its initial position before velocity solve. call li_restore_calving_front(domain, err_tmp) err = ior(err, err_tmp) endif From 13f0908c0336495ec14580f3c735edf882d2dab1 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Mon, 30 Oct 2023 12:55:20 -0600 Subject: [PATCH 27/36] Remove support for config_rk_order = 1 Remove support for config_rk_order = 1 to avoid unnecessary redundancy with config_time_integration = forward_euler while maintaining backward compatibility. --- components/mpas-albany-landice/src/Registry.xml | 4 ++-- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 4 +--- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index ac119be2e7f3..9eae3fb81f27 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -489,8 +489,8 @@ possible_values="'forward_euler' or 'runge_kutta'" /> Date: Mon, 30 Oct 2023 19:25:01 -0600 Subject: [PATCH 28/36] Include calving in Runge-Kutta loop Include calving in Runge-Kutta loop, rather than implementing it afterwards. --- .../mpas_li_time_integration_fe_rk.F | 52 +++++++------------ 1 file changed, 20 insertions(+), 32 deletions(-) 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 5dda3f792766..97e50ba9ee9c 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 @@ -127,6 +127,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & + calvingThickness, & fluxAcrossGroundingLine, & fluxAcrossGroundingLineOnCells, & groundedToFloatingThickness @@ -135,6 +136,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum, & groundedBasalMassBalAccum, & floatingBasalMassBalAccum, & + calvingThicknessAccum, & fluxAcrossGroundingLineAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection @@ -201,6 +203,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(groundedSfcMassBalAccum(nCells+1)) allocate(groundedBasalMassBalAccum(nCells+1)) allocate(floatingBasalMassBalAccum(nCells+1)) + allocate(calvingThicknessAccum(nCells+1)) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) temperaturePrev(:,:) = 0.0_RKIND @@ -218,6 +221,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum(:) = 0.0_RKIND groundedBasalMassBalAccum(:) = 0.0_RKIND floatingBasalMassBalAccum(:) = 0.0_RKIND + calvingThicknessAccum(:) = 0.0_RKIND fluxAcrossGroundingLineAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== @@ -375,6 +379,18 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advect thickness and tracers") + ! ==== Apply iceberg calving ==== + call mpas_timer_start("calve_ice") + call li_calve_ice(domain, err_tmp) + err = ior(err, err_tmp) + + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif + call mpas_timer_stop("calve_ice") + ! If using SSP RK, then update thickness and tracers incrementally. ! For first RK stage, thickness and tracer updates above are sufficient. ! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update. @@ -453,6 +469,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) call mpas_pool_get_array(geometryPool, 'thickness', thickness) @@ -462,6 +479,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied + calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine ! Halo updates @@ -475,12 +493,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("halo updates") - if (config_restore_calving_front) then - ! restore the calving front to its initial position before velocity solve. - call li_restore_calving_front(domain, err_tmp) - err = ior(err, err_tmp) - endif - ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) @@ -506,6 +518,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) + calvingThickness(:) = calvingThicknessAccum(:) fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND @@ -529,32 +542,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Reset time step to full length after RK loop deltat = deltatFull -! === Calve ice ======================== - call mpas_timer_start("calve_ice") - - ! ice calving - call li_calve_ice(domain, err_tmp) - err = ior(err, err_tmp) - - if (config_restore_calving_front .and. config_restore_calving_front_prevent_retreat) then - ! restore the calving front to its initial position before velocity solve. - call li_restore_calving_front(domain, err_tmp) - err = ior(err, err_tmp) - endif - - call mpas_timer_stop("calve_ice") - - call mpas_timer_start("halo updates") - call mpas_dmpar_field_halo_exch(domain, 'cellMask') - call mpas_dmpar_field_halo_exch(domain, 'edgeMask') - call mpas_dmpar_field_halo_exch(domain, 'vertexMask') - call mpas_timer_stop("halo updates") - - ! Update stresses etc - ! === Solve Velocity ===================== - call li_velocity_solve(domain, solveVelo=.false., err=err_tmp) - err = ior(err, err_tmp) - ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice @@ -589,6 +576,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(groundedSfcMassBalAccum) deallocate(groundedBasalMassBalAccum) deallocate(floatingBasalMassBalAccum) + deallocate(calvingThicknessAccum) deallocate(fluxAcrossGroundingLineAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From add0633d37bb0c4a799b37bcc1ff41f476e21189 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 1 Nov 2023 19:26:49 -0700 Subject: [PATCH 29/36] Calculate layerThickness after calving in RK loop Calculate layerThickness after calving in RK loop, which is necessary when calculating weighted averages for RK integration. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 2 ++ 1 file changed, 2 insertions(+) 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 97e50ba9ee9c..4b5ee90956c0 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 @@ -403,6 +403,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(thermalPool, 'temperature', temperature) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + ! get layerThickness after calving + call li_calculate_layerThickness(meshPool, thickness, layerThickness) layerThicknessTmp(:,:) = layerThickness(:,:) layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) From 7073b89680ca181c1746d509f1df8d998d151358 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 6 Nov 2023 10:51:05 -0800 Subject: [PATCH 30/36] Revert "Calculate layerThickness after calving in RK loop" This reverts commit 5e23b98c14fdf9cb2f1277f4da0b0e6de8bdb9cb. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 2 -- 1 file changed, 2 deletions(-) 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 4b5ee90956c0..97e50ba9ee9c 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 @@ -403,8 +403,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(thermalPool, 'temperature', temperature) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) - ! get layerThickness after calving - call li_calculate_layerThickness(meshPool, thickness, layerThickness) layerThicknessTmp(:,:) = layerThickness(:,:) layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) From fb0315ff4e4c8b0e5bfd7f6d0c2598cd77aeea97 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 6 Nov 2023 10:51:23 -0800 Subject: [PATCH 31/36] Revert "Include calving in Runge-Kutta loop" This reverts commit 813ee0555ec6239f6f0f13dd1a8e353e1a10721b. After testing and discussion, it is apparent that including calving within the Runge-Kutta loop is not desirable because it changes the domain between stages of the Runge-Kutta integration. --- .../mpas_li_time_integration_fe_rk.F | 52 ++++++++++++------- 1 file changed, 32 insertions(+), 20 deletions(-) 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 97e50ba9ee9c..5dda3f792766 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 @@ -127,7 +127,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalApplied, & groundedBasalMassBalApplied, & floatingBasalMassBalApplied, & - calvingThickness, & fluxAcrossGroundingLine, & fluxAcrossGroundingLineOnCells, & groundedToFloatingThickness @@ -136,7 +135,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum, & groundedBasalMassBalAccum, & floatingBasalMassBalAccum, & - calvingThicknessAccum, & fluxAcrossGroundingLineAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection @@ -203,7 +201,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(groundedSfcMassBalAccum(nCells+1)) allocate(groundedBasalMassBalAccum(nCells+1)) allocate(floatingBasalMassBalAccum(nCells+1)) - allocate(calvingThicknessAccum(nCells+1)) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) temperaturePrev(:,:) = 0.0_RKIND @@ -221,7 +218,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum(:) = 0.0_RKIND groundedBasalMassBalAccum(:) = 0.0_RKIND floatingBasalMassBalAccum(:) = 0.0_RKIND - calvingThicknessAccum(:) = 0.0_RKIND fluxAcrossGroundingLineAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== @@ -379,18 +375,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advect thickness and tracers") - ! ==== Apply iceberg calving ==== - call mpas_timer_start("calve_ice") - call li_calve_ice(domain, err_tmp) - err = ior(err, err_tmp) - - if (config_restore_calving_front) then - ! restore the calving front to its initial position before velocity solve. - call li_restore_calving_front(domain, err_tmp) - err = ior(err, err_tmp) - endif - call mpas_timer_stop("calve_ice") - ! If using SSP RK, then update thickness and tracers incrementally. ! For first RK stage, thickness and tracer updates above are sufficient. ! Likewise, for the 4-stage SSP RK3 the last stage is just a forward euler update. @@ -469,7 +453,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'basalMassBalApplied', basalMassBalApplied) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) - call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) call mpas_pool_get_array(geometryPool, 'thickness', thickness) @@ -479,7 +462,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) basalMassBalAccum = basalMassBalAccum + rkTendWeights(rkStage) * basalMassBalApplied groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied - calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine ! Halo updates @@ -493,6 +475,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("halo updates") + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif + ! Update velocity for each RK step ! === Solve Velocity ===================== call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) @@ -518,7 +506,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) - calvingThickness(:) = calvingThicknessAccum(:) fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND @@ -542,6 +529,32 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Reset time step to full length after RK loop deltat = deltatFull +! === Calve ice ======================== + call mpas_timer_start("calve_ice") + + ! ice calving + call li_calve_ice(domain, err_tmp) + err = ior(err, err_tmp) + + if (config_restore_calving_front .and. config_restore_calving_front_prevent_retreat) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif + + call mpas_timer_stop("calve_ice") + + call mpas_timer_start("halo updates") + call mpas_dmpar_field_halo_exch(domain, 'cellMask') + call mpas_dmpar_field_halo_exch(domain, 'edgeMask') + call mpas_dmpar_field_halo_exch(domain, 'vertexMask') + call mpas_timer_stop("halo updates") + + ! Update stresses etc + ! === Solve Velocity ===================== + call li_velocity_solve(domain, solveVelo=.false., err=err_tmp) + err = ior(err, err_tmp) + ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice @@ -576,7 +589,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(groundedSfcMassBalAccum) deallocate(groundedBasalMassBalAccum) deallocate(floatingBasalMassBalAccum) - deallocate(calvingThicknessAccum) deallocate(fluxAcrossGroundingLineAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta From 51802e1bca96f70da0a72ab6786116b88dfa9c67 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 6 Nov 2023 11:47:15 -0800 Subject: [PATCH 32/36] Average fluxAcrossGroundingLineOnCells in time integration module Calculate fluxAcrossGroundingLineOnCells in advection module for each RK stage and then take weighted average after RK integration instead of calcuating fluxAcrossGroundingLineOnCells on final state. This is more consistent with how fluxAcrossGroundingLine (on edges) is treated. --- .../src/mode_forward/mpas_li_advection.F | 4 -- .../mpas_li_time_integration_fe_rk.F | 37 ++++++++----------- 2 files changed, 15 insertions(+), 26 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 5e0b90846735..56b018cdd399 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 @@ -630,10 +630,6 @@ subroutine li_advection_thickness_tracers(& call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state - call li_grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells) - ! Calculate flux across grounding line ! Do this after new thickness & mask have been calculated, including SMB/BMB fluxAcrossGroundingLine(:) = 0.0_RKIND 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 5dda3f792766..6c9874f26486 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 @@ -135,7 +135,8 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedSfcMassBalAccum, & groundedBasalMassBalAccum, & floatingBasalMassBalAccum, & - fluxAcrossGroundingLineAccum + fluxAcrossGroundingLineAccum, & + fluxAcrossGroundingLineOnCellsAccum real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & @@ -202,6 +203,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(groundedBasalMassBalAccum(nCells+1)) allocate(floatingBasalMassBalAccum(nCells+1)) allocate(fluxAcrossGroundingLineAccum(nEdges+1)) + allocate(fluxAcrossGroundingLineOnCellsAccum(nCells+1)) temperaturePrev(:,:) = 0.0_RKIND waterFracPrev(:,:) = 0.0_RKIND @@ -219,6 +221,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedBasalMassBalAccum(:) = 0.0_RKIND floatingBasalMassBalAccum(:) = 0.0_RKIND fluxAcrossGroundingLineAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCellsAccum(:) = 0.0_RKIND ! === Prepare for advection (including CFL checks) =========== ! This has to come first currently, because it sets the time step! @@ -454,6 +457,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBalApplied', floatingBasalMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBalApplied', groundedBasalMassBalApplied) call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) call mpas_pool_get_array(geometryPool, 'thickness', thickness) ! update budgets @@ -463,6 +467,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) groundedBasalMassBalAccum = groundedBasalMassBalAccum + rkTendWeights(rkStage) * groundedBasalMassBalApplied floatingBasalMassBalAccum = floatingBasalMassBalAccum + rkTendWeights(rkStage) * floatingBasalMassBalApplied fluxAcrossGroundingLineAccum = fluxAcrossGroundingLineAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLine + fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells ! Halo updates call mpas_timer_start("halo updates") @@ -490,8 +495,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) enddo ! Finalize budget updates - ! Calculate volume converted from grounded to floating - ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state + ! Update masks after RK integration + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + err = ior(err, err_tmp) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) @@ -500,32 +507,18 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'edgeMask') call mpas_timer_stop("halo updates") + ! Calculate volume converted from grounded to floating + ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state call li_grounded_to_floating(cellMaskPrev, cellMask, thickness, groundedToFloatingThickness, nCells) + sfcMassBalApplied(:) = sfcMassBalAccum(:) groundedSfcMassBalApplied(:) = groundedSfcMassBalAccum(:) basalMassBalApplied(:) = basalMassBalAccum(:) groundedBasalMassBalApplied(:) = groundedBasalMassBalAccum(:) floatingBasalMassBalApplied(:) = floatingBasalMassBalAccum(:) fluxAcrossGroundingLine(:) = fluxAcrossGroundingLineAccum(:) - - fluxAcrossGroundingLineOnCells(:) = 0.0_RKIND - do iEdge = 1, nEdges - 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 - theGroundedCell = iCell1 - else - theGroundedCell = iCell2 - endif - if ( thickness(theGroundedCell) > 0.0_RKIND ) then - fluxAcrossGroundingLineOnCells(theGroundedCell) = fluxAcrossGroundingLineOnCells(theGroundedCell) + & - fluxAcrossGroundingLine(iEdge) / thickness(theGroundedCell) * config_ice_density ! adjust to correct units - else - fluxAcrossGroundingLineOnCells(theGroundedCell) = 0.0_RKIND - endif - endif - enddo ! edges + fluxAcrossGroundingLineOnCells(:) = fluxAcrossGroundingLineOnCellsAccum(:) + ! Reset time step to full length after RK loop deltat = deltatFull From ca28f5bc3fece5ce49d39d425ed12f66bfe32310 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Tue, 7 Nov 2023 19:44:48 -0800 Subject: [PATCH 33/36] Add logic to control whether velocity is solved before and/or after calving Add logic to control whether velocity is solved before and/or after calving. This allows both for backward compatibility and for the abilility to use self-consistent velocty, stress, and strain-rate fields when calculating calving. --- .../mpas-albany-landice/src/Registry.xml | 4 ++ .../src/mode_forward/mpas_li_calving.F | 44 ++++++++++++++++--- .../mpas_li_time_integration_fe_rk.F | 41 +++++++++-------- 3 files changed, 64 insertions(+), 25 deletions(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 9eae3fb81f27..d5507295fe0a 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -292,6 +292,10 @@ description="If true, then distribute unablatedVolumeDynCell among dynamic neighbors when converting ablation velocity to ablation thickness. This should only be used as a clean-up measure, while limiting the timestep based on ablation velocity should be used as the primary method of getting accurate ablation thickness from ablation velocity. If you choose to set config_adaptive_timestep_calvingCFL_fraction much larger than 1.0 (which is not recommended), setting this option to true usually results in more accurate calving behavior. " possible_values=".true. or .false." /> + diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 36302a1a78c3..304ed0b0349c 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F @@ -74,7 +74,7 @@ module li_calving !> (3) Calve ice based on an ice thickness threshold !----------------------------------------------------------------------- - subroutine li_calve_ice(domain, err) + subroutine li_calve_ice(domain, err, solveVeloAfterCalving) use li_advection @@ -92,7 +92,7 @@ subroutine li_calve_ice(domain, err) ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - + logical, intent(out), optional :: solveVeloAfterCalving !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -111,13 +111,14 @@ subroutine li_calve_ice(domain, err) logical, pointer :: config_apply_calving_mask real(kind=RKIND), pointer :: config_calving_timescale - integer, pointer :: nCells + integer, pointer :: nCells, nCellsSolve integer, pointer :: config_number_of_blocks real (kind=RKIND), pointer :: deltat !< time step (s) integer, dimension(:), pointer :: & - indexToCellID ! list of global cell IDs + indexToCellID, & ! list of global cell IDs + cellMask real (kind=RKIND) :: & calvingFraction ! fraction of ice that calves in each column; depends on calving_timescale @@ -136,14 +137,22 @@ subroutine li_calve_ice(domain, err) real (kind=RKIND), dimension(:), pointer :: calvingVelocity + real (kind=RKIND), dimension(:), pointer :: areaCell + type (field1dReal), pointer :: originalThicknessField real (kind=RKIND), dimension(:), pointer :: originalThickness + type (field1dInteger), pointer :: cellMaskTemporaryField + integer :: iCell + real (kind=RKIND), parameter :: smallNumber = 1.0e-6_RKIND + integer :: err_tmp + real (kind=RKIND), dimension(1) :: calvingSumLocal, calvingSumGlobal + err = 0 err_tmp = 0 @@ -155,6 +164,15 @@ subroutine li_calve_ice(domain, err) call mpas_pool_get_config(liConfigs, 'config_data_calving', config_data_calving) call mpas_pool_get_config(liConfigs, 'config_number_of_blocks', config_number_of_blocks) + if (present(solveVeloAfterCalving)) then + call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) + call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField) + call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + ! Store cellMask prior to calving + cellMaskTemporaryField % array(:) = cellMask(:) + endif + ! Zero calvingThickness here instead of or in addition to in individual subroutines. ! This is necessary when using damage threshold calving with other calving ! routines. Some individual routines still set calvingThickness to zero, but @@ -206,7 +224,7 @@ subroutine li_calve_ice(domain, err) ! However, the eigencalving method requires multiple applications of the calvingThickness ! to the thickness. So the simplest method to apply data calving is to store the old ! thickness and then set it back when we are done. - if (config_data_calving) then + if ( config_data_calving .or. present(solveVeloAfterCalving) ) then call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool) call mpas_pool_get_field(scratchPool, 'workCell2', originalThicknessField) call mpas_allocate_scratch_field(originalThicknessField, single_block_in = .false.) @@ -337,8 +355,22 @@ subroutine li_calve_ice(domain, err) block => block % next end do + if (present(solveVeloAfterCalving)) then + call mpas_pool_get_array(meshPool, 'areaCell', areaCell) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + calvingSumLocal(1) = sum(calvingThickness(1:nCellsSolve) * areaCell(1:nCellsSolve) * & + real(li_mask_is_dynamic_ice_int(cellMaskTemporaryField % array(1:nCellsSolve)), RKIND)) + call mpas_dmpar_sum_real_array(domain % dminfo, 1, calvingSumLocal(1), calvingSumGlobal(1)) + if (calvingSumGlobal(1) > smallNumber) then + solveVeloAfterCalving = .true. + else + solveVeloAfterCalving = .false. + endif + call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.) + endif - if (config_data_calving) then + if ( config_data_calving .or. present(solveVeloAfterCalving) ) then call mpas_deallocate_scratch_field(originalThicknessField, single_block_in=.false.) endif 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 6c9874f26486..4858ef0abbd2 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 @@ -96,12 +96,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) type (domain_type), intent(inout) :: & domain !< Input/Output: domain object - !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- integer, intent(out) :: err !< Output: error flag - !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- @@ -113,6 +111,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) config_restore_calving_front_prevent_retreat logical, pointer :: config_calculate_damage logical, pointer :: config_finalize_damage_after_advection + logical, pointer :: config_update_velocity_before_calving character (len=StrKIND), pointer :: config_thickness_advection character (len=StrKIND), pointer :: config_thermal_solver character (len=StrKIND), pointer :: config_time_integration @@ -157,9 +156,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) real (kind=RKIND), dimension(4) :: rkSSPweights real (kind=RKIND), dimension(4) :: rkTendWeights ! Weights used for calculating budget terms integer :: nRKstages + logical :: solveVeloAfterCalving err = 0 err_tmp = 0 + solveVeloAfterCalving = .false. call mpas_pool_get_config(liConfigs, 'config_restore_calving_front', config_restore_calving_front) call mpas_pool_get_config(liConfigs, 'config_restore_calving_front_prevent_retreat', config_restore_calving_front_prevent_retreat) @@ -171,6 +172,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_pool_get_config(liConfigs, 'config_rk3_stages', config_rk3_stages) call mpas_pool_get_config(liConfigs, 'config_time_integration', config_time_integration) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) + call mpas_pool_get_config(liConfigs, 'config_update_velocity_before_calving', config_update_velocity_before_calving) call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(domain % blocklist % structs, 'thermal', thermalPool) @@ -480,17 +482,20 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_stop("halo updates") - if (config_restore_calving_front) then - ! restore the calving front to its initial position before velocity solve. - call li_restore_calving_front(domain, err_tmp) - err = ior(err, err_tmp) - endif - ! Update velocity for each RK step ! === Solve Velocity ===================== - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) - err = ior(err, err_tmp) + if ( config_update_velocity_before_calving .or. ( (.not. config_update_velocity_before_calving) & + .and. (rkStage < nRKstages) ) ) then + + if (config_restore_calving_front) then + ! restore the calving front to its initial position before velocity solve. + call li_restore_calving_front(domain, err_tmp) + err = ior(err, err_tmp) + endif + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + endif ! *** end RK loop *** enddo @@ -526,15 +531,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_start("calve_ice") ! ice calving - call li_calve_ice(domain, err_tmp) + call li_calve_ice(domain, err_tmp, solveVeloAfterCalving) err = ior(err, err_tmp) - - if (config_restore_calving_front .and. config_restore_calving_front_prevent_retreat) then + if (config_restore_calving_front) then ! restore the calving front to its initial position before velocity solve. call li_restore_calving_front(domain, err_tmp) err = ior(err, err_tmp) endif - call mpas_timer_stop("calve_ice") call mpas_timer_start("halo updates") @@ -543,11 +546,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") - ! Update stresses etc - ! === Solve Velocity ===================== - call li_velocity_solve(domain, solveVelo=.false., err=err_tmp) - err = ior(err, err_tmp) - +! === Solve velocity for final state ===================== + if (solveVeloAfterCalving) then + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + endif ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice From 4cda6487b7f1a516359639e5a7be97ffa9d27b22 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 8 Nov 2023 12:02:35 -0700 Subject: [PATCH 34/36] Make config_update_velocity_before_calving = .false. the default Make config_update_velocity_before_calving = .false. the default. This is not necessarily the desired behavior for production runs, but it is necessary for baseline comparison when using forward Euler time integration. --- components/mpas-albany-landice/src/Registry.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index d5507295fe0a..279d7fb3c0a0 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -292,7 +292,7 @@ description="If true, then distribute unablatedVolumeDynCell among dynamic neighbors when converting ablation velocity to ablation thickness. This should only be used as a clean-up measure, while limiting the timestep based on ablation velocity should be used as the primary method of getting accurate ablation thickness from ablation velocity. If you choose to set config_adaptive_timestep_calvingCFL_fraction much larger than 1.0 (which is not recommended), setting this option to true usually results in more accurate calving behavior. " possible_values=".true. or .false." /> - From 9ac9af25e56c2b637f7337eab3f456a92e55b14d Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 8 Nov 2023 12:05:09 -0700 Subject: [PATCH 35/36] Always update velocity after calving if it was not updated before calving. Always update velocity after calving if it was not updated before calving. --- .../src/mode_forward/mpas_li_time_integration_fe_rk.F | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) 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 4858ef0abbd2..1ae17ed878db 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 @@ -531,7 +531,13 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_start("calve_ice") ! ice calving - call li_calve_ice(domain, err_tmp, solveVeloAfterCalving) + if ( config_update_velocity_before_calving ) then + call li_calve_ice(domain, err_tmp, solveVeloAfterCalving) + else + call li_calve_ice(domain, err_tmp) + solveVeloAfterCalving = .true. + endif + err = ior(err, err_tmp) if (config_restore_calving_front) then ! restore the calving front to its initial position before velocity solve. From e2a2f0b9389ed3810d0332664a25596a9e316441 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Wed, 8 Nov 2023 12:06:18 -0700 Subject: [PATCH 36/36] Move subglacial hydro and bed topography updates back to their former positions Move subglacial hydro and bed topography updates back to their original positions. These were initially moved due to the complexity of interaction with RK integration, but I think subsequent changes have made that move unnecessary. --- .../mpas_li_time_integration_fe_rk.F | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) 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 1ae17ed878db..6ca13f8deb0d 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 @@ -259,17 +259,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("vertical therm") -! === Update subglacial hydrology =========== -! It's not clear where the best place to call this should be. -! Seems sensible to put it after thermal evolution is complete to get updated basal melting source term. -! Also seems (might be?) better to put it after geometry evolution. -! We want it before the velocity solve since the hydro model can control the velo basal b.c. -! This was moved to be before geometry evolution to simplify RK time integration. - call mpas_timer_start("subglacial hydro") - call li_SGH_solve(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("subglacial hydro") - ! *** 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, 'layerThickness', layerThickness) @@ -510,6 +499,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_timer_start("halo updates") call mpas_dmpar_field_halo_exch(domain, 'cellMask') call mpas_dmpar_field_halo_exch(domain, 'edgeMask') + call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") ! Calculate volume converted from grounded to floating @@ -527,6 +517,16 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! Reset time step to full length after RK loop deltat = deltatFull +! === Update subglacial hydrology =========== +! It's not clear where the best place to call this should be. +! Seems sensible to put it after thermal evolution is complete to get updated basal melting source term. +! Also seems (might be?) better to put it after geometry evolution. +! We want it before the velocity solve since the hydro model can control the velo basal b.c. + call mpas_timer_start("subglacial hydro") + call li_SGH_solve(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("subglacial hydro") + ! === Calve ice ======================== call mpas_timer_start("calve_ice") @@ -552,11 +552,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call mpas_dmpar_field_halo_exch(domain, 'vertexMask') call mpas_timer_stop("halo updates") -! === Solve velocity for final state ===================== - if (solveVeloAfterCalving) then - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) - err = ior(err, err_tmp) - endif ! === Update bed topo ===================== ! It's not clear when the best time to do this is. ! Seems cleaner to do it either before or after all of the time evolution of the ice @@ -564,6 +559,12 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) call li_bedtopo_solve(domain, err=err_tmp) err = ior(err, err_tmp) +! === Solve velocity for final state ===================== + if (solveVeloAfterCalving) then + call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + err = ior(err, err_tmp) + endif + ! === Calculate diagnostic variables for new state ===================== call li_calculate_diagnostic_vars(domain, err=err_tmp)