diff --git a/src/icebergs.F90 b/src/icebergs.F90
index b10b4f8..a50b58f 100644
--- a/src/icebergs.F90
+++ b/src/icebergs.F90
@@ -58,7 +58,7 @@ module ice_bergs
 use ice_bergs_io,        only: ice_bergs_io_init, write_restart_bergs, write_trajectory, write_bond_trajectory
 use ice_bergs_io,        only: read_restart_bergs, read_restart_calving
 use ice_bergs_io,        only: read_restart_bonds
-use ice_bergs_io,        only: read_ocean_depth
+use ice_bergs_io,        only: read_ocean_depth, read_ice_sheet_basins
 
 implicit none ; private
 
@@ -117,8 +117,8 @@ subroutine icebergs_init(bergs, &
   logical, intent(in), optional :: fractional_area !< If true, ice_area contains cell area as fraction of entire spherical surface
   ! Local variables
   type(icebergs_gridded), pointer :: grd => null()
-  integer :: nbonds
-  logical :: check_bond_quality
+  integer :: nbonds, nbasins, id_class
+  logical :: check_bond_quality, lerr
   integer :: stdlogunit, stderrunit
 
   ! Get the stderr and stdlog unit numbers
@@ -131,6 +131,8 @@ subroutine icebergs_init(bergs, &
              dt, Time, ice_lon, ice_lat, ice_wet, ice_dx, ice_dy, ice_area, &
              cos_rot, sin_rot, ocean_depth=ocean_depth, maskmap=maskmap, fractional_area=fractional_area)
 
+  grd=>bergs%grd
+
   call unit_testing(bergs)
 
   call mpp_clock_begin(bergs%clock_ior)
@@ -149,6 +151,18 @@ subroutine icebergs_init(bergs, &
 
   !Reading ocean depth from a file
   if (bergs%read_ocean_depth_from_file) call read_ocean_depth(bergs%grd)
+  ! Reading the ice-sheet basins of origin for the bergs
+  if (bergs%use_berg_origin_basins) then
+    call read_ice_sheet_basins(bergs%grd)
+    id_class=register_static_field('icebergs', 'ice_sheet_basins_static', axes, &
+                 'Static ice-sheet basins of origin for icebergs', 'none')
+    if (id_class>0) lerr=send_data(id_class, grd%ice_sheet_basins(grd%isc:grd%iec,grd%jsc:grd%jec))
+  else
+    bergs%nbasins=1
+    grd%ice_sheet_basins(:,:)=0.0
+  endif
+  allocate( grd%melt_by_ice_sheet_basin(grd%isd:grd%ied, grd%jsd:grd%jed, bergs%nbasins) )
+  grd%melt_by_ice_sheet_basin(:,:,:)=0.
 
   if (bergs%iceberg_bonds_on) then
     if (bergs%manually_initialize_bonds) then
@@ -3132,6 +3146,13 @@ subroutine thermodynamics(bergs)
           grd%melt_by_class(i,j,k)=grd%melt_by_class(i,j,k)+melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s
         endif
 
+        if (bergs%use_berg_origin_basins) then
+          if (this%basin>0) then
+            grd%melt_by_ice_sheet_basin(i,j,this%basin)=grd%melt_by_ice_sheet_basin(i,j,this%basin)+&
+                                                        melt/grd%area(i,j)*this%mass_scaling ! kg/m2/s
+          endif
+        endif
+
         melt=melt*this%heat_density ! kg/s x J/kg = J/s
         grd%calving_hflx(i,j)=grd%calving_hflx(i,j)+melt/grd%area(i,j)*this%mass_scaling ! W/m2
         bergs%net_heat_to_ocean=bergs%net_heat_to_ocean+melt*this%mass_scaling*bergs%dt ! J
@@ -5152,6 +5173,7 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh,
   grd%mass(:,:)=0.
   grd%virtual_area(:,:)=0.
   grd%melt_by_class(:,:,:)=0.
+  grd%melt_by_ice_sheet_basin(:,:,:) = 0.
   grd%melt_buoy_fl(:,:)=0.
   grd%melt_eros_fl(:,:)=0.
   grd%melt_conv_fl(:,:)=0.
@@ -5613,6 +5635,10 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh,
     lerr=send_data(grd%id_fay, tauya(:,:), Time)
   if (grd%id_melt_by_class>0) &
     lerr=send_data(grd%id_melt_by_class, grd%melt_by_class(grd%isc:grd%iec,grd%jsc:grd%jec,:), Time)
+  if (grd%id_melt_by_ice_sheet_basin>0) &
+    lerr=send_data(grd%id_melt_by_ice_sheet_basin, grd%melt_by_ice_sheet_basin(grd%isc:grd%iec,grd%jsc:grd%jec,:), Time)
+  if (grd%id_ice_sheet_basins>0) &
+    lerr=send_data(grd%id_ice_sheet_basins, grd%ice_sheet_basins(grd%isc:grd%iec,grd%jsc:grd%jec), Time)
   if (grd%id_melt_buoy_fl>0) &
     lerr=send_data(grd%id_melt_buoy_fl, grd%melt_buoy_fl(grd%isc:grd%iec,grd%jsc:grd%jec), Time)
   if (grd%id_melt_eros_fl>0) &
@@ -6362,6 +6388,13 @@ subroutine calve_icebergs(bergs)
             newberg%ang_vel=0.; newberg%ang_accel=0.; newberg%rot=0.
           endif
 
+          if (bergs%use_berg_origin_basins) then
+            if (.not. allocations_done) then
+              if (.not. allocated(newberg%basin)) allocate(newberg%basin)
+            endif
+            newberg%basin=int(grd%ice_sheet_basins(i,j))
+          endif
+
           if (.not. bergs%old_interp_flds_order) then
             !interpolate gridded variables to new iceberg
             if (grd%tidal_drift>0.) then
@@ -6572,6 +6605,11 @@ subroutine calve_fl_icebergs(bergs,pberg,k,l_b,fl_disp_x,fl_disp_y,berg_from_bit
     cberg%ang_vel=0.; cberg%ang_accel=0.; cberg%rot=0.
   endif
 
+  if (bergs%use_berg_origin_basins) then
+    allocate(cberg%basin)
+    cberg%basin=pberg%basin
+  endif
+
   call add_new_berg_to_list(bergs%list(cberg%ine,cberg%jne)%first, cberg)
 end subroutine calve_fl_icebergs
 
@@ -8230,6 +8268,8 @@ subroutine icebergs_end(bergs)
   deallocate(bergs%grd%cn)
   deallocate(bergs%grd%hi)
   deallocate(bergs%grd%melt_by_class)
+  deallocate(bergs%grd%melt_by_ice_sheet_basin)
+  deallocate(bergs%grd%ice_sheet_basins)
   deallocate(bergs%grd%melt_buoy_fl)
   deallocate(bergs%grd%melt_eros_fl)
   deallocate(bergs%grd%melt_conv_fl)
diff --git a/src/icebergs_fms2io.F90 b/src/icebergs_fms2io.F90
index ddc97f2..b41c03b 100644
--- a/src/icebergs_fms2io.F90
+++ b/src/icebergs_fms2io.F90
@@ -44,12 +44,12 @@ module ice_bergs_fms2io
 use ice_bergs_framework, only: force_all_pes_traj
 use ice_bergs_framework, only: check_for_duplicates_in_parallel
 use ice_bergs_framework, only: split_id, id_from_2_ints, generate_id
-! for MTS/DEM/fracture/footloose:
+! for MTS/DEM/fracture/footloose/basins:
 use ice_bergs_framework, only: mts,save_bond_traj
 use ice_bergs_framework, only: push_bond_posn, append_bond_posn
 use ice_bergs_framework, only: pack_bond_traj_into_buffer2,unpack_bond_traj_from_buffer2
 use ice_bergs_framework, only: dem, iceberg_bonds_on
-use ice_bergs_framework, only: footloose
+use ice_bergs_framework, only: footloose, use_berg_origin_basins
 
 
 implicit none ; private
@@ -59,7 +59,7 @@ module ice_bergs_fms2io
 public ice_bergs_io_init
 public read_restart_bergs, write_restart_bergs, write_trajectory, write_bond_trajectory
 public read_restart_calving, read_restart_bonds
-public read_ocean_depth
+public read_ocean_depth, read_ice_sheet_basins
 
 !Local Vars
 integer, parameter :: file_format_major_version=0
@@ -187,7 +187,8 @@ subroutine write_restart_bergs(bergs, time_stamp)
                                       first_berg_ine,         &
                                       other_berg_jne,         &
                                       other_berg_ine,         &
-                                      broken
+                                      broken,           &
+                                      basin
 
 
 integer :: grdi, grdj
@@ -258,6 +259,7 @@ subroutine write_restart_bergs(bergs, time_stamp)
      allocate(ang_accel(nbergs))
      allocate(rot(nbergs))
    endif
+   if (use_berg_origin_basins) allocate(basin(nbergs))
 
   i = 0
   do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
@@ -296,6 +298,7 @@ subroutine write_restart_bergs(bergs, time_stamp)
         ang_accel(i) = this%ang_accel
         rot(i) = this%rot
       endif
+      if (use_berg_origin_basins) basin(i) = this%basin
       this=>this%next
     enddo
   enddo ; enddo
@@ -393,6 +396,10 @@ subroutine write_restart_bergs(bergs, time_stamp)
                                      dim_names_1d,longname='dem accumulated rotation',units='rad')
   endif
 
+  if (use_berg_origin_basins) then
+    call register_restart_field_wrap(fileobj,'basin',basin,&
+                                     dim_names_1d,longname='ice-sheet basin of origin',units='none')
+  endif
   !Checking if any icebergs are static in order to decide whether to save static_berg
   n_static_bergs = 0
   do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec
@@ -457,6 +464,8 @@ subroutine write_restart_bergs(bergs, time_stamp)
              rot)
   endif
 
+  if (use_berg_origin_basins) deallocate(basin)
+
   deallocate(           &
              ine,       &
              jne,       &
@@ -711,7 +720,8 @@ subroutine read_restart_bergs(bergs,Time)
                                       iceberg_num,&
                                       id_cnt,     &
                                       id_ij,      &
-                                      start_year
+                                      start_year, &
+                                      basin
 
 type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2_io fileobj
 character(len=1), dimension(1) :: dim_names_1d
@@ -809,6 +819,10 @@ subroutine read_restart_bergs(bergs,Time)
        allocate(ang_accel(nbergs_in_file))
        allocate(rot(nbergs_in_file))
      endif
+     if (use_berg_origin_basins) then
+       allocate(localberg%basin)
+       allocate(basin(nbergs_in_file))
+     endif
 
      call register_restart_field(fileobj,'lon',lon,dim_names_1d)
      call register_restart_field(fileobj,'lat',lat,dim_names_1d)
@@ -858,6 +872,11 @@ subroutine read_restart_bergs(bergs,Time)
        call register_restart_field(fileobj,'ang_accel',ang_accel,dim_names_1d,is_optional=.true.)
        call register_restart_field(fileobj,'rot'      ,rot      ,dim_names_1d,is_optional=.true.)
      endif
+
+     if (use_berg_origin_basins) then
+       basin = 0
+       call register_restart_field(fileobj,'basin'    ,basin    ,dim_names_1d,is_optional=.true.)
+     endif
       call read_restart(fileobj)
       call close_file(fileobj)
   elseif (bergs%require_restart) then
@@ -943,6 +962,10 @@ subroutine read_restart_bergs(bergs,Time)
        localberg%rot      =rot(k)
      endif
 
+     if (use_berg_origin_basins) then
+       localberg%basin    =basin(k)
+     endif
+
       if (really_debug) lres=is_point_in_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, explain=.true.)
       lres=pos_within_cell(grd, localberg%lon, localberg%lat, localberg%ine, localberg%jne, localberg%xi, localberg%yj)
       !call add_new_berg_to_list(bergs%first, localberg)
@@ -1001,6 +1024,7 @@ subroutine read_restart_bergs(bergs,Time)
                ang_accel,    &
                rot)
     endif
