Skip to content

Commit

Permalink
Fixed the berg deletion error, which was caused by nullifying a cell'…
Browse files Browse the repository at this point in the history 
…s entire berg list when some non-fully-calved bergs should have still been present in the cell. Cleaned up most of the debugging print statements.
  • Loading branch information
@alex-huth
alex-huth committed May 21, 2024
1 parent f7d2665 commit af6f97e
Show file tree
Hide file tree
Showing 3 changed files with 24 additions and 95 deletions.
6 changes: 1 addition & 5 deletions src/ice_shelf_calving.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ module ice_shelf_tabular_calving
use ice_bergs_framework, only : square_into_quadrants_using_triangles
use ice_bergs_framework, only : initialize_iceberg_bonds, count_bonds
use ice_bergs_framework, only : find_cell, pos_within_cell, generate_id
use ice_bergs_framework, only : debug, footloose, connect_all_bonds, delete_all_bonds, berg_exists
use ice_bergs_framework, only : debug, footloose, connect_all_bonds, delete_all_bonds
use ice_bergs_framework, only : update_halo_calved_tabular_icebergs, assign_n_bonds,transfer_mts_bergs
use fms_mod, only : error_mesg, FATAL, WARNING, stderr

Expand Down Expand Up @@ -734,7 +734,6 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC)
! call transfer_mts_bergs(bergs)
! else
!call update_halo_icebergs(bergs)
if (berg_exists(bergs)) print *,'BE: pre connect CALVING'

!includes all particles that just initialized
call connect_all_bonds(bergs, match_bond_pairs=.true.,tabular_calving_only=.true.)
Expand All @@ -748,9 +747,6 @@ subroutine ice_shelf_to_bonded_bergs(bergs, TC)

!only includes particles that just initialized
call assign_n_bonds(bergs,tabular_calving_only=.true.)

if (berg_exists(bergs)) print *,'BE: post bonding CALVING'

endif

