Made the Glen flow-law exponent 'n' a config variable

Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90
as an integer parameter, gn = 3.

With this commit, n is renamed n_glen (a more greppable name) for use in Glissade.
It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0.

Since n_glen is no longer a parameter, it can now be set in the config file like other
physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can
take different values in different models or benchmarking experiments.

To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90.
This parameter is now used only in the Glide dycore (e.g., glide_velo.F90).

In Glissade, I replaced gn with n_glen in several places:
(1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen.
    Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs.
(2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity.
(3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A.
(4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used
    to compute effective_viscosity for BP, DIVA, or L1L2.
(5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to
    tau**((n_glen-1.)/2.) in the vertical integral.

For (1) and (2), n_glen was previously assumed to be an integer.  Switching it to real(dp)
is answer-changing at the machine roundoff level for runs using the local SIA solver
in glissade_velo_sia.F90.

For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp).
For this reason, answers are BFB when using the BP, DIVA or L1L2 solver.

In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option,
HO_BABC_COULOMB_FRICTION.  Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default)
to be consistent with the expressions for beta in the Schoof and Tsai laws.

In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale,
Here, I defined a local integer parameter gn = 3 so that the scales are unchanged.

It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test.

Another minor change: I added some logic to the subroutine that computes L1L2 velocities.
For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from
the effective strain rate using an equation from Perego et al. (2012).
But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor),
this strain rate equation in the code does not apply. For these options, I added an
alternative that computes velocity in terms of the strain-rate-independent efvs.
This allows us to use L1L2 for problems with constant efvs (e.g., the slab test).

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

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)
 

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 )

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<SUP>-3</SUP>) 
   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<SUP>-1</SUP>)
   real(dp),parameter :: actenl = 60.0d3          !< Activation energy in Glen's flow law for \f$T^{*}<263\f$K. (J mol<SUP>-1</SUP>)  
   real(dp),parameter :: arrmlh = 1.733d3         !< Constant of proportionality in Arrhenius relation

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
 

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

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

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

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