+    if (use_berg_origin_basins) deallocate(basin)
 
     if (replace_iceberg_num) then
       deallocate(iceberg_num)
@@ -1076,6 +1100,7 @@ subroutine generate_bergs(bergs,Time)
     allocate(localberg%ang_accel)
     allocate(localberg%rot)
   endif
+  if (use_berg_origin_basins) allocate(localberg%basin)
 
   do j=grd%jsc,grd%jec; do i=grd%isc,grd%iec
     if (grd%msk(i,j)>0. .and. abs(grd%latc(i,j))>80.0) then
@@ -1125,6 +1150,9 @@ subroutine generate_bergs(bergs,Time)
         localberg%ang_accel=0.
         localberg%rot=0.
       endif
+      if (use_berg_origin_basins) then
+        localberg%basin=0
+      endif
 
       !Berg A
       call loc_set_berg_pos(grd, 0.9, 0.5, 1., 0., localberg)
@@ -1603,7 +1631,7 @@ subroutine read_ocean_depth(grd)
 ! Local variables
 character(len=37) :: filename
 type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2_io fileobj
-  ! Read stored ice
+  ! Read depth
   filename=trim(restart_input_dir)//'topog.nc'
   if (open_file(fileobj, filename, "read", grd%domain)) then
     if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') &
