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..13295397 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 @@ -661,12 +704,7 @@ 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. - ! Adjust the thinning based on the estimated calving-front length in the cell. call apply_calving_dthck(& nx, ny, & @@ -680,15 +718,17 @@ 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) 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 +755,105 @@ 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) - ! 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 + + 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. + ! 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 @@ -2264,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 @@ -2384,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) & @@ -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)/(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) + 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)/(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) + 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)/(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) + 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)/(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) + 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..0c3f5f3a 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 :: & @@ -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/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..51fd1eaa 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,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 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 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) @@ -236,10 +239,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) @@ -328,19 +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,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 +# 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::]) -#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, 0, dx,..., R +y0[:] = y[:-1] + dy/2. # x = -R,..., -dy, 0, dy,..., R # Set SMB acab[:,:,:] = accum @@ -349,16 +353,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 +394,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 +490,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 +509,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 +545,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..27f8a6a9 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 @@ -61,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. @@ -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