diff --git a/libglide/felix_dycore_interface.F90 b/libglide/felix_dycore_interface.F90 index f702d9a1..ce1ed0dc 100644 --- a/libglide/felix_dycore_interface.F90 +++ b/libglide/felix_dycore_interface.F90 @@ -146,7 +146,7 @@ end subroutine felix_velo_init subroutine felix_velo_driver(model) use glimmer_global, only : dp - use glimmer_physcon, only: gn, scyr + use glimmer_physcon, only: scyr use glimmer_paramets, only: thk0, len0, vel0, vis0 use glimmer_log use glide_types diff --git a/libglide/glide.F90 b/libglide/glide.F90 index 9683d6a8..5f7d6174 100644 --- a/libglide/glide.F90 +++ b/libglide/glide.F90 @@ -602,8 +602,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice) call calcbwat(model, & model%options%whichbwat, & model%basal_melt%bmlt, & - model%temper%bwat, & - model%temper%bwatflx, & + model%basal_hydro%bwat, & + model%basal_hydro%bwatflx, & model%geometry%thck, & model%geometry%topg, & model%temper%temp(model%general%upn,:,:), & @@ -612,8 +612,8 @@ subroutine glide_init_state_diagnostic(model, evolve_ice) ! This call is redundant for now, but is needed if the call to calcbwat is removed - call stagvarb(model%temper%bwat, & - model%temper%stagbwat ,& + call stagvarb(model%basal_hydro%bwat, & + model%basal_hydro%stagbwat ,& model%general%ewn, & model%general%nsn) @@ -867,8 +867,8 @@ subroutine glide_tstep_p1(model,time) call calcbwat(model, & model%options%whichbwat, & model%basal_melt%bmlt, & - model%temper%bwat, & - model%temper%bwatflx, & + model%basal_hydro%bwat, & + model%basal_hydro%bwatflx, & model%geometry%thck, & model%geometry%topg, & model%temper%temp(model%general%upn,:,:), & diff --git a/libglide/glide_bwater.F90 b/libglide/glide_bwater.F90 index 0126d1b7..35efa72b 100644 --- a/libglide/glide_bwater.F90 +++ b/libglide/glide_bwater.F90 @@ -195,8 +195,8 @@ subroutine calcbwat(model, which, bmlt, bwat, bwatflx, thck, topg, btem, floater end select ! now also calculate basal water in velocity (staggered) coord system - call stagvarb(model%temper%bwat, & - model%temper%stagbwat ,& + call stagvarb(model%basal_hydro%bwat, & + model%basal_hydro%stagbwat ,& model%general%ewn, & model%general%nsn) @@ -429,7 +429,7 @@ subroutine pressure_wphi(thck,topg,N,wphi,thicklim,floater) !> Compute the pressure wphi at the base of the ice sheet according to !> ice overburden plus bed height minus effective pressure. !> - !> whpi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g) + !> wphi/(rhow*g) = topg + bwat * rhoi / rhow * thick - N / (rhow * g) use glimmer_physcon, only : rhoi,rhow,grav implicit none diff --git a/libglide/glide_diagnostics.F90 b/libglide/glide_diagnostics.F90 index 8da14b2b..f8bb3718 100644 --- a/libglide/glide_diagnostics.F90 +++ b/libglide/glide_diagnostics.F90 @@ -200,9 +200,12 @@ subroutine glide_write_diag (model, time) max_thck, max_thck_global, & ! max ice thickness (m) max_temp, max_temp_global, & ! max ice temperature (deg C) min_temp, min_temp_global, & ! min ice temperature (deg C) - max_bmlt, max_bmlt_global, & ! max basal melt rate (m/yr) - max_spd_sfc, max_spd_sfc_global, & ! max surface ice speed (m/yr) - max_spd_bas, max_spd_bas_global, & ! max basal ice speed (m/yr) + max_bmlt, & ! max basal melt rate (m/yr) + max_bmlt_global, & + max_bmlt_ground, & ! max basal melt rate, grounded ice (m/yr) + max_bmlt_ground_global, & + max_spd_sfc, max_spd_sfc_global, & ! max surface ice speed (m/yr) + max_spd_bas, max_spd_bas_global, & ! max basal ice speed (m/yr) spd, & ! speed thck_diag, usrf_diag, & ! local column diagnostics topg_diag, relx_diag, & @@ -768,7 +771,8 @@ subroutine glide_write_diag (model, time) min_temp_global, imin_global, jmin_global, kmin_global call write_log(trim(message), type = GM_DIAGNOSTIC) - ! max basal melt rate + ! max applied basal melt rate + ! Usually, this will be for floating ice, if floating ice is present imax = 0 jmax = 0 max_bmlt = unphys_val @@ -791,11 +795,39 @@ subroutine glide_write_diag (model, time) ! Write to diagnostics only if nonzero if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then - write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/yr), i, j ', & + write(message,'(a25,f24.16,2i6)') 'Max bmlt (m/y), i, j ', & max_bmlt_global*thk0*scyr/tim0, imax_global, jmax_global call write_log(trim(message), type = GM_DIAGNOSTIC) endif + ! max basal melt rate for grounded ice + imax = 0 + jmax = 0 + max_bmlt_ground = unphys_val + + do j = lhalo+1, nsn-uhalo + do i = lhalo+1, ewn-uhalo + if (model%basal_melt%bmlt_ground(i,j) > max_bmlt_ground) then + max_bmlt_ground = model%basal_melt%bmlt_ground(i,j) + imax = i + jmax = j + endif + enddo + enddo + + call parallel_reduce_maxloc(xin=max_bmlt_ground, xout=max_bmlt_ground_global, xprocout=procnum) + call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel) + call broadcast(imax_global, proc = procnum) + call broadcast(jmax_global, proc = procnum) + + ! Write to diagnostics only if nonzero + + if (abs(max_bmlt_global*thk0*scyr/tim0) > eps) then + write(message,'(a25,f24.16,2i6)') 'Max bmlt grnd (m/y), i, j', & + max_bmlt_ground_global*thk0*scyr/tim0, imax_global, jmax_global + call write_log(trim(message), type = GM_DIAGNOSTIC) + endif + ! max surface speed imax = 0 jmax = 0 @@ -818,7 +850,7 @@ subroutine glide_write_diag (model, time) call broadcast(imax_global, proc = procnum) call broadcast(jmax_global, proc = procnum) - write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/yr), i, j ', & + write(message,'(a25,f24.16,2i6)') 'Max sfc spd (m/y), i, j ', & max_spd_sfc_global*vel0*scyr, imax_global, jmax_global call write_log(trim(message), type = GM_DIAGNOSTIC) @@ -843,7 +875,7 @@ subroutine glide_write_diag (model, time) call parallel_globalindex(imax, jmax, imax_global, jmax_global, parallel) call broadcast(imax_global, proc = procnum) call broadcast(jmax_global, proc = procnum) - write(message,'(a25,f24.16,2i6)') 'Max base spd (m/yr), i, j', & + write(message,'(a25,f24.16,2i6)') 'Max base spd (m/y), i, j ', & max_spd_bas_global*vel0*scyr, imax_global, jmax_global call write_log(trim(message), type = GM_DIAGNOSTIC) @@ -884,7 +916,7 @@ subroutine glide_write_diag (model, time) artm_diag = model%climate%artm_corrected(i,j) ! artm_corrected = artm + artm_anomaly acab_diag = model%climate%acab_applied(i,j) * thk0*scyr/tim0 bmlt_diag = model%basal_melt%bmlt_applied(i,j) * thk0*scyr/tim0 - bwat_diag = model%temper%bwat(i,j) * thk0 + bwat_diag = model%basal_hydro%bwat(i,j) * thk0 bheatflx_diag = model%temper%bheatflx(i,j) temp_diag(:) = model%temper%temp(1:upn,i,j) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 3b967d05..7bfe06de 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -160,7 +160,7 @@ end subroutine glide_printconfig subroutine glide_scale_params(model) !> scale parameters use glide_types - use glimmer_physcon, only: scyr, gn + use glimmer_physcon, only: scyr use glimmer_paramets, only: thk0, tim0, len0, vel0, vis0, acc0, tau0 implicit none @@ -745,6 +745,7 @@ subroutine handle_options(section, model) call GetValue(section,'cull_calving_front', model%options%cull_calving_front) call GetValue(section,'adjust_input_thickness', model%options%adjust_input_thickness) call GetValue(section,'smooth_input_topography', model%options%smooth_input_topography) + call GetValue(section,'smooth_input_usrf', model%options%smooth_input_usrf) call GetValue(section,'adjust_input_topography', model%options%adjust_input_topography) call GetValue(section,'read_lat_lon',model%options%read_lat_lon) call GetValue(section,'dm_dt_diag',model%options%dm_dt_diag) @@ -786,6 +787,7 @@ subroutine handle_ho_options(section, model) call GetValue(section, 'which_ho_bmlt_inversion', model%options%which_ho_bmlt_inversion) call GetValue(section, 'which_ho_bmlt_basin_inversion', model%options%which_ho_bmlt_basin_inversion) call GetValue(section, 'which_ho_bwat', model%options%which_ho_bwat) + call GetValue(section, 'ho_flux_routing_scheme', model%options%ho_flux_routing_scheme) call GetValue(section, 'which_ho_effecpress', model%options%which_ho_effecpress) call GetValue(section, 'which_ho_resid', model%options%which_ho_resid) call GetValue(section, 'which_ho_nonlinear', model%options%which_ho_nonlinear) @@ -883,11 +885,10 @@ subroutine print_options(model) 'advective-diffusive balance ',& 'temp from external file ' /) - character(len=*), dimension(0:3), parameter :: flow_law = (/ & - 'const 1e-16 Pa^-n a^-1 ', & + character(len=*), dimension(0:2), parameter :: flow_law = (/ & + 'uniform factor flwa ', & 'Paterson and Budd (T = -5 C)', & - 'Paterson and Budd ', & - 'read flwa/flwastag from file' /) + 'Paterson and Budd ' /) !TODO - Rename slip_coeff to which_btrc? character(len=*), dimension(0:5), parameter :: slip_coeff = (/ & @@ -959,10 +960,11 @@ subroutine print_options(model) 'artm and d(artm)/dz input as function of (x,y)', & 'artm input as function of (x,y,z) ' /) - character(len=*), dimension(0:2), parameter :: overwrite_acab = (/ & + character(len=*), dimension(0:3), parameter :: overwrite_acab = (/ & 'do not overwrite acab anywhere ', & 'overwrite acab where input acab = 0 ', & - 'overwrite acab where input thck <= minthck' /) + 'overwrite acab where input thck <= minthck', & + 'overwrite acab based on input mask ' /) ! NOTE: Set gthf = 1 in the config file to read the geothermal heat flux from an input file. ! Otherwise it will be overwritten, even if the 'bheatflx' field is present. @@ -1063,17 +1065,24 @@ subroutine print_options(model) 'invert for basin-based basal melting parameters ', & 'apply basin basal melting parameters from earlier inversion' /) - character(len=*), dimension(0:2), parameter :: ho_whichbwat = (/ & + character(len=*), dimension(0:3), parameter :: ho_whichbwat = (/ & 'zero basal water depth ', & 'constant basal water depth ', & - 'basal water depth computed from local till model' /) + 'basal water depth computed from local till model', & + 'steady-state water routing with flux calculation' /) - character(len=*), dimension(0:4), parameter :: ho_whicheffecpress = (/ & + character(len=*), dimension(0:2), parameter :: ho_flux_routing_scheme = (/ & + 'D8; route flux to lowest-elevation neighbor ', & + 'Dinf; route flux to two lower-elevation neighbors', & + 'FD8; route flux to all lower-elevation neighbors ' /) + + character(len=*), dimension(0:5), parameter :: ho_whicheffecpress = (/ & 'full overburden pressure ', & 'reduced effecpress near pressure melting point ', & 'reduced effecpress where there is melting at the bed ', & 'reduced effecpress where bed is connected to ocean ', & - 'reduced effecpress with increasing basal water '/) + 'reduced effecpress with increasing basal water (B/vP)', & + 'reduced effecpress with increasing basal water (ramp)'/) character(len=*), dimension(0:1), parameter :: which_ho_nonlinear = (/ & 'use standard Picard iteration ', & @@ -1234,7 +1243,7 @@ subroutine print_options(model) end if if (tasks > 1 .and. model%options%whichbwat==BWATER_FLUX) then - call write_log('Error, flux-based basal water option not supported for more than one processor', GM_FATAL) + call write_log('Error, flux-based basal water option not yet supported for more than one processor', GM_FATAL) endif ! Forbidden options associated with Glissade dycore @@ -1407,6 +1416,11 @@ subroutine print_options(model) call write_log(message) endif + if (model%options%smooth_input_usrf) then + write(message,*) ' Input usrf will be smoothed' + call write_log(message) + endif + if (model%options%smooth_input_topography) then write(message,*) ' Input topography will be smoothed' call write_log(message) @@ -1766,6 +1780,16 @@ subroutine print_options(model) call write_log('Error, HO basal water input out of range', GM_FATAL) end if + if (model%options%which_ho_bwat == HO_BWAT_FLUX_ROUTING) then + write(message,*) 'ho_flux_routing_scheme : ',model%options%ho_flux_routing_scheme, & + ho_flux_routing_scheme(model%options%ho_flux_routing_scheme) + call write_log(message) + if (model%options%ho_flux_routing_scheme < 0.or. & + model%options%ho_flux_routing_scheme >= size(ho_flux_routing_scheme)) then + call write_log('Error, HO flux routing scheme out of range', GM_FATAL) + end if + end if + write(message,*) 'ho_whicheffecpress : ',model%options%which_ho_effecpress, & ho_whicheffecpress(model%options%which_ho_effecpress) call write_log(message) @@ -1996,7 +2020,7 @@ subroutine handle_parameters(section, model) use glimmer_config use glide_types use glimmer_log - use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt + use glimmer_physcon, only: rhoi, rhoo, grav, shci, lhci, trpt, n_glen implicit none type(ConfigSection), pointer :: section @@ -2021,6 +2045,7 @@ subroutine handle_parameters(section, model) call GetValue(section,'lhci', lhci) call GetValue(section,'trpt', trpt) #endif + call GetValue(section,'n_glen', n_glen) loglevel = GM_levels-GM_ERROR call GetValue(section,'log_level',loglevel) @@ -2033,9 +2058,9 @@ subroutine handle_parameters(section, model) call GetValue(section,'pmp_offset', model%temper%pmp_offset) call GetValue(section,'pmp_threshold', model%temper%pmp_threshold) call GetValue(section,'geothermal', model%paramets%geot) - !TODO - Change default_flwa to flwa_constant? Would have to change config files. call GetValue(section,'flow_factor', model%paramets%flow_enhancement_factor) call GetValue(section,'flow_factor_float', model%paramets%flow_enhancement_factor_float) + !TODO - Change default_flwa to flwa_constant? Would have to change config files. call GetValue(section,'default_flwa', model%paramets%default_flwa) call GetValue(section,'efvs_constant', model%paramets%efvs_constant) call GetValue(section,'effstrain_min', model%paramets%effstrain_min) @@ -2103,9 +2128,9 @@ subroutine handle_parameters(section, model) call GetValue(section, 'effecpress_bmlt_threshold', model%basal_physics%effecpress_bmlt_threshold) ! basal water parameters - call GetValue(section, 'const_bwat', model%basal_physics%const_bwat) - call GetValue(section, 'bwat_till_max', model%basal_physics%bwat_till_max) - call GetValue(section, 'c_drainage', model%basal_physics%c_drainage) + call GetValue(section, 'const_bwat', model%basal_hydro%const_bwat) + call GetValue(section, 'bwat_till_max', model%basal_hydro%bwat_till_max) + call GetValue(section, 'c_drainage', model%basal_hydro%c_drainage) ! pseudo-plastic parameters !TODO - Put pseudo-plastic and other basal sliding parameters in a separate section @@ -2206,7 +2231,7 @@ end subroutine handle_parameters subroutine print_parameters(model) - use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav + use glimmer_physcon, only: rhoi, rhoo, lhci, shci, trpt, grav, n_glen use glide_types use glimmer_log implicit none @@ -2371,6 +2396,9 @@ subroutine print_parameters(model) write(message,*) 'triple point of water (K) : ', trpt call write_log(message) + write(message,*) 'Glen flow law exponent : ', n_glen + call write_log(message) + write(message,*) 'geothermal flux (W/m^2) : ', model%paramets%geot call write_log(message) @@ -2627,12 +2655,12 @@ subroutine print_parameters(model) endif if (model%options%which_ho_bwat == HO_BWAT_CONSTANT) then - write(message,*) 'constant basal water depth (m): ', model%basal_physics%const_bwat + write(message,*) 'constant basal water depth (m): ', model%basal_hydro%const_bwat call write_log(message) elseif (model%options%which_ho_bwat == HO_BWAT_LOCAL_TILL) then - write(message,*) 'maximum till water depth (m) : ', model%basal_physics%bwat_till_max + write(message,*) 'maximum till water depth (m) : ', model%basal_hydro%bwat_till_max call write_log(message) - write(message,*) 'till drainage rate (m/yr) : ', model%basal_physics%c_drainage + write(message,*) 'till drainage rate (m/yr) : ', model%basal_hydro%c_drainage call write_log(message) endif @@ -3323,7 +3351,7 @@ subroutine define_glide_restart_variables(options) ! other Glissade options ! If overwriting acab in certain grid cells, than overwrite_acab_mask needs to be in the restart file. - ! This mask is set at model initialization based on the input acab or ice thickness. + ! This mask is read in at model initialization, or is set based on the input acab or ice thickness. if (options%overwrite_acab /= 0) then call glide_add_to_restart_variable_list('overwrite_acab_mask') endif diff --git a/libglide/glide_temp.F90 b/libglide/glide_temp.F90 index 23cb5293..7fb47060 100644 --- a/libglide/glide_temp.F90 +++ b/libglide/glide_temp.F90 @@ -512,7 +512,7 @@ subroutine glide_temp_driver(model,whichtemp) call corrpmpt(model%temper%temp(:,ew,ns), & model%geometry%thck(ew,ns), & - model%temper%bwat(ew,ns), & + model%basal_hydro%bwat(ew,ns), & model%numerics%sigma, & model%general%upn) @@ -560,7 +560,7 @@ subroutine glide_temp_driver(model,whichtemp) call corrpmpt(model%temper%temp(:,ew,ns), & model%geometry%thck(ew,ns), & - model%temper%bwat(ew,ns), & + model%basal_hydro%bwat(ew,ns), & model%numerics%sigma, & model%general%upn) diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index 147608e9..35402833 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -104,7 +104,6 @@ module glide_types integer, parameter :: FLWA_CONST_FLWA = 0 integer, parameter :: FLWA_PATERSON_BUDD_CONST_TEMP = 1 integer, parameter :: FLWA_PATERSON_BUDD = 2 - integer, parameter :: FLWA_INPUT = 3 integer, parameter :: BTRC_ZERO = 0 integer, parameter :: BTRC_CONSTANT = 1 @@ -161,6 +160,7 @@ module glide_types integer, parameter :: OVERWRITE_ACAB_NONE = 0 integer, parameter :: OVERWRITE_ACAB_ZERO_ACAB = 1 integer, parameter :: OVERWRITE_ACAB_THCKMIN = 2 + integer, parameter :: OVERWRITE_ACAB_INPUT_MASK = 3 integer, parameter :: GTHF_UNIFORM = 0 integer, parameter :: GTHF_PRESCRIBED_2D = 1 @@ -283,6 +283,11 @@ module glide_types integer, parameter :: HO_BWAT_NONE = 0 integer, parameter :: HO_BWAT_CONSTANT = 1 integer, parameter :: HO_BWAT_LOCAL_TILL = 2 + integer, parameter :: HO_BWAT_FLUX_ROUTING = 3 + + integer, parameter :: HO_FLUX_ROUTING_D8 = 0 + integer, parameter :: HO_FLUX_ROUTING_DINF = 1 + integer, parameter :: HO_FLUX_ROUTING_FD8 = 2 !TODO - Remove option 2? Rarely used integer, parameter :: HO_EFFECPRESS_OVERBURDEN = 0 @@ -290,6 +295,7 @@ module glide_types integer, parameter :: HO_EFFECPRESS_BMLT = 2 integer, parameter :: HO_EFFECPRESS_OCEAN_PENETRATION = 3 integer, parameter :: HO_EFFECPRESS_BWAT = 4 + integer, parameter :: HO_EFFECPRESS_BWAT_RAMP = 5 !WHL - added Picard acceleration option integer, parameter :: HO_NONLIN_PICARD = 0 @@ -398,9 +404,12 @@ module glide_types integer :: global_bc = 0 ! 0 for periodic, 1 for outflow, 2 for no_penetration, 3 for no_ice - !WHL - added to handle the active-blocks option + ! added to handle the active-blocks option integer, dimension(:,:), pointer :: ice_domain_mask => null() ! = 1 for cells that are potentially active + ! mask to identify cells at the edge of the global domain + integer, dimension(:,:), pointer :: global_edge_mask => null() ! = 1 for cells at edge of global domain + integer :: nx_block = 0 ! user-specified block sizes integer :: ny_block = 0 ! one task per block; optionally, tasks not assigned to inactive blocks @@ -470,7 +479,6 @@ module glide_types !> \item[1] \emph{Paterson and Budd} relationship, !> with temperature set to $-5^{\circ}\mathrm{C}$ !> \item[2] \emph{Paterson and Budd} relationship - !> \item[3] Read flwa/flwastag from file !> \end{description} integer :: whichbtrc = 0 @@ -589,6 +597,7 @@ module glide_types !> \item[0] Do not overwrite acab anywhere !> \item[1] Overwrite acab where input acab = 0 !> \item[2] Overwrite acab where input thickness <= threshold value + !> \item[3] Overwrite acab where input mask = 1 !> \end{description} integer :: gthf = 0 @@ -667,6 +676,9 @@ module glide_types logical :: adjust_input_thickness = .false. !> if true, then adjust thck to maintain usrf, instead of deriving usrf from topg and thck + logical :: smooth_input_usrf = .false. + !> if true, then apply Laplacian smoothing to usrf at initialization + logical :: smooth_input_topography = .false. !> if true, then apply Laplacian smoothing to the topography at initialization @@ -836,6 +848,15 @@ module glide_types !> \item[0] Set to zero everywhere !> \item[1] Set to constant everywhere, to force T = Tpmp. !> \item[2] Local basal till model with constant drainage + !> \item[3] Steady-state water routing with flux calculation + !> \end{description} + + integer :: ho_flux_routing_scheme = 0 + !> Flux routing scheme for basal water: + !> \begin{description} + !> \item[0] D8; send flux to lowest-elevation neighbor + !> \item[1] Dinf; divide flux between two lower-elevation neighbors + !> \item[2] FD8; divide flux amond all lower-elevation neighbors !> \end{description} integer :: which_ho_effecpress = 0 @@ -845,7 +866,8 @@ module glide_types !> \item[1] N is reduced where the bed is at or near the pressure melting point !> \item[2] N is reduced where there is melting at the bed !> \item[3] N is reduced due to connection of subglacial water to the ocean - !> \item[4] N is reduced where basal water is present + !> \item[4] N is reduced where basal water is present, following Bueler/van Pelt + !> \item[5] N is reduced where basal water is present, with a ramp function !> \end{description} integer :: which_ho_nonlinear = 0 @@ -1525,12 +1547,6 @@ module glide_types real(dp),dimension(:,:), pointer :: lcondflx => null() !> conductive heat flux (W/m^2) at lower sfc (positive down) real(dp),dimension(:,:), pointer :: dissipcol => null() !> total heat dissipation rate (W/m^2) in column (>= 0) - ! fields related to basal water - !TODO - Move these fields to the basal_physics type? - real(dp),dimension(:,:), pointer :: bwat => null() !> Basal water depth - real(dp),dimension(:,:), pointer :: bwatflx => null() !> Basal water flux - real(dp),dimension(:,:), pointer :: stagbwat => null() !> Basal water depth on velo grid - real(dp) :: pmp_offset = 5.0d0 ! offset of initial Tbed from pressure melting point temperature (deg C) real(dp) :: pmp_threshold = 1.0d-3 ! bed is assumed thawed where Tbed >= pmptemp - pmp_threshold (deg C) @@ -1819,9 +1835,46 @@ module glide_types !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + type glide_basal_hydro + + !> Holds variables related to basal hydrology + !> See glissade_bwater.F90 for usage details + + ! fields related to basal water + ! Note: Ideally, bwat should have MKS units (m), but currently is scaled. + real(dp),dimension(:,:), pointer :: bwat => null() !> Basal water depth + real(dp),dimension(:,:), pointer :: stagbwat => null() !> Basal water depth on velo grid + real(dp),dimension(:,:), pointer :: bwatflx => null() !> Basal water flux (m^3/s) + real(dp),dimension(:,:), pointer :: head => null() !> Hydraulic head (m) + + ! parameter for constant basal water + ! Note: This parameter applies to teh case HO_BWAT_CONSTANT. + ! For Glide's BWATER_CONST, the constant value is hardwired in subroutine calcbwat. + real(dp) :: const_bwat = 10.d0 !> constant basal water depth (m) + + ! parameters for local till model + ! These parameters apply to the case HO_BWAT_LOCAL_TILL. + ! The default values are from Aschwanden et al. (2016) and Bueler and van Pelt (2015). + real(dp) :: bwat_till_max = 2.0d0 !> maximum water depth in till (m) + real(dp) :: c_drainage = 1.0d-3 !> uniform drainage rate (m/yr) + real(dp) :: N_0 = 1000.d0 !> reference effective pressure (Pa) + real(dp) :: e_0 = 0.69d0 !> reference void ratio (dimensionless) + real(dp) :: C_c = 0.12d0 !> till compressibility (dimensionless) + !> Note: The ratio (e_0/C_c) is the key parameter + + ! parameters for steady-state flux-routing model + ! Could add visc_water and omega_hydro here; currently set in glissade_bwater module + ! Some of these parameters might apply to more general models like SHAKTI + + end type glide_basal_hydro + + !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + type glide_basal_physics + !> Holds variables related to basal physics associated with ice dynamics !> See glissade_basal_traction.F90 for usage details + !> TODO: Divide into separate types for basal friction/sliding and basal hydrology? !Note: By default, beta_grounded_min is set to a small nonzero value. ! Larger values (~10 to 100 Pa yr/m) might be needed for stability in realistic simulations. @@ -1891,22 +1944,8 @@ module glide_types real(dp) :: beta_powerlaw_umax = 0.0d0 !> upper limit of ice speed (m/yr) when evaluating powerlaw beta !> Where u > umax, let u = umax when evaluating beta(u) - ! parameter for constant basal water - ! Note: This parameter applies to HO_BWAT_CONSTANT only. - ! For Glide's BWATER_CONST, the constant value is hardwired in subroutine calcbwat. - real(dp) :: const_bwat = 10.d0 !> constant basal water depth (m) - - ! parameters for local till model - ! The default values are from Aschwanden et al. (2016) and Bueler and van Pelt (2015). - real(dp) :: bwat_till_max = 2.0d0 !> maximum water depth in till (m) - real(dp) :: c_drainage = 1.0d-3 !> uniform drainage rate (m/yr) - real(dp) :: N_0 = 1000.d0 !> reference effective pressure (Pa) - real(dp) :: e_0 = 0.69d0 !> reference void ratio (dimensionless) - real(dp) :: C_c = 0.12d0 !> till compressibility (dimensionless) - !> Note: The ratio (e_0/C_c) is the key parameter - ! Note: A basal process model is not currently supported, but a specified mintauf can be passed to subroutine calcbeta - ! to simulate a plastic bed.. + ! to simulate a plastic bed. real(dp),dimension(:,:) ,pointer :: mintauf => null() ! Bed strength (yield stress) calculated with basal process model end type glide_basal_physics @@ -2148,6 +2187,7 @@ module glide_types !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !TODO - Move these parameters to types associated with a certain kind of physics + !TODO - Set default geot = 0, so that idealized tests by default have no mass loss type glide_paramets real(dp),dimension(5) :: bpar = (/ 0.2d0, 0.5d0, 0.0d0 ,1.0d-2, 1.0d0/) real(dp) :: btrac_const = 0.d0 ! m yr^{-1} Pa^{-1} (gets scaled during init) @@ -2278,6 +2318,7 @@ module glide_types type(eismint_climate_type) :: eismint_climate type(glide_calving) :: calving type(glide_temper) :: temper + type(glide_basal_hydro) :: basal_hydro type(glide_basal_physics):: basal_physics type(glide_basal_melt) :: basal_melt type(glide_ocean_data) :: ocean_data @@ -2316,7 +2357,6 @@ subroutine glide_allocarr(model) !> \item \texttt{bheatflx(ewn,nsn)} !> \item \texttt{flwa(upn,ewn,nsn)} !WHL - 2 choices !> \item \texttt{dissip(upn,ewn,nsn)} !WHL - 2 choices - !> \item \texttt{bwat(ewn,nsn)} !> \item \texttt{bfricflx(ewn,nsn)} !> \item \texttt{ucondflx(ewn,nsn)} !> \item \texttt{lcondflx(ewn,nsn)} @@ -2496,6 +2536,9 @@ subroutine glide_allocarr(model) ! ice domain mask (to identify active blocks) call coordsystem_allocate(model%general%ice_grid, model%general%ice_domain_mask) + ! mask to identify cells at global domain edge + call coordsystem_allocate(model%general%ice_grid, model%general%global_edge_mask) + ! temperature arrays !NOTE: In the glide dycore (whichdycore = DYCORE_GLIDE), the temperature and @@ -2528,8 +2571,6 @@ subroutine glide_allocarr(model) model%temper%tempunstag(:,:,:) = unphys_val call coordsystem_allocate(model%general%ice_grid, model%temper%bheatflx) - call coordsystem_allocate(model%general%ice_grid, model%temper%bwat) - call coordsystem_allocate(model%general%velo_grid, model%temper%stagbwat) call coordsystem_allocate(model%general%ice_grid, model%temper%bpmp) call coordsystem_allocate(model%general%velo_grid, model%temper%stagbpmp) call coordsystem_allocate(model%general%ice_grid, model%temper%btemp) @@ -2538,9 +2579,14 @@ subroutine glide_allocarr(model) call coordsystem_allocate(model%general%velo_grid, model%temper%stagbtemp) call coordsystem_allocate(model%general%ice_grid, model%temper%ucondflx) - if (model%options%whichdycore == DYCORE_GLIDE) then ! glide only - call coordsystem_allocate(model%general%ice_grid, model%temper%bwatflx) - else ! glissade only + call coordsystem_allocate(model%general%ice_grid, model%basal_hydro%bwat) + call coordsystem_allocate(model%general%velo_grid, model%basal_hydro%stagbwat) + call coordsystem_allocate(model%general%ice_grid, model%basal_hydro%bwatflx) + if (model%options%which_ho_bwat == HO_BWAT_FLUX_ROUTING) then + call coordsystem_allocate(model%general%ice_grid, model%basal_hydro%head) + endif + + if (model%options%whichdycore == DYCORE_GLISSADE) then ! glissade only call coordsystem_allocate(model%general%ice_grid, model%temper%bfricflx) call coordsystem_allocate(model%general%ice_grid, model%temper%lcondflx) call coordsystem_allocate(model%general%ice_grid, model%temper%dissipcol) @@ -2939,6 +2985,8 @@ subroutine glide_deallocarr(model) deallocate(model%general%lon) if (associated(model%general%ice_domain_mask)) & deallocate(model%general%ice_domain_mask) + if (associated(model%general%global_edge_mask)) & + deallocate(model%general%global_edge_mask) ! vertical sigma coordinates @@ -2957,12 +3005,6 @@ subroutine glide_deallocarr(model) deallocate(model%temper%tempunstag) if (associated(model%temper%bheatflx)) & deallocate(model%temper%bheatflx) - if (associated(model%temper%bwat)) & - deallocate(model%temper%bwat) - if (associated(model%temper%bwatflx)) & - deallocate(model%temper%bwatflx) - if (associated(model%temper%stagbwat)) & - deallocate(model%temper%stagbwat) if (associated(model%temper%bpmp)) & deallocate(model%temper%bpmp) if (associated(model%temper%stagbpmp)) & @@ -3119,6 +3161,16 @@ subroutine glide_deallocarr(model) if (associated(model%stress%taudy)) & deallocate(model%stress%taudy) + ! basal hydrology arrays + if (associated(model%basal_hydro%bwat)) & + deallocate(model%basal_hydro%bwat) + if (associated(model%basal_hydro%stagbwat)) & + deallocate(model%basal_hydro%stagbwat) + if (associated(model%basal_hydro%bwatflx)) & + deallocate(model%basal_hydro%bwatflx) + if (associated(model%basal_hydro%head)) & + deallocate(model%basal_hydro%head) + ! basal physics arrays if (associated(model%basal_physics%bpmp_mask)) & deallocate(model%basal_physics%bpmp_mask) diff --git a/libglide/glide_vars.def b/libglide/glide_vars.def index 086fa520..a5847bde 100644 --- a/libglide/glide_vars.def +++ b/libglide/glide_vars.def @@ -1122,7 +1122,7 @@ load: 1 dimensions: time, y1, x1 units: meter long_name: basal water depth -data: data%temper%bwat +data: data%basal_hydro%bwat factor: thk0 load: 1 @@ -1130,8 +1130,15 @@ load: 1 dimensions: time, y1, x1 units: meter3/year long_name: basal water flux -data: data%temper%bwatflx -factor: thk0 +data: data%basal_hydro%bwatflx +factor: scyr + +[head] +dimensions: time, y1, x1 +units: meter +long_name: hydraulic head +data: data%basal_hydro%head +factor: 1 [effecpress] dimensions: time, y1, x1 diff --git a/libglide/glide_velo.F90 b/libglide/glide_velo.F90 index 10cdd120..675a0f56 100644 --- a/libglide/glide_velo.F90 +++ b/libglide/glide_velo.F90 @@ -1033,7 +1033,7 @@ subroutine calc_btrc(model,flag,btrc) do ns = 1,nsn-1 do ew = 1,ewn-1 - if (0.0d0 < model%temper%stagbwat(ew,ns)) then + if (0.0d0 < model%basal_hydro%stagbwat(ew,ns)) then btrc(ew,ns) = model%velocity%bed_softness(ew,ns) else btrc(ew,ns) = 0.0d0 @@ -1078,10 +1078,10 @@ subroutine calc_btrc(model,flag,btrc) do ns = 1,nsn-1 do ew = 1,ewn-1 - if (0.0d0 < model%temper%stagbwat(ew,ns)) then + if (0.0d0 < model%basal_hydro%stagbwat(ew,ns)) then btrc(ew,ns) = model%velowk%c(1) + model%velowk%c(2) * tanh(model%velowk%c(3) * & - model%temper%stagbwat(ew,ns) - model%velowk%c(4)) + model%basal_hydro%stagbwat(ew,ns) - model%velowk%c(4)) if (0.0d0 > sum(model%isostasy%relx(ew:ew+1,ns:ns+1))) then btrc(ew,ns) = btrc(ew,ns) * model%velowk%marine diff --git a/libglimmer/glimmer_paramets.F90 b/libglimmer/glimmer_paramets.F90 index 019f44e6..aa8b595d 100644 --- a/libglimmer/glimmer_paramets.F90 +++ b/libglimmer/glimmer_paramets.F90 @@ -33,7 +33,7 @@ module glimmer_paramets use glimmer_global, only : dp - use glimmer_physcon, only : scyr, gn + use glimmer_physcon, only : scyr implicit none save @@ -118,6 +118,7 @@ module glimmer_paramets real(dp), parameter :: grav_glam = 9.81d0 ! m s^{-2} ! GLAM scaling parameters; units are correct if thk0 has units of meters + integer, parameter :: gn = 3 ! Glen flow exponent; fixed at 3 for purposes of setting vis0 real(dp), parameter :: tau0 = rhoi_glam*grav_glam*thk0 ! stress scale in GLAM ( Pa ) real(dp), parameter :: evs0 = tau0 / (vel0/len0) ! eff. visc. scale in GLAM ( Pa s ) real(dp), parameter :: vis0 = tau0**(-gn) * (vel0/len0) ! rate factor scale in GLAM ( Pa^-3 s^-1 ) diff --git a/libglimmer/glimmer_physcon.F90 b/libglimmer/glimmer_physcon.F90 index 0c990d29..f697bf3e 100644 --- a/libglimmer/glimmer_physcon.F90 +++ b/libglimmer/glimmer_physcon.F90 @@ -79,11 +79,17 @@ module glimmer_physcon ! real(dp) :: trpt = 273.16d0 !< Triple point of water (K) #endif + ! The Glen flow-law exponent, n_glen, is used in Glissade. + ! It is not a parameter, because the default can be overridden in the config file. + ! TODO: Allow setting n_glen independently for each ice sheet instance? + ! Note: Earlier code used an integer parameter, gn = 3, for all flow-law calculations. + ! For backward compatiblity, gn = 3 is retained for Glide. + real(dp) :: n_glen = 3.0d0 !< Exponent in Glen's flow law; user-configurable real(dp) in Glissade + integer, parameter :: gn = 3 !< Exponent in Glen's flow law; fixed integer parameter in Glide real(dp),parameter :: celsius_to_kelvin = 273.15d0 !< Note: Not quite equal to trpt real(dp),parameter :: scyr = 31536000.d0 !< Number of seconds in a year of exactly 365 days real(dp),parameter :: rhom = 3300.0d0 !< The density of magma(?) (kg m-3) real(dp),parameter :: rhos = 2600.0d0 !< The density of solid till (kg m$^{-3}$) - integer, parameter :: gn = 3 !< The power dependency of Glen's flow law. real(dp),parameter :: actenh = 139.0d3 !< Activation energy in Glen's flow law for \f$T^{*}\geq263\f$K. (J mol-1) real(dp),parameter :: actenl = 60.0d3 !< Activation energy in Glen's flow law for \f$T^{*}<263\f$K. (J mol-1) real(dp),parameter :: arrmlh = 1.733d3 !< Constant of proportionality in Arrhenius relation diff --git a/libglimmer/glimmer_scales.F90 b/libglimmer/glimmer_scales.F90 index f95c6a86..380ff2b8 100644 --- a/libglimmer/glimmer_scales.F90 +++ b/libglimmer/glimmer_scales.F90 @@ -45,7 +45,7 @@ subroutine glimmer_init_scales ! set scale factors for I/O (can't have non-integer powers) - use glimmer_physcon, only : scyr, gn + use glimmer_physcon, only : scyr use glimmer_paramets, only : thk0, tim0, vel0, vis0, len0, acc0, tau0, evs0 implicit none diff --git a/libglimmer/parallel_mpi.F90 b/libglimmer/parallel_mpi.F90 index fa8ed875..0b7d6d29 100644 --- a/libglimmer/parallel_mpi.F90 +++ b/libglimmer/parallel_mpi.F90 @@ -287,12 +287,19 @@ module cism_parallel module procedure parallel_get_var_real8_2d end interface + interface parallel_global_sum + module procedure parallel_global_sum_integer_2d + module procedure parallel_global_sum_real4_2d + module procedure parallel_global_sum_real8_2d + end interface + interface parallel_halo module procedure parallel_halo_integer_2d module procedure parallel_halo_logical_2d module procedure parallel_halo_real4_2d module procedure parallel_halo_real8_2d module procedure parallel_halo_real8_3d + module procedure parallel_halo_real8_4d end interface interface parallel_halo_extrapolate @@ -2576,7 +2583,7 @@ subroutine distributed_grid(ewn, nsn, & ! of the global domain, so staggered_ilo = staggered_jlo = staggered_lhalo on ! processors that include these rows. ! Note: For no_ice BC, we assume (uvel,vvel) = 0 along the global boundary. - ! In this case, vertices along the southern and western rows of the global boundary + ! In this case, vertices along the southern and western edges of the global boundary ! are not considered to be locally owned by any task. if (outflow_bc .and. this_rank <= west) then ! on west edge of global domain @@ -5668,6 +5675,59 @@ function parallel_get_var_real8_2d(ncid, varid, values) end function parallel_get_var_real8_2d +!======================================================================= + + subroutine parallel_global_edge_mask(global_edge_mask, parallel) + + ! Create a mask = 1 in locally owned cells at the edge of the global domain, + ! = 0 elsewhere + + integer, dimension(:,:), intent(out) :: global_edge_mask + type(parallel_type) :: parallel + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn, & + east => parallel%east, & + west => parallel%west, & + north => parallel%north, & + south => parallel%south) + + ! Check array dimensions + + ! unknown grid + if (size(global_edge_mask,1)/=local_ewn .or. size(global_edge_mask,2)/=local_nsn) then + write(*,*) "Unknown Grid: Size a=(", size(global_edge_mask,1), ",", size(global_edge_mask,2), & + ") and local_ewn and local_nsn = ", local_ewn, ",", local_nsn + call parallel_stop(__FILE__,__LINE__) + endif + + ! Identify cells at the edge of the global domain + + global_edge_mask = 0 + + if (this_rank >= east) then ! at east edge of global domain + global_edge_mask(local_ewn-uhalo,:) = 1 + endif + + if (this_rank <= west) then ! at west edge of global domain + global_edge_mask(lhalo+1,:) = 1 + endif + + if (this_rank >= north) then ! at north edge of global domain + global_edge_mask(:,local_nsn-uhalo) = 1 + endif + + if (this_rank <= south) then ! at south edge of global domain + global_edge_mask(:,lhalo+1) = 1 + endif + + call parallel_halo(global_edge_mask, parallel) + + end associate + + end subroutine parallel_global_edge_mask + !======================================================================= !TODO - Is function parallel_globalID still needed? No longer called except from glissade_test_halo. @@ -5815,6 +5875,93 @@ subroutine parallel_globalindex(ilocal, jlocal, iglobal, jglobal, parallel) end subroutine parallel_globalindex +!======================================================================= + + function parallel_global_sum_integer_2d(a, parallel) + + ! Calculates the global sum of a 2D integer field + + integer,dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + integer :: local_sum + integer :: parallel_global_sum_integer_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0 + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_integer_2d = parallel_reduce_sum(local_sum) + + end associate + + end function parallel_global_sum_integer_2d + +!======================================================================= + + function parallel_global_sum_real4_2d(a, parallel) + + ! Calculates the global sum of a 2D single-precision field + + real(sp),dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + real(sp) :: local_sum + real(sp) :: parallel_global_sum_real4_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0.0 + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_real4_2d = parallel_reduce_sum(local_sum) + + end associate + + end function parallel_global_sum_real4_2d + +!======================================================================= + + function parallel_global_sum_real8_2d(a, parallel) + + ! Calculates the global sum of a 2D double-precision field + + real(dp),dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + real(dp) :: local_sum + real(dp) :: parallel_global_sum_real8_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0.0d0 + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_real8_2d = parallel_reduce_sum(local_sum) + + end associate + + end function parallel_global_sum_real8_2d + !======================================================================= subroutine parallel_localindex(iglobal, jglobal, ilocal, jlocal, rlocal, parallel) @@ -6519,6 +6666,130 @@ subroutine parallel_halo_real8_3d(a, parallel) end subroutine parallel_halo_real8_3d + + subroutine parallel_halo_real8_4d(a, parallel) + + use mpi_mod + implicit none + real(dp),dimension(:,:,:,:) :: a + type(parallel_type) :: parallel + + integer :: erequest,ierror,one,nrequest,srequest,wrequest + real(dp),dimension(size(a,1), size(a,2), lhalo, parallel%local_nsn-lhalo-uhalo) :: esend,wrecv + real(dp),dimension(size(a,1), size(a,2), uhalo, parallel%local_nsn-lhalo-uhalo) :: erecv,wsend + real(dp),dimension(size(a,1), size(a,2), parallel%local_ewn, lhalo) :: nsend,srecv + real(dp),dimension(size(a,1), size(a,2), parallel%local_ewn, uhalo) :: nrecv,ssend + + ! begin + associate( & + outflow_bc => parallel%outflow_bc, & + no_ice_bc => parallel%no_ice_bc, & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn, & + east => parallel%east, & + west => parallel%west, & + north => parallel%north, & + south => parallel%south, & + southwest_corner => parallel%southwest_corner, & + southeast_corner => parallel%southeast_corner, & + northeast_corner => parallel%northeast_corner, & + northwest_corner => parallel%northwest_corner & + ) + + ! staggered grid + if (size(a,3)==local_ewn-1.and.size(a,4)==local_nsn-1) return + + ! unknown grid + if (size(a,3)/=local_ewn.or.size(a,4)/=local_nsn) then + write(*,*) "Unknown Grid: Size a=(", size(a,1), ",", size(a,2), ",", size(a,3), ",", size(a,4), ") & + &and local_ewn and local_nsn = ", local_ewn, ",", local_nsn + call parallel_stop(__FILE__,__LINE__) + endif + + ! unstaggered grid + call mpi_irecv(wrecv,size(wrecv),mpi_real8,west,west,& + comm,wrequest,ierror) + call mpi_irecv(erecv,size(erecv),mpi_real8,east,east,& + comm,erequest,ierror) + call mpi_irecv(srecv,size(srecv),mpi_real8,south,south,& + comm,srequest,ierror) + call mpi_irecv(nrecv,size(nrecv),mpi_real8,north,north,& + comm,nrequest,ierror) + + esend(:,:,:,:) = & + a(:,:,local_ewn-uhalo-lhalo+1:local_ewn-uhalo,1+lhalo:local_nsn-uhalo) + call mpi_send(esend,size(esend),mpi_real8,east,this_rank,comm,ierror) + wsend(:,:,:,:) = a(:,:,1+lhalo:1+lhalo+uhalo-1,1+lhalo:local_nsn-uhalo) + call mpi_send(wsend,size(wsend),mpi_real8,west,this_rank,comm,ierror) + + call mpi_wait(wrequest,mpi_status_ignore,ierror) + a(:,:,:lhalo,1+lhalo:local_nsn-uhalo) = wrecv(:,:,:,:) + call mpi_wait(erequest,mpi_status_ignore,ierror) + a(:,:,local_ewn-uhalo+1:,1+lhalo:local_nsn-uhalo) = erecv(:,:,:,:) + + nsend(:,:,:,:) = a(:,:,:,local_nsn-uhalo-lhalo+1:local_nsn-uhalo) + call mpi_send(nsend,size(nsend),mpi_real8,north,this_rank,comm,ierror) + ssend(:,:,:,:) = a(:,:,:,1+lhalo:1+lhalo+uhalo-1) + call mpi_send(ssend,size(ssend),mpi_real8,south,this_rank,comm,ierror) + + call mpi_wait(srequest,mpi_status_ignore,ierror) + a(:,:,:,:lhalo) = srecv(:,:,:,:) + call mpi_wait(nrequest,mpi_status_ignore,ierror) + a(:,:,:,local_nsn-uhalo+1:) = nrecv(:,:,:,:) + + if (outflow_bc) then ! set values in global halo to zero + ! interior halo cells should not be affected + + if (this_rank >= east) then ! at east edge of global domain + a(:,:,local_ewn-uhalo+1:,:) = 0.d0 + endif + + if (this_rank <= west) then ! at west edge of global domain + a(:,:,:lhalo,:) = 0.d0 + endif + + if (this_rank >= north) then ! at north edge of global domain + a(:,:,:,local_nsn-uhalo+1:) = 0.d0 + endif + + if (this_rank <= south) then ! at south edge of global domain + a(:,:,:,:lhalo) = 0.d0 + endif + + elseif (no_ice_bc) then + + ! Set values to zero in cells adjacent to the global boundary; + ! includes halo cells and one row of locally owned cells + + if (this_rank >= east) then ! at east edge of global domain + a(:,:,local_ewn-uhalo:,:) = 0.d0 + endif + + if (this_rank <= west) then ! at west edge of global domain + a(:,:,:lhalo+1,:) = 0.d0 + endif + + if (this_rank >= north) then ! at north edge of global domain + a(:,:,:,local_nsn-uhalo:) = 0.d0 + endif + + if (this_rank <= south) then ! at south edge of global domain + a(:,:,:,:lhalo+1) = 0.d0 + endif + + ! Some interior blocks have a single cell at a corner of the global boundary. + ! Set values in corner cells to zero, along with adjacent halo cells. + if (southwest_corner) a(:,:,:lhalo+1,:lhalo+1) = 0.d0 + if (southeast_corner) a(:,:,local_ewn-lhalo:,:lhalo+1) = 0.d0 + if (northeast_corner) a(:,:,local_ewn-lhalo:,local_nsn-lhalo:) = 0.d0 + if (northwest_corner) a(:,:,:lhalo+1,local_nsn-lhalo:) = 0.d0 + + endif ! outflow or no_ice bc + + end associate + + end subroutine parallel_halo_real8_4d + !======================================================================= ! subroutines belonging to the parallel_halo_extrapolate interface diff --git a/libglimmer/parallel_slap.F90 b/libglimmer/parallel_slap.F90 index df5991cf..f0ac86b9 100644 --- a/libglimmer/parallel_slap.F90 +++ b/libglimmer/parallel_slap.F90 @@ -254,12 +254,19 @@ module cism_parallel module procedure parallel_get_var_real8_2d end interface + interface parallel_global_sum + module procedure parallel_global_sum_integer_2d + module procedure parallel_global_sum_real4_2d + module procedure parallel_global_sum_real8_2d + end interface + interface parallel_halo module procedure parallel_halo_integer_2d module procedure parallel_halo_logical_2d module procedure parallel_halo_real4_2d module procedure parallel_halo_real8_2d module procedure parallel_halo_real8_3d + module procedure parallel_halo_real8_4d end interface interface parallel_halo_extrapolate @@ -2394,6 +2401,42 @@ function parallel_get_var_real8_2d(ncid, varid, values) end function parallel_get_var_real8_2d +!======================================================================= + + subroutine parallel_global_edge_mask(global_edge_mask, parallel) + + ! Create a mask = 1 in locally owned cells at the edge of the global domain, + ! = 0 elsewhere + + integer, dimension(:,:), intent(out) :: global_edge_mask + type(parallel_type) :: parallel + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + ! Check array dimensions + + ! unknown grid + if (size(global_edge_mask,1)/=local_ewn .or. size(global_edge_mask,2)/=local_nsn) then + write(*,*) "Unknown Grid: Size a=(", size(global_edge_mask,1), ",", size(global_edge_mask,2), & + ") and local_ewn and local_nsn = ", local_ewn, ",", local_nsn + call parallel_stop(__FILE__,__LINE__) + endif + + ! Identify cells at the edge of the global domain + + global_edge_mask = 0 + + global_edge_mask(lhalo+1,:) = 1 + global_edge_mask(local_ewn-uhalo,:) = 1 + global_edge_mask(:,lhalo+1) = 1 + global_edge_mask(:,local_nsn-uhalo) = 1 + + end associate + + end subroutine parallel_global_edge_mask + !======================================================================= !TODO - Is function parallel_globalID still needed? No longer called except from glissade_test_halo. @@ -2468,6 +2511,91 @@ function parallel_globalID_scalar(locew, locns, upstride, parallel) end function parallel_globalID_scalar +!======================================================================= + + function parallel_global_sum_integer_2d(a, parallel) + + ! Calculates the global sum of a 2D integer field + + integer,dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + integer :: local_sum + integer :: parallel_global_sum_integer_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0 + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_integer_2d = local_sum + + end associate + + end function parallel_global_sum_integer_2d + + + function parallel_global_sum_real4_2d(a, parallel) + + ! Calculates the global sum of a 2D single-precision field + + real(sp),dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + real(sp) :: local_sum + real(sp) :: parallel_global_sum_real4_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0. + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_real4_2d = local_sum + + end associate + + end function parallel_global_sum_real4_2d + + + function parallel_global_sum_real8_2d(a, parallel) + + ! Calculates the global sum of a 2D integer field + + real(dp),dimension(:,:),intent(in) :: a + type(parallel_type) :: parallel + + integer :: i, j + real(dp) :: local_sum + real(dp) :: parallel_global_sum_real8_2d + + associate( & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + local_sum = 0.0d0 + do j = nhalo+1, local_nsn-nhalo + do i = nhalo+1, local_ewn-nhalo + local_sum = local_sum + a(i,j) + enddo + enddo + parallel_global_sum_real8_2d = local_sum + + end associate + + end function parallel_global_sum_real8_2d + !======================================================================= subroutine parallel_globalindex(ilocal, jlocal, iglobal, jglobal, parallel) @@ -2852,6 +2980,68 @@ subroutine parallel_halo_real8_3d(a, parallel) end subroutine parallel_halo_real8_3d + + subroutine parallel_halo_real8_4d(a, parallel) + + implicit none + real(dp),dimension(:,:,:,:) :: a + type(parallel_type) :: parallel + + real(dp),dimension(size(a,1),size(a,2),lhalo,parallel%local_nsn-lhalo-uhalo) :: ecopy + real(dp),dimension(size(a,1),size(a,2),uhalo,parallel%local_nsn-lhalo-uhalo) :: wcopy + real(dp),dimension(size(a,1),size(a,2),parallel%local_ewn,lhalo) :: ncopy + real(dp),dimension(size(a,1),size(a,2),parallel%local_ewn,uhalo) :: scopy + + ! begin + + associate( & + outflow_bc => parallel%outflow_bc, & + no_ice_bc => parallel%no_ice_bc, & + local_ewn => parallel%local_ewn, & + local_nsn => parallel%local_nsn) + + ! staggered grid + if (size(a,3)==local_ewn-1 .and. size(a,4)==local_nsn-1) return + + ! unknown grid + if (size(a,3)/=local_ewn .or. size(a,4)/=local_nsn) then + write(*,*) "Unknown Grid: Size a=(", size(a,2), ",", size(a,3), ",", size(a,4), ") & + &and local_ewn and local_nsn = ", local_ewn, ",", local_nsn + call parallel_stop(__FILE__,__LINE__) + endif + + if (outflow_bc) then + + a(:,:,:lhalo,1+lhalo:local_nsn-uhalo) = 0.d0 + a(:,:,local_ewn-uhalo+1:,1+lhalo:local_nsn-uhalo) = 0.d0 + a(:,:,:,:lhalo) = 0.d0 + a(:,:,:,local_nsn-uhalo+1:) = 0.d0 + + elseif (no_ice_bc) then + + a(:,:,:lhalo+1,1+lhalo:local_nsn-uhalo) = 0.d0 + a(:,:,local_ewn-uhalo:,1+lhalo:local_nsn-uhalo) = 0.d0 + a(:,:,:,:lhalo+1) = 0.d0 + a(:,:,:,local_nsn-uhalo:) = 0.d0 + + else ! periodic BC + + ecopy(:,:,:,:) = a(:,:,local_ewn-uhalo-lhalo+1:local_ewn-uhalo,1+lhalo:local_nsn-uhalo) + wcopy(:,:,:,:) = a(:,:,1+lhalo:1+lhalo+uhalo-1,1+lhalo:local_nsn-uhalo) + a(:,:,:lhalo,1+lhalo:local_nsn-uhalo) = ecopy(:,:,:,:) + a(:,:,local_ewn-uhalo+1:,1+lhalo:local_nsn-uhalo) = wcopy(:,:,:,:) + + ncopy(:,:,:,:) = a(:,:,:,local_nsn-uhalo-lhalo+1:local_nsn-uhalo) + scopy(:,:,:,:) = a(:,:,:,1+lhalo:1+lhalo+uhalo-1) + a(:,:,:,:lhalo) = ncopy(:,:,:,:) + a(:,:,:,local_nsn-uhalo+1:) = scopy(:,:,:,:) + + endif + + end associate + + end subroutine parallel_halo_real8_4d + !======================================================================= ! subroutines belonging to the parallel_halo_extrapolate interface diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index ef0dc5dd..614effae 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -90,7 +90,7 @@ subroutine glissade_initialise(model, evolve_ice) use cism_parallel, only: parallel_type, distributed_gather_var, & distributed_scatter_var, parallel_finalise, & - distributed_grid, distributed_grid_active_blocks, & + distributed_grid, distributed_grid_active_blocks, parallel_global_edge_mask, & parallel_halo, parallel_halo_extrapolate, parallel_reduce_max, & staggered_parallel_halo_extrapolate, staggered_no_penetration_mask, & parallel_create_comm_row, parallel_create_comm_col, not_parallel @@ -104,7 +104,7 @@ subroutine glissade_initialise(model, evolve_ice) use glissade_basal_water, only: glissade_basal_water_init use glissade_masks, only: glissade_get_masks, glissade_marine_connection_mask use glimmer_scales - use glimmer_paramets, only: eps11, thk0, len0, tim0 + use glimmer_paramets, only: eps11, thk0, len0, tim0, scyr use glimmer_physcon, only: rhow, rhoi use glide_mask use isostasy, only: init_isostasy, isos_relaxed @@ -117,8 +117,8 @@ subroutine glissade_initialise(model, evolve_ice) use glissade_inversion, only: glissade_init_inversion, verbose_inversion use glissade_bmlt_float, only: glissade_bmlt_float_thermal_forcing_init, verbose_bmlt_float use glissade_grounding_line, only: glissade_grounded_fraction - use glissade_utils, only: & - glissade_adjust_thickness, glissade_smooth_topography, glissade_adjust_topography + use glissade_utils, only: glissade_adjust_thickness, glissade_smooth_usrf, & + glissade_smooth_topography, glissade_adjust_topography use glissade_utils, only: glissade_stdev use felix_dycore_interface, only: felix_velo_init @@ -204,7 +204,7 @@ subroutine glissade_initialise(model, evolve_ice) model%general%ice_domain_mask = 0 ! Read ice_domain_mask from the input or restart file - ! Note: In generaly, input arrays are read from subroutine glide_io_readall (called below) in glide_io.F90. + ! Note: In general, input arrays are read from subroutine glide_io_readall (called below) in glide_io.F90. ! However, ice_domain_mask is needed now to identify active blocks. infile => model%funits%in_first ! assume ice_domain_mask is in the input or restart file @@ -218,11 +218,6 @@ subroutine glissade_initialise(model, evolve_ice) ! The subroutine will report how many tasks are needed to compute on all active blocks, and then abort. ! The user can then resubmit (on an optimal number of processors) with model%options%compute_blocks = ACTIVE_BLOCKS. -! call distributed_grid_active_blocks(model%general%ewn, model%general%nsn, & -! model%general%nx_block, model%general%ny_block, & -! model%general%ice_domain_mask, & -! inquire_only = .true.) - call distributed_grid_active_blocks(model%general%ewn, model%general%nsn, & model%general%nx_block, model%general%ny_block, & model%general%ice_domain_mask, & @@ -240,10 +235,6 @@ subroutine glissade_initialise(model, evolve_ice) model%general%global_bc = GLOBAL_BC_NO_ICE endif -! call distributed_grid_active_blocks(model%general%ewn, model%general%nsn, & -! model%general%nx_block, model%general%ny_block, & -! model%general%ice_domain_mask) - call distributed_grid_active_blocks(model%general%ewn, model%general%nsn, & model%general%nx_block, model%general%ny_block, & model%general%ice_domain_mask, & @@ -256,17 +247,11 @@ subroutine glissade_initialise(model, evolve_ice) elseif (model%general%global_bc == GLOBAL_BC_OUTFLOW) then -! call distributed_grid(model%general%ewn, model%general%nsn, global_bc_in = 'outflow') - - !WHL - temporary call to fill the parallel derived type call distributed_grid(model%general%ewn, model%general%nsn, & model%parallel, global_bc_in = 'outflow') - elseif (model%general%global_bc == GLOBAL_BC_NO_ICE) then -! call distributed_grid(model%general%ewn, model%general%nsn, global_bc_in = 'no_ice') - call distributed_grid(model%general%ewn, model%general%nsn, & model%parallel, global_bc_in = 'no_ice') @@ -276,8 +261,6 @@ subroutine glissade_initialise(model, evolve_ice) ! The difference is that we also use no-penetration masks for (uvel,vvel) at the global boundary ! (computed by calling staggered_no_penetration_mask below). -! call distributed_grid(model%general%ewn, model%general%nsn, global_bc_in = 'no_penetration') - call distributed_grid(model%general%ewn, model%general%nsn, & model%parallel, global_bc_in = 'no_penetration') @@ -319,8 +302,13 @@ subroutine glissade_initialise(model, evolve_ice) ! allocate arrays call glide_allocarr(model) - ! set masks at global boundary for no-penetration boundary conditions - ! this subroutine includes a halo update + ! Compute a mask to identify cells at the edge of the global domain + ! (Currently used only to compute bwat_mask for basal water routing) + ! Includes a halo update for global_edge_mask + call parallel_global_edge_mask(model%general%global_edge_mask, parallel) + + ! Set masks at global boundary for no-penetration boundary conditions + ! Includes a halo update for the masks if (model%general%global_bc == GLOBAL_BC_NO_PENETRATION) then call staggered_no_penetration_mask(model%velocity%umask_no_penetration, & model%velocity%vmask_no_penetration, & @@ -408,7 +396,15 @@ subroutine glissade_initialise(model, evolve_ice) call glissade_adjust_thickness(model) endif - ! Optionally, smooth the input topography with a 9-point Laplacian smoother. + ! Optionally, smooth the input surface elevation with a Laplacian smoother. + ! This subroutine does not change the topg, but returns thck consistent with the new usrf. + ! If the initial usrf is rough, then multiple smoothing passes may be needed to stabilize the flow. + + if (model%options%smooth_input_usrf .and. model%options%is_restart == RESTART_FALSE) then + call glissade_smooth_usrf(model, nsmooth = 5) + endif ! smooth_input_usrf + + ! Optionally, smooth the input topography with a Laplacian smoother. if (model%options%smooth_input_topography .and. model%options%is_restart == RESTART_FALSE) then call glissade_smooth_topography(model) @@ -457,7 +453,11 @@ subroutine glissade_initialise(model, evolve_ice) ! handle relaxed/equilibrium topo ! Initialise isostasy first - call init_isostasy(model) + if (model%options%isostasy == ISOSTASY_COMPUTE) then + + call init_isostasy(model) + + endif select case(model%isostasy%whichrelaxed) @@ -633,17 +633,17 @@ subroutine glissade_initialise(model, evolve_ice) if (make_ice_domain_mask) then where (model%geometry%thck > 0.0d0 .or. model%geometry%topg > 0.0d0) +!! where (model%geometry%thck > 0.0d0) ! uncomment for terrestrial margins model%general%ice_domain_mask = 1 elsewhere model%general%ice_domain_mask = 0 endwhere - ! Extend the mask a couple of cells in each direction to be on the safe side. + ! Extend the mask a few cells in each direction to be on the safe side. ! The number of buffer layers could be made a config parameter. allocate(ice_domain_mask(model%general%ewn,model%general%nsn)) -!! do k = 1, 2 do k = 1, 3 call parallel_halo(model%general%ice_domain_mask, parallel) ice_domain_mask = model%general%ice_domain_mask ! temporary copy @@ -737,8 +737,6 @@ subroutine glissade_initialise(model, evolve_ice) if (model%climate%overwrite_acab_value /= 0 .and. model%options%is_restart == RESTART_FALSE) then -!! print*, 'Setting acab = overwrite value (m/yr):', model%climate%overwrite_acab_value * scyr*thk0/tim0 - call glissade_overwrite_acab_mask(model%options%overwrite_acab, & model%climate%acab, & model%geometry%thck, & @@ -1172,6 +1170,14 @@ subroutine glissade_tstep(model, time) enddo write(6,*) ' ' enddo + print*, ' ' + print*, 'bmlt_ground (m/yr):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') model%basal_melt%bmlt_ground(i,j)*scyr + enddo + write(6,*) ' ' + enddo endif ! ------------------------------------------------------------------------ @@ -1305,8 +1311,8 @@ subroutine glissade_bmlt_float_solve(model) integer, dimension(model%general%ewn, model%general%nsn) :: & ice_mask, & ! = 1 if ice is present (thck > 0, else = 0 - floating_mask, & ! = 1 if ice is present (thck > 0) and floating - ocean_mask, & ! = 0 if ice is absent (thck = 0) and topg < 0 + floating_mask, & ! = 1 if ice is present (thck > 0) and floating, else = 0 + ocean_mask, & ! = 1 if topg is below sea level and ice is absent, else = 0 land_mask ! = 1 if topg - eus >= 0 real(dp), dimension(model%general%ewn, model%general%nsn) :: & @@ -1761,11 +1767,15 @@ subroutine glissade_thermal_solve(model, dt) use cism_parallel, only: parallel_type, parallel_halo use glimmer_paramets, only: tim0, thk0, len0 - use glimmer_physcon, only: scyr + use glimmer_physcon, only: rhow, rhoi, scyr use glissade_therm, only: glissade_therm_driver - use glissade_basal_water, only: glissade_calcbwat + use glissade_basal_water, only: glissade_calcbwat, glissade_bwat_flux_routing use glissade_transport, only: glissade_add_2d_anomaly use glissade_grid_operators, only: glissade_vertical_interpolate + use glissade_masks, only: glissade_get_masks + + !WHL - debug + use cism_parallel, only: parallel_reduce_max implicit none @@ -1784,6 +1794,15 @@ subroutine glissade_thermal_solve(model, dt) integer :: i, j, up integer :: itest, jtest, rtest + integer, dimension(model%general%ewn, model%general%nsn) :: & + ice_mask, & ! = 1 if ice is present (thck > thklim_temp), else = 0 + floating_mask, & ! = 1 if ice is present (thck > thklim_temp) and floating, else = 0 + ocean_mask, & ! = 1 if topg is below sea level and ice is absent, else = 0 + bwat_mask ! = 1 for cells through which basal water is routed, else = 0 + + !WHL - debug + real(dp) :: head_max + type(parallel_type) :: parallel ! info for parallel communication rtest = -999 @@ -1888,33 +1907,122 @@ subroutine glissade_thermal_solve(model, dt) model%temper%bheatflx, model%temper%bfricflx, & ! W/m2 model%temper%dissip, & ! deg/s model%temper%pmp_threshold, & ! deg C - model%temper%bwat*thk0, & ! m + model%basal_hydro%bwat*thk0, & ! m model%temper%temp, & ! deg C model%temper%waterfrac, & ! unitless model%temper%bpmp, & ! deg C model%temper%btemp_ground, & ! deg C model%temper%btemp_float, & ! deg C bmlt_ground_unscaled) ! m/s - + ! Update basal hydrology, if needed ! Note: glissade_calcbwat uses SI units if (main_task .and. verbose_glissade) print*, 'Call glissade_calcbwat' ! convert bwat to SI units for input to glissade_calcbwat - bwat_unscaled(:,:) = model%temper%bwat(:,:) * thk0 + bwat_unscaled(:,:) = model%basal_hydro%bwat(:,:) * thk0 + + !TODO - Move the following calls to a new basal hydrology solver? + + if (model%options%which_ho_bwat == HO_BWAT_FLUX_ROUTING) then + + !WHL - Temporary code for debugging: Make up a simple basal melt field. +! model%basal_hydro%head(:,:) = & +! model%geometry%thck(:,:)*thk0 + (rhow/rhoi)*model%geometry%topg(:,:)*thk0 +! head_max = maxval(model%basal_hydro%head) ! max on local processor +! head_max = parallel_reduce_max(head_max) ! global max +! do j = 1, model%general%nsn +! do i = 1, model%general%ewn +! if (head_max - model%basal_hydro%head(i,j) < 1000.d0) then +!! if (head_max - model%basal_hydro%head(i,j) < 200.d0) then +! bmlt_ground_unscaled(i,j) = 1.0d0/scyr ! units are m/s +! else +! bmlt_ground_unscaled(i,j) = 0.0d0 +! endif +! enddo +! enddo + + ! Compute some masks needed below + + call glissade_get_masks(& + model%general%ewn, model%general%nsn, & + model%parallel, & + model%geometry%thck, model%geometry%topg, & + model%climate%eus, model%numerics%thklim, & + ice_mask, & + floating_mask = floating_mask, & + ocean_mask = ocean_mask) + + ! Compute a mask that sets the domain for flux routing. + ! Cells excluded from the domain are: + ! (1) floating or ocean cells + ! (2) cells at the edge of the global domain + ! (3) ice-free cells in the region where the SMB is overwritten + ! by a prescribed negative value (on the assumption that + ! such cells are supposed to be beyond the ice margin) + ! + ! Note: Cells with bwat_mask = 0 can have bwat_flux > 0 if they receive water + ! from adjacent cells with bwat_mask = 1. + ! But once the flux reaches a cell with bwat_mask = 0, it is not routed further. + ! Thus, the total flux in cells with bwat_mask = 0 should be equal to the + ! total input flux of basal meltwater. + + bwat_mask = 1 ! initially, include the entire domain - call glissade_calcbwat(model%options%which_ho_bwat, & - model%basal_physics, & - dt, & ! s - model%geometry%thck*thk0, & ! m - model%numerics%thklim_temp*thk0, & ! m - bmlt_ground_unscaled, & ! m/s - bwat_unscaled) ! m + where (floating_mask == 1 .or. ocean_mask == 1 .or. & + model%general%global_edge_mask == 1) + bwat_mask = 0 + endwhere + + if (model%options%overwrite_acab /= OVERWRITE_ACAB_NONE .and. & + model%climate%overwrite_acab_value < 0.0d0) then + where (model%climate%overwrite_acab_mask == 1 .and. & + model%geometry%thck < model%numerics%thklim) + bwat_mask = 0 + endwhere + endif + + !WHL - debug - Set mask = 0 where thck = 0 for dome test +! where (model%geometry%thck == 0) +! bwat_mask = 0 +! endwhere + + call parallel_halo(bwat_mask, parallel) + + ! Compute bwat based on a steady-state flux routing scheme + + call glissade_bwat_flux_routing(& + model%general%ewn, model%general%nsn, & + model%numerics%dew*len0, model%numerics%dns*len0, & ! m + model%parallel, & + itest, jtest, rtest, & + model%options%ho_flux_routing_scheme, & + model%geometry%thck*thk0, & ! m + model%geometry%topg*thk0, & ! m + model%numerics%thklim_temp*thk0, & ! m + bwat_mask, & + floating_mask, & + bmlt_ground_unscaled, & ! m/s + bwat_unscaled, & ! m + model%basal_hydro%bwatflx, & ! m^3/s + model%basal_hydro%head) ! m + + else ! simpler basal water options + + call glissade_calcbwat(model%options%which_ho_bwat, & + model%basal_hydro, & + dt, & ! s + model%geometry%thck*thk0, & ! m + model%numerics%thklim_temp*thk0, & ! m + bmlt_ground_unscaled, & ! m/s + bwat_unscaled) ! m + + endif ! convert bmlt and bwat from SI units (m/s and m) to scaled model units model%basal_melt%bmlt_ground(:,:) = bmlt_ground_unscaled(:,:) * tim0/thk0 - model%temper%bwat(:,:) = bwat_unscaled(:,:) / thk0 + model%basal_hydro%bwat(:,:) = bwat_unscaled(:,:) / thk0 ! Update tempunstag as sigma weighted interpolation from temp to layer interfaces do up = 2, model%general%upn-1 @@ -1932,7 +2040,7 @@ subroutine glissade_thermal_solve(model, dt) !------------------------------------------------------------------------ ! Note: bwat is needed in halos to compute effective pressure if which_ho_effecpress = HO_EFFECPRESS_BWAT - call parallel_halo(model%temper%bwat, parallel) + call parallel_halo(model%basal_hydro%bwat, parallel) call t_stopf('glissade_thermal_solve') @@ -1973,6 +2081,7 @@ subroutine glissade_thickness_tracer_solve(model) use glissade_bmlt_float, only: verbose_bmlt_float use glissade_calving, only: verbose_calving use glissade_grid_operators, only: glissade_vertical_interpolate + use glide_stop, only: glide_finalise implicit none @@ -2161,21 +2270,25 @@ subroutine glissade_thickness_tracer_solve(model) model%geomderv%dusrfdew*thk0/len0, model%geomderv%dusrfdns*thk0/len0, & model%velocity%uvel * scyr * vel0, model%velocity%vvel * scyr * vel0, & model%numerics%dt_transport * tim0 / scyr, & + model%numerics%adaptive_cfl_threshold, & model%numerics%adv_cfl_dt, model%numerics%diff_cfl_dt) ! Set the transport timestep. ! The timestep is model%numerics%dt by default, but optionally can be reduced for subcycling + !WHL - debug +! if (main_task) then +! print*, 'Checked advective CFL threshold' +! print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr +! print*, 'adv_cfl_dt =', model%numerics%adv_cfl_dt +! endif + + advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt + if (model%numerics%adaptive_cfl_threshold > 0.0d0) then - !WHL - debug -! if (main_task) then -! print*, 'Check advective CFL threshold' -! print*, 'model dt (yr) =', model%numerics%dt * tim0/scyr -! print*, 'adv_cfl_dt =', model%numerics%adv_cfl_dt -! endif + ! subcycle the transport when advective_cfl exceeds the threshold - advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt if (advective_cfl > model%numerics%adaptive_cfl_threshold) then ! compute the number of subcycles @@ -2188,14 +2301,29 @@ subroutine glissade_thickness_tracer_solve(model) print*, 'Ratio =', advective_cfl / model%numerics%adaptive_cfl_threshold print*, 'nsubcyc =', nsubcyc endif + else nsubcyc = 1 endif dt_transport = model%numerics%dt * tim0 / real(nsubcyc,dp) ! convert to s else ! no adaptive subcycling - nsubcyc = model%numerics%subcyc - dt_transport = model%numerics%dt_transport * tim0 ! convert to s + + advective_cfl = model%numerics%dt*(tim0/scyr) / model%numerics%adv_cfl_dt + + ! If advective_cfl exceeds 1.0, then abort cleanly. Otherwise, set dt_transport and proceed. + ! Note: Usually, it would be enough to write a fatal abort message. + ! The call to glide_finalise was added to allow CISM to finish cleanly when running + ! a suite of automated stability tests, e.g. with the stabilitySlab.py script. + if (advective_cfl > 1.0d0) then + if (main_task) print*, 'advective CFL violation; call glide_finalise and exit cleanly' + call glide_finalise(model, crash=.true.) + stop + else + nsubcyc = model%numerics%subcyc + dt_transport = model%numerics%dt_transport * tim0 ! convert to s + endif + endif !------------------------------------------------------------------------- @@ -3054,6 +3182,14 @@ subroutine glissade_calving_solve(model, init_calving) ! Note: Currently hardwired to include 13 of the 16 ISMIP6 basins. ! Does not include the three largest shelves (Ross, Filchner-Ronne, Amery) + call glissade_get_masks(nx, ny, & + parallel, & + model%geometry%thck*thk0, model%geometry%topg*thk0, & + model%climate%eus*thk0, 0.0d0, & ! thklim = 0 + ice_mask, & + floating_mask = floating_mask, & + land_mask = land_mask) + if (init_calving .and. model%options%expand_calving_mask) then ! Identify basins whose floating ice will be added to the calving mask @@ -3071,14 +3207,6 @@ subroutine glissade_calving_solve(model, init_calving) enddo endif - call glissade_get_masks(nx, ny, & - parallel, & - model%geometry%thck*thk0, model%geometry%topg*thk0, & - model%climate%eus*thk0, 0.0d0, & ! thklim = 0 - ice_mask, & - floating_mask = floating_mask, & - land_mask = land_mask) - if (verbose_calving .and. this_rank==rtest) then print*, ' ' print*, 'initial calving_mask, itest, jtest, rank =', itest, jtest, rtest @@ -3475,6 +3603,12 @@ subroutine glissade_calving_solve(model, init_calving) ! halo updates call parallel_halo(model%geometry%thck, parallel) ! Updated halo values of thck are needed below in calclsrf + ! update the upper and lower surfaces + + call glide_calclsrf(model%geometry%thck, model%geometry%topg, & + model%climate%eus, model%geometry%lsrf) + model%geometry%usrf(:,:) = max(0.d0, model%geometry%thck(:,:) + model%geometry%lsrf(:,:)) + if (verbose_calving .and. this_rank == rtest) then print*, ' ' print*, 'Final calving_thck (m), itest, jtest, rank =', itest, jtest, rtest @@ -3494,14 +3628,17 @@ subroutine glissade_calving_solve(model, init_calving) enddo write(6,*) ' ' enddo + print*, ' ' + print*, 'Final usrf (m):' + do j = jtest+3, jtest-3, -1 + write(6,'(i6)',advance='no') j + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') model%geometry%usrf(i,j) * thk0 + enddo + write(6,*) ' ' + enddo endif ! verbose_calving - ! update the upper and lower surfaces - - call glide_calclsrf(model%geometry%thck, model%geometry%topg, & - model%climate%eus, model%geometry%lsrf) - model%geometry%usrf(:,:) = max(0.d0, model%geometry%thck(:,:) + model%geometry%lsrf(:,:)) - end subroutine glissade_calving_solve !======================================================================= @@ -4852,7 +4989,7 @@ subroutine glissade_cleanup_icefree_cells(model) if (model%geometry%thck_old(i,j) > 0.0d0 .and. model%geometry%thck(i,j) == 0.0d0) then ! basal water - model%temper%bwat(i,j) = 0.0d0 + model%basal_hydro%bwat(i,j) = 0.0d0 ! thermal variables if (model%options%whichtemp == TEMP_INIT_ZERO) then diff --git a/libglissade/glissade_basal_traction.F90 b/libglissade/glissade_basal_traction.F90 index 69f1fd75..805a9110 100644 --- a/libglissade/glissade_basal_traction.F90 +++ b/libglissade/glissade_basal_traction.F90 @@ -440,9 +440,12 @@ subroutine calcbeta (whichbabc, & ! If this factor is not present in the input file, it is set to 1 everywhere. ! Compute beta - ! gn = Glen's n from physcon module - beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/gn - 1.0d0) * & - (speed(:,:) + basal_physics%effecpress_stag(:,:)**gn * big_lambda)**(-1.0d0/gn) + ! Note: Where this equation has powerlaw_m, we used to have Glen's flow exponent n, + ! following the notation of Leguy et al. (2014). + ! Changed to powerlaw_m to be consistent with the Schoof and Tsai laws. + m = basal_physics%powerlaw_m + beta(:,:) = coulomb_c * basal_physics%effecpress_stag(:,:) * speed(:,:)**(1.0d0/m - 1.0d0) * & + (speed(:,:) + basal_physics%effecpress_stag(:,:)**m * big_lambda)**(-1.0d0/m) ! If c_space_factor /= 1.0 everywhere, then multiply beta by c_space_factor if (maxval(abs(basal_physics%c_space_factor_stag(:,:) - 1.0d0)) > tiny(0.0d0)) then @@ -696,7 +699,7 @@ end subroutine calcbeta subroutine calc_effective_pressure (which_effecpress, & ewn, nsn, & - basal_physics, & + basal_physics, basal_hydro, & ice_mask, floating_mask, & thck, topg, & eus, & @@ -723,6 +726,10 @@ subroutine calc_effective_pressure (which_effecpress, & basal_physics ! basal physics object ! includes effecpress, effecpress_stag and various parameters + type(glide_basal_hydro), intent(inout) :: & + basal_hydro ! basal hydro object + ! includes bwat and various parameters + integer, dimension(:,:), intent(in) :: & ice_mask, & ! = 1 where ice is present (thk > thklim), else = 0 floating_mask ! = 1 where ice is present and floating, else = 0 @@ -843,7 +850,7 @@ subroutine calc_effective_pressure (which_effecpress, & if (present(bwat)) then - ! Reduce N where basal water is present. + ! Reduce N where basal water is present, following Bueler % van Pelt (2015). ! The effective pressure decreases from overburden P_0 for bwat = 0 to a small value for bwat = bwat_till_max. ! Note: Instead of using a linear ramp for the variation between overburden and the small value ! (as for the BPMP and BMLT options above), we use the published formulation of Bueler & van Pelt (2015). @@ -855,24 +862,19 @@ subroutine calc_effective_pressure (which_effecpress, & if (bwat(i,j) > 0.0d0) then - relative_bwat = max(0.0d0, min(bwat(i,j)/basal_physics%bwat_till_max, 1.0d0)) + relative_bwat = max(0.0d0, min(bwat(i,j)/basal_hydro%bwat_till_max, 1.0d0)) ! Eq. 23 from Bueler & van Pelt (2015) - basal_physics%effecpress(i,j) = basal_physics%N_0 & - * (basal_physics%effecpress_delta * overburden(i,j) / basal_physics%N_0)**relative_bwat & - * 10.d0**((basal_physics%e_0/basal_physics%C_c) * (1.0d0 - relative_bwat)) + basal_physics%effecpress(i,j) = basal_hydro%N_0 & + * (basal_physics%effecpress_delta * overburden(i,j) / basal_hydro%N_0)**relative_bwat & + * 10.d0**((basal_hydro%e_0/basal_hydro%C_c) * (1.0d0 - relative_bwat)) ! The following line (if uncommented) would implement Eq. 5 of Aschwanden et al. (2016). ! Results are similar to Bueler & van Pelt, but the dropoff in N from P_0 to delta*P_0 begins ! with a larger value of bwat (~0.7*bwat_till_max instead of 0.6*bwat_till_max). !! basal_physics%effecpress(i,j) = basal_physics%effecpress_delta * overburden(i,j) & -!! * 10.d0**((basal_physics%e_0/basal_physics%C_c) * (1.0d0 - relative_bwat)) - - !WHL - Uncomment to try a linear ramp in place of the Bueler & van Pelt relationship. - ! This might lead to smoother variations in N with spatial variation in bwat. -!! basal_physics%effecpress(i,j) = overburden(i,j) * & -!! (basal_physics%effecpress_delta + (1.0d0 - relative_bwat) * (1.0d0 - basal_physics%effecpress_delta)) +!! * 10.d0**((basal_hydro%e_0/basal_hydro%C_c) * (1.0d0 - relative_bwat)) ! limit so as not to exceed overburden basal_physics%effecpress(i,j) = min(basal_physics%effecpress(i,j), overburden(i,j)) @@ -887,6 +889,37 @@ subroutine calc_effective_pressure (which_effecpress, & basal_physics%effecpress = 0.0d0 end where + case(HO_EFFECPRESS_BWAT_RAMP) ! Similar to HO_EFFECPRESS_BWAT, but with a ramp function + + ! Initialize for the case where bwat isn't present, and also for points with bwat == 0 + + basal_physics%effecpress(:,:) = overburden(:,:) + + if (present(bwat)) then + + ! Reduce N where basal water is present. + ! The effective pressure decreases from overburden P_0 for bwat = 0 to a small value for bwat = bwat_till_max. + + do j = 1, nsn + do i = 1, ewn + if (bwat(i,j) > 0.0d0) then + + relative_bwat = max(0.0d0, min(bwat(i,j)/basal_hydro%bwat_till_max, 1.0d0)) + + basal_physics%effecpress(i,j) = overburden(i,j) * & + (basal_physics%effecpress_delta + (1.0d0 - relative_bwat) * (1.0d0 - basal_physics%effecpress_delta)) + + end if + enddo + enddo + + endif ! present(bwat) + + where (floating_mask == 1) + ! set to zero for floating ice + basal_physics%effecpress = 0.0d0 + end where + case(HO_EFFECPRESS_OCEAN_PENETRATION) ! Reduce N for ice grounded below sea level based on connectivity of subglacial water to the ocean diff --git a/libglissade/glissade_basal_water.F90 b/libglissade/glissade_basal_water.F90 index c7f4e476..68c3efd8 100644 --- a/libglissade/glissade_basal_water.F90 +++ b/libglissade/glissade_basal_water.F90 @@ -24,21 +24,33 @@ ! !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -!TODO - Support Jesse's water-routing code (or something similar) in parallel? +!TODO - Test and parallelize Jesse's water-routing code. ! Currently supported only for serial Glide runs, in module glide_bwater.F90 module glissade_basal_water use glimmer_global, only: dp + use glimmer_paramets, only: eps11, eps08 + use glimmer_physcon, only: rhoi, rhow, grav, scyr + use glimmer_log use glide_types + use cism_parallel, only: main_task, this_rank, nhalo, parallel_type, parallel_halo implicit none private - public :: glissade_basal_water_init, glissade_calcbwat + public :: glissade_basal_water_init, glissade_calcbwat, glissade_bwat_flux_routing + +!! logical, parameter :: verbose_bwat = .false. + logical, parameter :: verbose_bwat = .true. + + integer, parameter :: pdiag = 4 ! range for diagnostic prints +!! integer, parameter :: pdiag = 3 ! range for diagnostic prints contains +!============================================================== + subroutine glissade_basal_water_init(model) use glimmer_paramets, only: thk0 @@ -50,14 +62,15 @@ subroutine glissade_basal_water_init(model) ! HO_BWAT_NONE: basal water depth = 0 ! HO_BWAT_CONSTANT: basal water depth = prescribed constant ! HO_BWAT_LOCAL_TILL: local basal till model with prescribed drainage rate + ! HO_BWAT_FLUX_ROUTING: steady-state water routing with flux calculation case(HO_BWAT_CONSTANT) ! Set a constant water thickness where ice is present where (model%geometry%thck > model%numerics%thklim) - model%temper%bwat(:,:) = model%basal_physics%const_bwat / thk0 + model%basal_hydro%bwat(:,:) = model%basal_hydro%const_bwat / thk0 elsewhere - model%temper%bwat(:,:) = 0.0d0 + model%basal_hydro%bwat(:,:) = 0.0d0 endwhere case default @@ -68,9 +81,10 @@ subroutine glissade_basal_water_init(model) end subroutine glissade_basal_water_init +!============================================================== subroutine glissade_calcbwat(which_ho_bwat, & - basal_physics, & + basal_hydro, & dt, & thck, & thklim, & @@ -81,13 +95,12 @@ subroutine glissade_calcbwat(which_ho_bwat, & ! Note: This subroutine assumes SI units. ! Currently, only a few simple options are supported. - use glimmer_physcon, only: rhow, scyr use glide_types integer, intent(in) :: & which_ho_bwat !> basal water options - type(glide_basal_physics), intent(in) :: basal_physics ! basal physics object + type(glide_basal_hydro), intent(inout) :: basal_hydro ! basal hydro object real(dp), intent(in) :: & dt, & !> time step (s) @@ -113,6 +126,7 @@ subroutine glissade_calcbwat(which_ho_bwat, & ! HO_BWAT_NONE: basal water depth = 0 ! HO_BWAT_CONSTANT: basal water depth = prescribed constant ! HO_BWAT_LOCAL_TILL: local basal till model with prescribed drainage rate + ! HO_BWAT_FLUX_ROUTING: steady-state water routing with flux calculation (handled in a different subroutine) case(HO_BWAT_NONE) @@ -122,7 +136,7 @@ subroutine glissade_calcbwat(which_ho_bwat, & ! Use a constant water thickness where ice is present, to force Tbed = Tpmp where (thck > thklim) - bwat(:,:) = basal_physics%const_bwat + bwat(:,:) = basal_hydro%const_bwat elsewhere bwat(:,:) = 0.0d0 endwhere @@ -137,11 +151,11 @@ subroutine glissade_calcbwat(which_ho_bwat, & ! compute new bwat, given source (bmlt) and sink (drainage) ! Note: bmlt > 0 for ice melting. Freeze-on will reduce bwat. - dbwat_dt = bmlt(i,j)*rhoi/rhow - basal_physics%c_drainage/scyr ! convert c_drainage from m/yr to m/s + dbwat_dt = bmlt(i,j)*rhoi/rhow - basal_hydro%c_drainage/scyr ! convert c_drainage from m/yr to m/s bwat(i,j) = bwat(i,j) + dbwat_dt*dt ! limit to the range [0, bwat_till_max] - bwat(i,j) = min(bwat(i,j), basal_physics%bwat_till_max) + bwat(i,j) = min(bwat(i,j), basal_hydro%bwat_till_max) bwat(i,j) = max(bwat(i,j), 0.0d0) enddo @@ -151,4 +165,2225 @@ subroutine glissade_calcbwat(which_ho_bwat, & end subroutine glissade_calcbwat +!============================================================== + + subroutine glissade_bwat_flux_routing(& + nx, ny, & + dx, dy, & + parallel, & + itest, jtest, rtest, & + flux_routing_scheme, & + thck, topg, & + thklim, & + bwat_mask, floating_mask, & + bmlt, bwat, & + bwatflx, head) + + ! This subroutine is a recoding of Jesse Johnson's steady-state water routing scheme in Glide. + ! It has been parallelized for Glissade. + + use cism_parallel, only: tasks ! while code is serial only + + ! Input/output arguments + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid cell size (m) + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + integer, intent(in) :: & + flux_routing_scheme ! flux routing scheme: D8, Dinf or FD8; see subroutine route_basal_water + + real(dp), dimension(nx,ny), intent(in) :: & + thck, & ! ice thickness (m) + topg, & ! bed topography (m) + bmlt ! basal melt rate (m/s) + + real(dp), intent(in) :: & + thklim ! minimum ice thickness for basal melt and hydropotential calculations (m) + ! Note: This is typically model%geometry%thklim_temp + + integer, dimension(nx,ny), intent(in) :: & + bwat_mask, & ! mask to identify cells through which basal water is routed; + ! = 0 for floating and ocean cells; cells at global domain edge; + ! and cells with thck = 0 and forced negative SMB + floating_mask ! = 1 if ice is present (thck > thklim) and floating, else = 0 + + + real(dp), dimension(nx,ny), intent(inout) :: & + bwat ! basal water depth (m) + + real(dp), dimension(nx,ny), intent(out) :: & + bwatflx, & ! basal water flux (m^3/s) + head ! hydraulic head (m) + + ! Local variables + + integer :: i, j, p + + !TODO - Make effecpress in/out? + real(dp), dimension(nx, ny) :: & + effecpress, & ! effective pressure + lakes ! difference between filled head and original head (m) + + ! parameters related to effective pressure + real(dp), parameter :: & + c_effective_pressure = 0.0d0 ! for now estimated as N = c/bwat + + ! parameters related to subglacial fluxes + ! The basal water flux is given by Sommers et al. (2018), Eq. 5: + ! + ! q = (b^3*g)/[(12*nu)*(1 + omega*Re)] * (-grad(h)) + ! + ! where q = basal water flux per unit width (m^2/s) = bwatflx/dx + ! b = water depth (m) = bwat + ! g = gravitational constant (m/s^2) = grad + ! nu = kinematic viscosity of water (m^2/s)= visc_water + ! omega = parameter controlling transition between laminar and turbulent flow + ! Re = Reynolds number (large for turbulent flow) + ! h = hydraulic head (m) + ! + ! By default, we set Re = 0, which means the flow is purely laminar, as in Sommers et al. (2018), Eq. 6. + + ! Optionally, one or more of these parameters could be made a config parameter in the basal_hydro type + real(dp), parameter :: & + visc_water = 1.787e-6, & ! kinematic viscosity of water (m^2/s); Sommers et al. (2018), Table 2 + omega_hydro = 1.0d-3, & ! omega (unitless) in Sommers et al (2018), Eq. 6 + Reynolds = 0.0d0 ! Reynolds number (unitless), = 0 for pure laminar flow + + real(dp), parameter :: & + c_flux_to_depth = 1.0d0/((12.0d0*visc_water)*(1.0d0 + omega_hydro*Reynolds)), & ! proportionality coefficient in Eq. 6 + p_flux_to_depth = 2.0d0, & ! exponent for water depth; = 2 if q is proportional to b^3 + q_flux_to_depth = 1.0d0 ! exponent for potential gradient; = 1 if q is linearly proportional to grad(h) + + + ! WHL - debug fix_flats subroutine + logical :: test_fix_flats = .false. +!! logical :: test_fix_flats = .true. + integer :: nx_test, ny_test + real(dp), dimension(:,:), allocatable :: phi_test + integer, dimension(:,:), allocatable :: mask_test + + !WHL - debug + !Note: This test works in serial, but does not work with parallel updates. + ! To use it again, would need to comment out parallel calls in fix_flats. + if (test_fix_flats) then + + ! Solve the example problem of Garbrecht & Martz (1997) + ! Their problem is 7x7, but easier to solve if padded with a row of low numbers. + + nx_test = 9 + ny_test = 9 + allocate (phi_test(nx_test,ny_test)) + allocate (mask_test(nx_test,ny_test)) + + mask_test = 1 + do j = 1, ny_test + do i = 1, nx_test + if (i == 1 .or. i == nx_test .or. j == 1 .or. j == ny_test) then + mask_test(i,j) = 0 + endif + enddo + enddo + + phi_test(:,9) = (/ 1.d0, 1.d0,1.d0,1.d0,1.d0,1.d0,1.d0,1.d0, 1.d0 /) + phi_test(:,8) = (/ 1.d0, 9.d0,9.d0,9.d0,9.d0,9.d0,9.d0,9.d0, 1.d0 /) + phi_test(:,7) = (/ 1.d0, 9.d0,6.d0,6.d0,6.d0,6.d0,6.d0,9.d0, 1.d0 /) + phi_test(:,6) = (/ 1.d0, 8.d0,6.d0,6.d0,6.d0,6.d0,6.d0,9.d0, 1.d0 /) + phi_test(:,5) = (/ 1.d0, 8.d0,6.d0,6.d0,6.d0,6.d0,6.d0,9.d0, 1.d0 /) + phi_test(:,4) = (/ 1.d0, 7.d0,6.d0,6.d0,6.d0,6.d0,6.d0,8.d0, 1.d0 /) + phi_test(:,3) = (/ 1.d0, 7.d0,6.d0,6.d0,6.d0,6.d0,6.d0,8.d0, 1.d0 /) + phi_test(:,2) = (/ 1.d0, 7.d0,7.d0,5.d0,7.d0,7.d0,8.d0,8.d0, 1.d0 /) + phi_test(:,1) = (/ 1.d0, 1.d0,1.d0,1.d0,1.d0,1.d0,1.d0,1.d0, 1.d0 /) + + call fix_flats(& + nx_test, ny_test, & + parallel, & + 5, 5, rtest, & + phi_test, & + mask_test) + + deallocate(phi_test, mask_test) + + endif + + !WHL - debug + if (this_rank == rtest) then + print*, 'In glissade_bwat_flux_routing: rtest, itest, jtest =', rtest, itest, jtest + endif + + ! Uncomment if the following fields are not already up to date in halo cells +! call parallel_halo(thk, parallel) +! call parallel_halo(topg, parallel) + call parallel_halo(bwat, parallel) + call parallel_halo(bmlt, parallel) + + ! Compute effective pressure N as a function of water depth + + call effective_pressure(& + bwat, & + c_effective_pressure, & + effecpress) + + ! Compute the hydraulic head + + call compute_head(& + nx, ny, & + thck, & + topg, & + effecpress, & + thklim, & + floating_mask, & + head) + + p = pdiag + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'thck (m):' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') thck(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'topg (m):' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') topg(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'effecpress (Pa):' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') effecpress(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'bmlt (m/yr):' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') bmlt(i,j) * scyr + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'bwat_mask:' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') bwat_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'Before fill: head (m):' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') head(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! Route basal water down the gradient of hydraulic head, giving a water flux + + call route_basal_water(& + nx, ny, & + dx, dy, & + parallel, & + itest, jtest, rtest, & + flux_routing_scheme, & + head, & + bmlt, & + bwat_mask, & + bwatflx, & + lakes) + + ! Convert the water flux to a basal water depth + + call flux_to_depth(& + nx, ny, & + dx, dy, & + itest, jtest, rtest, & + bwatflx, & + head, & + c_flux_to_depth, & + p_flux_to_depth, & + q_flux_to_depth, & + bwat_mask, & + bwat) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + write(6,*) 'bwatflx (m^3/s):' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') bwatflx(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'bwat (mm):' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') bwat(i,j) * 1000.d0 + enddo + write(6,*) ' ' + enddo + endif + + end subroutine glissade_bwat_flux_routing + +!============================================================== + + subroutine effective_pressure(& + bwat, & + c_effective_pressure, & + effecpress) + + ! Compute the effective pressure: the part of ice overburden not balanced by water pressure + ! TODO: Try c_effective_pressure > 0, or call calc_effecpress instead + + real(dp),dimension(:,:),intent(in) :: bwat ! water depth + real(dp) ,intent(in) :: c_effective_pressure ! constant of proportionality + real(dp),dimension(:,:),intent(out) :: effecpress ! effective pressure + + ! Note: By default, c_effective_pressure = 0 + ! This implies N = 0; full support of the ice by water at the bed + ! Alternatively, could call the standard glissade subroutine, calc_effective_pressure + + where (bwat > 0.d0) + effecpress = c_effective_pressure / bwat + elsewhere + effecpress = 0.d0 + endwhere + + end subroutine effective_pressure + +!============================================================== + + subroutine compute_head(& + nx, ny, & + thck, & + topg, & + effecpress, & + thklim, & + floating_mask, & + head) + + ! Compute the hydraulic head as the bed elevation plus the scaled water pressure: + ! + ! head = z_b + p_w / (rhow*g) + ! + ! where z_b = bed elevation (m) = topg + ! p_w = water pressure (Pa) = p_i - N + ! p_i = ice overburden pressure = rhoi*g*H + ! N = effective pressure (Pa) = part of overburden not supported by water + ! H = ice thickness (m) + ! + ! If we make the approximation p_w =~ p_i, then + ! + ! head =~ z_b + (rhoi/rhow) * H + + implicit none + + ! Input/output variables + + integer, intent(in) :: & + nx, ny ! number of grid cells in each direction + + real(dp), dimension(nx,ny), intent(in) :: & + thck, & ! ice thickness (m) + topg, & ! bed elevation (m) + effecpress ! effective pressure (Pa) + + real(dp), intent(in) :: & + thklim ! minimum ice thickness for bmlt and head calculations + + integer, dimension(nx,ny), intent(in) :: & + floating_mask ! = 1 if ice is present (thck > thklim) and floating, else = 0 + + real(dp), dimension(nx,ny), intent(out) :: & + head ! hydraulic head (m) + + where (thck > thklim .and. floating_mask /= 1) + head = topg + (rhoi/rhow)*thck - effecpress/(rhow*grav) + elsewhere + head = max(topg, 0.0d0) + endwhere + + end subroutine compute_head + +!============================================================== + + subroutine route_basal_water(& + nx, ny, & + dx, dy, & + parallel, & + itest, jtest, rtest, & + flux_routing_scheme, & + head, & + bmlt, & + bwat_mask, & + bwatflx, & + lakes) + + ! Route water from the basal melt field to its destination, recording the water flux along the way. + ! Water flow direction is determined according to the gradient of the hydraulic head. + ! For the algorithm to work correctly, surface depressions must be filled, + ! so that all cells have an outlet to the ice sheet margin. + ! This results in the lakes field, which is the difference between the filled head and the original head. + ! The method used is by Quinn et. al. (1991). + ! + ! Based on code by Jesse Johnson (2005), adapted from the glimmer_routing file by Ian Rutt. + + use cism_parallel, only: parallel_global_sum + + !WHL - debug + use cism_parallel, only: parallel_globalindex, parallel_reduce_max + + implicit none + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid cell size (m) + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + integer, intent(in) :: & + flux_routing_scheme ! flux routing scheme: D8, Dinf or FD8 + + real(dp), dimension(nx,ny), intent(in) :: & + bmlt ! basal melt beneath grounded ice (m/s) + + real(dp), dimension(nx,ny), intent(inout) :: & + head ! hydraulic head (m) + ! intent inout because it can be filled and adjusted below + + integer, dimension(nx,ny), intent(in) :: & + bwat_mask ! mask to identify cells through which basal water is routed; + ! = 1 where ice is present and not floating + + real(dp), dimension(nx,ny), intent(out) :: & + bwatflx, & ! water flux through a grid cell (m^3/s) + lakes ! lakes field, difference between filled and original head + + ! Local variables + + integer :: nlocal ! number of locally owned cells + integer :: count, count_max ! iteration counters + integer :: i, j, k, ii, jj, ip, jp, p + integer :: i1, j1, i2, j2, itmp, jtmp + + logical :: finished ! true when an iterative loop has finished + + integer, dimension(:,:), allocatable :: & + sorted_ij ! i and j indices of all cells, sorted from low to high values of head + + real(dp), dimension(-1:1,-1:1,nx,ny) :: & + flux_fraction, & ! fraction of flux from each cell that flows downhill to each of 8 neighbors + bwatflx_halo ! water flux (m^3/s) routed to a neighboring halo cell; routed further in next iteration + + real(dp), dimension(nx,ny) :: & + head_filled, & ! head after depressions are filled (m) + bwatflx_accum, & ! water flux (m^3/s) accumulated over multiple iterations + sum_bwatflx_halo ! bwatflx summed over the first 2 dimensions in each grid cell + + integer, dimension(nx,ny) :: & + local_mask, & ! = 1 for cells owned by the local processor, else = 0 + halo_mask, & ! = 1 for the layer of halo cells adjacent to locally owned cells, else = 0 + margin_mask ! = 1 for cells at the grounded ice margin, as defined by bwat_mask, else = 0 + + real(dp) :: & + total_flux_in, & ! total input flux (m^3/s), computed as sum of bmlt*dx*dy + total_flux_out, & ! total output flux (m^3/s), computed as sum of bwatflx at ice margin + err, & ! water conservation error + global_flux_sum ! flux sum over all cells in global domain + + character(len=100) :: message + + !WHL - debug + real(dp) :: bmlt_max, bmlt_max_global + integer :: imax, jmax, rmax, iglobal, jglobal + ! Allocate the sorted_ij array + + nlocal = parallel%own_ewn * parallel%own_nsn + allocate(sorted_ij(nlocal,2)) + + ! Compute mask of locally owned and halo cells. + ! These masks are used to transfer fluxes between processors on subsequent iterations. + + local_mask = 0 + halo_mask = 0 + do j = nhalo, ny-nhalo+1 + do i = nhalo, nx-nhalo+1 + if (j == nhalo .or. j == ny-nhalo+1 .or. i == nhalo .or. i == nx-nhalo+1) then + halo_mask(i,j) = 1 + elseif (j > nhalo .or. j <= ny-nhalo .or. i > nhalo .or. i <= nx-nhalo+1) then + local_mask(i,j) = 1 + endif + enddo + enddo + + ! Fill depressions in the head field, so that no interior cells are sinks + + call fill_depressions(& + nx, ny, & + parallel, & + itest, jtest, rtest, & + head, & + bwat_mask, & + head_filled) + + ! Note: In an earlier code version, fix_flats was called here. + ! It is not needed, however, if the fill_depressions scheme is run with epsilon > 0. + + ! Compute the lake depth + lakes = head_filled - head + + ! Update head with the filled values + head = head_filled + + p = pdiag + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'After fill: head (m):' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') head(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'lakes (m):' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') lakes(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! Sort heights. + ! The sorted_ij array stores the i and j index for each locally owned cell, from lowest to highest value. + + call sort_heights(& + nx, ny, nlocal, & + itest, jtest, rtest, & + head, sorted_ij) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'sorted, from the top:' + do k = nlocal, nlocal-10, -1 + i = sorted_ij(k,1) + j = sorted_ij(k,2) + print*, k, i, j, head(i,j) + enddo + endif + + call get_flux_fraction(& + nx, ny, nlocal, & + dx, dy, & + itest, jtest, rtest, & + flux_routing_scheme, & + sorted_ij, & + head, & + bwat_mask, & + flux_fraction) + + ! Initialize bwatflx in locally owned cells with the basal melt, which will be routed downslope. + ! Multiply by area, so units are m^3/s. + ! The halo water flux, bwatflx_halo, holds water routed to halo cells; + ! it will be routed downhill on the next iteration. + ! The accumulated flux, bwatflx_accum, holds the total flux over multiple iterations. + ! Note: This subroutine conserves water only if bmlt >= 0 everywhere. + ! One way to account for refreezing would be to do the thermal calculation after + ! computing bwat in this subroutine. At that point, refreezing would take away + ! from the bwat computed here. In the next time step, positive values of bmlt + ! would provide a new source for bwat. + ! In other words, the sequence would be: + ! (1) Ice transport and calving + ! (2) Basal water routing: apply bmlt and diagnose bwat + ! (3) Vertical heat flow: + ! (a) compute bmlt + ! (b) use bmlt < 0 to reduce bwat + ! (c) save bmlt > 0 for the next time step (and write to restart) + ! (4) Diagnose velocity + + bwatflx = 0.0d0 + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + bwatflx(i,j) = bmlt(i,j) * dx * dy + bwatflx(i,j) = max(bwatflx(i,j), 0.0d0) ! not conservative unless refreezing is handled elsewhere + enddo + enddo + bwatflx_halo = 0.0d0 + bwatflx_accum = 0.0d0 + + ! Compute total input of meltwater (m^3/s) + total_flux_in = parallel_global_sum(bwatflx, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'Total input basal melt flux (m^3/s):', total_flux_in + endif + + ! Loop over locally owned cells, from highest to lowest. + ! During each iteration, there are two possible outcomes for routing: + ! (1) Routed to the ice sheet margin, to a cell with bwat_mask = 0. + ! In this case, the routing of that flux is done. + ! (2) Routed to a halo cell, i.e. a downslope cell on a neighboring processor. + ! In this case, the flux will be routed further downhill on the next iteration. + ! When all the water has been routed to the margin, we are done. + + count = 0 + ! Note: It is hard to predict how many iterations will be sufficient. + ! With Dinf or FD8, we can have flow back and forth across processor boundaries, + ! requiring many iterations to reach the margin. + ! For Greenland 4 km, Dinf requires ~20 iterations on 4 cores, and FD8 can require > 40. + count_max = 50 + finished = .false. + + do while (.not.finished) + + count = count + 1 + if (verbose_bwat .and. this_rank == rtest) then + print*, 'flux routing, count =', count + endif + + do k = nlocal, 1, -1 + + ! Get i and j indices of current point + i = sorted_ij(k,1) + j = sorted_ij(k,2) + + ! Apportion the flux among downslope neighbors + if (bwat_mask(i,j) == 1 .and. bwatflx(i,j) > 0.0d0) then + do jj = -1,1 + do ii = -1,1 + ip = i + ii + jp = j + jj + if (flux_fraction(ii,jj,i,j) > 0.0d0) then + if (halo_mask(ip,jp) == 1) then + bwatflx_halo(ii,jj,i,j) = bwatflx(i,j)*flux_fraction(ii,jj,i,j) + elseif (local_mask(ip,jp) == 1) then + bwatflx(ip,jp) = bwatflx(ip,jp) + bwatflx(i,j)*flux_fraction(ii,jj,i,j) + endif + endif ! flux_fraction > 0 + enddo + enddo + endif + + enddo ! loop from high to low + + ! Accumulate bwatflx from the latest iteration. + ! Reset to zero for the next iteration, if needed. + + bwatflx_accum = bwatflx_accum + bwatflx + bwatflx = 0.0d0 + + ! If bwatflx_halo = 0 everywhere, then we are done. + ! (If the remaining flux is very small (< eps11), discard it to avoid + ! unnecessary extra iterations.) + ! If bwatflx_halo remains, then communicate it to neighboring tasks and + ! continue routing on the next iteration. + + do j = 1, ny + do i = 1, nx + sum_bwatflx_halo(i,j) = sum(bwatflx_halo(:,:,i,j)) +!! if (verbose_bwat .and. sum_bwatflx_halo(i,j) > 0.0d0) then + if (verbose_bwat .and. sum_bwatflx_halo(i,j) > eps11 .and. count > 10) then + print*, 'Nonzero bwatflx_halo, count, rank, i, j, sum_bwatflx_halo:', & + count, this_rank, i, j, sum_bwatflx_halo(i,j) + call parallel_globalindex(i, j, iglobal, jglobal, parallel) + print*, ' iglobal, jglobal:', iglobal, jglobal + endif + enddo + enddo + global_flux_sum = parallel_global_sum(sum_bwatflx_halo, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, 'Before halo update, sum of bwatflx_halo:', global_flux_sum + print*, ' ' + print*, 'sum_bwatflx_halo:' + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(e10.3)',advance='no') sum_bwatflx_halo(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'rank, i, j, bwatflx_halo:' + do j = jtest+1, jtest + do i = itest-4, itest + 4 + write(6, '(3i5,9e10.3)') this_rank, i, j, bwatflx_halo(:,:,i,j) + enddo + enddo + endif + + if (global_flux_sum > eps11) then + + finished = .false. + + ! Communicate bmltflx_halo to the halo cells of neighboring processors + call parallel_halo(bwatflx_halo(:,:,:,:), parallel) + + ! bmltflx_halo is now available in the halo cells of the local processor. + ! Route downslope to the adjacent locally owned cells. + ! These fluxes will be routed further downslope during the next iteration. + + do j = 2, ny-1 + do i = 2, nx-1 + if (halo_mask(i,j) == 1 .and. sum(bwatflx_halo(:,:,i,j)) > 0.0d0) then + do jj = -1,1 + do ii = -1,1 + if (bwatflx_halo(ii,jj,i,j) > 0.0d0) then + ip = i + ii + jp = j + jj + if (local_mask(ip,jp) == 1) then + bwatflx(ip,jp) = bwatflx(ip,jp) + bwatflx_halo(ii,jj,i,j) + if (verbose_bwat) then +!!! print*, 'Nonzero bwatflx, rank, i, j:', this_rank, ip, jp, bwatflx(ip,jp) + endif + endif + endif ! bwatflx_halo > 0 to this local cell + enddo ! ii + enddo ! jj + endif ! bwatflx_halo > 0 from this halo cell + enddo ! i + enddo ! j + + ! Reset bwatflx_halo for the next iteration + bwatflx_halo = 0.0d0 + + global_flux_sum = parallel_global_sum(bwatflx, parallel) + if (verbose_bwat .and. this_rank == rtest) then + ! Should be equal to the global sum of bwatflx_halo computed above + print*, 'After halo update, sum(bwatflx) =', global_flux_sum + endif + + else ! bwatflx_halo = 0 everywhere; no fluxes to route to adjacent processors + if (verbose_bwat .and. this_rank == rtest) print*, 'Done routing fluxes' + finished = .true. + bwatflx = bwatflx_accum + endif + + if (count > count_max) then + call write_log('Hydrology error: too many iterations in route_basal_water', GM_FATAL) + endif + + enddo ! finished routing + + ! Identify cells just beyond the ice sheet margin, which can receive from upstream but not send downstream + where (bwat_mask == 0 .and. bwatflx > 0.0d0) + margin_mask = 1 + elsewhere + margin_mask = 0 + endwhere + + ! Compute total output of meltwater (m^3/s) and check that input = output, within roundoff. + + total_flux_out = parallel_global_sum(bwatflx*margin_mask, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, 'Total output basal melt flux (m^3/s):', total_flux_out + print*, 'Difference between input and output =', total_flux_in - total_flux_out + endif + + ! Not sure if a threshold of eps11 is large enough. Increase if needed. + if (total_flux_in > 0.0d0) then + err = abs(total_flux_in - total_flux_out) + if (err > eps11) then +! write(message,*) 'Hydrology error: total water not conserved, error =', err +! call write_log(message, GM_FATAL) + write(message,*) 'WARNING: Hydrology error: total water not conserved, error =', err + call write_log(message, GM_WARNING) + endif + endif + + ! clean up + deallocate(sorted_ij) + + end subroutine route_basal_water + +!============================================================== + + subroutine flux_to_depth(& + nx, ny, & + dx, dy, & + itest, jtest, rtest, & + bwatflx, & + head, & + c_flux_to_depth, & + p_flux_to_depth, & + q_flux_to_depth, & + bwat_mask, & + bwat) + + ! Assuming that the flow is steady state, this function simply solves + ! flux = depth * velocity + ! for the depth, assuming that the velocity is a function of depth, + ! and pressure potential. This amounts to assuming a Weertman film, + ! or Manning flow, both of which take the form of a constant times water + ! depth to a power, times grad(head) to a power. + + use glissade_grid_operators, only: glissade_gradient_at_edges + + ! Input/ouput variables + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid spacing in each direction + + real(dp), dimension(nx,ny), intent(in) :: & + bwatflx, & ! basal water flux (m^3/s) + head ! hydraulic head (m) + + real(dp), intent(in) :: & + c_flux_to_depth, & ! constant of proportionality + p_flux_to_depth, & ! exponent for water depth + q_flux_to_depth ! exponent for potential_gradient + + integer, dimension(nx,ny), intent(in) :: & + bwat_mask ! mask to identify cells through which basal water is routed; + ! = 1 where ice is present and not floating + + real(dp), dimension(nx,ny), intent(out):: & + bwat ! water depth + + ! Local variables + + integer :: i, j, p + + real(dp), dimension(nx,ny) :: & + grad_head ! gradient of hydraulic head (m/m), averaged to cell centers + + real(dp), dimension(nx-1,ny) :: & + dhead_dx ! gradient component on E edges + + real(dp), dimension(nx,ny-1) :: & + dhead_dy ! gradient component on N edges + + real(dp) :: & + dhead_dx_ctr, dhead_dy_ctr, & ! gradient components averaged to cell centers + p_exponent ! p-dependent exponent in bwat expression + + integer, dimension(nx,ny) :: & + ice_mask ! mask passed to glissade_gradient_at edges; = 1 everywhere + + ice_mask = 1 + + ! Compute gradient components at cell edges + ! HO_GRADIENT_MARGIN_LAND: Use all field values when computing the gradient, including values in ice-free cells. + + call glissade_gradient_at_edges(& + nx, ny, & + dx, dy, & + head, & + dhead_dx, dhead_dy, & + ice_mask, & + gradient_margin_in = HO_GRADIENT_MARGIN_LAND) + + grad_head = 0.0d0 ! will remain 0 in outer row of halo cells + do j = 2, ny-1 + do i = 2, nx-1 + dhead_dx_ctr = 0.5d0 * (dhead_dx(i-1,j) + dhead_dx(i,j)) + dhead_dy_ctr = 0.5d0 * (dhead_dy(i,j-1) + dhead_dy(i,j)) + grad_head(i,j) = sqrt(dhead_dx_ctr**2 + dhead_dy_ctr**2) + enddo + enddo + + !TODO - If a halo update is needed for grad_head, then pass in 'parallel'. But may not be needed. +!! call parallel_halo(grad_head, parallel) + + !WHL - debug + p = 5 + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'grad(head):' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(e10.3)',advance='no') grad_head(i,j) + enddo + write(6,*) ' ' + enddo + endif + + p_exponent = 1.d0 / (p_flux_to_depth + 1.d0) + + ! Note: In Sommers et al. (2018), Eq. 5, the basal water flux q (m^2/s) is + ! q = (b^3 * g) / [(12*nu)(1 + omega*Re)] * (-grad(h)) + ! where nu = kinematic viscosity of water = 1.787d-06 m^2/s + ! omega = 0.001 + ! Re = Reynolds number + ! + ! Following Aleah's formulation: + ! F = b^3 * c * g * dx * -grad(h) where c = 1/[(12*nu)(1 + omega*Re)] + ! b^3 = F / [c * g * dx * -grad(h)] + ! b = { F / [c * g * dx * -grad(h)] }^(1/3) + ! + ! In the context of a formulation with general exponents, + ! we have q_flux_to_depth = 1 and p_flux_to_depth = 2 (so p_exponent = 1/3) + ! + ! Jesse's Glimmer code has this: + ! bwat = ( bwatflx / (c_flux_to_depth * scyr * dy * grad_wphi**q_flux_to_depth) ) ** exponent + ! which is missing the grav term and seems to have an extra scyr term. + ! Also, c_flux_to_depth = 1 / (12 * 1.6d-3) in Jesse's code. Note exponent of d-3 instead of d-6 for nu. + ! + ! Note: Assuming dx = dy + ! TODO: Modify for the case dx /= dy? + + where (grad_head /= 0.d0 .and. bwat_mask == 1) + bwat = ( bwatflx / (c_flux_to_depth * grav * dy * grad_head**q_flux_to_depth) ) ** p_exponent + elsewhere + bwat = 0.d0 + endwhere + + end subroutine flux_to_depth + +!============================================================== + + subroutine fill_depressions(& + nx, ny, & + parallel, & + itest, jtest, rtest, & + phi_in, & + phi_mask, & + phi) + + ! Fill depressions in the input field, phi_in. + ! The requirements for the output field, phi_out, are: + ! (1) phi_out >= phi_in everywhere + ! (2) For each cell with phi_mask = 1, there is a descending path to the boundary. + ! That is, phi1 >= phi2 for any two adjacent cells along the path, where the flow + ! is from cell 1 to cell 2. + ! (3) phi_out is the lowest surface consistent with properties (1) and (2). + ! + ! The algorithm is based on this paper: + ! Planchon, O., and F. Darboux (2001): A fast, simple and versatile algorithm + ! to fill the depressions of digital elevation models, Catena (46), 159-176. + ! + ! The basic idea is: + ! Let phi = the current best guess for phi_out. + ! Initially, set phi = phi_in on the boundary, and set phi = a large number elsewhere. + ! Loop through the domain. For each cell c, with value phi(c) not yet fixed as a known value, + ! compute phi_min8(n), the current minimum of phi in the 8 neighbor cells. + ! If phi_in(c) > phi_min8(n) + eps, then set phi(c) = phi_in(c) and mark that cell as having a known value, + ! since phi(c) cannot go any lower. Here, eps is a small number greater than zero. + ! If phi_in(c) < phi_min8(n) + eps, but phi(c) > phi_min8(c) + eps, set phi(c) = phi_min8(n) + eps. + ! Do not mark the cell as having a known value, because it might be lowered further. + ! Continue until no further lowering of phi is possible. At that point, phi = phi_out. + ! Note: Setting eps = 0 would result in flat surfaces that would need to be fixed later. + + use cism_parallel, only: parallel_reduce_sum + use cism_parallel, only: parallel_globalindex + + implicit none + + ! Input/output variables + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + real(dp), dimension(nx,ny), intent(in) :: & + phi_in ! input field with depressions to be filled + + integer, dimension(nx,ny), intent(in) :: & + phi_mask ! = 1 in the domain where depressions need to be filled, else = 0 + + real(dp), dimension(nx,ny), intent(out) :: & + phi ! output field with depressions filled + + ! Local variables -------------------------------------- + + logical, dimension(nx,ny) :: & + known ! = true for cells where the final phi(i,j) is known + + integer :: & + local_lowered, & ! local sum of cells where phi is lowered + global_lowered ! global sum of cells where phi is lowered + + real(dp) :: & + phi_min8 ! current minval of phi in a cell's 8 neighbors, + ! considering only cells with phi_mask = 1 + + real(dp) :: epsilon ! small increment in phi, either epsilon_edge or epsilon_diag + + logical :: finished ! true when an iterative loop has finished + + integer :: count ! iteration counter + + integer :: i, j, ii, jj, ip, jp, p + integer :: iglobal, jglobal + integer :: i1, i2, istep, j1, j2, jstep + + real(dp), parameter :: big_number = 1.d+20 ! initial large value for phi + + ! According to Planchon & Darboux (2001), there should be one value of epsilon for edge neighbors + ! and another value for corner neighbors. + real(dp), parameter :: & + epsilon_edge = 1.d-4, & ! small increment in phi to avoid flat regions, applied to edge neighbors + epsilon_diag = 1.d-4*sqrt(2.d0) ! small increment in phi to avoid flat regions, applied to diagonal neighbors + + !WHL - Typically, it takes ~10 iterations to fill all depressions on a large domain. + integer, parameter :: count_max = 100 + +!! logical, parameter :: verbose_depression = .false. + logical, parameter :: verbose_depression = .true. + + ! Initial halo update, in case phi_in is not up to date in halo cells + call parallel_halo(phi_in, parallel) + + ! Initialize phi to a large value + where (phi_mask == 1) + phi = big_number + known = .false. + elsewhere + phi = 0.0d0 + known = .true. + endwhere + + ! Set phi = phi_in for boundary cells. + ! A boundary cell is a cell with phi_mask = 1, adjacent to one or more cells with phi_mask = 0. + do j = 2, ny-1 + do i = 2, nx-1 + if (phi_mask(i,j) == 1) then + if (phi_mask(i-1,j+1)==0 .or. phi_mask(i,j+1)==0 .or. phi_mask(i+1,j+1)==0 .or. & + phi_mask(i-1,j) ==0 .or. phi_mask(i+1,j) ==0 .or. & + phi_mask(i-1,j-1)==0 .or. phi_mask(i,j-1)==0 .or. phi_mask(i+1,j-1)==0) then + phi(i,j) = phi_in(i,j) + known(i,j) = .true. + endif + endif + enddo + enddo + + ! The resulting mask applies to locally owned cells and one layer of halo cells. + ! A halo update brings it up to date in all halo cells. + + call parallel_halo(phi, parallel) + + p = pdiag + + if (verbose_depression .and. this_rank == rtest) then + print*, ' ' + print*, 'Initial phi:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(e11.4)',advance='no') phi(i,j) + enddo + write(6,*) ' ' + enddo + endif + + count = 0 + finished = .false. + + do while (.not.finished) + + count = count + 1 + local_lowered = 0 + + if (verbose_depression .and. this_rank == rtest) then + write(6,*) ' ' + print*, 'fill_depressions, count =', count + endif + + ! Loop through cells + ! Iterate until phi cannot be lowered further. + ! + ! To vary the route through the cells and reduce the required number of iterations, + ! we alternate between four possible sequences: + ! (1) j lo to hi, i lo to hi + ! (2) j hi to lo, i hi to lo + ! (3) j lo to hi, i hi to lo + ! (4) j hi to lo, i lo to hi + ! Other sequences would be possible with i before j, but these are not Fortran-friendly. + + if (mod(count,4) == 1) then + j1 = 2; j2 = ny-1; jstep = 1 + i1 = 2; i2 = nx-1; istep = 1 + elseif (mod(count,4) == 2) then + j1 = ny-1; j2 = 2; jstep = -1 + i1 = nx-1; i2 = 2; istep = -1 + elseif (mod(count,4) == 3) then + j1 = 2; j2 = ny-1; jstep = 1 + i1 = nx-1; i2 = 2; istep = -1 + elseif (mod(count,4) == 0) then + j1 = ny-1; j2 = 2; jstep = -1 + i1 = 2; i2 = nx-1; istep = 1 + endif + + do j = j1, j2, jstep + do i = i1, i2, istep + if (phi_mask(i,j) == 1 .and. .not.known(i,j)) then + + ! In each cell, compute the min value of phi in the 8 neighbors + phi_min8 = big_number + do jj = -1,1 + do ii = -1,1 + ! If this is the centre point, ignore + if (ii == 0 .and. jj == 0) then + continue + else ! check whether this neighbor has the minimum phi value + ip = i + ii + jp = j + jj + if (phi(ip,jp) < phi_min8) phi_min8 = phi(ip,jp) + if (mod(ii+jj,2) == 0) then ! diagonal neighbor + epsilon = epsilon_diag + else ! edge neighbor + epsilon = epsilon_edge + endif + endif + enddo + enddo + + ! If phi_in(i,j) > phi_min8 + eps, set phi(i,j) = phi_in(i,j); mark cell as known. + ! Else if phi(i,j) > phi_min8 + eps, set phi(i,j) = phi_min8 + eps; do not mark as known. + ! Note: epsilon could be either epsilon_edge or epsilon_diag. + + if (phi_in(i,j) > phi_min8 + epsilon) then + + !WHL - debug + if (verbose_depression .and. count >= 20) then + call parallel_globalindex(i, j, iglobal, jglobal, parallel) + print*, ' ' + print*, 'rank, i, j, ig, jg:', this_rank, i, j, iglobal, jglobal + print*, ' phi_in, phi:', phi_in(i,j), phi(i,j) + print*, ' phi_min8 =', phi_min8 + print*, ' new phi = phi_in' + endif + + phi(i,j) = phi_in(i,j) + known(i,j) = .true. + local_lowered = local_lowered + 1 + + elseif (phi(i,j) > phi_min8 + epsilon) then + + !WHL - debug + if (verbose_depression .and. count >= 20) then + call parallel_globalindex(i, j, iglobal, jglobal, parallel) + print*, ' ' + print*, 'rank, i, j, ig, jg:', this_rank, i, j, iglobal, jglobal + print*, ' phi_in, phi:', phi_in(i,j), phi(i,j) + print*, ' phi_min8 =', phi_min8 + print*, ' new phi = phi_min8' + endif + + phi(i,j) = phi_min8 + epsilon + local_lowered = local_lowered + 1 + + endif ! phi_in > phi_min8 + eps, phi > phi_min8 + eps + + end if ! phi_mask = 1 and .not.known + end do ! i + end do ! j + + if (verbose_depression .and. this_rank == rtest) then + print*, ' ' + print*, 'New phi:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f11.4)',advance='no') phi(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! If one or more cells was lowered, then repeat; else exit the local loop. + + global_lowered = parallel_reduce_sum(local_lowered) + + if (global_lowered == 0) then + finished = .true. + if (verbose_depression .and. this_rank == rtest) then + print*, 'finished lowering' + endif + else + finished = .false. + if (verbose_depression .and. this_rank == rtest) then + print*, 'cells lowered on this iteration:', global_lowered + endif + call parallel_halo(phi, parallel) + endif + + if (count > count_max) then + call write_log('Hydrology error, exceeded max number of global iterations', GM_FATAL) + endif + + enddo ! finished + + if (verbose_bwat .and. this_rank == rtest) then + print*, 'Filled depressions, count =', count + endif + + end subroutine fill_depressions + +!============================================================== + + subroutine fix_flats(& + nx, ny, & + parallel, & + itest, jtest, rtest, & + phi, & + phi_mask) + + ! Add a small increment to flat regions in the input field phi, + ! so that all cells have a downhill outlet. + ! + ! Use the algorithm of Garbrecht & Mertz: + ! Garbrecht, J., and L. W. Mertz (1997), The assignment of drainage direction + ! over flat surfaces in raster digital elevation models, J. Hydrol., 193, + ! 204-213. + ! + ! Note: This subroutine is not currently called, because the depression-filling algorithm + ! above is configured not to leave any flats. + ! I am leaving it here in case it is useful for debugging. + + use cism_parallel, only: parallel_global_sum + + implicit none + + ! Input/output variables + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + real(dp), dimension(nx,ny), intent(inout) :: & + phi ! input field with flat regions to be fixed + + integer, dimension(nx,ny), intent(in) :: & + phi_mask ! = 1 where any flat regions of phi will need to be fixed, else = 0 + ! corresponds to the grounded ice sheet (bmlt_mask) for the flux-routing problem + + ! Local variables -------------------------------------- + + real(dp), dimension(nx,ny) :: & + phi_input, & ! input value of phi, before any corrections + phi_new, & ! new value of phi, after incremental corrections + dphi1, & ! sum of increments applied in step 1 + dphi2 ! sum of increments applied in step 2 + + integer, dimension(nx,ny) :: & + flat_mask, & ! = 1 for cells with phi_mask = 1 and without a downslope gradient, else = 0 + flat_mask_input, & ! flat_mask as computed from phi_input + n_uphill, & ! number of uphill neighbors for each cell, as computed from input phi + n_downhill, & ! number of downhill neighbors for each cell, as computed from input phi + incremented_mask, & ! = 1 for cells that have already been incremented (in step 2) + unincremented_mask, & ! = 1 for cells in input flat regions, not yet incremented + incremented_neighbor_mask ! = 1 for cells that have not been incremented, but have an incremented neighbor + + integer :: & + global_sum ! global sum of cells meeting a mask criterion + + logical :: finished ! true when an iterative loop has finished + + real(dp), parameter :: & + phi_increment = 2.0d-5 ! fractional increment in phi (Garbrecht & Martz use 2.0e-5) + + integer :: i, j, ii, jj, ip, jp, p + integer :: count + integer, parameter :: count_max = 100 + + ! Uncomment if the input fields are not up to date in halos +! call parallel_halo(phi, parallel) +! call parallel_halo(phi_mask, parallel) + + p = pdiag + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'In fix_flats, rtest, itest, jtest =', rtest, itest, jtest + print*, ' ' + print*, 'input phi:' + write(6,'(a3)',advance='no') ' ' + do i = itest-p, itest+p + write(6,'(i10)',advance='no') i + enddo + write(6,*) ' ' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') phi(i,j) + enddo + write(6,*) ' ' + enddo + write(6,*) ' ' + print*, 'phi_mask:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') phi_mask(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! initialize + + phi_input = phi + + ! For use in Step 2, count the uphill and downhill neighbors of each cell. + + n_uphill = 0 + n_downhill = 0 + + do j = 2, ny-1 + do i = 2, nx-1 + if (phi_mask(i,j) == 1) then + do jj = -1,1 + do ii = -1,1 + ! If this is the centre point, ignore + if (ii == 0 .and. jj == 0) then + continue + else ! check for nonzero elevation gradients + ip = i + ii + jp = j + jj + if (phi(ip,jp) - phi(i,j) > eps11) then ! uphill neighbor + n_uphill(i,j) = n_uphill(i,j) + 1 + elseif (phi(i,j) - phi(ip,jp) > eps11) then ! downhill neighbor + n_downhill(i,j) = n_downhill(i,j) + 1 + endif + endif + enddo ! ii + enddo ! jj + endif ! phi_mask = 1 + enddo ! i + enddo ! j + + ! Identify the flat regions in the input field. + ! This includes all cells with phi_mask = 1 and without downslope neighbors. + ! The resulting flat_mask is valid in all cells except the outer halo. + + call find_flats(& + nx, ny, & + itest, jtest, rtest, & + phi_input, & + phi_mask, & + flat_mask_input) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'n_uphill:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') n_uphill(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'n_downhill:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') n_downhill(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'input flat_mask:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') flat_mask_input(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! Step 1: Gradient toward lower terrain + + dphi1 = 0.0d0 + flat_mask = flat_mask_input + finished = .false. + count = 0 + + ! Increment phi in all cells with flat_mask = 1 (no downslope gradient). + ! Repeat until all cells on the global grid have a downslope gradient. + + do while(.not.finished) + + count = count + 1 + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'step 1, count =', count + endif + + where (flat_mask == 1) + dphi1 = dphi1 + phi_increment + endwhere + + call parallel_halo(dphi1, parallel) + + phi_new = phi_input + dphi1 + + ! From the original flat region, identify cells that still have no downslope gradient. + ! The resulting flat_mask is valid in all cells except the outer halo. + + call find_flats(& + nx, ny, & + itest, jtest, rtest, & + phi_new, & + flat_mask_input, & + flat_mask) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'Updated dphi1/phi_increment:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.1)',advance='no') dphi1(i,j)/ phi_increment + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'Updated flat_mask:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') flat_mask(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! Compute the number of cells in the remaining flat regions on the global grid. + ! If there are no such cells, then exit the loop. + + global_sum = parallel_global_sum(flat_mask, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, 'global sum of flat_mask =', global_sum + endif + + if (global_sum > 0) then + finished = .false. + else + finished = .true. + endif + + if (count > count_max) then + call write_log('Hydrology error: abort in step 1 of fix_flats', GM_FATAL) + endif + + enddo ! step 1 finished + + ! Step 2: Gradient away from higher terrain + + dphi2 = 0.0d0 + incremented_mask = 0 + finished = .false. + count = 0 + + ! In the first pass, increment the elevation in all cells of the input flat region that are + ! adjacent to higher terrain and have no adjacent downhill cell. + ! The resulting dphi2 and incremented_mask are valid in all cells except the outer halo. + ! Need a halo update for incremented_mask to compute incremented_neighbor_mask below. + + do j = 2, ny-1 + do i = 2, nx-1 + if (flat_mask_input(i,j) == 1) then + if (n_uphill(i,j) >= 1 .and. n_downhill(i,j) == 0) then + dphi2(i,j) = dphi2(i,j) + phi_increment + incremented_mask(i,j) = 1 + endif + endif + enddo + enddo + + call parallel_halo(dphi2, parallel) + call parallel_halo(incremented_mask, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'step 2, input flat_mask:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') flat_mask_input(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'Updated dphi2/phi_increment' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.1)',advance='no') dphi2(i,j)/phi_increment + enddo + write(6,*) ' ' + enddo + endif + + ! Compute the number of cells incremented in the first pass. + ! If no cells are incremented, then skip step 2. + ! This will be the case if the flat region lies at the highest elevation in the domain. + + global_sum = parallel_global_sum(incremented_mask, parallel) + + if (global_sum == 0) then + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'No cells to increment; skip step 2' + endif + finished = .true. + endif + + ! In subsequent passes, increment the elevation in the following cells: + ! (1) all cells that have been previously incremented; and + ! (2) all cells in the input flat region that have not been previously incremented, + ! are adjacent to an incremented cell, and are not adjacent to a cell downhill + ! from the input flat region. + ! Repeat until no unincremented cells remain on the input flat region. + ! Note: This iterated loop uses flat_mask_input, which is not incremented. + + do while(.not.finished) + + count = count + 1 + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'step 2, count =', count + endif + + ! Identify cells that have not been incremented, but are adjacent to incremented cells + ! The resulting incremented_neighbor mask is valid in all cells except the outer halo. + + incremented_neighbor_mask = 0 + do j = 2, ny-1 + do i = 2, nx-1 + if (flat_mask_input(i,j) == 1 .and. incremented_mask(i,j) == 0) then + do jj = -1,1 + do ii = -1,1 + ! If this is the centre point, ignore + if (ii == 0 .and. jj == 0) then + continue + else ! check for an incremented neighbor + ip = i + ii + jp = j + jj + if (incremented_mask(ip,jp) == 1) then + incremented_neighbor_mask(i,j) = 1 + endif + endif + enddo ! ii + enddo ! jj + endif ! flat_mask = 1 and incremented = F + enddo ! i + enddo ! j + + ! Increment cells of type (1) and (2) + ! Note: n_downhill was computed before step 1. + + do j = 2, ny-1 + do i = 2, nx-1 + if (incremented_mask(i,j) == 1) then + dphi2(i,j) = dphi2(i,j) + phi_increment + elseif (n_downhill(i,j) == 0 .and. incremented_neighbor_mask(i,j) == 1) then + dphi2(i,j) = dphi2(i,j) + phi_increment + incremented_mask(i,j) = 1 + endif + enddo + enddo + + call parallel_halo(dphi2, parallel) + call parallel_halo(incremented_mask, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'incremented_neighbor_mask:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(i10)',advance='no') incremented_neighbor_mask(i,j) + enddo + write(6,*) ' ' + enddo + print*, 'Updated dphi2/phi_increment' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.1)',advance='no') dphi2(i,j)/phi_increment + enddo + write(6,*) ' ' + enddo + endif + + ! Compute the number of cells in the input flat region that have not been incremented. + ! If all the flat cells have been incremented, then exit the loop. + + where (flat_mask_input == 1 .and. incremented_mask == 0) + unincremented_mask = 1 + elsewhere + unincremented_mask = 0 + endwhere + global_sum = parallel_global_sum(unincremented_mask, parallel) + + + if (global_sum > 0) then + if (verbose_bwat .and. this_rank == rtest) then + print*, 'number of flat cells not yet incremented =', global_sum + endif + finished = .false. + else + finished = .true. + endif + + if (count > count_max) then + call write_log('Hydrology error: abort in step 2 of fix_flats', GM_FATAL) + endif + + enddo ! step 2 finished + + + ! Step 3: + + ! Add the increments from steps 1 and 2 + ! The result is a surface with gradients both toward lower terrain and away from higher terrain. + ! No halo update is needed here, since dphi1 and dphi2 have been updated in halos. + + phi = phi_input + dphi1 + dphi2 + + ! Check for cells with flat_mask = 1 (no downslope gradient). + ! Such cells are possible because of cancelling dphi1 and dphi2. + + count = 0 + finished = .false. + + do while (.not.finished) + + count = count + 1 + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'step 3, count =', count + endif + + ! Identify cells without downslope neighbors. + ! The resulting flat_mask is valid in all cells except the outer halo. + + call find_flats(& + nx, ny, & + itest, jtest, rtest, & + phi, & + phi_mask, & + flat_mask) + + ! Add a half increment to any cells without downslope neighbors + where (flat_mask == 1) + phi = phi + 0.5d0 * phi_increment + endwhere + + call parallel_halo(phi, parallel) + + ! Compute the number of cells without downslope neighbors. + ! If there are no such cells, then exit the loop. + + global_sum = parallel_global_sum(flat_mask, parallel) + + if (verbose_bwat .and. this_rank == rtest) then + print*, 'global sum of flat_mask =', global_sum + endif + + if (global_sum > 0) then + finished = .false. + else + finished = .true. + endif + + if (count > count_max) then + call write_log('Hydrology error: abort in step 3 of fix_flats', GM_FATAL) + endif + + enddo ! step 3 finished + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'Final phi:' + do j = jtest+p, jtest-p, -1 + write(6,'(i6)',advance='no') j + do i = itest-p, itest+p + write(6,'(f10.3)',advance='no') phi(i,j) + enddo + write(6,*) ' ' + enddo + endif + + end subroutine fix_flats + +!============================================================== + + subroutine find_flats(& + nx, ny, & + itest, jtest, rtest, & + phi, phi_mask, & + flat_mask) + + ! Compute a mask that = 1 for cells in flat regions. + ! These are defined as cells with phi_mask = 1 and without a downslope gradient. + ! Note: This subroutine is not currently called. See the comment in fix_flats. + + ! Input/output arguments + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), dimension(nx,ny), intent(inout) :: & + phi ! elevation field with potential flat regions + + integer, dimension(nx,ny), intent(in) :: & + phi_mask ! = 1 for cells in the region where flats need to be identified + + integer, dimension(nx,ny), intent(out) :: & + flat_mask ! = 1 for cells with phi_mask = 1 and without a downslope gradient + + ! Local variables + + integer :: i, j, ii, jj, ip, jp + + where (phi_mask == 1) + flat_mask = 1 ! assume flat until shown otherwise + elsewhere + flat_mask = 0 + endwhere + + ! Identify cells that have no downslope neighbors, and mark them as flat. + + do j = 2, ny-1 + do i = 2, nx-1 + if (phi_mask(i,j) == 1) then + !TODO - Add an exit statement? + do jj = -1,1 + do ii = -1,1 + ! If this is the centre point, ignore + if (ii == 0 .and. jj == 0) then + continue + else ! check for a downslope gradient + ip = i + ii + jp = j + jj + if (phi(i,j) - phi(ip,jp) > eps11) then + flat_mask(i,j) = 0 + endif + endif + enddo ! ii + enddo ! jj + endif ! phi_mask = 1 + enddo ! i + enddo ! j + +! call parallel_halo(flat_mask, parallel) + + end subroutine find_flats + +!============================================================== + + subroutine sort_heights(& + nx, ny, nlocal, & + itest, jtest, rtest, & + phi, sorted_ij) + + ! Create an array with the x and y location of each cell, sorted from from low to high values of head. + ! Note: This subroutine sorts locally owned cells and excludes halo cells. + + ! Input/output arguments + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + nlocal, & ! number of locally owned cells + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), dimension(nx,ny), intent(in) :: & + phi ! input field, to be sorted from low to high + + integer, dimension(nlocal,2), intent(inout) :: & + sorted_ij ! i and j indices of each cell, sorted from from low phi to high phi + + ! Local variables + + integer :: i, j, k + integer :: ilo, ihi, jlo, jhi + integer :: nx_local, ny_local + + real(dp), dimension(nlocal) :: vect + integer, dimension(nlocal) :: ind + + ! Set array bounds for locally owned cells + ilo = nhalo+1 + ihi = nx - nhalo + jlo = nhalo+1 + jhi = ny - nhalo + nx_local = ihi-ilo+1 + ny_local = jhi-jlo+1 + + ! Fill a work vector with head values of locally owned cells + k = 1 + do i = ilo, ihi + do j = jlo, jhi + vect(k) = phi(i,j) + k = k + 1 + enddo + enddo + + ! Sort the vector from low to high values + ! The resulting 'ind' vector contains the k index for each cell, arranged from lowest to highest. + ! E.g., if the lowest-ranking cell has k = 5 and the highest-ranking cell has k = 50, + ! then ind(1) = 5 and ind(nlocal) = 50. + ! Note: For a large problem with a small number of processors, the code can fail here + ! because of too much recursion. + + call indexx(vect, ind) + + if (verbose_bwat .and. this_rank == rtest) then + print*, ' ' + print*, 'Sort from low to high, nlocal =', nlocal + print*, 'k, local i and j, ind(k), phi:' + do k = nlocal, nlocal-10, -1 + i = floor(real(ind(k)-1)/real(ny_local)) + 1 + nhalo + j = mod(ind(k)-1,ny_local) + 1 + nhalo + print*, k, i, j, ind(k), phi(i,j) + enddo + endif + + ! Fill the sorted_ij array with the i and j values of each cell. + ! Note: These are the i and j values we would have if there were no halo cells. + do k = 1, nlocal + sorted_ij(k,1) = floor(real(ind(k)-1)/real(ny_local)) + 1 + sorted_ij(k,2) = mod(ind(k)-1,ny_local) + 1 + enddo + + ! Correct the i and j values in the sorted array for halo offsets + sorted_ij(:,:) = sorted_ij(:,:) + nhalo + + end subroutine sort_heights + +!============================================================== + + subroutine get_flux_fraction(& + nx, ny, nlocal, & + dx, dy, & + itest, jtest, rtest, & + flux_routing_scheme, & + sorted_ij, & + head, & + bwat_mask, & + flux_fraction) + + ! For each cell, compute the flux fraction sent to each of the 8 neighbors, + ! based on the chosen flux routing scheme (D8, Dinf or FD8). + + ! Input/output arguments + + integer, intent(in) :: & + nx, ny, & ! number of grid cells in each direction + nlocal, & ! number of locally owned cells + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid spacing in each direction (m) + + integer, intent(in) :: & + flux_routing_scheme ! flux routing scheme: D8, Dinf or FD8 + ! D8: Flow is downhill toward the single cell with the lowest elevation. + ! Dinf: Flow is downhill toward the two cells with the lowest elevations. + ! FD8: Flow is downhill toward all cells with lower elevation. + ! D8 scheme gives the narrowest flow, and FD8 gives the most diffuse flow. + + integer, dimension(nlocal,2), intent(in) :: & + sorted_ij ! i and j indices of each cell, sorted from from low phi to high phi + + real(dp), dimension(nx,ny), intent(in) :: & + head ! hydraulic head (m) + + integer, dimension(nx,ny), intent(in) :: & + bwat_mask ! = 1 for cells in the region where basal water fluxes can be nonzero + + real(dp), dimension(-1:1,-1:1,nx,ny), intent(out) :: & + flux_fraction ! fraction of flux from a cell that flows downhill to each of 8 neighbors + + ! Local variables + + integer :: i, j, k, ii, jj, ip, jp, i1, i2, j1, j2, itmp, jtmp + + real(dp), dimension(-1:1,-1:1) :: & + dists, & ! distance (m) to adjacent grid cell + slope ! slope of head between adjacent grid cells, positive downward + + real(dp) :: & + slope1, & ! largest value of slope array + slope2, & ! second largest value of slope array + sum_slope, & ! sum of positive downward slopes + slope_tmp ! temporary slope value + + !WHL - debug + real(dp) :: sum_frac + + ! Compute distances to adjacent grid cells for slope determination + + dists(-1,:) = (/ sqrt(dx**2 + dy**2), dy, sqrt(dx**2 + dy**2) /) + dists(0,:) = (/ dx, 0.0d0, dx /) + dists(1,:) = dists(-1,:) + + ! Loop through locally owned cells and compute the flux fraction sent to each neighbor cell. + ! This fraction is stored in an array of dimension (-1:1,-1:1,nx,ny). + ! The (0,0) element refers to the cell itself and is equal to 0 for each i and j. + + flux_fraction = 0.0d0 + + do k = nlocal, 1, -1 + + ! Get i and j indices of current point + i = sorted_ij(k,1) + j = sorted_ij(k,2) + + if (bwat_mask(i,j) == 1) then + + ! Compute the slope between this cell and each neighbor. + ! Slopes are defined as positive for downhill neighbors, and zero otherwise. + + slope = 0.0d0 + + ! Loop over adjacent points and calculate slope + do jj = -1,1 + do ii = -1,1 + ! If this is the centre point, ignore + if (ii == 0 .and. jj == 0) then + continue + else ! compute slope + ip = i + ii + jp = j + jj + if (ip >= 1 .and. ip <= nx .and. jp > 1 .and. jp <= ny) then + if (head(ip,jp) < head(i,j)) then + slope(ii,jj) = (head(i,j) - head(ip,jp)) / dists(ii,jj) + endif + endif + endif + enddo + enddo + + sum_slope = sum(slope) + + !WHL - debug + if (this_rank == rtest .and. i == itest .and. j == jtest) then + print*, ' ' + print*, 'slope: task, i, j =', rtest, i, j + print*, slope(:,1) + print*, slope(:,0) + print*, slope(:,-1) + print*, 'sum(slope) =', sum(slope) + endif + + ! Distribute the downslope flux according to the flux-routing scheme: + ! to the lowest-elevation neighbor (D8), the two lowest-elevation neighbors (Dinf), or + ! all lower-elevation neighbors (FD8). + ! The D8 and FD8 schemes have been tested with a simple dome problem. + ! Dinf is less suited for the dome problem because there are many ties for 2nd greatest slope, + ! so i2 and j2 for slope2 are not well defined. + + if (flux_routing_scheme == HO_FLUX_ROUTING_D8) then + + ! route to the adjacent cell with the lowest elevation + slope1 = 0.0d0 + if (sum_slope > 0.d0) then + i1 = 0; j1 = 0 + do jj = -1,1 + do ii = -1,1 + ip = i + ii + jp = j + jj + if (slope(ii,jj) > slope1) then + slope1 = slope(ii,jj) + i1 = ip + j1 = jp + endif + enddo + enddo + endif + + if (slope1 > 0.0d0) then + ii = i1 - i + jj = j1 - j + flux_fraction(ii,jj,i,j) = 1.0d0 ! route the entire flux to one downhill cell + else + ! Do a fatal abort? + print*, 'Warning: Cell with no downhill neighbors, i, j =', i, j + endif + + if (this_rank == rtest .and. i == itest .and. j == jtest) then + print*, 'i1, j1, slope1 =', i1, j1, slope1 + endif + + elseif (flux_routing_scheme == HO_FLUX_ROUTING_DINF) then + + ! route to the two adjacent cells with the lowest elevation + i1 = 0; j1 = 0 + i2 = 0; j2 = 0 + slope1 = 0.0d0 + slope2 = 0.0d0 + do jj = -1,1 + do ii = -1,1 + ip = i + ii + jp = j + jj + if (slope(ii,jj) > slope1) then + slope_tmp = slope1 + itmp = i1 + jtmp = j1 + slope1 = slope(ii,jj) + i1 = ip + j1 = jp + slope2 = slope_tmp + i2 = itmp + j2 = jtmp + elseif (slope(ii,jj) > slope2) then + slope2 = slope(ii,jj) + i2 = ip + j2 = jp + endif + enddo + enddo + + sum_slope = slope1 + slope2 ! divide the flux between cells (i1,j1) and (i2,j2) + if (sum_slope > 0.0d0) then + if (slope1 > 0.0d0) then + ii = i1 - i + jj = j1 - j + flux_fraction(ii,jj,i,j) = slope1/sum_slope + endif + if (slope2 > 0.0d0) then + ii = i2 - i + jj = j2 - j + flux_fraction(ii,jj,i,j) = slope2/sum_slope + endif + else + print*, 'Warning: Cell with no downhill neighbors, i, j =', i, j + endif + + if (this_rank == rtest .and. i == itest .and. j == jtest) then + print*, 'i1, j1, slope1:', i1, j1, slope1 + print*, 'i2, j2, slope2:', i2, j2, slope2 + print*, 'sum_slope:', sum_slope + print*, 'slope(:, 1):', slope(:, 1) + print*, 'slope(:, 0):', slope(:, 0) + print*, 'slope(:,-1):', slope(:,-1) + print*, 'flux_fraction(:, 1,i,j):', flux_fraction(:, 1,i,j) + print*, 'flux_fraction(:, 0,i,j):', flux_fraction(:, 0,i,j) + print*, 'flux_fraction(:,-1,i,j):', flux_fraction(:,-1,i,j) + endif + + !WHL - bug check - make sure fractions add to 1 + sum_frac = 0.0d0 + do jj = -1,1 + do ii = -1,1 + sum_frac = sum_frac + flux_fraction(ii,jj,i,j) + enddo + enddo + if (abs(sum_frac - 1.0d0) > eps11) then +!! print*, 'sum_frac error: r, i, j, sum:', this_rank, i, j, sum_frac + endif + + elseif (flux_routing_scheme == HO_FLUX_ROUTING_FD8) then + + ! route to all adjacent downhill cells in proportion to grad(head) + if (sum_slope > 0.d0) then + do jj = -1,1 + do ii = -1,1 + ip = i + ii + jp = j + jj + if (slope(ii,jj) > 0.d0) then + flux_fraction(ii,jj,i,j) = slope(ii,jj)/sum_slope + endif + enddo + enddo + endif ! sum(slope) > 0 + + endif ! flux_routing_scheme: D8, Dinf, FD8 + + endif ! bwat_mask = 1 + + enddo ! loop from high to low + + end subroutine get_flux_fraction + +!============================================================== + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ! + ! The following two subroutines perform an index-sort of an array. + ! They are a GPL-licenced replacement for the Numerical Recipes routine indexx. + ! They are not derived from any NR code, but are based on a quicksort routine by + ! Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + ! in C, and issued under the GNU General Public License. The conversion to + ! Fortran 90, and modification to do an index sort was done by Ian Rutt. + ! + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine indexx(array, index) + + !> Performs an index sort of \texttt{array} and returns the result in + !> \texttt{index}. The order of elements in \texttt{array} is unchanged. + !> + !> This is a GPL-licenced replacement for the Numerical Recipes routine indexx. + !> It is not derived from any NR code, but are based on a quicksort routine by + !> Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + !> in C, and issued under the GNU General Public License. The conversion to + !> Fortran 90, and modification to do an index sort was done by Ian Rutt. + + real(dp),dimension(:) :: array !> Array to be indexed. + integer, dimension(:) :: index !> Index of elements of \texttt{array}. + integer :: i + + if (size(array) /= size(index)) then + call write_log('ERROR: INDEXX size mismatch.',GM_FATAL,__FILE__,__LINE__) + endif + + do i = 1,size(index) + index(i) = i + enddo + + call q_sort_index(array, index, 1, size(array)) + + end subroutine indexx + +!============================================================== + + recursive subroutine q_sort_index(numbers, index, left, right) + + !> This is the recursive subroutine actually used by \texttt{indexx}. + !> + !> This is a GPL-licenced replacement for the Numerical Recipes routine indexx. + !> It is not derived from any NR code, but are based on a quicksort routine by + !> Michael Lamont (http://linux.wku.edu/~lamonml/kb.html), originally written + !> in C, and issued under the GNU General Public License. The conversion to + !> Fortran 90, and modification to do an index sort was done by Ian Rutt. + !> + !> Note: For a large problem with a small number of processors, the code can + ! fail here with a seg fault because there is too much recursion. + + implicit none + + real(dp),dimension(:) :: numbers !> Numbers being sorted + integer, dimension(:) :: index !> Returned index + integer :: left, right !> Limit of sort region + + integer :: ll, rr + integer :: pv_int, l_hold, r_hold, pivpos + real(dp) :: pivot + + ll = left + rr = right + + l_hold = ll + r_hold = rr + pivot = numbers(index(ll)) + pivpos = index(ll) + + do + if (.not.(ll < rr)) exit + + do + if (.not.((numbers(index(rr)) >= pivot) .and. (ll < rr))) exit + rr = rr - 1 + enddo + + if (ll /= rr) then + index(ll) = index(rr) + ll = ll + 1 + endif + + do + if (.not.((numbers(index(ll)) <= pivot) .and. (ll < rr))) exit + ll = ll + 1 + enddo + + if (ll /= rr) then + index(rr) = index(ll) + rr = rr - 1 + endif + enddo + + index(ll) = pivpos + pv_int = ll + ll = l_hold + rr = r_hold + if (ll < pv_int) call q_sort_index(numbers, index,ll, pv_int-1) + if (rr > pv_int) call q_sort_index(numbers, index,pv_int+1, rr) + + end subroutine q_sort_index + +!============================================================== + end module glissade_basal_water + +!============================================================== diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 96925892..12fadf6a 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -1392,6 +1392,7 @@ subroutine glissade_remove_icebergs(& real(dp), dimension(nx,ny) :: & thck_calving_front ! effective ice thickness at the calving front + !TODO - Make this a config parameter? real(dp), parameter :: & ! threshold for counting cells as grounded f_ground_threshold = 0.10d0 @@ -1632,7 +1633,7 @@ subroutine glissade_remove_isthmuses(& ocean_plus_thin_ice_mask ! = 1 for ocean cells and cells with thin floating ice ! Both floating and weakly grounded cells can be identified as isthmuses and removed; - ! isthmuses_f_ground_threshold is used to identify weakly grounded cells. + ! isthmus_f_ground_threshold is used to identify weakly grounded cells. real(dp), parameter :: & ! threshold for counting cells as grounded isthmus_f_ground_threshold = 0.50d0 diff --git a/libglissade/glissade_grid_operators.F90 b/libglissade/glissade_grid_operators.F90 index 09757411..e6354769 100644 --- a/libglissade/glissade_grid_operators.F90 +++ b/libglissade/glissade_grid_operators.F90 @@ -578,6 +578,9 @@ subroutine glissade_gradient_at_edges(nx, ny, & ! !-------------------------------------------------------- + ! TODO - Make HO_GRADIENT_MARGIN_LAND the default, since it is simple and requires no optional arguments? + ! TODO - Make ice_mask an optional argument, = 1 everywhere by default. + if (present(gradient_margin_in)) then gradient_margin = gradient_margin_in else @@ -585,7 +588,6 @@ subroutine glissade_gradient_at_edges(nx, ny, & endif ! Set logical edge mask based on gradient_margin. - edge_mask_x(:,:) = .false. edge_mask_y(:,:) = .false. diff --git a/libglissade/glissade_therm.F90 b/libglissade/glissade_therm.F90 index 10b1fcca..f155521b 100644 --- a/libglissade/glissade_therm.F90 +++ b/libglissade/glissade_therm.F90 @@ -1115,11 +1115,14 @@ subroutine glissade_therm_driver(whichtemp, & dissipcol(ew,ns) = dissipcol(ew,ns) * thck(ew,ns)*rhoi*shci ! Verify that the net input of energy into the column is equal to the change in - ! internal energy. + ! internal energy. delta_e = (ucondflx(ew,ns) - lcondflx(ew,ns) + dissipcol(ew,ns)) * dttem - if (abs((efinal-einit-delta_e)/dttem) > 1.0d-8) then + ! Note: For very small dttem (e.g., 1.0d-6 year or less), this error can be triggered + ! by roundoff error. In that case, the user may need to increase the threshold. + ! July 2021: Increased from 1.0d-8 to 1.0d-7 to allow smaller dttem. + if (abs((efinal-einit-delta_e)/dttem) > 1.0d-7) then if (verbose_column) then print*, 'Ice thickness:', thck(ew,ns) @@ -1640,10 +1643,11 @@ subroutine glissade_enthalpy_matrix_elements(dttem, & ! At each temperature point, compute the temperature part of the enthalpy. ! enth_T = enth for cold ice, enth_T < enth for temperate ice - enth_T(0) = rhoi*shci*temp(0) !WHL - not sure enth_T(0) is needed - do up = 1, upn + do up = 1, upn-1 enth_T(up) = (1.d0 - waterfrac(up)) * rhoi*shci*temp(up) enddo + enth_T(0) = rhoi*shci*temp(0) + enth_T(up) = rhoi*shci*temp(up) !WHL - debug if (verbose_column) then @@ -2250,8 +2254,7 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, & ! Compute the dissipation source term associated with strain heating, ! based on the shallow-ice approximation. - use glimmer_physcon, only : gn ! Glen's n - use glimmer_physcon, only: rhoi, shci, grav + use glimmer_physcon, only: rhoi, shci, grav, n_glen integer, intent(in) :: ewn, nsn, upn ! grid dimensions @@ -2267,12 +2270,14 @@ subroutine glissade_interior_dissipation_sia(ewn, nsn, & real(dp), dimension(:,:,:), intent(out) :: & dissip ! interior heat dissipation (deg/s) - integer, parameter :: p1 = gn + 1 - integer :: ew, ns real(dp), dimension(upn-1) :: sia_dissip_fact ! factor in SIA dissipation calculation real(dp) :: geom_fact ! geometric factor + real(dp) :: p1 ! exponent = n_glen + 1 + + p1 = n_glen + 1.0d0 + ! Two methods of doing this calculation: ! 1. find dissipation at u-pts and then average ! 2. find dissipation at H-pts by averaging quantities from u-pts @@ -2414,7 +2419,7 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & integer, intent(in) :: whichflwa !> which method of calculating A integer, intent(in) :: whichtemp !> which method of calculating temperature; - !> include waterfrac in calculation if using enthalpy method + !> include waterfrac in calculation if using enthalpy method real(dp),dimension(:), intent(in) :: stagsigma !> vertical coordinate at layer midpoints real(dp),dimension(:,:), intent(in) :: thck !> ice thickness (m) real(dp),dimension(:,:,:), intent(in) :: temp !> 3D temperature field (deg C) @@ -2488,17 +2493,16 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & endif ! Multiply the default rate factor by the enhancement factor if applicable - ! Note: Here, default_flwa is assumed to have units of Pa^{-n} s^{-1}, + ! Note: Here, the input default_flwa is assumed to have units of Pa^{-n} s^{-1}, ! whereas model%paramets%default_flwa has units of Pa^{-n} yr^{-1}. ! initialize - if (whichflwa /= FLWA_INPUT) then - do ns = 1, nsn - do ew = 1, ewn - flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa - enddo + !TODO - Move the next few lines inside the select case construct. + do ns = 1, nsn + do ew = 1, ewn + flwa(:,ew,ns) = enhancement_factor(ew,ns) * default_flwa enddo - endif + enddo select case(whichflwa) @@ -2558,21 +2562,12 @@ subroutine glissade_flow_factor(whichflwa, whichtemp, & end do case(FLWA_CONST_FLWA) - - ! do nothing (flwa is initialized to default_flwa above) - - case(FLWA_INPUT) - ! do nothing - use flwa from input or forcing file - print *, 'FLWA', minval(flwa), maxval(flwa) + ! do nothing (flwa is set above, with units Pa^{-n} s^{-1}) end select - ! This logic assumes that the input flwa is already in dimensionless model units. - ! TODO: Make a different assumption about input units? - if (whichflwa /= FLWA_INPUT) then - ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1}) - flwa(:,:,:) = flwa(:,:,:) / vis0 - endif + ! Change flwa to model units (glissade_flow_factor assumes SI units of Pa{-n} s^{-1}) + flwa(:,:,:) = flwa(:,:,:) / vis0 deallocate(enhancement_factor) diff --git a/libglissade/glissade_transport.F90 b/libglissade/glissade_transport.F90 index a1c57219..583ccb84 100644 --- a/libglissade/glissade_transport.F90 +++ b/libglissade/glissade_transport.F90 @@ -979,6 +979,7 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & parallel, & stagthk, dusrfdew, dusrfdns, & uvel, vvel, deltat, & + adaptive_cfl_threshold, & allowable_dt_adv, allowable_dt_diff) ! Calculate maximum allowable time step based on both @@ -1015,6 +1016,10 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & real(dp), intent(in) :: & deltat ! model deltat (yrs) + real(dp), intent(in) :: & + adaptive_cfl_threshold ! threshold for adaptive subcycling + ! if = 0, there is no adaptive subcycling; code aborts when CFL > 1 + real(dp), intent(out) :: & allowable_dt_adv ! maximum allowable dt (yrs) based on advective CFL @@ -1087,7 +1092,8 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & maxvel = maxvvel indices_adv = maxloc(abs(vvel_layer(:,xs:xe,ys:ye))) endif - indices_adv(2:3) = indices_adv(2:3) + staggered_lhalo ! want the i,j coordinates WITH the halo present - we got indices into the slice of owned cells + indices_adv(2:3) = indices_adv(2:3) + staggered_lhalo ! want the i,j coordinates WITH the halo present - + ! we got indices into the slice of owned cells ! Finally, determine maximum allowable time step based on advectice CFL condition. my_allowable_dt_adv = dew / (maxvel + 1.0d-20) @@ -1159,13 +1165,16 @@ subroutine glissade_check_cfl(ewn, nsn, nlyr, & write(message,*) 'Advective CFL violation! Maximum allowable time step for advective CFL condition is ' & // trim(adjustl(dt_string)) // ' yr, limited by global position i=' & // trim(adjustl(xpos_string)) // ' j=' //trim(adjustl(ypos_string)) + call write_log(trim(message),GM_WARNING) + + ! If adaptive subcyling is allowed, then make the code abort for egregious CFL violations, + ! (defined as deltat > 10 * allowable_dt_adv), to prevent excessive subcycling. - ! If the violation is egregious (defined as deltat > 10 * allowable_dt_adv), then abort. - ! Otherwise, write a warning and proceed. - if (deltat > 10.d0 * allowable_dt_adv) then - call write_log(trim(message),GM_FATAL) - else - call write_log(trim(message),GM_WARNING) + if (main_task .and. adaptive_cfl_threshold > 0.0d0) then + if (deltat > 10.d0 * allowable_dt_adv) then + print*, 'deltat, allowable_dt_adv, ratio =', deltat, allowable_dt_adv, deltat/allowable_dt_adv + call write_log('Aborting with CFL violation', GM_FATAL) + endif endif endif @@ -1662,9 +1671,10 @@ subroutine glissade_overwrite_acab_mask(overwrite_acab, & use glide_types ! If overwrite_acab /=0 , then set overwrite_acab_mask = 1 for grid cells - ! where acab is to be overwritten. Currently, two options are supported: + ! where acab is to be overwritten. Currently, three options are supported: ! (1) Overwrite acab where the input acab = 0 at initialization ! (2) Overwrite acab where the input thck <= overwrite_acab_minthck at initialization + ! (3) Overwrite acab based on an input mask ! ! Note: This subroutine should be called only on initialization, not on restart. @@ -1683,6 +1693,7 @@ subroutine glissade_overwrite_acab_mask(overwrite_acab, & integer :: ewn, nsn integer :: i, j + integer :: max_mask_local, max_mask_global ewn = size(overwrite_acab_mask,1) nsn = size(overwrite_acab_mask,2) @@ -1716,6 +1727,25 @@ subroutine glissade_overwrite_acab_mask(overwrite_acab, & enddo enddo + elseif (overwrite_acab == OVERWRITE_ACAB_INPUT_MASK) then + + ! Make sure a mask was read in with some nonzero values + ! If not, then write a warning + + max_mask_local = maxval(overwrite_acab_mask) + max_mask_global = parallel_reduce_max(max_mask_local) + if (main_task) then + print*, 'rank, max_mask_local, max_mask_global:', & + this_rank, max_mask_local, max_mask_global + endif + if (max_mask_global == 1) then + ! continue + elseif (max_mask_global == 0) then + call write_log('Using overwrite_acab_mask without any values > 0', GM_WARNING) + else + call write_log('Using overwrite_acab_mask with values other than 0 and 1', GM_FATAL) + endif + endif ! overwrite_acab end subroutine glissade_overwrite_acab_mask diff --git a/libglissade/glissade_utils.F90 b/libglissade/glissade_utils.F90 index 0c4de494..36c2bece 100644 --- a/libglissade/glissade_utils.F90 +++ b/libglissade/glissade_utils.F90 @@ -38,7 +38,8 @@ module glissade_utils implicit none private - public :: glissade_adjust_thickness, glissade_smooth_topography, glissade_adjust_topography + public :: glissade_adjust_thickness, glissade_smooth_usrf, & + glissade_smooth_topography, glissade_adjust_topography public :: glissade_stdev, verbose_stdev logical, parameter :: verbose_stdev = .true. @@ -216,6 +217,187 @@ subroutine glissade_adjust_thickness(model) end subroutine glissade_adjust_thickness +!**************************************************************************** + + subroutine glissade_smooth_usrf(model, nsmooth) + + ! Use a Laplacian smoother to smooth the upper surface elevation, + ! and compute a thickness consistent with this new elevation. + ! This can be useful if the input thickness and topography are inconsistent, + ! such that their sum has large gradients. + + use glimmer_paramets, only: thk0 + use glide_thck, only: glide_calclsrf + use glissade_masks, only: glissade_get_masks + use glissade_grid_operators, only: glissade_laplacian_smoother + use cism_parallel, only: parallel_halo + + !---------------------------------------------------------------- + ! Input-output arguments + !---------------------------------------------------------------- + + type(glide_global_type), intent(inout) :: model ! derived type holding ice-sheet info + + integer, intent(in), optional :: nsmooth ! number of smoothing passes + + ! local variables + + real(dp), dimension(model%general%ewn, model%general%nsn) :: & + topg, & ! bed topography (m) + thck, & ! thickness (m) + usrf, & ! surface elevation (m) + usrf_smoothed ! surface elevation after smoothing + + integer, dimension(model%general%ewn, model%general%nsn) :: & + ice_mask, & ! = 1 if ice is present (thck > 0, else = 0 + floating_mask, & ! = 1 if ice is present (thck > 0) and floating, else = 0 + ocean_mask ! = 1 if topg < 0 and ice is absent, else = 0 + + integer :: n_smoothing_passes ! local version of nsmooth + integer :: i, j, n + integer :: nx, ny + integer :: itest, jtest, rtest + +! logical, parameter :: verbose_smooth_usrf = .false. + logical, parameter :: verbose_smooth_usrf = .true. + + ! Initialize + + if (present(nsmooth)) then + n_smoothing_passes = nsmooth + else + n_smoothing_passes = 1 + endif + + ! Copy some model variables to local variables + + nx = model%general%ewn + ny = model%general%nsn + + rtest = -999 + itest = 1 + jtest = 1 + if (this_rank == model%numerics%rdiag_local) then + rtest = model%numerics%rdiag_local + itest = model%numerics%idiag_local + jtest = model%numerics%jdiag_local + endif + + ! compute the initial upper surface elevation + call glide_calclsrf(model%geometry%thck, model%geometry%topg, model%climate%eus, model%geometry%lsrf) + model%geometry%usrf = max(0.d0, model%geometry%thck + model%geometry%lsrf) + + ! Save input fields + topg = (model%geometry%topg - model%climate%eus) * thk0 + thck = model%geometry%thck * thk0 + usrf = model%geometry%usrf * thk0 + + if (verbose_smooth_usrf .and. this_rank == rtest) then + i = itest + j = jtest + print*, ' ' + print*, 'itest, jtest, rank =', itest, jtest, rtest + print*, ' ' + print*, 'Before Laplacian smoother, topg (m):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') topg(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'Before Laplacian smoother, usrf (m):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') usrf(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'Before Laplacian smoother, thck (m):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') thck(i,j) + enddo + write(6,*) ' ' + enddo + endif + + ! compute initial masks + call glissade_get_masks(nx, ny, & + model%parallel, & + model%geometry%thck, model%geometry%topg, & + model%climate%eus, 0.0d0, & ! thklim = 0 + ice_mask, & + floating_mask = floating_mask, & + ocean_mask = ocean_mask) + + do n = 1, n_smoothing_passes + + call glissade_laplacian_smoother(nx, ny, & + usrf, usrf_smoothed, & + npoints_stencil = 9) + + ! Force usrf = topg on ice-free land + where (topg > 0.0d0 .and. ice_mask == 0) usrf_smoothed = topg + + ! Force usrf = unsmoothed value for floating ice and ice-free ocean, to avoid advancing the calving front + where (floating_mask == 1 .or. ocean_mask == 1) + usrf_smoothed = usrf + endwhere + + ! Force usrf >= topg + usrf_smoothed = max(usrf_smoothed, topg) + + usrf = usrf_smoothed + call parallel_halo(usrf, model%parallel) + + enddo + + ! Given the smoothed usrf, adjust the input thickness such that topg is unchanged. + ! Do this only where ice is present. Elsewhere, usrf = topg. + + where (usrf > topg) ! ice is present + where (topg < 0.0d0) ! marine-based ice + where (topg*(1.0d0 - rhoo/rhoi) > usrf) ! ice is floating + thck = usrf / (1.0d0 - rhoi/rhoo) + elsewhere ! ice is grounded + thck = usrf - topg + endwhere + elsewhere ! land-based ice + thck = usrf - topg + endwhere + endwhere + + ! Copy the new thickness and usrf to the model derived type + model%geometry%thck = thck/thk0 + model%geometry%usrf = usrf/thk0 + + if (verbose_smooth_usrf .and. this_rank == rtest) then + i = itest + j = jtest + print*, ' ' + print*, 'itest, jtest, rank =', itest, jtest, rtest + print*, ' ' + print*, 'After Laplacian smoother, usrf (m):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') usrf(i,j) + enddo + write(6,*) ' ' + enddo + print*, ' ' + print*, 'After Laplacian smoother, thck (m):' + do j = jtest+3, jtest-3, -1 + do i = itest-3, itest+3 + write(6,'(f10.3)',advance='no') thck(i,j) + enddo + write(6,*) ' ' + enddo + endif + + end subroutine glissade_smooth_usrf + !**************************************************************************** subroutine glissade_smooth_topography(model) diff --git a/libglissade/glissade_velo.F90 b/libglissade/glissade_velo.F90 index bb6e6e38..a9dadd17 100644 --- a/libglissade/glissade_velo.F90 +++ b/libglissade/glissade_velo.F90 @@ -43,9 +43,6 @@ subroutine glissade_velo_driver(model) ! Glissade higher-order velocity driver - use glimmer_global, only : dp - use glimmer_physcon, only: gn, scyr - use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0, evs0 use glimmer_log use glide_types use glissade_velo_higher, only: glissade_velo_higher_solve diff --git a/libglissade/glissade_velo_higher.F90 b/libglissade/glissade_velo_higher.F90 index 033ca254..a58f4a22 100644 --- a/libglissade/glissade_velo_higher.F90 +++ b/libglissade/glissade_velo_higher.F90 @@ -57,7 +57,7 @@ module glissade_velo_higher use glimmer_global, only: dp - use glimmer_physcon, only: gn, rhoi, rhoo, grav, scyr, pi + use glimmer_physcon, only: n_glen, rhoi, rhoo, grav, scyr, pi use glimmer_paramets, only: eps08, eps10, thk0, len0, tim0, tau0, vel0, vis0, evs0 use glimmer_paramets, only: vel_scale, len_scale ! used for whichefvs = HO_EFVS_FLOWFACT use glimmer_log @@ -200,8 +200,8 @@ module glissade_velo_higher ! logical :: verbose = .true. logical :: verbose_init = .false. ! logical :: verbose_init = .true. -! logical :: verbose_solver = .false. - logical :: verbose_solver = .true. + logical :: verbose_solver = .false. +! logical :: verbose_solver = .true. logical :: verbose_Jac = .false. ! logical :: verbose_Jac = .true. logical :: verbose_residual = .false. @@ -1108,7 +1108,7 @@ subroutine glissade_velo_higher_solve(model, & beta_internal => model%velocity%beta_internal(:,:) bfricflx => model%temper%bfricflx(:,:) bpmp => model%temper%bpmp(:,:) - bwat => model%temper%bwat(:,:) + bwat => model%basal_hydro%bwat(:,:) bmlt => model%basal_melt%bmlt(:,:) uvel => model%velocity%uvel(:,:,:) @@ -2016,12 +2016,14 @@ subroutine glissade_velo_higher_solve(model, & ! Note: Ideally, bpmp and temp(nz) are computed after the transport solve, ! just before the velocity solve. Then they will be consistent with the ! current thickness field. + ! TODO: Move this call to a higher level. Does not need any velocity information. !------------------------------------------------------------------------------ !TODO - Use btemp_ground instead of temp(nz)? call calc_effective_pressure(whicheffecpress, & nx, ny, & model%basal_physics, & + model%basal_hydro, & ice_mask, floating_mask, & thck, topg, & eus, & @@ -2082,7 +2084,7 @@ subroutine glissade_velo_higher_solve(model, & ! gn = exponent in Glen's flow law (= 3 by default) do k = 1, nz-1 if (flwa(k,i,j) > 0.0d0) then - flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/real(gn,dp)) + flwafact(k,i,j) = 0.5d0 * flwa(k,i,j)**(-1.d0/n_glen) endif enddo enddo @@ -4222,6 +4224,7 @@ subroutine glissade_velo_higher_solve(model, & usrf, & dusrf_dx, dusrf_dy, & flwa, efvs, & + whichefvs, efvs_constant, & whichgradient_margin, & max_slope, & uvel, vvel) @@ -6426,6 +6429,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & usrf, & dusrf_dx, dusrf_dy, & flwa, efvs, & + whichefvs, efvs_constant, & whichgradient_margin, & max_slope, & uvel, vvel) @@ -6486,6 +6490,12 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & flwa, & ! temperature-based flow factor A, Pa^{-n} yr^{-1} efvs ! effective viscosity, Pa yr + integer, intent(in) :: & + whichefvs ! option for effective viscosity calculation + + real(dp), intent(in) :: & + efvs_constant ! constant value of effective viscosity (Pa yr) + integer, intent(in) :: & whichgradient_margin ! option for computing gradient at ice margin ! 0 = include all neighbor cells in gradient calculation @@ -6840,7 +6850,7 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & ! Compute vertical integration factor at each active vertex ! This is int_b_to_z{-2 * A * tau^2 * rho*g*(s-z) * dz}, - ! similar to the factor computed in Glide and glissade_velo_sia.. + ! similar to the factor computed in Glide and glissade_velo_sia. ! Note: tau_xz ~ rho*g*(s-z)*ds_dx; ds_dx term is computed on edges below do j = 1, ny-1 @@ -6921,9 +6931,27 @@ subroutine compute_3d_velocity_L1L2(nx, ny, & tau_eff_sq = stagtau_parallel_sq(i,j) & + tau_xz(k,i,j)**2 + tau_yz(k,i,j)**2 - ! Note: This formula is correct for any value of Glen's n, but currently efvs is computed - ! only for gn = 3 (in which case (n-1)/2 = 1). - fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((gn-1.d0)/2.d0) * (sigma(k+1) - sigma(k))*stagthck(i,j) + ! Note: The first formula below is correct for whichefvs = 2 (efvs computed from effective strain rate), + ! but not for whichefvs = 0 (constant efvs) or whichefvs = 1 (multiple of flow factor). + ! For these options we need a modified formula. + ! + ! Recall: efvs = 1/2 * A^(-1/n) * eps_e^[(1-n)/n] + ! = 1/2 * A^(-1/n) * [A tau_e^n]^[(1-n)/n] + ! = 1/2 * A^(-1) * tau_e^(1-n) + ! => 1/efvs = 2 * A * tau_e(n-1) + ! + ! Thus, for options 0 and 1, we can replace 2 * A * tau_e^(n-1) below with 1/efvs. + + if (whichefvs == HO_EFVS_NONLINEAR) then + fact = 2.d0 * stagflwa(i,j) * tau_eff_sq**((n_glen-1.d0)/2.d0) & + * (sigma(k+1) - sigma(k))*stagthck(i,j) + else ! HO_EFVS_CONSTANT, HO_EFVS_FLOWFACT + if (efvs(k,i,j) > 0.0d0) then + fact = (sigma(k+1) - sigma(k))*stagthck(i,j) / efvs(k,i,j) + else + fact = 0.0d0 + endif + endif ! reset velocity to prescribed basal value if Dirichlet condition applies ! else compute velocity at this level @@ -7875,15 +7903,6 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & integer, intent(in) :: i, j, k, p - !---------------------------------------------------------------- - ! Local parameters - !---------------------------------------------------------------- - - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation - p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation - - !---------------------------------------------------------------- ! Local variables !---------------------------------------------------------------- @@ -7896,8 +7915,14 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & integer :: n + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation + real(dp), parameter :: p2 = -1.d0/3.d0 + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen)/n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -7988,11 +8013,11 @@ subroutine compute_effective_viscosity (whichefvs, whichapprox, & ! Compute effective viscosity (PGB 2012, eq. 4) ! Units: flwafact has units Pa yr^{1/n} ! effstrain has units yr^{-1} - ! p2_effstr = (1-n)/(2n) - ! = -1/3 for n=3 + ! p_effstr = (1-n)/n + ! = -2/3 for n=3 ! Thus efvs has units Pa yr - efvs = flwafact * effstrainsq**p2_effstr + efvs = flwafact * effstrainsq**(p_effstr/2.d0) if (verbose_efvs .and. this_rank==rtest .and. i==itest .and. j==jtest .and. k==ktest .and. p==ptest) then print*, ' ' @@ -8081,8 +8106,8 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, ! Local parameters !---------------------------------------------------------------- - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp)) / real(gn,dp) ! exponent (1-n)/n in effective viscosity relation + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation !---------------------------------------------------------------- ! Local variables @@ -8107,6 +8132,9 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, integer :: n, k + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen) / n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -8125,7 +8153,7 @@ subroutine compute_effective_viscosity_L1L2(whichefvs, ! ! Units: flwafact has units Pa yr^{1/n} ! effstrain has units yr^{-1} - ! p_effstr = (1-n)/n + ! p_effstr = (1-n)/n ! = -2/3 for n=3 ! Thus efvs has units Pa yr @@ -8320,14 +8348,6 @@ subroutine compute_effective_viscosity_diva(whichefvs, integer, intent(in) :: i, j, p - !---------------------------------------------------------------- - ! Local parameters - !---------------------------------------------------------------- - - real(dp), parameter :: & - p_effstr = (1.d0 - real(gn,dp))/real(gn,dp), &! exponent (1-n)/n in effective viscosity relation - p2_effstr = p_effstr/2 ! exponent (1-n)/(2n) in effective viscosity relation - !---------------------------------------------------------------- ! Local variables !---------------------------------------------------------------- @@ -8346,11 +8366,17 @@ subroutine compute_effective_viscosity_diva(whichefvs, integer :: n, k real(dp) :: du_dz, dv_dz + real(dp) :: & + p_effstr ! exponent (1-n)/n in effective viscosity relation + !WHL - For ISMIP-HOM, the cubic solve is not robust. It leads to oscillations ! in successive iterations between uvel_2d/vvel_2d and btractx/btracty !TODO - Remove the cubic solve for efvs, unless we find a way to make it robust? logical, parameter :: cubic = .false. + ! Set exponent that depends on Glen's exponent + p_effstr = (1.d0 - n_glen)/n_glen + select case(whichefvs) case(HO_EFVS_CONSTANT) @@ -8493,7 +8519,8 @@ subroutine compute_effective_viscosity_diva(whichefvs, effstrainsq = effstrain_min**2 & + du_dx**2 + dv_dy**2 + du_dx*dv_dy + 0.25d0*(dv_dx + du_dy)**2 & + 0.25d0 * (du_dz**2 + dv_dz**2) - efvs(k) = flwafact(k) * effstrainsq**p2_effstr + efvs(k) = flwafact(k) * effstrainsq**(p_effstr/2.d0) + enddo endif ! cubic diff --git a/libglissade/glissade_velo_higher_pcg.F90 b/libglissade/glissade_velo_higher_pcg.F90 index 7bac76c9..1c2003ea 100644 --- a/libglissade/glissade_velo_higher_pcg.F90 +++ b/libglissade/glissade_velo_higher_pcg.F90 @@ -63,14 +63,6 @@ module glissade_velo_higher_pcg module procedure global_sum_staggered_2d_real8_nvar end interface - ! linear solver settings - !TODO - Pass in these solver settings as arguments? -! integer, parameter :: & -! maxiters = 200 ! max number of linear iterations before quitting - ! TODO - change to maxiters_default? -! real(dp), parameter :: & -! tolerance = 1.d-08 ! tolerance for linear solver - logical, parameter :: verbose_pcg = .false. logical, parameter :: verbose_tridiag = .false. !! logical, parameter :: verbose_pcg = .true. diff --git a/libglissade/glissade_velo_sia.F90 b/libglissade/glissade_velo_sia.F90 index 66884aaf..3949f9f4 100644 --- a/libglissade/glissade_velo_sia.F90 +++ b/libglissade/glissade_velo_sia.F90 @@ -55,7 +55,7 @@ module glissade_velo_sia use glimmer_global, only: dp - use glimmer_physcon, only: gn, rhoi, grav, scyr + use glimmer_physcon, only: n_glen, rhoi, grav, scyr use glimmer_paramets, only: thk0, len0, vel0, vis0, tau0 ! use glimmer_log, only: write_log @@ -205,7 +205,7 @@ subroutine glissade_velo_sia_solve(model, & usrf => model%geometry%usrf(:,:) topg => model%geometry%topg(:,:) - bwat => model%temper%bwat(:,:) + bwat => model%basal_hydro%bwat(:,:) btrc => model%velocity%btrc(:,:) bfricflx => model%temper%bfricflx(:,:) temp => model%temper%temp(:,:,:) @@ -881,16 +881,15 @@ subroutine glissade_velo_sia_interior(nx, ny, nz, & if (stagthck(i,j) > thklim) then - siafact = 2.d0 * (rhoi*grav)**gn * stagthck(i,j)**(gn+1) & - * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((gn-1)/2) - + siafact = 2.d0 * (rhoi*grav)**n_glen * stagthck(i,j)**(n_glen+1) & + * (dusrf_dx(i,j)**2 + dusrf_dy(i,j)**2) ** ((n_glen-1)/2) vintfact(nz,i,j) = 0.d0 do k = nz-1, 1, -1 - vintfact(k,i,j) = vintfact(k+1,i,j) - & - siafact * stagflwa(k,i,j) & - * ((sigma(k) + sigma(k+1))/2.d0) ** gn & + vintfact(k,i,j) = vintfact(k+1,i,j) - & + siafact * stagflwa(k,i,j) & + * ((sigma(k) + sigma(k+1))/2.d0) ** n_glen & * (sigma(k+1) - sigma(k)) enddo ! k diff --git a/tests/slab/README.md b/tests/slab/README.md index c3767f86..71b95feb 100644 --- a/tests/slab/README.md +++ b/tests/slab/README.md @@ -1,18 +1,88 @@ Slab test case ============== -WARNING: THIS TEST CASE IS IN DEVELOPMENT AND HAS NOT BEEN SCIENTIFICALLY VALIDATED. -USE AT YOUR OWN RISK! +This directory contains python scripts for running an experiment involving a +uniform, infinite ice sheet ("slab") on an inclined plane. +The test case is described in sections 5.1-5.2 of: + Dukowicz, J. K., 2012, Reformulating the full-Stokes ice sheet model for a + more efficient computational solution. The Cryosphere, 6, 21-34, + doi:10.5194/tc-6-21-2012. -This directory contains python scripts for running an experiment involving a -uniform and infinite ice sheet ("slab") on an inclined plane. +Some results from this test case are described in Sect. 3.4 of: + Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance + of depth-integrated ice-dynamics solvers. Submitted to The Cryosphere, Aug. 2021. + +The test case consists of an ice slab of uniform thickness moving down an +inclined plane by a combination of sliding and shearing. +Analytic Stokes and first-order velocity solutions exist for all values of Glen's n >= 1. +The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1 +are derived in an unpublished manuscript by Dukowicz (2013). + +The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman +with support for Glens' n = 1. They came with warnings that the test is not supported. +The test is now supported, and the scripts include some new features: + +* The user may specify any n >= 1 (not necessarily an integer). + The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant A). +* Physics parameters are no longer hard-coded. The user can enter the ice thickness, + beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line. +* The user can specify time parameters dt (the dynamic time step) and nt (number of steps). + The previous version did not support transient runs. +* The user can specify a small thickness perturbation dh, which is added to the initial + uniform thickness via random sampling from a Gaussian distribution. + The perturbation will grow or decay, depending on the solver stability for given dx and dt. + +The run script is executed by a command like the following: + +> python runSlab.py -n 4 -a DIVA -theta 0.0375 -thk 1000. -mu 1.e5 -beta 1000. + +In this case, the user runs on 4 processors with the DIVA solver, a slope angle of 0.0375 degrees, +Glen's n = 1 (the default), slab thickness H = 1000 m, sliding coefficient beta = 1000 Pa (m/yr)^{-1}, +and viscosity coefficient 1.e5 Pa yr. +These parameters correspond to the thick shearing test case described by Robinson et al. (2021). + +To see the full set of command-line options, type 'python runSlab.py -h'. + +Notes on effective viscosity: + * For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation + mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate. + * For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n + such that the basal and surface speeds are nearly the same as for an n = 1 case with the + mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta. + * There is a subtle difference between the Dukowicz and CISM definitions of the + effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful + to make the Dukowicz convention consistent with CISM.) + +The plotting script, plotSlab.py, is run by typing 'python plotSlab.py'. It creates two plots. +The first plot shows the vertical velocity profile in nondimensional units and in units of m/yr. +There is excellent agreement between higher-order CISM solutions and the analytic solution +for small values of the slope angle theta. For steep slopes, the answers diverge as expected. + +For the second plot, the extent of the y-axis is wrong. This remains to be fixed. + +This directory also includes a new script, stabilitySlab.py, to carry out the stability tests +described in Robinson et al. (2021). + +For a given set of physics parameters and stress-balance approximation (DIVA, L1L2, etc.), +the script launches multiple CISM runs at a range of grid resolutions. +At each grid resolution, the script determines the maximum stable time step. +A run is deemed stable when the standard deviation of an initial small thickness perturbation +is reduced over the course of 100 time steps. A run is unstable if the standard deviation +increases or if the model aborts (usually with a CFL violation). + +To run the stability script, type a command like the following: + +> python stabilitySlab.py -n 4 -a DIVA -theta 0.0375 -thk 1000. -mu 1.e5 -beta 1000. \ + -dh 0.1 -nt 100 -nr 12 -rmin 10. -rmax 40000. + +Here, the first few commands correspond to the thick shearing test case and are passed repeatedly +to the run script. The remaining commands specify that each run will be initialized +with a Gaussian perturbation of amplitude 0.1 m and run for 100 timesteps. +The maximum stable timestep will be determined at 12 resolutions ranging from 10m to 40 km. +This test takes several minutes to complete on a Macbook Pro with 4 cores. -The test case is described in sections 5.1-2 of: - J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a - more efficient computational solution. The Cryosphere, 6, 21-34. - www.the-cryosphere.net/6/21/2012/ +To see the full set of commmand line options, type 'python stabilitySlab.py -h'. -Blatter-Pattyn First-order solution is described in J.K. Dukowicz, manuscript -in preparation. +For questions, please contact Willian Lipscomb (lipscomb@ucar.edu) or Gunter Leguy (gunterl@ucar.edu). diff --git a/tests/slab/plotSlab.py b/tests/slab/plotSlab.py index 214c6531..6bfa7663 100755 --- a/tests/slab/plotSlab.py +++ b/tests/slab/plotSlab.py @@ -1,10 +1,9 @@ #!/usr/bin/env python2 - """ This script plots the results of an experiment with an ice "slab" on an inclined plane. Test case is described in sections 5.1-2 of: - J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a + J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/ @@ -12,13 +11,12 @@ in preparation. """ #FIXME: Manuscript likely not in prep anymore -- JHK, 08/07/2015 +# Not published as of July 2021 -- WHL # Written by Matt Hoffman, Dec. 16, 2013 # Reconfigured by Joseph H Kennedy at ORNL on August 7, 2015 to work with the regression testing # NOTE: Did not adjust inner workings except where needed. - - -#NOTE: this script is assuming n=3, but more general solutions are available. +# Revised by William Lipscomb in 2021 to support more options, including general values of Glen's n. import os import sys @@ -28,8 +26,12 @@ import matplotlib.pyplot as plt from netCDF import * -from math import tan, pi, sin, cos -from runSlab import n, rhoi, grav, theta, beta, efvs, thickness # Get the values used to run the experiment +from math import tan, pi, sin, cos, atan + +# Get hard-coded parameters from the run script. +from runSlab import rhoi, grav + +import ConfigParser import argparse parser = argparse.ArgumentParser(description=__doc__, @@ -46,16 +48,15 @@ help="The tests output file you would like to plot. If a path is" \ +"passed via this option, the -o/--output-dir option will be ignored.") +parser.add_argument('-c','--config-file', + help="The configure file used to set up the test case and run CISM") + # =========================================================== # Define some variables and functions used in the main script # =========================================================== -# Calculate scales from Ducowicz unpub. man. -eta = beta * thickness * efvs**-n * (rhoi * grav * thickness)**(n-1) -velscale = (rhoi * grav * thickness / efvs)**n * thickness -thetar = theta * pi/180.0 # theta in radians - +#WHL args.output-file with a hyphen? def get_in_file(): if args.output_file: out_d, out_f = os.path.split(args.output_file) @@ -76,7 +77,7 @@ def get_in_file(): newest = max(matching, key=os.path.getmtime) print("\nWARNING: MULTIPLE *.out.nc FILES DETECTED!") print( "==========================================") - print( "Ploting the most recently modified file in the output directory:") + print( "Plotting the most recently modified file in the output directory:") print( " "+newest) print( "To plot another file, specify it with the -f/--outfile option.\n") @@ -94,6 +95,25 @@ def get_in_file(): return filein +def split_file_name(file_name): + """ + Get the root name, size, and number of processors from an out.nc filename. + #WHL - Adapted from plotISMIP_HOM.py + """ + root = '' + size = '' + proc = '' + + file_details = file_name.replace('.out.nc','') .split('.') +# print(file_details) +# print('len = ' + str(len(file_details))) + + if len(file_details) > 2: + proc = '.'+file_details[2] + size = '.'+file_details[1] + root = file_details[0] + + return (root, size, proc) # ========================= # Actual script starts here @@ -103,10 +123,7 @@ def main(): Plot the slab test results. """ - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - - - filein = get_in_file() + filein = get_in_file() # Get needed variables from the output file x1 = filein.variables['x1'][:] @@ -120,28 +137,96 @@ def main(): # use integer floor division operator to get an index close to the center xp = len(x0)//2 yp = len(y0)//2 - #yp = 15 - #xp = 15 # ===================================================================== - print 'Using x index of '+str(xp)+'='+str(x0[xp]) - print 'Using y index of '+str(yp)+'='+str(y0[yp]) + + print('Using x index of '+str(xp)+'='+str(x0[xp])) + print('Using y index of '+str(yp)+'='+str(y0[yp])) thk = filein.variables['thk'][:] if netCDF_module == 'Scientific.IO.NetCDF': - thk = thk * filein.variables['thk'].scale_factor + thk = thk * filein.variables['thk'].scale_factor topg = filein.variables['topg'][:] if netCDF_module == 'Scientific.IO.NetCDF': - topg = topg * filein.variables['topg'].scale_factor + topg = topg * filein.variables['topg'].scale_factor uvel = filein.variables['uvel'][:] if netCDF_module == 'Scientific.IO.NetCDF': - uvel = uvel * filein.variables['uvel'].scale_factor - + uvel = uvel * filein.variables['uvel'].scale_factor + beta_2d = filein.variables['beta'][:] + if netCDF_module == 'Scientific.IO.NetCDF': + beta_2d = beta_2d * filein.variables['beta'].scale_factor + + # Get the name of the config file + # If not entered on the command line, then construct from the output filename + + if not args.config_file: + root, size, proc = split_file_name(args.output_file) + args.config_file = root + size + proc + '.config' + + configpath = os.path.join(args.output_dir, args.config_file) + print('configpath = ' + configpath) + + # Get gn and default_flwa from the config file + + try: + config_parser = ConfigParser.SafeConfigParser() + config_parser.read( configpath ) + + gn = float(config_parser.get('parameters','n_glen')) + flwa = float(config_parser.get('parameters', 'default_flwa')) + + except ConfigParser.Error as error: + print("Error parsing " + args.config ) + print(" "), + print(error) + sys.exit(1) + + # Derive the viscosity constant mu_n from flwa + # This expression is derived in the comments on flwa in runSlab.py. + mu_n = 1.0 / (2.0**((1.0+gn)/(2.0*gn)) * flwa**(1.0/gn)) + + # Get the ice thickness from the output file. + # If thickness = constant (i.e., the optional perturbation dh = 0), it does not matter where we sample. + # Note: In general, this thickness will differ from the baseline 'thk' that is used in runSlab.py + # to create the input file. + # This is because the baseline value is measured perpendicular to the sloped bed, + # whereas the CISM value is in the vertical direction, which is not perpendicular to the bed. + thickness = thk[0,yp,xp] + + # Get beta from the output file. + # Since beta = constant, it does not matter where we sample. + beta = beta_2d[0,yp,xp] + + # Derive theta from the output file as atan(slope(topg)) + # Since the slope is constant, it does not matter where we sample. + slope = (topg[0,yp,xp] - topg[0,yp,xp+1]) / (x0[xp+1] - x0[xp]) + thetar = atan(slope) + theta = thetar * 180.0/pi + + # Compute the dimensionless parameter eta and the velocity scale, + # which appear in the scaled velocity solution. + eta = (beta * thickness / mu_n**gn) * (rhoi * grav * thickness)**(gn-1) + velscale = (rhoi * grav * thickness / mu_n)**gn * thickness + + print('gn = ' + str(gn)) + print('rhoi = ' + str(rhoi)) + print('grav = ' + str(grav)) + print('thck = ' + str(thickness)) + print('mu_n = ' + str(mu_n)) + print('flwa = ' + str(flwa)) + print('beta = ' + str(beta)) + print('eta = ' + str(eta)) + print('theta= ' + str(theta)) + print('velscale = ' + str(velscale)) # === Plot the results at the given location === # Note we are not plotting like in Fig 3 of paper. # That figure plotted a profile against zprime. # It seemed more accurate to plot a profile against z to avoid interpolating model results (analytic solution can be calculated anywhere). - # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid, so we need to transform the analytic solution to the CISM coordinate system. + # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid, + # so we need to transform the analytic solution to the CISM coordinate system. + + #WHL - I think the analytic solution is actually for u(z'), which is not bed-parallel. + # The bed-parallel solution would be u'(z'), with w'(z') = 0. fig = plt.figure(1, facecolor='w', figsize=(12, 6)) @@ -151,24 +236,23 @@ def main(): x = (x0-x0[xp]) / thickness # calculate rotated zprime coordinates for this column (we assume the solution truly is spatially uniform) zprime = x[xp] * sin(thetar) + z * cos(thetar) - #print 'zprime', zprime # Calculate analytic solution for x-component of velocity (eq. 39 in paper) for the CISM-column - #uvelStokesAnalyticScaled = sin(theta * pi/180.0) * cos(theta * pi/180.0) * (0.5 * zprime**2 - zprime - 1.0/eta) - uvelStokesAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * sin(thetar)**n * cos(thetar) / (n+1) \ - * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + sin(thetar) * cos(thetar) / eta + uvelStokesAnalyticScaled = sin(thetar) * cos(thetar) / eta \ + - 2**((1.0-gn)/2.0) * sin(thetar)**gn * cos(thetar) / (gn+1) * ( (1.0 - zprime)**(gn+1) - 1.0 ) - # Calculate the BP FO solution for x-component of velocity (Ducowicz, in prep. paper, Eq.30, n=3) - #uvelFOAnalyticScaled = (tan(theta * pi/180.0))**3 / (8.0 * (1.0 + 3.0 * (sin(theta * pi/180.0)**2))**2) \ - uvelFOAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * tan(thetar)**n / \ - ( (n + 1) * (1.0 + 3.0 * sin(thetar)**2)**((n+1.0)/2.0) ) \ - * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + tan(thetar) / eta + # Calculate the BP FO solution for x-component of velocity (Dukowicz, in prep. paper, Eq.30, n=3) + uvelFOAnalyticScaled = + tan(thetar) / eta \ + - 2**((1.0-gn)/2.0) * tan(thetar)**gn / \ + ( (gn + 1) * (1.0 + 3.0 * sin(thetar)**2)**((gn+1.0)/2.0) ) \ + * ( (1.0 - zprime)**(gn+1) - 1.0 ) ### 1. Plot as nondimensional variables # Plot analytic solution fig.add_subplot(1,2,1) plt.plot(uvelStokesAnalyticScaled, z, '-kx', label='Analytic Stokes') plt.plot(uvelFOAnalyticScaled, z, '-ko', label='Analytic FO') + # Plot model results plt.plot(uvel[0,:,yp,xp] / velscale, z, '--ro', label='CISM') plt.ylim((-0.05, 1.05)) @@ -191,7 +275,16 @@ def main(): plt.title('Velocity profile at x=' + str(x0[xp]) + ' m, y=' + str(y0[yp]) + ' m\n(Unscaled coordinates)') ################# +# print('y0_min:') +# print(y0.min()) +# print('y0_max:') +# print(y0.max()) + # Now plot maps to show if the velocities vary over the domain (they should not) + # For some reason, the y-axis has a greater extent than the range (y0.min, y0.max). + #TODO - Fix the y-axis extent. Currently, the extent is too large for small values of ny. + #TODO - Plot the thickness relative to the initial thickness. + fig = plt.figure(2, facecolor='w', figsize=(12, 6)) fig.add_subplot(1,2,1) uvelDiff = uvel[0,0,:,:] - uvel[0,0,yp,xp] @@ -224,14 +317,11 @@ def main(): #plt.plot(level, tan(thetar)**3 / (8.0 * (1.0 + 3.0 * sin(thetar)**2)**2) * (1.0 - (level-1.0)**4 ) + tan(thetar)/eta, 'b--' , label='nonlinear fo') #plt.ylim((0.0, 0.04)); plt.xlabel("z'"); plt.ylabel('u'); plt.legend() - plt.draw() plt.show() filein.close() - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - # Run only if this is being run as a script. if __name__=='__main__': @@ -240,4 +330,3 @@ def main(): # run the script sys.exit(main()) - diff --git a/tests/slab/runSlab.py b/tests/slab/runSlab.py index 2fc0217a..b6009ed5 100755 --- a/tests/slab/runSlab.py +++ b/tests/slab/runSlab.py @@ -1,6 +1,5 @@ #!/usr/bin/env python2 -#FIXME: More detailed description of this test case!!! """ Run an experiment with an ice "slab". """ @@ -8,10 +7,12 @@ # Authors # ------- # Modified from dome.py by Matt Hoffman, Dec. 16, 2013 -# Test case described in sections 5.1-2 of: -# J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a -# more efficient computational solution. The Cryosphere, 6, 21-34. www.the-cryosphere.net/6/21/2012/ -# Reconfigured by Joseph H Kennedy at ORNL on April 27, 2015 to work with the regression testing +# Test case described in sections 5.1- 5.2 of: +# J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a +# more efficient computational solution. The Cryosphere, 6, 21-34, +# https://doi.org/10.5194/tc-6-21-2012. +# Reconfigured by Joseph H Kennedy at ORNL on April 27, 2015 to work with the regression testing. +# Revised by William Lipscomb in 2021 to support more options. import os import sys @@ -19,10 +20,10 @@ import subprocess import ConfigParser -import numpy +import numpy as np import netCDF -from math import sqrt, tan, pi, cos +from math import sqrt, sin, cos, tan, pi # Parse the command line options # ------------------------------ @@ -56,11 +57,36 @@ def unsigned_int(x): parser.add_argument('-s','--setup-only', action='store_true', help="Set up the test, but don't actually run it.") - -# Additional test specific options: -#parser.add_argument('--scale', type=unsigned_int, default=0, -# help="Scales the problem size by 2**SCALE. SCALE=0 creates a 31x31 grid, SCALE=1 " -# +"creates a 62x62 grid, and SCALE=2 creates a 124x124 grid.") +# Additional options to set grid, solver, physics parameters, etc.: +#Note: The default mu_n = 0.0 is not actually used. +# Rather, mu_n is computed below, unless mu_n > 0 is specified in the command line. +# For n = 1, the default is mu_1 = 1.0e6 Pa yr. +parser.add_argument('-a','--approx', default='BP', + help="Stokes approximation (SIA, SSA, BP, L1L2, DIVA)") +parser.add_argument('-beta','--beta', default=2000.0, + help="Friction parameter beta (Pa (m/yr)^{-1})") +parser.add_argument('-dh','--delta_thck', default=0.0, + help="Thickness perturbation (m)") +parser.add_argument('-dt','--tstep', default=0.01, + help="Time step (yr)") +parser.add_argument('-gn','--glen_exponent', default=1, + help="Exponent in Glen flow law") +parser.add_argument('-l','--levels', default=10, + help="Number of vertical levels") +parser.add_argument('-mu','--mu_n', default=0.0, + help="Viscosity parameter mu_n (Pa yr^{1/n})") +parser.add_argument('-nt','--n_tsteps', default=0, + help="Number of timesteps") +parser.add_argument('-nx','--nx_grid', default=50, + help="Number of grid cells in x direction") +parser.add_argument('-ny','--ny_grid', default=5, + help="Number of grid cells in y direction") +parser.add_argument('-r','--resolution', default=100.0, + help="Grid resolution (m)") +parser.add_argument('-theta','--theta', default=5.0, + help="Slope angle (deg)") +parser.add_argument('-thk','--thickness', default=1000.0, + help="Ice thickness (m)") # Some useful functions @@ -112,28 +138,11 @@ def prep_commands(args, config_name): return commands - -# Hard coded test specific parameters # ----------------------------------- -#FIXME: Some of these could just be options! - -# Physical parameters -n = 1 # flow law parameter - only the n=1 case is currently supported -# (implementing the n=3 case would probably require implementing a new efvs option in CISM) -rhoi = 910.0 # kg/m3 -grav = 9.1801 # m^2/s - -# Test case parameters -theta = 18 # basal inclination angle (degrees) unpub. man. uses example with theta=18 -thickness = 1000.0 # m thickness in the rotated coordinate system, not in CISM coordinates +# Hard-cosed test case parameters +rhoi = 917.0 # kg/m^3 +grav = 9.81 # m^2/s baseElevation = 1000.0 # arbitrary height to keep us well away from sea level - -efvs = 2336041.42829 # hardcoded in CISM for constant viscosity setting (2336041.42829 Pa yr) - -eta = 10.0 # unpub. man. uses example with eta=10.0 -beta = eta / thickness / efvs**-n / (rhoi * grav * thickness)**(n-1) # Pa yr m^-1 -# Note: Fig. 3 in Ducowicz (2013) uses eta=18, where eta=beta*H/efvs - # the main script function # ------------------------ @@ -142,24 +151,24 @@ def main(): Run the slab test. """ - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") - # check that file name modifier, if it exists, starts with a '-' if not (args.modifier == '') and not args.modifier.startswith('-') : args.modifier = '-'+args.modifier # get the configuration # --------------------- + + dx = float(args.resolution) + dy = dx + nx = int(args.nx_grid) + ny = int(args.ny_grid) + nz = int(args.levels) + dt = float(args.tstep) + tend = float(args.n_tsteps) * dt + try: config_parser = ConfigParser.SafeConfigParser() config_parser.read( args.config ) - - nz = int(config_parser.get('grid','upn')) - nx = int(config_parser.get('grid','ewn')) - ny = int(config_parser.get('grid','nsn')) - dx = float(config_parser.get('grid','dew')) - dy = float(config_parser.get('grid','dns')) - file_name = config_parser.get('CF input', 'name') root, ext = os.path.splitext(file_name) @@ -169,7 +178,8 @@ def main(): print(error) sys.exit(1) - res = str(nx).zfill(4) + res=str(int(dx)).zfill(5) # 00100 for 100m, 01000 for 1000m, etc. + if args.parallel > 0: mod = args.modifier+'.'+res+'.p'+str(args.parallel).zfill(3) else: @@ -180,32 +190,146 @@ def main(): out_name = root+mod+'.out'+ext - # Setup the domain + # Set up the domain # ---------------- - offset = 1.0 * float(nx)*dx * tan(theta * pi/180.0) - - # create the new config file + # Create the new config file # -------------------------- if not args.quiet: print("\nCreating config file: "+config_name) + # Set grid and time parameters config_parser.set('grid', 'upn', str(nz)) config_parser.set('grid', 'ewn', str(nx)) config_parser.set('grid', 'nsn', str(ny)) config_parser.set('grid', 'dew', str(dx)) config_parser.set('grid', 'dns', str(dy)) + config_parser.set('time', 'dt', str(dt)) + config_parser.set('time', 'tend',str(tend)) + + # Set physics parameters that are needed to create the config file and the input netCDF file. + # Note: rhoi and grav are hardwired above. + + # ice thickness + thickness = float(args.thickness) + + # friction parameter beta (Pa (m/yr)^{-1}) + beta = float(args.beta) + + # basal inclination angle (degrees) + theta = float(args.theta) + theta_rad = theta * pi/180.0 # convert to radians + + # flow law exponent + # Any value n >= 1 is supported. + gn = float(args.glen_exponent) + + # Note: Fig. 3 in Dukowicz (2012) uses eta = 18 and theta = 18 deg. + # This gives u(1) = 10.0 * u(0), where u(1) = usfc and u(0) = ubas. + + # viscosity coefficient mu_n, dependent on n (Pa yr^{1/n}) + # The nominal default is mu_n = 0.0, but this value is never used. + # If a nonzero value is specified on the command line, it is used; + # else, mu_n is computed here. The goal is to choose a value mu_n(n) that + # will result in vertical shear similar to a default case with n = 1 and mu_1, + # provided we have similar values of H and theta. + # In the Dukowicz unpublished manuscript, the viscosity mu is given by + # mu = mu_n * eps_e^[(1-n)/n], where eps_e is the effective strain rate. + # For n = 1, we choose a default value of 1.0e6 Pa yr. + # For n > 1, we choose mu_n (units of Pa yr^{1/n}) to match the surface velocity + # we would have with n = 1 and the same values of H and theta. + # The general velocity solution is + # u(z') = u_b + du(z') + # where u_b = rhoi * grav * sin(theta) * cos(theta) / beta + # and du(z') = 2^{(1-n)/2}/(n+1) * sin^n(theta) * cos(theta) + # * (rhoi*grav*H/mu_n)^n * H * [1 - (1 - z'/H)^{n+1}] + # For z' = H and general n, we have + # du_n(H) = 2^{(1-n)/2}/(n+1) * sin^n(theta) * cos(theta) + # * (rhoi*grav*H/mu_n)^n * H + # For z' = H and n = 1, we have + # du_1(H) = (1/2) * sin(theta) * cos(theta) * (rhoi*grav*H/mu_1) * H + # If we equate du_n(H) with du_1(H), we can solve for mu_n: + # mu_n = [ 2^{(3-n)/(2n)}/(n+1) * sin^{n-1}(theta) * (rhoi*grav*H)^{n-1} * mu_1 ]^{1/n} + # with units Pa yr^{1/n} + # This value should give nearly the same shearing velocity du(H) for exponent n > 1 + # as we would get for n = 1, given mu_1 and the same values of H and theta. + + if float(args.mu_n) > 0.0: + mu_n = float(args.mu_n) + mu_n_pwrn = mu_n**gn + else: + mu_1 = 1.0e6 # default value for mu_1 (Pa yr) + mu_n_pwrn = 2.0**((3.0-gn)/2.0)/(gn+1.0) * sin(theta_rad)**(gn-1.0) \ + * (rhoi*grav*thickness)**(gn-1.0) * mu_1 # (mu_n)^n + mu_n = mu_n_pwrn**(1.0/gn) + + # Given mu_n, compute the temperature-independent flow factor A = 1 / [2^((1+n)/2) * mu_n^n]. + # This is how CISM incorporates a prescribed mu_n (with flow_law = 0, i.e. constant flwa). + # Note: The complicated exponent of 2 in the denominator results from CISM and the Dukowicz papers + # having different conventions for the squared effective strain rate, eps_sq. + # In CISM: mu = 1/2 * A^(-1/n) * eps_sq_c^((1-n)/(2n)) + # where eps_sq_c = 1/2 * eps_ij * eps_ij + # eps_ij = strain rate tensor + # In Dukowicz: mu = mu_n * eps_sq_d^((1-n)/(2n)) + # where eps_sq_d = eps_ij * eps_ij = 2 * eps_sq_c + # Equating the two values of mu, we get mu_n * 2^((1-n)/(2n)) = 1/2 * A^(-1/n), + # which we solve to find A = 1 / [2^((1+n)/2) * mu_n^n] + # Conversely, we have mu_n = 1 / [2^((1+n)/(2n)) * A^(1/n)] + #TODO: Modify the Dukowicz derivations to use the same convention as CISM. + flwa = 1.0 / (2.0**((1.0+gn)/2.0) * mu_n_pwrn) + + # Find the dimensionless parameter eta + # This is diagnostic only; not used directly by CISM + eta = (beta * thickness / mu_n**gn) * (rhoi * grav * thickness)**(gn-1) + + # periodic offset; depends on theta and dx + offset = 1.0 * float(nx)*dx * tan(theta_rad) + + # Print some values + print('nx = ' + str(nx)) + print('ny = ' + str(ny)) + print('nz = ' + str(nz)) + print('dt = ' + str(dt)) + print('tend = ' + str(tend)) + print('rhoi = ' + str(rhoi)) + print('grav = ' + str(grav)) + print('thck = ' + str(thickness)) + print('beta = ' + str(beta)) + print('gn = ' + str(gn)) + print('mu_n = ' + str(mu_n)) + print('flwa = ' + str(flwa)) + print('eta = ' + str(eta)) + print('theta = ' + str(theta)) + print('offset = ' + str(offset)) + + # Write some options and parameters to the config file config_parser.set('parameters', 'periodic_offset_ew', str(offset)) - + config_parser.set('parameters', 'rhoi', str(rhoi)) + config_parser.set('parameters', 'grav', str(grav)) + config_parser.set('parameters', 'n_glen', str(gn)) + config_parser.set('parameters', 'default_flwa', str(flwa)) + + if (args.approx == 'SIA'): + approx = 0 + elif (args.approx == 'SSA'): + approx = 1 + elif (args.approx == 'BP'): + approx = 2 + elif (args.approx == 'L1L2'): + approx = 3 + elif (args.approx == 'DIVA'): + approx = 4 + config_parser.set('ho_options', 'which_ho_approx', str(approx)) + config_parser.set('CF input', 'name', file_name) config_parser.set('CF output', 'name', out_name) config_parser.set('CF output', 'xtype', 'double') - + config_parser.set('CF output', 'frequency', str(tend)) # write output at start and end of run + with open(config_name, 'wb') as config_file: config_parser.write(config_file) - # create the input netCDF file # ---------------------------- if not args.quiet: @@ -222,8 +346,8 @@ def main(): nc_file.createDimension('x0',nx-1) # staggered grid nc_file.createDimension('y0',ny-1) - x = dx*numpy.arange(nx,dtype='float32') - y = dx*numpy.arange(ny,dtype='float32') + x = dx*np.arange(nx,dtype='float32') + y = dx*np.arange(ny,dtype='float32') nc_file.createVariable('time','f',('time',))[:] = [0] nc_file.createVariable('x1','f',('x1',))[:] = x @@ -231,20 +355,49 @@ def main(): nc_file.createVariable('x0','f',('x0',))[:] = dx/2 + x[:-1] # staggered grid nc_file.createVariable('y0','f',('y0',))[:] = dy/2 + y[:-1] - # Calculate values for the required variables. - thk = numpy.zeros([1,ny,nx],dtype='float32') - topg = numpy.zeros([1,ny,nx],dtype='float32') - artm = numpy.zeros([1,ny,nx],dtype='float32') - unstagbeta = numpy.zeros([1,ny,nx],dtype='float32') + thk = np.zeros([1,ny,nx],dtype='float32') + topg = np.zeros([1,ny,nx],dtype='float32') + artm = np.zeros([1,ny,nx],dtype='float32') + unstagbeta = np.zeros([1,ny,nx],dtype='float32') # Calculate the geometry of the slab of ice - thk[:] = thickness / cos(theta * pi/180.0) + # Note: Thickness is divided by cos(theta), since thck in CISM is the vertical distance + # from bed to surface. On a slanted bed, this is greater than the distance measured + # in the direction perpendicular to the bed. + # Why does topg use a tan function? Is the bed slanted? + # Do we need unstagbeta instead of beta? Compare to ISMIP-HOM tests. + + thk[:] = thickness / cos(theta_rad) xmax = x[:].max() for i in range(nx): - topg[0,:,i] = (xmax - x[i]) * tan(theta * pi/180.0) + baseElevation + topg[0,:,i] = (xmax - x[i]) * tan(theta_rad) + baseElevation unstagbeta[:] = beta + # Optionally, add a small perturbation to the thickness field + + if args.delta_thck: + dh = float(args.delta_thck) + for i in range(nx): + + # Apply a Gaussian perturbation, using the Box-Mueller algorithm: + # https://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution + + mu = 0.0 # mean of normal distribution + sigma = 1.0 # stdev of normal distribution + + rnd1 = np.random.random() # two random numbers between 0 and 1 + rnd2 = np.random.random() + + # Either of the next two lines will sample a number at random from a normal distribution + rnd_normal = mu + sigma * sqrt(-2.0 * np.log(rnd1)) * cos(2.0*pi*rnd2) +# rnd_normal = mu + sigma * sqrt(-2.0 * np.log(rnd2)) * sin(2.0*pi*rnd1) + + dthk = dh * rnd_normal + thk[0,:,i] = thk[0,:,i] + dthk + print(i, dthk, thk[0,ny/2,i]) + thk_in = thk # for comparing later to final thk + # Create the required variables in the netCDF file. nc_file.createVariable('thk', 'f',('time','y1','x1'))[:] = thk nc_file.createVariable('topg','f',('time','y1','x1'))[:] = topg @@ -274,6 +427,8 @@ def main(): print("\nRunning CISM slab test") print( "======================\n") + print('command_list =' + str(command_list)) + process = subprocess.check_call(str.join("; ",command_list), shell=True) try: @@ -289,6 +444,7 @@ def main(): if not args.quiet: print("\nFinished running the CISM slab test") print( "===================================\n") + else: run_script = args.output_dir+os.sep+root+mod+".run" @@ -304,7 +460,6 @@ def main(): print( "======================================") print( " To run the test, use: "+run_script) - print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!") # Run only if this is being run as a script. if __name__=='__main__': @@ -314,4 +469,3 @@ def main(): # run the script sys.exit(main()) - diff --git a/tests/slab/slab.config b/tests/slab/slab.config index d9ffcd61..fbba9139 100644 --- a/tests/slab/slab.config +++ b/tests/slab/slab.config @@ -1,30 +1,34 @@ [grid] -upn = 50 +upn = 20 ewn = 30 -nsn = 20 +nsn = 5 dew = 50 dns = 50 [time] tstart = 0. tend = 0. -dt = 1. +dt = 0.01 +dt_diag = 0.01 +idiag = 15 +jdiag = 5 [options] -dycore = 2 # 1 = glam, 2 = glissade -flow_law = 0 # 0 = constant +dycore = 2 # 2 = glissade +flow_law = 0 # 0 = constant flwa (default = 1.e-16 Pa-n yr-1) evolution = 3 # 3 = remapping -temperature = 1 # 1 = prognostic, 3 = enthalpy +temperature = 1 # 1 = prognostic +basal_mass_balance = 0 # 0 = basal mbal not in continuity eqn [ho_options] which_ho_babc = 5 # 5 = externally-supplied beta(required by test case) -which_ho_efvs = 0 # 0 = constant (required by test case - makes n effectively 1) -which_ho_sparse = 3 # 1 = SLAP GMRES, 3 = glissade parallel PCG, 4 = Trilinos for linear solver +which_ho_sparse = 3 # 1 = SLAP GMRES, 3 = glissade parallel PCG which_ho_nonlinear = 0 # 0 = Picard, 1 = accelerated Picard +which_ho_approx = 4 # 2 = BP, 3 = L1L2, 4 = DIVA [parameters] ice_limit = 1. # min thickness (m) for dynamics -periodic_offset_ew = 487.379544349 +geothermal = 0. [CF default] comment = created with slab.py diff --git a/tests/slab/stabilitySlab.py b/tests/slab/stabilitySlab.py new file mode 100644 index 00000000..5529896a --- /dev/null +++ b/tests/slab/stabilitySlab.py @@ -0,0 +1,387 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- + +""" +This script runs a series of CISM experiments at different resolutions. +At each resolution, it determines the maximum stable time step. +A run is deemed to be stable if the standard deviation of a small thickness perturbation +decreases during a transient run (100 timesteps by default). + +Used to obtain the CISM stability results described in: +Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance +of depth-integrated ice-dynamics solvers, to be submitted. +""" + +# Authors +# ------- +# Created by William Lipscomb, July 2021 + +import os +import sys +import errno +import subprocess +import ConfigParser + +import numpy as np +import netCDF +from math import sqrt, log10 + +# Parse the command line options +# ------------------------------ +import argparse +parser = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.ArgumentDefaultsHelpFormatter) + +# small helper function so argparse will understand unsigned integers +def unsigned_int(x): + x = int(x) + if x < 1: + raise argparse.ArgumentTypeError("This argument is an unsigned int type! Should be an integer greater than zero.") + return x + +# The following command line arguments determine the set of resolutions to run the slab test. +# At each resolution, we aim to find the maximum stable time step. +# Note: If args.n_resolution > 1, then args.resolution (see below) is ignored. + +parser.add_argument('-nr','--n_resolution', default=1, + help="number of resolutions") +parser.add_argument('-rmin','--min_resolution', default=10.0, + help="minimum resolution (m)") +parser.add_argument('-rmax','--max_resolution', default=40000.0, + help="minimum resolution (m)") + +# The following command line arguments are the same as in runSlab.py. +# Not sure how to avoid code repetition. + +parser.add_argument('-c','--config', default='./slab.config', + help="The configure file used to setup the test case and run CISM") +parser.add_argument('-e','--executable', default='./cism_driver', + help="The CISM driver") +parser.add_argument('-m', '--modifier', metavar='MOD', default='', + help="Add a modifier to file names. FILE.EX will become FILE.MOD.EX") +parser.add_argument('-n','--parallel', metavar='N', type=unsigned_int, default=0, + help="Run in parallel using N processors.") +parser.add_argument('-o', '--output_dir',default='./output', + help="Write all created files here.") +parser.add_argument('-a','--approx', default='BP', + help="Stokes approximation (SIA, SSA, BP, L1L2, DIVA)") +parser.add_argument('-beta','--beta', default=2000.0, + help="Friction parameter beta (Pa (m/yr)^{-1})") +parser.add_argument('-dh','--delta_thck', default=0.0, + help="Thickness perturbation (m)") +parser.add_argument('-dt','--tstep', default=0.01, + help="Time step (yr)") +parser.add_argument('-gn','--glen_exponent', default=1, + help="Exponent in Glen flow law") +parser.add_argument('-l','--levels', default=10, + help="Number of vertical levels") +parser.add_argument('-mu','--mu_n', default=0.0, + help="Viscosity parameter mu_n (Pa yr^{1/n})") +parser.add_argument('-nt','--n_tsteps', default=0, + help="Number of timesteps") +parser.add_argument('-nx','--nx_grid', default=50, + help="Number of grid cells in x direction") +parser.add_argument('-ny','--ny_grid', default=5, + help="Number of grid cells in y direction") +parser.add_argument('-r','--resolution', default=100.0, + help="Grid resolution (m)") +parser.add_argument('-theta','--theta', default=5.0, + help="Slope angle (deg)") +parser.add_argument('-thk','--thickness', default=1000.0, + help="Ice thickness") + + ############ + # Functions # + ############ + +def reading_file(inputFile): + + #Check whether a netCDF file exists, and return a list of times in the file + + ReadVarFile = True + try: + filein = netCDF.NetCDFFile(inputFile,'r') + time = filein.variables['time'][:] + filein.close() + print('Was able to read file ' + inputFile) + print(time) + except: + ReadVarFile = False + time = [0.] + print('Was not able to read file' + inputFile) + + return time, ReadVarFile + + +def check_output_file(outputFile, time_end): + + # Check that the output file exists with the expected final time slice + + # Path to experiment + path_to_slab_output = './output/' + + # File to check + filename = path_to_slab_output + outputFile + + # Read the output file + print('Reading file ' + str(filename)) + time_var, VarRead = reading_file(filename) + +# print(time_var) + + # Checking that the last time entry is the same as we expect from time_end + # Allow for a small roundoff difference. + if abs(time_var[-1] - time_end) < 1.0e-7: + check_time_var = True + else: + check_time_var = False + + print('time_end = ' + str(time_end)) + print('last time in file = ' + str(time_var[-1])) + + # Creating the status of both checks + check_passed = check_time_var and VarRead + + if check_passed: + print('Found output file with expected file time slice') + else: + if (not VarRead): + print('Output file cannot be read') + else: + if not check_time_var: + print('Output file is missing time slices') + + return check_passed + + +def main(): + + print('In main') + + """ + For each of several values of the horizontal grid resolution, determine the maximum + stable time step for a given configuration of the slab test. + """ + + resolution = [] + + # Based on the input arguments, make a list of resolutions at which to run the test. + # The formula and the default values of rmin and rmax give resolutions agreeing with + # those used by Alex Robinson for Yelmo, for the case nres = 12: + # resolution = [10., 21., 45., 96., 204., 434., 922., 1960., 4170., 8850., 18800., 40000.] + + print('Computing resolutions') + print(args.n_resolution) + if int(args.n_resolution) > 1: + nres = int(args.n_resolution) + resolution = [0. for n in range(nres)] + rmin = float(args.min_resolution) + rmax = float(args.max_resolution) + for n in range(nres): + res = 10.0**(log10(rmin) + (log10(rmax) - log10(rmin))*float(n)/float(nres-1)) + # Round to 3 significant figures (works for log10(res) < 5) + if log10(res) > 4.: + resolution[n] = round(res, -2) + elif log10(res) > 3.: + resolution[n] = round(res, -1) + else: + resolution[n] = round(res) + else: + nres = 1 + resolution.append(float(args.resolution)) + + print('nres = ' + str(nres)) + print(resolution) + + # Create an array to store max time step for each resolution + rows, cols = (nres, 2) + res_tstep = [[0. for i in range(cols)] for j in range(rows)] + for n in range(nres): + res_tstep[n][0] = resolution[n] + + for n in range(nres): + + print('output_dir: ' + args.output_dir) + + # Construct the command for calling the main runSlab script + run_command = 'python runSlab.py' + run_command = run_command + ' -c ' + args.config + run_command = run_command + ' -e ' + args.executable + if args.parallel > 0: + run_command = run_command + ' -n ' + str(args.parallel) + run_command = run_command + ' -o ' + args.output_dir + run_command = run_command + ' -a ' + args.approx + run_command = run_command + ' -beta ' + str(args.beta) + run_command = run_command + ' -dh ' + str(args.delta_thck) + run_command = run_command + ' -gn ' + str(args.glen_exponent) + run_command = run_command + ' -l ' + str(args.levels) + run_command = run_command + ' -mu ' + str(args.mu_n) + run_command = run_command + ' -nt ' + str(args.n_tsteps) + run_command = run_command + ' -nx ' + str(args.nx_grid) + run_command = run_command + ' -ny ' + str(args.ny_grid) + run_command = run_command + ' -theta '+ str(args.theta) + run_command = run_command + ' -thk '+ str(args.thickness) + + tend = float(args.n_tsteps) * args.tstep + + res = resolution[n] + run_command = run_command + ' -r ' + str(res) + + # Choose the time step. + # Start by choosing a very small timestep that can be assumed stable + # and a large step that can be assumed unstable. + # Note: SIA-type solvers at 10m resolution can require dt <~ 1.e-6 yr. + + tstep_lo = 1.0e-7 + tstep_hi = 1.0e+5 + tstep_log_precision = 1.0e-4 + print('Initial tstep_lo = ' + str(tstep_lo)) + print('Initial tstep_hi = ' + str(tstep_hi)) + print('Log precision = ' + str(tstep_log_precision)) + + while (log10(tstep_hi) - log10(tstep_lo)) > tstep_log_precision: + + # Compute the time step as the geometric mean of the tstep_lo and tstep_hi. + # tstep_lo is the largest time step known to be stable. + # tstep_hi is the smallest time step known to be unstable. + + tstep = sqrt(tstep_lo*tstep_hi) + + run_command_full = run_command + ' -dt ' + str(tstep) + + print("\nRunning CISM slab test...") + print('resolution = ' + str(res)) + print('tstep = ' + str(tstep)) + print('run_command = ' + run_command_full) + + process = subprocess.check_call(run_command_full, shell=True) + + print("\nFinished running the CISM slab test") + + # Determine the name of the output file. + # Must agree with naming conventions in runSlab.py + + file_name = args.config + root, ext = os.path.splitext(file_name) + + res=str(int(res)).zfill(5) # 00100 for 100m, 01000 for 1000m, etc. + + if args.parallel > 0: + mod = args.modifier + '.' + res + '.p' + str(args.parallel).zfill(3) + else: + mod = args.modifier + '.' + res + + outputFile = root + mod + '.out.nc' + + # Check whether the output file exists with the desired final time slice. + + time_end = float(args.n_tsteps) * tstep + + print('outputFile = ' + str(outputFile)) + print('n_tsteps = ' + str(float(args.n_tsteps))) + print('tstep = ' + str(tstep)) + print('time_end = ' + str(time_end)) + + check_passed = check_output_file(outputFile, time_end) + + if check_passed: + + print('Compute stdev of initial and final thickness for j = ny/2') + nx = int(args.nx_grid) + ny = int(args.ny_grid) + + # Read initial and final thickness from output file + outpath = os.path.join(args.output_dir, outputFile) + print('outpath = ' + outpath) + filein = netCDF.NetCDFFile(outpath,'r') + thk = filein.variables['thk'][:] + + j = ny/2 + thk_in = thk[0,j,:] + thk_out = thk[1,j,:] + + # Compute + Hav_in = 0.0 + Hav_out = 0.0 + for i in range(nx): + Hav_in = Hav_in + thk_in[i] + Hav_out = Hav_out + thk_out[i] + Hav_in = Hav_in / nx + Hav_out = Hav_out / nx + + # Compute + H2av_in = 0.0 + H2av_out = 0.0 + for i in range(nx): + H2av_in = H2av_in + thk_in[i]**2 + H2av_out = H2av_out + thk_out[i]**2 + H2av_in = H2av_in / nx + H2av_out = H2av_out / nx + + print('H2av_out =' + str(H2av_out)) + print('Hav_out^2 =' + str(Hav_out**2)) + + # Compute stdev = sqrt( - ^2) + var_in = H2av_in - Hav_in**2 + var_out = H2av_out - Hav_out**2 + + if var_in > 0.: + stdev_in = sqrt(H2av_in - Hav_in**2) + else: + stdev_in = 0. + + if var_out > 0.: + stdev_out = sqrt(H2av_out - Hav_out**2) + else: + stdev_out = 0. + + if stdev_in > 0.: + ratio = stdev_out/stdev_in + else: + ratio = 0. + + print('stdev_in = ' + str(stdev_in)) + print('stdev_out = ' + str(stdev_out)) + print('ratio = ' + str(ratio)) + + # Determine whether the run was stable. + # A run is defined to be stable if the final standard deviation of thickness + # is less than the initial standard deviation + + if ratio < 1.: + tstep_lo = max(tstep_lo, tstep) + print('Stable, new tstep_lo =' + str(tstep_lo)) + else: + tstep_hi = min(tstep_hi, tstep) + print('Unstable, new tstep_hi =' + str(tstep_hi)) + + else: # check_passed = F; not stable + tstep_hi = min(tstep_hi, tstep) + print('Unstable, new tstep_hi =' + str(tstep_hi)) + + print('Latest tstep_lo = ' + str(tstep_lo)) + print('Latest tstep_hi = ' + str(tstep_hi)) + + # Add to the array containing the max stable timestep at each resolution. + # Take the max stable timestep to be the average of tstep_lo and tstep_hi. + res_tstep[n][1] = 0.5 * (tstep_lo + tstep_hi) + + print('New res_tstep, res #' + str(n)) + print(res_tstep) + + # Print a table containing the max timestep for each resolution + for n in range(nres): + float_res = res_tstep[n][0] + float_dt = res_tstep[n][1] + formatted_float_res = "{:8.1f}".format(float_res) + formatted_float_dt = "{:.3e}".format(float_dt) # exponential notation with 3 decimal places + print(formatted_float_res + ' ' + formatted_float_dt) + +# Run only if this is being run as a script. +if __name__=='__main__': + + # get the command line arguments + args = parser.parse_args() + + # run the script + sys.exit(main())