Skip to content

Commit

Permalink
Basal hydrology changes to support N. Hemisphere ice sheets
Browse files Browse the repository at this point in the history 
This commit contains several minor changes to support the flux-routing hydrology scheme
in runs with Northern Hemisphere paleo ice sheets

* Added a new effective pressure option: which_ho_effecpress = HO_EFFECPRESS_BWAT_RAMP = 5.
  This is similar to the existing option HO_EFFECPRESS_BWAT = 4.
  The only difference is that as bwat increases from 0 to bwat_till_max,
  the effective pressure ramps down linearly, unlike the more complex formulation
  of Bueler and van Pelt (2015).

* Added a new logical option, smooth_input_usrf.  When this option is set to true,
  the initial upper surface elevation field (usrf) is smoothed, and the thickness
  is then adjusted to be consistent.  This is helpful for a 4-km N. Hemisphere simulation
  using the input file cism_USGS-huy3-S1D_4km.nc.
  In this file, parts of the GrIS have rough topography and fairly smooth thck,
  leading to rough usrf and large surface gradients that crash the solver in the
  first few time steps.  With several smoothing passes, the code starts cleanly.

  The smoothing is done in subroutine glissade_smooth_usrf, in glissade_utils.F90.
  Note that usrf is smoothed only for grounded ice, and not for floating ice,
  ice-free ocean, or ice-free land.

* Made the fill_depression routine more efficient by restructuring with an outer and inner loop.
  Halo updates are called only from the outer loop.
  Although more efficient than before, it is still very expensive to fill depressions with the
  current algorithm.  The cost of the entire code more than doubles on a 4-km N. Hemisphere grid.
  In a future commit, I will try to implement a different algorithm that scales better.

* Added the max value of bmlt_ground in the diagnostic log file
  • Loading branch information
@whlipscomb
whlipscomb committed Aug 5, 2021
1 parent 929ca40 commit be4e749
Show file tree
Hide file tree
Showing 10 changed files with 457 additions and 177 deletions.
46 changes: 39 additions & 7 deletions libglide/glide_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)

Expand All @@ -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)

Expand Down
11 changes: 9 additions & 2 deletions libglide/glide_setup.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -1075,12 +1076,13 @@ subroutine print_options(model)
'Dinf; route flux to two lower-elevation neighbors', &
'FD8; route flux to all lower-elevation neighbors ' /)

character(len=*), dimension(0:4), parameter :: ho_whicheffecpress = (/ &
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 ', &
Expand Down Expand Up @@ -1414,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)
Expand Down
7 changes: 6 additions & 1 deletion libglide/glide_types.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,7 @@ module glide_types
integer, parameter :: HO_EFFECPRESS_BMLT = 2
integer, parameter :: HO_EFFECPRESS_OCEAN_PENETRATION = 3
integer, parameter :: HO_EFFECPRESS_BWAT = 4
integer, parameter :: HO_EFFECPRESS_BWAT_RAMP = 5

!WHL - added Picard acceleration option
integer, parameter :: HO_NONLIN_PICARD = 0
Expand Down Expand Up @@ -675,6 +676,9 @@ module glide_types
logical :: adjust_input_thickness = .false.
!> if true, then adjust thck to maintain usrf, instead of deriving usrf from topg and thck

logical :: smooth_input_usrf = .false.
!> if true, then apply Laplacian smoothing to usrf at initialization

logical :: smooth_input_topography = .false.
!> if true, then apply Laplacian smoothing to the topography at initialization

Expand Down Expand Up @@ -862,7 +866,8 @@ module glide_types
!> \item[1] N is reduced where the bed is at or near the pressure melting point
!> \item[2] N is reduced where there is melting at the bed
!> \item[3] N is reduced due to connection of subglacial water to the ocean
!> \item[4] N is reduced where basal water is present
!> \item[4] N is reduced where basal water is present, following Bueler/van Pelt
!> \item[5] N is reduced where basal water is present, with a ramp function
!> \end{description}

integer :: which_ho_nonlinear = 0
Expand Down
110 changes: 49 additions & 61 deletions libglissade/glissade.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,8 @@ subroutine glissade_initialise(model, evolve_ice)
use glissade_inversion, only: glissade_init_inversion, verbose_inversion
use glissade_bmlt_float, only: glissade_bmlt_float_thermal_forcing_init, verbose_bmlt_float
use glissade_grounding_line, only: glissade_grounded_fraction
use glissade_utils, only: &
glissade_adjust_thickness, glissade_smooth_topography, glissade_adjust_topography
use glissade_utils, only: glissade_adjust_thickness, glissade_smooth_usrf, &
glissade_smooth_topography, glissade_adjust_topography
use glissade_utils, only: glissade_stdev
use felix_dycore_interface, only: felix_velo_init

Expand Down Expand Up @@ -396,7 +396,15 @@ subroutine glissade_initialise(model, evolve_ice)
call glissade_adjust_thickness(model)
endif

