diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index a14205e0d783..0295c2ac9279 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -558,6 +558,38 @@ description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2." possible_values="3 or 4" /> + + + + + + + + domain % blocklist do while (associated(block)) @@ -1405,7 +1416,7 @@ end subroutine iceshelf_melt_ismip6 subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & - applyToFloating, applyToGroundingLine, err) + applyToFloating, applyToGroundingLine, err, computeFaceMeltSpeedNow, applyFrontAblationNow) use li_calving use li_constants, only: oceanFreezingTempDepthDependence @@ -1416,6 +1427,8 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & type (domain_type), intent(inout) :: & domain !< Input/Output: domain object logical, intent(in) :: applyToFloating, applyToGrounded, applyToGroundingLine + logical, intent(in), optional :: computeFaceMeltSpeedNow + logical, intent(in), optional :: applyFrontAblationNow !----------------------------------------------------------------- ! output variables !----------------------------------------------------------------- @@ -1464,11 +1477,16 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & logical, pointer :: config_calculate_thermal_forcing logical, pointer :: config_ocean_data_extrapolation integer, pointer :: config_ocean_data_extrap_ncells_extra + logical :: doComputeFaceMeltSpeed, doApplyFrontAblation real (kind=RKIND) :: dzAccumulated, dz integer :: err_tmp, kk, iLayer err = 0 call mpas_log_write('Starting face melt routine') + doComputeFaceMeltSpeed = .true. + doApplyFrontAblation = .true. + if (present(computeFaceMeltSpeedNow)) doComputeFaceMeltSpeed = computeFaceMeltSpeedNow + if (present(applyFrontAblationNow)) doApplyFrontAblation = applyFrontAblationNow ! Get sea level, bedTopography, ice density call mpas_pool_get_config(liConfigs, 'config_ice_density', rhoi) @@ -1511,7 +1529,7 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & call mpas_pool_get_array(geometryPool, 'faceMeltingCFLdt', faceMeltingCFLdt) call mpas_pool_get_array(geometryPool, 'dtFaceMeltingCFLratio', dtFaceMeltingCFLratio) - if ( config_use_3d_thermal_forcing_for_face_melt ) then + if (doComputeFaceMeltSpeed .and. config_use_3d_thermal_forcing_for_face_melt ) then call mpas_pool_get_dimension(meshPool, 'nISMIP6OceanLayers', nISMIP6OceanLayers) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_zOcean', zOcean) @@ -1624,29 +1642,37 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & TFocean = TFocean + ismip6shelfMelt_deltaT !apply dT correction endif - faceMeltSpeed(:) = 0.0_RKIND - allocate(faceMeltSpeedVertAvg(nCells+1)) - faceMeltSpeedVertAvg(:) = 0.0_RKIND + if (doComputeFaceMeltSpeed) then + faceMeltSpeed(:) = 0.0_RKIND - ! Calculate face melting for each cell - do iCell = 1, nCells + ! Calculate face-melt speed for each cell + do iCell = 1, nCells - if ( bedTopography(iCell) < 0.0_RKIND ) then - waterDepth = seaLevel - bedTopography(iCell) - else - waterDepth = 0.0_RKIND - endif + if ( bedTopography(iCell) < 0.0_RKIND ) then + waterDepth = seaLevel - bedTopography(iCell) + else + waterDepth = 0.0_RKIND + endif - if (trim(config_front_mass_bal_grounded) == 'ismip6') then - ! Calculate ice front melt rate at each cell - faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1 - (ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial & - + B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay - elseif (trim(config_front_mass_bal_grounded) == 'uniform') then - faceMeltSpeed(iCell) = config_uniform_face_melt_rate - endif + if (trim(config_front_mass_bal_grounded) == 'ismip6') then + ! Calculate ice front melt rate at each cell + faceMeltSpeed(iCell) = (aSubglacial * waterDepth * & ! m s^-1 + (ismip6Runoff(iCell) / rhow * secPerDay)**alphaSubglacial & + + B) * max(0.0_RKIND, TFocean(iCell) + addTFocean)**betaTF / secPerDay + elseif (trim(config_front_mass_bal_grounded) == 'uniform') then + faceMeltSpeed(iCell) = config_uniform_face_melt_rate + endif - enddo + enddo + endif + + if (.not. doApplyFrontAblation) then + block => block % next + cycle + endif + + allocate(faceMeltSpeedVertAvg(nCells+1)) + faceMeltSpeedVertAvg(:) = 0.0_RKIND where ( (thickness > 0.0_RKIND) .and. (lowerSurface < 0.0_RKIND) ) faceMeltSpeedVertAvg(:) = faceMeltSpeed(:) * abs(lowerSurface(:)) / thickness(:) @@ -1660,12 +1686,12 @@ subroutine grounded_face_melt_ismip6(domain, applyToGrounded, & maxDt=faceMeltingCFLdt, CFLratio=dtFaceMeltingCFLratio, err=err_tmp) err = ior(err, err_tmp) + deallocate(faceMeltSpeedVertAvg) + block => block % next enddo ! associated(block) - deallocate(faceMeltSpeedVertAvg) - end subroutine grounded_face_melt_ismip6 !----------------------------------------------------------------------- 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 9581e32adc2d..5aca473bb557 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 @@ -110,16 +110,26 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) type (block_type), pointer :: block integer :: err_tmp type (mpas_pool_type), pointer :: geometryPool, thermalPool, meshPool, velocityPool + type (mpas_pool_type), pointer :: meshPoolTmp logical, pointer :: config_restore_calving_front, & 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 + logical, pointer :: config_use_embedded_ssprk character (len=StrKIND), pointer :: config_thickness_advection, & config_tracer_advection character (len=StrKIND), pointer :: config_thermal_solver character (len=StrKIND), pointer :: config_time_integration + real (kind=RKIND), pointer :: config_embedded_rk_rtol, & + config_embedded_rk_atol_thickness, & + config_embedded_rk_atol_tracer, & + config_embedded_rk_safety_factor, & + config_embedded_rk_max_growth, & + config_embedded_rk_min_shrink, & + config_min_adaptive_timestep + integer, pointer :: config_embedded_rk_max_retries integer, pointer :: config_rk_order, config_rk3_stages integer :: rkStage, iCell, iTracer, k real (kind=RKIND), dimension(:,:), pointer :: layerThickness, & @@ -145,32 +155,73 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) fluxAcrossGroundingLineAccum, & fluxAcrossGroundingLineOnCellsAccum, & calvingThicknessAccum - real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions - integer, dimension(:), pointer :: cellMaskPrev ! cell mask before advection + real (kind=RKIND), dimension(:), pointer :: layerCenterSigma, layerThicknessFractions, areaCell + integer, dimension(:), allocatable :: cellMaskPrev, edgeMaskPrev, vertexMaskPrev ! cell mask before advection real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, & layerThicknessTmp, & + layerThicknessEmbedded, & + layerThicknessStage1, & + layerThicknessStage2, & temperaturePrev, & + temperatureEmbedded, & + temperatureStage1, & + temperatureStage2, & waterFracPrev, & + waterFracEmbedded, & + waterFracStage1, & + waterFracStage2, & passiveTracer3d, & passiveTracer3dPrev, & damage3d, & damage3dPrev + real (kind=RKIND), dimension(:), allocatable :: damagePrev, & + passiveTracer2dPrev, & + damageEmbedded, & + passiveTracer2dEmbedded, & + damageStage1, & + damageStage2, & + passiveTracer2dStage1, & + passiveTracer2dStage2 integer, pointer :: nVertLevels - integer, pointer :: nCells, nEdges + integer, pointer :: nCells, nEdges, nVertices integer, pointer :: albanyVelocityError integer :: iCell1, iCell2, iEdge, theGroundedCell - integer, dimension(:), pointer :: edgeMask, cellMask + integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask real (kind=RKIND), pointer :: deltat, config_ice_density + real (kind=RKIND), pointer :: deltatTmp real (kind=RKIND) :: deltatFull + real (kind=RKIND) :: deltatStage 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 logical :: solveVeloAfterCalving + logical :: useEmbeddedSSPRK, includeThermalForError, includeDamageForError, stepAccepted + integer :: stepAttempt, maxStepAttempts, embeddedHighOrder + real (kind=RKIND) :: rkErrorNorm, rkErrorNormGlobal, rkErrorSumSqLocal, rkErrorSumSqGlobal, dtScale + real (kind=RKIND) :: rkErrorWeightLocal, rkErrorWeightGlobal + real (kind=RKIND) :: rk2EmbeddedW1, rk2EmbeddedW2 + real (kind=RKIND), save :: previousEmbeddedAcceptedDt = -1.0_RKIND err = 0 err_tmp = 0 solveVeloAfterCalving = .false. + useEmbeddedSSPRK = .false. + includeThermalForError = .false. + includeDamageForError = .false. + stepAccepted = .false. + stepAttempt = 1 + maxStepAttempts = 1 + embeddedHighOrder = 3 + rkErrorNorm = 0.0_RKIND + rkErrorNormGlobal = 0.0_RKIND + rkErrorSumSqLocal = 0.0_RKIND + rkErrorSumSqGlobal = 0.0_RKIND + rkErrorWeightLocal = 0.0_RKIND + rkErrorWeightGlobal = 0.0_RKIND + rk2EmbeddedW1 = 0.694021459207626_RKIND + rk2EmbeddedW2 = 0.305978540792374_RKIND + dtScale = 1.0_RKIND 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) @@ -184,6 +235,15 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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_config(liConfigs, 'config_use_embedded_ssprk', config_use_embedded_ssprk) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_rtol', config_embedded_rk_rtol) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_atol_thickness', config_embedded_rk_atol_thickness) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_atol_tracer', config_embedded_rk_atol_tracer) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_safety_factor', config_embedded_rk_safety_factor) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_max_growth', config_embedded_rk_max_growth) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_min_shrink', config_embedded_rk_min_shrink) + call mpas_pool_get_config(liConfigs, 'config_embedded_rk_max_retries', config_embedded_rk_max_retries) + call mpas_pool_get_config(liConfigs, 'config_min_adaptive_timestep', config_min_adaptive_timestep) call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(domain % blocklist % structs, 'thermal', thermalPool) @@ -193,9 +253,11 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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_array(meshPool, 'areaCell', areaCell) 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_dimension(meshPool, 'nVertices', nVertices) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) call mpas_pool_get_array(velocityPool, 'albanyVelocityError', albanyVelocityError) @@ -220,6 +282,28 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) allocate(passiveTracer3dPrev(nVertLevels, nCells+1)) allocate(passiveTracer3d(nVertLevels, nCells+1)) allocate(cellMaskPrev(nCells+1)) + allocate(edgeMaskPrev(nEdges+1)) + allocate(vertexMaskPrev(nVertices+1)) + allocate(damagePrev(nCells+1)) + allocate(passiveTracer2dPrev(nCells+1)) + + if (config_use_embedded_ssprk) then + allocate(layerThicknessEmbedded(nVertLevels, nCells+1)) + allocate(layerThicknessStage1(nVertLevels, nCells+1)) + allocate(layerThicknessStage2(nVertLevels, nCells+1)) + allocate(temperatureEmbedded(nVertLevels, nCells+1)) + allocate(waterFracEmbedded(nVertLevels, nCells+1)) + allocate(temperatureStage1(nVertLevels, nCells+1)) + allocate(temperatureStage2(nVertLevels, nCells+1)) + allocate(waterFracStage1(nVertLevels, nCells+1)) + allocate(waterFracStage2(nVertLevels, nCells+1)) + allocate(damageEmbedded(nCells+1)) + allocate(passiveTracer2dEmbedded(nCells+1)) + allocate(damageStage1(nCells+1)) + allocate(damageStage2(nCells+1)) + allocate(passiveTracer2dStage1(nCells+1)) + allocate(passiveTracer2dStage2(nCells+1)) + endif allocate(sfcMassBalAccum(nCells+1)) allocate(basalMassBalAccum(nCells+1)) @@ -238,7 +322,29 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) damage3d(:,:) = 0.0_RKIND passiveTracer3dPrev(:,:) = 0.0_RKIND cellMaskPrev(:) = 0 + edgeMaskPrev(:) = 0 + vertexMaskPrev(:) = 0 passiveTracer3d(:,:) = 0.0_RKIND + damagePrev(:) = 0.0_RKIND + passiveTracer2dPrev(:) = 0.0_RKIND + + if (config_use_embedded_ssprk) then + layerThicknessEmbedded(:,:) = 0.0_RKIND + layerThicknessStage1(:,:) = 0.0_RKIND + layerThicknessStage2(:,:) = 0.0_RKIND + temperatureEmbedded(:,:) = 0.0_RKIND + waterFracEmbedded(:,:) = 0.0_RKIND + temperatureStage1(:,:) = 0.0_RKIND + temperatureStage2(:,:) = 0.0_RKIND + waterFracStage1(:,:) = 0.0_RKIND + waterFracStage2(:,:) = 0.0_RKIND + damageEmbedded(:) = 0.0_RKIND + passiveTracer2dEmbedded(:) = 0.0_RKIND + damageStage1(:) = 0.0_RKIND + damageStage2(:) = 0.0_RKIND + passiveTracer2dStage1(:) = 0.0_RKIND + passiveTracer2dStage2(:) = 0.0_RKIND + endif sfcMassBalAccum(:) = 0.0_RKIND basalMassBalAccum(:) = 0.0_RKIND @@ -256,12 +362,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) err = ior(err, err_tmp) call mpas_timer_stop("advection prep") -! === Advance the clock before all other physics happen =========== - call mpas_timer_start("advancing clock") - call advance_clock(domain, err_tmp) - err = ior(err, err_tmp) - call mpas_timer_stop("advancing clock") - !TODO: Determine whether grounded melting should in fact be called first ! === Ocean forcing extrapolation into ice-shelf cavities =========== call mpas_timer_start("ocean forcing extrapolation") @@ -271,7 +371,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! === Face melting for grounded ice =========== call mpas_timer_start("face melting for grounded ice") - call li_face_melt_grounded_ice(domain, err_tmp) + ! Compute face-melt tendencies/CFL diagnostics now, but defer applying + ! the thickness update until after RK accepts the timestep. + call li_face_melt_grounded_ice(domain, err_tmp, applyFaceMeltNow=.false.) err = ior(err, err_tmp) call mpas_timer_stop("face melting for grounded ice") @@ -293,6 +395,9 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) 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(geometryPool, 'edgeMask', edgeMask) + call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) + call mpas_pool_get_array(thermalPool, 'temperature', temperature) call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) @@ -301,6 +406,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) waterFracPrev(:,:) = waterFrac(:,:) layerThicknessPrev(:,:) = layerThickness(:,:) cellMaskPrev(:) = cellMask(:) + edgeMaskPrev(:) = edgeMask(:) + vertexMaskPrev(:) = vertexMask(:) + damagePrev(:) = damage(:) + passiveTracer2dPrev(:) = passiveTracer2d(:) do k = 1, nVertLevels damage3dPrev(k,:) = damage(:) passiveTracer3dPrev(k,:) = passiveTracer2d(:) @@ -393,151 +502,376 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! to preserve the accuracy of time integration. call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) err = ior(err, err_tmp) - -! *** Start RK loop *** - do rkStage = 1, nRKstages - call mpas_log_write('beginning rk stage $i of $i', & - intArgs=(/rkStage, nRKstages/)) - deltat = deltatFull * rkSubstepWeights(rkStage) - - ! === 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") + useEmbeddedSSPRK = config_use_embedded_ssprk + if (useEmbeddedSSPRK) then + if (.not. ((trim(config_time_integration) == 'runge_kutta') .and. & + ((config_rk_order == 2) .or. (config_rk_order == 3)))) then + err = 1 + call mpas_log_write('Embedded SSPRK currently supports only runge_kutta with config_rk_order=2 or 3', & + messageType=MPAS_LOG_ERR) + return + endif + if ((config_rk_order == 3) .and. (config_rk3_stages /= 3 .and. config_rk3_stages /= 4)) then + err = 1 + call mpas_log_write('Embedded SSPRK currently supports config_rk3_stages = 3 or 4', & + messageType=MPAS_LOG_ERR) + return endif + 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") + embeddedHighOrder = config_rk_order - ! 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. - if ( (rkStage > 1) .and. (rkStage < 4) ) 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) + includeThermalForError = trim(config_thermal_solver) .ne. 'none' + includeDamageForError = .false. + if (useEmbeddedSSPRK .and. config_calculate_damage) then + call mpas_log_write('Embedded SSPRK error control excludes damage field to avoid non-smooth limiter dominance.') + endif - call mpas_pool_get_array(thermalPool, 'temperature', temperature) - call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac) + if (useEmbeddedSSPRK .and. (previousEmbeddedAcceptedDt > 0.0_RKIND)) then + deltatFull = min(deltatFull, previousEmbeddedAcceptedDt * config_embedded_rk_max_growth) + endif + maxStepAttempts = 1 + if (useEmbeddedSSPRK) then + maxStepAttempts = max(1, config_embedded_rk_max_retries + 1) + endif - layerThicknessTmp(:,:) = layerThickness(:,:) - layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & - (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) - thickness = sum(layerThickness, 1) - ! Do not calculate masks after updating thickness! We need to keep masks - ! constant for now to preserve accuracy of time integration + stepAccepted = .false. + do stepAttempt = 1, maxStepAttempts + ! Reset prognostic state to pre-RK values before each attempt. + layerThickness(:,:) = layerThicknessPrev(:,:) + thickness(:) = sum(layerThickness, 1) + temperature(:,:) = temperaturePrev(:,:) + waterFrac(:,:) = waterFracPrev(:,:) + damage(:) = damagePrev(:) + passiveTracer2d(:) = passiveTracer2dPrev(:) + cellMask(:) = cellMaskPrev(:) + edgeMask(:) = edgeMaskPrev(:) + vertexMask(:) = vertexMaskPrev(:) + + sfcMassBalAccum(:) = 0.0_RKIND + basalMassBalAccum(:) = 0.0_RKIND + groundedSfcMassBalAccum(:) = 0.0_RKIND + groundedBasalMassBalAccum(:) = 0.0_RKIND + floatingBasalMassBalAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineAccum(:) = 0.0_RKIND + fluxAcrossGroundingLineOnCellsAccum(:) = 0.0_RKIND + calvingThicknessAccum(:) = 0.0_RKIND + calvingThickness(:) = 0.0_RKIND + calvingThicknessFromThreshold(:) = 0.0_RKIND + calvingVelocity(:) = 0.0_RKIND + + if (useEmbeddedSSPRK) then + layerThicknessEmbedded(:,:) = layerThicknessPrev(:,:) + temperatureEmbedded(:,:) = temperaturePrev(:,:) + waterFracEmbedded(:,:) = waterFracPrev(:,:) + damageEmbedded(:) = damagePrev(:) + passiveTracer2dEmbedded(:) = passiveTracer2dPrev(:) + endif - if (trim(config_thermal_solver) .ne. 'none') then - do iCell = 1, nCells - 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 - endif +! *** Start RK loop *** + do rkStage = 1, nRKstages + call mpas_log_write('beginning rk stage $i of $i', & + intArgs=(/rkStage, nRKstages/)) + deltatStage = deltatFull * rkSubstepWeights(rkStage) + ! Keep deltat consistent on all blocks for retry attempts. + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) + call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) + deltatTmp = deltatStage + block => block % next + enddo + deltat = deltatStage + ! === 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") + + ! 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. + if ( (rkStage > 1) .and. (rkStage < 4) ) 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, 'waterFrac', waterFrac) + + layerThicknessTmp(:,:) = layerThickness(:,:) + layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + & + (1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:) + thickness = sum(layerThickness, 1) + ! Do not calculate masks after updating thickness! We need to keep masks + ! constant for now to preserve accuracy of time integration + + if (trim(config_thermal_solver) .ne. 'none') then + do iCell = 1, nCells + 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 + 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) * damage3dPrev(k,iCell) * layerThicknessPrev(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 - damage3d(k,iCell) = ( rkSSPweights(rkStage) * damage3dPrev(k,iCell) * layerThicknessPrev(k,iCell) + & - (1.0_RKIND - rkSSPweights(rkStage)) * damage(iCell) * & - layerThicknessTmp(k,iCell) ) / layerThickness(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 - damage3d(k,iCell) = 0.0_RKIND + passiveTracer3d(k,iCell) = 0.0_RKIND endif enddo - damage(iCell) = sum(damage3d(:, iCell) * layerThicknessFractions) + passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) enddo + endif + + ! === 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) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) + call mpas_pooL_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) + call mpas_pool_get_array(geometryPool, 'thickness', thickness) - do iCell = 1, nCells - do k = 1, nVertLevels - if (layerThickness(k,iCell) > 0.0_RKIND) then - 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 - passiveTracer3d(k,iCell) = 0.0_RKIND + ! 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 + fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells + + ! 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, 'damage') + call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + + call mpas_timer_stop("halo updates") + + if (useEmbeddedSSPRK) then + if (rkStage == 1) then + layerThicknessStage1(:,:) = layerThickness(:,:) + if (includeThermalForError) then + temperatureStage1(:,:) = temperature(:,:) + waterFracStage1(:,:) = waterFrac(:,:) endif - enddo - passiveTracer2d(iCell) = sum(passiveTracer3d(:,iCell) * layerThicknessFractions) - enddo + damageStage1(:) = damage(:) + passiveTracer2dStage1(:) = passiveTracer2d(:) + elseif (rkStage == 2) then + layerThicknessStage2(:,:) = layerThickness(:,:) + if (includeThermalForError) then + temperatureStage2(:,:) = temperature(:,:) + waterFracStage2(:,:) = waterFrac(:,:) + endif + damageStage2(:) = damage(:) + passiveTracer2dStage2(:) = passiveTracer2d(:) + endif + endif + + ! Update velocity for each RK step + ! === Solve Velocity ===================== + 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) + call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) + calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness + endif + ! We need to remove icebergs between RK stages because the + ! main calving routine is not called until after the RK loop. + ! This frequently results in icebergs that causes intermediate + ! RK stage velocity solves to fail. + call li_remove_icebergs(domain) + call mpas_log_write('Resetting masks') + + cellMask(:) = cellMaskPrev(:) + edgeMask(:) = edgeMaskPrev(:) + vertexMask(:) = vertexMaskPrev(:) + call li_velocity_solve(domain, solveVelo=.true., updateMask=.false., err=err_tmp) + err = ior(err, err_tmp) + endif +! *** end RK loop *** + enddo + if (.not. useEmbeddedSSPRK) then + stepAccepted = .true. + exit endif - - ! === 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") + + if (config_rk_order == 2) then + ! SSPRK2(2,2) embedded companion from Fekete et al. (2022), + ! using convex stage-state weights. + layerThicknessEmbedded(:,:) = rk2EmbeddedW1 * layerThicknessStage1(:,:) + & + rk2EmbeddedW2 * layerThickness(:,:) + if (includeThermalForError) then + temperatureEmbedded(:,:) = rk2EmbeddedW1 * temperatureStage1(:,:) + & + rk2EmbeddedW2 * temperature(:,:) + waterFracEmbedded(:,:) = rk2EmbeddedW1 * waterFracStage1(:,:) + & + rk2EmbeddedW2 * waterFrac(:,:) + endif + damageEmbedded(:) = rk2EmbeddedW1 * damageStage1(:) + & + rk2EmbeddedW2 * damage(:) + passiveTracer2dEmbedded(:) = rk2EmbeddedW1 * passiveTracer2dStage1(:) + & + rk2EmbeddedW2 * passiveTracer2d(:) + elseif (config_rk3_stages == 3) then + ! SSPRK3(3,3): use stage-2 endpoint as lower-order companion. + layerThicknessEmbedded(:,:) = layerThicknessStage2(:,:) + if (includeThermalForError) then + temperatureEmbedded(:,:) = temperatureStage2(:,:) + waterFracEmbedded(:,:) = waterFracStage2(:,:) + endif + damageEmbedded(:) = damageStage2(:) + passiveTracer2dEmbedded(:) = passiveTracer2dStage2(:) + else + ! SSPRK3(4,3): use stage-2 endpoint as SSP-consistent lower-order companion. + layerThicknessEmbedded(:,:) = layerThicknessStage2(:,:) + if (includeThermalForError) then + temperatureEmbedded(:,:) = temperatureStage2(:,:) + waterFracEmbedded(:,:) = waterFracStage2(:,:) + endif + damageEmbedded(:) = damageStage2(:) + passiveTracer2dEmbedded(:) = passiveTracer2dStage2(:) 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) - 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 - 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 - fluxAcrossGroundingLineOnCellsAccum = fluxAcrossGroundingLineOnCellsAccum + rkTendWeights(rkStage) * fluxAcrossGroundingLineOnCells - - ! Halo updates - call mpas_timer_start("halo updates") + call li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & + waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & + damage, damageEmbedded, areaCell, nVertLevels, nCells, includeThermalForError, & + includeDamageForError, config_embedded_rk_rtol, & + config_embedded_rk_atol_thickness, config_embedded_rk_atol_tracer, & + rkErrorSumSqLocal, rkErrorWeightLocal) + + call mpas_dmpar_sum_real(domain % dminfo, rkErrorSumSqLocal, rkErrorSumSqGlobal) + call mpas_dmpar_sum_real(domain % dminfo, rkErrorWeightLocal, rkErrorWeightGlobal) + if (rkErrorWeightGlobal > 0.0_RKIND) then + rkErrorNormGlobal = sqrt(rkErrorSumSqGlobal / rkErrorWeightGlobal) + else + rkErrorNormGlobal = 0.0_RKIND + endif + rkErrorNorm = rkErrorNormGlobal - 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, 'damage') - call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d') + if (rkErrorNorm <= 1.0_RKIND) then + stepAccepted = .true. + call mpas_log_write('Embedded SSPRK accepted attempt $i with normalized error $r', & + intArgs=(/stepAttempt/), realArgs=(/rkErrorNorm/)) + exit + endif - call mpas_timer_stop("halo updates") + if (stepAttempt == maxStepAttempts) then + err = 1 + call mpas_log_write('Embedded SSPRK failed to satisfy tolerance after $i attempts (error=$r)', & + intArgs=(/maxStepAttempts/), realArgs=(/rkErrorNorm/), & + messageType=MPAS_LOG_ERR) + return + endif - ! Update velocity for each RK step - ! === Solve Velocity ===================== - if ( config_update_velocity_before_calving .or. ( (.not. config_update_velocity_before_calving) & - .and. (rkStage < nRKstages) ) ) then + dtScale = config_embedded_rk_safety_factor * rkErrorNorm ** (-1.0_RKIND / real(embeddedHighOrder, RKIND)) + dtScale = min(config_embedded_rk_max_growth, dtScale) + dtScale = max(config_embedded_rk_min_shrink, dtScale) + deltatFull = deltatFull * dtScale - 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) - call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - calvingThicknessAccum = calvingThicknessAccum + rkTendWeights(rkStage) * calvingThickness - endif - ! We need to remove icebergs between RK stages because the - ! main calving routine is not called until after the RK loop. - ! This frequently results in icebergs that causes intermediate - ! RK stage velocity solves to fail. - call li_remove_icebergs(domain) + if (deltatFull < config_min_adaptive_timestep) then + err = 1 + call mpas_log_write('Embedded SSPRK requested deltat below config_min_adaptive_timestep', & + messageType=MPAS_LOG_ERR) + return + endif - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) - err = ior(err, err_tmp) + call mpas_log_write('Embedded SSPRK rejected attempt $i; reducing deltat by factor $r to $r seconds', & + intArgs=(/stepAttempt/), realArgs=(/dtScale, deltatFull/)) + enddo + + if (.not. stepAccepted) then + err = 1 + call mpas_log_write('Embedded SSPRK did not accept a step', messageType=MPAS_LOG_ERR) + return + endif + + if (useEmbeddedSSPRK) then + deltatFull = real(floor(deltatFull, KIND=8), RKIND) + if (deltatFull < config_min_adaptive_timestep) then + err = 1 + call mpas_log_write('Embedded SSPRK accepted deltat falls below config_min_adaptive_timestep after flooring', & + messageType=MPAS_LOG_ERR) + return endif -! *** end RK loop *** + endif + + if (useEmbeddedSSPRK) previousEmbeddedAcceptedDt = deltatFull + + ! Restore full accepted timestep on all blocks before advancing the clock. + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPoolTmp) + call mpas_pool_get_array(meshPoolTmp, 'deltat', deltatTmp) + deltatTmp = deltatFull + block => block % next enddo + deltat = deltatFull + + ! Apply grounded face melt once, using the accepted timestep. + call mpas_timer_start("face melting for grounded ice") + call li_face_melt_grounded_ice(domain, err_tmp, computeFaceMeltNow=.false.) + err = ior(err, err_tmp) + call mpas_timer_stop("face melting for grounded ice") + + ! Advance model time only once, after this step has been accepted. + call mpas_timer_start("advancing clock") + call advance_clock(domain, err_tmp) + err = ior(err, err_tmp) + call mpas_timer_stop("advancing clock") ! Finalize budget updates ! Update masks after RK integration @@ -576,9 +910,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) advMask2ndOrder) endif -! 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. @@ -623,7 +954,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) ! === Solve velocity for final state ===================== if (solveVeloAfterCalving) then - call li_velocity_solve(domain, solveVelo=.true., err=err_tmp) + call li_velocity_solve(domain, solveVelo=.true., updateMask=.true., err=err_tmp) err = ior(err, err_tmp) endif @@ -649,12 +980,32 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err) deallocate(passiveTracer3dPrev) deallocate(passiveTracer3d) deallocate(cellMaskPrev) + deallocate(edgeMaskPrev) + deallocate(vertexMaskPrev) + deallocate(damagePrev) + deallocate(passiveTracer2dPrev) + if (allocated(layerThicknessEmbedded)) deallocate(layerThicknessEmbedded) + if (allocated(layerThicknessStage1)) deallocate(layerThicknessStage1) + if (allocated(layerThicknessStage2)) deallocate(layerThicknessStage2) + if (allocated(temperatureEmbedded)) deallocate(temperatureEmbedded) + if (allocated(waterFracEmbedded)) deallocate(waterFracEmbedded) + if (allocated(temperatureStage1)) deallocate(temperatureStage1) + if (allocated(temperatureStage2)) deallocate(temperatureStage2) + if (allocated(waterFracStage1)) deallocate(waterFracStage1) + if (allocated(waterFracStage2)) deallocate(waterFracStage2) + if (allocated(damageEmbedded)) deallocate(damageEmbedded) + if (allocated(passiveTracer2dEmbedded)) deallocate(passiveTracer2dEmbedded) + if (allocated(damageStage1)) deallocate(damageStage1) + if (allocated(damageStage2)) deallocate(damageStage2) + if (allocated(passiveTracer2dStage1)) deallocate(passiveTracer2dStage1) + if (allocated(passiveTracer2dStage2)) deallocate(passiveTracer2dStage2) deallocate(sfcMassBalAccum) deallocate(basalMassBalAccum) deallocate(groundedSfcMassBalAccum) deallocate(groundedBasalMassBalAccum) deallocate(floatingBasalMassBalAccum) deallocate(fluxAcrossGroundingLineAccum) + deallocate(fluxAcrossGroundingLineOnCellsAccum) deallocate(calvingThicknessAccum) !-------------------------------------------------------------------- end subroutine li_time_integrator_forwardeuler_rungekutta @@ -1281,7 +1632,6 @@ subroutine advection_solver(domain, err) call li_calculate_layerThickness(meshPool, thickness, layerThickness) endif - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) call li_update_geometry(geometryPool) @@ -1557,6 +1907,92 @@ subroutine calculate_layerThicknessEdge(meshPool, geometryPool, velocityPool, er !-------------------------------------------------------------------- end subroutine calculate_layerThicknessEdge +!*********************************************************************** +! +! routine li_embedded_ssprk_error_norm +! +!> \brief Compute normalized error estimate for embedded SSPRK control +!> \details +!> Uses the current RK endpoint and a formally reconstructed embedded +!> lower-order endpoint to compute an area-weighted normalized error over +!> thickness and tracer-like fields. +!----------------------------------------------------------------------- + subroutine li_embedded_ssprk_error_norm(layerThickness, layerThicknessEmbedded, temperature, temperatureEmbedded, & + waterFrac, waterFracEmbedded, passiveTracer2d, passiveTracer2dEmbedded, & + damage, damageEmbedded, areaCell, nVertLevels, nCells, includeThermal, & + includeDamage, relTol, absTolThickness, absTolTracer, errorSumSq, errorWeight) + + real (kind=RKIND), dimension(:,:), intent(in) :: layerThickness, layerThicknessEmbedded + real (kind=RKIND), dimension(:,:), intent(in) :: temperature, temperatureEmbedded + real (kind=RKIND), dimension(:,:), intent(in) :: waterFrac, waterFracEmbedded + real (kind=RKIND), dimension(:), intent(in) :: passiveTracer2d, passiveTracer2dEmbedded + real (kind=RKIND), dimension(:), intent(in) :: damage, damageEmbedded + real (kind=RKIND), dimension(:), intent(in) :: areaCell + integer, intent(in) :: nVertLevels, nCells + logical, intent(in) :: includeThermal, includeDamage + real (kind=RKIND), intent(in) :: relTol, absTolThickness, absTolTracer + real (kind=RKIND), intent(out) :: errorSumSq + real (kind=RKIND), intent(out) :: errorWeight + + integer :: iCell, k + real (kind=RKIND) :: yPrimary, yEmbedded, scaleVal, thisErr + real (kind=RKIND), parameter :: minIceForError = 1.0e-12_RKIND + + errorSumSq = 0.0_RKIND + errorWeight = 0.0_RKIND + + do iCell = 1, nCells + yPrimary = sum(layerThickness(:, iCell)) + yEmbedded = sum(layerThicknessEmbedded(:, iCell)) + scaleVal = absTolThickness + relTol * max(abs(yPrimary), abs(yEmbedded)) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(yPrimary - yEmbedded) / scaleVal + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) + endif + + if (max(yPrimary, yEmbedded) <= minIceForError) then + cycle + endif + + scaleVal = absTolTracer + relTol * max(abs(passiveTracer2d(iCell)), abs(passiveTracer2dEmbedded(iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(passiveTracer2d(iCell) - passiveTracer2dEmbedded(iCell)) / scaleVal + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) + endif + + if (includeDamage) then + scaleVal = absTolTracer + relTol * max(abs(damage(iCell)), abs(damageEmbedded(iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(damage(iCell) - damageEmbedded(iCell)) / scaleVal + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) + endif + endif + + if (includeThermal) then + do k = 1, nVertLevels + if (max(layerThickness(k,iCell), layerThicknessEmbedded(k,iCell)) <= minIceForError) cycle + scaleVal = absTolTracer + relTol * max(abs(temperature(k,iCell)), abs(temperatureEmbedded(k,iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(temperature(k,iCell) - temperatureEmbedded(k,iCell)) / scaleVal + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) + endif + + scaleVal = absTolTracer + relTol * max(abs(waterFrac(k,iCell)), abs(waterFracEmbedded(k,iCell))) + if (scaleVal > 0.0_RKIND) then + thisErr = abs(waterFrac(k,iCell) - waterFracEmbedded(k,iCell)) / scaleVal + errorSumSq = errorSumSq + areaCell(iCell) * thisErr * thisErr + errorWeight = errorWeight + areaCell(iCell) + endif + enddo + endif + enddo + + end subroutine li_embedded_ssprk_error_norm + !*********************************************************************** ! ! routine advance_clock @@ -1599,8 +2035,6 @@ subroutine advance_clock(domain, err) type (MPAS_TimeInterval_type) :: timeStepInterval !< the current time step as an interval type (MPAS_Time_type) :: simulationStartTime_timeType - logical, pointer :: config_adaptive_timestep - character (len=StrKIND), pointer :: xtime character (len=StrKIND), pointer :: simulationStartTime character(len=StrKIND) :: timeStamp !< current time as a string @@ -1612,24 +2046,19 @@ subroutine advance_clock(domain, err) err = 0 - call mpas_pool_get_config(liConfigs, 'config_adaptive_timestep', config_adaptive_timestep) - ! === - ! === Adaptive timestep: update clock information + ! === Update clock timestep from mesh deltat ! === - ! Set time step in clock object since the time step could have changed - ! Need to get value out of a block, but all blocks should have the same value, so just use first block - if (config_adaptive_timestep) then - block => domain % blocklist - call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) - call mpas_pool_get_array(meshPool, 'deltat', deltat) - ! convert dt in seconds to timeInterval type - call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) - err = ior(err,err_tmp) - ! update the clock with the timeInterval - call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) - err = ior(err,err_tmp) - endif + ! Need to get value out of a block, but all blocks should have the same value, so just use first block. + block => domain % blocklist + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_array(meshPool, 'deltat', deltat) + ! convert dt in seconds to timeInterval type + call mpas_set_timeInterval(timeStepInterval, dt=deltat, ierr=err_tmp) + err = ior(err,err_tmp) + ! update the clock with the timeInterval + call mpas_set_clock_timestep(domain % clock, timeStepInterval, err_tmp) + err = ior(err,err_tmp) ! === ! === Update clock information diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F index 4d07dbbfdb84..6b01e5124ba4 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity.F @@ -222,7 +222,7 @@ end subroutine li_velocity_block_init !> This routine calls velocity solvers. ! !----------------------------------------------------------------------- - subroutine li_velocity_solve(domain, solveVelo, err) + subroutine li_velocity_solve(domain, solveVelo, updateMask, err) use mpas_vector_reconstruction use li_mask @@ -236,6 +236,8 @@ subroutine li_velocity_solve(domain, solveVelo, err) !< velocity should be solved, or if an existing solution should be !< used to calculate the diagnostic fields related to velocty !< (e.g. on a restart) + logical, intent(in) :: updateMask !< Input: update mask before + !< velocity solve only if outside of RK loop. !----------------------------------------------------------------- ! input/output variables @@ -326,14 +328,15 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_dmpar_field_halo_exch(domain, 'thickness') call mpas_timer_stop("halo updates") - ! Update mask just to be safe (might be redundant) block => domain % blocklist do while (associated(block)) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - - call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + ! Only update mask outside of RK loop + if ( updateMask ) then + call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp) + endif call li_update_geometry(geometryPool) block => block % next diff --git a/components/mpas-albany-landice/src/shared/mpas_li_mask.F b/components/mpas-albany-landice/src/shared/mpas_li_mask.F index 64050cb83adf..6f97e5ed874f 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_mask.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_mask.F @@ -354,7 +354,7 @@ subroutine li_calculate_mask(meshPool, velocityPool, geometryPool, err) ! ==== ! Calculate cellMask values=========================== ! ==== - + call mpas_log_write("Updating masks") !call mpas_timer_start('calculate mask cell') ! Set mask to 0 everywhere, but need to preserve bits the initial ice extent bit do i=1, nCells