Added a steady-state basal water routing scheme

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.F90

@@ -602,8 +602,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
           call calcbwat(model, &
                         model%options%whichbwat, &
                         model%basal_melt%bmlt, &
-                        model%temper%bwat, &
-                        model%temper%bwatflx, &
+                        model%basal_hydro%bwat, &
+                        model%basal_hydro%bwatflx, &
                         model%geometry%thck, &
                         model%geometry%topg, &
                         model%temper%temp(model%general%upn,:,:), &
@@ -612,8 +612,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice)
 
 
           ! This call is redundant for now, but is needed if the call to calcbwat is removed
-          call stagvarb(model%temper%bwat, &
-                        model%temper%stagbwat ,&
+          call stagvarb(model%basal_hydro%bwat, &
+                        model%basal_hydro%stagbwat ,&
                         model%general%ewn, &
                         model%general%nsn)
 
@@ -867,8 +867,8 @@ subroutine glide_tstep_p1(model,time)
        call calcbwat(model, &
                      model%options%whichbwat, &
                      model%basal_melt%bmlt, &
-                     model%temper%bwat, &
-                     model%temper%bwatflx, &
+                     model%basal_hydro%bwat, &
+                     model%basal_hydro%bwatflx, &
                      model%geometry%thck, &
                      model%geometry%topg, &
                      model%temper%temp(model%general%upn,:,:), &

libglide/glide_bwater.F90

@@ -195,8 +195,8 @@ subroutine calcbwat(model, which, bmlt, bwat, bwatflx, thck, topg, btem, floater
     end select
 
     ! now also calculate basal water in velocity (staggered) coord system
-    call stagvarb(model%temper%bwat, &
-                  model%temper%stagbwat ,&
+    call stagvarb(model%basal_hydro%bwat, &
+                  model%basal_hydro%stagbwat ,&
                   model%general%ewn, &
                   model%general%nsn)
 
@@ -429,7 +429,7 @@ subroutine pressure_wphi(thck,topg,N,wphi,thicklim,floater)
   !> Compute the pressure wphi at the base of the ice sheet according to
   !> ice overburden plus bed height minus effective pressure.
   !>
-  !> whpi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g)
+  !> wphi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g)
 
     use glimmer_physcon, only : rhoi,rhow,grav
     implicit none

libglide/glide_diagnostics.F90

@@ -200,9 +200,12 @@ subroutine glide_write_diag (model,       time)
          max_thck, max_thck_global,     &    ! max ice thickness (m)
          max_temp, max_temp_global,     &    ! max ice temperature (deg C)
          min_temp, min_temp_global,     &    ! min ice temperature (deg C)
-         max_bmlt, max_bmlt_global,     &    ! max basal melt rate (m/yr)
-         max_spd_sfc, max_spd_sfc_global,   &    ! max surface ice speed (m/yr)
-         max_spd_bas, max_spd_bas_global,   &    ! max basal ice speed (m/yr)
+         max_bmlt,                      &    ! max basal melt rate (m/yr)
+         max_bmlt_global,               &
+         max_bmlt_ground,               &    ! max basal melt rate, grounded ice (m/yr)
+         max_bmlt_ground_global,        &
+         max_spd_sfc, max_spd_sfc_global, &  ! max surface ice speed (m/yr)
+         max_spd_bas, max_spd_bas_global, &  ! max basal ice speed (m/yr)
          spd,                           &    ! speed
          thck_diag, usrf_diag,          &    ! local column diagnostics
          topg_diag, relx_diag,          &    
@@ -768,7 +771,8 @@ subroutine glide_write_diag (model,       time)
                     min_temp_global, imin_global, jmin_global, kmin_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
-    ! max basal melt rate
+    ! max applied basal melt rate
+    ! Usually, this will be for floating ice, if floating ice is present
     imax = 0
     jmax = 0
     max_bmlt = unphys_val
@@ -791,11 +795,39 @@ subroutine glide_write_diag (model,       time)
     ! Write to diagnostics only if nonzero
 
     if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then
-       write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/yr), i, j    ',   &
+       write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/y), i, j     ',   &
             max_bmlt_global*thk0*scyr/tim0, imax_global, jmax_global
        call write_log(trim(message), type = GM_DIAGNOSTIC)
     endif
 