@@ -1627,6 +1655,34 @@ subroutine read_ocean_depth(grd)
   !call grd_chksum2(bergs%grd, bergs%grd%ocean_depth, 'read_ocean_depth, ocean_depth')
 end subroutine read_ocean_depth
 
+!> Read ice-sheet basins from file
+subroutine read_ice_sheet_basins(grd)
+! Arguments
+type(icebergs_gridded), pointer :: grd !< Container for gridded fields
+! Local variables
+character(len=37) :: filename, actual_filename
+type(FmsNetcdfDomainFile_t) :: fileobj !< Fms2_io fileobj
+  ! Read sub_basin
+  filename=trim(restart_input_dir)//'ice_sheet_basins.nc'
+  if (open_file(fileobj, filename, "read", grd%domain)) then
+    if (mpp_pe().eq.mpp_root_pe()) write(*,'(3a)') &
+     'KID, read_ice_sheet_basins: reading ',actual_filename, filename
+    call register_axis_wrapper(fileobj)
+    if (variable_exists(fileobj, "sub_basin")) then
+      if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(a)') &
+       'KID, read_ice_sheet_basins: reading sub_basins from ice-shelf basins file.'
+      call read_data(fileobj, 'sub_basin', grd%ice_sheet_basins)
+    else
+      call error_mesg('KID, read_ice_sheet_basins', &
+        'variable sub_basin ice_sheet_basins.nc not found in ice_sheet_basins.nc!', FATAL)
+    endif
+    call close_file(fileobj)
+  else
+    call error_mesg('KID, read_ice_sheet_basins', 'ice_sheet_basins.nc not found!', FATAL)
+  endif
+
+end subroutine read_ice_sheet_basins
+
 !> Write a trajectory-based diagnostics file
 subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name)
 ! Arguments
