From e5ce4362fc3aec0dd1472f3c47a4be030d001f39 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Tue, 22 Aug 2023 17:53:01 -0600 Subject: [PATCH 1/7] Add damage coupling in velocity solver Pass max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor to velocity solver instead of just stiffnessFactor. --- .../src/mode_forward/mpas_li_calving.F | 12 ------------ .../src/mode_forward/mpas_li_velocity.F | 7 ++++++- .../src/mode_forward/mpas_li_velocity_external.F | 6 +++++- 3 files changed, 11 insertions(+), 14 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 1df2391cdb2c..a18d0bb7bad2 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 @@ -3769,18 +3769,6 @@ subroutine li_finalize_damage_after_advection(domain, err) damage = 1.0_RKIND end where - if (config_damage_rheology_coupling) then - do iCell = 1, nCells - if (li_mask_is_floating_ice(cellMask(iCell))) then - stiffnessFactor(iCell) = 1.0_RKIND - damage(iCell) - if (stiffnessFactor(iCell) < config_damage_stiffness_min) then - stiffnessFactor(iCell) = config_damage_stiffness_min - end if - end if - end do - end if - - block => block % next enddo 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..0d40a31736ec 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,14 @@ 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 + real (kind=RKIND), pointer :: config_damage_stiffness_min integer :: iCell @@ -926,6 +928,7 @@ 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_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -949,6 +952,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) @@ -985,7 +989,8 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo effectiveViscosity(iCell) = 0.0_RKIND else effectiveViscosity(iCell) = & - 0.5_RKIND*stiffnessFactor(iCell)*meanFlowParamA(iCell)**(-1.0_RKIND/config_flowLawExponent)* & + 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) 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..a1458ba7fa90 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 @@ -461,6 +462,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro logical, pointer :: config_nonconvergence_error 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 +487,7 @@ 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) ! Mesh variables call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) @@ -502,6 +505,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) @@ -596,7 +600,7 @@ 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, & + 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 From b8a07aa1d1edbedb5d3f8a61ec4090ae9c6711d9 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Thu, 24 Aug 2023 10:09:06 -0600 Subject: [PATCH 2/7] Fix inadvertent coupling when config_damage_rheology_coupling is false Fix inadvertent coupling when config_damage_rheology_coupling is false --- .../mode_forward/mpas_li_velocity_external.F | 23 ++++++++++++++----- 1 file changed, 17 insertions(+), 6 deletions(-) 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 a1458ba7fa90..24bc98103b8d 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 @@ -460,6 +460,7 @@ 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 @@ -488,6 +489,7 @@ subroutine li_velocity_external_solve(meshPool, geometryPool, thermalPool, hydro 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) @@ -599,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, 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 + 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 From ba9c5a9d9676d9870ba6d6384ea704abcc7585ee Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 25 Aug 2023 08:45:26 -0600 Subject: [PATCH 3/7] Don't include damage in uncoupled effectiveViscosity calculation Remove damage factor from calculation for uncoupled effectiveViscosity --- .../src/mode_forward/mpas_li_velocity.F | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) 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 0d40a31736ec..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 @@ -916,6 +916,7 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo integer, pointer :: nVertLevels integer, pointer :: nCells + logical, pointer :: config_damage_rheology_coupling real (kind=RKIND), pointer :: config_damage_stiffness_min integer :: iCell @@ -929,6 +930,7 @@ subroutine calculate_strain_rates_and_stresses(meshPool, geometryPool, thermalPo 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) @@ -988,10 +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 * 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) + 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 From 9a94b05da4607baafd68979adc1ed9ab1a2ffc21 Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 13 Oct 2023 16:08:55 -0600 Subject: [PATCH 4/7] Some cleanup Remove config_damage_stiffness_min as a variable in li_calculate_damage, where it is no longer used. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 2 -- 1 file changed, 2 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index a18d0bb7bad2..8e05d4e0edcf 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 @@ -3645,7 +3645,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 +3666,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) From 27e07256cf1d0916226fb4b35f7a7dff852f7beb Mon Sep 17 00:00:00 2001 From: Trevor Ray Hillebrand Date: Fri, 29 Sep 2023 08:29:38 -0600 Subject: [PATCH 5/7] Do not zero out damage for grounded ice Do not zero out damage for grounded ice. Damage should remain small anyway, but forcing damage to be zero on grounded ice can cause unrealistically large gradients in damage at the grounding line. --- .../src/mode_forward/mpas_li_calving.F | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 8e05d4e0edcf..2f2db3803149 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 @@ -3713,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 @@ -3749,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) Date: Sun, 15 Oct 2023 15:48:10 -0600 Subject: [PATCH 6/7] Disable damage threshold calving for grounded ice Disable damage threshold calving for grounded ice now that damage is no longer set to zero for grounded ice. --- .../mpas-albany-landice/src/mode_forward/mpas_li_calving.F | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F index 2f2db3803149..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 @@ -3843,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 From 305eb6a37ddfb82bb26c2a3b9a5d22c49d99804a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Mon, 13 Nov 2023 18:40:56 -0800 Subject: [PATCH 7/7] Write effective stiffness to ascii mesh when coupling damage to viscosity Write effective stiffness (defined as max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor) to ascii mesh when coupling damage to viscosity. --- .../mode_forward/mpas_li_velocity_external.F | 44 +++++++++++++------ 1 file changed, 30 insertions(+), 14 deletions(-) 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 24bc98103b8d..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 @@ -785,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 @@ -799,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 @@ -825,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.", & @@ -851,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) @@ -911,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()