Skip to content

UM Damage Evolution #28

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
morgwhitcomb wants to merge 10 commits into
base: master
Choose a base branch
from
Open

UM Damage Evolution #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
d9ab613
Implementing UM damage evolution
morgwhitcomb Feb 20, 2019
6b3c8c8
Damage input options, new diagnostic variables, and calving
morgwhitcomb Apr 5, 2019
7dbba49
Ice tongue test case
morgwhitcomb Apr 8, 2019
98d5e93
Minor bugfixes
morgwhitcomb Apr 24, 2019
c825590
Merge branch 'master' into morgwhit/damage
morgwhitcomb Apr 24, 2019
0e666fa
Slight doc fix
morgwhitcomb Apr 24, 2019
0b33606
Significant bug fixes in the physics/discretization + support for sub…
morgwhitcomb Aug 13, 2019
79c4e30
Bugfixes
morgwhitcomb Aug 14, 2019
50d2f2c
Commit reversion + more bugfixes
morgwhitcomb Oct 2, 2019
9f9c1aa
Further reversion + bugfixes
morgwhitcomb Oct 16, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 60 additions & 1 deletion libglide/glide_setup.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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'/)
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
42 changes: 42 additions & 0 deletions libglide/glide_types.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)) &
Expand Down Expand Up @@ -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

Expand Down
17 changes: 17 additions & 0 deletions libglide/glide_vars.def
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
64 changes: 44 additions & 20 deletions libglissade/glissade.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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. &
Expand Down Expand Up @@ -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

!------------------------------------------------------------------------
Expand Down Expand Up @@ -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

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

Expand All @@ -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
Expand Down
Loading