@@ -1642,7 +1698,7 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name
 integer :: cnid, hiid, hsid
 integer :: mid, smid, did, wid, lid, mbid, mflbid, mflbbid, hdid, nbid, odid, flkid
 integer :: axnid,aynid,bxnid,bynid,axnfid,aynfid,bxnfid,bynfid, msid
-integer :: avid, aaid, rid
+integer :: avid, aaid, rid, baid
 character(len=70) :: filename
 character(len=7) :: pe_name
 type(xyt), pointer :: this, next
@@ -1819,6 +1875,9 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name
           rid  = inq_varid(ncid, 'rot')
         endif
 
+        if (use_berg_origin_basins) then
+          baid = inq_varid(ncid, 'basin')
+        endif
       endif
     else
       ! Dimensions
@@ -1889,6 +1948,10 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name
           aaid = def_var(ncid, 'ang_accel', NF_DOUBLE, i_dim)
           rid  = def_var(ncid, 'rot', NF_DOUBLE, i_dim)
         endif
+
+        if (use_berg_origin_basins) then
+          baid = def_var(ncid, 'basin', NF_INT, i_dim)
+        endif
       endif
 
       ! Attributes
@@ -2006,6 +2069,10 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name
           call put_att(ncid, rid, 'long_name', 'accumulated rotation')
           call put_att(ncid, rid, 'units', 'rad')
         endif
+        if (use_berg_origin_basins) then
+          call put_att(ncid, baid, 'long_name', 'ice-sheet basin of origin')
+          call put_att(ncid, baid, 'units', 'none')
+        endif
       endif
     endif
 
@@ -2087,6 +2154,9 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name
           call put_double(ncid, rid,  i, this%rot)
         endif
 
+        if (use_berg_origin_basins) then
+          call put_int(ncid, baid, i, this%basin)
+        endif
       endif
       next=>this%next
       deallocate(this)
diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90
index c327135..ed7f1e0 100644
--- a/src/icebergs_framework.F90
+++ b/src/icebergs_framework.F90
@@ -62,6 +62,7 @@ module ice_bergs_framework
 logical :: skip_first_outer_mts_step=.false.
 logical :: no_frac_first_ts=.false.
 logical :: footloose=.false. !< Turn footloose calving on/off
