From 9149dd1a13c67a57b5c7c41ff08f4f5c50e5b4e1 Mon Sep 17 00:00:00 2001 From: Mauro Perego Date: Fri, 8 Nov 2019 14:01:53 -0700 Subject: [PATCH] Modified logic to determine when FEM needs to be recomputed The mesh needs to be recomputed non only of the vertices ice mask has changed, but also if the floating edges mask and the cell ice mask has changed. Note: in order to compare the clee mask with the previous time steps, we would need the cell mask to have 2 time levels. Howver, if we set the cell mask to have 2 time levels we get different results. This needs to be fixed. --- src/core_landice/Registry.xml | 11 +++- .../analysis_members/mpas_li_global_stats.F | 2 +- .../analysis_members/mpas_li_regional_stats.F | 2 +- .../mode_forward/mpas_li_advection.F | 4 +- .../mode_forward/mpas_li_calving.F | 20 +++---- .../mode_forward/mpas_li_diagnostic_vars.F | 3 +- .../mode_forward/mpas_li_iceshelf_melt.F | 4 +- src/core_landice/mode_forward/mpas_li_sia.F | 2 +- .../mode_forward/mpas_li_statistics.F | 2 +- .../mode_forward/mpas_li_subglacial_hydro.F | 16 +++--- .../mode_forward/mpas_li_thermal.F | 2 +- .../mode_forward/mpas_li_velocity.F | 57 ++++++++++++++++--- .../mode_forward/mpas_li_velocity_external.F | 14 +++-- 13 files changed, 95 insertions(+), 44 deletions(-) diff --git a/src/core_landice/Registry.xml b/src/core_landice/Registry.xml index ec1c2a015a..625a3e5ce5 100644 --- a/src/core_landice/Registry.xml +++ b/src/core_landice/Registry.xml @@ -941,7 +941,7 @@ is the value of that variable from the *previous* time level! - + + diff --git a/src/core_landice/analysis_members/mpas_li_global_stats.F b/src/core_landice/analysis_members/mpas_li_global_stats.F index 41ddc9aebd..b367875979 100644 --- a/src/core_landice/analysis_members/mpas_li_global_stats.F +++ b/src/core_landice/analysis_members/mpas_li_global_stats.F @@ -283,7 +283,7 @@ subroutine li_compute_global_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'basalMassBal', basalMassBal) diff --git a/src/core_landice/analysis_members/mpas_li_regional_stats.F b/src/core_landice/analysis_members/mpas_li_regional_stats.F index 0e0543df57..4c6df7def9 100644 --- a/src/core_landice/analysis_members/mpas_li_regional_stats.F +++ b/src/core_landice/analysis_members/mpas_li_regional_stats.F @@ -275,7 +275,7 @@ subroutine li_compute_regional_stats(domain, memberName, timeLevel, err) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'sfcMassBalApplied', sfcMassBalApplied) call mpas_pool_get_array(geometryPool, 'groundedSfcMassBalApplied', groundedSfcMassBalApplied) diff --git a/src/core_landice/mode_forward/mpas_li_advection.F b/src/core_landice/mode_forward/mpas_li_advection.F index 7dad2e6cb3..fbd80d0afd 100644 --- a/src/core_landice/mode_forward/mpas_li_advection.F +++ b/src/core_landice/mode_forward/mpas_li_advection.F @@ -257,7 +257,7 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'layerThicknessEdge', layerThicknessEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'dynamicThickening', dynamicThickening) call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness) @@ -630,7 +630,7 @@ subroutine li_update_geometry(geometryPool) call mpas_pool_get_dimension(geometryPool, 'nCells', nCells) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) diff --git a/src/core_landice/mode_forward/mpas_li_calving.F b/src/core_landice/mode_forward/mpas_li_calving.F index 91146b7c05..6dff0a04dc 100644 --- a/src/core_landice/mode_forward/mpas_li_calving.F +++ b/src/core_landice/mode_forward/mpas_li_calving.F @@ -402,7 +402,7 @@ subroutine li_restore_calving_front(domain, err) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'restoreThickness', restoreThickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) ! get required fields from the thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -870,7 +870,7 @@ subroutine floating_calving(domain, calvingFraction, err) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) calvingThickness = 0.0_RKIND @@ -929,7 +929,7 @@ subroutine remove_small_islands(meshPool, geometryPool) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) do iCell = 1, nCellsSolve @@ -1039,7 +1039,7 @@ subroutine topographic_calving(domain, calvingFraction, err) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) calvingThickness = 0.0_RKIND @@ -1152,7 +1152,7 @@ subroutine eigencalving(domain, err) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'dvEdge', dvEdge) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'eigencalvingParameter', eigencalvingParameter) @@ -1336,7 +1336,7 @@ subroutine distribute_calving_flux(meshPool, geometryPool, scratchPool, calvingF ! get fields call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'deltat', deltat) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'requiredCalvingVolumeRate', requiredCalvingVolumeRate) @@ -1566,7 +1566,7 @@ subroutine calculate_calving_front_mask(meshPool, geometryPool, calvingFrontMask call mpas_pool_get_dimension(meshPool, 'nCells', nCells) 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, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) calvingFrontMask = 0 !initialize @@ -1656,7 +1656,7 @@ subroutine remove_icebergs(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(geometryPool, 'nCellsSolve', nCellsSolve) @@ -1710,7 +1710,7 @@ subroutine remove_icebergs(domain) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block % structs, 'velocity', velocityPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIceOld', contiguousDynamicIceOld) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -1765,7 +1765,7 @@ subroutine remove_icebergs(domain) block => domain % blocklist do while (associated(block)) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'contiguousDynamicIce', contiguousDynamicIce) diff --git a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F index ab47f0d226..a8f2664137 100644 --- a/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F +++ b/src/core_landice/mode_forward/mpas_li_diagnostic_vars.F @@ -258,10 +258,9 @@ subroutine li_calculate_apparent_diffusivity(meshPool, velocityPool, scratchPool call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel=1) call mpas_pool_get_array(geometryPool, 'apparentDiffusivity', apparentDiffusivity) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'uReconstructZonal', uReconstructAxis1) call mpas_pool_get_array(velocityPool, 'uReconstructMeridional', uReconstructAxis2) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_field(scratchPool, 'workCell', slopeReconstructXField) call mpas_allocate_scratch_field(slopeReconstructXField, .true.) call mpas_pool_get_field(scratchPool, 'workCell2', slopeReconstructYField) diff --git a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F index 95f6716fa8..0b99e9d843 100644 --- a/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F +++ b/src/core_landice/mode_forward/mpas_li_iceshelf_melt.F @@ -247,7 +247,7 @@ subroutine li_basal_melt_floating_ice(domain, err) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) ! get fields from the scratch pool @@ -853,7 +853,7 @@ subroutine calc_iceshelf_draft_info(domain, GLdepth, CFdepth, err) call mpas_pool_get_subpool(block % structs, 'geometry', geometryPool) call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnCell', cellsOnCell) call mpas_pool_get_array(meshPool, 'areaCell', areaCell) diff --git a/src/core_landice/mode_forward/mpas_li_sia.F b/src/core_landice/mode_forward/mpas_li_sia.F index b22001c4f2..d400c12e35 100644 --- a/src/core_landice/mode_forward/mpas_li_sia.F +++ b/src/core_landice/mode_forward/mpas_li_sia.F @@ -316,7 +316,7 @@ subroutine li_sia_solve(meshPool, geometryPool, thermalPool, velocityPool, err) call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'slopeEdge', slopeEdge) call mpas_pool_get_array(geometryPool, 'normalSlopeEdge', normalSlopeEdge) call mpas_pool_get_array(thermalPool, 'temperature', temperature) diff --git a/src/core_landice/mode_forward/mpas_li_statistics.F b/src/core_landice/mode_forward/mpas_li_statistics.F index 47a7be44e8..99de7c4939 100644 --- a/src/core_landice/mode_forward/mpas_li_statistics.F +++ b/src/core_landice/mode_forward/mpas_li_statistics.F @@ -235,7 +235,7 @@ subroutine li_compute_statistics(domain, itimestep) ! Geometry arrays call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness, timeLevel = 1) diff --git a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F index 9e2e4f79ca..596119d01c 100644 --- a/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F +++ b/src/core_landice/mode_forward/mpas_li_subglacial_hydro.F @@ -163,7 +163,7 @@ subroutine li_SGH_init(domain, err) waterPressure = max(0.0_RKIND, waterPressure) waterPressure = min(waterPressure, rhoi * gravity * thickness) ! set pressure correctly under floating ice and open ocean - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) where ( (li_mask_is_floating_ice(cellMask)) .or. & ( (.not. li_mask_is_ice(cellMask)) .and. (bedTopography < config_sea_level) ) ) @@ -379,7 +379,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'deltatSGH', deltatSGH) call mpas_pool_get_array(hydroPool, 'basalMeltInput', basalMeltInput) call mpas_pool_get_array(hydroPool, 'externalWaterInput', externalWaterInput) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) WtillOld = Wtill Wtill = Wtill + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - Cd) @@ -420,7 +420,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'waterVelocity', waterVelocity) call mpas_pool_get_array(hydroPool, 'waterVelocityCellX', waterVelocityCellX) call mpas_pool_get_array(hydroPool, 'waterVelocityCellY', waterVelocityCellY) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) allocate(cellJunk(nCells+1)) call mpas_reconstruct(meshPool, waterVelocity, waterVelocityCellX, waterVelocityCellY, & cellJunk, cellJunk, cellJunk) @@ -476,7 +476,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) divergence(:) = 0.0_RKIND ! zero div before starting ! loop over locally owned cells @@ -551,7 +551,7 @@ subroutine li_SGH_solve(domain, err) call mpas_pool_get_array(hydroPool, 'divergence', divergence) call mpas_pool_get_array(hydroPool, 'divergenceChannel', divergenceChannel) call mpas_pool_get_array(hydroPool, 'channelAreaChangeCell', channelAreaChangeCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) waterThicknessOld = waterThickness waterThickness = waterThicknessOld + deltatSGH * ( (basalMeltInput + externalWaterInput) / rho_water - divergence & @@ -755,7 +755,7 @@ subroutine calc_edge_quantities(block, err) call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeNormal', hydropotentialBaseSlopeNormal) call mpas_pool_get_array(hydroPool, 'waterPressureSlopeNormal', waterPressureSlopeNormal) call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'edgeSignOnCell', edgeSignOnCell) call mpas_pool_get_array(hydroPool, 'hydropotentialBaseSlopeTangent', hydropotentialBaseSlopeTangent) @@ -1217,7 +1217,7 @@ subroutine calc_pressure(block, err) call mpas_pool_get_array(velocityPool, 'basalSpeed', basalSpeed) call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA) call mpas_pool_get_array(geometryPool, 'thickness', thickness) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) openingRate(:) = bedRough * basalSpeed(:) * (bedRoughMax - waterThickness(:)) @@ -1338,7 +1338,7 @@ subroutine calc_pressure_diag_vars(block, err) call mpas_pool_get_array(hydroPool, 'hydropotentialBase', hydropotentialBase) call mpas_pool_get_array(hydroPool, 'waterThickness', waterThickness) call mpas_pool_get_array(hydroPool, 'hydropotential', hydropotential) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) effectivePressure = rhoi * gravity * thickness - waterPressure ! < this should evalute to 0 for floating ice if Pw set correctly there. diff --git a/src/core_landice/mode_forward/mpas_li_thermal.F b/src/core_landice/mode_forward/mpas_li_thermal.F index 5bceb3c9d2..296d5070c6 100644 --- a/src/core_landice/mode_forward/mpas_li_thermal.F +++ b/src/core_landice/mode_forward/mpas_li_thermal.F @@ -820,7 +820,7 @@ subroutine li_thermal_solver(domain, err) ! get fields from the geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'groundedBasalMassBal', groundedBasalMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'basalWaterThickness', basalWaterThickness) diff --git a/src/core_landice/mode_forward/mpas_li_velocity.F b/src/core_landice/mode_forward/mpas_li_velocity.F index e1e543e929..3bd6541abe 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity.F +++ b/src/core_landice/mode_forward/mpas_li_velocity.F @@ -269,8 +269,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) integer, pointer :: nEdges integer, pointer :: nVertInterfaces integer, pointer :: nCells - integer, dimension(:), pointer :: edgeMask, cellMask, vertexMask, vertexMaskOld + integer, dimension(:), pointer :: edgeMask, cellMask, cellMaskOld, vertexMask, vertexMaskOld integer, dimension(:,:), pointer :: dirichletVelocityMaskOld, dirichletVelocityMaskNew + integer, dimension(:), pointer :: floatingEdgesOld, floatingEdgesNew real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, normalVelocityInitial real (kind=RKIND), dimension(:,:), pointer :: uReconstructX, uReconstructY, uReconstructZ, & uReconstructZonal, uReconstructMeridional @@ -288,7 +289,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) real (kind=RKIND), dimension(:), pointer :: upperSurface integer, dimension(:), pointer :: indexToEdgeID integer, pointer :: anyDynamicVertexMaskChanged + integer, pointer :: cellIceMaskChanged integer, pointer :: dirichletMaskChanged + integer, pointer :: floatingEdgesChanged ! truly local variables integer :: cell1, cell2 integer :: cell3, cell4, thisCell @@ -301,7 +304,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) real (kind=RKIND) :: maxThicknessOnProc, maxThicknessAllProcs real (kind=RKIND) :: xVelEdge, yVelEdge integer :: blockDynamicVertexMaskChanged, procDynamicVertexMaskChanged + integer :: blockCellIceMaskChanged, procCellIceMaskChanged integer :: blockDirichletMaskChanged, procDirichletMaskChanged + integer :: blockFloatingEdgesChanged, procFloatingEdgesChanged call mpas_timer_start("velocity solve") @@ -386,7 +391,9 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Don't bother checking this with SIA because it requires an extra global reduce maxThicknessOnProc = 0.0_RKIND ! initialize to procDynamicVertexMaskChanged = 0 + procCellIceMaskChanged = 0 procDirichletMaskChanged = 0 + procFloatingEdgesChanged = 0 block => domain % blocklist do while (associated(block)) @@ -398,7 +405,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) maxThicknessOnProc = max(maxThicknessOnProc, maxval(thickness)) ! The interface expects an array where 1's are floating edges and 0's are non-floating edges. - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) floatingEdges = li_mask_is_floating_ice_int(edgeMask) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel=1) @@ -408,7 +415,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! TODO: There may be some aspects of the mask that are ok change for external dycores, ! but for now just check the whole thing. call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMaskOld, timeLevel=2) - if ( sum(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld)) /= 0 ) then + if ( sum(abs(li_mask_is_dynamic_ice_int(vertexMask) - li_mask_is_dynamic_ice_int(vertexMaskOld))) /= 0 ) then blockDynamicVertexMaskChanged = 1 else blockDynamicVertexMaskChanged = 0 @@ -417,11 +424,26 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Determine if any blocks on this processor had a change to the vertex mask procDynamicVertexMaskChanged = max(procDynamicVertexMaskChanged, blockDynamicVertexMaskChanged) !print *,'procVertexMaskChanged', procVertexMaskChanged + + ! Determine if the cell mask changed during this time step for this block (needed for external dycores) + ! TODO: There may be some aspects of the mask that are ok change for external dycores, + ! but for now just check the whole thing. + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel=1) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMaskOld, timeLevel=2) + if (( sum(abs(li_mask_is_dynamic_ice_int(cellMask) - li_mask_is_dynamic_ice_int(cellMaskOld))) /= 0 ) .or. & + ( sum(abs(li_mask_is_ice_int(cellMask) - li_mask_is_ice_int(cellMaskOld))) /= 0 )) then + blockCellIceMaskChanged = 1 + else + blockCellIceMaskChanged = 0 + endif + + ! Determine if any blocks on this processor had a change to the cell mask + procCellIceMaskChanged = max(procCellIceMaskChanged, blockCellIceMaskChanged) - ! Also check to see if the Dirichlet b.c. mask has changed + ! Check to see if the Dirichlet b.c. mask has changed call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskOld, timeLevel=2) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMaskNew, timeLevel=1) - if ( sum(dirichletVelocityMaskNew - dirichletVelocityMaskOld) /= 0 ) then + if ( sum(abs(dirichletVelocityMaskNew - dirichletVelocityMaskOld)) /= 0 ) then blockDirichletMaskChanged = 1 else blockDirichletMaskChanged = 0 @@ -429,6 +451,17 @@ subroutine li_velocity_solve(domain, solveVelo, err) ! Determine if any blocks on this processor had a change to the vertex mask procDirichletMaskChanged = max(procDirichletMaskChanged, blockDirichletMaskChanged) + ! Check to see if the Floating Edges Ids have changed + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesOld, timeLevel=2) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdgesNew, timeLevel=1) + if ( sum(abs(floatingEdgesNew - floatingEdgesOld)) /= 0 ) then + blockFloatingEdgesChanged = 1 + else + blockFloatingEdgesChanged = 0 + endif + ! Determine if any blocks on this processor had a change to the vertex mask + procFloatingEdgesChanged = max(procFloatingEdgesChanged, blockFloatingEdgesChanged) + ! Set beta for solver to use ! (this could potentially be applied to an SIA solver, so it is calculated in this module) call calculate_beta(block, err_tmp) @@ -452,10 +485,16 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) call mpas_dmpar_max_int(domain % dminfo, procDynamicVertexMaskChanged, anyDynamicVertexMaskChanged) !print *,'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged + ! Do the same for the cell mask + call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged) + call mpas_dmpar_max_int(domain % dminfo, procCellIceMaskChanged, cellIceMaskChanged) ! Do the same for the Dirichlet b.c. mask call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_dmpar_max_int(domain % dminfo, procDirichletMaskChanged, dirichletMaskChanged) - !print *,'dirichletMaskChanged', dirichletMaskChanged + ! Do the same for the Floating edges + call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged) + call mpas_dmpar_max_int(domain % dminfo, procFloatingEdgesChanged, floatingEdgesChanged) + !print *,'floatingEdgesChanged', floatingEdgesChanged endif @@ -479,7 +518,7 @@ subroutine li_velocity_solve(domain, solveVelo, err) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) call mpas_pool_get_array(meshPool, 'cellsOnVertex', cellsOnVertex) call mpas_pool_get_array(meshPool, 'verticesOnEdge', verticesOnEdge) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) @@ -844,7 +883,7 @@ subroutine calculate_beta(block, err) call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve) call mpas_pool_get_array(velocityPool, 'beta', beta) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) betaSolve = beta @@ -1028,7 +1067,7 @@ subroutine remove_fast_ice(meshPool, uReconstructX, uReconstructY, normalVelocit call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) diff --git a/src/core_landice/mode_forward/mpas_li_velocity_external.F b/src/core_landice/mode_forward/mpas_li_velocity_external.F index 53a8e32e17..f1fa8bd529 100644 --- a/src/core_landice/mode_forward/mpas_li_velocity_external.F +++ b/src/core_landice/mode_forward/mpas_li_velocity_external.F @@ -452,7 +452,9 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro logical, pointer :: config_output_external_velocity_solver_data real (kind=RKIND), pointer :: config_ice_density integer, pointer :: anyDynamicVertexMaskChanged + integer, pointer :: cellIceMaskChanged integer, pointer :: dirichletMaskChanged + integer, pointer :: floatingEdgesChanged integer, pointer :: nEdges integer, pointer :: nVertLevels integer, pointer :: timestepNumber @@ -486,7 +488,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) @@ -501,9 +503,11 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(velocityPool, 'uReconstructZ', uReconstructZ) call mpas_pool_get_array(velocityPool, 'betaSolve', betaSolve) call mpas_pool_get_array(velocityPool, 'anyDynamicVertexMaskChanged', anyDynamicVertexMaskChanged) + call mpas_pool_get_array(velocityPool, 'cellIceMaskChanged', cellIceMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletMaskChanged', dirichletMaskChanged) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) + call mpas_pool_get_array(velocityPool, 'floatingEdgesChanged', floatingEdgesChanged) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) ! Hydro variables @@ -521,7 +525,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro ! Note these functions will always be called on the first solve because we ! initialize vertexMask to garbage which sets anyDynamicVertexMaskChanged to 1. if ((anyDynamicVertexMaskChanged == 1) .or. (config_always_compute_fem_grid) .or. & - (dirichletMaskChanged == 1) ) then + (dirichletMaskChanged == 1) .or. (floatingEdgesChanged == 1) .or. (cellIceMaskChanged == 1)) then call mpas_log_write("Generating new external velocity solver FEM grid.", flushNow=.true.) call generate_fem_grid(config_velocity_solver, vertexMask, cellMask, dirichletVelocityMask, & floatingEdges, layerThicknessFractions, lowerSurface, thickness, err) @@ -805,13 +809,13 @@ subroutine li_velocity_external_write_albany_mesh(domain) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal', floatingBasalMassBal) call mpas_pool_get_array(geometryPool, 'vertexMask', vertexMask, timeLevel = 1) - call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask, timeLevel = 1) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) ! Velocity variables call mpas_pool_get_array(velocityPool, 'beta', beta) - call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges) + call mpas_pool_get_array(velocityPool, 'floatingEdges', floatingEdges, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'dirichletVelocityMask', dirichletVelocityMask, timeLevel = 1) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor)