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 1df2391cdb2c..cef805e8d965 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 @@ -3517,11 +3517,10 @@ subroutine li_calculate_damage(domain, err) ! damage is always larger than the Nye value for initialization of damage evolution. do iCell = 1, nCells - if ((li_mask_is_grounded_ice(cellMask(iCell))) .or. (thickness(iCell) .eq. 0.0_RKIND)) then + if ( thickness(iCell) .eq. 0.0_RKIND ) then damage(iCell) = 0.0_RKIND end if end do - ! the damage of grounded ice is kept as 0, as the strain rate calculation is only valid for ice shelf ! Options for initializing damage where ice goes afloat: if (trim(config_damage_gl_setting) == 'extrapolate') then @@ -3554,7 +3553,7 @@ subroutine li_calculate_damage(domain, err) if (li_mask_is_grounding_line(cellMask(iCell))) then do iNeighbor = 1, nEdgesOnCell(iCell) jCell = cellsOnCell(iNeighbor, iCell) - if (li_mask_is_floating_ice(cellMask(jCell))) then + if ( li_mask_is_floating_ice(cellMask(jCell)) .and. (damage(jCell) < damageNye(jCell)) ) then damage(jCell) = damageNye(jCell) end if end do @@ -3645,7 +3644,6 @@ subroutine li_finalize_damage_after_advection(domain, err) type (mpas_pool_type), pointer :: meshPool type (mpas_pool_type), pointer :: velocityPool type (mpas_pool_type), pointer :: scratchPool - real(kind=RKIND), pointer :: config_damage_stiffness_min logical, pointer :: config_damage_rheology_coupling logical, pointer :: config_preserve_damage logical, pointer :: config_print_calving_info @@ -3667,7 +3665,6 @@ subroutine li_finalize_damage_after_advection(domain, err) err = 0 - call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min) call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling) call mpas_pool_get_config(liConfigs, 'config_damage_gl_setting', config_damage_gl_setting) call mpas_pool_get_config(liConfigs, 'config_preserve_damage', config_preserve_damage) @@ -3715,7 +3712,7 @@ subroutine li_finalize_damage_after_advection(domain, err) do iCell = 1, nCells - if (li_mask_is_grounded_ice(cellMask(iCell)) .or. .not. li_mask_is_ice(cellMask(iCell))) then + if (.not. li_mask_is_ice(cellMask(iCell))) then damage(iCell) = 0.0_RKIND end if end do @@ -3751,7 +3748,7 @@ subroutine li_finalize_damage_after_advection(domain, err) if (li_mask_is_grounding_line(cellMask(iCell))) then do iNeighbor = 1, nEdgesOnCell(iCell) jCell = cellsOnCell(iNeighbor, iCell) - if (li_mask_is_floating_ice(cellMask(jCell))) then + if ( li_mask_is_floating_ice(cellMask(jCell)) .and. (damage(jCell) block % next enddo @@ -3858,11 +3843,11 @@ subroutine apply_calving_damage_threshold(meshPool, geometryPool, scratchPool, d growMask(:) = 0 ! define seed and grow masks for flood fill. - where ( li_mask_is_dynamic_margin(cellMask) .and. (damage .ge. config_damage_calving_threshold) ) + where ( li_mask_is_dynamic_margin(cellMask) .and. li_mask_is_floating_ice(cellMask) .and. (damage .ge. config_damage_calving_threshold) ) seedMask = 1 end where - where ( seedMask == 0 .and. (damage .ge. config_damage_calving_threshold) ) + where ( seedMask == 0 .and. li_mask_is_floating_ice(cellMask) .and. (damage .ge. config_damage_calving_threshold) ) growMask = 1 end where 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 cf80b2497ebe..51d1cb92d7bc 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 @@ -909,12 +909,15 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo real (kind=RKIND), dimension(:,:), pointer :: flowParamA real (kind=RKIND), dimension(:), pointer :: effectiveViscosity real (kind=RKIND), dimension(:), pointer :: stiffnessFactor + real (kind=RKIND), dimension(:), pointer :: damage type (field1dReal), pointer :: meanFlowParamAVar real (kind=RKIND), dimension(:), pointer :: meanFlowParamA integer, pointer :: nVertLevels integer, pointer :: nCells + logical, pointer :: config_damage_rheology_coupling + real (kind=RKIND), pointer :: config_damage_stiffness_min integer :: iCell @@ -926,6 +929,8 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo call mpas_pool_get_array(velocityPool, 'yvelmean', yvelmean) call mpas_pool_get_config(liConfigs, 'config_flowLawExponent', config_flowLawExponent) + call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min) + call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -949,6 +954,7 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo call mpas_pool_get_array(velocityPool, 'flowParamA', flowParamA) call mpas_pool_get_array(velocityPool, 'stiffnessFactor', stiffnessFactor) + call mpas_pool_get_array(geometryPool, 'damage', damage) call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'effectiveViscosity', effectiveViscosity) @@ -984,9 +990,18 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo if ( (eEff == 0.0_RKIND) .or. (meanFlowParamA(iCell) == 0.0_RKIND) ) then effectiveViscosity(iCell) = 0.0_RKIND else - effectiveViscosity(iCell) = & - 0.5_RKIND*stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* & - eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent) + if ( config_damage_rheology_coupling ) then + effectiveViscosity(iCell) = & + 0.5_RKIND * max(1.0_RKIND - damage(iCell), config_damage_stiffness_min) * & + stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* & + eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent) + else + effectiveViscosity(iCell) = & + 0.5_RKIND * stiffnessFactor(iCell)*meanFlowParamA(iCell)** & + (-1.0_RKIND/config_flowLawExponent)* & + eEff**((1.0_RKIND-config_flowLawExponent)/config_flowLawExponent) + endif + endif enddo diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F index 077ff4c1773f..c87738181269 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_velocity_external.F @@ -449,6 +449,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro normalVelocity, uReconstructX, uReconstructY, uReconstructZ real (kind=RKIND), dimension(:,:), pointer :: temperature real (kind=RKIND), dimension(:), pointer :: stiffnessFactor + real (kind=RKIND), dimension(:), pointer :: damage real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: effectivePressureLimited real (kind=RKIND), dimension(:), pointer :: muFriction @@ -459,8 +460,10 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro logical, pointer :: config_always_compute_fem_grid logical, pointer :: config_output_external_velocity_solver_data logical, pointer :: config_nonconvergence_error + logical, pointer :: config_damage_rheology_coupling real (kind=RKIND), pointer :: config_ice_density real (kind=RKIND), pointer :: config_effective_pressure_max + real (kind=RKIND), pointer :: config_damage_stiffness_min integer, pointer :: anyDynamicVertexMaskChanged integer, pointer :: dirichletMaskChanged integer, pointer :: nEdges @@ -485,6 +488,8 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_config(liConfigs, 'config_nonconvergence_error', config_nonconvergence_error) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_effective_pressure_max', config_effective_pressure_max) + call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min) + call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling) ! Mesh variables call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) @@ -502,6 +507,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) call mpas_pool_get_array(geometryPool, 'edgeMask', edgeMask) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) + call mpas_pool_get_array(geometryPool, 'damage', damage) ! Thermal variables call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -595,12 +601,21 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro call mpas_timer_start("velocity_solver_solve_FO") - call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & - betaSolve, sfcMassBal, temperature, stiffnessFactor, & - effectivePressureLimited, muFriction, & - uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 - normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values - deltat, albanyVelocityError) ! return values + if ( config_damage_rheology_coupling ) then + call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & + betaSolve, sfcMassBal, temperature, max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor, & + effectivePressureLimited, muFriction, & + uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 + normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values + deltat, albanyVelocityError) ! return values + else + call velocity_solver_solve_FO(bedTopography, lowerSurface, thickness, & + betaSolve, sfcMassBal, temperature, stiffnessFactor, & + effectivePressureLimited, muFriction, & + uReconstructX, uReconstructY, & ! Dirichlet boundary values to apply where dirichletVelocityMask=1 + normalVelocity, drivingStressVert, dissipationVertexField % array, uReconstructX, uReconstructY, & ! return values + deltat, albanyVelocityError) ! return values + endif call mpas_timer_stop("velocity_solver_solve_FO") if (albanyVelocityError == 1) then @@ -770,7 +785,8 @@ subroutine li_velocity_external_write_albany_mesh(domain) !----------------------------------------------------------------- ! local variables !----------------------------------------------------------------- - logical, pointer :: config_write_albany_ascii_mesh + logical, pointer :: config_write_albany_ascii_mesh, & + config_damage_rheology_coupling character (len=StrKIND), pointer :: config_velocity_solver real (kind=RKIND), dimension(:), pointer :: & bedTopography, lowerSurface, upperSurface, layerThicknessFractions, beta @@ -784,12 +800,12 @@ subroutine li_velocity_external_write_albany_mesh(domain) real (kind=RKIND), dimension(:), pointer :: surfaceAirTemperature, basalHeatFlux integer, dimension(:), pointer :: vertexMask, cellMask, edgeMask, indexToCellID real (kind=RKIND), dimension(:,:), pointer :: layerThickness - real (kind=RKIND), dimension(:), pointer :: stiffnessFactor + real (kind=RKIND), dimension(:), pointer :: stiffnessFactor, damage real (kind=RKIND), dimension(:), pointer :: effectivePressure real (kind=RKIND), dimension(:), pointer :: muFriction integer, dimension(:,:), pointer :: dirichletVelocityMask type (mpas_pool_type), pointer :: meshPool, geometryPool, thermalPool, observationsPool, velocityPool, scratchPool, hydroPool - real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density + real (kind=RKIND), pointer :: config_sea_level, config_ice_density, config_ocean_density, config_damage_stiffness_min integer :: iCell integer :: err @@ -810,6 +826,8 @@ subroutine li_velocity_external_write_albany_mesh(domain) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) + call mpas_pool_get_config(liConfigs, 'config_damage_stiffness_min', config_damage_stiffness_min) + call mpas_pool_get_config(liConfigs, 'config_damage_rheology_coupling', config_damage_rheology_coupling) if (trim(config_velocity_solver) /= 'FO') then call mpas_log_write("config_velocity solver needs to be set to 'FO' for config_write_albany_ascii_mesh to work.", & @@ -836,6 +854,7 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! Geometry variables call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'damage', damage) call mpas_pool_get_array(geometryPool, 'bedTopography', bedTopography) call mpas_pool_get_array(geometryPool, 'lowerSurface', lowerSurface) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) @@ -896,17 +915,29 @@ subroutine li_velocity_external_write_albany_mesh(domain) ! call the C++ routine to write the mesh call mpas_log_write("Writing Albany ASCII mesh.", flushNow=.true.) - call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, & - beta, temperature, & - surfaceAirTemperature, basalHeatFlux, & - stiffnessFactor, & - effectivePressure, muFriction, & - thickness, thicknessUncertainty, & - sfcMassBal, sfcMassBalUncertainty, & - floatingBasalMassBal, floatingBasalMassBalUncertainty, & - observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, & - observedThicknessTendency, observedThicknessTendencyUncertainty) - + if ( config_damage_rheology_coupling ) then + call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, & + beta, temperature, & + surfaceAirTemperature, basalHeatFlux, & + max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor, & + effectivePressure, muFriction, & + thickness, thicknessUncertainty, & + sfcMassBal, sfcMassBalUncertainty, & + floatingBasalMassBal, floatingBasalMassBalUncertainty, & + observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, & + observedThicknessTendency, observedThicknessTendencyUncertainty) + else + call write_ascii_mesh(indexToCellID, bedTopography, lowerSurface, & + beta, temperature, & + surfaceAirTemperature, basalHeatFlux, & + stiffnessFactor, & + effectivePressure, muFriction, & + thickness, thicknessUncertainty, & + sfcMassBal, sfcMassBalUncertainty, & + floatingBasalMassBal, floatingBasalMassBalUncertainty, & + observedSurfaceVelocityX, observedSurfaceVelocityY, observedSurfaceVelocityUncertainty, & + observedThicknessTendency, observedThicknessTendencyUncertainty) + endif !---- call interface_reset_stdout()