+logical :: use_berg_origin_basins=.false. !< If T, save berg melt associated with the each ice-sheet basin of origin for the bergs
 
 !Public params !Niki: write a subroutine to expose these
 public nclasses,buffer_width,buffer_width_traj,buffer_width_bond_traj
@@ -72,7 +73,7 @@ module ice_bergs_framework
 public dem, save_bond_forces, orig_dem_moment_of_inertia
 public short_step_mts_grounding, radius_based_drag
 public A68_test, A68_xdisp, A68_ydisp
-public footloose
+public footloose, use_berg_origin_basins
 
 !Public types
 public icebergs_gridded, xyt, iceberg, icebergs, buffer, bond, bond_xyt
@@ -145,6 +146,8 @@ module ice_bergs_framework
   real, dimension(:,:), pointer :: cos=>null() !< Cosine from rotation matrix to lat-lon coords
   real, dimension(:,:), pointer :: sin=>null() !< Sine from rotation matrix to lat-lon coords
   real, dimension(:,:), pointer :: ocean_depth=>NULL() !< Depth of ocean (m)
+  real, dimension(:,:), pointer :: ice_sheet_basins=>NULL() !< Ice-sheet basins of origin for bergs (e.g. IMBIE basins)
+  real, dimension(:,:,:), pointer :: melt_by_ice_sheet_basin=>null() !< Total melt rate for bergs from each ice-sheet basin (kg/s/m^2)
   real, dimension(:,:), pointer :: uo=>null() !< Ocean zonal flow (m/s)
   real, dimension(:,:), pointer :: vo=>null() !< Ocean meridional flow (m/s)
   real, dimension(:,:), pointer :: ui=>null() !< Ice zonal flow (m/s)
@@ -219,7 +222,7 @@ module ice_bergs_framework
   integer :: id_count=-1, id_chksum=-1, id_u_iceberg=-1, id_v_iceberg=-1, id_sss=-1, id_ustar_iceberg
   integer :: id_spread_uvel=-1, id_spread_vvel=-1
   integer :: id_melt_m_per_year=-1
-  integer :: id_ocean_depth=-1
+  integer :: id_ocean_depth=-1, id_ice_sheet_basins=-1, id_melt_by_ice_sheet_basin=-1
   integer :: id_melt_by_class=-1, id_melt_buoy_fl=-1, id_melt_eros_fl=-1, id_melt_conv_fl=-1
   integer :: id_fl_parent_melt=-1, id_fl_child_melt=-1
   !>@}
@@ -284,6 +287,8 @@ module ice_bergs_framework
   real, allocatable :: ang_vel !< Angular velocity
   real, allocatable :: ang_accel !< Angular acceleration
   real, allocatable :: rot !< Accumulated rotation
+  ! If use_berg_origin_basins
+  integer, allocatable :: basin
 end type xyt
 
 !> An iceberg object, used as a link in a linked list
@@ -356,6 +361,8 @@ module ice_bergs_framework
   real, allocatable :: ang_vel !< Angular velocity
   real, allocatable :: ang_accel !< Angular acceleration
   real, allocatable :: rot !< Accumulated rotation
+  ! If use_berg_origin_basins
+  integer, allocatable :: basin
 end type iceberg
 
 !> A bond object connecting two bergs, used as a link in a linked list
@@ -518,6 +525,8 @@ module ice_bergs_framework
   logical :: tang_crit_int_damp_on=.true. !<crit interaction damping for tangential component?
   logical :: use_old_spreading=.true. !< If true, spreads iceberg mass as if the berg is one grid cell wide
   logical :: read_ocean_depth_from_file=.false. !< If true, ocean depth is read from a file.
+  logical :: use_berg_origin_basins=.false. !< If T, save berg melt associated with the each ice-sheet basin of origin for the bergs
+  integer :: nbasins=1 !< Number of ice-sheet basins of origin for the bergs
   integer(kind=8) :: debug_iceberg_with_id = -1 !< If positive, monitors a berg with this id
 
   real :: length_for_manually_initialize_bonds=1000.0 !< If manually initialing bonds, only  bond if dist between particles is < this length -Added by Alex
@@ -603,7 +612,8 @@ module ice_bergs_framework
   real :: ocean_drag_scale=1. !< Scaling factor for the ocean drag coefficients
   ! Footloose calving parameters
   logical :: footloose=.false. !< Turn footloose calving on/off
