Write effective stiffness to ascii mesh when coupling damage to visco…

…sity Write effective stiffness (defined as max(1.0_RKIND - damage, config_damage_stiffness_min) * stiffnessFactor) to ascii mesh when coupling damage to viscosity.
MALI-Dev · Nov 14, 2023 · c9051a7 · c9051a7
1 parent e6689f6
commit c9051a7
Showing 1 changed file with 30 additions and 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
@@ -781,7 +781,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
@@ -795,12 +796,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
 
@@ -821,6 +822,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.", &
@@ -847,6 +850,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)
@@ -907,17 +911,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()