From bac0a543b34e6c7456811cf577dc1b4334d229d7 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Sat, 9 Nov 2024 14:26:30 -0700 Subject: [PATCH 1/4] Changes to improve CalvingMIP accuracy Comparing to results from other groups, we found that CISM has too much retreat during the retreat phase of Experiments 2 and 4, mainly because of the method used to compute cf_length. I was computing cf_length by reconstructing the CF along a path connecting cell centers. This path can become noisy (i.e., with spurious projections and indentations), giving excessive retreat. I tried two new approaches: (1) Compute the CF length in a cell based on the number of cell edges bordering the ocean. Cells with two ocean edges have their thinning rate multiplied by sqrt(2), reflecting the longer CF (but not by a factor of 2, since the calving direction is not normal to the two edges). (2) Compute the CF length based on the secant of the angle between the calving direction and the x and y axes. The thinning rate is multiplied by 1 for the compass directions, by sqrt(2) for the diagonals, and by intermediate values in between. Method (1) does not fix the problem of excessive retreat. Method (2), however, gives more accurate retreat and is now the default for CalvingMIP. For realistic ice sheet geometries, the calving rate is not radially symmetric, so method (2) is not appropriate. Method (1) is now the default for subgrid CF schemes other than CalvingMIP. It is a bit cleaner than the previous method. Other changes: (1) Modified the method of estimating the CF location along the diagonal trajectories in the circular domain. This method uses the x0 and y0 coordinates and is a bit more accurate than the old method; it is less prone to getting stuck for several years at the same value. (2) Added a subroutine to estimate the CF location along the Halbrane and Caprona trajectories for the Thule geometry. The Caprona method is approximate and might not be as accurate as offline methods using 2D interpolation schemes. (3) The 1D fields x0 and y0 are now loadable at startup, like x1 and y1. (4) To deal with I/O involving x0 and y0, I added some code in the distributed get_var and put_var functions in parallel_mpi.F90. (5) Added a config option, which_ho_calvingmip_domain, to specify whether the user is running on the circular or Thule domain. This parameter is set automatically by the CalvingMip setup script. Note: This commit will be followed by a commit that includes a modifed CalvingMIP setup script. That script creates a grid with the origin at a cell corner instead of a cell center. --- libglide/glide_setup.F90 | 18 +- libglide/glide_types.F90 | 12 + libglide/glide_vars.def | 2 + libglimmer/parallel_mpi.F90 | 80 +- libglissade/glissade.F90 | 13 +- libglissade/glissade_calving.F90 | 1325 ++++++++++++++++++++++++------ libglissade/glissade_masks.F90 | 2 +- 7 files changed, 1162 insertions(+), 290 deletions(-) diff --git a/libglide/glide_setup.F90 b/libglide/glide_setup.F90 index 1fcfb99d..54b586fc 100644 --- a/libglide/glide_setup.F90 +++ b/libglide/glide_setup.F90 @@ -809,6 +809,7 @@ subroutine handle_ho_options(section, model) call GetValue(section, 'which_ho_assemble_bfric', model%options%which_ho_assemble_bfric) call GetValue(section, 'which_ho_assemble_lateral', model%options%which_ho_assemble_lateral) call GetValue(section, 'which_ho_calving_front', model%options%which_ho_calving_front) + call GetValue(section, 'which_ho_calvingmip_domain', model%options%which_ho_calvingmip_domain) call GetValue(section, 'which_ho_ground', model%options%which_ho_ground) call GetValue(section, 'which_ho_fground_no_glp', model%options%which_ho_fground_no_glp) call GetValue(section, 'which_ho_ground_bmlt', model%options%which_ho_ground_bmlt) @@ -1169,6 +1170,11 @@ subroutine print_options(model) 'no subgrid calving front parameterization ', & 'subgrid calving front parameterization ' /) + character(len=*), dimension(0:2), parameter :: ho_calvingmip_domain = (/ & + 'none ', & + 'circular ', & + 'Thule' /) + character(len=*), dimension(0:2), parameter :: ho_whichground = (/ & 'f_ground = 0 or 1; no GLP (glissade dycore) ', & 'GLP for basal friction with 0 <= f_ground <= 1 for vertices', & @@ -2026,7 +2032,17 @@ subroutine print_options(model) call write_log('Error, calving front option out of range for glissade dycore', GM_FATAL) end if - write(message,*) 'ho_whichground : ',model%options%which_ho_ground, & + if (model%options%which_ho_calvingmip_domain /= HO_CALVINGMIP_DOMAIN_NONE) then + write(message,*) 'ho_calvingmip_domain : ',model%options%which_ho_calvingmip_domain, & + ho_calvingmip_domain(model%options%which_ho_calvingmip_domain) + call write_log(message) + if (model%options%which_ho_calvingmip_domain < 0 .or. & + model%options%which_ho_calvingmip_domain >= size(ho_calvingmip_domain)) then + call write_log('Error, calvingMIP domain option out of range for glissade dycore', GM_FATAL) + end if + end if + + write(message,*) 'ho_whichground : ', model%options%which_ho_ground, & ho_whichground(model%options%which_ho_ground) call write_log(message) if (model%options%which_ho_ground < 0 .or. model%options%which_ho_ground >= size(ho_whichground)) then diff --git a/libglide/glide_types.F90 b/libglide/glide_types.F90 index 63e7aa82..a5a4b68e 100644 --- a/libglide/glide_types.F90 +++ b/libglide/glide_types.F90 @@ -366,6 +366,10 @@ module glide_types integer, parameter :: HO_CALVING_FRONT_NO_SUBGRID = 0 integer, parameter :: HO_CALVING_FRONT_SUBGRID = 1 + integer, parameter :: HO_CALVINGMIP_DOMAIN_NONE = 0 + integer, parameter :: HO_CALVINGMIP_DOMAIN_CIRCULAR = 1 + integer, parameter :: HO_CALVINGMIP_DOMAIN_THULE = 2 + integer, parameter :: HO_GROUND_NO_GLP = 0 integer, parameter :: HO_GROUND_GLP_BASAL_FRICTION = 1 integer, parameter :: HO_GROUND_GLP_DELUXE = 2 @@ -1090,6 +1094,14 @@ module glide_types !> \item[1] subgrid parameterization with partially filled cells at the calving front !> \end{description} + integer :: which_ho_calvingmip_domain = 0 + !> Flag that indicates the desired domain for CalvingMIP experiments + !> \begin{description} + !> \item[0] none + !> \item[1] circular (radially symmetric) + !> \item[1] Thule (complex topography) + !> \end{description} + integer :: which_ho_ground = 0 !> Flag that indicates how to compute the grounded fraction of each gridcell in the glissade dycore. !> Not valid for other dycores diff --git a/libglide/glide_vars.def b/libglide/glide_vars.def index bd776c6c..115943aa 100644 --- a/libglide/glide_vars.def +++ b/libglide/glide_vars.def @@ -23,6 +23,7 @@ long_name: Cartesian x-coordinate, velocity grid axis: X data: data%general%x0 dimlen: model%parallel%global_ewn-1 +load: 1 [y0] dimensions: y0 @@ -31,6 +32,7 @@ long_name: Cartesian y-coordinate, velocity grid axis: Y data: data%general%y0 dimlen: model%parallel%global_nsn-1 +load: 1 [x1] dimensions: x1 diff --git a/libglimmer/parallel_mpi.F90 b/libglimmer/parallel_mpi.F90 index 6df851f5..f4c878a2 100644 --- a/libglimmer/parallel_mpi.F90 +++ b/libglimmer/parallel_mpi.F90 @@ -1920,7 +1920,7 @@ function distributed_get_var_real4_1d(ncid, varid, values, parallel, start) real(sp),dimension(:) :: values type(parallel_type) :: parallel - integer :: i,ierror,myn,status,x1id,y1id + integer :: i,ierror,myn,status,x0id,y0id,x1id,y1id integer,dimension(2) :: mybounds integer,dimension(:),allocatable :: displs,sendcounts integer,dimension(:,:),allocatable :: bounds @@ -1945,14 +1945,26 @@ function distributed_get_var_real4_1d(ncid, varid, values, parallel, start) if (main_task) then allocate(bounds(2,tasks)) + status = nf90_inq_varid(ncid,"x0",x0id) + status = nf90_inq_varid(ncid,"y0",y0id) status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y1",y1id) else allocate(bounds(1,1)) end if + call broadcast(x0id) + call broadcast(y0id) call broadcast(x1id) call broadcast(y1id) - if (varid==x1id) then + if (varid==x0id) then + mybounds(1) = ewlb + mybounds(2) = ewub-1 + myn = global_ewn-1 + else if (varid==y0id) then + mybounds(1) = nslb + mybounds(2) = nsub-1 + myn = global_nsn-1 + else if (varid==x1id) then mybounds(1) = ewlb mybounds(2) = ewub myn = global_ewn @@ -1968,7 +1980,11 @@ function distributed_get_var_real4_1d(ncid, varid, values, parallel, start) if (main_task) then !WHL - See comments above on allocating the global_values array !! allocate(global_values(minval(bounds(1,:)):maxval(bounds(2,:)))) - if (varid==x1id) then + if (varid==x0id) then + allocate(global_values(global_minval_ewlb:global_maxval_ewub-1)) + elseif (varid==y0id) then + allocate(global_values(global_minval_nslb:global_maxval_nsub-1)) + elseif (varid==x1id) then allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y1id) then allocate(global_values(global_minval_nslb:global_maxval_nsub)) @@ -2105,7 +2121,7 @@ function distributed_get_var_real8_1d(ncid, varid, values, parallel, start) real(dp),dimension(:) :: values type(parallel_type) :: parallel - integer :: i,ierror,myn,status,x1id,y1id + integer :: i,ierror,myn,status,x0id,y0id,x1id,y1id integer,dimension(2) :: mybounds integer,dimension(:),allocatable :: displs,sendcounts integer,dimension(:,:),allocatable :: bounds @@ -2130,14 +2146,26 @@ function distributed_get_var_real8_1d(ncid, varid, values, parallel, start) if (main_task) then allocate(bounds(2,tasks)) + status = nf90_inq_varid(ncid,"x0",x0id) + status = nf90_inq_varid(ncid,"y0",y0id) status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y1",y1id) else allocate(bounds(1,1)) end if + call broadcast(x0id) + call broadcast(y0id) call broadcast(x1id) call broadcast(y1id) - if (varid==x1id) then + if (varid==x0id) then + mybounds(1) = ewlb + mybounds(2) = ewub-1 + myn = global_ewn-1 + else if (varid==y0id) then + mybounds(1) = nslb + mybounds(2) = nsub-1 + myn = global_nsn-1 + else if (varid==x1id) then mybounds(1) = ewlb mybounds(2) = ewub myn = global_ewn @@ -2153,7 +2181,11 @@ function distributed_get_var_real8_1d(ncid, varid, values, parallel, start) if (main_task) then !WHL - See comments above on allocating the global_values array !! allocate(global_values(minval(bounds(1,:)):maxval(bounds(2,:)))) - if (varid==x1id) then + if (varid==x0id) then + allocate(global_values(global_minval_ewlb:global_maxval_ewub-1)) + elseif (varid==y0id) then + allocate(global_values(global_minval_nslb:global_maxval_nsub-1)) + elseif (varid==x1id) then allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y1id) then allocate(global_values(global_minval_nslb:global_maxval_nsub)) @@ -3753,7 +3785,7 @@ function distributed_put_var_real4_1d(ncid, varid, values, parallel, start) type(parallel_type) :: parallel integer,dimension(:),optional :: start - integer :: i,ierror,myn,status,x0id,x1id,y0id,y1id + integer :: i,ierror,myn,status,x0id,y0id,x1id,y1id integer,dimension(2) :: mybounds integer,dimension(:),allocatable :: displs,recvcounts integer,dimension(:,:),allocatable :: bounds @@ -3777,28 +3809,28 @@ function distributed_put_var_real4_1d(ncid, varid, values, parallel, start) if (main_task) then allocate(bounds(2,tasks)) status = nf90_inq_varid(ncid,"x0",x0id) - status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y0",y0id) + status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y1",y1id) else allocate(bounds(1,1)) end if call broadcast(x0id) - call broadcast(x1id) call broadcast(y0id) + call broadcast(x1id) call broadcast(y1id) if (varid==x0id) then mybounds(1) = ewlb mybounds(2) = ewub-1 myn = global_ewn-1 - else if (varid==x1id) then - mybounds(1) = ewlb - mybounds(2) = ewub - myn = global_ewn else if (varid==y0id) then mybounds(1) = nslb mybounds(2) = nsub-1 myn = global_nsn-1 + else if (varid==x1id) then + mybounds(1) = ewlb + mybounds(2) = ewub + myn = global_ewn else if (varid==y1id) then mybounds(1) = nslb mybounds(2) = nsub @@ -3813,10 +3845,10 @@ function distributed_put_var_real4_1d(ncid, varid, values, parallel, start) !! allocate(global_values(minval(bounds(1,:)):maxval(bounds(2,:)))) if (varid==x0id) then allocate(global_values(global_minval_ewlb:global_maxval_ewub-1)) - elseif (varid==x1id) then - allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y0id) then allocate(global_values(global_minval_nslb:global_maxval_nsub-1)) + elseif (varid==x1id) then + allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y1id) then allocate(global_values(global_minval_nslb:global_maxval_nsub)) endif @@ -3948,7 +3980,7 @@ function distributed_put_var_real8_1d(ncid, varid, values, parallel, start) type(parallel_type) :: parallel integer,dimension(:),optional :: start - integer :: i,ierror,myn,status,x0id,x1id,y0id,y1id + integer :: i,ierror,myn,status,x0id,y0id,x1id,y1id integer,dimension(2) :: mybounds integer,dimension(:),allocatable :: displs,recvcounts @@ -3973,28 +4005,28 @@ function distributed_put_var_real8_1d(ncid, varid, values, parallel, start) if (main_task) then allocate(bounds(2,tasks)) status = nf90_inq_varid(ncid,"x0",x0id) - status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y0",y0id) + status = nf90_inq_varid(ncid,"x1",x1id) status = nf90_inq_varid(ncid,"y1",y1id) else allocate(bounds(1,1)) end if call broadcast(x0id) - call broadcast(x1id) call broadcast(y0id) + call broadcast(x1id) call broadcast(y1id) if (varid==x0id) then mybounds(1) = ewlb mybounds(2) = ewub-1 myn = global_ewn-1 - else if (varid==x1id) then - mybounds(1) = ewlb - mybounds(2) = ewub - myn = global_ewn else if (varid==y0id) then mybounds(1) = nslb mybounds(2) = nsub-1 myn = global_nsn-1 + else if (varid==x1id) then + mybounds(1) = ewlb + mybounds(2) = ewub + myn = global_ewn else if (varid==y1id) then mybounds(1) = nslb mybounds(2) = nsub @@ -4009,10 +4041,10 @@ function distributed_put_var_real8_1d(ncid, varid, values, parallel, start) !! allocate(global_values(minval(bounds(1,:)):maxval(bounds(2,:)))) if (varid==x0id) then allocate(global_values(global_minval_ewlb:global_maxval_ewub-1)) - elseif (varid==x1id) then - allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y0id) then allocate(global_values(global_minval_nslb:global_maxval_nsub-1)) + elseif (varid==x1id) then + allocate(global_values(global_minval_ewlb:global_maxval_ewub)) elseif (varid==y1id) then allocate(global_values(global_minval_nslb:global_maxval_nsub)) endif diff --git a/libglissade/glissade.F90 b/libglissade/glissade.F90 index 6a72e1f2..af1b80c7 100644 --- a/libglissade/glissade.F90 +++ b/libglissade/glissade.F90 @@ -2419,11 +2419,10 @@ subroutine glissade_thickness_tracer_solve(model) effective_areafrac = model%calving%effective_areafrac) if (verbose_calving) then - call point_diag(calving_front_mask, 'Before transport, calving_front_mask', itest, jtest, rtest, 7, 7) + call point_diag(calving_front_mask, 'calving_front_mask', itest, jtest, rtest, 7, 7) call point_diag(partial_cf_mask, 'partial_cf_mask', itest, jtest, rtest, 7, 7) call point_diag(full_mask, 'full_mask', itest, jtest, rtest, 7, 7) - call point_diag(ocean_mask, 'ocean_mask', itest, jtest, rtest, 7, 7) - call point_diag(model%geometry%thck*thk0, 'thck', itest, jtest, rtest, 7, 7) +! call point_diag(ocean_mask, 'ocean_mask', itest, jtest, rtest, 7, 7) call point_diag(model%calving%thck_effective, 'thck_effective', itest, jtest, rtest, 7, 7) call point_diag(model%calving%effective_areafrac, & 'effective_areafrac', itest, jtest, rtest, 7, 7, '(f10.6)') @@ -3525,7 +3524,8 @@ subroutine glissade_calving_solve(model, init_calving) nx, ny, & model%options%whichcalving, & model%options%calving_domain, & - model%options%which_ho_calving_front, & + model%options%which_ho_calving_front, & + model%options%which_ho_calvingmip_domain, & parallel, & model%calving, & ! calving object; includes calving_thck (m) itest, jtest, rtest, & @@ -3533,6 +3533,8 @@ subroutine glissade_calving_solve(model, init_calving) model%numerics%time*scyr, & ! s model%numerics%dew*len0, & ! m model%numerics%dns*len0, & ! m + model%general%x0, & ! m + model%general%y0, & ! m model%general%x1, & ! m model%general%y1, & ! m model%numerics%sigma, & @@ -4124,7 +4126,8 @@ subroutine glissade_diagnostic_variable_solve(model) thck_effective = model%calving%thck_effective, & thck_effective_min = model%calving%thck_effective_min, & partial_cf_mask = partial_cf_mask, & - full_mask = full_mask) + full_mask = full_mask, & + effective_areafrac = model%calving%effective_areafrac) ! ------------------------------------------------------------------------ ! Compute the fraction of grounded ice in each cell and at each vertex. diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 64db14df..2fa9e245 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -36,7 +36,7 @@ module glissade_calving parallel_halo, parallel_globalindex, & parallel_reduce_sum, parallel_reduce_max, parallel_reduce_log_or - use glimmer_paramets, only: eps08, thk0 + use glimmer_paramets, only: eps11, thk0 use glimmer_physcon, only: rhoi, rhoo, grav, scyr use glide_diagnostics, only: point_diag @@ -49,7 +49,8 @@ module glissade_calving glissade_stress_tensor_eigenvalues, glissade_strain_rate_tensor_eigenvalues public :: verbose_calving - logical, parameter :: verbose_calving = .false. +!! logical, parameter :: verbose_calving = .false. + logical, parameter :: verbose_calving = .true. contains @@ -234,11 +235,13 @@ subroutine glissade_calve_ice(nx, ny, & which_calving, & calving_domain, & which_ho_calving_front, & + which_ho_calvingmip_domain, & parallel, & calving, & ! calving derived type itest, jtest, rtest, & dt, time, & ! s dx, dy, & ! m + x0, y0, & ! m x1, y1, & ! m sigma, & thklim, & ! m @@ -262,6 +265,7 @@ subroutine glissade_calve_ice(nx, ny, & integer, intent(in) :: nx, ny !> horizontal grid dimensions + !TODO: Move these options to the calving derived type integer, intent(in) :: which_calving !> option for calving law integer, intent(in) :: calving_domain !> option for where calving can occur !> = 0 if calving occurs at the ocean edge only @@ -269,6 +273,7 @@ subroutine glissade_calve_ice(nx, ny, & !> = 2 if calving occurs where criterion is met and there is a connected path !> to the ocean through other cells where the criterion is met integer, intent(in) :: which_ho_calving_front !> = 1 for subgrid calving-front scheme, else = 0 + integer, intent(in) :: which_ho_calvingmip_domain !> = 1 for circular, 2 for Thule; otherwise = 0 type(parallel_type), intent(in) :: parallel !> info for parallel communication type(glide_calving), intent(inout) :: calving !> calving object @@ -301,6 +306,8 @@ subroutine glissade_calve_ice(nx, ny, & real(dp), intent(in) :: dt !> model timestep (s) real(dp), intent(in) :: time !> model time (s) real(dp), intent(in) :: dx, dy !> grid cell size in x and y directions (m) + real(dp), dimension(nx-1), intent(in) :: x0 !> x coordinates of NE cell corners (m) + real(dp), dimension(ny-1), intent(in) :: y0 !> y coordinates of NE cell corners (m) real(dp), dimension(nx), intent(in) :: x1 !> x coordinates of cell centers (m) real(dp), dimension(ny), intent(in) :: y1 !> y coordinates of cell centers (m) real(dp), dimension(:), intent(in) :: sigma !> vertical sigma coordinate @@ -370,6 +377,12 @@ subroutine glissade_calve_ice(nx, ny, & real(dp), dimension(2,8) :: & cf_location ! x and y components of calving front location ! first index is (x,y); second corresponds to 8 CalvingMIP axes + + ! some optional diagnostics + real(dp) :: & + total_ice_area, & ! total effective ice area (with weighting by effective_areafrac) + total_cf_length ! total length of the calving front + character(len=100) :: message ! initialize @@ -400,7 +413,7 @@ subroutine glissade_calve_ice(nx, ny, & !WHL - Changed definition of calving fraction; now it is the fraction lost ! rather than the fraction remaining float_fraction_calve = calving%calving_fraction - + else ! other calving options if (calving%timescale == 0.0d0) then ! calve the entire column for eligible columns (this is the default) @@ -472,6 +485,10 @@ subroutine glissade_calve_ice(nx, ny, & flux_in, & ! m^3/s parallel) + if (verbose_calving) then + call point_diag(thck, 'Before redistribution, thck (m)', itest, jtest, rtest, 7, 7) + endif + ! Gather ice that has flowed to unprotected cells and move it back upstream call redistribute_unprotected_ice(& @@ -483,6 +500,10 @@ subroutine glissade_calve_ice(nx, ny, & thck, & ! m calving%calving_thck) ! m + if (verbose_calving) then + call point_diag(thck, 'After redistribution, thck (m)', itest, jtest, rtest, 7, 7) + endif + ! Compute masks for calving. call glissade_get_masks(& @@ -513,43 +534,67 @@ subroutine glissade_calve_ice(nx, ny, & effective_areafrac = calving%effective_areafrac) if (verbose_calving) then - call point_diag(thck, 'thck (m)', itest, jtest, rtest, 7, 7) - call point_diag(calving%thck_effective, 'thck_effective (m)', itest, jtest, rtest, 7, 7) - call point_diag(calving%effective_areafrac, 'effective_areafrac', itest, jtest, rtest, 7, 7) call point_diag(calving_front_mask, 'calving_front_mask', itest, jtest, rtest, 7, 7) call point_diag(partial_cf_mask, 'partial_cf_mask', itest, jtest, rtest, 7, 7) call point_diag(full_mask, 'full_mask', itest, jtest, rtest, 7, 7) + call point_diag(calving%thck_effective, 'thck_effective (m)', itest, jtest, rtest, 7, 7) + call point_diag(calving%effective_areafrac, 'effective_areafrac', itest, jtest, rtest, 7, 7) endif ! Compute the effective length of the calving front in each grid cell - call compute_calving_front_length(& - nx, ny, & - dx, dy, & - itest, jtest, rtest, & - parallel, & - calving_front_mask, & - cf_length) + if (which_calving == CF_ADVANCE_RETREAT_RATE) then + + ! compute the CF length as a function of a cell's location on the unit circle + ! surrounding the origin (assuming a radially symmetric calving rate). + + call compute_calving_front_length_radial(& + nx, ny, & + dx, dy, & + x1, y1, & + itest, jtest, rtest, & + calving_front_mask, & + ocean_mask, & + cf_length) + + else + + ! compute the CF length for each cell based on its number of ocean neighbors + ! (i.e., enhanced calving for cells with 2 or 3 ocean neighbors). + + call compute_calving_front_length(& + nx, ny, & + dx, dy, & + itest, jtest, rtest, & + calving_front_mask, & + ocean_mask, & + cf_length) + + endif ! which_calving call parallel_halo(cf_length, parallel) if (verbose_calving) then - call point_diag(cf_length, 'cf_length(m)', itest, jtest, rtest, 7, 7) + call point_diag(cf_length, 'cf_length (m)', itest, jtest, rtest, 7, 7) + ! Diagnose the total CF length + total_cf_length = 0.0d0 + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + if (calving_front_mask(i,j) == 1) then + total_cf_length = total_cf_length + cf_length(i,j) + endif + enddo + enddo + total_cf_length = parallel_reduce_sum(total_cf_length) + if (this_rank == rtest) then + print*, 'Total CF length (km)', total_cf_length/1000.d0 + endif endif - if (which_calving == CF_ADVANCE_RETREAT_RATE) then - - ! Diagnose the current calving front location, based on the CalvingMIP domain. - ! (The advance/retreat option is used for several CalvingMIP experiments.) + ! Depending on the calving method, compute the calving rate for each grid cell + ! and convert to an equivalent thinning rate. - call locate_calving_front_calvingMIP(& - nx, ny, & - dx, dy, & ! m - x1, y1, & ! m - parallel, & - itest, jtest, rtest, & - calving%effective_areafrac, & - cf_location) ! m + if (which_calving == CF_ADVANCE_RETREAT_RATE) then call calving_front_advance_retreat(& nx, ny, & @@ -567,8 +612,6 @@ subroutine glissade_calve_ice(nx, ny, & elseif (which_calving == CALVING_THCK_THRESHOLD) then - !TODO - Pass in cf_length? - call thickness_based_calving(& nx, ny, & dx, dy, & ! m @@ -666,7 +709,6 @@ subroutine glissade_calve_ice(nx, ny, & endif ! Apply calving_dthck as computed above. - ! Adjust the thinning based on the estimated calving-front length in the cell. call apply_calving_dthck(& nx, ny, & @@ -686,9 +728,6 @@ subroutine glissade_calve_ice(nx, ny, & call parallel_halo(thck, parallel) call parallel_halo(calving%calving_thck, parallel) - !TODO - Is it possible to skip the call to advance_calving_front - ! and the two mask calls? The mask calls are needed to compute thck_effective. - ! Recompute the calving masks call glissade_get_masks(& @@ -715,22 +754,106 @@ subroutine glissade_calve_ice(nx, ny, & thck_effective = calving%thck_effective, & thck_effective_min = calving%thck_effective_min, & partial_cf_mask = partial_cf_mask, & - full_mask = full_mask) + full_mask = full_mask, & + effective_areafrac = calving%effective_areafrac) + + if (verbose_calving) then + call point_diag(thck, 'After calving dthck: thck (m)', itest, jtest, rtest, 7, 7) +! call point_diag(calving%thck_effective, 'thck_effective (m)', itest, jtest, rtest, 7, 7) +! call point_diag(calving%effective_areafrac, 'effective_areafrac', itest, jtest, rtest, 7, 7) +! call point_diag(thck/calving%effective_areafrac, 'CF thickness', itest, jtest, rtest, 7, 7) + endif - ! Where thck > thck_effective, allow the CF to advance by distributing ice downstream + if (which_calving == CF_ADVANCE_RETREAT_RATE) then - call advance_calving_front(& - nx, ny, & - itest, jtest, rtest, & - ocean_mask, & - calving_front_mask, & - flux_in, & - calving%thck_effective, & - thck) + ! Compute some CalvingMIP diagnostics. + ! Note: If running with a prescribed advance/retreat rate on a grid other than + ! the CalvingMIP circular and Thule domains, we would need some additional + ! logic to identify the domain. + ! Note: With a prescribed advance/retreat rate, it isn't necessary to call advance_calving_front. + + if (which_ho_calvingmip_domain == HO_CALVINGMIP_DOMAIN_CIRCULAR) then + + call locate_calving_front_circular(& + nx, ny, & + dx, dy, & ! m + x0, y0, & ! m + x1, y1, & ! m + parallel, & + itest, jtest, rtest, & + calving%effective_areafrac, & + cf_location) ! m + + elseif (which_ho_calvingmip_domain == HO_CALVINGMIP_DOMAIN_THULE) then + + call locate_calving_front_thule(& + nx, ny, & + dx, dy, & ! m + x0, y0, & ! m + x1, y1, & ! m + parallel, & + itest, jtest, rtest, & + calving%effective_areafrac, & + cf_location) ! m + + endif + + ! Compute the total ice area and the area of each quadrant + total_ice_area = 0.0d0 + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + total_ice_area = total_ice_area + dx*dy*calving%effective_areafrac(i,j) + enddo + enddo + total_ice_area = parallel_reduce_sum(total_ice_area) + if (this_rank == rtest) then + print*, 'Total ice area (km^2)=', total_ice_area/1.0d6 + print*, 'Quadrant area (km^2)=', total_ice_area/4.0d6 + endif + + else ! other calving schemes + + ! Where thck > thck_effective, allow the CF to advance by distributing ice downstream. + ! Note: This is not necessary when running with a prescribed advance/retreat rate. + + call advance_calving_front(& + nx, ny, & + itest, jtest, rtest, & + ocean_mask, & + calving_front_mask, & + flux_in, & + calving%thck_effective, & + thck) + + ! Recompute the calving masks + + call glissade_get_masks(& + nx, ny, & + parallel, & + thck, topg, & + eus, thklim, & + ice_mask, & + floating_mask = floating_mask, & + ocean_mask = ocean_mask, & + land_mask = land_mask) + + call glissade_calving_front_mask(& + nx, ny, & + which_ho_calving_front, & + parallel, & + thck, topg, & + eus, & + ice_mask, floating_mask, & + ocean_mask, land_mask, & + calving_front_mask, & + dthck_dx_cf = calving%dthck_dx_cf, & + dx = dx, dy = dy, & + thck_effective = calving%thck_effective, & + thck_effective_min = calving%thck_effective_min, & + partial_cf_mask = partial_cf_mask, & + full_mask = full_mask, & + effective_areafrac = calving%effective_areafrac) - if (verbose_calving) then - call point_diag(thck, 'After advance: thck (m)', itest, jtest, rtest, 7, 7) - call point_diag(calving%thck_effective, 'thck_effective (m)', itest, jtest, rtest, 7, 7) endif else ! other calving options (no subgrid calving front) @@ -900,11 +1023,6 @@ subroutine redistribute_unprotected_ice(& ! I think the halo update is needed only to get the right halo values for diagnostics call parallel_halo(thck, parallel) - if (verbose_calving) then - call point_diag(thck, 'Before redistribution, thck (m)', itest, jtest, rtest, 7, 7) - call point_diag(protected_mask, 'protected_mask', itest, jtest, rtest, 7, 7) - endif - ! Identify unprotected ice with nonzero thickness. ! Instead of calving this ice, move it to one or more protected upstream CF cell ! (from which most or all of the ice likely arrived during transport). @@ -953,221 +1071,200 @@ subroutine redistribute_unprotected_ice(& call parallel_halo(thck, parallel) - ! Write the final thicknesses - if (verbose_calving) then - call point_diag(thck, 'After redistribution, thck (m)', itest, jtest, rtest, 7, 7) - endif - end subroutine redistribute_unprotected_ice !--------------------------------------------------------------------------- - subroutine locate_calving_front_calvingMIP(& - nx, ny, & - dx, dy, & - x1, y1, & - parallel, & - itest, jtest, rtest, & - areafrac, & - cf_location) + subroutine compute_calving_front_length(& + nx, ny, & + dx, dy, & + itest, jtest, rtest, & + calving_front_mask, & + ocean_mask, & + cf_length) - use cism_parallel, only: parallel_reduce_maxloc, parallel_reduce_minloc, broadcast + ! Compute the effective length of the calving front in each grid cell, + ! based on the number of ocean neighbors. + ! Cells with a single ocean edge receive a length of dx or dy. + ! Cells with two or three adjacent edges receive a length of sqrt(dx^2 + dy^2). + ! + ! input/output arguments integer, intent(in) :: & nx, ny, & ! grid dimensions itest, jtest, rtest ! coordinates of diagnostic point real(dp), intent(in) :: & - dx, dy ! grid cell size (m) - - real(dp), dimension(nx), intent(in) :: x1 - real(dp), dimension(ny), intent(in) :: y1 + dx, dy ! grid cell size (m) - type(parallel_type), intent(in) :: & - parallel ! info for parallel communication + integer, dimension(nx,ny), intent(in) :: & + calving_front_mask, & ! = 1 for floating cells with an ice-free ocean neighber + ocean_mask ! = 1 for ice-free ocean cells - real(dp), dimension(nx,ny), intent(in) :: areafrac - real(dp), dimension(2,8), intent(out) :: cf_location + real(dp), dimension(nx,ny), intent(out) :: & + cf_length ! calving front length (m) through a given grid cell ! local variables - integer :: i, j, iglobal, jglobal - integer :: axis - integer :: procnum - real(dp) :: cf_location_xmax, cf_location_ymax, cf_location_xmin, cf_location_ymin, radius - - ! Find the x and y coordinates of the calving front along the different axes - ! specified in CalvingMIP. - ! The code assumes a circular domain with center at (0,0). - !TODO: Do this processing offline with a Python script + integer :: i, j, ip, jp, ii, jj + integer :: count, count_ew, count_ns ! number of neighbors of a CF cell - if (this_rank == rtest) print*, 'Locate_calving_front for calvingMIP, rtest =', rtest + ! Compute the CF length based on the number of ocean neighbors; + ! the CF is longer for cells with two ocean neighbors. - cf_location(:,:) = 0.0d0 + cf_length = 0.0d0 - ! Find the CF location along each of 8 axes - ! The CF lies in the last cell along a given axis with areafrac > 0 - ! All loops are over locally owned cells + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + if (calving_front_mask(i,j) == 1) then - axis = 1 ! index for the positive y-axis - do i = nhalo+1, nx-nhalo - if (x1(i) == 0.0d0) then - cf_location(1,axis) = 0.0d0 - do j = nhalo+1, ny-nhalo - if (areafrac(i,j) > 0.0d0 .and. areafrac(i,j+1) == 0.0d0) then - cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)*dy + ! Count the number of ocean edges in each direction + count_ns = ocean_mask(i,j-1) + ocean_mask(i,j+1) ! N/S neighbors; edge of length dx + count_ew = ocean_mask(i-1,j) + ocean_mask(i+1,j) ! E/W neighbors; edge of length dy + count = count_ns + count_ew + if (count == 1) then ! one ocean neighbor + cf_length(i,j) = count_ns*dx + count_ew*dy + elseif (count == 2) then ! two ocean neighbors; usually a corner in the CF + if (count_ew == 1 .and. count_ns == 1) then + cf_length(i,j) = sqrt(dx*dx + dy*dy) + else ! rare case of an isthmus with ocean on opposite edges + cf_length(i,j) = count_ns*dx + count_ew*dy + endif + elseif (count == 3) then + ! This is a bit arbitrary; assume the same length as a CF bordering two edges + cf_length(i,j) = sqrt(dx*dx + dy*dy) endif - enddo - endif - enddo - - ! If this proc has a positive value of y, then broadcast the coordinates to all procs - call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) - axis = 2 ! index for the line y = x in the positive x and y direction - do i = nhalo+1, nx-nhalo - do j = nhalo+1, ny-nhalo - if (x1(i) == y1(j)) then ! on the line y = x - if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j+1) == 0.0d0) then - cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)*dx - cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)*dy - endif - endif + endif ! calving_mask enddo ! i enddo ! j - ! If this proc has a positive value of x, then broadcast the coordinates to all procs - call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + end subroutine compute_calving_front_length - axis = 3 ! index for the positive x-axis - do j = nhalo+1, ny-nhalo - if (y1(j) == 0.0d0) then - cf_location(2,axis) = 0.0d0 - do i = nhalo+1, nx-nhalo - if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j) == 0.0d0) then - cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)*dx - endif - enddo - endif - enddo +!--------------------------------------------------------------------------- - ! If this proc has a positive value of x, then broadcast the coordinates to all procs - call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_xmax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + subroutine compute_calving_front_length_radial(& + nx, ny, & + dx, dy, & + x1, y1, & + itest, jtest, rtest, & + calving_front_mask, & + ocean_mask, & + cf_length) - axis = 4 ! index for the line y = -x in the positive x and negative y direction - do i = nhalo+1, nx-nhalo - do j = nhalo+1, ny-nhalo - if (x1(i) == -y1(j)) then ! on the line y = -x - if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j-1) == 0.0d0) then - cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)*dx - cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))*dy - endif - endif - enddo ! i - enddo ! j + ! This is an alternate method of computing the calving front length in idealized problems + ! (e.g., CalvingMIP) in which the calving rate has radial symmetry. + ! For this method, the CF length for a given gridcell has a minimum value = dx) when a line drawn + ! from the origin to the cell center is parallel to the x or y axis. + ! The CF length has a maximum value when a line drawn from the origin to the cell center + ! makes an angle of pi/4 = 45 degrees with the x and y axes. + ! At intermediate angles, the CF length has an intermediate value. - ! If this proc has a positive value of x, then broadcast the coordinates to all procs - call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + ! input/output arguments - axis = 5 ! index for the negative y-axis - do i = nhalo+1, nx-nhalo - if (x1(i) == 0.0d0) then - cf_location(1,axis) = 0.0d0 - do j = ny-nhalo, nhalo+1, -1 - if (areafrac(i,j) > 0.0d0 .and. areafrac(i,j-1) == 0.0d0) then - cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))*dy - endif - enddo - endif - enddo + integer, intent(in) :: & + nx, ny, & ! grid dimensions + itest, jtest, rtest ! coordinates of diagnostic point - ! If this proc has a negative value of y, then broadcast the coordinates to all procs - call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + real(dp), intent(in) :: & + dx, dy ! grid cell size (m) - axis = 6 ! index for the line y = x in the negative x and y direction - do i = nhalo+1, nx-nhalo - do j = nhalo+1, ny-nhalo - if (x1(i) == y1(j)) then ! on the line y = x - if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j-1) == 0.0d0) then - cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))*dx - cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))*dy - endif - endif - enddo ! i - enddo ! j + real(dp), dimension(nx), intent(in) :: x1 ! x coordinate of cell centers + real(dp), dimension(ny), intent(in) :: y1 ! y coordinate of cell centers - ! If this proc has a negative value of x, then broadcast the coordinates to all procs - call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + integer, dimension(nx,ny), intent(in) :: & + calving_front_mask, & ! = 1 for floating cells with an ice-free ocean neighber + ocean_mask ! = 1 for ice-free ocean cells - axis = 7 ! index for the negative x-axis - do j = nhalo+1, ny-nhalo - if (y1(j) == 0.0d0) then - cf_location(2,axis) = 0.0d0 - do i = nx-nhalo, nhalo+1, -1 - if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j) == 0.0d0) then - cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))*dx - endif - enddo - endif - enddo + real(dp), dimension(nx,ny), intent(out) :: & + cf_length ! calving front length (m) through a given grid cell - ! If this proc has a negative value of x, then broadcast the coordinates to all procs - call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_xmin, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + ! local variables - axis = 8 ! index for the line y = -x in the negative x and positive y direction - do i = nhalo+1, nx-nhalo - do j = nhalo+1, ny-nhalo - if (x1(i) == -y1(j)) then ! on the line y = -x - if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j+1) == 0.0d0) then - cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))*dx - cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)*dy - endif - endif - enddo ! i - enddo ! j + integer :: i, j - ! If this proc has a negative value of x, then broadcast the coordinates to all procs - call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) - call broadcast(cf_location(:,axis), proc=procnum) + real(dp) :: absx, absy ! absolute value of (x1,y1) + real(dp) :: theta ! angle relative to x or y axis - if (verbose_calving .and. main_task) then - print*, ' ' - print*, 'axis, CF location, radius (km):' - do axis = 1, 8 - radius = sqrt(cf_location(1,axis)**2 + cf_location(2,axis)**2) - write(6,'(i4,3f10.3)') axis, cf_location(:,axis)/1000.d0, radius/1000.d0 - enddo + if (abs(dx - dy) > eps11) then + write(6,*) 'Error, compute_calving_front_length_radial, this_rank, i, j,', this_rank, i, j + call write_log('Must have dx = dy for the cf_length method', GM_FATAL) endif - end subroutine locate_calving_front_calvingMIP + ! Compute the CF length for each grid cell: + ! * equal to dx when the calving direction is parallel to the x or y axis + ! * equal to sqrt(2)*dx when the calving direction is along a diagonal + ! * of intermediate length for intermediate angles + ! Note: The method assumes dx = dy. + ! It fails if the center of a CF cell lies at the origin (0,0). + + do j = nhalo+1, ny-nhalo + do i = nhalo+1, nx-nhalo + if (calving_front_mask(i,j) == 1) then + absx = abs(x1(i)) + absy = abs(y1(j)) + if (absx == 0.0d0 .and. absy == 0.0d0) then + write(6,*) 'Error, compute_calving_front_length_radial, this_rank, i, j,', this_rank, i, j + call write_log('Cannot compute angle for CF cell at the origin', GM_FATAL) + else + if (absx >= absy) then + theta = atan(absy/absx) + else ! absx < absy + theta = atan(absx/absy) + endif + endif + ! Note: theta lies in the range [0, pi/4], so cos(theta) is in the range [sqrt(2)/2, 1] + cf_length(i,j) = dx / cos(theta) + endif ! calving front cell + enddo ! i + enddo ! j + + end subroutine compute_calving_front_length_radial !--------------------------------------------------------------------------- - subroutine compute_calving_front_length(& + subroutine compute_calving_front_length_advance(& nx, ny, & dx, dy, & itest, jtest, rtest, & - parallel, & calving_front_mask, & cf_length) - use glimmer_physcon, only: pi - - ! Compute the effective length of the calving front in each grid cell. - ! The method is as follows: - ! * For each CF cell, identify its CF neighbors (usually 2 neighbors). - ! * Connect the center of the cell with the center of each neighbor. - ! * Compute the distance associated with these connections: - ! dx/2 or dy/2 for edge neighbors and sqrt(dx^2+dy^2)/2 for corner neighbors - - ! input/output arguments + ! Compute the calving front length by a method slightly different from the method + ! in subroutine compute_calving_front_length. + ! This method is based on connecting all the adjacent cell centers on the calving front, + ! such that cells with diagonal CF neighbors thin more than cells with edge neighbors. + ! + ! This was originally the method used for CalvingMIP experiments, + ! before switching to the more accurate radial method. + ! When applied to CalvingMIP problems, it is more accurate than the method + ! in subroutine compute_calving_front_length during the advance phases, + ! but less accurate during the retreat phase. + ! + ! A picture can be helpful. Let 'x' denote ice-filled cells and let 'o' denote ocean cells. + ! + ! | o | o | o | o | + ! |--------|--------|--------|--------| + ! | | | | | + ! | x | x | o | o | + ! | | (i,j+1)| | | + ! |--------|--------|--------|--------| + ! | | | | | + ! | x | x | x | x | + ! | | (i,j) | (i+1,j)| | + ! |------- |--------|--------|--------| + ! + ! Note the corner in the CF at the NE vertex of cell (i,j). + ! In general, we want calving to smooth out the CF rather than created teeth and indentations. + ! During retreat, we can smooth the corner by enhancing the calving in cell (i,j+1), + ! which borders two ocean cells. We do this by giving this cell an increased CF length. + ! But when the CF is prescribed to advance, it's better to smooth the corner by reducing the + ! amount of ice entering ocean cell (i+1,j+1). We do this by increasing the CF length in both + ! (i,j+1) and (i+1,j), to reflect the diagonal connection between the two cells. + ! This subroutine is not currently called but is here for now in case it's useful later. + !--------------------------------------------------------------------------- + + ! input/output variables integer, intent(in) :: & nx, ny, & ! grid dimensions @@ -1176,9 +1273,6 @@ subroutine compute_calving_front_length(& real(dp), intent(in) :: & dx, dy ! grid cell size (m) - type(parallel_type), intent(in) :: & - parallel ! info for parallel communication - integer, dimension(nx,ny), intent(in) :: & calving_front_mask ! = 1 for floating cells with an ice-free ocean neighber @@ -1187,19 +1281,13 @@ subroutine compute_calving_front_length(& ! local variables - integer :: i, j, ip, jp, ii, jj, ig, jg + integer :: i, j, ip, jp, ii, jj integer :: count - real(dp) :: & - local_cf_length, & ! length of CF on local domain - global_cf_length ! length of CF on global domain - - ! Initialize - cf_length = 0.0d0 - do j = nhalo+1, ny-nhalo do i = nhalo+1, nx-nhalo if (calving_front_mask(i,j) == 1) then + count = 0 ! start with edge CF neighbors do jp = -1,1 @@ -1216,6 +1304,7 @@ subroutine compute_calving_front_length(& endif enddo ! ip enddo ! jp + ! check for corner CF neighbors ! Note: Occasionally, a cell can have 3 CF neighbors. ! If so, then count corner neighbors only if there are not already 2 edge neighbors @@ -1225,45 +1314,28 @@ subroutine compute_calving_front_length(& ii = i + ip; jj = j + jp if (calving_front_mask(ii,jj) == 1) then count = count + 1 - if (count <= 2) cf_length(i,j) = cf_length(i,j) + 0.5d0*sqrt(dx*dx + dy*dy) + if (count <= 2) then + cf_length(i,j) = cf_length(i,j) + 0.5d0*sqrt(dx*dx + dy*dy) + endif endif enddo enddo endif + ! if we have not found 2 neighbors, then add a CF length of sqrt(dx*dy)/2 for each missing neighbor if (count == 0) then cf_length(i,j) = sqrt(dx*dy) elseif (count == 1) then cf_length(i,j) = cf_length(i,j) + 0.5d0*sqrt(dx*dy) endif - !WHL - debug - if (cf_length(i,j) > 1.01d0*sqrt(dx*dx + dy*dy)) then - call parallel_globalindex(i, j, ig, jg, parallel) - print*, 'Warning: ig, jg, CF length:', ig, jg, cf_length(i,j) - endif + endif ! calving_mask enddo ! i enddo ! j - ! Diagnose the total CF length - if (verbose_calving) then - local_cf_length = 0.0d0 - do j = nhalo+1, ny-nhalo - do i = nhalo+1, nx-nhalo - if (calving_front_mask(i,j) == 1) then - local_cf_length = local_cf_length + cf_length(i,j) - endif - enddo - enddo - global_cf_length = parallel_reduce_sum(local_cf_length) - if (this_rank == rtest) then - print*, 'Total CF length (km)', global_cf_length/1000.d0 - endif - endif + end subroutine compute_calving_front_length_advance - end subroutine compute_calving_front_length - - !--------------------------------------------------------------------------- +!--------------------------------------------------------------------------- subroutine thickness_based_calving(& nx, ny, & @@ -1855,20 +1927,13 @@ subroutine calving_front_advance_retreat(& ! If so, it is possible to undo this calving so the CF can advance. ! Decrease the calving if the CF is advancing, increase if the CF is retreating + calving_dthck(i,j) = calving_dthck(i,j) - & (cf_advance_retreat_rate*dt * thck_effective(i,j) * cf_length(i,j)) / (dx*dy) ! Make sure dthck >= 0. This prevents unphysical CF advance via negative calving calving_dthck(i,j) = max(calving_dthck(i,j), 0.0d0) - if (verbose_calving .and. i==itest .and. j==jtest .and. this_rank==rtest) then - write(6,*) ' ' - write(6,*) 'thck:', thck(i,j) - write(6,*) 'thck_effective (m) =', thck_effective(i,j) - write(6,*) 'dthck transport (m) =', dthck_dt_transport(i,j) * dt - write(6,*) 'dthck calving (m) =', calving_dthck(i,j) - endif - endif ! calving_front cell enddo ! i enddo ! j @@ -2876,6 +2941,748 @@ subroutine extrapolate_to_calving_front(& end subroutine extrapolate_to_calving_front +!--------------------------------------------------------------------------- +! The next two subroutines are diagnostic subroutines for CalvingMIP. +! They estimate the calving front location along 8 prescribed axes +! for the circular and Thule domains. +! They are not necessary if we have offline tools for locating the CF, +! but are left here for reference. +!--------------------------------------------------------------------------- + + subroutine locate_calving_front_circular(& + nx, ny, & + dx, dy, & + x0, y0, & + x1, y1, & + parallel, & + itest, jtest, rtest, & + areafrac, & + cf_location) + + use cism_parallel, only: parallel_reduce_maxloc, parallel_reduce_minloc, broadcast + use glissade_grid_operators, only: glissade_stagger + use glide_diagnostics, only: point_diag + + ! Find the calving front location along eight profiles on the circular domain. + ! These profiles are the four cardinal directions (N, S, E, W) along with the diagonals + ! that form 45-degree angles with the cardinal directions. + + integer, intent(in) :: & + nx, ny, & ! grid dimensions + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid cell size (m) + + real(dp), dimension(nx-1), intent(in) :: x0 ! x coordinate of NE cell corners + real(dp), dimension(ny-1), intent(in) :: y0 ! y coordinate of NE cell corners + real(dp), dimension(nx), intent(in) :: x1 ! x coordinate of cell centers + real(dp), dimension(ny), intent(in) :: y1 ! y coordinate of cell centers + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + real(dp), dimension(nx,ny), intent(in) :: areafrac + real(dp), dimension(2,8), intent(out) :: cf_location + + ! local variables + + integer :: i, j, iglobal, jglobal + integer :: axis + integer :: procnum + real(dp) :: cf_location_xmax, cf_location_ymax, cf_location_xmin, cf_location_ymin, radius + real(dp) :: & + this_areafrac_avg, next_areafrac_avg ! average of areafrac in two adjacent cells + real(dp) :: areafrac_ne, areafrac_nw, areafrac_se, areafrac_sw + + ! Note: For the original CalvingMIP grid, the origin was located at a cell center, + ! so both axes passed through cell centers. + ! For the new CalvingMIP grid (as of Nov. 2024), the origin is located at a cell corner, + ! so both axes lie along cell edges. + ! The following code is written generally to find the CF along the x-axis and y-axis + ! in either case. + ! The code aborts if one of the other is not true. + + logical :: & + x_axis_thru_centers, & ! true if the x-axis passes through cell centers + x_axis_thru_edges, & ! true if the x-axis passes through cell edges + y_axis_thru_centers, & ! true if the y-axis passes through cell centers + y_axis_thru_edges ! true if the y-axis passes through cell edges + + ! Find the x and y coordinates of the calving front along the different axes + ! specified in CalvingMIP. + ! The code assumes a circular domain with center at (0,0). + ! The logic depends on whether the N, S, E and W axes pass through cell centers or edges. + + if (this_rank == rtest) print*, 'Locate_calving_front for calvingMIP, rtest =', rtest + + ! Determine whether the x and y axes passes through cell centers, or through cell edges. + ! They should pass through one or the other. + x_axis_thru_centers = .false. + do j = nhalo+1, ny-nhalo + if (y1(j) == 0.0d0) then + x_axis_thru_centers = .true. + endif + enddo + + x_axis_thru_edges = .false. + do j = nhalo+1, ny-nhalo + if (y0(j) == 0.0d0) then + x_axis_thru_edges = .true. + endif + enddo + + y_axis_thru_centers = .false. + do i = nhalo+1, nx-nhalo + if (x1(i) == 0.0d0) then + y_axis_thru_centers = .true. + endif + enddo + + y_axis_thru_edges = .false. + do i = nhalo+1, nx-nhalo + if (x0(i) == 0.0d0) then + y_axis_thru_edges = .true. + endif + enddo + +! if (x_axis_thru_centers) then +! print*, this_rank, 'x_axis_thru_centers', x_axis_thru_centers +! endif +! if (y_axis_thru_centers) then +! print*, this_rank, 'y_axis_thru_centers', y_axis_thru_centers +! endif +! if (x_axis_thru_edges) then +! print*, this_rank, 'x_axis_thru_edges', x_axis_thru_edges +! endif +! if (y_axis_thru_edges) then +! print*, this_rank, 'y_axis_thru_edges', y_axis_thru_edges +! endif + + cf_location(:,:) = 0.0d0 + + ! Find the CF location along each of 8 axes + ! The CF lies in the last cell along a given axis with areafrac > 0 + ! All loops are over locally owned cells + + axis = 1 ! index for the positive y-axis (profile A) + if (y_axis_thru_centers) then + do i = nhalo+1, nx-nhalo + if (x1(i) == 0.0d0) then + cf_location(1,axis) = 0.0d0 + do j = nhalo+1, ny-nhalo + if (areafrac(i,j) > 0.0d0 .and. areafrac(i,j+1) == 0.0d0) then + cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)*dy + endif + enddo + endif + enddo + elseif (y_axis_thru_edges) then + do i = nhalo+1, nx-nhalo + if (x0(i) == 0.0d0) then ! E edge of cell lies on the y-axis + cf_location(1,axis) = 0.0d0 + do j = nhalo+1, ny-nhalo + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j+1) + areafrac(i+1,j+1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (this_areafrac_avg - 0.5d0)*dy + endif + enddo + endif + enddo + endif ! y_axis_thru_centers + + ! If this proc has a positive value of y, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 2 ! index for the line y = x in the positive x and y direction (profile B) + do i = nhalo+1, nx-nhalo + do j = nhalo+1, ny-nhalo + if (x1(i) == y1(j)) then ! on the line y = x + if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j+1) == 0.0d0) then + areafrac_ne = 0.5d0 * (areafrac(i+1,j) + areafrac(i,j+1)) + areafrac_sw = 0.5d0 * (areafrac(i,j-1) + areafrac(i-1,j)) + if (areafrac_ne >= 0.5d0) then ! CF in cell (i+1,j+1) + cf_location(1,axis) = x0(i) + (areafrac_ne - 0.5d0)/areafrac_ne * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (areafrac_ne - 0.5d0)/areafrac_ne * (0.5d0*dy) + elseif (areafrac_sw < 0.5d0) then ! CF in cell (i-1,j-1) + cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_sw)/(areafrac(i-1,j-1) - areafrac_sw) * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_sw)/(areafrac(i-1,j-1) - areafrac_sw) * (0.5d0*dy) + else ! CF in cell (i,j) + if (areafrac(i,j) >= 0.5d0) then ! CF in upper right of cell + cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_ne) * (0.5d0*dx) + cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_ne) * (0.5d0*dy) + else ! areafrac(i,j) < 0.5; CF in lower left of cell + cf_location(1,axis) = x1(i) - (0.5d0 - areafrac(i,j))/(areafrac_sw - areafrac(i,j)) * (0.5d0*dx) + cf_location(2,axis) = y1(j) - (0.5d0 - areafrac(i,j))/(areafrac_sw - areafrac(i,j)) * (0.5d0*dy) + endif + endif + endif + endif ! on the line y = x + enddo ! i + enddo ! j + + ! If this proc has a positive value of x, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 3 ! index for the positive x-axis (profile C) + if (x_axis_thru_centers) then + do j = nhalo+1, ny-nhalo + if (y1(j) == 0.0d0) then + cf_location(2,axis) = 0.0d0 + do i = nhalo+1, nx-nhalo + if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j) == 0.0d0) then + cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)*dx + endif + enddo + endif + enddo + elseif (x_axis_thru_edges) then + do j = nhalo+1, ny-nhalo + if (y0(j) == 0.0d0) then + cf_location(2,axis) = 0.0d0 + do i = nhalo+1, nx-nhalo + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i,j+1)) + next_areafrac_avg = 0.5d0 * (areafrac(i+1,j) + areafrac(i+1,j+1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(1,axis) = x1(i) + (this_areafrac_avg - 0.5d0)*dx + endif + enddo + endif + enddo + endif ! x_axis_thru_centers + + ! If this proc has a positive value of x, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_xmax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 4 ! index for the line y = -x in the positive x and negative y direction (profile D) + do i = nhalo+1, nx-nhalo + do j = nhalo+1, ny-nhalo + if (x1(i) == -y1(j)) then ! on the line y = -x + if (areafrac(i,j) > 0.0d0 .and. areafrac(i+1,j-1) == 0.0d0) then + areafrac_se = 0.5d0 * (areafrac(i+1,j) + areafrac(i,j-1)) + areafrac_nw = 0.5d0 * (areafrac(i-1,j) + areafrac(i,j+1)) + if (areafrac_se >= 0.5d0) then ! CF in cell (i+1,j-1) + cf_location(1,axis) = x0(i) + (areafrac_se - 0.5d0)/areafrac_se * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (areafrac_se - 0.5d0)/areafrac_se * (0.5d0*dy) + elseif (areafrac_nw < 0.5d0) then ! CF in cell (i-1,j+1) + cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_nw)/(areafrac(i-1,j-1) - areafrac_nw) * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_nw)/(areafrac(i-1,j-1) - areafrac_nw) * (0.5d0*dy) + else ! CF in cell (i,j) + if (areafrac(i,j) >= 0.5d0) then ! CF in lower right of cell + cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_se) * (0.5d0*dx) + cf_location(2,axis) = y1(j) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_se) * (0.5d0*dy) + else ! areafrac(i,j) < 0.5; CF in upper left of cell + cf_location(1,axis) = x1(i) - (0.5d0 - areafrac(i,j))/(areafrac_nw - areafrac(i,j)) * (0.5d0*dx) + cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))/(areafrac_nw - areafrac(i,j)) * (0.5d0*dy) + endif + endif + endif + endif + enddo ! i + enddo ! j + + ! If this proc has a positive value of x, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 5 ! index for the negative y-axis (profile E) + if (y_axis_thru_centers) then + do i = nhalo+1, nx-nhalo + if (x1(i) == 0.0d0) then + cf_location(1,axis) = 0.0d0 + do j = ny-nhalo, nhalo+1, -1 + if (areafrac(i,j) > 0.0d0 .and. areafrac(i,j-1) == 0.0d0) then + cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))*dy + endif + enddo + endif + enddo + elseif (y_axis_thru_edges) then + do i = nhalo+1, nx-nhalo + if (x0(i) == 0.0d0) then ! E edge of cell lies on the y-axis + cf_location(1,axis) = 0.0d0 + do j = ny-nhalo, nhalo+1, -1 + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j-1) + areafrac(i+1,j-1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (0.5d0 - this_areafrac_avg)*dy + endif + enddo + endif + enddo + endif ! y_axis_thru_centers + + ! If this proc has a negative value of y, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 6 ! index for the line y = x in the negative x and y direction (profile F) + do i = nhalo+1, nx-nhalo + do j = nhalo+1, ny-nhalo + if (x1(i) == y1(j)) then ! on the line y = x + if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j-1) == 0.0d0) then + areafrac_sw = 0.5d0 * (areafrac(i,j-1) + areafrac(i-1,j)) + areafrac_ne = 0.5d0 * (areafrac(i,j+1) + areafrac(i+1,j)) + if (areafrac_sw >= 0.5d0) then ! CF in cell (i-1,j-1) + cf_location(1,axis) = x0(i-1) - (areafrac_sw - 0.5d0)/areafrac_sw * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (areafrac_sw - 0.5d0)/areafrac_sw * (0.5d0*dy) + elseif (areafrac_ne < 0.5d0) then ! CF in cell (i+1,j+1) + cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_ne)/(areafrac(i-1,j-1) - areafrac_ne) * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_ne)/(areafrac(i-1,j-1) - areafrac_ne) * (0.5d0*dy) + else ! CF in cell (i,j) + if (areafrac(i,j) >= 0.5d0) then ! CF in lower left of cell + cf_location(1,axis) = x1(i) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_sw) * (0.5d0*dx) + cf_location(2,axis) = y1(j) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_sw) * (0.5d0*dy) + else ! areafrac(i,j) < 0.5; CF in upper right of cell + cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))/(areafrac_ne - areafrac(i,j)) * (0.5d0*dx) + cf_location(2,axis) = y1(j) + (0.5d0 - areafrac(i,j))/(areafrac_ne - areafrac(i,j)) * (0.5d0*dy) + endif + endif + endif + endif + enddo ! i + enddo ! j + + ! If this proc has a negative value of x, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 7 ! index for the negative x-axis (profile g) + if (x_axis_thru_centers) then + do j = nhalo+1, ny-nhalo + if (y1(j) == 0.0d0) then + cf_location(2,axis) = 0.0d0 + do i = nx-nhalo, nhalo+1, -1 + if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j) == 0.0d0) then + cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))*dx + endif + enddo + endif + enddo + elseif (x_axis_thru_edges) then + do j = nhalo+1, ny-nhalo + if (y0(j) == 0.0d0) then + cf_location(2,axis) = 0.0d0 + do i = nx-nhalo, nhalo+1, -1 + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i,j+1)) + next_areafrac_avg = 0.5d0 * (areafrac(i-1,j) + areafrac(i-1,j+1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(1,axis) = x1(i) + (0.5d0 - this_areafrac_avg)*dx + endif + enddo + endif + enddo + endif ! x_axis_thru_centers + + ! If this proc has a negative value of x, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_xmin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 8 ! index for the line y = -x in the negative x and positive y direction (profile H) + do i = nhalo+1, nx-nhalo + do j = nhalo+1, ny-nhalo + if (x1(i) == -y1(j)) then ! on the line y = -x + if (areafrac(i,j) > 0.0d0 .and. areafrac(i-1,j+1) == 0.0d0) then + areafrac_nw = 0.5d0 * (areafrac(i-1,j) + areafrac(i,j+1)) + areafrac_se = 0.5d0 * (areafrac(i+1,j) + areafrac(i,j-1)) + if (areafrac_nw >= 0.5d0) then ! CF in cell (i-1,j+1) + cf_location(1,axis) = x0(i-1) - (areafrac_nw - 0.5d0)/areafrac_nw * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (areafrac_nw - 0.5d0)/areafrac_nw * (0.5d0*dy) + elseif (areafrac_se < 0.5d0) then ! CF in cell (i+1,j-1) + cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_se)/(areafrac(i+1,j-1) - areafrac_se) * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_se)/(areafrac(i+1,j-1) - areafrac_se) * (0.5d0*dy) + else ! CF in cell (i,j) + if (areafrac(i,j) >= 0.5d0) then ! CF in upper left of cell + cf_location(1,axis) = x1(i) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_nw) * (0.5d0*dx) + cf_location(2,axis) = y1(j) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_nw) * (0.5d0*dy) + else ! areafrac(i,j) < 0.5; CF in lower right of cell + cf_location(1,axis) = x1(i) + (0.5d0 - areafrac(i,j))/(areafrac_se - areafrac(i,j)) * (0.5d0*dx) + cf_location(2,axis) = y1(j) - (0.5d0 - areafrac(i,j))/(areafrac_se - areafrac(i,j)) * (0.5d0*dy) + endif + endif + endif + endif + enddo ! i + enddo ! j + + ! If this proc has a negative value of x, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(1,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + if (verbose_calving .and. main_task) then + print*, ' ' + print*, 'Circular domain: axis, CF location, radius (km)' + do axis = 1, 8 + radius = sqrt(cf_location(1,axis)**2 + cf_location(2,axis)**2) + write(6,'(i4,3f10.3)') axis, cf_location(:,axis)/1000.d0, radius/1000.d0 + enddo + endif + + end subroutine locate_calving_front_circular + +!--------------------------------------------------------------------------- + + subroutine locate_calving_front_thule(& + nx, ny, & + dx, dy, & + x0, y0, & + x1, y1, & + parallel, & + itest, jtest, rtest, & + areafrac, & + cf_location) + + use cism_parallel, only: parallel_reduce_maxloc, parallel_reduce_minloc, broadcast + use glissade_grid_operators, only: glissade_stagger + use glide_diagnostics, only: point_diag + + ! Find the calving front location along eight profiles on the Thule domain. + ! These profiles are defined as follows: + ! + ! Halbrane profiles: + ! A: (-150,0) to (-150, 740) + ! B: (150, 0) to ( 150, 740) + ! C: (-150,0) to (-150,-740) + ! D: (150, 0) to ( 150,-740) + ! + ! Caprona profiles: + ! A: (-390,0) to (-590, 450) + ! B: (390,0) to ( 590, 450) + ! C: (-390,0) to (-590,-450) + ! D: (390,0) to ( 590,-450) + ! + ! The Halbrane profiles are easier. Not sure if I'm going to implement the Caprona profiles. + + integer, intent(in) :: & + nx, ny, & ! grid dimensions + itest, jtest, rtest ! coordinates of diagnostic point + + real(dp), intent(in) :: & + dx, dy ! grid cell size (m) + + real(dp), dimension(nx-1), intent(in) :: x0 ! x coordinate of NE cell corners + real(dp), dimension(ny-1), intent(in) :: y0 ! y coordinate of NE cell corners + real(dp), dimension(nx), intent(in) :: x1 ! x coordinate of cell centers + real(dp), dimension(ny), intent(in) :: y1 ! y coordinate of cell centers + + type(parallel_type), intent(in) :: & + parallel ! info for parallel communication + + real(dp), dimension(nx,ny), intent(in) :: areafrac + + real(dp), dimension(2,8), intent(out) :: & + cf_location ! x and y locations of CF along the Halbrane and Caprona profiles + ! first index: x and y + ! second index: 1 to 4 for Halbrane, 5 to 8 for Caprona + + ! local variables + + integer :: i, j, jj, iglobal, jglobal + integer :: axis + integer :: procnum + real(dp) :: cf_location_xmax, cf_location_ymax, cf_location_xmin, cf_location_ymin, radius + real(dp) :: this_areafrac, next_areafrac + real(dp) :: & + this_areafrac_avg, next_areafrac_avg ! average of areafrac in two adjacent cells + + real(dp), dimension(nx) :: y_int, areafrac_int + real(dp) :: x_intercept, y_intercept, slope ! properties of the profile + real(dp) :: x_lim, y_lim ! outer limits of the profile + real(dp) :: dist_y, frac_dist + + ! Find the x and y coordinates of the calving front for the Thule domain + ! along the different profiles specified in CalvingMIP. + + cf_location(:,:) = 0.0d0 + + ! Find the CF location along Halbrane profiles A, B, C and D. + ! All loops are over locally owned cells. + ! Assume that the x value of each profile coincides with a cell edge + ! (not a cell center). + + axis = 1 ! index for Halbrane A + x_intercept = -150.d3 + + do i = nhalo+1, nx-nhalo + if (x0(i) == x_intercept) then ! E edge of cell lies on the vertical Halbrane profile + cf_location(1,axis) = x_intercept + do j = nhalo+1, ny-nhalo + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j+1) + areafrac(i+1,j+1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (this_areafrac_avg - 0.5d0)*dy + endif + enddo + endif + enddo + + ! If this proc has a positive value of y, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 2 ! index for Halbrane B (same as A except for positive x_intercept) + x_intercept = 150.d3 + + do i = nhalo+1, nx-nhalo + if (x0(i) == x_intercept) then ! E edge of cell lies on the vertical Halbrane profile + cf_location(1,axis) = x_intercept + do j = nhalo+1, ny-nhalo + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j+1) + areafrac(i+1,j+1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (this_areafrac_avg - 0.5d0)*dy + endif + enddo + endif + enddo + + ! If this proc has a positive value of y, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymax, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 3 ! index for Halbrane C (same as A except in the negative y direction) + x_intercept = -150.d3 + + do i = nhalo+1, nx-nhalo + if (x0(i) == x_intercept) then ! E edge of cell lies on the vertical Halbrane profile + cf_location(1,axis) = x_intercept + do j = ny-nhalo, nhalo+1, -1 + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j-1) + areafrac(i+1,j-1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (0.5d0 - this_areafrac_avg)*dy + endif + enddo + endif + enddo + + ! If this proc has a negative value of y, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 4 ! index for Halbrane D (same as C except for positive x_intercept) + x_intercept = 150.d3 + + do i = nhalo+1, nx-nhalo + if (x0(i) == x_intercept) then ! E edge of cell lies on the vertical Halbrane profile + cf_location(1,axis) = x_intercept + do j = ny-nhalo, nhalo+1, -1 + this_areafrac_avg = 0.5d0 * (areafrac(i,j) + areafrac(i+1,j)) + next_areafrac_avg = 0.5d0 * (areafrac(i,j-1) + areafrac(i+1,j-1)) + if (this_areafrac_avg > 0.0d0 .and. next_areafrac_avg == 0.0d0) then + cf_location(2,axis) = y1(j) + (0.5d0 - this_areafrac_avg)*dy + endif + enddo + endif + enddo + + ! If this proc has a negative value of y, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + ! Find the CF location along Caprona profiles A, B, C and D. + ! The Caprona profiles cut across cells without passing through centers or corners. + ! As a result, the logic below is more complicated than for the Halbrane profiles, + ! and more approximate. Results might be better with offline interpolation of areafrac. + + axis = 5 ! index for Caprona A + x_intercept = -390.d3 + x_lim = -590.d3 + y_lim = 450.d3 + slope = y_lim/(x_lim - x_intercept) ! rise over run = 450/(-200) = -2.25 + y_intercept = -x_intercept * slope + + ! Adjust x_lim to allow the CF to be a little out of bounds + x_lim = -650.d3 + + y_int = 0.0d0 + areafrac_int = 0.0d0 + + ! Estimate areafrac at each point where the Caprona profile intersects the x1 grid + do i = nx-nhalo, nhalo+1, -1 + if (x1(i) < x_intercept .and. x1(i) >= x_lim) then ! x1 in range + y_int(i) = slope*x1(i) + y_intercept ! profile intersects x1 grid at (x1(i),y) + do j = nhalo+1, ny-nhalo + if (y_int(i) >= y1(j) .and. y_int(i) < y1(j+1)) then + ! Interpolate to estimate a_eff at (x1(i),y_int) + areafrac_int(i) = areafrac(i,j) + (y_int(i) - y1(j))/dy * (areafrac(i,j+1) - areafrac(i,j)) + exit + endif + enddo + endif + enddo + + ! Find a point along the profile where the interpolated areafrac = 0.5 + do i = nx-nhalo, nhalo+1, -1 + if (areafrac_int(i) > 0.5d0 .and. areafrac_int(i-1) < 0.5d0) then + dist_y = y_int(i-1) - y_int(i) ! y distance between neighboring intersection points + frac_dist = (areafrac_int(i) - 0.5d0) / (areafrac_int(i) - areafrac_int(i-1)) + cf_location(1,axis) = x1(i) - frac_dist*dx + cf_location(2,axis) = y_int(i) + frac_dist*dist_y + exit + endif + enddo + + ! If this proc has a positive value of y, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 6 ! index for Caprona B + x_intercept = 390.d3 + x_lim = 590.d3 + y_lim = 450.d3 + slope = y_lim/(x_lim - x_intercept) ! rise over run = 450/200 = 2.25 + y_intercept = -x_intercept * slope + + ! Adjust x_lim to allow the CF to be a little out of bounds + x_lim = 650.d3 + + y_int = 0.0d0 + areafrac_int = 0.0d0 + + ! Estimate areafrac at each point where the Caprona profile intersects the x1 grid + do i = nhalo+1, nx-nhalo + if (x1(i) >= x_intercept .and. x1(i) < x_lim) then ! x1 in range + y_int(i) = slope*x1(i) + y_intercept ! profile intersects x1 grid at (x1(i),y) + do j = nhalo+1, ny-nhalo + if (y_int(i) >= y1(j) .and. y_int(i) < y1(j+1)) then + ! Interpolate to estimate a_eff at (x1(i),y_int) + areafrac_int(i) = areafrac(i,j) + (y_int(i) - y1(j))/dy * (areafrac(i,j+1) - areafrac(i,j)) + exit + endif + enddo + endif + enddo + + if (verbose_calving .and. this_rank ==rtest) then +! print*, 'Caprona B intersection points: i, x, y, areafrac' +! do i = nhalo+1, nx-nhalo +! if (y_int(i) /= 0.0d0) then +! write(6,'(i4,3f10.3)'), i, x1(i)/1000.d0, y_int(i)/1000.d0, areafrac_int(i) +! endif +! enddo + endif + + ! Find a point along the profile where the interpolated areafrac = 0.5 + do i = nhalo+1, nx-nhalo + if (areafrac_int(i) > 0.5d0 .and. areafrac_int(i+1) < 0.5d0) then + dist_y = y_int(i+1) - y_int(i) ! y distance between neighboring intersection points + frac_dist = (areafrac_int(i) - 0.5d0) / (areafrac_int(i) - areafrac_int(i+1)) + cf_location(1,axis) = x1(i) + frac_dist*dx + cf_location(2,axis) = y_int(i) + frac_dist*dist_y + if (verbose_calving .and. this_rank == rtest) then +! print*, '1st IP: x, y, a_eff =', x1(i)/1000.d0, y_int(i)/1000.d0, areafrac_int(i) +! print*, '2nd IP: x, y, a_eff =', x1(i+1)/1000.d0, y_int(i+1)/1000.d0, areafrac_int(i+1) +! print*, 'dist_y, frac_dist =', dist_y/1000.d0, frac_dist +! print*, 'CF location =', cf_location(1,axis)/1000.d0, cf_location(2,axis)/1000.d0 +! print*, 'residual y - (mx + b):', cf_location(2,axis) - slope*cf_location(1,axis) - y_intercept + endif + exit + endif + enddo + + ! If this proc has a positive value of y, then broadcast the coordinates to all procs + call parallel_reduce_maxloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 7 ! index for Caprona C + x_intercept = -390.d3 + x_lim = -590.d3 + y_lim = -450.d3 + slope = y_lim/(x_lim - x_intercept) ! rise over run = -450/(-200) = 9/4 + y_intercept = -x_intercept * slope + + ! Adjust x_lim to allow the CF to be a little out of bounds + x_lim = -650.d3 + + y_int = 0.0d0 + areafrac_int = 0.0d0 + + ! Estimate areafrac at each point where the Caprona profile intersects the x1 grid + do i = nx-nhalo, nhalo+1, -1 + if (x1(i) < x_intercept .and. x1(i) >= x_lim) then ! x1 in range + y_int(i) = slope*x1(i) + y_intercept ! profile intersects x1 grid at (x1(i),y) + do j = ny-nhalo, nhalo+1, -1 + if (y_int(i) <= y1(j) .and. y_int(i) > y1(j-1)) then + ! Interpolate to estimate a_eff at (x1(i),y_int) + areafrac_int(i) = areafrac(i,j) + (y1(j) - y_int(i))/dy * (areafrac(i,j-1) - areafrac(i,j)) + exit + endif + enddo + endif + enddo + + ! Find a point along the profile where the interpolated areafrac = 0.5 + do i = nx-nhalo, nhalo+1, -1 + if (areafrac_int(i) > 0.5d0 .and. areafrac_int(i-1) < 0.5d0) then + dist_y = y_int(i-1) - y_int(i) ! y distance between neighboring intersection points + frac_dist = (areafrac_int(i) - 0.5d0) / (areafrac_int(i) - areafrac_int(i-1)) + cf_location(1,axis) = x1(i) - frac_dist*dx + cf_location(2,axis) = y_int(i) + frac_dist*dist_y + exit + endif + enddo + + ! If this proc has a negative value of y, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + axis = 8 ! index for Caprona D + x_intercept = 390.d3 + x_lim = 590.d3 + y_lim = -450.d3 + slope = y_lim/(x_lim - x_intercept) ! rise over run = -450/200 = -9/4 + y_intercept = -x_intercept * slope + + ! Adjust x_lim to allow the CF to be a little out of bounds + x_lim = 650.d3 + + y_int = 0.0d0 + areafrac_int = 0.0d0 + + ! Estimate areafrac at each point where the Caprona profile intersects the x1 grid + do i = nhalo+1, nx-nhalo + if (x1(i) >= x_intercept .and. x1(i) < x_lim) then ! x1 in range + y_int(i) = slope*x1(i) + y_intercept ! profile intersects x1 grid at (x1(i),y) + do j = ny-nhalo, nhalo+1, -1 + if (y_int(i) <= y1(j) .and. y_int(i) > y1(j-1)) then + ! Interpolate to estimate a_eff at (x1(i),y_int) + areafrac_int(i) = areafrac(i,j) + (y1(j) - y_int(i))/dy * (areafrac(i,j-1) - areafrac(i,j)) + exit + endif + enddo + endif + enddo + + ! Find a point along the profile where the interpolated areafrac = 0.5 + do i = nhalo+1, nx-nhalo + if (areafrac_int(i) > 0.5d0 .and. areafrac_int(i+1) < 0.5d0) then + dist_y = y_int(i+1) - y_int(i) ! y distance between neighboring intersection points + frac_dist = (areafrac_int(i) - 0.5d0) / (areafrac_int(i) - areafrac_int(i+1)) + cf_location(1,axis) = x1(i) + frac_dist*dx + cf_location(2,axis) = y_int(i) + frac_dist*dist_y + endif + enddo + + ! If this proc has a negative value of y, then broadcast the coordinates to all procs + call parallel_reduce_minloc(xin=cf_location(2,axis), xout=cf_location_ymin, xprocout=procnum) + call broadcast(cf_location(:,axis), proc=procnum) + + if (verbose_calving .and. main_task) then + print*, ' ' + print*, 'Thule domain: axis, CF location, radius (km)' + do axis = 1, 8 + radius = sqrt(cf_location(1,axis)**2 + cf_location(2,axis)**2) + write(6,'(i4,3f10.3)') axis, cf_location(:,axis)/1000.d0, radius/1000.d0 + enddo + endif + + end subroutine locate_calving_front_thule + !--------------------------------------------------------------------------- end module glissade_calving diff --git a/libglissade/glissade_masks.F90 b/libglissade/glissade_masks.F90 index 914bd986..d9990829 100644 --- a/libglissade/glissade_masks.F90 +++ b/libglissade/glissade_masks.F90 @@ -321,7 +321,7 @@ subroutine glissade_calving_front_mask(& real(dp), dimension(nx,ny), intent(out), optional :: & thck_effective, & ! effective ice thickness (m) for calving ! Generally, H_eff > H at the CF, with H_eff = H elsewhere - effective_areafrac ! effective ice-covered fraction, in range [0,1] + effective_areafrac ! effective ice-covered fraction, in range [0,1] ! 0 < f < 1 for partial calving-front cells integer, dimension(nx,ny), intent(out), optional :: & From db8e06b679ce2b5cbd9f93197551ca5116706c6f Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Sat, 9 Nov 2024 14:30:48 -0700 Subject: [PATCH 2/4] Modified CalvingMIP setup script and template This commit modifies the CalvingMIP setup script. The main change is that the grid origin is now located at a cell corner instead of a cell center. This means that on a 5-km grid measuring 1600 km along each edge, there are now 320 cells in each direction instead of 321. This change makes it easier to initialize the calving front at a distance of 750 km from the origin, instead of being extended or falling short by half a grid cell. Other changes: * Fixed a bug that gave the wrong retreat rate for Experiment 4 * Added the which_ho_calvingmip_domain variable to the config files * Edited the README file for clarity and accuracy --- tests/calvingMIP/README.calvingMIP | 25 ++++++--- tests/calvingMIP/calvingMIP.Setup.py | 60 ++++++++++----------- tests/calvingMIP/calvingMIP.config.template | 9 ++-- 3 files changed, 50 insertions(+), 44 deletions(-) diff --git a/tests/calvingMIP/README.calvingMIP b/tests/calvingMIP/README.calvingMIP index 1fcf46d6..38e9170c 100644 --- a/tests/calvingMIP/README.calvingMIP +++ b/tests/calvingMIP/README.calvingMIP @@ -4,6 +4,7 @@ Note: For setting up the experiments in an NCAR computing environment, follow the steps in the README.NCAR_HPC file in the tests directory. See this wiki for details on MISMIP+: +https://github.com/JRowanJordan/CalvingMIP/wiki There are two experimental domains for CalvingMIP: * a circular domain, which is radially symmetric @@ -17,6 +18,9 @@ Experiment3 1000 years on the Thule domain with a fixed calving front Experiment4 1000 years on the Thule domain with a prescribed rate of retreat and advance Experiment5 1000 years with a thickness-based calving front +Most groups ran the first four experiments. There were some issues with Experiment 5, +so it wasn't included in the model intercomparison. + To initialize these experiments, we also need a long spin-up on each domain. This directory contains two files needed to set up the experiments: @@ -31,10 +35,10 @@ for the different config files. These setup files are located in directory ../tests/calvingMIP. You might want to keep this directory clean and set up your tests in a duplicate directory. -For example, from ../tests/calvingMIP: +For example, from ../tests/ (i.e., one directory up from calvingMIP): -> cp -rf calvingMIP ../calvingMIP.test -> cd ../calvingMIP.test +> cp -rf calvingMIP calvingMIP.test +> cd calvingMIP.test The following instructions are for running calvingMIP on derecho: @@ -78,18 +82,23 @@ Notes on optional arguments: - If the resolution is relatively high (e.g., 5000 m or finer), it may take several minutes to create the input files. +To write output more frequently (e.g., annual output for Experiments 1 to 4), +you can set outputfreqhi and outputfreqmd when running the setup script:: + +> python calvingMIP.Setup.py -fm 1. -fh 1. + The config template file includes sensible default values for other parameters. -To change any of these parameters, you need to edit the template file. +To change any of these parameters, you will need to edit the template file. After setup, you should first run the SpinupCircular and SpinupThule experiments. To do this, you will need a run script in each directory. The test directory does not include a run script, -but the following is a sample run script for Derecho with a Land Ice Working Group account code: +The following is a sample run script for Derecho with a Land Ice Working Group account code: ***** #!/bin/bash # #PBS -N calvingMIP -#PBS -A P93300601 +#PBS -A P93300301 #PBS -l walltime=01:00:00 #PBS -q main #PBS -l job_priority=regular @@ -101,11 +110,11 @@ but the following is a sample run script for Derecho with a Land Ice Working Gro . derecho-intel-modules -mpiexec -n 128 -ppn 128 ./cism_driver calvingMIP.SpinupCircular.config > cism.out.$PBS_JOBID 2>&1 +mpiexec -n 128 -ppn 128 ./cism_driver calvingMIP.ExperimentName.config > cism.out.$PBS_JOBID 2>&1 ***** On 128 Derecho cores, the required walltime is typically about 30 minutes for a 10kyr spin-up, -and less than 5 minutes for the 1kry experiments. +and less than 5 minutes for the 1kyr experiments. You can launch the job by typing diff --git a/tests/calvingMIP/calvingMIP.Setup.py b/tests/calvingMIP/calvingMIP.Setup.py index 0ae6b547..f2aac7b7 100644 --- a/tests/calvingMIP/calvingMIP.Setup.py +++ b/tests/calvingMIP/calvingMIP.Setup.py @@ -26,11 +26,11 @@ def unsigned_int(x): # Constants # ############# -R = 800.e3 # (m) size of the circular and Thule domains; x and y range over (-R, R) +R = 800000. # (m) size of the circular and Thule domains; x and y range over (-R, R) # For the circular domain (Experiments 1 and 2) -Bc_circ = 900. # (m) maximum bed topography elevation at the center of the domain -Bl_circ = -2000. # (m) minimum bed topography elevation at the edge of the domain +Bc_circ = 900. # (m) maximum bed topography elevation at the center of the domain +Bl_circ = -2000. # (m) minimum bed topography elevation at the edge of the domain # For the Thule domain (Experiments 3, 4 and 5) Bc_thule = 900. # (m) maximum bed topography elevation at the center of the domain @@ -200,9 +200,10 @@ def createFileFromSource(src_file, trg_file): print( 'Number of vertical levels =', nz) # Set number of grid cells in each direction. -# E.g., if R = 800 km with dx = dy = 10 km, we have 161 cells in each direction -nx = 2*int(R/dx) + 1 -ny = 2*int(R/dy) + 1 +# E.g., if R = 800 km with dx = dy = 10 km, we have 160 cells in each direction +# Note: The origin lies at a cell corner +nx = 2*int(R/dx) +ny = 2*int(R/dy) print('grid dimension in x:',nx) print('grid dimension in y:',ny) @@ -236,10 +237,10 @@ def createFileFromSource(src_file, trg_file): config.set('time', 'dt_diag', str(args.timestep)) # Set the diagnostic cell -# By default, this cell lies at the center of the x-axis, near the +# By default, this cell lies just to the right of the x-axis, near the # northern boundary of the domain. idiag = int(nx/2) + 1 -jdiag = int(idiag*31./16.) - 1 +jdiag = int(nx/2) + int(rcalve/dy) config.set('time', 'idiag', str(idiag)) config.set('time', 'jdiag', str(jdiag)) print('idiag:',idiag) @@ -329,18 +330,16 @@ def createFileFromSource(src_file, trg_file): # Compute x and y on each grid. # The origin (x = y = 0) is placed at the center of the domain. -x = np.arange(-R,R+dx,dx,dtype='float32') # x = -R, -R+dx,..., -dx, 0, dx,..., R-dx, R -y = np.arange(-R,R+dy,dy,dtype='float32') # y = -R, -R+dx,..., -dx, 0, dx,..., R-dx, R +x = np.arange(-R+dx/2.,R+dx/2.,dx,dtype='float32') # x = -R+dx/2, -R+3*dx/2,..., -dx/2, dx/2.,..., R-dx/2 +y = np.arange(-R+dy/2.,R+dy/2.,dy,dtype='float32') # y = -R+dy/2, -R+3*dy/2,..., -dy/2, dy/2.,..., R-dy/2 #print('x=',x[-4::]) #print('y=',y[-4::]) -#sys.exit('coucou') x1[:] = x[:] y1[:] = y[:] - -x0[:] = x[:-1] + dx/2. # x = -R+dx/2,..., -dx/2, dx/2,..., R-dx/2 -y0[:] = y[:-1] + dy/2. # x = -R+dy/2,..., -dy/2, dy/2,..., R-dy/2 +x0[:] = x[:-1] + dx/2. # x = -R+dx,..., -dx, 0, dx,..., R-dx +y0[:] = y[:-1] + dy/2. # x = -R+dy,..., -dy, 0, dy,..., R-dy # Set SMB acab[:,:,:] = accum @@ -349,16 +348,12 @@ def createFileFromSource(src_file, trg_file): thk[:,:,:] = 0. calving_mask[:,:,:] = 1 -for i in range(1,nx-1): - for j in range(1,ny-1): - # Find the distance r from the origin to the most distant cell vertex. - # If r < rcalve, then put ice in the cell and set calving_mask = 0; if not, then set H = 0 and calving_mask = 1. - rsw = np.sqrt(x0[i-1]**2 + y0[j-1]**2) - rnw = np.sqrt(x0[i-1]**2 + y0[j]**2) - rne = np.sqrt(x0[i]**2 + y0[j]**2) - rse = np.sqrt(x0[i]**2 + y0[j-1]**2) - rmax = max(rsw,rnw,rne,rse) # distance from the origin to the most distant vertex - if rmax <= rcalve: +for i in range(0,nx): + for j in range(0,ny): + # Find the distance d from the origin to the cell center. + # If d < rcalve, then put ice in the cell and set calving_mask = 0; if not, then set H = 0 and calving_mask = 1. + d = np.sqrt(x1[i]**2 + y1[j]**2) + if d <= rcalve: thk[:,j,i] = initThickness calving_mask[:,j,i] = 0 @@ -394,7 +389,7 @@ def createFileFromSource(src_file, trg_file): ncfile['topg'][:,j,i] = computeBedThule(ncfile['x1'][i], ncfile['y1'][j], R, Bc_thule, Bl_thule, Ba_thule) -# Close the file +# Close the file ncfile.close() print('Experiments:', experiments) @@ -490,15 +485,10 @@ def createFileFromSource(src_file, trg_file): cf_advance_retreat_amplitude = -300. cf_advance_retreat_period = 1000. elif expt == 'Experiment4': - # for retreat500 + # for retreat_500 (will be changed below for advance_1000) cf_advance_retreat_amplitude = -750. cf_advance_retreat_period = 1000. - # for advance_500_1000 - cf_advance_retreat_amplitude = 5000. - cf_advance_retreat_period = 0. - #TODO: Modify for Experiment 5 - # Set other parameters specific to certain experiments # TODO: Do we need to read in the input temperature? Or do we always want temp_init = 1? @@ -514,6 +504,12 @@ def createFileFromSource(src_file, trg_file): calvingMinthck = 325. config.set('parameters', 'calving_minthck', str(calvingMinthck)) + # Set the calvingMIP domain (circular or Thule) + if expt == 'SpinupCircular' or expt == 'Experiment1' or expt == 'Experiment2': + config.set('ho_options', 'which_ho_calvingmip_domain', '1') + elif expt == 'SpinupCircular' or expt == 'Experiment3' or expt == 'Experiment4' or expt == 'Experiment5': + config.set('ho_options', 'which_ho_calvingmip_domain', '2') + # Set the start and end times config.set('time', 'tstart', str(tstart)) config.set('time', 'tend', str(tend)) @@ -544,7 +540,7 @@ def createFileFromSource(src_file, trg_file): # Set the output filename in the section '[CF output]'. outputfile = 'calvingMIP.' + expt + '.out.nc' config.set('CF output', 'name', outputfile) - config.set('CF output', 'freq', str(outputfreq)) + config.set('CF output', 'frequency', str(outputfreq)) # print('Output file:', outputfile) # Specify additional output files for CalvingMIP experiments. diff --git a/tests/calvingMIP/calvingMIP.config.template b/tests/calvingMIP/calvingMIP.config.template index 48969bff..acc99490 100644 --- a/tests/calvingMIP/calvingMIP.config.template +++ b/tests/calvingMIP/calvingMIP.config.template @@ -1,10 +1,11 @@ [grid] upn = 5 -ewn = 321 -nsn = 321 +ewn = 320 +nsn = 320 dew = 5000 dns = 5000 global_bc = 1 # 1 = outflow BCs + [time] tstart = 0. tend = 1000.0 @@ -75,13 +76,13 @@ time = 1 [CF output] variables = thk topg usurf uvel vvel beta_internal effecpress f_ground f_flotation floating_mask grounded_mask usfc vsfc ubas vbas uvel_mean vvel_mean stagthk ivol imass_above_flotation iareag calving_thck calving_mask -frequency = 1000 +frequency = 1000.0 name = calvingMIP.out.nc [CF restart] variables = restart xtype = double -frequency = 1000 +frequency = 1000.0 name = calvingMIP.restart.nc write_init = F From a52d3a5fc3ff54326c088c9a1921425534ff2bb7 Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Thu, 14 Nov 2024 16:41:18 -0700 Subject: [PATCH 3/4] Extended the CalvingMIP grid This commit extends the Circular and Thule grids by one cell in each direction. E.g., for a 5km run, where there were previously 160 cells in each direction, there are now 162. This extension makes it easier to interpolate scalars to the CalvingMIP results grid, which corresponds to the CISM (x0,y0) grid. I re-ran the experiments, which should provide our final submitted results. Compare to the old grid, answers change at the machine roundoff level, presumably because the grid decomposition is different. --- tests/calvingMIP/calvingMIP.Setup.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/tests/calvingMIP/calvingMIP.Setup.py b/tests/calvingMIP/calvingMIP.Setup.py index f2aac7b7..51fd1eaa 100644 --- a/tests/calvingMIP/calvingMIP.Setup.py +++ b/tests/calvingMIP/calvingMIP.Setup.py @@ -200,10 +200,12 @@ def createFileFromSource(src_file, trg_file): print( 'Number of vertical levels =', nz) # Set number of grid cells in each direction. -# E.g., if R = 800 km with dx = dy = 10 km, we have 160 cells in each direction -# Note: The origin lies at a cell corner -nx = 2*int(R/dx) -ny = 2*int(R/dy) +# E.g., if R = 800 km with dx = dy = 10 km, we have 162 cells in each direction. +# The origin lies at a cell corner. +# The outer row of cells lies outside the CalvingMIP domain but is +# useful for interpolating data to the CalvingMIP grid. +nx = 2*int(R/dx) + 2 +ny = 2*int(R/dy) + 2 print('grid dimension in x:',nx) print('grid dimension in y:',ny) @@ -329,17 +331,20 @@ def createFileFromSource(src_file, trg_file): # Compute x and y on each grid. # The origin (x = y = 0) is placed at the center of the domain. - -x = np.arange(-R+dx/2.,R+dx/2.,dx,dtype='float32') # x = -R+dx/2, -R+3*dx/2,..., -dx/2, dx/2.,..., R-dx/2 -y = np.arange(-R+dy/2.,R+dy/2.,dy,dtype='float32') # y = -R+dy/2, -R+3*dy/2,..., -dy/2, dy/2.,..., R-dy/2 +# The (x0,y0) grid extends a distance R along the x and y axes. +# The (x1,y1) grid extends by one additional grid cell in each direction. +# The extension makes it easier to interpolate scalars onto the CalvingMIP grid, +# which corresponds to (x0,y0). +x = np.arange(-R-dx/2.,R+3.*dx/2.,dx,dtype='float32') # x = -R-dx/2, -R+3*dx/2,..., -dx/2, dx/2.,..., R+dx/2 +y = np.arange(-R-dy/2.,R+3.*dy/2.,dy,dtype='float32') # y = -R-dy/2, -R+3*dy/2,..., -dy/2, dy/2.,..., R+dy/2 #print('x=',x[-4::]) #print('y=',y[-4::]) x1[:] = x[:] y1[:] = y[:] -x0[:] = x[:-1] + dx/2. # x = -R+dx,..., -dx, 0, dx,..., R-dx -y0[:] = y[:-1] + dy/2. # x = -R+dy,..., -dy, 0, dy,..., R-dy +x0[:] = x[:-1] + dx/2. # x = -R,..., -dx, 0, dx,..., R +y0[:] = y[:-1] + dy/2. # x = -R,..., -dy, 0, dy,..., R # Set SMB acab[:,:,:] = accum From 02f23cfa4c0fbd58b455429c177d8a02e5d9e9db Mon Sep 17 00:00:00 2001 From: William Lipscomb Date: Mon, 18 Nov 2024 17:49:41 -0700 Subject: [PATCH 4/4] Changed the thck_effective calculation at the calving front This is the version used for the CISM CalvingMIP submission in Nov. 2024. This commit changes the way thck_effective (the effective ice thickness) is computed in calving front cells. For any CF cell with an interior neighbor, the calculation is unchanged; thck_effective is set to the maximum thickness of the cell's interior edge neighbors. For CF cells without interior edge neighbors (i.e., cells whose only edge neighbors are other CF cells), thck_effective is now set to the maximum thickness of the CF neighbors. Previously, we set thck_effective = thck, giving values that are too low. I also fixed a minor bug in the way the CF location is computed along the diagonal axes for the circular domain. I set dthck_dx_cf = 0.0005 in the config template file for CalvingMIP. This value gives slightly more accurate results than the default value of 0. Setting dthck_dx_cf = 0 gives a little too much retreat. If desired, we could modify the python setup script to control this parameter at the command line. We ran the final CalvingMIP runs with dthck_dx_cf = 0.0005. --- libglissade/glissade_calving.F90 | 46 ++++++++++----------- libglissade/glissade_masks.F90 | 18 +++++--- tests/calvingMIP/calvingMIP.config.template | 2 +- 3 files changed, 37 insertions(+), 29 deletions(-) diff --git a/libglissade/glissade_calving.F90 b/libglissade/glissade_calving.F90 index 2fa9e245..13295397 100644 --- a/libglissade/glissade_calving.F90 +++ b/libglissade/glissade_calving.F90 @@ -704,10 +704,6 @@ subroutine glissade_calve_ice(nx, ny, & call parallel_halo(calving_dthck, parallel) - if (verbose_calving) then - call point_diag(calving_dthck, 'dthck (m)', itest, jtest, rtest, 7, 7) - endif - ! Apply calving_dthck as computed above. call apply_calving_dthck(& @@ -722,6 +718,11 @@ subroutine glissade_calve_ice(nx, ny, & thck, & calving%calving_thck) + if (verbose_calving) then + call point_diag(calving_dthck, 'dthck (m)', itest, jtest, rtest, 7, 7) + call point_diag(thck, 'New thck (m)', itest, jtest, rtest, 7, 7) + endif + !TODO - Add a bug check for negative thicknesses? thck = max(thck, 0.0d0) @@ -757,13 +758,6 @@ subroutine glissade_calve_ice(nx, ny, & full_mask = full_mask, & effective_areafrac = calving%effective_areafrac) - if (verbose_calving) then - call point_diag(thck, 'After calving dthck: thck (m)', itest, jtest, rtest, 7, 7) -! call point_diag(calving%thck_effective, 'thck_effective (m)', itest, jtest, rtest, 7, 7) -! call point_diag(calving%effective_areafrac, 'effective_areafrac', itest, jtest, rtest, 7, 7) -! call point_diag(thck/calving%effective_areafrac, 'CF thickness', itest, jtest, rtest, 7, 7) - endif - if (which_calving == CF_ADVANCE_RETREAT_RATE) then ! Compute some CalvingMIP diagnostics. @@ -811,6 +805,12 @@ subroutine glissade_calve_ice(nx, ny, & print*, 'Quadrant area (km^2)=', total_ice_area/4.0d6 endif + if (verbose_calving) then + call point_diag(calving%thck_effective, 'New thck_effective (m)', itest, jtest, rtest, 7, 7) + call point_diag(calving%effective_areafrac, 'New effective_areafrac', itest, jtest, rtest, 7, 7) +!! call point_diag(thck/calving%effective_areafrac, ' CF thickness', itest, jtest, rtest, 7, 7) + endif + else ! other calving schemes ! Where thck > thck_effective, allow the CF to advance by distributing ice downstream. @@ -2329,8 +2329,8 @@ subroutine glissade_remove_icebergs(& if (verbose_calving) then call point_diag(thck, 'Remove icebergs, thck (m)', itest, jtest, rtest, 7, 7) - call point_diag(ice_mask, 'ice_mask', itest, jtest, rtest, 7, 7) - call point_diag(f_ground_cell, 'f_ground_cell', itest, jtest, rtest, 7, 7) +!! call point_diag(ice_mask, 'ice_mask', itest, jtest, rtest, 7, 7) +!! call point_diag(f_ground_cell, 'f_ground_cell', itest, jtest, rtest, 7, 7) endif ! Initialize iceberg removal @@ -2449,8 +2449,8 @@ subroutine glissade_remove_icebergs(& global_count = parallel_reduce_sum(local_count) if (global_count == global_count_save) then - if (verbose_calving .and. main_task) & - write(6,*) 'Fill converged: iter, global_count =', iter, global_count +!! if (verbose_calving .and. main_task) & +!! write(6,*) 'Fill converged: iter, global_count =', iter, global_count exit else !! if (verbose_calving .and. main_task) & @@ -3107,8 +3107,8 @@ subroutine locate_calving_front_circular(& cf_location(1,axis) = x0(i) + (areafrac_ne - 0.5d0)/areafrac_ne * (0.5d0*dx) cf_location(2,axis) = y0(j) + (areafrac_ne - 0.5d0)/areafrac_ne * (0.5d0*dy) elseif (areafrac_sw < 0.5d0) then ! CF in cell (i-1,j-1) - cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_sw)/(areafrac(i-1,j-1) - areafrac_sw) * (0.5d0*dx) - cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_sw)/(areafrac(i-1,j-1) - areafrac_sw) * (0.5d0*dy) + cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_sw)/(1.0d0 - areafrac_sw) * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_sw)/(1.0d0 - areafrac_sw) * (0.5d0*dy) else ! CF in cell (i,j) if (areafrac(i,j) >= 0.5d0) then ! CF in upper right of cell cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_ne) * (0.5d0*dx) @@ -3169,8 +3169,8 @@ subroutine locate_calving_front_circular(& cf_location(1,axis) = x0(i) + (areafrac_se - 0.5d0)/areafrac_se * (0.5d0*dx) cf_location(2,axis) = y0(j-1) - (areafrac_se - 0.5d0)/areafrac_se * (0.5d0*dy) elseif (areafrac_nw < 0.5d0) then ! CF in cell (i-1,j+1) - cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_nw)/(areafrac(i-1,j-1) - areafrac_nw) * (0.5d0*dx) - cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_nw)/(areafrac(i-1,j-1) - areafrac_nw) * (0.5d0*dy) + cf_location(1,axis) = x0(i-1) - (0.5d0 - areafrac_nw)/(1.0d0 - areafrac_nw) * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_nw)/(1.0d0 - areafrac_nw) * (0.5d0*dy) else ! CF in cell (i,j) if (areafrac(i,j) >= 0.5d0) then ! CF in lower right of cell cf_location(1,axis) = x1(i) + (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_se) * (0.5d0*dx) @@ -3231,8 +3231,8 @@ subroutine locate_calving_front_circular(& cf_location(1,axis) = x0(i-1) - (areafrac_sw - 0.5d0)/areafrac_sw * (0.5d0*dx) cf_location(2,axis) = y0(j-1) - (areafrac_sw - 0.5d0)/areafrac_sw * (0.5d0*dy) elseif (areafrac_ne < 0.5d0) then ! CF in cell (i+1,j+1) - cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_ne)/(areafrac(i-1,j-1) - areafrac_ne) * (0.5d0*dx) - cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_ne)/(areafrac(i-1,j-1) - areafrac_ne) * (0.5d0*dy) + cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_ne)/(1.0d0 - areafrac_ne) * (0.5d0*dx) + cf_location(2,axis) = y0(j) + (0.5d0 - areafrac_ne)/(1.0d0 - areafrac_ne) * (0.5d0*dy) else ! CF in cell (i,j) if (areafrac(i,j) >= 0.5d0) then ! CF in lower left of cell cf_location(1,axis) = x1(i) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_sw) * (0.5d0*dx) @@ -3293,8 +3293,8 @@ subroutine locate_calving_front_circular(& cf_location(1,axis) = x0(i-1) - (areafrac_nw - 0.5d0)/areafrac_nw * (0.5d0*dx) cf_location(2,axis) = y0(j) + (areafrac_nw - 0.5d0)/areafrac_nw * (0.5d0*dy) elseif (areafrac_se < 0.5d0) then ! CF in cell (i+1,j-1) - cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_se)/(areafrac(i+1,j-1) - areafrac_se) * (0.5d0*dx) - cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_se)/(areafrac(i+1,j-1) - areafrac_se) * (0.5d0*dy) + cf_location(1,axis) = x0(i) + (0.5d0 - areafrac_se)/(1.0d0 - areafrac_se) * (0.5d0*dx) + cf_location(2,axis) = y0(j-1) - (0.5d0 - areafrac_se)/(1.0d0 - areafrac_se) * (0.5d0*dy) else ! CF in cell (i,j) if (areafrac(i,j) >= 0.5d0) then ! CF in upper left of cell cf_location(1,axis) = x1(i) - (areafrac(i,j) - 0.5d0)/(areafrac(i,j) - areafrac_nw) * (0.5d0*dx) diff --git a/libglissade/glissade_masks.F90 b/libglissade/glissade_masks.F90 index d9990829..0c3f5f3a 100644 --- a/libglissade/glissade_masks.F90 +++ b/libglissade/glissade_masks.F90 @@ -421,12 +421,20 @@ subroutine glissade_calving_front_mask(& else full_mask(i,j) = 1 endif ! dthck_dx > dthck_dx_cf - else ! no interior neighbors; mark as a partial cell + else ! no interior neighbors + ! Mark as a partial cell, and compute thck_effective from a CF neighbor partial_cf_mask(i,j) = 1 - call parallel_globalindex(i, j, ig, jg, parallel) -!! if (this_rank == rtest) then -!! write(6,*) 'No interior neighbor:', ig, jg, thck(i,j) -!! endif + max_neighbor_thck = max(& + calving_front_mask(i-1,j)*thck(i-1,j), calving_front_mask(i+1,j)*thck(i+1,j), & + calving_front_mask(i,j-1)*thck(i,j-1), calving_front_mask(i,j+1)*thck(i,j+1)) + distance = sqrt(dx*dy) + dthck_dx = (max_neighbor_thck - thck(i,j)) / distance + if (dthck_dx > dthck_dx_cf) then + thck_effective(i,j) = max_neighbor_thck - dthck_dx_cf*distance + endif +!! call parallel_globalindex(i, j, ig, jg, parallel) +!! write(6,*) 'No interior neighbor:', ig, jg, thck(i,j) +!! write(6,*) ' New H_eff:', thck_effective(i,j) endif ! max_neighbor_thck > 0 else ! not a CF cell; thck_effective = thck diff --git a/tests/calvingMIP/calvingMIP.config.template b/tests/calvingMIP/calvingMIP.config.template index acc99490..27f8a6a9 100644 --- a/tests/calvingMIP/calvingMIP.config.template +++ b/tests/calvingMIP/calvingMIP.config.template @@ -62,7 +62,7 @@ powerlaw_m = 3. # exponent for the power law stress geothermal = 0. beta_grounded_min = 1. thck_effective_min = 50. -dthck_dx_cf = 0. +dthck_dx_cf = 0.0005 # surface slope at the CF cf_advance_retreat_amplitude = 0. cf_advance_retreat_period = 0.