+    ! max basal melt rate for grounded ice
+    imax = 0
+    jmax = 0
+    max_bmlt_ground = unphys_val
+
+    do j = lhalo+1, nsn-uhalo
+       do i = lhalo+1, ewn-uhalo
+          if (model%basal_melt%bmlt_ground(i,j) > max_bmlt_ground) then
+             max_bmlt_ground = model%basal_melt%bmlt_ground(i,j)
+             imax = i
+             jmax = j
+          endif
+       enddo
+    enddo
+
+    call parallel_reduce_maxloc(xin=max_bmlt_ground, xout=max_bmlt_ground_global, xprocout=procnum)
+    call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel)
+    call broadcast(imax_global, proc = procnum)
+    call broadcast(jmax_global, proc = procnum)
+
+    ! Write to diagnostics only if nonzero
+
+    if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then
+       write(message,'(a25,f24.16,2i6)') 'Max bmlt grnd (m/y), i, j',   &
+            max_bmlt_ground_global*thk0*scyr/tim0, imax_global, jmax_global
+       call write_log(trim(message), type = GM_DIAGNOSTIC)
+    endif
+
     ! max surface speed
     imax = 0
     jmax = 0
@@ -818,7 +850,7 @@ subroutine glide_write_diag (model,       time)
     call broadcast(imax_global, proc = procnum)
     call broadcast(jmax_global, proc = procnum)
 
-    write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/yr), i, j ',   &
+    write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/y), i, j  ',   &
                     max_spd_sfc_global*vel0*scyr, imax_global, jmax_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
@@ -843,7 +875,7 @@ subroutine glide_write_diag (model,       time)
     call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel)
     call broadcast(imax_global, proc = procnum)
     call broadcast(jmax_global, proc = procnum)
-    write(message,'(a25,f24.16,2i6)') 'Max base spd (m/yr), i, j',   &
+    write(message,'(a25,f24.16,2i6)') 'Max base spd (m/y), i, j ',   &
                     max_spd_bas_global*vel0*scyr, imax_global, jmax_global
     call write_log(trim(message), type = GM_DIAGNOSTIC)
 
@@ -884,7 +916,7 @@ subroutine glide_write_diag (model,       time)
           artm_diag = model%climate%artm_corrected(i,j)  ! artm_corrected = artm + artm_anomaly
           acab_diag = model%climate%acab_applied(i,j) * thk0*scyr/tim0
           bmlt_diag = model%basal_melt%bmlt_applied(i,j) * thk0*scyr/tim0
-          bwat_diag = model%temper%bwat(i,j) * thk0
+          bwat_diag = model%basal_hydro%bwat(i,j) * thk0
           bheatflx_diag = model%temper%bheatflx(i,j)
 
           temp_diag(:) = model%temper%temp(1:upn,i,j)          

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
@@ -745,6 +745,7 @@ subroutine handle_options(section, model)
     call GetValue(section,'cull_calving_front', model%options%cull_calving_front)
     call GetValue(section,'adjust_input_thickness', model%options%adjust_input_thickness)
     call GetValue(section,'smooth_input_topography', model%options%smooth_input_topography)
+    call GetValue(section,'smooth_input_usrf', model%options%smooth_input_usrf)
     call GetValue(section,'adjust_input_topography', model%options%adjust_input_topography)
     call GetValue(section,'read_lat_lon',model%options%read_lat_lon)
     call GetValue(section,'dm_dt_diag',model%options%dm_dt_diag)
@@ -786,6 +787,7 @@ subroutine handle_ho_options(section, model)
     call GetValue(section, 'which_ho_bmlt_inversion',     model%options%which_ho_bmlt_inversion)
     call GetValue(section, 'which_ho_bmlt_basin_inversion', model%options%which_ho_bmlt_basin_inversion)
     call GetValue(section, 'which_ho_bwat',               model%options%which_ho_bwat)
+    call GetValue(section, 'ho_flux_routing_scheme',      model%options%ho_flux_routing_scheme)
     call GetValue(section, 'which_ho_effecpress',         model%options%which_ho_effecpress)
     call GetValue(section, 'which_ho_resid',              model%options%which_ho_resid)
     call GetValue(section, 'which_ho_nonlinear',          model%options%which_ho_nonlinear)
@@ -883,11 +885,10 @@ subroutine print_options(model)
          'advective-diffusive balance ',&
          'temp from external file     ' /)
 
