Merge pull request #70 from alex-huth/basins

Added capability to track melt-by-ice-sheet-basin. Fixed consistency under different PE layouts for positioning of new footloose bergs.

Author: alex-huth

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
@@ -2545,14 +2559,21 @@ subroutine footloose_calving(bergs, time)
     l_c  = pi/(2.*sqrt(2.)) !for length-scale of child iceberg
     lw_c = 1./(gravity*rho_seawater) !for buoyancy length
     B_c  = youngs/(12.*(1.-poisson**2.)) !for bending stiffness
-    seed = constructSeed(mpp_pe(),mpp_pe(),time) !Seed random numbers for Poisson distribution
-    rns = initializeRandomNumberStream(seed)
-    call getRandomNumbers(rns, rn)
+    if (bergs%fl_init_child_xy_by_pe) then !old bug that does not reproduce on different PE layouts
+      seed = constructSeed(mpp_pe(),mpp_pe(),time) !Seed random numbers for Poisson distribution
+      rns = initializeRandomNumberStream(seed)
+      call getRandomNumbers(rns, rn)
+    endif
     Visited=.true.
   endif
 
   do grdj = grd%jsc,grd%jec ; do grdi = grd%isc,grd%iec !computational domain only
     this=>bergs%list(grdi,grdj)%first
+    if ((.not. bergs%fl_init_child_xy_by_pe) .and. associated(this)) then
+      ! Seed random numbers based on space and "time"
+      rns = initializeRandomNumberStream( grdi + 10000*grdj + &
+                                          int( 16384.*abs( sin(262144.*grd%ssh(grdi,grdj)) ) ) )
+    endif
     do while(associated(this))
 
       !only non-static, non-footloose-child bergs are eligible for footloose calving:
@@ -3125,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
@@ -5145,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.
@@ -5606,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) &
@@ -6355,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
@@ -6565,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
 
@@ -8223,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)

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)