! Optionally, smooth the input topography with a 9-point Laplacian smoother.
! Optionally, smooth the input surface elevation with a Laplacian smoother.
! This subroutine does not change the topg, but returns thck consistent with the new usrf.
! If the initial usrf is rough, then multiple smoothing passes may be needed to stabilize the flow.

if (model%options%smooth_input_usrf .and. model%options%is_restart == RESTART_FALSE) then
call glissade_smooth_usrf(model, nsmooth = 5)
endif ! smooth_input_usrf

! Optionally, smooth the input topography with a Laplacian smoother.

if (model%options%smooth_input_topography .and. model%options%is_restart == RESTART_FALSE) then
call glissade_smooth_topography(model)
Expand Down Expand Up @@ -625,17 +633,17 @@ subroutine glissade_initialise(model, evolve_ice)
if (make_ice_domain_mask) then

where (model%geometry%thck > 0.0d0 .or. model%geometry%topg > 0.0d0)
!! where (model%geometry%thck > 0.0d0) ! uncomment for terrestrial margins
model%general%ice_domain_mask = 1
elsewhere
model%general%ice_domain_mask = 0
endwhere

! Extend the mask a couple of cells in each direction to be on the safe side.
! Extend the mask a few cells in each direction to be on the safe side.
! The number of buffer layers could be made a config parameter.

allocate(ice_domain_mask(model%general%ewn,model%general%nsn))

!! do k = 1, 2
do k = 1, 3
call parallel_halo(model%general%ice_domain_mask, parallel)
ice_domain_mask = model%general%ice_domain_mask ! temporary copy
Expand Down Expand Up @@ -1924,24 +1932,23 @@ subroutine glissade_thermal_solve(model, dt)
if (model%options%which_ho_bwat == HO_BWAT_FLUX_ROUTING) then

!WHL - Temporary code for debugging: Make up a simple basal melt field.
model%basal_hydro%head(:,:) = &
model%geometry%thck(:,:)*thk0 + (rhow/rhoi)*model%geometry%topg(:,:)*thk0
head_max = maxval(model%basal_hydro%head) ! max on local processor
head_max = parallel_reduce_max(head_max) ! global max
do j = 1, model%general%nsn
do i = 1, model%general%ewn
if (head_max - model%basal_hydro%head(i,j) < 1000.d0) then
! model%basal_hydro%head(:,:) = &
! model%geometry%thck(:,:)*thk0 + (rhow/rhoi)*model%geometry%topg(:,:)*thk0
! head_max = maxval(model%basal_hydro%head) ! max on local processor
! head_max = parallel_reduce_max(head_max) ! global max
! do j = 1, model%general%nsn
! do i = 1, model%general%ewn
! if (head_max - model%basal_hydro%head(i,j) < 1000.d0) then
!! if (head_max - model%basal_hydro%head(i,j) < 200.d0) then
bmlt_ground_unscaled(i,j) = 1.0d0/scyr ! units are m/s
else
bmlt_ground_unscaled(i,j) = 0.0d0
endif
enddo
enddo
! bmlt_ground_unscaled(i,j) = 1.0d0/scyr ! units are m/s
! else
! bmlt_ground_unscaled(i,j) = 0.0d0
! endif
! enddo
! enddo

! Compute some masks needed below