!This can be done with the rest of the bergs?
Expand Down
18 changes: 6 additions & 12 deletions src/icebergs.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,6 @@ module ice_bergs
use ice_bergs_framework, only: point_in_triangle, point_in_interval, point_is_on_the_line
use ice_bergs_framework, only: convert_from_grid_to_meters, convert_from_meters_to_grid
use ice_bergs_framework, only: spread_variable_across_cells, find_orientation_using_iceberg_bonds
use ice_bergs_framework, only: berg_exists

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
Expand Down Expand Up @@ -162,6 +161,7 @@ subroutine icebergs_init(bergs, &

!Reading ocean depth from a file
if (bergs%read_ocean_depth_from_file) call read_ocean_depth(bergs%grd)
call mpp_update_domains(bergs%grd%ocean_depth, bergs%grd%domain)

if (bergs%iceberg_bonds_on) then
if (bergs%manually_initialize_bonds) then
Expand Down Expand Up @@ -4907,7 +4907,6 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh,

if (.not.bergs%Static_icebergs) then
! call assign_n_bonds(bergs) ! for debugging tabular calving!
if (berg_exists(bergs)) print *,'BE: pre-evolve-icebergs'
if (bergs%mts) then
call evolve_icebergs_mts(bergs)
else
Expand All @@ -4934,20 +4933,15 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh,
if (bergs%footloose) call footloose_calving(bergs, time)
call mpp_clock_end(bergs%clock_fl1)

if (berg_exists(bergs)) print *,'BE: post_send_bergs_to_other_pes'

! Calving of tabular bonded bergs
! TODO: does this make sense here?
if (bergs%tabular_calving) then
call process_tabular_calving(bergs)
if (bergs%iceberg_bonds_on) call bond_address_update(bergs)
!if (bergs%iceberg_bonds_on) call bond_address_update(bergs)
!return gridded variables associated with tabular calving
frac_cberg_calved(:,:)=TC%frac_cberg_calved(grd%isc:grd%iec,grd%jsc:grd%jec)
frac_cberg(:,:) =TC%frac_cberg(grd%isc:grd%iec,grd%jsc:grd%jec)
endif

if (berg_exists(bergs)) print *,'BE: post_process_tabular_calving'

call mpp_clock_begin(bergs%clock_com2)
if (bergs%mts) then
call interp_gridded_fields_to_bergs(bergs)
Expand Down Expand Up @@ -6499,10 +6493,10 @@ subroutine evolve_icebergs_mts(bergs)
xi=berg%xi ; yj=berg%yj
! finalize new iceberg positions and index
call adjust_index_and_ground(grd, lonn, latn, uveln, vveln, i, j, xi, yj, bounced, error_flag, berg%id)
if (xi .ne. xi) then
print *,'id',berg%id,'EIM xi0',berg%xi,'xi1',xi,'u,v',uveln,vveln,'lonlat',berg%lon,berg%lat,&
'msk0',grd%msk(berg%ine,berg%jne),'msk1',grd%msk(i,j),'i,j 0',berg%ine,berg%jne,'i,j 1',i,j,'lonn,latn',lonn,latn
endif
!!$ if (xi .ne. xi) then
!!$ print *,'id',berg%id,'EIM xi0',berg%xi,'xi1',xi,'u,v',uveln,vveln,'lonlat',berg%lon,berg%lat,&
!!$ 'msk0',grd%msk(berg%ine,berg%jne),'msk1',grd%msk(i,j),'i,j 0',berg%ine,berg%jne,'i,j 1',i,j,'lonn,latn',lonn,latn
!!$ endif
berg%lon=lonn ; berg%lat=latn
berg%lon_old=lonn ; berg%lat_old=latn
berg%ine=i ; berg%jne=j
Expand Down
95 changes: 17 additions & 78 deletions src/icebergs_framework.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1582,12 +1582,6 @@ subroutine ice_bergs_framework_init(bergs, &

if (present(tabular_calving)) then
bergs%tabular_calving=tabular_calving
! if (bergs%tabular_calving) then
! if (.not. (bergs%mts .and. bergs%dem .and. (.not. bergs%old_interp_flds_order))) then
! call error_mesg('KID, ice_bergs_framework_init', &
! 'tabular calving requires (mts .and. dem .and. (.not. old_interp_flds_order))!', FATAL)
! endif
! endif
else
bergs%tabular_calving=.false.
endif
Expand Down Expand Up @@ -2490,8 +2484,6 @@ Subroutine transfer_mts_bergs(bergs)
grd=>bergs%grd ! for convenience
stderrunit = stderr() ! Get the stderr unit number

if (berg_exists(bergs)) print *,'BE: pre clear halo'

! Step 1: Clear the current halos
do grdj = grd%jsd,grd%jsc-1 ; do grdi = grd%isd,grd%ied
call delete_all_bergs_in_list(bergs,grdj,grdi)
Expand Down Expand Up @@ -2524,41 +2516,30 @@ Subroutine transfer_mts_bergs(bergs)
enddo
enddo;enddo

if (berg_exists(bergs)) print *,'BE: pre connect'

!Copy bergs between PEs
do i = 1,2 !run twice to account for diagonal transfers and guarantee robust transfers of conglomerates
if (mpp_pe()==6) print *,'connect a',i
call connect_all_bonds(bergs,ignore_unmatched=.true.)
nbergs_to_send_e=0; nbergs_to_send_w=0; nbergs_to_send_n=0; nbergs_to_send_s=0

if (berg_exists(bergs)) print *,'BE: pre mts_pack_in_dir',i

call mpp_sync_self()

call mts_pack_in_dir(bergs,nbergs_to_send_e,"e")
call mts_pack_in_dir(bergs,nbergs_to_send_w,"w")
call mts_pack_in_dir(bergs,nbergs_to_send_n,"n")
call mts_pack_in_dir(bergs,nbergs_to_send_s,"s")

if (berg_exists(bergs)) print *,'BE: post mts_pack_in_dir',i

call mts_send_and_receive(bergs,nbergs_to_send_e,nbergs_to_send_w,nbergs_to_send_n,nbergs_to_send_s)

call mpp_sync_self()
if (berg_exists(bergs)) print *,'BE: post_mts_send_and_receive',i
enddo

if (debug) then
call connect_all_bonds(bergs,ignore_unmatched=.false.,match_bond_pairs=.true.)
else
if (mpp_pe()==6) print *,'connect b',i
call connect_all_bonds(bergs,ignore_unmatched=.true.,match_bond_pairs=.true.)
! call connect_all_bonds(bergs,ignore_unmatched=.false.,match_bond_pairs=.true.)
endif

if (berg_exists(bergs)) print *,'BE: post-connect'

if (bergs%tabular_calving) then
count_d=0
count_k=0
Expand All @@ -2567,22 +2548,19 @@ Subroutine transfer_mts_bergs(bergs)
do while (associated(this))
if (this%mass_scaling == -1) then
! print *,'berg deleted',mpp_pe(),this%id
if (this%id==4294972978) print *,'4294972978 deleted on PE',mpp_pe(),grdi,grdj,this%halo_berg,this%mass_scaling
kick_the_bucket=>this
this=>this%next
call delete_all_bonds(kick_the_bucket)
call delete_iceberg_from_list(bergs%list(grdi,grdj)%first,kick_the_bucket)
count_d=count_d+1
else
! print *,'berg not deleted',mpp_pe(),this%id
if (this%id==4294972978) print *,'4294972978 exists on PE',mpp_pe(),grdi,grdj,this%halo_berg,this%mass_scaling
count_k=count_k+1
this=>this%next
endif
enddo
enddo;enddo
if (count_k+count_d>0) then
print *,'PE',mpp_pe(),'bergs kept', count_k, 'bergs deleted',count_d
endif
endif

Expand All @@ -2601,12 +2579,13 @@ Subroutine transfer_mts_bergs(bergs)
call show_all_bonds(bergs)
endif

if (berg_exists(bergs)) print *,'BE: transfer_mts_bergs'
end subroutine transfer_mts_bergs

!> Debugging function that detects if a berg exists on the current PE within a cell on within any bond
!! Edit the berg id (default id==4294972978) accordingly within the function
function berg_exists(bergs)
! Arguments
type(logical) :: berg_exists
type(logical) :: berg_exists !> True if berg is detected
type(icebergs), pointer :: bergs !< Container for all types and memory
type(icebergs_gridded), pointer :: grd
integer :: grdi,grdj,bond_count
Expand All @@ -2629,7 +2608,14 @@ function berg_exists(bergs)
endif
bond=>this%first_bond
do while (associated(bond))
if (abs(bond%other_id) .eq. 4294972978) bond_count=bond_count+1
if (abs(bond%other_id) .eq. 4294972978) then
bond_count=bond_count+1
elseif (associated(bond%other_berg)) then
if (abs(bond%other_berg%id) .eq. 4294972978) then
bond_count=bond_count+1
print *,'BERG 4294972978 exists as bond%other_berg%id, but not as bond%other_id'
endif
endif
bond=>bond%next_bond
enddo
this=>this%next
Expand Down Expand Up @@ -2806,22 +2792,18 @@ subroutine mts_pack_in_dir(bergs, nbergs_to_send, dir)
nbergs_to_send=nbergs_to_send+1
select case (dir)
case ("s")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send S from edgecontact',mpp_pe()
berg%conglom_id=2
call pack_berg_into_buffer2(berg, bergs%obuffer_s, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+1
case ("n")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send N from edgecontact',mpp_pe()
berg%conglom_id=1
call pack_berg_into_buffer2(berg, bergs%obuffer_n, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+2
case ("w")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send W from edgecontact',mpp_pe()
berg%lon=berg%lon-pfix; berg%conglom_id=4
call pack_berg_into_buffer2(berg, bergs%obuffer_w, nbergs_to_send, bergs%max_bonds)
berg%lon=berg%lon+pfix; berg%conglom_id=current_conglom_id+8
case ("e")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send E from edgecontact',mpp_pe()
berg%lon=berg%lon-pfix; berg%conglom_id=8
call pack_berg_into_buffer2(berg, bergs%obuffer_e, nbergs_to_send, bergs%max_bonds)
berg%lon=berg%lon+pfix; berg%conglom_id=current_conglom_id+4
Expand Down Expand Up @@ -2851,7 +2833,7 @@ recursive subroutine mts_mark_and_pack_halo_and_congloms(bergs, berg, dir, nberg
integer :: k !<bond counter
integer :: current_conglom_id
real :: current_halo_id
logical :: tf

!pack the berg for transfer if it has not been packed already
if (.not. mts_berg_sent(berg%conglom_id,dir)) then
current_halo_id=berg%halo_berg; current_conglom_id=berg%conglom_id
Expand All @@ -2871,22 +2853,18 @@ recursive subroutine mts_mark_and_pack_halo_and_congloms(bergs, berg, dir, nberg

select case (dir)
case ("e")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send E',mpp_pe()
berg%conglom_id=8; berg%lon=berg%lon-pfix;
call pack_berg_into_buffer2(berg,bergs%obuffer_e, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+4; berg%lon=berg%lon+pfix;
case ("w")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send W',mpp_pe()
berg%conglom_id=4; berg%lon=berg%lon-pfix;
call pack_berg_into_buffer2(berg,bergs%obuffer_w, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+8; berg%lon=berg%lon+pfix;
case ("n")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send N',mpp_pe()
berg%conglom_id=1
call pack_berg_into_buffer2(berg,bergs%obuffer_n, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+2
case ("s")
if (abs(berg%id)==4294972978) print *,'BERG 4294972978 packed to send S',mpp_pe()
berg%conglom_id=2
call pack_berg_into_buffer2(berg,bergs%obuffer_s, nbergs_to_send, bergs%max_bonds)
berg%conglom_id=current_conglom_id+1
Expand All @@ -2905,25 +2883,6 @@ recursive subroutine mts_mark_and_pack_halo_and_congloms(bergs, berg, dir, nberg
if (associated(current_bond%other_berg)) then
other_berg=>current_bond%other_berg
if (other_berg%id>0) then
if (other_berg%id==4294972978) then
print *,'BERG 4294972978 marked to send from connection',mpp_pe(),'by berg',berg%id,'. CID',other_berg%conglom_id
print *,'sending berg stats',mpp_pe(),berg%id,berg%lon,berg%lat,berg%ine,berg%jne

if (.not. berg_exists(bergs)) then
print *,'but berg exists still claims its not here'
this=>bergs%list(berg%ine,berg%jne)%first
if (.not. associated(this)) print *,'there arent even any bergs in the cell....'
do while (associated(this))
if (this%id==berg%id) then
print *,'sending berg found',mpp_pe(),this%id
this=>this%next
else
print *,'berg in cell',mpp_pe(),this%id
this=>this%next
endif
enddo
endif
endif
call mts_mark_and_pack_halo_and_congloms(bergs,other_berg,dir,nbergs_to_send,pfix,x,y)!,rhc)
endif
endif
Expand Down Expand Up @@ -3017,22 +2976,18 @@ recursive subroutine mts_pack_contact_bergs(bergs, berg, dir, pfix, nbergs_to_se
nbergs_to_send=nbergs_to_send+1
select case (dir)
case ("s")
if (abs(other_berg%id)==4294972978) print *,'BERG 4294972978 packed to send S from contact',mpp_pe()
other_berg%conglom_id=2
call pack_berg_into_buffer2(other_berg,bergs%obuffer_s, nbergs_to_send, bergs%max_bonds)
other_berg%conglom_id=current_conglom_id+1
case ("n")
if (abs(other_berg%id)==4294972978) print *,'BERG 4294972978 packed to send N from contact',mpp_pe()
other_berg%conglom_id=1
call pack_berg_into_buffer2(other_berg,bergs%obuffer_n, nbergs_to_send, bergs%max_bonds)
other_berg%conglom_id=current_conglom_id+2
case ("w")
if (abs(other_berg%id)==4294972978) print *,'BERG 4294972978 packed to send W from contact',mpp_pe()
other_berg%lon=other_berg%lon-pfix; other_berg%conglom_id=4
call pack_berg_into_buffer2(other_berg,bergs%obuffer_w, nbergs_to_send, bergs%max_bonds)
other_berg%lon=other_berg%lon+pfix; other_berg%conglom_id=current_conglom_id+8
case ("e")
if (abs(other_berg%id)==4294972978) print *,'BERG 4294972978 packed to send E from contact',mpp_pe()
other_berg%lon=other_berg%lon-pfix; other_berg%conglom_id=8
call pack_berg_into_buffer2(other_berg,bergs%obuffer_e, nbergs_to_send, bergs%max_bonds)
other_berg%lon=other_berg%lon+pfix; other_berg%conglom_id=current_conglom_id+4
Expand Down Expand Up @@ -3472,7 +3427,7 @@ subroutine delete_all_bergs_in_list(bergs, grdj, grdi, tabular_calving_only)
logical, optional :: tabular_calving_only !< true to only delete tabular bergs that are calving from the ice shelf
! Local variables
type(iceberg), pointer :: kick_the_bucket, this
logical :: new_tab_only
logical :: new_tab_only !< local version of tabular_calving_only

new_tab_only=.false.
if (present(tabular_calving_only)) then
Expand All @@ -3481,16 +3436,17 @@ subroutine delete_all_bergs_in_list(bergs, grdj, grdi, tabular_calving_only)

this=>bergs%list(grdi,grdj)%first
do while (associated(this))
if (new_tab_only .and. this%static_berg>=0 .and. this%halo_berg<=1) then! .and. this%static_berg<2) then
if (new_tab_only .and. this%static_berg>=0 .and. this%halo_berg<=1) then
this=>this%next
else
kick_the_bucket=>this
this=>this%next
call destroy_iceberg(kick_the_bucket)
call delete_iceberg_from_list(bergs%list(grdi,grdj)%first,kick_the_bucket)
! call destroy_iceberg(kick_the_bucket)
endif
!call delete_iceberg_from_list(bergs%list(grdi,grdj)%first,kick_the_bucket)
enddo
bergs%list(grdi,grdj)%first=>null()
! if (.not. new_tab_only) bergs%list(grdi,grdj)%first=>null()
end subroutine delete_all_bergs_in_list


Expand Down Expand Up @@ -5400,45 +5356,29 @@ subroutine delete_all_bonds(berg)
type(iceberg), intent(inout), pointer :: berg !<parent berg to delete associated bonds
type(iceberg), pointer :: other_berg
type(bond), pointer :: current_bond, matching_bond, kick_the_bucket
integer :: obbd

obbd=0
current_bond=>berg%first_bond
do while (associated(current_bond))
if (associated(current_bond%other_berg)) then
other_berg=>current_bond%other_berg
matching_bond=>other_berg%first_bond
if (berg%id==4294972978) then
print *,'other_berg%id',other_berg%id,'on PE',mpp_pe()
if (.not. associated(matching_bond)) print *,'matching_bond not associated!'
endif
do while (associated(matching_bond)) ! Looping over possible matching bonds in other_berg
if (matching_bond%other_id .eq. berg%id) then
if (berg%id==4294972978) then
obbd=obbd+1
print *,'matched bond for other_berg%id',other_berg%id,'on PE',mpp_pe()
endif
kick_the_bucket=>matching_bond
matching_bond=>matching_bond%next_bond
call delete_bond_from_list(other_berg,kick_the_bucket)
! matching_bond=>null()
other_berg%n_bonds=other_berg%n_bonds-1
else
if (berg%id==4294972978) print *,'matching_bond%other_id',matching_bond%other_id
matching_bond=>matching_bond%next_bond
endif
enddo
else
if (berg%id==4294972978) print *,'other berg not associated for the bond!',mpp_pe()
endif
kick_the_bucket=>current_bond
current_bond=>current_bond%next_bond
call delete_bond_from_list(berg,kick_the_bucket)
berg%n_bonds=berg%n_bonds-1
enddo
if (berg%id==4294972978) then
print *,'remaining bonds of deleting berg',berg%id,'=',berg%n_bonds, 'matched bonds deleted = ',obbd
endif
end subroutine delete_all_bonds

!> Bond two bergs together
Expand Down Expand Up @@ -5749,7 +5689,6 @@ subroutine connect_all_bonds(bergs, ignore_unmatched, match_bond_pairs, tabular_
endif
enddo
else
print *,'pe',mpp_pe(),'berg id,latlon', berg%id, berg%lat, berg%lon, 'missing berg',current_bond%other_id
call error_mesg('KID, connect_all_bonds', 'A bond is missing its second berg !!!', WARNING)
endif
endif
Expand Down

0 comments on commit af6f97e

Please sign in to comment.