diff --git a/libglide/felix_dycore_interface.F90 b/libglide/felix_dycore_interface.F90
index f702d9a1..ce1ed0dc 100644
--- a/libglide/felix_dycore_interface.F90
+++ b/libglide/felix_dycore_interface.F90
@@ -146,7 +146,7 @@ end subroutine felix_velo_init
subroutine felix_velo_driver(model)
use glimmer_global, only : dp
- use glimmer_physcon, only: gn, scyr
+ use glimmer_physcon, only: scyr
use glimmer_paramets, only: thk0, len0, vel0, vis0
use glimmer_log
use glide_types
diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90
index 3b967d05..909244fb 100644
--- a/libglide/glide_setup.F90
+++ b/libglide/glide_setup.F90
@@ -160,7 +160,7 @@ end subroutine glide_printconfig
subroutine glide_scale_params(model)
!> scale parameters
use glide_types
- use glimmer_physcon, only: scyr, gn
+ use glimmer_physcon, only: scyr
use glimmer_paramets, only: thk0, tim0, len0, vel0, vis0, acc0, tau0
implicit none
@@ -1996,7 +1996,7 @@ subroutine handle_parameters(section, model)
use glimmer_config
use glide_types
use glimmer_log
- use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt
+ use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt, n_glen
implicit none
type(ConfigSection), pointer :: section
@@ -2021,6 +2021,7 @@ subroutine handle_parameters(section, model)
call GetValue(section,'lhci', lhci)
call GetValue(section,'trpt', trpt)
#endif
+ call GetValue(section,'n_glen', n_glen)
loglevel = GM_levels-GM_ERROR
call GetValue(section,'log_level',loglevel)
@@ -2206,7 +2207,7 @@ end subroutine handle_parameters
subroutine print_parameters(model)
- use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav
+ use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav, n_glen
use glide_types
use glimmer_log
implicit none
@@ -2371,6 +2372,9 @@ subroutine print_parameters(model)
write(message,*) 'triple point of water (K) : ', trpt
call write_log(message)
+ write(message,*) 'Glen flow law exponent : ', n_glen
+ call write_log(message)
+
write(message,*) 'geothermal flux (W/m^2) : ', model%paramets%geot
call write_log(message)
diff --git a/libglimmer/glimmer_paramets.F90 b/libglimmer/glimmer_paramets.F90
index 019f44e6..aa8b595d 100644
--- a/libglimmer/glimmer_paramets.F90
+++ b/libglimmer/glimmer_paramets.F90
@@ -33,7 +33,7 @@
module glimmer_paramets
use glimmer_global, only : dp
- use glimmer_physcon, only : scyr, gn
+ use glimmer_physcon, only : scyr
implicit none
save
@@ -118,6 +118,7 @@ module glimmer_paramets
real(dp), parameter :: grav_glam = 9.81d0 ! m s^{-2}
! GLAM scaling parameters; units are correct if thk0 has units of meters
+ integer, parameter :: gn = 3 ! Glen flow exponent; fixed at 3 for purposes of setting vis0
real(dp), parameter :: tau0 = rhoi_glam*grav_glam*thk0 ! stress scale in GLAM ( Pa )
real(dp), parameter :: evs0 = tau0 / (vel0/len0) ! eff. visc. scale in GLAM ( Pa s )
real(dp), parameter :: vis0 = tau0**(-gn) * (vel0/len0) ! rate factor scale in GLAM ( Pa^-3 s^-1 )
diff --git a/libglimmer/glimmer_physcon.F90 b/libglimmer/glimmer_physcon.F90
index 0c990d29..f697bf3e 100644
--- a/libglimmer/glimmer_physcon.F90
+++ b/libglimmer/glimmer_physcon.F90
@@ -79,11 +79,17 @@ module glimmer_physcon
! real(dp) :: trpt = 273.16d0 !< Triple point of water (K)
#endif
+ ! The Glen flow-law exponent, n_glen, is used in Glissade.
+ ! It is not a parameter, because the default can be overridden in the config file.
+ ! TODO: Allow setting n_glen independently for each ice sheet instance?
+ ! Note: Earlier code used an integer parameter, gn = 3, for all flow-law calculations.
+ ! For backward compatiblity, gn = 3 is retained for Glide.
+ real(dp) :: n_glen = 3.0d0 !< Exponent in Glen's flow law; user-configurable real(dp) in Glissade
+ integer, parameter :: gn = 3 !< Exponent in Glen's flow law; fixed integer parameter in Glide
real(dp),parameter :: celsius_to_kelvin = 273.15d0 !< Note: Not quite equal to trpt
real(dp),parameter :: scyr = 31536000.d0 !< Number of seconds in a year of exactly 365 days
real(dp),parameter :: rhom = 3300.0d0 !< The density of magma(?) (kg m-3)
real(dp),parameter :: rhos = 2600.0d0 !< The density of solid till (kg m$^{-3}$)
- integer, parameter :: gn = 3 !< The power dependency of Glen's flow law.
real(dp),parameter :: actenh = 139.0d3 !< Activation energy in Glen's flow law for \f$T^{*}\geq263\f$K. (J mol-1)
real(dp),parameter :: actenl = 60.0d3 !< Activation energy in Glen's flow law for \f$T^{*}<263\f$K. (J mol-1)
real(dp),parameter :: arrmlh = 1.733d3 !< Constant of proportionality in Arrhenius relation
diff --git a/libglimmer/glimmer_scales.F90 b/libglimmer/glimmer_scales.F90
index f95c6a86..380ff2b8 100644
--- a/libglimmer/glimmer_scales.F90
+++ b/libglimmer/glimmer_scales.F90
@@ -45,7 +45,7 @@ subroutine glimmer_init_scales
! set scale factors for I/O (can't have non-integer powers)
- use glimmer_physcon, only : scyr, gn
+ use glimmer_physcon, only : scyr
use glimmer_paramets, only : thk0, tim0, vel0, vis0, len0, acc0, tau0, evs0
implicit none
diff --git a/libglissade/glissade_basal_traction.F90 b/libglissade/glissade_basal_traction.F90
index 69f1fd75..5dc1e2cb 100644
--- a/libglissade/glissade_basal_traction.F90
+++ b/libglissade/glissade_basal_traction.F90
@@ -440,9 +440,12 @@ subroutine calcbeta (whichbabc, &
! If this factor is not present in the input file, it is set to 1 everywhere.
! Compute beta
- ! gn = Glen's n from physcon module
- beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/gn - 1.0d0) * &
- (speed(:,:) + basal_physics%effecpress_stag(:,:)**gn * big_lambda)**(-1.0d0/gn)
+ ! Note: Where this equation has powerlaw_m, we used to have Glen's flow exponent n,
+ ! following the notation of Leguy et al. (2014).
+ ! Changed to powerlaw_m to be consistent with the Schoof and Tsai laws.
+ m = basal_physics%powerlaw_m
+ beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/m - 1.0d0) * &
+ (speed(:,:) + basal_physics%effecpress_stag(:,:)**m * big_lambda)**(-1.0d0/m)
! If c_space_factor /= 1.0 everywhere, then multiply beta by c_space_factor
if (maxval(abs(basal_physics%c_space_factor_stag(:,:) - 1.0d0)) > tiny(0.0d0)) then
diff --git a/libglissade/glissade_therm.F90 b/libglissade/glissade_therm.F90
index 10b1fcca..f6364650 100644
--- a/libglissade/glissade_therm.F90
+++ b/libglissade/glissade_therm.F90
@@ -1640,10 +1640,11 @@ subroutine glissade_enthalpy_matrix_elements(dttem, &
! At each temperature point, compute the temperature part of the enthalpy.
! enth_T = enth for cold ice, enth_T < enth for temperate ice
- enth_T(0) = rhoi*shci*temp(0) !WHL - not sure enth_T(0) is needed
- do up = 1, upn
+ do up = 1, upn-1
enth_T(up) = (1.d0 - waterfrac(up)) * rhoi*shci*temp(up)
enddo
+ enth_T(0) = rhoi*shci*temp(0)
+ enth_T(up) = rhoi*shci*temp(up)
!WHL - debug
if (verbose_column) then
@@ -2250,8 +2251,7 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, &
! Compute the dissipation source term associated with strain heating,
! based on the shallow-ice approximation.
- use glimmer_physcon, only : gn ! Glen's n
- use glimmer_physcon, only: rhoi, shci, grav
+ use glimmer_physcon, only: rhoi, shci, grav, n_glen
integer, intent(in) :: ewn, nsn, upn ! grid dimensions
@@ -2267,12 +2267,14 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, &
real(dp), dimension(:,:,:), intent(out) :: &
dissip ! interior heat dissipation (deg/s)
- integer, parameter :: p1 = gn + 1
-
integer :: ew, ns
real(dp), dimension(upn-1) :: sia_dissip_fact ! factor in SIA dissipation calculation
real(dp) :: geom_fact ! geometric factor
+ real(dp) :: p1 ! exponent = n_glen + 1
+
+ p1 = n_glen + 1.0d0
+
! Two methods of doing this calculation:
! 1. find dissipation at u-pts and then average
! 2. find dissipation at H-pts by averaging quantities from u-pts
diff --git a/libglissade/glissade_velo.F90 b/libglissade/glissade_velo.F90
index bb6e6e38..a9dadd17 100644
--- a/libglissade/glissade_velo.F90
+++ b/libglissade/glissade_velo.F90
@@ -43,9 +43,6 @@ subroutine glissade_velo_driver(model)
! Glissade higher-order velocity driver
- use glimmer_global, only : dp
- use glimmer_physcon, only: gn, scyr
- use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0, evs0
use glimmer_log
use glide_types
use glissade_velo_higher, only: glissade_velo_higher_solve
diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90
index 033ca254..e9845437 100644
--- a/libglissade/glissade_velo_higher.F90
+++ b/libglissade/glissade_velo_higher.F90
@@ -57,7 +57,7 @@
module glissade_velo_higher
use glimmer_global, only: dp
- use glimmer_physcon, only: gn, rhoi, rhoo, grav, scyr, pi
+ use glimmer_physcon, only: n_glen, rhoi, rhoo, grav, scyr, pi
use glimmer_paramets, only: eps08, eps10, thk0, len0, tim0, tau0, vel0, vis0, evs0
use glimmer_paramets, only: vel_scale, len_scale ! used for whichefvs = HO_EFVS_FLOWFACT
use glimmer_log
@@ -2082,7 +2082,7 @@ subroutine glissade_velo_higher_solve(model, &
! gn = exponent in Glen's flow law (= 3 by default)
do k = 1, nz-1
if (flwa(k,i,j) > 0.0d0) then
- flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/real(gn,dp))
+ flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/n_glen)
endif
enddo
enddo
@@ -4222,6 +4222,7 @@ subroutine glissade_velo_higher_solve(model, &
usrf, &
dusrf_dx, dusrf_dy, &
flwa, efvs, &
+ whichefvs, efvs_constant, &
whichgradient_margin, &
max_slope, &
uvel, vvel)
@@ -6426,6 +6427,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, &
usrf, &
dusrf_dx, dusrf_dy, &
flwa, efvs, &
+ whichefvs, efvs_constant, &
whichgradient_margin, &
max_slope, &
uvel, vvel)
@@ -6486,6 +6488,12 @@ subroutine compute_3d_velocity_L1L2(nx, ny, &
flwa, & ! temperature-based flow factor A, Pa^{-n} yr^{-1}
efvs ! effective viscosity, Pa yr
+ integer, intent(in) :: &
+ whichefvs ! option for effective viscosity calculation
+
+ real(dp), intent(in) :: &
+ efvs_constant ! constant value of effective viscosity (Pa yr)
+
integer, intent(in) :: &
whichgradient_margin ! option for computing gradient at ice margin
! 0 = include all neighbor cells in gradient calculation
@@ -6840,7 +6848,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, &
! Compute vertical integration factor at each active vertex
! This is int_b_to_z{-2 * A * tau^2 * rho*g*(s-z) * dz},
- ! similar to the factor computed in Glide and glissade_velo_sia..
+ ! similar to the factor computed in Glide and glissade_velo_sia.
! Note: tau_xz ~ rho*g*(s-z)*ds_dx; ds_dx term is computed on edges below
do j = 1, ny-1
@@ -6921,9 +6929,27 @@ subroutine compute_3d_velocity_L1L2(nx, ny, &
tau_eff_sq = stagtau_parallel_sq(i,j) &
+ tau_xz(k,i,j)**2 + tau_yz(k,i,j)**2
- ! Note: This formula is correct for any value of Glen's n, but currently efvs is computed
- ! only for gn = 3 (in which case (n-1)/2 = 1).
- fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((gn-1.d0)/2.d0) * (sigma(k+1) - sigma(k))*stagthck(i,j)
+ ! Note: The first formula below is correct for whichefvs = 2 (efvs computed from effective strain rate),
+ ! but not for whichefvs = 0 (constant efvs) or whichefvs = 1 (multiple of flow factor).
+ ! For these options we need a modified formula.
+ !
+ ! Recall: efvs = 1/2 * A^(-1/n) * eps_e^[(1-n)/n]
+ ! = 1/2 * A^(-1/n) * [A tau_e^n]^[(1-n)/n]
+ ! = 1/2 * A^(-1) * tau_e^(1-n)
+ ! => 1/efvs = 2 * A * tau_e(n-1)
+ !
+ ! Thus, for options 0 and 1, we can replace 2 * A * tau_e^(n-1) below with 1/efvs.
+
+ if (whichefvs == HO_EFVS_NONLINEAR) then
+ fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((n_glen-1.d0)/2.d0) &
+ * (sigma(k+1) - sigma(k))*stagthck(i,j)
+ else ! HO_EFVS_CONSTANT, HO_EFVS_FLOWFACT
+ if (efvs(k,i,j) > 0.0d0) then
+ fact = (sigma(k+1) - sigma(k))*stagthck(i,j) / efvs(k,i,j)
+ else
+ fact = 0.0d0
+ endif
+ endif
! reset velocity to prescribed basal value if Dirichlet condition applies
! else compute velocity at this level
@@ -7875,15 +7901,6 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, &
integer, intent(in) :: i, j, k, p
- !----------------------------------------------------------------
- ! Local parameters
- !----------------------------------------------------------------
-
- real(dp), parameter :: &
- p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation
- p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation
-
-
!----------------------------------------------------------------
! Local variables
!----------------------------------------------------------------
@@ -7896,8 +7913,14 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, &
integer :: n
+ real(dp) :: &
+ p_effstr ! exponent (1-n)/n in effective viscosity relation
+
real(dp), parameter :: p2 = -1.d0/3.d0
+ ! Set exponent that depends on Glen's exponent
+ p_effstr = (1.d0 - n_glen)/n_glen
+
select case(whichefvs)
case(HO_EFVS_CONSTANT)
@@ -7988,11 +8011,11 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, &
! Compute effective viscosity (PGB 2012, eq. 4)
! Units: flwafact has units Pa yr^{1/n}
! effstrain has units yr^{-1}
- ! p2_effstr = (1-n)/(2n)
- ! = -1/3 for n=3
+ ! p_effstr = (1-n)/n
+ ! = -2/3 for n=3
! Thus efvs has units Pa yr
- efvs = flwafact * effstrainsq**p2_effstr
+ efvs = flwafact * effstrainsq**(p_effstr/2.d0)
if (verbose_efvs .and. this_rank==rtest .and. i==itest .and. j==jtest .and. k==ktest .and. p==ptest) then
print*, ' '
@@ -8081,8 +8104,8 @@ subroutine compute_effective_viscosity_L1L2(whichefvs,
! Local parameters
!----------------------------------------------------------------
- real(dp), parameter :: &
- p_effstr = (1.d0 - real(gn,dp)) / real(gn,dp) ! exponent (1-n)/n in effective viscosity relation
+ real(dp) :: &
+ p_effstr ! exponent (1-n)/n in effective viscosity relation
!----------------------------------------------------------------
! Local variables
@@ -8107,6 +8130,9 @@ subroutine compute_effective_viscosity_L1L2(whichefvs,
integer :: n, k
+ ! Set exponent that depends on Glen's exponent
+ p_effstr = (1.d0 - n_glen) / n_glen
+
select case(whichefvs)
case(HO_EFVS_CONSTANT)
@@ -8125,7 +8151,7 @@ subroutine compute_effective_viscosity_L1L2(whichefvs,
!
! Units: flwafact has units Pa yr^{1/n}
! effstrain has units yr^{-1}
- ! p_effstr = (1-n)/n
+ ! p_effstr = (1-n)/n
! = -2/3 for n=3
! Thus efvs has units Pa yr
@@ -8320,14 +8346,6 @@ subroutine compute_effective_viscosity_diva(whichefvs,
integer, intent(in) :: i, j, p
- !----------------------------------------------------------------
- ! Local parameters
- !----------------------------------------------------------------
-
- real(dp), parameter :: &
- p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation
- p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation
-
!----------------------------------------------------------------
! Local variables
!----------------------------------------------------------------
@@ -8346,11 +8364,17 @@ subroutine compute_effective_viscosity_diva(whichefvs,
integer :: n, k
real(dp) :: du_dz, dv_dz
+ real(dp) :: &
+ p_effstr ! exponent (1-n)/n in effective viscosity relation
+
!WHL - For ISMIP-HOM, the cubic solve is not robust. It leads to oscillations
! in successive iterations between uvel_2d/vvel_2d and btractx/btracty
!TODO - Remove the cubic solve for efvs, unless we find a way to make it robust?
logical, parameter :: cubic = .false.
+ ! Set exponent that depends on Glen's exponent
+ p_effstr = (1.d0 - n_glen)/n_glen
+
select case(whichefvs)
case(HO_EFVS_CONSTANT)
@@ -8493,7 +8517,8 @@ subroutine compute_effective_viscosity_diva(whichefvs,
effstrainsq = effstrain_min**2 &
+ du_dx**2 + dv_dy**2 + du_dx*dv_dy + 0.25d0*(dv_dx + du_dy)**2 &
+ 0.25d0 * (du_dz**2 + dv_dz**2)
- efvs(k) = flwafact(k) * effstrainsq**p2_effstr
+ efvs(k) = flwafact(k) * effstrainsq**(p_effstr/2.d0)
+
enddo
endif ! cubic
diff --git a/libglissade/glissade_velo_sia.F90 b/libglissade/glissade_velo_sia.F90
index 66884aaf..1efcb554 100644
--- a/libglissade/glissade_velo_sia.F90
+++ b/libglissade/glissade_velo_sia.F90
@@ -55,7 +55,7 @@
module glissade_velo_sia
use glimmer_global, only: dp
- use glimmer_physcon, only: gn, rhoi, grav, scyr
+ use glimmer_physcon, only: n_glen, rhoi, grav, scyr
use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0
! use glimmer_log, only: write_log
@@ -881,16 +881,15 @@ subroutine glissade_velo_sia_interior(nx, ny, nz, &
if (stagthck(i,j) > thklim) then
- siafact = 2.d0 * (rhoi*grav)**gn * stagthck(i,j)**(gn+1) &
- * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((gn-1)/2)
-
+ siafact = 2.d0 * (rhoi*grav)**n_glen * stagthck(i,j)**(n_glen+1) &
+ * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((n_glen-1)/2)
vintfact(nz,i,j) = 0.d0
do k = nz-1, 1, -1
- vintfact(k,i,j) = vintfact(k+1,i,j) - &
- siafact * stagflwa(k,i,j) &
- * ((sigma(k) + sigma(k+1))/2.d0) ** gn &
+ vintfact(k,i,j) = vintfact(k+1,i,j) - &
+ siafact * stagflwa(k,i,j) &
+ * ((sigma(k) + sigma(k+1))/2.d0) ** n_glen &
* (sigma(k+1) - sigma(k))
enddo ! k