Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Iceberg basin-of-origin tracking when using FMS1 #72

Open
wants to merge 1 commit into
base: dev/gfdl
Choose a base branch
from
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
83 changes: 76 additions & 7 deletions src/icebergs_fmsio.F90
Original file line number Diff line number Diff line change
Expand Up @@ -41,12 +41,12 @@ module ice_bergs_fmsio
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
Expand All @@ -56,7 +56,7 @@ module ice_bergs_fmsio
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
Expand Down Expand Up @@ -181,7 +181,8 @@ subroutine write_restart_bergs(bergs, time_stamp)
first_berg_ine, &
other_berg_jne, &
other_berg_ine, &
broken
broken, &
basin


integer :: grdi, grdj
Expand Down Expand Up @@ -247,6 +248,7 @@ subroutine write_restart_bergs(bergs, time_stamp)
allocate(ang_accel(nbergs))
allocate(rot(nbergs))
endif
if (use_berg_origin_basins) allocate(basin(nbergs))

filename = trim("icebergs.res.nc")
call set_domain(bergs%grd%domain)
Expand Down Expand Up @@ -322,6 +324,10 @@ subroutine write_restart_bergs(bergs, time_stamp)
id = register_restart_field(bergs_restart,filename,'rot',rot,&
longname='dem accumulated rotation',units='rad')
endif
if (use_berg_origin_basins) then
id = register_restart_field(bergs_restart,filename,'basin',basin,&
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
Expand Down Expand Up @@ -376,6 +382,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
Expand Down Expand Up @@ -426,6 +433,8 @@ subroutine write_restart_bergs(bergs, time_stamp)
rot)
endif

if (use_berg_origin_basins) deallocate(basin)

deallocate( &
ine, &
jne, &
Expand Down Expand Up @@ -655,7 +664,8 @@ subroutine read_restart_bergs(bergs,Time)
iceberg_num,&
id_cnt, &
id_ij, &
start_year
start_year, &
basin

!integer, allocatable, dimension(:,:) :: iceberg_counter_grd

Expand Down Expand Up @@ -738,6 +748,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 read_unlimited_axis(filename,'lon',lon,domain=grd%domain)
call read_unlimited_axis(filename,'lat',lat,domain=grd%domain)
Expand Down Expand Up @@ -784,6 +798,11 @@ subroutine read_restart_bergs(bergs,Time)
call read_real_vector(filename,'ang_accel',ang_accel,grd%domain,value_if_not_in_file=0.)
call read_real_vector(filename,'rot' ,rot ,grd%domain,value_if_not_in_file=0.)
endif

if (use_berg_origin_basins) then
call read_int_vector(filename,'basin',basin,grd%domain,value_if_not_in_file=0)
endif

elseif (bergs%require_restart) then
stop 'read_restart_bergs, RESTART NOT FOUND!'
endif
Expand Down Expand Up @@ -867,6 +886,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)
Expand Down Expand Up @@ -927,6 +950,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)
Expand Down Expand Up @@ -1032,6 +1056,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
Expand Down Expand Up @@ -1081,6 +1106,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)
Expand Down Expand Up @@ -1549,7 +1577,7 @@ subroutine read_ocean_depth(grd)
! Local variables
character(len=37) :: filename

! Read stored ice
! Read depth
filename=trim(restart_input_dir)//'topog.nc'
if (file_exist(filename)) then
if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') &
Expand All @@ -1571,6 +1599,33 @@ 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

! Read sub_basin
filename=trim(restart_input_dir)//'ice_sheet_basins.nc'
if (file_exist(filename)) then
if (mpp_pe().eq.mpp_root_pe()) write(*,'(2a)') &
'KID, read_ice_sheet_basins: reading ',filename
if (field_exist(filename, '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(filename, 'sub_basin', grd%ice_sheet_basins, grd%domain)
else
if (verbose.and.mpp_pe().eq.mpp_root_pe()) write(*,'(a)') &
'KID, read_ice_sheet_basins: sub_basin WAS NOT FOUND in the file. Setting to 0.'
! grd%ice_sheet_basins(:,:)=0.
endif
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
Expand All @@ -1586,7 +1641,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
Expand Down Expand Up @@ -1763,6 +1818,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
Expand Down Expand Up @@ -1833,6 +1891,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
Expand Down Expand Up @@ -1950,6 +2012,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

Expand Down Expand Up @@ -2031,6 +2097,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)
Expand Down
Loading