From 291bdcbff2f280e96fc9c95e0170286df7f2efc0 Mon Sep 17 00:00:00 2001 From: alex-huth Date: Thu, 23 May 2024 12:08:12 -0400 Subject: [PATCH] Fixed error where particles on the edge of recently-calved iceberg conglomerates were moving unphysically: this was related to bouncing occuring when these bergs overlapped masked ice-shelf cells. Fixed by removing bouncing when using tabular_calving and by unmasking any cell that contains a berg particle. The latter unmasking fix is a bit of a hack, but is required for grid_to_berg interpolations and tidal_drift/coastal_drift to work as intended (note that at calving, it is possible for a cell to be partially-full with both icebergs and ice shelf). Overall, this hack works OK as long as tidal_drift or coastal_drift are set such that berg particles will not advect into cells that are full of ice shelf or land. An alterative solution would involve pushing bergs away from the ice front (partially-full or not) using some type of contact force instead of tidal_drift or coastal_drift, and preferably, skipping the unmasking step. To this end, and helpful for debugging, a mask_status variable was also added to the berg and berg trajectory types. While unused at this time, this mask_status could be used in the future, e.g. to track whether a berg can coexist with ice shelf in a masked cell because it is recently calved --- src/ice_shelf_calving.F90 | 20 ++++++- src/icebergs.F90 | 120 +++++++++++++++++-------------------- src/icebergs_fms2io.F90 | 64 ++++++++++++++++++-- src/icebergs_framework.F90 | 52 +++++++++++++++- 4 files changed, 181 insertions(+), 75 deletions(-) diff --git a/src/ice_shelf_calving.F90 b/src/ice_shelf_calving.F90 index 54235c4..f55e0d1 100644 --- a/src/ice_shelf_calving.F90 +++ b/src/ice_shelf_calving.F90 @@ -1172,6 +1172,18 @@ subroutine new_tabular_bergs_thickness_and_pressure(bergs) call mpp_update_domains(frac_cberg, grd%domain, complete=.false.) call mpp_update_domains(grd%msk, grd%domain, complete=.true.) if (allocated(pf_area)) deallocate(pf_area) + + if (.not. bergs%mts) then + !if bergs%mts, then this is done immediately following process_tabular_calving, within + !interp_gridded_fields_to_bergs, so no need to do it here... + do grdj = grd%jsc,grd%jec ; do grdi = grd%isc,grd%iec + berg=>bergs%list(grdi,grdj)%first + do while (associated(berg)) + berg%mask_status=grd%msk(grdi,grdj) + berg=>berg%next + enddo + enddo; enddo + endif end subroutine new_tabular_bergs_thickness_and_pressure !> Initialize (begin calving) a tabular iceberg particle from an ice shelf at the given lat/lon coordinates. @@ -1236,7 +1248,8 @@ subroutine begin_calving_tabular_iceberg_from_shelf(bergs, grd, lon, lat, calve_ call error_mesg('KID, calve_icebergs', 'berg xi,yj is not correct!', FATAL) endif - !Ignore bergs on the N and E boundry of the PE, as they will be included in the PEs to the N or E, respectively + !Ignore bergs on the N and E boundary of the PE, as they will be included in the PEs to the N or E, respectively + !But the find_cell call should not allow bergs to be found on these boundaries, anyway. if ((i==grd%iec .and. xi==1) .or. (j==grd%jec .and. yj==1)) return ! if (grd%msk(i,j)<0.5) then @@ -1439,6 +1452,11 @@ subroutine begin_calving_tabular_iceberg_from_shelf(bergs, grd, lon, lat, calve_ newberg%ang_vel=0.; newberg%ang_accel=0.; newberg%rot=0. endif + if (bergs%tabular_calving) then !(obviously this is true) + if (.not. allocated(newberg%mask_status)) allocate(newberg%mask_status) + newberg%mask_status=grd%msk(i,j) + endif + call add_new_berg_to_list(bergs%list(i,j)%first, newberg) ! calved_to_berg=initial_mass*mass_scaling ! Units of kg ! Heat content TODO diff --git a/src/icebergs.F90 b/src/icebergs.F90 index 3197b82..1720414 100644 --- a/src/icebergs.F90 +++ b/src/icebergs.F90 @@ -1924,6 +1924,7 @@ subroutine accel(bergs, berg, i, j, xi, yj, lat, uvel, vvel, uvel0, vvel0, dt, r ! !end if if (bergs%old_interp_flds_order) then + if (bergs%tabular_calving) berg%mask_status=grd%msk(i,j) call interp_flds(grd, berg%lon, berg%lat, i, j, xi, yj, rx, ry, uo, vo, ui, vi, ua, va, ssh_x, & ssh_y, sst, sss, cn, hi, od) else @@ -2784,11 +2785,12 @@ subroutine thermodynamics(bergs) this=>next else - if (bergs%old_interp_flds_order .or. (.not. mts .and. .not. dem .and. this%halo_berg.ge.0.5)) then - call interp_flds(grd, this%lon, this%lat, this%ine, this%jne, this%xi, this%yj, 0., 0., & - this%uo, this%vo, this%ui, this%vi, this%ua, this%va, this%ssh_x, & - this%ssh_y, this%sst, this%sss,this%cn, this%hi) - end if + if (bergs%old_interp_flds_order .or. (.not. mts .and. .not. dem .and. this%halo_berg.ge.0.5)) then + if (bergs%tabular_calving) this%mask_status=grd%msk(grdi,grdj) + call interp_flds(grd, this%lon, this%lat, this%ine, this%jne, this%xi, this%yj, 0., 0., & + this%uo, this%vo, this%ui, this%vi, this%ua, this%va, this%ssh_x, & + this%ssh_y, this%sst, this%sss,this%cn, this%hi) + end if SST=this%sst SSS=this%sss @@ -4019,51 +4021,6 @@ subroutine spread_mass_across_ocean_cells(bergs, berg, i, j, x, y, Mberg, Mbits, end subroutine spread_mass_across_ocean_cells -!!$!> Distribute a quantity among nine cells on a grid centered at cell i,j -!!$subroutine spread_variable_across_cells(grd, variable_on_ocean, Var, i, j, & -!!$ yDxL, yDxC,yDxR, yCxL, yCxC, yCxR, yUxL, yUxC, yUxR, I_fraction_used) -!!$ ! Arguments -!!$ type(icebergs_gridded), pointer, intent(in) :: grd !< Container for gridded fields -!!$ real, dimension(grd%isd:grd%ied, grd%jsd:grd%jed, 9), intent(inout) :: variable_on_ocean !< Gridded field to augment -!!$ real, intent(in) :: Var !< Variable to be spread accross cell -!!$ real, intent(in) :: yDxL !< Weight for the cell at i-1,j-1 -!!$ real, intent(in) :: yDxC !< Weight for the cell at i-1,j -!!$ real, intent(in) :: yDxR !< Weight for the cell at i-1,j+1 -!!$ real, intent(in) :: yCxL !< Weight for the cell at i,j-1 -!!$ real, intent(in) :: yCxC !< Weight for the cell at i,j -!!$ real, intent(in) :: yCxR !< Weight for the cell at i,j-1 -!!$ real, intent(in) :: yUxL !< Weight for the cell at i+1,j-1 -!!$ real, intent(in) :: yUxC !< Weight for the cell at i+1,j -!!$ real, intent(in) :: yUxR !< Weight for the cell at i+1,j+1 -!!$ real, intent(in) :: I_fraction_used !< Amount of iceberg used (inverse) -!!$ integer, intent(in) :: i !< i-index of cell containing center of berg -!!$ integer, intent(in) :: j !< j-index of cell containing center of berg -!!$ -!!$ !Spreading the iceberg mass onto the ocean -!!$ variable_on_ocean(i,j,1)=variable_on_ocean(i,j,1)+(yDxL*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,2)=variable_on_ocean(i,j,2)+(yDxC*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,3)=variable_on_ocean(i,j,3)+(yDxR*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,4)=variable_on_ocean(i,j,4)+(yCxL*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,5)=variable_on_ocean(i,j,5)+(yCxC*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,6)=variable_on_ocean(i,j,6)+(yCxR*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,7)=variable_on_ocean(i,j,7)+(yUxL*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,8)=variable_on_ocean(i,j,8)+(yUxC*Var*I_fraction_used) -!!$ variable_on_ocean(i,j,9)=variable_on_ocean(i,j,9)+(yUxR*Var*I_fraction_used) -!!$ -!!$end subroutine spread_variable_across_cells - -!!$!> Returns area of a triangle -!!$real function Area_of_triangle(Ax, Ay, Bx, By, Cx, Cy) -!!$ ! Arguments -!!$ real, intent(in) :: Ax !< x-position of corner A -!!$ real, intent(in) :: Ay !< y-position of corner A -!!$ real, intent(in) :: Bx !< x-position of corner B -!!$ real, intent(in) :: By !< y-position of corner B -!!$ real, intent(in) :: Cx !< x-position of corner C -!!$ real, intent(in) :: Cy !< y-position of corner C -!!$ Area_of_triangle = abs( 0.5*((Ax*(By-Cy))+(Bx*(Cy-Ay))+(Cx*(Ay-By))) ) -!!$end function Area_of_triangle - !> Returns x rounded of to sig_fig !! \todo What the heck is this for? -AJA real function roundoff(x,sig_fig) @@ -4111,6 +4068,7 @@ subroutine interp_gridded_fields_to_bergs(bergs) call getRandomNumbers(rns, ry) ry = 2.*ry - 1. endif + if (bergs%tabular_calving) grd%msk(grdi,grdj)=1 !berg%mask_status=grd%msk(grdi,grdj) call interp_flds(grd, berg%lon, berg%lat, berg%ine, berg%jne, berg%xi, berg%yj, rx, ry, berg%uo, berg%vo, & berg%ui, berg%vi, berg%ua, berg%va, berg%ssh_x, berg%ssh_y, berg%sst, berg%sss, berg%cn, berg%hi, berg%od) endif @@ -4510,6 +4468,7 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh, real, dimension(:,:), optional, pointer :: frac_cberg_calved !< Cell fraction of fully-calved bonded bergs from !! the ice sheet [nondim] ! Local variables + type(iceberg), pointer :: berg type(tabular_calving_state), pointer :: TC integer :: iyr, imon, iday, ihr, imin, isec, k type(icebergs_gridded), pointer :: grd @@ -4619,12 +4578,28 @@ subroutine icebergs_run(bergs, time, calving, uo, vo, ui, vi, tauxa, tauya, ssh, else grd%msk(i,j)=1 endif + +!!$ A TEMPORARY WORK-AROUND? +!!$ Adjust mask so that bergs are never in a masked cell + if (grd%msk(i,j).ne.1) then + if (associated(bergs%list(i,j)%first)) then + berg=>bergs%list(i,j)%first + do while (associated(berg)) + if (berg%static_berg==0.) then + grd%msk(i,j)=1 + berg=>NULL() + else + berg=>berg%next + endif + enddo + endif + endif enddo; enddo - call mpp_update_domains(TC%calve_mask, grd%domain)!, complete=.false.) - call mpp_update_domains(TC%h_shelf, grd%domain)!, complete=.false.) - call mpp_update_domains(TC%frac_shelf, grd%domain)!, complete=.false.) - call mpp_update_domains(grd%msk, grd%domain)!, complete=.true.) + call mpp_update_domains(TC%calve_mask, grd%domain, complete=.false.) + call mpp_update_domains(TC%h_shelf, grd%domain, complete=.false.) + call mpp_update_domains(TC%frac_shelf, grd%domain, complete=.false.) + call mpp_update_domains(grd%msk, grd%domain, complete=.true.) TC%frac_cberg_calved(:,:) = 0.0 TC%frac_cberg(:,:) = 0.0 @@ -5783,6 +5758,11 @@ subroutine calve_icebergs(bergs) newberg%ang_vel=0.; newberg%ang_accel=0.; newberg%rot=0. endif + if (bergs%tabular_calving) then + if (.not. allocated(newberg%mask_status)) allocate(newberg%mask_status) + newberg%mask_status=grd%msk(i,j) + endif + if (.not. bergs%old_interp_flds_order) then !interpolate gridded variables to new iceberg if (grd%tidal_drift>0.) then @@ -5993,6 +5973,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%tabular_calving) then + allocate(cberg%mask_status) + cberg%mask_status=grd%msk(cberg%ine,cberg%jne) + endif + call add_new_berg_to_list(bergs%list(cberg%ine,cberg%jne)%first, cberg) end subroutine calve_fl_icebergs @@ -6492,11 +6477,7 @@ subroutine evolve_icebergs_mts(bergs) i=berg%ine ; j=berg%jne 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 + call adjust_index_and_ground(grd, lonn, latn, uveln, vveln, i, j, xi, yj, bounced, error_flag, berg%id, ignore_bounce=.true.) berg%lon=lonn ; berg%lat=latn berg%lon_old=lonn ; berg%lat_old=latn berg%ine=i ; berg%jne=j @@ -7175,7 +7156,8 @@ subroutine update_verlet_position(bergs, berg) ! Adjusting mass... !MP3 i=berg%ine; j=berg%jne; xi = berg%xi; yj = berg%yj - call adjust_index_and_ground(grd, lonn, latn, uvel3, vvel3, i, j, xi, yj, bounced, error_flag, berg%id) !Alon:"unclear which velocity to use here?" + call adjust_index_and_ground(grd, lonn, latn, uvel3, vvel3, i, j, xi, yj, bounced, error_flag, berg%id, & + ignore_bounce=bergs%tabular_calving) !Alon:"unclear which velocity to use here?" !if (bounced) then ! print *, 'you have been bounce: big time!',mpp_pe(),berg%id,lonn, latn, uvel3, vvel3, i, j, xi, yj, bounced, error_flag @@ -7244,7 +7226,7 @@ subroutine rotvec_from_tang(lon, xdot, ydot, uvel, vvel) end subroutine rotvec_from_tang !> Moves berg's cell indexes,(i,j), checking for collisional with coasts -subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, bounced, error, id) +subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, bounced, error, id, ignore_bounce) ! Arguments type(icebergs_gridded), pointer :: grd !< Container for gridded fields real, intent(inout) :: lon !< Longitude (degree E) @@ -7257,6 +7239,7 @@ subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, boun integer, intent(inout) :: j !< j-index of cell logical, intent(out) :: bounced !< True if berg collided with coast logical, intent(out) :: error !< True if adjustments could not be made consistently + logical, intent(in), optional :: ignore_bounce !< Ignore bouncing if true integer(kind=8), intent(in) :: id !< Berg identifier ! Local variables logical lret, lpos @@ -7265,6 +7248,10 @@ subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, boun real :: xi0, yj0, lon0, lat0 integer :: stderrunit logical :: point_in_cell_using_xi_yj + logical :: no_bounce + + no_bounce=.false. + if (present(ignore_bounce)) no_bounce=ignore_bounce ! Get the stderr unit number stderrunit = stderr() @@ -7370,7 +7357,7 @@ subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, boun icount=icount+1 if (xi.lt.0.) then if (i>grd%isd) then - if (grd%msk(i-1,j)>0.) then + if (no_bounce .or. grd%msk(i-1,j)>0.) then if (i>grd%isd+1) i=i-1 else !write(stderr(),'(a,6f8.3,i)') 'KID, adjust: bouncing berg from west',lon,lat,xi,yj,uvel,vvel,mpp_pe() @@ -7380,7 +7367,7 @@ subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, boun elseif (xi.ge.1.) then !Alon!!!! ! elseif (xi.gt.1.) then if (i0.) then + if (no_bounce .or. grd%msk(i+1,j)>0.) then if (igrd%jsd) then - if (grd%msk(i,j-1)>0.) then + if (no_bounce .or. grd%msk(i,j-1)>0.) then if (j>grd%jsd+1) j=j-1 else !write(stderr(),'(a,6f8.3,i)') 'KID, adjust: bouncing berg from south',lon,lat,xi,yj,uvel,vvel,mpp_pe() @@ -7400,7 +7387,7 @@ subroutine adjust_index_and_ground(grd, lon, lat, uvel, vvel, i, j, xi, yj, boun elseif (yj.ge.1.) then !Alon. ! elseif (yj.gt.1.) then if (j0.) then + if (no_bounce .or. grd%msk(i,j+1)>0.) then if (j0.) return ! Landed in ocean without bouncing so all is well + if (.not.bounced.and.lret.and.(grd%msk(i,j)>0..or.no_bounce)) return ! Landed in ocean without bouncing so all is well if (.not.bounced.and..not.lret) then ! This implies the berg traveled many cells without getting far enough if (debug) then diff --git a/src/icebergs_fms2io.F90 b/src/icebergs_fms2io.F90 index bc5fcde..2fca6b9 100644 --- a/src/icebergs_fms2io.F90 +++ b/src/icebergs_fms2io.F90 @@ -50,6 +50,8 @@ module ice_bergs_fms2io 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 +! for tabular calving from ice shelves: +use ice_bergs_framework, only: tabular_calving_global implicit none ; private @@ -168,6 +170,7 @@ subroutine write_restart_bergs(bergs, time_stamp) ang_vel, & ang_accel, & rot, & + mask_status, & tangd1, & tangd2, & nstress, & @@ -258,6 +261,9 @@ subroutine write_restart_bergs(bergs, time_stamp) allocate(ang_accel(nbergs)) allocate(rot(nbergs)) endif + if (tabular_calving_global) then + allocate(mask_status(nbergs)) + endif i = 0 do grdj = bergs%grd%jsc,bergs%grd%jec ; do grdi = bergs%grd%isc,bergs%grd%iec @@ -296,6 +302,9 @@ subroutine write_restart_bergs(bergs, time_stamp) ang_accel(i) = this%ang_accel rot(i) = this%rot endif + if (tabular_calving_global) then + mask_status(i) = this%mask_status + endif this=>this%next enddo enddo ; enddo @@ -393,6 +402,11 @@ subroutine write_restart_bergs(bergs, time_stamp) dim_names_1d,longname='dem accumulated rotation',units='rad') endif + if (tabular_calving_global) then + call register_restart_field_wrap(fileobj,'mask_status',mask_status,& + dim_names_1d,longname='tabular calving mask status',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 @@ -456,7 +470,10 @@ subroutine write_restart_bergs(bergs, time_stamp) ang_accel, & rot) endif - + if (tabular_calving_global) then + deallocate( & + mask_status) + endif deallocate( & ine, & jne, & @@ -704,7 +721,8 @@ subroutine read_restart_bergs(bergs,Time) byn_fast, & ang_vel, & ang_accel, & - rot + rot, & + mask_status integer, allocatable, dimension(:) :: ine, & jne, & @@ -809,6 +827,9 @@ subroutine read_restart_bergs(bergs,Time) allocate(ang_accel(nbergs_in_file)) allocate(rot(nbergs_in_file)) endif + if (tabular_calving_global) then + allocate(localberg%mask_status) + endif call register_restart_field(fileobj,'lon',lon,dim_names_1d) call register_restart_field(fileobj,'lat',lat,dim_names_1d) @@ -858,6 +879,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 (tabular_calving_global) then + mask_status = 0 + call register_restart_field(fileobj,'mask_status',mask_status,dim_names_1d,is_optional=.true.) + endif call read_restart(fileobj) call close_file(fileobj) elseif (bergs%require_restart) then @@ -943,6 +969,10 @@ subroutine read_restart_bergs(bergs,Time) localberg%rot =rot(k) endif + if (tabular_calving_global) then + localberg%mask_status=grd%msk(ine(k),jne(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,7 +1031,9 @@ subroutine read_restart_bergs(bergs,Time) ang_accel, & rot) endif - + if (tabular_calving_global) then + deallocate(mask_status) + endif if (replace_iceberg_num) then deallocate(iceberg_num) else @@ -1076,6 +1108,9 @@ subroutine generate_bergs(bergs,Time) allocate(localberg%ang_accel) allocate(localberg%rot) endif + if (tabular_calving_global) then + allocate(localberg%mask_status) + endif 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,7 +1160,9 @@ subroutine generate_bergs(bergs,Time) localberg%ang_accel=0. localberg%rot=0. endif - + if (tabular_calving_global) then + localberg%mask_status=grd%msk(i,j) + endif !Berg A call loc_set_berg_pos(grd, 0.9, 0.5, 1., 0., localberg) localberg%id = generate_id(grd, i, j) @@ -1641,7 +1678,7 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name integer :: uoid, void, uiid, viid, uaid, vaid, sshxid, sshyid, sstid, sssid integer :: cnid, hiid, hsid, sbid integer :: mid, smid, did, wid, lid, mbid, mflbid, mflbbid, hdid, nbid, odid, flkid -integer :: axnid,aynid,bxnid,bynid,axnfid,aynfid,bxnfid,bynfid, msid +integer :: axnid,aynid,bxnid,bynid,axnfid,aynfid,bxnfid,bynfid, msid, tmid integer :: avid, aaid, rid character(len=70) :: filename character(len=7) :: pe_name @@ -1820,6 +1857,10 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name rid = inq_varid(ncid, 'rot') endif + if (tabular_calving_global) then + tmid = inq_varid(ncid, 'mask_status') + endif + endif else ! Dimensions @@ -1891,6 +1932,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 (tabular_calving_global) then + tmid = def_var(ncid, 'mask_status', NF_DOUBLE, i_dim) + endif endif ! Attributes @@ -2010,6 +2055,11 @@ 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 (tabular_calving_global) then + call put_att(ncid, tmid, 'long_name', 'mask status') + call put_att(ncid, tmid, 'units', 'none') + endif endif endif @@ -2092,6 +2142,10 @@ subroutine write_trajectory(trajectory, save_short_traj, save_fl_traj, traj_name call put_double(ncid, rid, i, this%rot) endif + if (tabular_calving_global) then + call put_double(ncid, tmid, i, this%mask_status) + endif + endif next=>this%next deallocate(this) diff --git a/src/icebergs_framework.F90 b/src/icebergs_framework.F90 index 4d1aab4..f4f78e3 100644 --- a/src/icebergs_framework.F90 +++ b/src/icebergs_framework.F90 @@ -50,6 +50,7 @@ module ice_bergs_framework logical :: ewsame=.false. !<(F) set T if periodic and 2 PEs along the x direction (zonal) (i.e. E/W PEs are the same) logical :: iceberg_bonds_on=.False. ! True=Allow icebergs to have bonds, False=don't allow. logical :: dem=.false. !< If T, run in DEM-mode with angular terms, variable stiffness, etc +logical :: tabular_calving_global=.false. !< True to allow tabular calving of icebergs from ice shelves logical :: save_bond_forces=.true. !< Saves forces on bonds so only 1 of 2 bonds in a pair need processing during DEM-MTS explicit sub-steps logical :: short_step_mts_grounding=.false. logical :: radius_based_drag=.false. !if T, hex bergs, and dem, 2r is used as the area of the vert face for drag/wave forces @@ -69,7 +70,7 @@ module ice_bergs_framework public ignore_ij_restart, use_slow_find,generate_test_icebergs,old_bug_rotated_weights,budget public orig_read, force_all_pes_traj public mts,save_bond_traj,ewsame,iceberg_bonds_on -public dem, save_bond_forces, orig_dem_moment_of_inertia +public dem, save_bond_forces, orig_dem_moment_of_inertia, tabular_calving_global public short_step_mts_grounding, radius_based_drag public A68_test, A68_xdisp, A68_ydisp public footloose @@ -298,6 +299,8 @@ module ice_bergs_framework real, allocatable :: ang_vel !< Angular velocity real, allocatable :: ang_accel !< Angular acceleration real, allocatable :: rot !< Accumulated rotation + ! For tabular calving from ice shelves + real, allocatable :: mask_status end type xyt !> An iceberg object, used as a link in a linked list @@ -372,6 +375,8 @@ module ice_bergs_framework real, allocatable :: ang_vel !< Angular velocity real, allocatable :: ang_accel !< Angular acceleration real, allocatable :: rot !< Accumulated rotation + ! For tabular calving from ice shelves + real, allocatable :: mask_status end type iceberg !> A bond object connecting two bergs, used as a link in a linked list @@ -638,7 +643,7 @@ module ice_bergs_framework logical :: snap_tabular_calving_to_bonded_grid=.true. !align tabular particles that calve from ice shelves with a constant cartesian grid !backwards compatibility logical :: old_interp_flds_order=.false. !< Use old order of when to interpolate grid variables to bergs. Will be false if MTS, DEM, or footloose - logical :: tabular_calving=.false. + logical :: tabular_calving=.false. !< True to allow tabular calving of icebergs from ice shelves real :: shelf_to_tabular_hours !< Time (hours) over which ice shelf is transitioned to bonded-particle tabular icebergs end type icebergs @@ -864,6 +869,7 @@ subroutine ice_bergs_framework_init(bergs, & real :: constant_length=0. ! If constant_interaction_LW, the constant length used. If zero in the nml, will be set to max initial L real :: constant_width=0. ! If constant_interaction_LW, the constant width used. If zero in the nml, will be set to max initial W real :: ocean_drag_scale=1. !< Scaling factor for the ocean drag coefficients +logical :: calculate_spring_from_dem_spring !< If constant_interaction_LW, will calculate spring_coef that corresponds to dem_spring_coef ! 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. @@ -907,7 +913,7 @@ subroutine ice_bergs_framework_init(bergs, & fl_youngs, fl_strength, save_all_traj_year, save_nonfl_traj_by_class,& save_traj_by_class_start_mass_thres_n, save_traj_by_class_start_mass_thres_s,traj_area_thres_sntbc,& traj_area_thres_fl,tau_is_velocity, ocean_drag_scale, A68_test, & - A68_xdisp,A68_ydisp,use_broken_bonds_for_substep_contact,print_fracture,& + A68_xdisp,A68_ydisp,use_broken_bonds_for_substep_contact,print_fracture,calculate_spring_from_dem_spring,& orig_dem_moment_of_inertia, break_bonds_on_sub_steps, skip_first_outer_mts_step, rev_mind, & no_frac_first_ts, use_grounding_torque, short_step_mts_grounding, radius_based_drag, save_bond_forces, shelf_to_tabular_hours @@ -1520,6 +1526,9 @@ subroutine ice_bergs_framework_init(bergs, & bergs%constant_radius=sqrt(bergs%constant_area/pi) ! Interaction radius of the iceberg (assuming circular icebergs) endif endif + if (calculate_spring_from_dem_spring .and. bergs%dem .and. bergs%constant_radius.ne.0) then + bergs%spring_coef = bergs%dem_spring_coef/(bergs%rho_bergs * 4.*(bergs%constant_radius**2)) + endif endif bergs%constant_radius_IS_berg=constant_radius_IS_berg bergs%snap_tabular_calving_to_bonded_grid=snap_tabular_calving_to_bonded_grid @@ -1585,6 +1594,11 @@ subroutine ice_bergs_framework_init(bergs, & else bergs%tabular_calving=.false. endif + if (bergs%tabular_calving) then + buffer_width=buffer_width+1 + buffer_width_traj=buffer_width_traj+1 + tabular_calving_global=.true. + endif !necessary? if (.not. mts) then if ((halo-1)berg%first_bond do k = 1,max_bonds @@ -3960,6 +3978,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 (tabular_calving_global) allocate(localberg%mask_status) + counter = 0 call pull_buffer_value(buff%data(:,n), counter, localberg%lon) call pull_buffer_value(buff%data(:,n), counter, localberg%lat) @@ -4027,6 +4047,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 (tabular_calving_global) then + call pull_buffer_value(buff%data(:,n), counter, localberg%mask_status) + endif + !These quantities no longer need to be passed between processors localberg%uvel_old=localberg%uvel localberg%vvel_old=localberg%vvel @@ -4285,6 +4309,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 (tabular_calving_global) then + call push_buffer_value(buff%data(:,n), counter, traj%mask_status) + endif endif end subroutine pack_traj_into_buffer2 @@ -4312,6 +4340,8 @@ subroutine unpack_traj_from_buffer2(first, buff, n, save_short_traj, save_fl_tra if (dem) allocate(traj%ang_vel,traj%ang_accel,traj%rot) + if (tabular_calving_global) allocate(traj%mask_status) + counter = 0 call pull_buffer_value(buff%data(:,n),counter,traj%lon) call pull_buffer_value(buff%data(:,n),counter,traj%lat) @@ -4378,6 +4408,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 (tabular_calving_global) then + call pull_buffer_value(buff%data(:,n), counter, traj%mask_status) + endif endif call append_posn(first, traj) @@ -4966,6 +5000,8 @@ subroutine create_iceberg(berg, bergvals) if (dem) allocate(berg%ang_vel,berg%ang_accel,berg%rot) + if (tabular_calving_global) allocate(berg%mask_status) + berg=bergvals berg%prev=>null() berg%next=>null() @@ -5945,6 +5981,8 @@ subroutine record_posn(bergs) if (dem) allocate(posn%ang_vel,posn%ang_accel,posn%rot) + if (tabular_calving_global) allocate(posn%mask_status) + 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) @@ -6036,6 +6074,10 @@ subroutine record_posn(bergs) posn%ang_accel=this%ang_accel posn%rot=this%rot endif + + if (tabular_calving_global) then + posn%mask_status=this%mask_status + endif endif call push_posn(this%trajectory, posn) @@ -6104,6 +6146,8 @@ subroutine push_posn(trajectory, posn_vals) if (dem) allocate(new_posn%ang_vel,new_posn%ang_accel,new_posn%rot) + if (tabular_calving_global) allocate(new_posn%mask_status) + new_posn=posn_vals new_posn%next=>trajectory trajectory=>new_posn @@ -6150,6 +6194,8 @@ subroutine append_posn(trajectory, posn_vals) if (dem) allocate(new_posn%ang_vel,new_posn%ang_accel,new_posn%rot) + if (tabular_calving_global) allocate(new_posn%mask_status) + new_posn=posn_vals new_posn%next=>null() if(.NOT. associated(trajectory)) then