From c1ca6fb25807968b5142c40b91c04e17f90deff0 Mon Sep 17 00:00:00 2001
From: alex-huth <alexander.huth@noaa.gov>
Date: Mon, 5 Aug 2024 16:51:31 -0400
Subject: [PATCH 1/3] Added capability to read in a gridded field of ice-sheet
 basin IDs (e.g. IMBIE basins), which are then associated with bergs as they
 calve from the basins. Note each basin must extend at least one ocean cell
 past the ice front. Then, the iceberg melt associated with each basin can
 also be saved

---
 src/icebergs.F90           | 41 ++++++++++++++++++-
 src/icebergs_fms2io.F90    | 84 ++++++++++++++++++++++++++++++++++----
 src/icebergs_framework.F90 | 68 ++++++++++++++++++++++++++++--
 3 files changed, 180 insertions(+), 13 deletions(-)

diff --git a/src/icebergs.F90 b/src/icebergs.F90
index b10b4f8..796762f 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,7 +117,7 @@ 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
+  integer :: nbonds, nbasins
   logical :: check_bond_quality
   integer :: stdlogunit, stderrunit
 
@@ -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,15 @@ 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)
+  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 +3143,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 +5170,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 +5632,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 +6385,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 +6602,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 +8265,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..4d85fca 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 true, save ice-sheet basin of origin for each berg and track berg melt associated with each basin
 
 !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() !< Basins of origin for bergs (e.g. IMBIE basins, extended to cover calving cells)
+  real, dimension(:,:,:), pointer :: melt_by_ice_sheet_basin=>null() !< Total icebergs melt rate associated with each ice-sheet basin
   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 true, save ice-sheet basin of origin for each berg and track berg melt associated with each basin
+  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
@@ -777,6 +786,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)
@@ -840,7 +851,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 +867,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 +1026,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 +1298,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 +1394,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 +1582,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 +1686,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 +3364,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 +3540,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 +3609,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 +3868,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 +3898,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 +3966,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 +4519,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 +5408,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 +5503,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 +5574,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 +5621,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()

From 734775dec60cf4d54afc769b7aa82f13863f263a Mon Sep 17 00:00:00 2001
From: alex-huth <alexander.huth@noaa.gov>
Date: Thu, 15 Aug 2024 11:09:31 -0400
Subject: [PATCH 2/3] Added ice_sheet_basins_static as a static diag field,
 rather than only having it as a dynamic field

---
 src/icebergs.F90 | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/icebergs.F90 b/src/icebergs.F90
index 796762f..a50b58f 100644
--- a/src/icebergs.F90
+++ b/src/icebergs.F90
@@ -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, nbasins
-  logical :: check_bond_quality
+  integer :: nbonds, nbasins, id_class
+  logical :: check_bond_quality, lerr
   integer :: stdlogunit, stderrunit
 
   ! Get the stderr and stdlog unit numbers
@@ -154,6 +154,9 @@ subroutine icebergs_init(bergs, &
   ! 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

From 23c41cc90a5ae4a165c9837fafaee0debeb80d8e Mon Sep 17 00:00:00 2001
From: alex-huth <alexander.huth@noaa.gov>
Date: Mon, 19 Aug 2024 10:54:13 -0400
Subject: [PATCH 3/3] Units on melt_by_ice_sheet_basin and some shortening of
 comments

---
 src/icebergs_framework.F90 | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90
index 4d85fca..ed7f1e0 100644
--- a/src/icebergs_framework.F90
+++ b/src/icebergs_framework.F90
@@ -62,7 +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 true, save ice-sheet basin of origin for each berg and track berg melt associated with each basin
+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
@@ -146,8 +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() !< Basins of origin for bergs (e.g. IMBIE basins, extended to cover calving cells)
-  real, dimension(:,:,:), pointer :: melt_by_ice_sheet_basin=>null() !< Total icebergs melt rate associated with each ice-sheet basin
+  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)
@@ -525,7 +525,7 @@ 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 true, save ice-sheet basin of origin for each berg and track berg melt associated with each basin
+  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
 
@@ -612,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
@@ -824,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