-  logical :: fl_init_child_xy_by_pe=.false. !< True: position of a new footloose child berg along the parent perimeter is random, yet constant for each PE (old bug). False: fully-random.
+  logical :: fl_init_child_xy_by_pe=.false. !< True: old bug that randomly positions a new footloose child berg along the parent perimeter acoording to PE
+                                            !! (not consistent under different PE layouts). False: Random positioning by element
   real :: fl_youngs=1.e7 !< Young's modulus for footloose calculations (Pa)
   real :: fl_strength=250. !< yield stress for footloose calculations (kPa)
   logical :: displace_fl_bergs=.true. !< footloose berg positions are randomly assigned along edges of parent berg
@@ -777,6 +787,8 @@ subroutine ice_bergs_framework_init(bergs, &
 logical :: read_old_restarts=.false. ! Legacy option that does nothing
 logical :: use_old_spreading=.true. ! If true, spreads iceberg mass as if the berg is one grid cell wide
 logical :: read_ocean_depth_from_file=.false. ! If true, ocean depth is read from a file.
+integer :: nbasins=1 !< Number of ice-sheet basins of origin for the bergs
+real, dimension(:), allocatable :: basin_arr !< array of basin IDs (/ I, I = 1, nbasins) /)
 integer :: mts_sub_steps=-1 ! If -1, the number of mts sub-steps will be automatically determined
 logical :: remove_unused_bergs=.true. ! Remove unneeded bergs after PEs transfers
 real :: contact_distance=0.0 ! For unbonded berg interactions, collision is assumed at max(contact_distance,sum of the 2 bergs radii)
