From d9ab6130158d1405bbe25a2f9487668f83190ae2 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Wed, 20 Feb 2019 10:34:12 -0500 Subject: [PATCH 1/9] Implementing UM damage evolution This commit adds the long-wavelength, dimensional damage evolution equation from Bassis and Ma 2015 to glissade_calving under the marine_margin=CALVING_DAMAGE option. Currently, this form of damage is a direct replacement of the simple effective stress scheme that existed there previously, but a later commit will add i/o options to allow for toggling between the two. The damage equation includes a term for the total melt rate, which I have implemented using -acab. I suspect that this will need to be updated to accommodate for the current state of basal melt parameterizations in the model, but I am not up-to-date on that and have coded what I know to exist. --- libglissade/glissade.F90 | 4 +++- libglissade/glissade_calving.F90 | 40 ++++++++++++++++++++++++++++---- 2 files changed, 39 insertions(+), 5 deletions(-) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 8c9f15e5..02877b3d 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -2007,7 +2007,7 @@ subroutine glissade_calving_solve(model, init_calving) use parallel - use glimmer_paramets, only: thk0, tim0, len0 + use glimmer_paramets, only: thk0, tim0, len0, evs0 use glissade_calving, only: glissade_calve_ice use glide_mask, only: glide_set_mask @@ -2057,6 +2057,8 @@ subroutine glissade_calving_solve(model, init_calving) model%numerics%dns*len0, & ! m model%numerics%sigma, & model%numerics%thklim*thk0, & ! m + model%stress%efvs*evs0, & ! Pa s + model%climate%acab*thk0/tim0, & ! m/s thck_unscaled, & ! m model%isostasy%relx*thk0, & ! m model%geometry%topg*thk0, & ! m diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index fa7c89c5..892b6a18 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -204,6 +204,8 @@ subroutine glissade_calve_ice(which_calving, & dx, dy, & ! m sigma, & thklim, & ! m + efvs, & ! Pa s + acab, & ! m/s thck, relx, & ! m topg, eus) ! m @@ -265,6 +267,8 @@ subroutine glissade_calve_ice(which_calving, & real(dp), intent(in) :: dx, dy !> grid cell size in x and y directions (m) real(dp), dimension(:), intent(in) :: sigma !> vertical sigma coordinate real(dp), intent(in) :: thklim !> minimum thickness for dynamically active grounded ice (m) + real(dp), dimension(:,:,:), intent(in) :: efvs !> effective viscosity (Pa s) + real(dp), dimension(:,:), intent(in) :: acab !> mass balance (m/s) real(dp), dimension(:,:), intent(inout) :: thck !> ice thickness (m) real(dp), dimension(:,:), intent(in) :: relx !> relaxed bedrock topography (m) real(dp), dimension(:,:), intent(in) :: topg !> present bedrock topography (m) @@ -278,6 +282,10 @@ subroutine glissade_calve_ice(which_calving, & integer :: i, j, k, n integer :: ii, jj + real(dp), dimension(:,:,:), allocatable :: & + eps_max, & ! maximum principal strain rate + d_damage_dt ! rate of change of damage scalar (1/s) + real(dp), dimension(:,:), allocatable :: & thck_calving_front, & ! effective ice thickness at the calving front thck_init, & ! value of thck before calving @@ -318,9 +326,11 @@ subroutine glissade_calve_ice(which_calving, & frac_lateral, & ! lateral_rate / lateral_rate_max areafrac, & ! fractional ice-covered area in a calving_front cell dthck, & ! thickness change (m) - d_damage_dt, & ! rate of change of damage scalar (1/s) thckmax_cliff, & ! max stable ice thickness in marine_cliff cells - factor ! factor in quadratic formula + factor, & ! factor in quadratic formula + beta, & ! flow law exponent anisotropy parameter + nstar, & ! effective flow law exponent + szero ! ratio of hydrostatic to tensile stress real(dp), parameter :: & thinning_limit = 0.99d0 ! When ice not originally on the calving front is allowed to thin, @@ -380,6 +390,7 @@ subroutine glissade_calve_ice(which_calving, & allocate (thck_calving_front(nx,ny)) allocate (thck_init(nx,ny)) allocate (marine_cliff_mask(nx,ny)) + allocate (d_damage_dt(nz,nx,ny)) !WHL - debug allocate(calving_thck_init(nx,ny)) @@ -591,11 +602,30 @@ subroutine glissade_calve_ice(which_calving, & ! Note: The damage is formally a 3D field, which makes it easier to advect, even though ! (in the current scheme) the damage source term is uniform in each column. + ! Allocate strain rate eigenvalue matrix + allocate(eps_max(nz,nx,ny)) + do j = 2, ny-1 do i = 2, nx-1 if (floating_mask(i,j) == 1) then - d_damage_dt = calving%damage_constant * calving%tau_eff(i,j) ! damage_constant has units of s^{-1}/(Pa) - calving%damage(:,i,j) = calving%damage(:,i,j) + d_damage_dt * dt + ! Compute the ratio of the principal stresses (equivalent to the ratio of the + ! principal strain rates) + beta = calving%tau_eigen2(i,j) / calving%tau_eigen1(i,j) + + ! Find the max principal strain rate using Glen's Flow Law + eps_max(:,i,j) = calving%tau_eigen1(i,j) / (2.0d0 * efvs(:,i,j)) + + ! Compute the effective flow law exponent + nstar = 12.0d0*(1.0d0+beta+beta**2) / (4.0d0*(1.0d0+beta+beta**2)+6.0d0*beta**2) + + ! Compute the hydrostatic to tensile stress ratio + szero = rhoi*(rhoo-rhoi)*grav*thck(i,j)/(2.0d0*calving%tau_eigen1(i,j)*rhoo) + + ! Prognose damage following Bassis & Ma 2015 + d_damage_dt(:,i,j) = (nstar*(1.0d0-szero)*eps_max(:,i,j) & + - acab(i,j)/thck(i,j)) * calving%damage(:,i,j) +! d_damage_dt = calving%damage_constant * calving%tau_eff(i,j) ! damage_constant has units of s^{-1}/(Pa) + calving%damage(:,i,j) = calving%damage(:,i,j) + d_damage_dt(:,i,j) * dt calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) else ! set damage to zero for grounded ice @@ -1375,11 +1405,13 @@ subroutine glissade_calve_ice(which_calving, & deallocate (thck_calving_front) deallocate (thck_init) deallocate (marine_cliff_mask) + deallocate (d_damage_dt) if (allocated(calving_thck_init)) deallocate(calving_thck_init) if (allocated(damage_column)) deallocate(damage_column) if (allocated(tau1)) deallocate(tau1) if (allocated(tau2)) deallocate(tau2) + if (allocated(eps_max)) deallocate(eps_max) end subroutine glissade_calve_ice From 6b3c8c8c024aa909b87ddf89cd131d58c41df838 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Fri, 5 Apr 2019 12:08:03 -0400 Subject: [PATCH 2/9] Damage input options, new diagnostic variables, and calving This commit adds several new input options associated with damage-based calving: * damage_src - Integer specifying the source term used to solve for damage - 0 = no damage source, 1 (default) = damage source based on the effective stress (what was previously already implemented), 2 = damage source from eq 27 in Bassis and Ma 2015 (what has been implemented here) * damage_floor - Integer used to set the lower-limit value to which damage is constrained - 0 (default) = zero damage floor, 1 = Nye zero-stress damage floor (currently due to only basal crevasses) * damage_advect - Boolean used to turn damage advection on (default) or off * damage_manufactured - Boolean used to toggle a manufactured solution for damage (default is off) - This is entirely a diagnostic option used for verification of the damage field in an ice tongue test case; the specific details of how this works will be added in a future commit, but in essence, when this is toggled on, damage is modified by a forcing term that allows the user to compare it to a known solution. In addition to the above input options, two new diagnostic variables have been added: * const_thk_mask - Similar to no_advance_mask, but it holds the thickness constant in time in the specified grid cells - This is used in the ice tongue test case (to be committed later) to enforce a constant grounding line flux. * damage_init - Holds the initial damage field specified (or not specified) in the input file - When damage_manufactured is true, we need this information to compute the forcing term. There may be a more elegant way to access this information in glissade_calving that I do not know about, in which case this variable would be unnecessary. Finally, this commit also adds damage calving support using the already-existing lateral calving rate conversion and subgrid calving front parameterization scheme. When the latter is active, damage leads to calving, and when it is inactive, damage does not calve (treating it entirely as a scalar field layered on top of the ice). --- libglide/glide_setup.F90 | 67 +++++- libglide/glide_types.F90 | 42 ++++ libglide/glide_vars.def | 17 ++ libglissade/glissade.F90 | 63 ++++-- libglissade/glissade_calving.F90 | 344 +++++++++++++++++++++-------- libglissade/glissade_transport.F90 | 6 +- 6 files changed, 418 insertions(+), 121 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 7acbfd24..06464e6d 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -570,9 +570,13 @@ subroutine handle_options(section, model) call GetValue(section,'marine_margin',model%options%whichcalving) call GetValue(section,'calving_init',model%options%calving_init) call GetValue(section,'calving_domain',model%options%calving_domain) + call GetValue(section,'damage_src', model%options%damage_src) + call GetValue(section,'damage_floor', model%options%damage_floor) call GetValue(section,'remove_icebergs', model%options%remove_icebergs) call GetValue(section,'limit_marine_cliffs', model%options%limit_marine_cliffs) call GetValue(section,'cull_calving_front', model%options%cull_calving_front) + call GetValue(section,'damage_advect', model%options%damage_advect) + call GetValue(section,'damage_manufactured', model%options%damage_manufactured) call GetValue(section,'dm_dt_diag',model%options%dm_dt_diag) call GetValue(section,'diag_minthck',model%options%diag_minthck) call GetValue(section,'vertical_integration',model%options%whichwvel) @@ -778,6 +782,15 @@ subroutine print_options(model) 'calving only at the ocean edge ', & 'calving in all cells where criterion is met'/) + character(len=*), dimension(0:2), parameter :: damage_src = (/ & + 'no damage source ', & + 'effective stress damage source ', & + 'bassis and ma damage source ' /) + + character(len=*), dimension(0:1), parameter :: damage_floor = (/ & + 'zero damage floor ', & + 'nye damage floor ' /) + character(len=*), dimension(0:1), parameter :: dm_dt_diag = (/ & 'write dmass/dt diagnostic in units of kg/s ', & 'write dmass/dt diagnostic in units of Gt/yr'/) @@ -1104,6 +1117,27 @@ subroutine print_options(model) write(message,*) 'calving_domain : ', model%options%calving_domain, domain_calving(model%options%calving_domain) call write_log(message) + if (model%options%damage_src < 0 .or. model%options%damage_src >= size(damage_src)) then + call write_log('Error, damage_src out of range',GM_FATAL) + elseif (model%options%damage_src /= EFF_STRESS_DAMAGE_SRC & + .and. model%options%whichcalving /= CALVING_DAMAGE) then + write(message,*) 'WARNING: damage source type can be selected for damage marine margin only; user selection ignored' + call write_log(message, GM_WARNING) + model%options%damage_src = EFF_STRESS_DAMAGE_SRC + end if + write(message,*) 'damage_src : ', model%options%damage_src, damage_src(model%options%damage_src) + call write_log(message) + + if (model%options%damage_floor < 0 .or. model%options%damage_floor >= size(damage_floor)) then + call write_log('Error, damage_floor out of range',GM_FATAL) + elseif (model%options%damage_floor /= ZERO_DAMAGE_FLOOR .and. model%options%whichcalving /= CALVING_DAMAGE) then + write(message,*) 'WARNING: damage floor type can be selected for damage marine margin only; user selection ignored' + call write_log(message, GM_WARNING) + model%options%damage_floor = ZERO_DAMAGE_FLOOR + end if + write(message,*) 'damage_floor : ', model%options%damage_floor, damage_floor(model%options%damage_floor) + call write_log(message) + ! dycore-dependent calving options if (model%options%whichdycore == DYCORE_GLISSADE) then @@ -1128,6 +1162,18 @@ subroutine print_options(model) call write_log(' Calving-front cells will not be culled at initialization') endif + if (model%options%damage_advect) then + call write_log(' Damage will advect') + else + call write_log(' Damage will not advect') + endif + + if (model%options%damage_manufactured) then + call write_log(' Damage will use a manufactured solution') + else + call write_log(' Damage will evolve normally') + endif + if (model%options%whichcalving == CALVING_FLOAT_FRACTION) then write(message,*) 'WARNING: calving float fraction option deprecated with Glissade_dycore; set calving_timescale instead' call write_log(message, GM_WARNING) @@ -1144,7 +1190,9 @@ subroutine print_options(model) if (model%options%whichcalving == CALVING_GRID_MASK) then call write_log('Error, calving grid mask option is supported for Glissade dycore only', GM_FATAL) endif - if (model%options%whichcalving == CALVING_DAMAGE) then + if (model%options%whichcalving == CALVING_DAMAGE & + .or. model%options%damage_src /= EFF_STRESS_DAMAGE_SRC & + .or. model%options%damage_floor /= ZERO_DAMAGE_FLOOR) then call write_log('Error, calving damage option is supported for Glissade dycore only', GM_FATAL) endif if (model%options%calving_domain /= CALVING_DOMAIN_OCEAN_EDGE) then @@ -1818,6 +1866,23 @@ subroutine print_parameters(model) write(message,*) 'Setting remove_icebergs = T for stability when using the calving_front subgrid scheme' call write_log(message) endif + + if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_manufactured) then + model%options%which_ho_calving_front = HO_CALVING_FRONT_NO_SUBGRID + write(message,*) 'Setting which_ho_calving_front =', HO_CALVING_FRONT_NO_SUBGRID + call write_log(message) + write(message,*) 'Manufactured damage calving option does not support subgrid calving front scheme' + call write_log(message) + elseif (model%options%whichcalving == CALVING_DAMAGE) then + write(message,*) 'Subgrid calving front scheme chosen: damage will be converted into a lateral calving rate' + call write_log(message) + endif + endif + + if (model%options%which_ho_calving_front == HO_CALVING_FRONT_NO_SUBGRID & + .and. model%options%whichcalving == CALVING_DAMAGE) then + write(message,*) 'Subgrid calving front scheme not chosen: ice will not calve based on damage' + call write_log(message) endif if (model%options%limit_marine_cliffs) then diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index ec6013e2..b9ee4413 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -170,6 +170,13 @@ module glide_types integer, parameter :: CALVING_DOMAIN_OCEAN_EDGE = 0 integer, parameter :: CALVING_DOMAIN_EVERYWHERE = 1 + integer, parameter :: NO_DAMAGE_SRC = 0 + integer, parameter :: EFF_STRESS_DAMAGE_SRC = 1 + integer, parameter :: BASSIS_MA_DAMAGE_SRC = 2 + + integer, parameter :: ZERO_DAMAGE_FLOOR = 0 + integer, parameter :: NYE_DAMAGE_FLOOR = 1 + integer, parameter :: VERTINT_STANDARD = 0 integer, parameter :: VERTINT_KINEMATIC_BC = 1 @@ -509,6 +516,21 @@ module glide_types !> \item[1] Calve wherever the calving criterion is met !> \end{description} + integer :: damage_src = 1 + !> \begin{description} + !> \item[0] No damage source + !> \item[1] Effective stress damage source; damage prognosed using a simple + !> scheme based on the effective tensile stress + !> \item[2] Bassis and Ma damage source; damage prognosed according to eq 27 + !> in Bassis and Ma 2015 + !> \end{description} + + integer :: damage_floor = 0 + !> \begin{description} + !> \item[0] Zero damage floor; damage cannot be reduced below zero + !> \item[1] Nye damage floor; damage cannot be reduced below the Nye zero stress value + !> \end{description} + logical :: remove_icebergs = .true. !> if true, then identify and remove icebergs after calving !> These are connected regions with zero basal traction and no connection to grounded ice. @@ -521,6 +543,12 @@ module glide_types !> if true, then cull calving_front cells at initialization !> This can make the run more stable by removing long, thin peninsulas + logical :: damage_advect = .true. + !> if true, then damage advects as a scalar tracer + + logical :: damage_manufactured = .false. + !> if true, then damage is forced to produce a manufactured solution + integer :: dm_dt_diag = 0 !> \begin{description} !> \item[0] Write dmass/dt diagnostic in units of kg/s @@ -897,6 +925,8 @@ module glide_types integer, dimension(:,:),pointer :: stagmask => null() !> see glide_mask.f90 for possible values + integer, dimension(:,:),pointer :: const_thk_mask => null() + ! mass fluxes at upper, lower and lateral boundaries ! TODO: Move to a flux derived type? ! Note: sfc_mbal_flux and basal_mbal_flux are not strictly needed, since they are equal to acab_applied and bmlt_applied @@ -1132,6 +1162,7 @@ module glide_types real(dp),dimension(:,:), pointer :: tau_eigen2 => null() !> second eigenvalue of 2D horizontal stress tensor (Pa) real(dp),dimension(:,:), pointer :: tau_eff => null() !> effective stress (Pa) for calving; derived from tau_eigen1, tau_eigen2 real(dp),dimension(:,:,:),pointer :: damage => null() !> 3D damage tracer, 0 > damage < 1 (whichcalving = CALVING_DAMAGE) + real(dp),dimension(:,:,:),pointer :: damage_init => null() !> initial damage field (whichcalving = CALVING_DAMAGE) real(dp) :: marine_limit = -200.d0 !> value of topg/relx at which floating ice calves (m) !> (whichcalving = CALVING_RELX_THRESHOLD, CALVING_TOPG_THRESHOLD) @@ -2185,6 +2216,7 @@ subroutine glide_allocarr(model) call coordsystem_allocate(model%general%ice_grid, model%geometry%thkmask) call coordsystem_allocate(model%general%velo_grid, model%geometry%stagmask) call coordsystem_allocate(model%general%ice_grid, model%geometry%cell_area) + call coordsystem_allocate(model%general%ice_grid, model%geometry%const_thk_mask) call coordsystem_allocate(model%general%velo_grid, model%geomderv%stagthck) call coordsystem_allocate(model%general%velo_grid, model%geomderv%dthckdew) @@ -2327,9 +2359,15 @@ subroutine glide_allocarr(model) call coordsystem_allocate(model%general%ice_grid, model%calving%tau_eff) if (model%options%whichcalving == CALVING_DAMAGE) then call coordsystem_allocate(model%general%ice_grid, upn-1, model%calving%damage) + if (model%options%damage_manufactured) then + call coordsystem_allocate(model%general%ice_grid, upn-1, model%calving%damage_init) + else + allocate(model%calving%damage_init(1,1,1)) + endif else ! allocate with size 1, since they need to be allocated to be passed to calving subroutine allocate(model%calving%damage(1,1,1)) + allocate(model%calving%damage_init(1,1,1)) endif ! matrix solver arrays @@ -2684,6 +2722,8 @@ subroutine glide_deallocarr(model) deallocate(model%geometry%thkmask) if (associated(model%geometry%stagmask)) & deallocate(model%geometry%stagmask) + if (associated(model%geometry%const_thk_mask)) & + deallocate(model%geometry%const_thk_mask) if (associated(model%geomderv%stagthck)) & deallocate(model%geomderv%stagthck) if (associated(model%geomderv%dthckdew)) & @@ -2818,6 +2858,8 @@ subroutine glide_deallocarr(model) deallocate(model%calving%tau_eff) if (associated(model%calving%damage)) & deallocate(model%calving%damage) + if (associated(model%calving%damage_init)) & + deallocate(model%calving%damage_init) ! matrix solver arrays diff --git a/libglide/glide_vars.def b/libglide/glide_vars.def index 57554459..49d6dda5 100644 --- a/libglide/glide_vars.def +++ b/libglide/glide_vars.def @@ -573,6 +573,14 @@ data: data%calving%damage(up,:,:) load: 1 coordinates: lon lat +[damage_init] +dimensions: time, staglevel, y1, x1 +units: unitless [0,1] +long_name: initial ice damage +data: data%calving%damage_init(up,:,:) +load: 1 +coordinates: lon lat + [area_factor] dimensions: time, y1, x1 units: unitless @@ -643,6 +651,15 @@ type: int coordinates: lon lat load: 1 +[const_thk_mask] +dimensions: time, y1, x1 +long_name: mask +units: 1 +data: data%geometry%const_thk_mask +type: int +coordinates: lon lat +load: 1 + [usurf] dimensions: time, y1, x1 units: meter diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 02877b3d..499d050a 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -510,6 +510,11 @@ subroutine glissade_initialise(model, evolve_ice) itest = model%numerics%idiag_local jtest = model%numerics%jdiag_local + ! Save the initial damage field + if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_manufactured) then + model%calving%damage_init(:,:,:) = model%calving%damage(:,:,:) + endif + ! initial calving, if desired ! Note: Do this only for a cold start with evolving ice, not for a restart if (l_evolve_ice .and. & @@ -1683,6 +1688,15 @@ subroutine glissade_thickness_tracer_solve(model) model%geometry%thck(:,:) = model%geometry%thck_old(:,:) endif + do i = 1, model%general%ewn + do j = 1, model%general%nsn + ! Enforce constant thickness where specified + if (model%geometry%const_thk_mask(i,j) == 1) then + model%geometry%thck(i,j) = model%geometry%thck_old(i,j) + endif + enddo + enddo + end select !------------------------------------------------------------------------ @@ -2007,7 +2021,8 @@ subroutine glissade_calving_solve(model, init_calving) use parallel - use glimmer_paramets, only: thk0, tim0, len0, evs0 + use glimmer_physcon, only: scyr + use glimmer_paramets, only: thk0, vel0, tim0, len0, evs0 use glissade_calving, only: glissade_calve_ice use glide_mask, only: glide_set_mask @@ -2042,27 +2057,33 @@ subroutine glissade_calving_solve(model, init_calving) thck_unscaled(:,:) = model%geometry%thck(:,:)*thk0 - call glissade_calve_ice(model%options%whichcalving, & - model%options%calving_domain, & + call glissade_calve_ice(model%options%whichcalving, & + model%options%calving_domain, & + model%options%damage_src, & + model%options%damage_floor, & model%options%which_ho_calving_front, & - model%options%remove_icebergs, & - model%options%limit_marine_cliffs, & - cull_calving_front, & - model%calving, & ! calving object; includes calving_thck (m) - model%numerics%idiag_local, & - model%numerics%jdiag_local, & - model%numerics%rdiag_local, & - model%numerics%dt*tim0, & ! s - model%numerics%dew*len0, & ! m - model%numerics%dns*len0, & ! m - model%numerics%sigma, & - model%numerics%thklim*thk0, & ! m - model%stress%efvs*evs0, & ! Pa s - model%climate%acab*thk0/tim0, & ! m/s - thck_unscaled, & ! m - model%isostasy%relx*thk0, & ! m - model%geometry%topg*thk0, & ! m - model%climate%eus*thk0) ! m + model%options%remove_icebergs, & + model%options%limit_marine_cliffs, & + cull_calving_front, & + model%options%damage_manufactured, & + model%options%damage_advect, & + model%calving, & ! calving object; includes calving_thck (m) + model%numerics%idiag_local, & + model%numerics%jdiag_local, & + model%numerics%rdiag_local, & + model%numerics%time*scyr, & ! s + model%numerics%dt*tim0, & ! s + model%numerics%dew*len0, & ! m + model%numerics%dns*len0, & ! m + model%numerics%sigma, & + model%numerics%thklim*thk0, & ! m + model%stress%efvs*evs0, & ! Pa s + model%velocity%uvel*vel0, & ! m/s + model%climate%acab*thk0/tim0, & ! m/s + thck_unscaled, & ! m + model%isostasy%relx*thk0, & ! m + model%geometry%topg*thk0, & ! m + model%climate%eus*thk0) ! m ! Convert geometry%thck and calving%calving_thck to scaled model units model%geometry%thck(:,:) = thck_unscaled(:,:)/thk0 diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 892b6a18..9717f1d8 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -194,18 +194,22 @@ end subroutine glissade_calving_mask_init subroutine glissade_calve_ice(which_calving, & calving_domain, & + damage_src, & + damage_floor, & which_ho_calving_front, & remove_icebergs, & limit_marine_cliffs, & cull_calving_front, & + damage_manufactured, & + damage_advect, & calving, & ! calving derived type itest, jtest, rtest, & - dt, & ! s + time, dt, & ! s dx, dy, & ! m sigma, & thklim, & ! m efvs, & ! Pa s - acab, & ! m/s + uvel, acab, & ! m/s thck, relx, & ! m topg, eus) ! m @@ -226,11 +230,20 @@ subroutine glissade_calve_ice(which_calving, & !> = 1 if calving occurs everywhere the calving criterion is met !> = 2 if calving occurs where criterion is met and there is a connected path !> to the ocean through other cells where the criterion is met + integer, intent(in) :: damage_src !> option for type of damage src + !> = 0 if no damage source + !> = 1 for effective stress-based damage source + !> = 2 for Bassis and Ma 2015 damage source + integer, intent(in) :: damage_floor !> option for type of damage floor + !> = 0 if zero floor + !> = 1 if nye zero stress floor integer, intent(in) :: which_ho_calving_front !> = 1 for subgrid calving-front scheme, else = 0 logical, intent(in) :: remove_icebergs !> if true, then remove icebergs after calving logical, intent(in) :: limit_marine_cliffs !> if true, then limit the thickness of marine-based ice cliffs logical, intent(in) :: cull_calving_front !> if true, then cull calving_front cells to improve model stability; !> generally applied only at initialization + logical, intent(in) :: damage_manufactured !> if true, use a manufactured solution for damage + logical, intent(in) :: damage_advect !> if true, damage advects type(glide_calving), intent(inout) :: calving !> calving object @@ -263,11 +276,13 @@ subroutine glissade_calve_ice(which_calving, & ! real(dp), dimension(:,:), intent(out) :: calving_thck !> thickness lost due to calving in each grid cell (m) integer, intent(in) :: itest, jtest, rtest !> coordinates of diagnostic point + real(dp), intent(in) :: time !> current model time (s) real(dp), intent(in) :: dt !> model timestep (s) real(dp), intent(in) :: dx, dy !> grid cell size in x and y directions (m) real(dp), dimension(:), intent(in) :: sigma !> vertical sigma coordinate real(dp), intent(in) :: thklim !> minimum thickness for dynamically active grounded ice (m) real(dp), dimension(:,:,:), intent(in) :: efvs !> effective viscosity (Pa s) + real(dp), dimension(:,:,:), intent(in) :: uvel !> ice velocity along x (m/s) real(dp), dimension(:,:), intent(in) :: acab !> mass balance (m/s) real(dp), dimension(:,:), intent(inout) :: thck !> ice thickness (m) real(dp), dimension(:,:), intent(in) :: relx !> relaxed bedrock topography (m) @@ -284,23 +299,34 @@ subroutine glissade_calve_ice(which_calving, & real(dp), dimension(:,:,:), allocatable :: & eps_max, & ! maximum principal strain rate - d_damage_dt ! rate of change of damage scalar (1/s) + source, & ! damage evolution source term + vel_exp, & ! exponential used in computing manufactured damage + force, & ! forcing term used in computing manufactured damage + d_damage_dt, & ! rate of change of damage scalar (1/s) + adv_term, & ! advective term in manufactured damage solution + uvel_unstag, & ! unstaggered x velocity (m/s) + efvs_eff ! effective ice viscosity everywhere (subgrid damage only) real(dp), dimension(:,:), allocatable :: & thck_calving_front, & ! effective ice thickness at the calving front thck_init, & ! value of thck before calving tau1, tau2, & ! tau_eigen1 and tau_eigen2 (Pa), modified for calving - damage_column ! 2D vertically integrated scalar damage parameter + damage_column, & ! 2D vertically integrated scalar damage parameter + thck_eff ! effective ice thickness everywhere (subgrid damage only) real(dp), dimension(:,:), allocatable :: & calving_thck_init ! debug diagnostic + real(dp), dimension(:), allocatable :: & + hgl ! grounding line thickness (m), assuming primary ice flow along x + ! basic masks integer, dimension(:,:), allocatable :: & ice_mask, & ! = 1 where ice is present (thck > thklim), else = 0 floating_mask, & ! = 1 where ice is present (thck > thklim) and floating, else = 0 ocean_mask, & ! = 1 where topg is below sea level and ice is absent, else = 0 land_mask, & ! = 1 where topg is at or above sea level, else = 0 + grounding_line_mask, & ! = 1 where ice is adjacent to the grounding line, else = 0 active_ice_mask, & ! = 1 for cells that are dynamically active, else = 0 calving_front_mask, & ! = 1 where ice is floating and borders at least one ocean cell, else = 0 marine_cliff_mask ! = 1 where ice is grounded and marine-based and borders at least @@ -328,9 +354,10 @@ subroutine glissade_calve_ice(which_calving, & dthck, & ! thickness change (m) thckmax_cliff, & ! max stable ice thickness in marine_cliff cells factor, & ! factor in quadratic formula - beta, & ! flow law exponent anisotropy parameter + alpha, & ! flow law exponent anisotropy parameter nstar, & ! effective flow law exponent - szero ! ratio of hydrostatic to tensile stress + szero, & ! ratio of hydrostatic to tensile stress + damage_nye ! nye's zero stress damage real(dp), parameter :: & thinning_limit = 0.99d0 ! When ice not originally on the calving front is allowed to thin, @@ -385,6 +412,7 @@ subroutine glissade_calve_ice(which_calving, & allocate (floating_mask(nx,ny)) allocate (ocean_mask(nx,ny)) allocate (land_mask(nx,ny)) + allocate (grounding_line_mask(nx,ny)) allocate (active_ice_mask(nx,ny)) allocate (calving_front_mask(nx,ny)) allocate (thck_calving_front(nx,ny)) @@ -457,6 +485,7 @@ subroutine glissade_calve_ice(which_calving, & ice_mask, & floating_mask = floating_mask, & ocean_mask = ocean_mask, & + grounding_line_mask = grounding_line_mask, & which_ho_calving_front = which_ho_calving_front, & calving_front_mask = calving_front_mask, & thck_calving_front = thck_calving_front) @@ -519,24 +548,49 @@ subroutine glissade_calve_ice(which_calving, & call parallel_halo(tau1) call parallel_halo(tau2) + ! Create effective thickness and viscosity matrices + if (which_calving == CALVING_DAMAGE) then + allocate(thck_eff(nx,ny)) + allocate(efvs_eff(nz,nx,ny)) + do j = 2, ny-1 + do i = 2, nx-1 + if (calving_front_mask(i,j) == 1) then + thck_eff(i,j) = thck_calving_front(i,j) + else + thck_eff(i,j) = thck(i,j) + endif + efvs_eff(:,i,j) = efvs(:,i,j) + enddo + enddo + endif + ! In inactive calving-front cells where both eigenvalues are zero (because a cell is dynamically inactive), ! extrapolate nonzero values in upstream cells. + ! It is also necessary to extrapolate the viscosity in the subgrid damage case. ! Note: A similar extrapolation is done in glissade_diagnostic_variable_solve, but an extra one ! may be useful here for cells where ice was just advected from upstream. do j = 2, ny-1 do i = 2, nx-1 if (calving_front_mask(i,j) == 1) then - if (tau1(i,j) == 0.0d0 .and. tau2(i,j) == 0) then - do jj = j-1, j+1 - do ii = i-1, i+1 - if (thck_calving_front(i,j) > 0.0d0 .and. thck_init(ii,jj) == thck_calving_front(i,j)) then + ! Look for nonzero values in an upstream cell + do jj = j-1, j+1 + do ii = i-1, i+1 + if (thck_calving_front(i,j) > 0.0d0 & + .and. thck_init(ii,jj) == thck_calving_front(i,j)) then + if (tau1(i,j) == 0.0d0 .and. tau2(i,j) == 0) then tau1(i,j) = tau1(ii,jj) tau2(i,j) = tau2(ii,jj) - endif - enddo - enddo - endif ! tau1 = tau2 = 0 + endif ! tau1 = tau2 = 0 + + do k = 1, nz-1 + if (allocated(efvs_eff) .and. efvs(k,i,j) == 0.0d0) then + efvs_eff(k,i,j) = efvs(k,ii,jj) + endif ! efvs = 0 + enddo ! k + endif ! calving front thickness + enddo ! ii + enddo ! jj endif ! calving_front_mask enddo ! i enddo ! j @@ -544,41 +598,43 @@ subroutine glissade_calve_ice(which_calving, & ! Compute the effective stress. ! Note: By setting eigen2_weight > 1, we can give greater weight to the second principle stress. ! This may be useful in calving unbuttressed shelves that are spreading in both directions. + if (which_calving == EIGENCALVING .or. damage_src == EFF_STRESS_DAMAGE_SRC) then + calving%tau_eff(:,:) = sqrt(tau1(:,:)**2 + (calving%eigen2_weight * tau2(:,:))**2) - calving%tau_eff(:,:) = sqrt(tau1(:,:)**2 + (calving%eigen2_weight * tau2(:,:))**2) - - if (verbose_calving .and. this_rank == rtest) then - print*, ' ' - print*, 'tau1 (Pa), itest, jtest, rank =', itest, jtest, rtest - do j = jtest+3, jtest-3, -1 - write(6,'(i6)',advance='no') j - do i = itest-3, itest+3 - write(6,'(f10.2)',advance='no') tau1(i,j) + if (verbose_calving .and. this_rank == rtest) then + print*, ' ' + print*, 'tau1 (Pa), itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.2)',advance='no') tau1(i,j) + enddo + write(6,*) ' ' enddo - write(6,*) ' ' - enddo - print*, ' ' - print*, 'tau2 (Pa), itest, jtest, rank =', itest, jtest, rtest - do j = jtest+3, jtest-3, -1 - write(6,'(i6)',advance='no') j - do i = itest-3, itest+3 - write(6,'(f10.2)',advance='no') tau2(i,j) + print*, ' ' + print*, 'tau2 (Pa), itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.2)',advance='no') tau2(i,j) + enddo + write(6,*) ' ' enddo - write(6,*) ' ' - enddo - print*, ' ' - print*, 'tau_eff (Pa), itest, jtest, rank =', itest, jtest, rtest - do j = jtest+3, jtest-3, -1 - write(6,'(i6)',advance='no') j - do i = itest-3, itest+3 - write(6,'(f10.2)',advance='no') calving%tau_eff(i,j) + print*, ' ' + print*, 'tau_eff (Pa), itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.2)',advance='no') calving%tau_eff(i,j) + enddo + write(6,*) ' ' enddo - write(6,*) ' ' - enddo - endif + endif + endif ! which_calving / damage_src ! Use the effective stress either to directly compute a lateral calving rate (for eigencalving), - ! or to accumulate damage which is then used to derive a lateral calving rate (for damage-based calving). + ! or to accumulate damage which is then used to derive a lateral calving rate (for damage-based + ! calving with the effective stress damage source option only). calving%lateral_rate(:,:) = 0.0d0 @@ -596,43 +652,127 @@ subroutine glissade_calve_ice(which_calving, & elseif (which_calving == CALVING_DAMAGE) then - ! Prognose changes in damage. - ! For now, this is done using a simple scheme based on the effective tensile stress, calving%tau_eff + ! Prognose changes in damage based on the specified source term. ! The damage is subsequently advected downstream. ! Note: The damage is formally a 3D field, which makes it easier to advect, even though ! (in the current scheme) the damage source term is uniform in each column. - ! Allocate strain rate eigenvalue matrix - allocate(eps_max(nz,nx,ny)) + ! Allocate relevant data matrices + if (damage_src == BASSIS_MA_DAMAGE_SRC) then + allocate(eps_max(nz,nx,ny)) + allocate(source(nz,nx,ny)) + allocate(force(nz,nx,ny)) + endif + if (damage_manufactured) then + allocate(vel_exp(nz,nx,ny)) + allocate(adv_term(nz,nx,ny)) + allocate(uvel_unstag(nz,nx,ny)) + allocate(hgl(ny)) + + ! Identify the grounding line thickness for manufactured damage on + ! an ice tongue (assumes only one grounding line point per j) + do j = 2, ny-1 + do i = 2, nx-1 + if (floating_mask(i,j) == 0 .and. grounding_line_mask(i,j) == 1) then + hgl(j) = thck(i,j) + endif + enddo + enddo + endif ! damage_manufactured do j = 2, ny-1 do i = 2, nx-1 if (floating_mask(i,j) == 1) then - ! Compute the ratio of the principal stresses (equivalent to the ratio of the - ! principal strain rates) - beta = calving%tau_eigen2(i,j) / calving%tau_eigen1(i,j) + ! Prognose damage using a simple scheme based on the + ! effective tensile stress, calving%tau_eff + if (damage_src == EFF_STRESS_DAMAGE_SRC) then + d_damage_dt = calving%damage_constant * calving%tau_eff(i,j) ! damage_constant has units of s^{-1}/(Pa) + + ! Prognose damage based on eq 27 in Bassis and Ma 2015 + else if (damage_src == BASSIS_MA_DAMAGE_SRC) then + ! Compute the ratio of the principal stresses (equivalent to the ratio of the + ! principal strain rates) + alpha = calving%tau_eigen2(i,j) / calving%tau_eigen1(i,j) + + ! Find the max principal strain rate using Glen's Flow Law + eps_max(:,i,j) = calving%tau_eigen1(i,j) / (2.0d0 * efvs_eff(:,i,j)) + + ! Compute the effective flow law exponent + nstar = 12.0d0*(1.0d0+alpha+alpha**2)/(4.0d0*(1.0d0+alpha+alpha**2)+6.0d0*alpha**2) + + ! Compute the hydrostatic to tensile stress ratio + szero = rhoi*(rhoo-rhoi)*grav*thck_eff(i,j)/(2.0d0*calving%tau_eigen1(i,j)*rhoo) + + ! Prognose damage following Bassis & Ma 2015 + source(:,i,j) = nstar*(1.0d0-szero)*eps_max(:,i,j)-acab(i,j)/thck_eff(i,j) + d_damage_dt(:,i,j) = source(:,i,j)*calving%damage(:,i,j) + + ! Compute the manufactured forcing term if requested; otherwise, set it to zero + if (damage_manufactured) then + ! Compute the advection term + adv_term(:,i,j) = 1.0d0 + if (damage_advect) then + adv_term(:,i,j) = adv_term(:,i,j)+eps_max(:,i,j)*time + endif - ! Find the max principal strain rate using Glen's Flow Law - eps_max(:,i,j) = calving%tau_eigen1(i,j) / (2.0d0 * efvs(:,i,j)) + ! Linearly interpolate uvel to the unstaggered grid + uvel_unstag(:,i,j) = 0.25d0*abs(uvel(:,i,j)+uvel(:,i-1,j)+uvel(:,i,j-1) & + +uvel(:,i-1,j-1)) + + ! Compute the manufactured forcing term + vel_exp(:,i,j) = exp(-uvel_unstag(:,i,j)*time/hgl(j)) + force(:,i,j) = -0.5d0*calving%damage_init(:,i,j)*(uvel_unstag(:,i,j)/hgl(j) & + *adv_term(:,i,j)*vel_exp(:,i,j) & + +source(:,i,j)*(1.0d0+vel_exp(:,i,j))) + else + force(:,i,j) = 0.0d0 + endif - ! Compute the effective flow law exponent - nstar = 12.0d0*(1.0d0+beta+beta**2) / (4.0d0*(1.0d0+beta+beta**2)+6.0d0*beta**2) + d_damage_dt(:,i,j) = d_damage_dt(:,i,j) + force(:,i,j) + calving%damage(:,i,j) = calving%damage(:,i,j) + d_damage_dt(:,i,j) * dt - ! Compute the hydrostatic to tensile stress ratio - szero = rhoi*(rhoo-rhoi)*grav*thck(i,j)/(2.0d0*calving%tau_eigen1(i,j)*rhoo) + endif ! damage_src - ! Prognose damage following Bassis & Ma 2015 - d_damage_dt(:,i,j) = (nstar*(1.0d0-szero)*eps_max(:,i,j) & - - acab(i,j)/thck(i,j)) * calving%damage(:,i,j) -! d_damage_dt = calving%damage_constant * calving%tau_eff(i,j) ! damage_constant has units of s^{-1}/(Pa) - calving%damage(:,i,j) = calving%damage(:,i,j) + d_damage_dt(:,i,j) * dt + ! Constrain damage to the specified range of values calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) - calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) - else ! set damage to zero for grounded ice + if (damage_floor == ZERO_DAMAGE_FLOOR) then + calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) + elseif (damage_floor == NYE_DAMAGE_FLOOR) then + ! Compute the Nye zero stress value + damage_nye = (2.0d0+alpha)*calving%tau_eigen1(i,j)/(grav*(rhoo-rhoi)*thck_eff(i,j)) + calving%damage(:,i,j) = max(calving%damage(:,i,j), damage_nye) + endif + endif ! floating_mask + enddo ! i + + do i = 2, nx-1 + ! Make sure that damage is zero where there isn't ice and where + ! the ice is grounded + if (thck(i,j) <= thklim) then + calving%damage(:,i,j) = 0.0d0 + elseif (floating_mask(i,j) == 0) then calving%damage(:,i,j) = 0.0d0 endif - enddo - enddo + + ! When the Nye damage floor is not enforced, we need to sweep + ! across i a second time to get the damage along grounding lines + ! right. This is important because damage advected into the + ! domain is set to match the damage of what's already there, + ! which for ice tongues happens to be the zero damage along the + ! grounding line. We want new damage to be as damaged as + ! existing ice, so we need the grounding line to also be as + ! damaged as the ice further downstream. + if (damage_floor == ZERO_DAMAGE_FLOOR) then + ! If i,j is the grounding line + if (floating_mask(i,j) == 0 .and. grounding_line_mask(i,j) == 1) then + ! If i+1,j is floating and ice exists there + if (floating_mask(i+1,j) == 1 .and. thck(i+1,j) > thklim) then + calving%damage(:,i,j) = calving%damage(:,i+1,j) + endif + endif + endif ! damage_floor + enddo ! i + enddo ! j ! Compute the vertically integrated damage in each column. allocate(damage_column(nx,ny)) @@ -646,41 +786,44 @@ subroutine glissade_calve_ice(which_calving, & enddo enddo - ! Convert damage in CF cells to a lateral calving rate (m/s). - ! Note: Although eigenprod = 0 in inactive calving-front cells, these cells can have significant damage - ! advected from upstream, so in general we should not have to interpolate damage from upstream. - !TODO - Verify this. - do j = 2, ny-1 - do i = 2, nx-1 - if (calving_front_mask(i,j) == 1) then - frac_lateral = (damage_column(i,j) - calving%damage_threshold) / (1.0d0 - calving%damage_threshold) - frac_lateral = max(0.0d0, min(1.0d0, frac_lateral)) - calving%lateral_rate(i,j) = calving%lateral_rate_max * frac_lateral ! m/s - endif + if (which_ho_calving_front == HO_CALVING_FRONT_SUBGRID) then + ! Convert damage in CF cells to a lateral calving rate (m/s). + ! Note: Although eigenprod = 0 in inactive calving-front cells, these cells can have significant damage + ! advected from upstream, so in general we should not have to interpolate damage from upstream. + !TODO - Verify this. + do j = 2, ny-1 + do i = 2, nx-1 + if (calving_front_mask(i,j) == 1) then + frac_lateral = (damage_column(i,j) - calving%damage_threshold) / (1.0d0 - calving%damage_threshold) + frac_lateral = max(0.0d0, min(1.0d0, frac_lateral)) + calving%lateral_rate(i,j) = calving%lateral_rate_max * frac_lateral ! m/s + endif + enddo enddo - enddo - if (verbose_calving .and. this_rank==rtest) then - print*, ' ' - print*, 'damage increment, itest, jtest, rank =', itest, jtest, rtest - do j = jtest+3, jtest-3, -1 - write(6,'(i6)',advance='no') j - do i = itest-3, itest+3 - write(6,'(f10.6)',advance='no') calving%damage_constant * calving%tau_eff(i,j) * dt + if (verbose_calving .and. this_rank==rtest) then + print*, ' ' + print*, 'damage increment, itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.6)',advance='no') calving%damage_constant * calving%tau_eff(i,j) * dt + enddo + write(6,*) ' ' enddo - write(6,*) ' ' - enddo - print*, ' ' - print*, 'new damage, itest, jtest, rank =', itest, jtest, rtest - do j = jtest+3, jtest-3, -1 - write(6,'(i6)',advance='no') j - do i = itest-3, itest+3 - write(6,'(f10.6)',advance='no') damage_column(i,j) + print*, ' ' + print*, 'new damage, itest, jtest, rank =', itest, jtest, rtest + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.6)',advance='no') damage_column(i,j) + enddo + write(6,*) ' ' enddo - write(6,*) ' ' - enddo - print*, ' ' - endif + print*, ' ' + endif + + endif ! which_ho_calving_front endif ! EIGENCALVING or CALVING_DAMAGE @@ -1400,6 +1543,7 @@ subroutine glissade_calve_ice(which_calving, & deallocate (floating_mask) deallocate (ocean_mask) deallocate (land_mask) + deallocate (grounding_line_mask) deallocate (active_ice_mask) deallocate (calving_front_mask) deallocate (thck_calving_front) @@ -1412,6 +1556,14 @@ subroutine glissade_calve_ice(which_calving, & if (allocated(tau1)) deallocate(tau1) if (allocated(tau2)) deallocate(tau2) if (allocated(eps_max)) deallocate(eps_max) + if (allocated(source)) deallocate(source) + if (allocated(vel_exp)) deallocate(vel_exp) + if (allocated(force)) deallocate(force) + if (allocated(hgl)) deallocate(hgl) + if (allocated(adv_term)) deallocate(adv_term) + if (allocated(uvel_unstag)) deallocate(uvel_unstag) + if (allocated(thck_eff)) deallocate(thck_eff) + if (allocated(efvs_eff)) deallocate(efvs_eff) end subroutine glissade_calve_ice diff --git a/libglissade/glissade_transport.F90 b/libglissade/glissade_transport.F90 index 394e7f65..a38e448b 100644 --- a/libglissade/glissade_transport.F90 +++ b/libglissade/glissade_transport.F90 @@ -107,7 +107,7 @@ subroutine glissade_transport_setup_tracers(model, & model%geometry%ntracers = model%geometry%ntracers + 1 endif - if (model%options%whichcalving == CALVING_DAMAGE) then + if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_advect) then model%geometry%ntracers = model%geometry%ntracers + 1 endif @@ -164,7 +164,7 @@ subroutine glissade_transport_setup_tracers(model, & endif ! damage parameter for prognostic calving scheme - if (model%options%whichcalving == CALVING_DAMAGE) then + if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_advect) then nt = nt + 1 do k = 1, nlyr @@ -246,7 +246,7 @@ subroutine glissade_transport_finish_tracers(model) endif ! damage parameter for prognostic calving scheme - if (model%options%whichcalving == CALVING_DAMAGE) then + if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_advect) then nt = nt + 1 do k = 1, nlyr model%calving%damage(k,:,:) = model%geometry%tracers(:,:,nt,k) From 7dbba496f98a333ef0d8428faf0acd63b5237c49 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Mon, 8 Apr 2019 15:42:35 -0400 Subject: [PATCH 3/9] Ice tongue test case This commit adds a test case designed for damage testing that constructs a one-dimensional, freely-floating ice tongue following the analytic equations developed in Chapter 5.5 of van der Veen [2013]. The test case itself is housed in tests/tongue and contains an executable python script (ice_tongue.py), a CISM configuration file (ice_tongue.config), and documentation in the tongue/doc folder (which also describes the details of the new damage_manufactured configuration option). In addition, the python script can be toggled to instead construct a two-dimensional ice tongue confined to a lateral embayment. The ice tongue geometries pre-loaded within this test case are: - Erebus - The unconfined portion of Drygalski - Idealized Amery - Idealized Ross --- tests/tongue/doc/README.tex | 130 ++++++++++ tests/tongue/doc/mybib.bib | 43 ++++ tests/tongue/ice_tongue.config | 57 ++++ tests/tongue/ice_tongue.py | 457 +++++++++++++++++++++++++++++++++ tests/tongue/netCDF.py | 124 +++++++++ 5 files changed, 811 insertions(+) create mode 100644 tests/tongue/doc/README.tex create mode 100644 tests/tongue/doc/mybib.bib create mode 100644 tests/tongue/ice_tongue.config create mode 100755 tests/tongue/ice_tongue.py create mode 100644 tests/tongue/netCDF.py diff --git a/tests/tongue/doc/README.tex b/tests/tongue/doc/README.tex new file mode 100644 index 00000000..5aba36e2 --- /dev/null +++ b/tests/tongue/doc/README.tex @@ -0,0 +1,130 @@ +\documentclass{article} +\title{README - Ice Tongue Test Case} +\author{Morgan Whitcomb, University of Michigan} +\date{\today} + +\usepackage{natbib} +\usepackage{booktabs} + +\begin{document} + +\maketitle + +\section{Introduction} + +This test case is designed to construct a one-dimensional, freely-floating ice tongue following Chapter 5.5 in \citet{van-der-Veen} for damage testing, while also providing an option to add embayment walls (making the ice tongue two-dimensional and no longer freely-floating). The script \texttt{ice\_tongue.py} feeds an input netCDF file to CISM, in which the ice tongue is constructed via a set of specified variables: +\begin{itemize} + \item \texttt{thk} - Analytic ice thickness (m) with a geometry-specific grounding line value + \item \texttt{uvel\_extend} - Extended analytic ice velocity (m/a) with a geometry-specific grounding line value; the extended version of the variable is necessary in order to avoid boundary condition errors where the ice tongue contacts the edge of the domain + \item \texttt{damage} - Scalar damage field, seeded as zero or a constant value (specified by the \texttt{const\_damage} variable in the python script) depending on the value of the \texttt{damage\_floor} configuration option + \item \texttt{acab} - Mass balance (m/a), set as a uniform geometry-specific basal melt rate + \item \texttt{topg} - Elevation of the basal topography (m), set to a uniform value beneath the ice tongue that grounds the row of grid cells furthest upstream as the grounding line but causes every other grid cell to float; when the ice tongue length is specified as variable, the elevation of the topography is increased to match the grounding line thickness along the lateral margins to simulate embayment walls + \item \texttt{beta} - Higher-order bed stress, set to zero + \item \texttt{kinbcmask} - Kinetic boundary condition mask, used to enforce a constant grounding line velocity + \item \texttt{no\_advance\_mask} - No advance mask, used to prevent the ice tongue from flowing upstream, and from advancing if specified by the input options + \item \texttt{const\_thk\_mask} - Constant thickness mask, used to enforce a constant grounding line thickness +\end{itemize} + +The user is also provided with a set of command-line input options that allow them to tweak some parameters of the simulation: +Input options: +\begin{itemize} + \item -c, --config : Name of the \texttt{.config} file to use for the simulation (default = \texttt{ice\_tongue.config}). + \item -m, --parallel : Execute the run in parallel with the specified number of processors (default = serial run). + \item -e, --exec : Path to the CISM executable file (default = \texttt{.\textbackslash cism\_driver}) + \item -s, --shelf : Ice shelf geometry to use for the simulation, which determines the grounding line flux, uniform basal melt rate, and length of embayment walls (confined ice tongues only). The current supported options are \texttt{erebus}, \texttt{drygalski}, \texttt{ross} (confined), and \texttt{amery} (confined); for the specific parameter values used to define each geometry, see Table~\ref{table:geometries} (default = \texttt{erebus}). + \item -l, --length : Toggle for variability in the length of the ice tongue; the \texttt{const} option contrains the length, while the \texttt{var} option allows the ice tongue to advance freely (default = \texttt{var}) + \item -r, --restart : Path to a netCDF restart file that will be used as the input file (default = none) +\end{itemize} + +\begin{table} + \centering + \begin{tabular}{l c c c c} + \toprule + & \multicolumn{2}{c}{Unconfined} & \multicolumn{2}{c}{Confined} \\ + \cmidrule(lr){2-3} + \cmidrule(lr){4-5} + Geometry-Specific Parameters & Erebus & Drygalski & Amery & Ross \\ + \cmidrule(lr){1-1} + \cmidrule(lr){2-2} + \cmidrule(lr){3-3} + \cmidrule(lr){4-4} + \cmidrule(lr){5-5} + \texttt{nx} & $73$ & $391$ & $301$ & $388$ \\ + \texttt{ny} & $50$ & $160$ & $60$ & $574$ \\ + \texttt{dx} and \texttt{dy} (m) & $500$ & $500$ & $5000$ & $5000$ \\ + \texttt{default\_flwa} ($10^{-17}$ Pa$^{-n}$ a$^{-1}$) & $1.294$ & $1.294$ & $1.105$ & $1.105$ \\ + $H_0$ (m) & $340$ & $575$ & $1090$ & $660$ \\ + $U_0$ (m/a) & $100$ & $940$ & $390$ & $340$ \\ + \texttt{acab} (m/a) & $2.2$ & $6.8$ & Variable & Variable \\ + Confinement length (km) & N/A & N/A & $505$ & $650$ \\ + \bottomrule + \end{tabular} + \caption{Characteristic parameters used to define the geometries specified by the -s input option. All of the entries except for \texttt{nx}, \texttt{ny}, \texttt{dx}, \texttt{dy}, and \texttt{default\_flwa} are specified in the \texttt{ice\_tongue.py} script; the exceptions must be set manually in the configuration file. For the confined ice tongues (Amery and Ross), the melt rates can be set by editing the \texttt{amery\_melt} and \texttt{ross\_melt} variables in the python script. See Whitcomb et al. [in press] for information about where these values come from. } + \label{table:geometries} +\end{table} + +\section{Ice Tongues} + +For a positive melt rate $\dot{m}$, the one-dimensional thickness $H\left(x\right)$ of the ice tongue is computed via Eq.~5.80 in \citet{van-der-Veen}: +\begin{equation} + H\left(x\right) = \left\{\frac{U_0^{n+1}\left[1+\left(C/\dot{m}\right)H_0^{n+1}\right]}{\left(H_0U_0 - \dot{m}x\right)^{n+1}} - \frac{C}{\dot{m}}\right\}^{-1/\left(n+1\right)} + \label{eq:thickness} +\end{equation} +where $H_0$ and $U_0$ are the grounding line thickness and velocity, respectively, and $n=3$ is the Glen's flow law exponent. The constant $C$ is defined by Eq.~5.65 in \citet{van-der-Veen} as +\begin{equation} + C = \left[\frac{\rho_ig\left(\rho_w-\rho_i\right)}{4B\rho_w}\right]^n, + \label{eq:C-constant} +\end{equation} +where $\rho_i$ is the ice density, $\rho_w$ is the density of seawater, $g$ is the acceleration due to gravity, and $B$ is the viscosity parameter in Glen's flow law. + +Given Eq.~\ref{eq:thickness}, the velocity $U\left(x\right)$ is computed directly from conservation of ice flux: +\begin{equation} + HU = H_0U_0 - \dot{m}x + \label{eq:conservation-of-flux} +\end{equation} + +\section{Damage} + +As stated above, this test case is specifically designed to verify damage $r$ as computed by Eq.~27 in \citet{Bassis-Ma}: +\begin{equation} + \frac{\mathrm{d}r}{\mathrm{d}t} = \left[n^* \left(1-S_0\right) \dot{\epsilon}_{xx} + \frac{\dot{m}}{H}\right] r + \label{eq:bassis-ma-damage} +\end{equation} +where $\dot{\epsilon}_{xx}$ is the largest principal strain rate and $n^*$ is an effective flow law exponent: +\begin{equation} + n^* = \frac{4n\left(1+\alpha+\alpha^2\right)}{4\left(1+\alpha+\alpha^2\right) + 3\left(n-1\right)\alpha^2} + \label{eq:nstar} +\end{equation} +The parameter $\alpha = \dot{\epsilon}_{yy}/\dot{\epsilon}_{xx}$ is the ratio of the principal strain rates, and $S_0$ is a dimensionless ratio between hydrostatic pressure and tensile stress: +\begin{equation} + S_0 = \frac{\rho_i \left(\rho_w-\rho_i\right)gH}{2\tau_{xx}\rho_w} + \label{eq:szero} +\end{equation} +where $\tau_{xx}$ is the largest principal deviatoric stress. + +Eq.~\ref{eq:bassis-ma-damage} takes in an initial damage field, evolving it in time and advecting it with the flow field. Depending on the value of the \texttt{damage\_floor} configuration option, the initial value $r_0$ is specified as either a uniform value (set as the \texttt{const\_damage} variable in the python script) or as the Nye's zero stress value, computed via \citep{Jezek,Nick-et-al,Nye} +\begin{equation} + r_0 = \frac{\rho_i}{\rho_w-\rho_i} \left[\frac{\left(2+\alpha\right) \tau_{xx}}{\rho_i gH}\right]. + \label{eq:nye-damage} +\end{equation} + +Verification of the damage solution can be performed by toggling the \texttt{damage\_} \texttt{manufactured} configuration option. When activated, this option adjusts Eq.~\ref{eq:bassis-ma-damage} by adding a forcing term $f$: +\begin{equation} + \frac{\mathrm{d}r}{\mathrm{d}t} = \left[n^* \left(1-S_0\right) \dot{\epsilon}_{xx} + \frac{\dot{m}}{H}\right] r + f + \label{eq:bassis-ma-damage-forced} +\end{equation} +This forcing term is defined as +\begin{equation} + f = -\frac{r_0}{2} \left\{\frac{U}{H_0} \left(1+\dot{\epsilon}_{xx}t\right) e^{-Ut/H_0} + \left[n^*\left(1-S_0\right)\dot{\epsilon}_{xx}+\frac{\dot{m}}{H}\right] \left(1+e^{-Ut/H_0}\right)\right\}, + \label{eq:damage-forcing} +\end{equation} +such that Eq.~\ref{eq:bassis-ma-damage-forced} produces a damage field given by +\begin{equation} + r = \frac{r_0}{2} \left(1+e^{-Ut/H_0}\right). + \label{eq:manufactured-damage} +\end{equation} + +\bibliographystyle{plainnat} +\bibliography{mybib} + +\end{document} diff --git a/tests/tongue/doc/mybib.bib b/tests/tongue/doc/mybib.bib new file mode 100644 index 00000000..dcf38814 --- /dev/null +++ b/tests/tongue/doc/mybib.bib @@ -0,0 +1,43 @@ +@Book{van-der-Veen, + Author = "van der Veen, C. J.", + Title = "Fundamentals of Glacier Dynamics", + Publisher = "CRC Press", + Year = 2013, + Edition = 2, +} + +@Article{Bassis-Ma, + Author = "Bassis, J. N. and Ma, Y.", + Title = "Evolution of basal crevasses links ice shelf stability to ocean forcing", + Journal = "Earth and Planetary Science Letters", + Volume = 409, + Pages = "203-211", + Year = 2015, +} + +@Article{Jezek, + Author = "Jezek, K. C.", + Title = "A modified theory of bottom crevasses used as a means for measuring the buttressing effect of ice shelves on inland ice sheets", + Journal = "Journal of Geophysical Research", + Volume = 89, + Pages = "1925-1931", + Year = 1984, +} + +@Article{Nick-et-al, + Author = "Nick, F. M. and van der Veen, C. J. and Vieli, A. and Benn, D. I.", + Title = "A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics", + Journal = "Journal of Glaciology", + Volume = 56, + Pages = "781-794", + Year = 2010, +} + +@Article{Nye, + Author = "Nye, J. F.", + Title = "The distribution of stress and velocity in glaciers and ice-sheets", + Journal = "Proceedings of the Royal Society A", + Volume = 239, + Pages="113-133", + Year = 1957, +} diff --git a/tests/tongue/ice_tongue.config b/tests/tongue/ice_tongue.config new file mode 100644 index 00000000..61d5928a --- /dev/null +++ b/tests/tongue/ice_tongue.config @@ -0,0 +1,57 @@ +[grid] +ewn = 73 +nsn = 50 +upn = 5 +dew = 500 +dns = 500 +global_bc = 2 # 0 = GLOBAL_BC_PERIODIC, 1 = outflow, 2 = GLOBAL_BC_NO_PENETRATION + +[time] +tstart = 0. +tend = 300. +dt = 0.1 +dt_diag = 10. +idiag = 5 +jdiag = 5 + +[options] +dycore = 2 # 1 = glam, 2 = glissade +flow_law = 0 # 0 = isothermal, 2 = temperature dependent +evolution = 3 # 3 = inc. remapping, 4 = FO upwind, 5 = none +temperature = 2 # 0 = set column to surf. air temp, 1 = prognostic, 2 = hold at init. values +marine_margin = 8 # 0 = do nothing, 8 = damage +calving_init = 0 # 0 = calve after first time step, 1 = calve at initialization +calving_domain = 0 # 0 = ocean edge, 1 = everywhere, 2 = ocean connect +damage_src = 2 # 0 = none, 1 = effective stress, 2 = bassis and ma +damage_floor = 1 # 0 = none, 1 = nye +damage_advect = true +damage_manufactured = false + +[ho_options] +which_ho_babc = 5 # 5 = take basal traction param value from .nc input +which_ho_efvs = 2 # 2 = nonlinear eff. visc. w/ n=3 +which_ho_sparse = 3 # 1 = SLAP GMRES, 3 = glissade parallel PCG, 4 = Trilinos for linear solver +which_ho_nonlinear = 0 # 0 = Picard, 1=JFNK +which_ho_approx = 1 # 1 = SSA, 2 = Blatter-Pattyn, 3 = L1L2 +which_ho_precond = 1 # 0 = none, 1 = diagonal, 2 = SIA based (ONLY use for which_ho_approx = 2) +which_ho_vertical_remap = 1 # 0 = first order, 1 = second order +which_ho_calving_front = 1 # 0 = no subgrid, 1 = subgrid +glissade_maxiter = 500 + +[parameters] +default_flwa = 1.294e-17 +damage_threshold = 0.99 # Calving criterion value +calving_timescale = 0.1 # Calving timescale + +[CF default] +title = Ice Tongue Test Case + +[CF input] +name = tongue_input.nc +time = 1 + +[CF output] +variables = thk uvel vvel damage kinbcmask const_thk_mask no_advance_mask efvs tau_xx tau_yy eps_xx eps_yy topg acab +frequency = 1.0 +name = results/erebus.out.nc + diff --git a/tests/tongue/ice_tongue.py b/tests/tongue/ice_tongue.py new file mode 100755 index 00000000..98d1a2f1 --- /dev/null +++ b/tests/tongue/ice_tongue.py @@ -0,0 +1,457 @@ +#!/usr/bin/env python +# This script runs a one-dimensional ice tongue test case used for verifying damage. +# Pre-loaded ice shelves: Erebus & Drygalski (unconfined); Amery & Ross (confined) +# Last edited by Morgan Whitcomb at the University of Michigan on April 8, 2019. + +# Imports +from optparse import OptionParser +from ConfigParser import ConfigParser +from netCDF import * +import numpy as np +import sys,os,datetime + +# User-defined inputs ## +const_damage = 0.1 # Constant initial damage value +#amery_melt = 0.6 # Melt rate across the Amery ice shelf (m/a) +amery_melt = 0.8 # +#amery_melt = 1. # +#ross_melt = 0.2 # Melt rate across the Ross ice shelf (m/a) +ross_melt = 0.25 # +#ross_melt = 0.5 # +## + +# Classes ##### +class constants: + # Structure for holding general constants + def __init__(s): + s.n = 3. # Glen's Flow Law exponent + s.rhoi = 910. # Ice density (kg/m^3) + s.rhoo = 1028. # Seawater density (kg/m^3) + s.grav = 9.81 # Gravitational acceleration (m/s^2) + +class geometry_vals: + # Structure for holding parameters specific to the geometry of interest + def __init__(s,options,c,nx,dx,ny,dy,nz,flwa,dinit,amery_melt,ross_melt): + # Initialize resolution values + s.nx = nx + s.nxst = nx-1 + s.dx = dx + s.ny = ny + s.nyst = ny-1 + s.dy = dy + s.nz = nz + s.nzst = nz-1 + + # Compute vectors and matrices of spatial grid points + s.x = s.dx*np.arange(s.nx,dtype='float32') + s.xst = s.x[:-1] + s.dx/2. + s.y = s.dy*np.arange(s.ny,dtype='float32') + s.yst = s.y[:-1] + s.dy/2. + s.xgrid,s.ygrid = np.meshgrid(s.x,s.y,indexing='xy') + s.xstgrid,s.ystgrid = np.meshgrid(s.xst,s.yst,indexing='xy') + + # Initialize geometry-specific parameters + if options.shelf == 'erebus': + s.set_shelf_params(340.,100.,2.2) + elif options.shelf == 'drygalski': + s.set_shelf_params(575.,940.,6.8) + elif options.shelf == 'amery': + s.set_shelf_params(1090.,390.,amery_melt,conflen=505.) + elif options.shelf == 'ross': + s.set_shelf_params(660.,340.,ross_melt,conflen=650.) + s.flwa = flwa # Ice viscosity (Pa^[-n] yr^[-1]) + s.C = s.flwa*(c.rhoi*c.grav*(c.rhoo-c.rhoi)/(4.*c.rhoo))**c.n + + # Set the initial damage field + s.dinit = dinit + + # Specify matrix indices + s.ptnum = int((s.nx+2)/3) # Number of points along the initial ice tongue + s.gl = 2 # Grounding line location + s.cf = s.ptnum+2 # Calving front location + + def set_shelf_params(s,hgl,ugl,mdot,conflen=0.): + # Load the ice shelf grounding line flux and melt rate into the struct + s.hgl = hgl # Grounding line thickness (m) + s.ugl = ugl # Grounding line velocity (m/a) + s.mdot = mdot # Melt rate (m/a) + s.conflen = conflen # Downstream extent of lateral embayment walls (km) + s.confend = int(s.conflen*1000./s.dx) # Matrix index corresponding to conflen + + def construct_ice_tongue(s,options,c): + # Compute the necessary fluid fields and establish boundary conditions + # Initialize matrices + s.thk = np.zeros([1,s.ny,s.nx],dtype='float32') # Thickness (m) + s.acab = np.zeros([1,s.ny,s.nx],dtype='float32') # Mass balance (m/a) + s.topg = np.zeros([1,s.ny,s.nx],dtype='float32') # Bed elevation (m) + s.beta = np.zeros([1,s.nyst,s.nxst],dtype='float32') # Higher-order bed stress + s.kbc = np.zeros([1,s.nyst,s.nxst],dtype='int') # Kinetic boundary condition mask + s.no_advance_mask = np.zeros([1,s.ny,s.nx],dtype='int') # No advance mask + s.const_thk_mask = np.zeros([1,s.ny,s.nx],dtype='int') # Constant thickness mask + # We need to use the extended staggered mesh for our velocities to avoid boundary + # condition problems - the fact that the standard staggered mesh is smaller than the + # unstaggered mesh by 1 point per dimension causes numerical issues when using + # periodic global boundaries with nonzero velocities at the walls. + s.uvel_extend = np.zeros([1,s.nz,s.ny,s.nx],dtype='float32') # X-velocity (m/a) + + # Set the climate forcing + s.acab[0,:,s.gl+1:] = -s.mdot # Mass balance is applied everywhere but the gl + + # Prevent the ice from flowing backwards from the gl and from advancing beyond the cf + s.no_advance_mask[0,:,:s.gl] = 1 # Backwards from gl + if options.length == 'const': + s.no_advance_mask[0,:,s.cf+1:] = 1 # Advancing from cf + + # Set the kinetic boundary conditions to hold the gl velocity constant in time + s.kbc[0,:,:s.gl+1] = 1 + + # Hold the gl thickness constant in time + s.const_thk_mask[0,:,s.gl] = 1 + + # Compute the ice thickness + s.thk[0,:,s.gl] = s.hgl + for j in range(0,s.ny): + s.thk[0,j,s.gl+1:s.cf+1] = s.compute_thk(c,s.x[s.gl+1:s.cf+1],s.x[s.gl]) + + # Compute the ice velocity + s.compute_uvel(c) + + # Lower the bedrock elevation everywhere to make the ice shelf float everywhere but the gl + # We have to raise topg by a small amount in order to account for floating point inaccuracies, + # which we've chosen here to use half the difference between the grounding line thickness and the + # thickness in the grid cell immediately downstream from it + for i in range(0,s.nx): + s.topg[0,:,i] = -(c.rhoi/c.rhoo)*s.hgl + 0.5*(s.hgl-s.thk[0,:,s.gl+1]) + + # Enforce lateral confinement + if s.conflen > 0.: # If laterally confined + s.enforce_confinement(s.acab[0],0.,'acab',s.confend) + s.enforce_confinement(s.thk[0],0.,'thk',min(s.cf,s.confend)) + s.enforce_confinement(s.uvel_extend[0],0.,'uvel',min(s.cf,s.confend)) + s.enforce_confinement(s.topg[0],s.hgl,'topg',s.confend) + s.enforce_confinement(s.kbc[0],1,'kbc',s.confend) + + def compute_thk(s,c,x,xgl): + # Compute the analytic ice thickness, via van der Veen Ch. 5 + if s.mdot == 0.: # Accumulation = ablation + thk = ((c.n+1.)*s.C*(x-xgl)/(s.hgl*s.ugl)+s.hgl**(-(c.n+1.)))**(-1./(c.n+1.)) + + else: # Nonzero accumulation + thk = (s.ugl**(c.n+1.)*(s.C*s.hgl**(c.n+1.)/s.mdot+1.) \ + /(s.hgl*s.ugl-s.mdot*(x-xgl))**(c.n+1.)-s.C/s.mdot)**(-1./(c.n+1.)) + + return thk + + def compute_uvel(s,c): + # Compute the analytic ice velocity from conservation of flux + # Calculate the thickness beyond the calving front for smooth interpolation + thktemp = s.thk[0,0,s.gl:s.cf+2].copy() + thktemp[-1] = s.compute_thk(c,s.x[s.cf+1],s.x[s.gl]) + + # Linearly interpolate the thickness onto the staggered grid + thk_stag = np.zeros(len(thktemp[:-1])) + thk_stag = 0.5*(thktemp[:-1]+thktemp[1:]) + + # Compute the velocity + for i in range(s.gl,s.cf+1): + s.uvel_extend[0,:,:,i] = (s.hgl*s.ugl-s.mdot*(s.xst[i]-s.x[s.gl]))/thk_stag[i-s.gl] + + def enforce_confinement(s,field,val,fieldstr,xend): + # Enforce lateral confinement in the specified fluid field + twod_fields = ['acab','thk','topg','kbc'] + threed_fields = ['uvel','damage'] + if fieldstr in twod_fields: + field[0,:xend+1] = val + field[-1,:xend+1] = val + elif fieldstr in threed_fields: + field[:,0,:xend+1] = val + field[:,-1,:xend+1] = val + + if fieldstr == 'uvel': # Need to include an extra row for uvel_extend + field[:,-2,:xend+1] = val + + # Save the altered field to the struct + if fieldstr == 'acab': + s.acab[0] = field + elif fieldstr == 'thk': + s.thk[0] = field + elif fieldstr == 'uvel': + s.uvel_extend[0] = field + elif fieldstr == 'topg': + s.topg[0] = field + elif fieldstr == 'damage': + s.damage[0] = field + +class restart_file: + # Structure for holding data from the restart file, when specified + def __init__(s,g,options): + message('Copying restart data from "'+options.restfile+'"') + + # Open the restart file + infile = open_netcdf_file(options.restfile,'r') + + # Extract the grid and time information from the restart file + s.x = np.array(infile.variables['x1'][:]) + s.xst = np.array(infile.variables['x0'][:]) # Staggered grid + s.y = np.array(infile.variables['y1'][:]) + s.yst = np.array(infile.variables['y0'][:]) # Staggered grid + s.z = np.array(infile.variables['level'][:]) + s.zst = np.array(infile.variables['staglevel'][:]) # Staggered grid + s.t = np.array(infile.variables['time'][:]) + + # Parse the spatial and temporal vectors to find the resolution + s.nx,s.dx = s.find_res(s.x) + s.nxst = len(s.xst) + s.ny,s.dy = s.find_res(s.y) + s.nyst = len(s.yst) + s.nz = len(s.z) + s.nzst = len(s.zst) + s.nt,s.dt = s.find_res(s.t) + + # Error checking + if s.nx != g.nx or s.dx != g.dx or s.ny != g.ny or s.dy != g.dy or s.nz != g.nz: + # If the resolutions don't match, quit + error('the resolution specified in the configuration file must match that of the restart file') + + # Extract the necessary variables from the restart file + s.thk = np.array(infile.variables['thk'][-1,:,:]) # Ice thickness (m) + s.uvel = np.array(infile.variables['uvel'][-1,:,:,:]) # X-velocity (m/a) + s.vvel = np.array(infile.variables['vvel'][-1,:,:,:]) # Y-velocity (m/a) + s.kbc = np.array(infile.variables['kinbcmask'][-1,:,:]) # Kinetic boundary condition mask + s.const_thk_mask = np.array(infile.variables['const_thk_mask'][-1,:,:]) # Constant thickness mask + s.no_advance_mask = np.array(infile.variables['no_advance_mask'][-1,:,:]) + s.efvs = np.array(infile.variables['efvs'][-1,:,:,:]) # Ice viscosity + s.tau_xx = np.array(infile.variables['tau_xx'][-1,:,:,:]) # Deviatoric X-stress (Pa) + s.tau_yy = np.array(infile.variables['tau_yy'][-1,:,:,:]) # Deviatoric Y-stress (Pa) + s.eps_xx = np.array(infile.variables['eps_xx'][-1,:,:,:]) # X-strain rate (1/a) + s.eps_yy = np.array(infile.variables['eps_yy'][-1,:,:,:]) # Y-strain rate (1/a) + s.topg = np.array(infile.variables['topg'][-1,:,:]) # Bed topography (m) + s.acab = np.array(infile.variables['acab'][-1,:,:]) # Mass balance (m/a) + + # Create the higher-order bed stress - zeros everywhere for floating ice + s.beta = np.zeros([1,s.nyst,s.nxst],dtype='float32') # Higher-order bed stress + + # Create extended versions of the velocity matrices by directly copying from nearest-neighbor + # interpolation + s.uvel_extend = s.create_extended_vel(s.uvel) + s.vvel_extend = s.create_extended_vel(s.vvel) + + def find_res(s,i): + # Pull the number of points and resolution from a vector + ni = len(i) + di = i[1]-i[0] + + return ni,di + + def create_extended_vel(s,vel): + # Create an extended version of a velocity matrix + vel_extend = np.zeros([1,s.nz,s.ny,s.nx],dtype='float32') + + # Copy the non-extended velocity into the extended matrix + vel_extend[0,:,:-1,:-1] = vel + + # Fill the remaining values by nearest-neighbor interpolation + vel_extend[0,:,-1,:-1] = vel[:,-1,:] + vel_extend[0,:,:-1,-1] = vel[:,:,-1] + vel_extend[0,:,-1,-1] = vel[:,-1,-1] + + return vel_extend +##### + +# Functions ##### +def cmdline_opts(cfgfile='ice_tongue.config'): + # Parse command-line options + optparser = OptionParser() + optparser.add_option('-c','--config',dest='configfile',type='string',default=cfgfile, \ + help='Name of .config file to use for the run', metavar='FILE') + optparser.add_option('-m','--parallel',dest='parallel',type='int', \ + help='Number of processors to run the model with; if specified, then execute'+\ + ' the run in parallel [default: perform a serial run]', metavar="NUMPROCS") + optparser.add_option('-e','--exec',dest='executable',type='string',default='./cism_driver', \ + help='Set path to the CISM executable',metavar='FILE') + optparser.add_option('-s','--shelf',dest='shelf',type='string',default='erebus',help='Ice shelf to'+\ + ' use representative parameters for; the supported options are: erebus,'+\ + ' drygalski, ross, amery') + optparser.add_option('-l','--length',dest='length',type='string',default='var',help="Mode for"+\ + " determining the final length of the ice tongue; const constrains the ice"+\ + " tongue's maximum length, and var allows the ice tongue to advance freely") + optparser.add_option('-r','--restart',dest='restfile',type='string', \ + help='Specify a netCDF to restart from with a fresh damage field') + for option in optparser.option_list: + if option.default != ("NO", "DEFAULT"): + option.help += (" " if option.help else "") + "[default: %default]" + options,_ = optparser.parse_args() + + # Error checking + if os.path.isfile(options.configfile) is False: # Config file doesn't exist + error('file "'+options.configfile+'" does not exist') + if os.path.isfile(options.executable) is False: # CISM executable doesn't exist + error('file "'+options.executable+'" does not exist') + # Restart file was specified but doesn't exist + if options.restfile is not None and os.path.isfile(options.restfile) is False: + error('file "'+options.restfile+'" does not exist') + shelf_list = ['erebus','drygalski','ross','amery'] + if options.shelf not in shelf_list: + error('Invalid input for shelf; accepted entries are: erebus, drygalski, ross, amery') + length_list = ['const','var'] + if options.length not in length_list: + error('Invalid input for length; accepted entries are: const, var') + + return options + +def error(string): + # Send an error to the user and exit the program + print 'ERROR: '+string + sys.exit() + +def message(string): + print string+' ... ', datetime.datetime.now().time() + +def parse_config_file(options,amery_melt,ross_melt,no_dam_floor=0,nye_dam_floor=1): + # Read resolution data from the config file + parser = ConfigParser() + parser.read(options.configfile) + nx = int(parser.get('grid','ewn')) + ny = int(parser.get('grid','nsn')) + nz = int(parser.get('grid','upn')) + dx = float(parser.get('grid','dew')) + dy = float(parser.get('grid','dns')) + filename = str(parser.get('CF input','name')) + flwa = float(parser.get('parameters','default_flwa')) + dam_floor = parser.get('options','damage_floor') + dam_floor = int(dam_floor.split(' ')[0]) # Remove post-option comments + if dam_floor == no_dam_floor: # No damage floor (constant input value) + dinit = const_damage + elif dam_floor == nye_dam_floor: # Nye damage floor + dinit = 0. + + # Create a geometry class instance and save the resolution data to it + c = constants() + geom = geometry_vals(options,c,nx,dx,ny,dy,nz,flwa,dinit,amery_melt,ross_melt) + + return geom,c,filename + +def open_netcdf_file(filename,action): + # Open a netCDF file for the specified action + try: + openfile = NetCDFFile(filename,action,format='NETCDF3_CLASSIC') + except TypeError: + openfile = NetCDFFile(filename,action) + + return openfile + +def init_netcdf(s,filename): + # Initialize input netCDF file for CISM, using information from the configuration file + # Create the netCDF file + message('Writing '+filename) + infile = open_netcdf_file(filename,'w') + + # Create netCDF dimensions for the resolution values + infile.createDimension('time',1) + infile.createDimension('x1',s.nx) + infile.createDimension('x0',s.nxst) # Staggered grid + infile.createDimension('y1',s.ny) + infile.createDimension('y0',s.nyst) # Staggered grid + infile.createDimension('level',s.nz) + infile.createDimension('staglevel',s.nzst) + + # Create netCDF vectors for the corresponding dimensions + infile.createVariable('time','f',('time',))[:] = [0] + infile.createVariable('x1','f',('x1',))[:] = s.x + infile.createVariable('x0','f',('x0',))[:] = s.xst # Staggered grid + infile.createVariable('y1','f',('y1',))[:] = s.y + infile.createVariable('y0','f',('y0',))[:] = s.yst # Staggered grid + + return infile + +def initialize_damage(s,g): + # Initialize the damage field across the ice tongue + # Create the damage matrix + s.damage = np.zeros([1,s.nzst,s.ny,s.nx],dtype='float32') # Damage + + # Set the initial damage + s.damage[0,:,:,g.gl:g.cf+1] = g.dinit + + # Enforce lateral confinement + if g.conflen > 0.: # If laterally confined + g.enforce_confinement(s.damage[0],0.,'damage',min(g.cf,g.confend)) + +def finalize_netcdf(s,infile,restart=False): + # Finalize the input netCDF file for outputting to CISM + message('Saving parameter matrices to input netCDF file') + + # Create the required variables in the netCDF file + infile.createVariable('thk','f',('time','y1','x1'))[:] = s.thk + infile.createVariable('uvel_extend','f',('time','level','y1','x1'))[:] = s.uvel_extend + infile.createVariable('damage','f',('time','staglevel','y1','x1'))[:] = s.damage + infile.createVariable('acab','f',('time','y1','x1'))[:] = s.acab + infile.createVariable('topg','f',('time','y1','x1'))[:] = s.topg + infile.createVariable('beta','f',('time','y0','x0'))[:] = s.beta + infile.createVariable('kinbcmask','i',('time','y0','x0'))[:] = s.kbc + infile.createVariable('no_advance_mask','i',('time','y1','x1'))[:] = s.no_advance_mask + infile.createVariable('const_thk_mask','i',('time','y1','x1'))[:] = s.const_thk_mask + if restart == True: + infile.createVariable('vvel_extend','f',('time','level','y1','x1'))[:] = s.vvel_extend + infile.createVariable('efvs','f',('time','staglevel','y1','x1'))[:] = s.efvs + infile.createVariable('tau_xx','f',('time','staglevel','y1','x1'))[:] = s.tau_xx + infile.createVariable('tau_yy','f',('time','staglevel','y1','x1'))[:] = s.tau_yy + infile.createVariable('eps_xx','f',('time','staglevel','y1','x1'))[:] = s.eps_xx + infile.createVariable('eps_yy','f',('time','staglevel','y1','x1'))[:] = s.eps_yy + + infile.close() + +def run_cism(options,test_type='unconfined'): + # Submit the input netCDF file for CISM to run + print 'Running CISM for the '+test_type+' ice-tongue experiment' + print '==============================================\n' + if options.parallel == None: + # Perform a serial run + runstring = options.executable + ' ' + options.configfile + print 'Executing serial run with: ' + runstring + '\n\n' + os.system(runstring) + else: + # Perform a parallel run + if options.parallel <= 0: + error('Number of processors specified for parallel run is <=0.') + else: + # These calls to os.system will return the exit status: + # 0 for success (the command exists), some other integer for failure + if os.system('which openmpirun > /dev/null') == 0: + mpiexec = 'openmpirun -np ' + str(options.parallel) + elif os.system('which mpirun > /dev/null') == 0: + mpiexec = 'mpirun -np ' + str(options.parallel) + elif os.system('which aprun > /dev/null') == 0: + mpiexec = 'aprun -n ' + str(options.parallel) + elif os.system('which mpirun.lsf > /dev/null') == 0: + # mpirun.lsf does NOT need the number of processors (options.parallel) + mpiexec = 'mpirun.lsf' + else: + error('Unable to execute parallel run. Please edit the script to use your MPI run'+\ + ' command, or run the model manually with something like: mpirun -np 4 ./cism_driver'+\ + ' ice-tongue.config') + runstring = mpiexec + ' ' + options.executable + ' ' + options.configfile + print 'Executing parallel run with: ' + runstring + '\n\n' + os.system(runstring) # Here is where the parallel run is actually executed! +##### + +# MAIN CODE ##### +options = cmdline_opts() +g,c,filename = parse_config_file(options,amery_melt,ross_melt) + +# When we're not restarting from a previous file, create an input netCDF file and construct an ice +# tongue; when we are restarting, copy the input data to a new input netCDF file and create a new +# damage field. +if options.restfile is None: # Not a restart + infile = init_netcdf(g,filename) + g.construct_ice_tongue(options,c) + initialize_damage(g,g) + finalize_netcdf(g,infile) +else: # Restarting + r = restart_file(g,options) + infile = init_netcdf(r,filename) + initialize_damage(r,g) + finalize_netcdf(r,infile,restart=True) + +run_cism(options) +##### + diff --git a/tests/tongue/netCDF.py b/tests/tongue/netCDF.py new file mode 100644 index 00000000..c9da3a7f --- /dev/null +++ b/tests/tongue/netCDF.py @@ -0,0 +1,124 @@ +# This script allows use of any of three python netCDF modules: +# Scientific.IO.NetCDF, netCDF4, or pycdf +# To use whichever netCDF module you might have installed put +# from netCDF import * +# in your script. +# Programs should use the Scientific.IO.NetCDF syntax; +# Generally, netCDF4 matches the Scientific.IO.NetCDF syntax and functionality. However there are some differences which require knowing which module is in use to properly call methods. The variable netCDF_module is provided to accomplish this. +# If the pycdf module is to be used, an appropriate "translation" of the method calls is provided. +# Written March 16, 2010 by Glen Granzow + +try: + from Scientific.IO.NetCDF import NetCDFFile + netCDF_module = 'Scientific.IO.NetCDF' +except ImportError: + try: + from netCDF4 import Dataset as NetCDFFile + netCDF_module = 'netCDF4' + except ImportError: + try: + import pycdf + netCDF_module = 'pycdf' + except ImportError: + print 'Unable to import any of the following python modules:' + print ' Scientific.IO.NetCDF \n netcdf4 \n pycdf' + print 'One of them must be installed.' + raise ImportError('No netCDF module found') + + def NCtype(value): + if isinstance(value,int): return pycdf.NC.INT + if isinstance(value,float): return pycdf.NC.FLOAT + if isinstance(value,str): return pycdf.NC.CHAR + + class NetCDFFile(object): + def __init__(self,filename,mode): + if mode == 'w': + self.FILE = pycdf.CDF(filename, pycdf.NC.WRITE | pycdf.NC.CREATE | pycdf.NC.TRUNC) + if mode == 'r': + self.FILE = pycdf.CDF(filename, pycdf.NC.NOWRITE) + if mode == 'a': + self.FILE = pycdf.CDF(filename, pycdf.NC.WRITE) + self.FILE.automode() + + def __setattr__(self,name,value): + if name == 'FILE': + object.__setattr__(self,name,value) + else: # Used to assign global attributes + self.FILE.attr(name).put(NCtype(value),value) + + def __getattr__(self,name): + if name == 'dimensions': + return self.FILE.dimensions() + if name == 'variables': + dictionary = dict() + for variable in self.FILE.variables().keys(): + dictionary[variable] = NetCDFvariable(self.FILE.var(variable),None,None,None) + return dictionary + global_attributes = self.FILE.attributes() + if name in global_attributes: + return global_attributes[name] + return object.__getattribute__(self,name) + + def __dir__(self): + return self.FILE.attributes().keys() + + def hasattr(self,name): + return name in dir() + + def createDimension(self,name,size): + self.FILE.def_dim(name,size) + + def createVariable(self,name,datatype,dimensions): + dictNC = {'f':pycdf.NC.FLOAT, 'd':pycdf.NC.DOUBLE, 'i':pycdf.NC.INT, 'c':pycdf.NC.CHAR, 'b':pycdf.NC.BYTE} + return NetCDFvariable(self.FILE,name,dictNC[datatype],dimensions) + + def sync(self): + self.FILE.sync() + + def close(self): + self.FILE.close() + + class NetCDFvariable(object): + def __init__(self,FILE,name,datatype,dimensions): + if isinstance(FILE,pycdf.pycdf.CDFVar): + # FILE is an already defined netCDF variable (not a file) + self.VARIABLE = FILE + else: + # Create a new variable in the netCDF file + self.VARIABLE = FILE.def_var(name,datatype,dimensions) + self.shape = self.VARIABLE.shape() + + def __setitem__(self,elem,data): + self.VARIABLE[elem] = data + + def __getitem__(self,elem): + return self.VARIABLE[elem] + + def __setattr__(self,name,value): + if name in ('VARIABLE','shape'): + object.__setattr__(self,name,value) + else: # Used to assign variable attributes + self.VARIABLE.attr(name).put(NCtype(value),value) + + def __getattr__(self,name): + if name == 'dimensions': + return self.VARIABLE.dimensions() + variable_attributes = self.VARIABLE.attributes() + if name in variable_attributes: + return variable_attributes[name] + return object.__getattribute__(self,name) + + def assignValue(self,value): + self.VARIABLE.put(value) + + def getValue(self): + return self.VARIABLE.get() + + def __dir__(self): + return self.VARIABLE.attributes().keys() + + def typecode(self): + NCdict = {pycdf.NC.FLOAT:'f', pycdf.NC.DOUBLE:'d', pycdf.NC.INT:'i', pycdf.NC.CHAR:'c', pycdf.NC.BYTE:'b'} + return NCdict[self.VARIABLE.inq_type()] + +print 'Using',netCDF_module,'for netCDF file I/O' From 98d5e9383bcb19cb15f8a55491fa96ab39a167b9 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Wed, 24 Apr 2019 12:54:11 -0400 Subject: [PATCH 4/9] Minor bugfixes General: - Damage should now never be larger than one In the ice tongue test case: - Damage confinement is now enforced for Amery and Ross when using a restart file - The bedrock elevation has been lowered slightly to fix pockets of zero damage downstream from the grounding line for Amery and Ross --- libglissade/glissade_calving.F90 | 3 +- tests/tongue/ice_tongue.py | 90 +++++++++++++++++++------------- 2 files changed, 55 insertions(+), 38 deletions(-) diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 9717f1d8..2dced985 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -734,7 +734,6 @@ subroutine glissade_calve_ice(which_calving, & endif ! damage_src ! Constrain damage to the specified range of values - calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) if (damage_floor == ZERO_DAMAGE_FLOOR) then calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) elseif (damage_floor == NYE_DAMAGE_FLOOR) then @@ -742,6 +741,8 @@ subroutine glissade_calve_ice(which_calving, & damage_nye = (2.0d0+alpha)*calving%tau_eigen1(i,j)/(grav*(rhoo-rhoi)*thck_eff(i,j)) calving%damage(:,i,j) = max(calving%damage(:,i,j), damage_nye) endif + calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) + endif ! floating_mask enddo ! i diff --git a/tests/tongue/ice_tongue.py b/tests/tongue/ice_tongue.py index 98d1a2f1..55baa7ac 100755 --- a/tests/tongue/ice_tongue.py +++ b/tests/tongue/ice_tongue.py @@ -117,19 +117,14 @@ def construct_ice_tongue(s,options,c): s.compute_uvel(c) # Lower the bedrock elevation everywhere to make the ice shelf float everywhere but the gl - # We have to raise topg by a small amount in order to account for floating point inaccuracies, - # which we've chosen here to use half the difference between the grounding line thickness and the - # thickness in the grid cell immediately downstream from it - for i in range(0,s.nx): - s.topg[0,:,i] = -(c.rhoi/c.rhoo)*s.hgl + 0.5*(s.hgl-s.thk[0,:,s.gl+1]) + s.topg[0] = compute_topg(s,c,s,s.topg[0],s.thk[0]) # Enforce lateral confinement if s.conflen > 0.: # If laterally confined - s.enforce_confinement(s.acab[0],0.,'acab',s.confend) - s.enforce_confinement(s.thk[0],0.,'thk',min(s.cf,s.confend)) - s.enforce_confinement(s.uvel_extend[0],0.,'uvel',min(s.cf,s.confend)) - s.enforce_confinement(s.topg[0],s.hgl,'topg',s.confend) - s.enforce_confinement(s.kbc[0],1,'kbc',s.confend) + enforce_confinement(s,s.acab[0],0.,'acab',s.confend) + enforce_confinement(s,s.thk[0],0.,'thk',min(s.cf,s.confend)) + enforce_confinement(s,s.uvel_extend[0],0.,'uvel',min(s.cf,s.confend)) + enforce_confinement(s,s.kbc[0],1,'kbc',s.confend) def compute_thk(s,c,x,xgl): # Compute the analytic ice thickness, via van der Veen Ch. 5 @@ -156,32 +151,6 @@ def compute_uvel(s,c): for i in range(s.gl,s.cf+1): s.uvel_extend[0,:,:,i] = (s.hgl*s.ugl-s.mdot*(s.xst[i]-s.x[s.gl]))/thk_stag[i-s.gl] - def enforce_confinement(s,field,val,fieldstr,xend): - # Enforce lateral confinement in the specified fluid field - twod_fields = ['acab','thk','topg','kbc'] - threed_fields = ['uvel','damage'] - if fieldstr in twod_fields: - field[0,:xend+1] = val - field[-1,:xend+1] = val - elif fieldstr in threed_fields: - field[:,0,:xend+1] = val - field[:,-1,:xend+1] = val - - if fieldstr == 'uvel': # Need to include an extra row for uvel_extend - field[:,-2,:xend+1] = val - - # Save the altered field to the struct - if fieldstr == 'acab': - s.acab[0] = field - elif fieldstr == 'thk': - s.thk[0] = field - elif fieldstr == 'uvel': - s.uvel_extend[0] = field - elif fieldstr == 'topg': - s.topg[0] = field - elif fieldstr == 'damage': - s.damage[0] = field - class restart_file: # Structure for holding data from the restart file, when specified def __init__(s,g,options): @@ -364,6 +333,26 @@ def init_netcdf(s,filename): return infile +def compute_topg(s,c,g,topg,thk): + # Compute the bedrock elevation + # + # We have to raise topg by a small amount in order to account for floating point inaccuracies, + # which we've chosen here to use quarter the difference between the grounding line thickness and the + # thickness in the grid cell immediately downstream from it + # + # This needs to be done again in the case of a restart file because the thickness is no longer + # laterally uniform, and that is how topg is initially defined with the guess thickness field. + # If it's not re-calculated, what this means is that parts of the ice shelf that do not lie along + # the grounding line will be treated as grounded rather than floating. + for i in range(0,s.nx): + topg[:,i] = -(c.rhoi/c.rhoo)*g.hgl + 0.25*(g.hgl-thk[:,g.gl+1]) + + # Enforce lateral confinement + if g.conflen > 0.: # If laterally confined + enforce_confinement(s,topg,g.hgl,'topg',g.confend) + + return topg + def initialize_damage(s,g): # Initialize the damage field across the ice tongue # Create the damage matrix @@ -374,7 +363,33 @@ def initialize_damage(s,g): # Enforce lateral confinement if g.conflen > 0.: # If laterally confined - g.enforce_confinement(s.damage[0],0.,'damage',min(g.cf,g.confend)) + enforce_confinement(s,s.damage[0],0.,'damage',min(g.cf,g.confend)) + +def enforce_confinement(s,field,val,fieldstr,xend): + # Enforce lateral confinement in the specified fluid field + twod_fields = ['acab','thk','topg','kbc'] + threed_fields = ['uvel','damage'] + if fieldstr in twod_fields: + field[0,:xend+1] = val + field[-1,:xend+1] = val + elif fieldstr in threed_fields: + field[:,0,:xend+1] = val + field[:,-1,:xend+1] = val + + if fieldstr == 'uvel': # Need to include an extra row for uvel_extend + field[:,-2,:xend+1] = val + + # Save the altered field to the struct + if fieldstr == 'acab': + s.acab[0] = field + elif fieldstr == 'thk': + s.thk[0] = field + elif fieldstr == 'uvel': + s.uvel_extend[0] = field + elif fieldstr == 'topg': + s.topg = field + elif fieldstr == 'damage': + s.damage[0] = field def finalize_netcdf(s,infile,restart=False): # Finalize the input netCDF file for outputting to CISM @@ -449,6 +464,7 @@ def run_cism(options,test_type='unconfined'): else: # Restarting r = restart_file(g,options) infile = init_netcdf(r,filename) + r.topg = compute_topg(r,c,g,r.topg,r.thk) initialize_damage(r,g) finalize_netcdf(r,infile,restart=True) From 0e666fa0c98df6912b6e3f59bb8974749a09a36f Mon Sep 17 00:00:00 2001 From: morgwhit Date: Wed, 24 Apr 2019 13:12:36 -0400 Subject: [PATCH 5/9] Slight doc fix Fixed the lateral sizes of Ross and Amery in the table of the ice tongue test case documentation --- tests/tongue/doc/README.tex | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tongue/doc/README.tex b/tests/tongue/doc/README.tex index 5aba36e2..55c67889 100644 --- a/tests/tongue/doc/README.tex +++ b/tests/tongue/doc/README.tex @@ -50,7 +50,7 @@ \section{Introduction} \cmidrule(lr){4-4} \cmidrule(lr){5-5} \texttt{nx} & $73$ & $391$ & $301$ & $388$ \\ - \texttt{ny} & $50$ & $160$ & $60$ & $574$ \\ + \texttt{ny} & $50$ & $160$ & $22$ & $193$ \\ \texttt{dx} and \texttt{dy} (m) & $500$ & $500$ & $5000$ & $5000$ \\ \texttt{default\_flwa} ($10^{-17}$ Pa$^{-n}$ a$^{-1}$) & $1.294$ & $1.294$ & $1.105$ & $1.105$ \\ $H_0$ (m) & $340$ & $575$ & $1090$ & $660$ \\ From 0b3360610f7a6e7ec87056e9dc26f59292fbf0e2 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Tue, 13 Aug 2019 16:33:13 -0400 Subject: [PATCH 6/9] Significant bug fixes in the physics/discretization + support for subgrid calving front scheme Fixed a variety of bugs in the Bassis and Ma damage evolution scheme, both in the physics and in the numerical discretization: - Damage is now converted to crevasse depth before being advected and evolved, to match the form of the advective equation solved by the incremental remapping scheme; it is then converted back to damage afterwards - Adjusted the damage source term, manufactured forcing, and Nye's zero stress criterion to match the aforementioned crevasse depth parameterization - Updated the ice tongue test case documentation to match the reformulated damage evolution parameterization - Replaced thck with thck_old in the damage evolution scheme to be consistent with the forward Euler discretization - Removed thck_eff, which only existed because of a misinterpretation of how the subgrid calving front scheme works - Fixed the linear extrapolation routine at the calving front for the velocities and viscosity, and added one for damage at the grounding line (which the model currently does not have the capacity to handle) - Identified and corrected loop indexing bugs in the damage evolution scheme and added the appropriate calls to parallel_halo In addition, support was added for the subgrid calving front scheme to the manufactured damage solver. To-do: Fix discrepancy in the damage field near the grounding line. --- libglide/glide_setup.F90 | 8 +- libglissade/glissade.F90 | 12 +- libglissade/glissade_calving.F90 | 295 ++++++++++++++++++++----------- tests/tongue/doc/README.tex | 22 ++- 4 files changed, 223 insertions(+), 114 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 06464e6d..913b07e6 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -1867,13 +1867,7 @@ subroutine print_parameters(model) call write_log(message) endif - if (model%options%whichcalving == CALVING_DAMAGE .and. model%options%damage_manufactured) then - model%options%which_ho_calving_front = HO_CALVING_FRONT_NO_SUBGRID - write(message,*) 'Setting which_ho_calving_front =', HO_CALVING_FRONT_NO_SUBGRID - call write_log(message) - write(message,*) 'Manufactured damage calving option does not support subgrid calving front scheme' - call write_log(message) - elseif (model%options%whichcalving == CALVING_DAMAGE) then + if (model%options%whichcalving == CALVING_DAMAGE) then write(message,*) 'Subgrid calving front scheme chosen: damage will be converted into a lateral calving rate' call write_log(message) endif diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 499d050a..a2605633 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -2039,7 +2039,7 @@ subroutine glissade_calving_solve(model, init_calving) logical :: cull_calving_front ! true iff init_calving = T and options%cull_calving_front = T - integer :: i, j + integer :: i, j, k !TODO - Make sure no additional halo updates are needed before glissade_calve_ice @@ -2055,7 +2055,14 @@ subroutine glissade_calving_solve(model, init_calving) ! Pass in thck, topg, etc. with units of meters. ! ------------------------------------------------------------------------ - thck_unscaled(:,:) = model%geometry%thck(:,:)*thk0 + ! Convert damage to crevasse depth + if (model%options%whichcalving == CALVING_DAMAGE) then + do k = 1, size(model%numerics%stagsigma) + model%calving%damage(k,:,:) = model%calving%damage(k,:,:)*model%geometry%thck_old(:,:)*thk0 + enddo + endif + + thck_unscaled(:,:) = model%geometry%thck_old(:,:)*thk0 call glissade_calve_ice(model%options%whichcalving, & model%options%calving_domain, & @@ -2079,6 +2086,7 @@ subroutine glissade_calving_solve(model, init_calving) model%numerics%thklim*thk0, & ! m model%stress%efvs*evs0, & ! Pa s model%velocity%uvel*vel0, & ! m/s + model%velocity%vvel*vel0, & ! m/s model%climate%acab*thk0/tim0, & ! m/s thck_unscaled, & ! m model%isostasy%relx*thk0, & ! m diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 2dced985..f74dcb3d 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -209,7 +209,8 @@ subroutine glissade_calve_ice(which_calving, & sigma, & thklim, & ! m efvs, & ! Pa s - uvel, acab, & ! m/s + uvel, vvel, & ! m/s + acab, & ! m/s thck, relx, & ! m topg, eus) ! m @@ -283,6 +284,7 @@ subroutine glissade_calve_ice(which_calving, & real(dp), intent(in) :: thklim !> minimum thickness for dynamically active grounded ice (m) real(dp), dimension(:,:,:), intent(in) :: efvs !> effective viscosity (Pa s) real(dp), dimension(:,:,:), intent(in) :: uvel !> ice velocity along x (m/s) + real(dp), dimension(:,:,:), intent(in) :: vvel !> ice velocity along y (m/s) real(dp), dimension(:,:), intent(in) :: acab !> mass balance (m/s) real(dp), dimension(:,:), intent(inout) :: thck !> ice thickness (m) real(dp), dimension(:,:), intent(in) :: relx !> relaxed bedrock topography (m) @@ -295,6 +297,7 @@ subroutine glissade_calve_ice(which_calving, & integer :: nz ! number of vertical levels ! Note: number of ice layers = nz-1 integer :: i, j, k, n + integer :: ud, vd ! Velocity sign index modulations integer :: ii, jj real(dp), dimension(:,:,:), allocatable :: & @@ -305,14 +308,15 @@ subroutine glissade_calve_ice(which_calving, & d_damage_dt, & ! rate of change of damage scalar (1/s) adv_term, & ! advective term in manufactured damage solution uvel_unstag, & ! unstaggered x velocity (m/s) - efvs_eff ! effective ice viscosity everywhere (subgrid damage only) + vvel_unstag, & ! unstaggered y velocity (m/s) + speed_unstag, & ! unstaggered horizontal speed (m/s) + efvs_eff ! effective viscosity including the calving front (Pa s) real(dp), dimension(:,:), allocatable :: & thck_calving_front, & ! effective ice thickness at the calving front thck_init, & ! value of thck before calving tau1, tau2, & ! tau_eigen1 and tau_eigen2 (Pa), modified for calving - damage_column, & ! 2D vertically integrated scalar damage parameter - thck_eff ! effective ice thickness everywhere (subgrid damage only) + damage_column ! 2D vertically integrated scalar damage parameter real(dp), dimension(:,:), allocatable :: & calving_thck_init ! debug diagnostic @@ -354,6 +358,8 @@ subroutine glissade_calve_ice(which_calving, & dthck, & ! thickness change (m) thckmax_cliff, & ! max stable ice thickness in marine_cliff cells factor, & ! factor in quadratic formula + xext, & ! extrapolated x value + yext, & ! extrapolated x value alpha, & ! flow law exponent anisotropy parameter nstar, & ! effective flow law exponent szero, & ! ratio of hydrostatic to tensile stress @@ -535,9 +541,14 @@ subroutine glissade_calve_ice(which_calving, & allocate(tau1(nx,ny)) allocate(tau2(nx,ny)) - ! Ignore negative eigenvalues corresponding to compressive stresses - tau1 = max(calving%tau_eigen1, 0.0d0) - tau2 = max(calving%tau_eigen2, 0.0d0) + if (which_calving == EIGENCALVING) then + ! Ignore negative eigenvalues corresponding to compressive stresses + tau1 = max(calving%tau_eigen1, 0.0d0) + tau2 = max(calving%tau_eigen2, 0.0d0) + else + tau1 = calving%tau_eigen1 + tau2 = calving%tau_eigen2 + endif ! Ignore values on grounded ice where (floating_mask == 0) @@ -548,53 +559,6 @@ subroutine glissade_calve_ice(which_calving, & call parallel_halo(tau1) call parallel_halo(tau2) - ! Create effective thickness and viscosity matrices - if (which_calving == CALVING_DAMAGE) then - allocate(thck_eff(nx,ny)) - allocate(efvs_eff(nz,nx,ny)) - do j = 2, ny-1 - do i = 2, nx-1 - if (calving_front_mask(i,j) == 1) then - thck_eff(i,j) = thck_calving_front(i,j) - else - thck_eff(i,j) = thck(i,j) - endif - efvs_eff(:,i,j) = efvs(:,i,j) - enddo - enddo - endif - - ! In inactive calving-front cells where both eigenvalues are zero (because a cell is dynamically inactive), - ! extrapolate nonzero values in upstream cells. - ! It is also necessary to extrapolate the viscosity in the subgrid damage case. - ! Note: A similar extrapolation is done in glissade_diagnostic_variable_solve, but an extra one - ! may be useful here for cells where ice was just advected from upstream. - - do j = 2, ny-1 - do i = 2, nx-1 - if (calving_front_mask(i,j) == 1) then - ! Look for nonzero values in an upstream cell - do jj = j-1, j+1 - do ii = i-1, i+1 - if (thck_calving_front(i,j) > 0.0d0 & - .and. thck_init(ii,jj) == thck_calving_front(i,j)) then - if (tau1(i,j) == 0.0d0 .and. tau2(i,j) == 0) then - tau1(i,j) = tau1(ii,jj) - tau2(i,j) = tau2(ii,jj) - endif ! tau1 = tau2 = 0 - - do k = 1, nz-1 - if (allocated(efvs_eff) .and. efvs(k,i,j) == 0.0d0) then - efvs_eff(k,i,j) = efvs(k,ii,jj) - endif ! efvs = 0 - enddo ! k - endif ! calving front thickness - enddo ! ii - enddo ! jj - endif ! calving_front_mask - enddo ! i - enddo ! j - ! Compute the effective stress. ! Note: By setting eigen2_weight > 1, we can give greater weight to the second principle stress. ! This may be useful in calving unbuttressed shelves that are spreading in both directions. @@ -657,22 +621,99 @@ subroutine glissade_calve_ice(which_calving, & ! Note: The damage is formally a 3D field, which makes it easier to advect, even though ! (in the current scheme) the damage source term is uniform in each column. + ! Set initial dummy values + d_damage_dt(:,:,:) = 0.0d0 + ! Allocate relevant data matrices if (damage_src == BASSIS_MA_DAMAGE_SRC) then + allocate(uvel_unstag(nz,nx,ny)) + allocate(vvel_unstag(nz,nx,ny)) + allocate(speed_unstag(nz,nx,ny)) + allocate(efvs_eff(nz-1,nx,ny)) allocate(eps_max(nz,nx,ny)) allocate(source(nz,nx,ny)) allocate(force(nz,nx,ny)) - endif + + ! Set initial dummy values + uvel_unstag(:,:,:) = 0.0d0 + vvel_unstag(:,:,:) = 0.0d0 + speed_unstag(:,:,:) = 0.0d0 + efvs_eff(:,:,:) = efvs(:,:,:) + eps_max(:,:,:) = 0.0d0 + source(:,:,:) = 0.0d0 + force(:,:,:) = 0.0d0 + + ! Construct unstaggered velocity matrices + do j = 3, ny-2 + do i = 3, nx-2 + ! Linearly interpolate velocities to the unstaggered grid + uvel_unstag(:,i,j) = 0.25d0*abs(uvel(:,i,j)+uvel(:,i-1,j)+uvel(:,i,j-1) & + +uvel(:,i-1,j-1)) + vvel_unstag(:,i,j) = 0.25d0*abs(vvel(:,i,j)+vvel(:,i-1,j)+vvel(:,i,j-1) & + +vvel(:,i-1,j-1)) + + ! Compute the magnitude of the horizontal velocity vectors + speed_unstag(:,i,j) = SQRT(uvel_unstag(:,i,j)**2+vvel_unstag(:,i,j)**2) + enddo ! i + enddo ! j + + call parallel_halo(uvel_unstag) + call parallel_halo(vvel_unstag) + + ! We need to restart the loops here to have the matrices fully computed for extrapolation + ! We want to only sweep over points that exist (so we need a buffer of at + ! least two in the halo). + do j = 3, ny-2 + do i = 3, nx-2 + ! To do this right we also need to extrapolate the velocities to + ! the calving front, since there are too few points on the + ! staggered grid to fully interpolate + if (calving_front_mask(i,j) == 1) then + do k = 1, nz + ! Store the velocity directions as matrix index modulations + ! e.g. if uvel is positive, set ud = -1 because ice will flow + ! from i-1 to i + ud = find_velocity_sign(uvel_unstag(k,i,j)) + vd = find_velocity_sign(vvel_unstag(k,i,j)) + + uvel_unstag(k,i,j) = linear_extrapolation(uvel_unstag(k,i+ud,j), & + uvel_unstag(k,i+2*ud,j)) + vvel_unstag(k,i,j) = linear_extrapolation(vvel_unstag(k,i,j+vd), & + vvel_unstag(k,i,j+2*vd)) + + ! Compute the magnitude of the horizontal velocity vectors + speed_unstag(k,i,j) = SQRT(uvel_unstag(k,i,j)**2+vvel_unstag(k,i,j)**2) + + ! We also need to extrapolate the ice viscosity to the + ! calving front, by weighting the components against the + ! magnitudes of the horizontal velocities + if (k < nz) then ! efvs is on staglevel + xext = linear_extrapolation(efvs(k,i+ud,j),efvs(k,i+2*ud,j)) + yext = linear_extrapolation(efvs(k,i,j+vd),efvs(k,i,j+2*vd)) + efvs_eff(k,i,j) = (uvel_unstag(k,i,j)*xext+vvel_unstag(k,i,j) & + *yext)/speed_unstag(k,i,j) + endif + enddo ! k + endif ! calving_front_mask + enddo ! i + enddo ! j + ! Finished unstaggered velocity computation !!! + endif ! Bassis & Ma damage + if (damage_manufactured) then allocate(vel_exp(nz,nx,ny)) allocate(adv_term(nz,nx,ny)) - allocate(uvel_unstag(nz,nx,ny)) allocate(hgl(ny)) + ! Set initial dummy values + vel_exp(:,:,:) = 0.0d0 + adv_term(:,:,:) = 0.0d0 + hgl(:) = 0.0d0 + ! Identify the grounding line thickness for manufactured damage on ! an ice tongue (assumes only one grounding line point per j) - do j = 2, ny-1 - do i = 2, nx-1 + do j = 3, ny-2 + do i = 3, nx-2 if (floating_mask(i,j) == 0 .and. grounding_line_mask(i,j) == 1) then hgl(j) = thck(i,j) endif @@ -680,8 +721,8 @@ subroutine glissade_calve_ice(which_calving, & enddo endif ! damage_manufactured - do j = 2, ny-1 - do i = 2, nx-1 + do j = 3, ny-2 + do i = 3, nx-2 if (floating_mask(i,j) == 1) then ! Prognose damage using a simple scheme based on the ! effective tensile stress, calving%tau_eff @@ -692,37 +733,38 @@ subroutine glissade_calve_ice(which_calving, & else if (damage_src == BASSIS_MA_DAMAGE_SRC) then ! Compute the ratio of the principal stresses (equivalent to the ratio of the ! principal strain rates) - alpha = calving%tau_eigen2(i,j) / calving%tau_eigen1(i,j) + alpha = tau2(i,j) / tau1(i,j) ! Find the max principal strain rate using Glen's Flow Law - eps_max(:,i,j) = calving%tau_eigen1(i,j) / (2.0d0 * efvs_eff(:,i,j)) + eps_max(:,i,j) = tau1(i,j) / (2.0d0 * efvs_eff(:,i,j)) ! Compute the effective flow law exponent nstar = 12.0d0*(1.0d0+alpha+alpha**2)/(4.0d0*(1.0d0+alpha+alpha**2)+6.0d0*alpha**2) ! Compute the hydrostatic to tensile stress ratio - szero = rhoi*(rhoo-rhoi)*grav*thck_eff(i,j)/(2.0d0*calving%tau_eigen1(i,j)*rhoo) + szero = rhoi*(rhoo-rhoi)*grav*thck(i,j)/(2.0d0*tau1(i,j)*rhoo) ! Prognose damage following Bassis & Ma 2015 - source(:,i,j) = nstar*(1.0d0-szero)*eps_max(:,i,j)-acab(i,j)/thck_eff(i,j) + source(:,i,j) = nstar*(1.0d0-szero)*eps_max(:,i,j) d_damage_dt(:,i,j) = source(:,i,j)*calving%damage(:,i,j) ! Compute the manufactured forcing term if requested; otherwise, set it to zero if (damage_manufactured) then + ! Compute the exponential + vel_exp(:,i,j) = exp(-uvel_unstag(:,i,j)*time/hgl(j)) + ! Compute the advection term adv_term(:,i,j) = 1.0d0 if (damage_advect) then adv_term(:,i,j) = adv_term(:,i,j)+eps_max(:,i,j)*time endif - - ! Linearly interpolate uvel to the unstaggered grid - uvel_unstag(:,i,j) = 0.25d0*abs(uvel(:,i,j)+uvel(:,i-1,j)+uvel(:,i,j-1) & - +uvel(:,i-1,j-1)) + adv_term(:,i,j) = adv_term(:,i,j)*uvel_unstag(:,i,j)*vel_exp(:,i,j)/hgl(j) + if (damage_advect) then + adv_term(:,i,j) = adv_term(:,i,j)-eps_max(:,i,j)*(1.0d0+vel_exp(:,i,j)) + endif ! Compute the manufactured forcing term - vel_exp(:,i,j) = exp(-uvel_unstag(:,i,j)*time/hgl(j)) - force(:,i,j) = -0.5d0*calving%damage_init(:,i,j)*(uvel_unstag(:,i,j)/hgl(j) & - *adv_term(:,i,j)*vel_exp(:,i,j) & + force(:,i,j) = -0.5d0*calving%damage_init(:,i,j)*thck(i,j)*(adv_term(:,i,j) & +source(:,i,j)*(1.0d0+vel_exp(:,i,j))) else force(:,i,j) = 0.0d0 @@ -733,21 +775,27 @@ subroutine glissade_calve_ice(which_calving, & endif ! damage_src + ! Convert crevasse depth back into damage + calving%damage(:,i,j) = calving%damage(:,i,j)/thck(i,j) + ! Constrain damage to the specified range of values if (damage_floor == ZERO_DAMAGE_FLOOR) then calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) elseif (damage_floor == NYE_DAMAGE_FLOOR) then ! Compute the Nye zero stress value - damage_nye = (2.0d0+alpha)*calving%tau_eigen1(i,j)/(grav*(rhoo-rhoi)*thck_eff(i,j)) + damage_nye = (2.0d0+alpha)*tau1(i,j)/(grav*(rhoo-rhoi)) calving%damage(:,i,j) = max(calving%damage(:,i,j), damage_nye) endif calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) - endif ! floating_mask enddo ! i + enddo ! j - do i = 2, nx-1 - ! Make sure that damage is zero where there isn't ice and where + call parallel_halo(calving%damage) + + do j = 3, ny-2 + do i = 3, nx-2 + ! Make sure that damage is zero where there isn't active ice and where ! the ice is grounded if (thck(i,j) <= thklim) then calving%damage(:,i,j) = 0.0d0 @@ -755,37 +803,50 @@ subroutine glissade_calve_ice(which_calving, & calving%damage(:,i,j) = 0.0d0 endif - ! When the Nye damage floor is not enforced, we need to sweep - ! across i a second time to get the damage along grounding lines - ! right. This is important because damage advected into the - ! domain is set to match the damage of what's already there, + ! Linearly extrapolate damage upstream to the grounding line + ! This is important because damage advected into the domain is + ! set to match the damage of what's already there, ! which for ice tongues happens to be the zero damage along the - ! grounding line. We want new damage to be as damaged as - ! existing ice, so we need the grounding line to also be as - ! damaged as the ice further downstream. - if (damage_floor == ZERO_DAMAGE_FLOOR) then - ! If i,j is the grounding line - if (floating_mask(i,j) == 0 .and. grounding_line_mask(i,j) == 1) then - ! If i+1,j is floating and ice exists there - if (floating_mask(i+1,j) == 1 .and. thck(i+1,j) > thklim) then - calving%damage(:,i,j) = calving%damage(:,i+1,j) + ! grounding line. + if (grounding_line_mask(i,j) == 1) then + do k = 1, nz + ! Store the velocity directions as matrix index modulations + ! e.g. if uvel is positive, set ud = -1 because ice will flow + ! from i-1 to i + ud = find_velocity_sign(uvel_unstag(k,i,j)) + vd = find_velocity_sign(vvel_unstag(k,i,j)) + + ! If the two points downstream along either horizontal axis + ! are floating, then linearly extrapolate upstream to + ! guess the grounding line damage + ! We use - instead of + here for ud/vd because we've + ! flipped the direction of extrapolation compared to what + ! we did for uvel/vvel/efvs_eff + if (floating_mask(i-ud,j) == 1 .and. floating_mask(i-2*ud,j) == 1) then + xext = linear_extrapolation(calving%damage(k,i-ud,j), & + calving%damage(k,i-2*ud,j)) + else + xext = 0.0d0 endif - endif - endif ! damage_floor + if (floating_mask(i,j-vd) == 1 .and. floating_mask(i,j-2*vd) == 1) then + yext = linear_extrapolation(calving%damage(k,i,j-vd), & + calving%damage(k,i,j-2*vd)) + else + yext = 0.0d0 + endif + + calving%damage(k,i,j) = (uvel_unstag(k,i,j)*xext+vvel_unstag(k,i,j) & + *yext)/speed_unstag(k,i,j) + enddo ! k + endif ! grounding_line_mask enddo ! i enddo ! j + call parallel_halo(calving%damage) + ! Compute the vertically integrated damage in each column. allocate(damage_column(nx,ny)) - damage_column(:,:) = 0.0d0 - - do j = 1, ny - do i = 1, nx - do k = 1, nz-1 - damage_column(i,j) = damage_column(i,j) + calving%damage(k,i,j) * (sigma(k+1) - sigma(k)) - enddo - enddo - enddo + damage_column(:,:) = calving%damage(1,:,:) if (which_ho_calving_front == HO_CALVING_FRONT_SUBGRID) then ! Convert damage in CF cells to a lateral calving rate (m/s). @@ -1563,7 +1624,8 @@ subroutine glissade_calve_ice(which_calving, & if (allocated(hgl)) deallocate(hgl) if (allocated(adv_term)) deallocate(adv_term) if (allocated(uvel_unstag)) deallocate(uvel_unstag) - if (allocated(thck_eff)) deallocate(thck_eff) + if (allocated(vvel_unstag)) deallocate(vvel_unstag) + if (allocated(speed_unstag)) deallocate(speed_unstag) if (allocated(efvs_eff)) deallocate(efvs_eff) end subroutine glissade_calve_ice @@ -2245,6 +2307,39 @@ recursive subroutine glissade_fill(nx, ny, & end subroutine glissade_fill +!------------------------------------------------------------------------------- + + function find_velocity_sign(vel_unstag) result(d) + ! Find the direction of the velocity at a given point, and report it as a + ! matrix index integer modulation (negative for positive velocities, + ! positive for negative velocities) + + real(dp), intent(in) :: vel_unstag ! Velocity on unstaggered grid (m/s) + integer :: d ! Output, sign of velocity as index modulation + + if (vel_unstag > 0.0d0) then + d = -1 + elseif (vel_unstag < 0.0d0) then + d = 1 + endif + end function + +!--------------------------------------------------------------------------- + + function linear_extrapolation(yup,yuptwo) result(ynew) + ! Extrapolate a quantity linearly in one dimension + + real(dp), intent(in) :: & + yup, & ! Function evaluated at xup + yuptwo ! Function evaluated at xuptwo + real(dp) :: ynew ! Output, extrapolated function value + + ! We're doing everything on a regular grid, so the absolute magnitudes of + ! the positions don't matter (only the relative magnitudes) + ynew = yuptwo + 2.0d0*(yup-yuptwo) + + end function + !--------------------------------------------------------------------------- end module glissade_calving diff --git a/tests/tongue/doc/README.tex b/tests/tongue/doc/README.tex index 55c67889..6873a39b 100644 --- a/tests/tongue/doc/README.tex +++ b/tests/tongue/doc/README.tex @@ -85,7 +85,7 @@ \section{Ice Tongues} \section{Damage} -As stated above, this test case is specifically designed to verify damage $r$ as computed by Eq.~27 in \citet{Bassis-Ma}: +As stated above, this test case is specifically designed to verify damage $r$ as governed by Eq.~27 in \citet{Bassis-Ma}: \begin{equation} \frac{\mathrm{d}r}{\mathrm{d}t} = \left[n^* \left(1-S_0\right) \dot{\epsilon}_{xx} + \frac{\dot{m}}{H}\right] r \label{eq:bassis-ma-damage} @@ -102,7 +102,19 @@ \section{Damage} \end{equation} where $\tau_{xx}$ is the largest principal deviatoric stress. -Eq.~\ref{eq:bassis-ma-damage} takes in an initial damage field, evolving it in time and advecting it with the flow field. Depending on the value of the \texttt{damage\_floor} configuration option, the initial value $r_0$ is specified as either a uniform value (set as the \texttt{const\_damage} variable in the python script) or as the Nye's zero stress value, computed via \citep{Jezek,Nick-et-al,Nye} +CISM solves Eq.~\ref{eq:bassis-ma-damage} by separating the equation into an advective component and an advection-less component. The former is handled by CISM's built-in incremental remapping scheme, which assumes an equation of the form +\begin{equation} + \frac{\delta \left( hr \right)}{\delta t} + \nabla\cdot\left( hr\mathbf{U} \right) = 0 + \label{eq:advection} +\end{equation} + +In order to make use of this scheme, we need to reformulate the damage evolution equation in terms of the crevasse depth $\tilde{r} \equiv rH$: +\begin{equation} + \frac{\mathrm{d}\tilde{r}}{\mathrm{d}t} = n^* \left(1-S_0\right) \dot{\epsilon}_{xx} \tilde{r} + \label{eq:bassis-ma-crevasse-depth} +\end{equation} + +Eq.~\ref{eq:bassis-ma-crevasse-depth} takes in an initial damage field, evolving it in time and advecting it with the flow field. Depending on the value of the \texttt{damage\_floor} configuration option, the initial value $r_0$ is specified as either a uniform value (set as the \texttt{const\_damage} variable in the python script) or as the Nye's zero stress value, computed via \citep{Jezek,Nick-et-al,Nye} \begin{equation} r_0 = \frac{\rho_i}{\rho_w-\rho_i} \left[\frac{\left(2+\alpha\right) \tau_{xx}}{\rho_i gH}\right]. \label{eq:nye-damage} @@ -110,17 +122,17 @@ \section{Damage} Verification of the damage solution can be performed by toggling the \texttt{damage\_} \texttt{manufactured} configuration option. When activated, this option adjusts Eq.~\ref{eq:bassis-ma-damage} by adding a forcing term $f$: \begin{equation} - \frac{\mathrm{d}r}{\mathrm{d}t} = \left[n^* \left(1-S_0\right) \dot{\epsilon}_{xx} + \frac{\dot{m}}{H}\right] r + f + \frac{\mathrm{d}\tilde{r}}{\mathrm{d}t} = n^* \left(1-S_0\right) \dot{\epsilon}_{xx} \tilde{r} + f \label{eq:bassis-ma-damage-forced} \end{equation} This forcing term is defined as \begin{equation} - f = -\frac{r_0}{2} \left\{\frac{U}{H_0} \left(1+\dot{\epsilon}_{xx}t\right) e^{-Ut/H_0} + \left[n^*\left(1-S_0\right)\dot{\epsilon}_{xx}+\frac{\dot{m}}{H}\right] \left(1+e^{-Ut/H_0}\right)\right\}, + f = -\frac{\tilde{r}_0}{2} \left\{\frac{U}{H_0} \left(1+\dot{\epsilon}_{xx}t\right) e^{-Ut/H_0} + \dot{\epsilon}_{xx} \left[ n^*\left(1-S_0\right)-1 \right] \left(1+e^{-Ut/H_0}\right)\right\}, \label{eq:damage-forcing} \end{equation} such that Eq.~\ref{eq:bassis-ma-damage-forced} produces a damage field given by \begin{equation} - r = \frac{r_0}{2} \left(1+e^{-Ut/H_0}\right). + \tilde{r} = \frac{\tilde{r}_0}{2} \left(1+e^{-Ut/H_0}\right). \label{eq:manufactured-damage} \end{equation} From 79c4e307f9ab802cd5b160d5577c812769c3e104 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Wed, 14 Aug 2019 14:30:15 -0400 Subject: [PATCH 7/9] Bugfixes - Fixed a bug that caused a segmentation fault when the Bassis and Ma damage source term was deactivated - Small tweak to the ice tongue test case documentation file for consistency --- libglissade/glissade_calving.F90 | 50 +++++++++++++++++--------------- tests/tongue/doc/README.tex | 2 +- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index f74dcb3d..af5374fa 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -621,45 +621,47 @@ subroutine glissade_calve_ice(which_calving, & ! Note: The damage is formally a 3D field, which makes it easier to advect, even though ! (in the current scheme) the damage source term is uniform in each column. + ! Allocate relevant data matrices + allocate(uvel_unstag(nz,nx,ny)) + allocate(vvel_unstag(nz,nx,ny)) + allocate(speed_unstag(nz,nx,ny)) + ! Set initial dummy values + uvel_unstag(:,:,:) = 0.0d0 + vvel_unstag(:,:,:) = 0.0d0 + speed_unstag(:,:,:) = 0.0d0 d_damage_dt(:,:,:) = 0.0d0 - ! Allocate relevant data matrices + ! Construct unstaggered velocity matrices + do j = 3, ny-2 + do i = 3, nx-2 + ! Linearly interpolate velocities to the unstaggered grid + uvel_unstag(:,i,j) = 0.25d0*abs(uvel(:,i,j)+uvel(:,i-1,j)+uvel(:,i,j-1) & + +uvel(:,i-1,j-1)) + vvel_unstag(:,i,j) = 0.25d0*abs(vvel(:,i,j)+vvel(:,i-1,j)+vvel(:,i,j-1) & + +vvel(:,i-1,j-1)) + + ! Compute the magnitude of the horizontal velocity vectors + speed_unstag(:,i,j) = SQRT(uvel_unstag(:,i,j)**2+vvel_unstag(:,i,j)**2) + enddo ! i + enddo ! j + + call parallel_halo(uvel_unstag) + call parallel_halo(vvel_unstag) + if (damage_src == BASSIS_MA_DAMAGE_SRC) then - allocate(uvel_unstag(nz,nx,ny)) - allocate(vvel_unstag(nz,nx,ny)) - allocate(speed_unstag(nz,nx,ny)) + ! Allocate relevant data matrices allocate(efvs_eff(nz-1,nx,ny)) allocate(eps_max(nz,nx,ny)) allocate(source(nz,nx,ny)) allocate(force(nz,nx,ny)) ! Set initial dummy values - uvel_unstag(:,:,:) = 0.0d0 - vvel_unstag(:,:,:) = 0.0d0 - speed_unstag(:,:,:) = 0.0d0 efvs_eff(:,:,:) = efvs(:,:,:) eps_max(:,:,:) = 0.0d0 source(:,:,:) = 0.0d0 force(:,:,:) = 0.0d0 - ! Construct unstaggered velocity matrices - do j = 3, ny-2 - do i = 3, nx-2 - ! Linearly interpolate velocities to the unstaggered grid - uvel_unstag(:,i,j) = 0.25d0*abs(uvel(:,i,j)+uvel(:,i-1,j)+uvel(:,i,j-1) & - +uvel(:,i-1,j-1)) - vvel_unstag(:,i,j) = 0.25d0*abs(vvel(:,i,j)+vvel(:,i-1,j)+vvel(:,i,j-1) & - +vvel(:,i-1,j-1)) - - ! Compute the magnitude of the horizontal velocity vectors - speed_unstag(:,i,j) = SQRT(uvel_unstag(:,i,j)**2+vvel_unstag(:,i,j)**2) - enddo ! i - enddo ! j - - call parallel_halo(uvel_unstag) - call parallel_halo(vvel_unstag) - ! We need to restart the loops here to have the matrices fully computed for extrapolation ! We want to only sweep over points that exist (so we need a buffer of at ! least two in the halo). diff --git a/tests/tongue/doc/README.tex b/tests/tongue/doc/README.tex index 6873a39b..436d8cfb 100644 --- a/tests/tongue/doc/README.tex +++ b/tests/tongue/doc/README.tex @@ -104,7 +104,7 @@ \section{Damage} CISM solves Eq.~\ref{eq:bassis-ma-damage} by separating the equation into an advective component and an advection-less component. The former is handled by CISM's built-in incremental remapping scheme, which assumes an equation of the form \begin{equation} - \frac{\delta \left( hr \right)}{\delta t} + \nabla\cdot\left( hr\mathbf{U} \right) = 0 + \frac{\delta \left( rH \right)}{\delta t} + \nabla\cdot\left( rH\mathbf{U} \right) = 0 \label{eq:advection} \end{equation} From 50d2f2c57a47d448068b74d2bf8c0f490b2d4588 Mon Sep 17 00:00:00 2001 From: morgwhit Date: Wed, 2 Oct 2019 12:23:00 -0400 Subject: [PATCH 8/9] Commit reversion + more bugfixes Reverting changes from the physics/discretization commit 0b33606, such that damage is no longer converted to crevasse depth during the course of the damage evolution routine. Also, this commit extends the unstaggered velocity extrapolation procedure to occur at the grounding line in addition to the calving front, where there are not enough staggered grid points for the velocity to be defined. Some of this procedural extension involves minor bugfixes in identifying the correct points within the domain. --- libglissade/glissade.F90 | 7 ------ libglissade/glissade_calving.F90 | 43 ++++++++++++++------------------ 2 files changed, 19 insertions(+), 31 deletions(-) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index a2605633..063521e3 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -2055,13 +2055,6 @@ subroutine glissade_calving_solve(model, init_calving) ! Pass in thck, topg, etc. with units of meters. ! ------------------------------------------------------------------------ - ! Convert damage to crevasse depth - if (model%options%whichcalving == CALVING_DAMAGE) then - do k = 1, size(model%numerics%stagsigma) - model%calving%damage(k,:,:) = model%calving%damage(k,:,:)*model%geometry%thck_old(:,:)*thk0 - enddo - endif - thck_unscaled(:,:) = model%geometry%thck_old(:,:)*thk0 call glissade_calve_ice(model%options%whichcalving, & diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index af5374fa..d0be983e 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -668,15 +668,21 @@ subroutine glissade_calve_ice(which_calving, & do j = 3, ny-2 do i = 3, nx-2 ! To do this right we also need to extrapolate the velocities to - ! the calving front, since there are too few points on the - ! staggered grid to fully interpolate - if (calving_front_mask(i,j) == 1) then + ! the calving front and the grounding line, since there are too few + ! points on the staggered grid to fully interpolate + if (calving_front_mask(i,j) == 1 .or. (grounding_line_mask(i,j) == 1 & + .and. floating_mask(i,j) == 0)) then do k = 1, nz ! Store the velocity directions as matrix index modulations ! e.g. if uvel is positive, set ud = -1 because ice will flow ! from i-1 to i - ud = find_velocity_sign(uvel_unstag(k,i,j)) - vd = find_velocity_sign(vvel_unstag(k,i,j)) + if (calving_front_mask(i,j) == 1) then + ud = find_velocity_sign(uvel_unstag(k,i,j),1) + vd = find_velocity_sign(vvel_unstag(k,i,j),1) + elseif (grounding_line_mask(i,j) == 1) then + ud = find_velocity_sign(uvel_unstag(k,i,j),-1) + vd = find_velocity_sign(vvel_unstag(k,i,j),-1) + endif ! calving front vs grounding line uvel_unstag(k,i,j) = linear_extrapolation(uvel_unstag(k,i+ud,j), & uvel_unstag(k,i+2*ud,j)) @@ -689,7 +695,7 @@ subroutine glissade_calve_ice(which_calving, & ! We also need to extrapolate the ice viscosity to the ! calving front, by weighting the components against the ! magnitudes of the horizontal velocities - if (k < nz) then ! efvs is on staglevel + if (calving_front_mask(i,j) == 1 .and. k < nz) then ! efvs is on staglevel xext = linear_extrapolation(efvs(k,i+ud,j),efvs(k,i+2*ud,j)) yext = linear_extrapolation(efvs(k,i,j+vd),efvs(k,i,j+2*vd)) efvs_eff(k,i,j) = (uvel_unstag(k,i,j)*xext+vvel_unstag(k,i,j) & @@ -747,7 +753,7 @@ subroutine glissade_calve_ice(which_calving, & szero = rhoi*(rhoo-rhoi)*grav*thck(i,j)/(2.0d0*tau1(i,j)*rhoo) ! Prognose damage following Bassis & Ma 2015 - source(:,i,j) = nstar*(1.0d0-szero)*eps_max(:,i,j) + source(:,i,j) = nstar*(1.0d0-szero)*eps_max(:,i,j)-acab(i,j)/thck(i,j) d_damage_dt(:,i,j) = source(:,i,j)*calving%damage(:,i,j) ! Compute the manufactured forcing term if requested; otherwise, set it to zero @@ -761,12 +767,9 @@ subroutine glissade_calve_ice(which_calving, & adv_term(:,i,j) = adv_term(:,i,j)+eps_max(:,i,j)*time endif adv_term(:,i,j) = adv_term(:,i,j)*uvel_unstag(:,i,j)*vel_exp(:,i,j)/hgl(j) - if (damage_advect) then - adv_term(:,i,j) = adv_term(:,i,j)-eps_max(:,i,j)*(1.0d0+vel_exp(:,i,j)) - endif ! Compute the manufactured forcing term - force(:,i,j) = -0.5d0*calving%damage_init(:,i,j)*thck(i,j)*(adv_term(:,i,j) & + force(:,i,j) = -0.5d0*calving%damage_init(:,i,j)*(adv_term(:,i,j) & +source(:,i,j)*(1.0d0+vel_exp(:,i,j))) else force(:,i,j) = 0.0d0 @@ -777,9 +780,6 @@ subroutine glissade_calve_ice(which_calving, & endif ! damage_src - ! Convert crevasse depth back into damage - calving%damage(:,i,j) = calving%damage(:,i,j)/thck(i,j) - ! Constrain damage to the specified range of values if (damage_floor == ZERO_DAMAGE_FLOOR) then calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) @@ -810,14 +810,8 @@ subroutine glissade_calve_ice(which_calving, & ! set to match the damage of what's already there, ! which for ice tongues happens to be the zero damage along the ! grounding line. - if (grounding_line_mask(i,j) == 1) then + if (grounding_line_mask(i,j) == 1 .and. floating_mask(i,j) == 0) then do k = 1, nz - ! Store the velocity directions as matrix index modulations - ! e.g. if uvel is positive, set ud = -1 because ice will flow - ! from i-1 to i - ud = find_velocity_sign(uvel_unstag(k,i,j)) - vd = find_velocity_sign(vvel_unstag(k,i,j)) - ! If the two points downstream along either horizontal axis ! are floating, then linearly extrapolate upstream to ! guess the grounding line damage @@ -2311,18 +2305,19 @@ end subroutine glissade_fill !------------------------------------------------------------------------------- - function find_velocity_sign(vel_unstag) result(d) + function find_velocity_sign(vel_unstag,dir) result(d) ! Find the direction of the velocity at a given point, and report it as a ! matrix index integer modulation (negative for positive velocities, ! positive for negative velocities) real(dp), intent(in) :: vel_unstag ! Velocity on unstaggered grid (m/s) + integer, intent(in) :: dir ! Direction of extrapolation integer :: d ! Output, sign of velocity as index modulation if (vel_unstag > 0.0d0) then - d = -1 + d = -dir elseif (vel_unstag < 0.0d0) then - d = 1 + d = dir endif end function From 9f9c1aa0b54b28b39ce371197f8f33367ccabd6b Mon Sep 17 00:00:00 2001 From: morgwhit Date: Tue, 15 Oct 2019 20:13:24 -0400 Subject: [PATCH 9/9] Further reversion + bugfixes Completed the commit reversion by fixing the Nye zero stress criterion for damage. Also fixed an issue with the thickness field never evolving, and an issue with the way the bedrock topography was initialized by ice_tongue.py for 2D ice shelves. --- libglissade/glissade.F90 | 2 +- libglissade/glissade_calving.F90 | 2 +- tests/tongue/ice_tongue.py | 5 ++++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 063521e3..993f8e8a 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -2055,7 +2055,7 @@ subroutine glissade_calving_solve(model, init_calving) ! Pass in thck, topg, etc. with units of meters. ! ------------------------------------------------------------------------ - thck_unscaled(:,:) = model%geometry%thck_old(:,:)*thk0 + thck_unscaled(:,:) = model%geometry%thck(:,:)*thk0 call glissade_calve_ice(model%options%whichcalving, & model%options%calving_domain, & diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index d0be983e..4badc1b3 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -785,7 +785,7 @@ subroutine glissade_calve_ice(which_calving, & calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) elseif (damage_floor == NYE_DAMAGE_FLOOR) then ! Compute the Nye zero stress value - damage_nye = (2.0d0+alpha)*tau1(i,j)/(grav*(rhoo-rhoi)) + damage_nye = (2.0d0+alpha)*tau1(i,j)/(grav*(rhoo-rhoi)*thck(i,j)) calving%damage(:,i,j) = max(calving%damage(:,i,j), damage_nye) endif calving%damage(:,i,j) = min(calving%damage(:,i,j), 1.0d0) diff --git a/tests/tongue/ice_tongue.py b/tests/tongue/ice_tongue.py index 55baa7ac..fbbe98e3 100755 --- a/tests/tongue/ice_tongue.py +++ b/tests/tongue/ice_tongue.py @@ -387,7 +387,10 @@ def enforce_confinement(s,field,val,fieldstr,xend): elif fieldstr == 'uvel': s.uvel_extend[0] = field elif fieldstr == 'topg': - s.topg = field + if np.shape(s.topg) == np.shape(field): # Restart files + s.topg = field + else: # Otherwise + s.topg[0] = field elif fieldstr == 'damage': s.damage[0] = field