Fix inadvertent coupling when config_damage_rheology_coupling is false

MALI-Dev · Sep 28, 2023 · db25426 · db25426
1 parent 3f9c17a
commit db25426
Showing 1 changed file with 17 additions and 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
@@ -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