From 165fb8a98de9ecb750124d6829a654e806cc5c59 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Tue, 20 Apr 2021 17:51:26 -0600 Subject: [PATCH] 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 | 2 +- libglide/glide_setup.F90 | 10 ++- libglimmer/glimmer_paramets.F90 | 3 +- libglimmer/glimmer_physcon.F90 | 8 ++- libglimmer/glimmer_scales.F90 | 2 +- libglissade/glissade_basal_traction.F90 | 9 ++- libglissade/glissade_therm.F90 | 14 ++-- libglissade/glissade_velo.F90 | 3 - libglissade/glissade_velo_higher.F90 | 85 ++++++++++++++++--------- libglissade/glissade_velo_sia.F90 | 13 ++-- 10 files changed, 93 insertions(+), 56 deletions(-) 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