diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 7acbfd24..913b07e6 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,17 @@ 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) 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 417e79ce..669b89d1 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 @@ -1135,6 +1165,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) @@ -2188,6 +2219,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) @@ -2330,9 +2362,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 @@ -2687,6 +2725,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)) & @@ -2821,6 +2861,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 ce1c2557..1255c836 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 @@ -646,6 +654,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 8c9f15e5..993f8e8a 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 + 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 @@ -2024,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 @@ -2042,25 +2057,34 @@ 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 - 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%velocity%vvel*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 fa7c89c5..4badc1b3 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -194,16 +194,23 @@ 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 + uvel, vvel, & ! m/s + acab, & ! m/s thck, relx, & ! m topg, eus) ! m @@ -224,11 +231,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 @@ -261,10 +277,15 @@ 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) :: 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) real(dp), dimension(:,:), intent(in) :: topg !> present bedrock topography (m) @@ -276,8 +297,21 @@ 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 :: & + eps_max, & ! maximum principal strain rate + 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) + 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 @@ -287,12 +321,16 @@ subroutine glissade_calve_ice(which_calving, & 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 @@ -318,9 +356,14 @@ 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 + 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 + 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, @@ -375,11 +418,13 @@ 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)) 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)) @@ -446,6 +491,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) @@ -495,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) @@ -508,66 +559,46 @@ subroutine glissade_calve_ice(which_calving, & call parallel_halo(tau1) call parallel_halo(tau2) - ! In inactive calving-front cells where both eigenvalues are zero (because a cell is dynamically inactive), - ! extrapolate nonzero values in upstream cells. - ! 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 - tau1(i,j) = tau1(ii,jj) - tau2(i,j) = tau2(ii,jj) - endif - enddo - enddo - endif ! tau1 = tau2 = 0 - 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. + 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 @@ -585,72 +616,272 @@ 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. - do j = 2, ny-1 - do i = 2, nx-1 + ! 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 + + ! 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 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 + efvs_eff(:,:,:) = efvs(:,:,:) + eps_max(:,:,:) = 0.0d0 + source(:,:,:) = 0.0d0 + force(:,:,:) = 0.0d0 + + ! 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 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 + 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)) + 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 (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) & + *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(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 = 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 + enddo + enddo + endif ! damage_manufactured + + do j = 3, ny-2 + do i = 3, nx-2 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 + ! 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 = tau2(i,j) / tau1(i,j) + + ! Find the max principal strain rate using Glen's Flow Law + 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(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(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 + adv_term(:,i,j) = adv_term(:,i,j)*uvel_unstag(:,i,j)*vel_exp(:,i,j)/hgl(j) + + ! Compute the manufactured forcing term + 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 + endif + + 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 + + endif ! damage_src + + ! 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)*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) - calving%damage(:,i,j) = max(calving%damage(:,i,j), 0.0d0) - else ! set damage to zero for grounded ice + endif ! floating_mask + enddo ! i + enddo ! j + + 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 + elseif (floating_mask(i,j) == 0) then calving%damage(:,i,j) = 0.0d0 endif - enddo - enddo + + ! 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. + if (grounding_line_mask(i,j) == 1 .and. floating_mask(i,j) == 0) then + do k = 1, nz + ! 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 + 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)) + 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). + ! 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 - - ! 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 - 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 @@ -1370,16 +1601,28 @@ 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) 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) + 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(vvel_unstag)) deallocate(vvel_unstag) + if (allocated(speed_unstag)) deallocate(speed_unstag) + if (allocated(efvs_eff)) deallocate(efvs_eff) end subroutine glissade_calve_ice @@ -2060,6 +2303,40 @@ recursive subroutine glissade_fill(nx, ny, & end subroutine glissade_fill +!------------------------------------------------------------------------------- + + 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 = -dir + elseif (vel_unstag < 0.0d0) then + d = dir + 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/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) diff --git a/tests/tongue/doc/README.tex b/tests/tongue/doc/README.tex new file mode 100644 index 00000000..436d8cfb --- /dev/null +++ b/tests/tongue/doc/README.tex @@ -0,0 +1,142 @@ +\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$ & $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$ \\ + $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 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} +\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. + +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( rH \right)}{\delta t} + \nabla\cdot\left( rH\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} +\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}\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{\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} + \tilde{r} = \frac{\tilde{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..fbbe98e3 --- /dev/null +++ b/tests/tongue/ice_tongue.py @@ -0,0 +1,476 @@ +#!/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 + s.topg[0] = compute_topg(s,c,s,s.topg[0],s.thk[0]) + + # Enforce lateral confinement + if s.conflen > 0.: # If laterally confined + 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 + 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] + +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 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 + 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 + 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': + 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 + +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) + r.topg = compute_topg(r,c,g,r.topg,r.thk) + 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'