call glissade_get_masks(&
model%general%ewn, model%general%nsn, &
model%parallel, &
Expand Down Expand Up @@ -1980,34 +1987,6 @@ subroutine glissade_thermal_solve(model, dt)
endwhere
endif

!WHL - debug
print*, ' '
print*, 'edge_mask:'
write(6,'(a6)',advance='no') ' '
do i = itest-5, itest+5
write(6,'(i5)',advance='no') i
enddo
write(6,*) ' '
do j = jtest+5, jtest-5, -1
write(6,'(i6)',advance='no') j
do i = itest-5, itest+5
write(6,'(i5)',advance='no') model%general%global_edge_mask(i,j)
enddo
write(6,*) ' '
enddo
print*, ' '
print*, ' '
print*, 'overwrite_acab_mask:'
write(6,*) ' '
do j = jtest+5, jtest-5, -1
write(6,'(i6)',advance='no') j
do i = itest-5, itest+5
write(6,'(i5)',advance='no') model%climate%overwrite_acab_mask(i,j)
enddo
write(6,*) ' '
enddo
print*, ' '

call parallel_halo(bwat_mask, parallel)

! Compute bwat based on a steady-state flux routing scheme
Expand Down Expand Up @@ -3202,6 +3181,14 @@ subroutine glissade_calving_solve(model, init_calving)
! Note: Currently hardwired to include 13 of the 16 ISMIP6 basins.
! Does not include the three largest shelves (Ross, Filchner-Ronne, Amery)

call glissade_get_masks(nx, ny, &
parallel, &
model%geometry%thck*thk0, model%geometry%topg*thk0, &
model%climate%eus*thk0, 0.0d0, & ! thklim = 0
ice_mask, &
floating_mask = floating_mask, &
land_mask = land_mask)

if (init_calving .and. model%options%expand_calving_mask) then

! Identify basins whose floating ice will be added to the calving mask
Expand All @@ -3219,14 +3206,6 @@ subroutine glissade_calving_solve(model, init_calving)
enddo
endif

call glissade_get_masks(nx, ny, &
parallel, &
model%geometry%thck*thk0, model%geometry%topg*thk0, &
model%climate%eus*thk0, 0.0d0, & ! thklim = 0
ice_mask, &
floating_mask = floating_mask, &
land_mask = land_mask)

if (verbose_calving .and. this_rank==rtest) then
print*, ' '
print*, 'initial calving_mask, itest, jtest, rank =', itest, jtest, rtest
Expand Down Expand Up @@ -3623,6 +3602,12 @@ subroutine glissade_calving_solve(model, init_calving)
! halo updates
call parallel_halo(model%geometry%thck, parallel) ! Updated halo values of thck are needed below in calclsrf

! update the upper and lower surfaces

call glide_calclsrf(model%geometry%thck, model%geometry%topg, &
model%climate%eus, model%geometry%lsrf)
model%geometry%usrf(:,:) = max(0.d0, model%geometry%thck(:,:) + model%geometry%lsrf(:,:))

if (verbose_calving .and. this_rank == rtest) then
print*, ' '
print*, 'Final calving_thck (m), itest, jtest, rank =', itest, jtest, rtest
Expand All @@ -3642,14 +3627,17 @@ subroutine glissade_calving_solve(model, init_calving)
enddo
write(6,*) ' '
enddo
print*, ' '
print*, 'Final usrf (m):'
do j = jtest+3, jtest-3, -1
write(6,'(i6)',advance='no') j
do i = itest-3, itest+3
write(6,'(f10.3)',advance='no') model%geometry%usrf(i,j) * thk0
enddo
write(6,*) ' '
enddo
endif ! verbose_calving

! update the upper and lower surfaces

call glide_calclsrf(model%geometry%thck, model%geometry%topg, &
model%climate%eus, model%geometry%lsrf)
model%geometry%usrf(:,:) = max(0.d0, model%geometry%thck(:,:) + model%geometry%lsrf(:,:))

end subroutine glissade_calving_solve

!=======================================================================
Expand Down
38 changes: 32 additions & 6 deletions libglissade/glissade_basal_traction.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -850,7 +850,7 @@ subroutine calc_effective_pressure (which_effecpress, &

if (present(bwat)) then

! Reduce N where basal water is present.
! Reduce N where basal water is present, following Bueler % van Pelt (2015).
! The effective pressure decreases from overburden P_0 for bwat = 0 to a small value for bwat = bwat_till_max.
! Note: Instead of using a linear ramp for the variation between overburden and the small value
! (as for the BPMP and BMLT options above), we use the published formulation of Bueler & van Pelt (2015).
Expand All @@ -876,11 +876,6 @@ subroutine calc_effective_pressure (which_effecpress, &
!! basal_physics%effecpress(i,j) = basal_physics%effecpress_delta * overburden(i,j) &
!! * 10.d0**((basal_hydro%e_0/basal_hydro%C_c) * (1.0d0 - relative_bwat))

!WHL - Uncomment to try a linear ramp in place of the Bueler & van Pelt relationship.
! This might lead to smoother variations in N with spatial variation in bwat.
!! basal_physics%effecpress(i,j) = overburden(i,j) * &
!! (basal_physics%effecpress_delta + (1.0d0 - relative_bwat) * (1.0d0 - basal_physics%effecpress_delta))

! limit so as not to exceed overburden
basal_physics%effecpress(i,j) = min(basal_physics%effecpress(i,j), overburden(i,j))
end if
Expand All @@ -894,6 +889,37 @@ subroutine calc_effective_pressure (which_effecpress, &
basal_physics%effecpress = 0.0d0
end where

case(HO_EFFECPRESS_BWAT_RAMP) ! Similar to HO_EFFECPRESS_BWAT, but with a ramp function

! Initialize for the case where bwat isn't present, and also for points with bwat == 0

basal_physics%effecpress(:,:) = overburden(:,:)

if (present(bwat)) then

! Reduce N where basal water is present.
! The effective pressure decreases from overburden P_0 for bwat = 0 to a small value for bwat = bwat_till_max.

do j = 1, nsn
do i = 1, ewn
if (bwat(i,j) > 0.0d0) then

relative_bwat = max(0.0d0, min(bwat(i,j)/basal_hydro%bwat_till_max, 1.0d0))

basal_physics%effecpress(i,j) = overburden(i,j) * &
(basal_physics%effecpress_delta + (1.0d0 - relative_bwat) * (1.0d0 - basal_physics%effecpress_delta))

end if
enddo
enddo

endif ! present(bwat)

where (floating_mask == 1)
! set to zero for floating ice
basal_physics%effecpress = 0.0d0
end where

case(HO_EFFECPRESS_OCEAN_PENETRATION)

! Reduce N for ice grounded below sea level based on connectivity of subglacial water to the ocean
Expand Down
Loading

0 comments on commit be4e749

Please sign in to comment.