-    character(len=*), dimension(0:3), parameter :: flow_law = (/ &
-         'const 1e-16 Pa^-n a^-1      ', &
+    character(len=*), dimension(0:2), parameter :: flow_law = (/ &
+         'uniform factor flwa         ', &
          'Paterson and Budd (T = -5 C)', &
-         'Paterson and Budd           ', &
-         'read flwa/flwastag from file' /)
+         'Paterson and Budd           ' /)
 
     !TODO - Rename slip_coeff to which_btrc?
     character(len=*), dimension(0:5), parameter :: slip_coeff = (/ &
@@ -959,10 +960,11 @@ subroutine print_options(model)
          'artm and d(artm)/dz input as function of (x,y)', &
          'artm input as function of (x,y,z)             ' /)
 
-    character(len=*), dimension(0:2), parameter :: overwrite_acab = (/ &
+    character(len=*), dimension(0:3), parameter :: overwrite_acab = (/ &
          'do not overwrite acab anywhere            ', &
          'overwrite acab where input acab = 0       ', &
-         'overwrite acab where input thck <= minthck' /)
+         'overwrite acab where input thck <= minthck', &
+         'overwrite acab based on input mask        ' /)
 
     ! NOTE: Set gthf = 1 in the config file to read the geothermal heat flux from an input file.
     !       Otherwise it will be overwritten, even if the 'bheatflx' field is present.
@@ -1063,17 +1065,24 @@ subroutine print_options(model)
          'invert for basin-based basal melting parameters            ', &
          'apply basin basal melting parameters from earlier inversion' /)
 
-    character(len=*), dimension(0:2), parameter :: ho_whichbwat = (/ &
+    character(len=*), dimension(0:3), parameter :: ho_whichbwat = (/ &
          'zero basal water depth                          ', &
          'constant basal water depth                      ', &
-         'basal water depth computed from local till model' /)
+         'basal water depth computed from local till model', &
+         'steady-state water routing with flux calculation' /)
 
-    character(len=*), dimension(0:4), parameter :: ho_whicheffecpress = (/ &
+    character(len=*), dimension(0:2), parameter :: ho_flux_routing_scheme = (/ &
+         'D8; route flux to lowest-elevation neighbor      ', &
+         'Dinf; route flux to two lower-elevation neighbors', &
+         'FD8; route flux to all lower-elevation neighbors ' /)
+
+    character(len=*), dimension(0:5), parameter :: ho_whicheffecpress = (/ &
          'full overburden pressure                             ', &
          'reduced effecpress near pressure melting point       ', &
          'reduced effecpress where there is melting at the bed ', &
          'reduced effecpress where bed is connected to ocean   ', &
-         'reduced effecpress with increasing basal water       '/)
+         'reduced effecpress with increasing basal water (B/vP)', &
+         'reduced effecpress with increasing basal water (ramp)'/)
 
     character(len=*), dimension(0:1), parameter :: which_ho_nonlinear = (/ &
          'use standard Picard iteration          ', &
@@ -1234,7 +1243,7 @@ subroutine print_options(model)
     end if
 
     if (tasks > 1 .and. model%options%whichbwat==BWATER_FLUX) then
-       call write_log('Error, flux-based basal water option not supported for more than one processor', GM_FATAL)
+       call write_log('Error, flux-based basal water option not yet supported for more than one processor', GM_FATAL)
     endif
 
     ! Forbidden options associated with Glissade dycore
@@ -1407,6 +1416,11 @@ subroutine print_options(model)
           call write_log(message)
        endif
 
+       if (model%options%smooth_input_usrf) then
+          write(message,*) ' Input usrf will be smoothed'
+          call write_log(message)
+       endif
+
        if (model%options%smooth_input_topography) then
           write(message,*) ' Input topography will be smoothed'
           call write_log(message)
@@ -1766,6 +1780,16 @@ subroutine print_options(model)
           call write_log('Error, HO basal water input out of range', GM_FATAL)
        end if
 
+       if (model%options%which_ho_bwat == HO_BWAT_FLUX_ROUTING) then
+          write(message,*) 'ho_flux_routing_scheme  : ',model%options%ho_flux_routing_scheme,  &
+               ho_flux_routing_scheme(model%options%ho_flux_routing_scheme)
+          call write_log(message)
+          if (model%options%ho_flux_routing_scheme < 0.or. &
+              model%options%ho_flux_routing_scheme >= size(ho_flux_routing_scheme)) then
+             call write_log('Error, HO flux routing scheme out of range', GM_FATAL)
+          end if
+       end if
+
        write(message,*) 'ho_whicheffecpress      : ',model%options%which_ho_effecpress,  &
                          ho_whicheffecpress(model%options%which_ho_effecpress)
        call write_log(message)
@@ -1996,7 +2020,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 +2045,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)
@@ -2033,9 +2058,9 @@ subroutine handle_parameters(section, model)
     call GetValue(section,'pmp_offset',         model%temper%pmp_offset)
     call GetValue(section,'pmp_threshold',      model%temper%pmp_threshold)
     call GetValue(section,'geothermal',         model%paramets%geot)