@@ -813,7 +825,8 @@ subroutine ice_bergs_framework_init(bergs, &
 real :: ocean_drag_scale=1. !< Scaling factor for the ocean drag coefficients
 ! Footloose calving parameters
 !logical :: footloose=.false. !< Turn footloose calving on/off
-logical :: fl_init_child_xy_by_pe=.false. !< True: position of a new footloose child berg along the parent perimeter is random, yet constant for each PE (old bug). False: fully-random.
+logical :: fl_init_child_xy_by_pe=.false. !< True: old bug that randomly positions a new footloose child berg along the parent perimeter acoording to PE
+                                          !! (not consistent under different PE layouts). False: Random positioning by element
 real :: fl_youngs=1.e7 !< Young's modulus for footloose calculations (Pa)
 real :: fl_strength=250. !< yield stress for footloose calculations (kPa)
 logical :: displace_fl_bergs=.true. ! footloose berg positions are randomly assigned along edges of parent berg
@@ -840,7 +853,7 @@ subroutine ice_bergs_framework_init(bergs, &
          read_old_restarts, tau_calving, read_ocean_depth_from_file, melt_cutoff,apply_thickness_cutoff_to_gridded_melt,&
          apply_thickness_cutoff_to_bergs_melt, use_mixed_melting, internal_bergs_for_drag, coastal_drift, tidal_drift,&
          mts,ewsame,mts_sub_steps,contact_distance,length_for_manually_initialize_bonds,&
-         manually_initialize_bonds_from_radii,contact_spring_coef,fracture_criterion, &
+         manually_initialize_bonds_from_radii,contact_spring_coef,fracture_criterion, use_berg_origin_basins, nbasins, &
          debug_write,cdrag_grounding,h_to_init_grounding,frac_thres_scaling,frac_thres_n,frac_thres_t,save_bond_traj,&
          remove_unused_bergs,force_convergence,explicit_inner_mts,convergence_tolerance,dem,ignore_tangential_force,poisson,&
          dem_spring_coef,dem_damping_coef,dem_beam_test,constant_interaction_LW,constant_length,constant_width,&
@@ -856,7 +869,7 @@ subroutine ice_bergs_framework_init(bergs, &
          no_frac_first_ts, use_grounding_torque, short_step_mts_grounding, radius_based_drag, save_bond_forces
 
 ! Local variables
-integer :: ierr, iunit, i, j, id_class, axes3d(3), is,ie,js,je,np
+integer :: ierr, iunit, i, j, id_class, axes3d(3), axes3d_b(3), is,ie,js,je,np
 type(icebergs_gridded), pointer :: grd
 real :: lon_mod, big_number
 logical :: lerr
@@ -1015,6 +1028,7 @@ subroutine ice_bergs_framework_init(bergs, &
   allocate( grd%parity_x(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%parity_x(:,:)=1.
   allocate( grd%parity_y(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%parity_y(:,:)=1.
   allocate( grd%iceberg_counter_grd(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%iceberg_counter_grd(:,:)=0
+  allocate( grd%ice_sheet_basins(grd%isd:grd%ied, grd%jsd:grd%jed) ); grd%ice_sheet_basins(:,:)=0.
 
  !write(stderrunit,*) 'KID: copying grid'
   ! Copy data declared on ice model computational domain
@@ -1286,6 +1300,11 @@ subroutine ice_bergs_framework_init(bergs, &
   buffer_width_traj=buffer_width_traj+1
 endif
 
+if (use_berg_origin_basins) then
+  buffer_width=buffer_width+1
+  buffer_width_traj=buffer_width_traj+1
+endif
+
 !must use verlet with mts - Alex
 if (mts) then
   buffer_width_traj=buffer_width_traj+4 !to accomodate n_bonds,fast-step accel terms, overall accel terms
@@ -1377,6 +1396,7 @@ subroutine ice_bergs_framework_init(bergs, &
   bergs%apply_thickness_cutoff_to_gridded_melt=apply_thickness_cutoff_to_gridded_melt
   bergs%melt_cutoff=melt_cutoff
   bergs%read_ocean_depth_from_file=read_ocean_depth_from_file
+  bergs%use_berg_origin_basins=use_berg_origin_basins
   bergs%const_gamma=const_gamma
   bergs%Gamma_T_3EQ=Gamma_T_3EQ
   bergs%pass_fields_to_ocean_model=pass_fields_to_ocean_model
@@ -1564,6 +1584,16 @@ subroutine ice_bergs_framework_init(bergs, &
   id_class = diag_axis_init('mass_class', initial_mass, 'kg','Z', 'iceberg mass')
   axes3d(1:2)=axes
   axes3d(3)=id_class
+
+  if (.not. bergs%use_berg_origin_basins) nbasins=1
+  nbasins=max(nbasins,1)
+  bergs%nbasins=nbasins
+  allocate(basin_arr(bergs%nbasins))
+  basin_arr = (/ (I, I = 1, bergs%nbasins) /)
+  id_class = diag_axis_init('basin_id', basin_arr, 'none','Z', 'iceberg source basins')
+  axes3d_b(1:2)=axes
+  axes3d_b(3)=id_class
+
   grd%id_calving=register_diag_field('icebergs', 'calving', axes, Time, &
      'Incoming Calving mass rate', 'kg/s')
   grd%id_calving_hflx_in=register_diag_field('icebergs', 'calving_hflx_in', axes, Time, &
@@ -1658,9 +1688,14 @@ subroutine ice_bergs_framework_init(bergs, &
      'Y-stress on ice from atmosphere', 'N m^-2')
   grd%id_ocean_depth=register_diag_field('icebergs', 'Depth', axes, Time, &
      'Ocean Depth', 'm')
+  grd%id_ice_sheet_basins=register_diag_field('icebergs', 'ice_sheet_basins', axes, Time, &
+     'Ice-sheet basins of origin for icebergs', 'none')
   grd%id_melt_by_class=register_diag_field('icebergs', 'melt_by_class', axes3d, Time, &
      'Total ice melt (bergs+bits+FL_bits) by class (z-axis labels correspond to Southern hemisphere classes)', &
      'kg/(m^2*s)')
+  grd%id_melt_by_ice_sheet_basin=register_diag_field('icebergs', 'melt_by_ice_sheet_basin', axes3d_b, Time, &
+     'Total ice melt (bergs+bits+FL_bits) by ice-sheet basin', &
+     'kg/(m^2*s)')
   grd%id_melt_buoy_fl=register_diag_field('icebergs', 'melt_buoy_fl', axes, Time, &
      'Buoyancy component of footloose iceberg melt rate', 'kg/(m^2*s)')
   grd%id_melt_eros_fl=register_diag_field('icebergs', 'melt_eros_fl', axes, Time, &
@@ -3331,6 +3366,10 @@ subroutine pack_berg_into_buffer2(berg, buff, n, max_bonds_in)
     call push_buffer_value(buff%data(:,n), counter, berg%rot)
   endif
 
+  if (use_berg_origin_basins) then
+    call push_buffer_value(buff%data(:,n), counter, berg%basin)
+  endif
+
   if (max_bonds .gt. 0) then
     current_bond=>berg%first_bond
     do k = 1,max_bonds
@@ -3503,6 +3542,8 @@ subroutine unpack_berg_from_buffer2(bergs, buff, n, grd, force_append, max_bonds
 
   if (dem) allocate(localberg%ang_vel,localberg%ang_accel,localberg%rot)
 
+  if (use_berg_origin_basins) allocate(localberg%basin)
+
   counter = 0
   call pull_buffer_value(buff%data(:,n), counter, localberg%lon)
   call pull_buffer_value(buff%data(:,n), counter, localberg%lat)
@@ -3570,6 +3611,10 @@ subroutine unpack_berg_from_buffer2(bergs, buff, n, grd, force_append, max_bonds
     call pull_buffer_value(buff%data(:,n), counter, localberg%rot)
   endif
 
+  if (use_berg_origin_basins) then
+    call pull_buffer_value(buff%data(:,n), counter, localberg%basin)
+  endif
+
   !These quantities no longer need to be passed between processors
   localberg%uvel_old=localberg%uvel
   localberg%vvel_old=localberg%vvel
@@ -3825,6 +3870,10 @@ subroutine pack_traj_into_buffer2(traj, buff, n, save_short_traj, save_fl_traj)
       call push_buffer_value(buff%data(:,n), counter, traj%ang_accel)
       call push_buffer_value(buff%data(:,n), counter, traj%rot)
     endif
+
+    if (use_berg_origin_basins) then
+      call push_buffer_value(buff%data(:,n), counter, traj%basin)
+    endif
   endif
 
 end subroutine pack_traj_into_buffer2
@@ -3851,6 +3900,7 @@ subroutine unpack_traj_from_buffer2(first, buff, n, save_short_traj, save_fl_tra
   endif
 
   if (dem) allocate(traj%ang_vel,traj%ang_accel,traj%rot)
+  if (use_berg_origin_basins) allocate(traj%basin)
 
   counter = 0
   call pull_buffer_value(buff%data(:,n),counter,traj%lon)
@@ -3918,6 +3968,10 @@ subroutine unpack_traj_from_buffer2(first, buff, n, save_short_traj, save_fl_tra
       call pull_buffer_value(buff%data(:,n), counter, traj%ang_accel)
       call pull_buffer_value(buff%data(:,n), counter, traj%rot)
     endif
+
+    if (use_berg_origin_basins) then
+      call pull_buffer_value(buff%data(:,n), counter, traj%basin)
+    endif
   endif
   call append_posn(first, traj)
 
@@ -4467,6 +4521,7 @@ subroutine create_iceberg(berg, bergvals)
   endif
 
   if (dem) allocate(berg%ang_vel,berg%ang_accel,berg%rot)
+  if (use_berg_origin_basins) allocate(berg%basin)
 
   berg=bergvals
   berg%prev=>null()
@@ -5355,6 +5410,7 @@ subroutine record_posn(bergs)
   endif
 
   if (dem) allocate(posn%ang_vel,posn%ang_accel,posn%rot)
+  if (use_berg_origin_basins) allocate(posn%basin)
 
   if (save_bond_traj .and. dem) allocate(bond_posn%tangd1,bond_posn%tangd2,bond_posn%nstress,&
                                          bond_posn%sstress,bond_posn%rel_rotation,bond_posn%broken)
@@ -5449,6 +5505,10 @@ subroutine record_posn(bergs)
             posn%ang_accel=this%ang_accel
             posn%rot=this%rot
           endif
+
+          if (use_berg_origin_basins) then
+            posn%basin=this%basin
+          endif
         endif
 
         call push_posn(this%trajectory, posn)
@@ -5516,6 +5576,7 @@ subroutine push_posn(trajectory, posn_vals)
   endif
 
   if (dem) allocate(new_posn%ang_vel,new_posn%ang_accel,new_posn%rot)
+  if (use_berg_origin_basins) allocate(new_posn%basin)
 
   new_posn=posn_vals
   new_posn%next=>trajectory
@@ -5562,6 +5623,7 @@ subroutine append_posn(trajectory, posn_vals)
   endif
 
   if (dem) allocate(new_posn%ang_vel,new_posn%ang_accel,new_posn%rot)
+  if (use_berg_origin_basins) allocate(new_posn%basin)
 
   new_posn=posn_vals
   new_posn%next=>null()