-    !TODO - Change default_flwa to flwa_constant?  Would have to change config files.
     call GetValue(section,'flow_factor',        model%paramets%flow_enhancement_factor)
     call GetValue(section,'flow_factor_float',  model%paramets%flow_enhancement_factor_float)
+    !TODO - Change default_flwa to flwa_constant?  Would have to change config files.
     call GetValue(section,'default_flwa',       model%paramets%default_flwa)
     call GetValue(section,'efvs_constant',      model%paramets%efvs_constant)
     call GetValue(section,'effstrain_min',      model%paramets%effstrain_min)
@@ -2103,9 +2128,9 @@ subroutine handle_parameters(section, model)
     call GetValue(section, 'effecpress_bmlt_threshold', model%basal_physics%effecpress_bmlt_threshold)
 
     ! basal water parameters
-    call GetValue(section, 'const_bwat', model%basal_physics%const_bwat)
-    call GetValue(section, 'bwat_till_max', model%basal_physics%bwat_till_max)
-    call GetValue(section, 'c_drainage', model%basal_physics%c_drainage)
+    call GetValue(section, 'const_bwat', model%basal_hydro%const_bwat)
+    call GetValue(section, 'bwat_till_max', model%basal_hydro%bwat_till_max)
+    call GetValue(section, 'c_drainage', model%basal_hydro%c_drainage)
 
     ! pseudo-plastic parameters
     !TODO - Put pseudo-plastic and other basal sliding parameters in a separate section
@@ -2206,7 +2231,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 +2396,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)
 
@@ -2627,12 +2655,12 @@ subroutine print_parameters(model)
     endif
 
     if (model%options%which_ho_bwat == HO_BWAT_CONSTANT) then
-       write(message,*) 'constant basal water depth (m): ', model%basal_physics%const_bwat
+       write(message,*) 'constant basal water depth (m): ', model%basal_hydro%const_bwat
        call write_log(message)
     elseif (model%options%which_ho_bwat == HO_BWAT_LOCAL_TILL) then
-       write(message,*) 'maximum till water depth (m)  : ', model%basal_physics%bwat_till_max
+       write(message,*) 'maximum till water depth (m)  : ', model%basal_hydro%bwat_till_max
        call write_log(message)
-       write(message,*) 'till drainage rate (m/yr)     : ', model%basal_physics%c_drainage
+       write(message,*) 'till drainage rate (m/yr)     : ', model%basal_hydro%c_drainage
        call write_log(message)
     endif
 
@@ -3323,7 +3351,7 @@ subroutine define_glide_restart_variables(options)
         ! other Glissade options
 
         ! If overwriting acab in certain grid cells, than overwrite_acab_mask needs to be in the restart file.
-        ! This mask is set at model initialization based on the input acab or ice thickness.
+        ! This mask is read in at model initialization, or is set based on the input acab or ice thickness.
         if (options%overwrite_acab /= 0) then
            call glide_add_to_restart_variable_list('overwrite_acab_mask')
         endif

libglide/glide_temp.F90

@@ -512,7 +512,7 @@ subroutine glide_temp_driver(model,whichtemp)
 
                 call corrpmpt(model%temper%temp(:,ew,ns),     &
                               model%geometry%thck(ew,ns),     &
-                              model%temper%bwat(ew,ns),       &
+                              model%basal_hydro%bwat(ew,ns),  &
                               model%numerics%sigma,           &
                               model%general%upn)
 
@@ -560,7 +560,7 @@ subroutine glide_temp_driver(model,whichtemp)
 
                    call corrpmpt(model%temper%temp(:,ew,ns),     &
                                  model%geometry%thck(ew,ns),     &
-                                 model%temper%bwat(ew,ns),       &
+                                 model%basal_hydro%bwat(ew,ns),  &
                                  model%numerics%sigma,           &
                                  model%general%upn)