Skip to content

Commit

Permalink
added test dem_test_contact, for varying L_ij, the critical-interacti…
Browse files Browse the repository at this point in the history
…ve-length scale, for MTS-DEM simulations. Artifically increasing L_ij is needed to prevent overlapping of conglomerates under contact if their particle radii are low. Also added long_run_no_inc_contact_dist.sh and .nml in tests/a68, to test A68a test without increaing L_ij (it has essentially no effect, given the large particle radii). Some edits to icebergs.F90 to increase computational efficiency.
  • Loading branch information
alex-huth committed Mar 9, 2022
1 parent 1a02384 commit 2686245
Show file tree
Hide file tree
Showing 15 changed files with 928 additions and 75 deletions.
38 changes: 34 additions & 4 deletions driver/icebergs_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ program icebergs_driver
use mpp_mod, only : mpp_pe
use mpp_mod, only : mpp_root_pe
use mpp_mod, only: mpp_sync
use mpp_mod, only: mpp_exit
use time_manager_mod, only : time_type
use time_manager_mod, only : set_date
use time_manager_mod, only : set_calendar_type
Expand Down Expand Up @@ -60,6 +61,7 @@ program icebergs_driver
logical :: chaotic_test=.false.
logical :: grounding_test=.false.
logical :: big_grounding_test=.false.
logical :: dem_test_contact=.false.
logical :: a68_test=.false.
logical :: vert_int_ocean_vel=.false.
logical :: tau_is_velocity=.true.
Expand Down Expand Up @@ -87,7 +89,7 @@ program icebergs_driver
saverestart,ibui,ibvi,collision_test,chaotic_test,grounding_test,&
gridres,write_time_inc,bump_depth, sst, a68_test, data_dir, vert_int_ocean_vel,&
reverse_a68_forcings,tau_is_velocity,transient_a68_data_start_ind,&
Old_wrong_Rearth_and_SSH,no_wind,REarth,big_grounding_test !mom_Rearth
Old_wrong_Rearth_and_SSH,no_wind,REarth,big_grounding_test,dem_test_contact !mom_Rearth
! For loops
integer :: isc !< Start of i-index for computational domain (used for loops)
integer :: iec !< End of i-index for computational domain (used for loops)
Expand Down Expand Up @@ -334,6 +336,32 @@ program icebergs_driver
enddo; enddo
print *,'min, max depth',minval(depth(:,:)),maxval(depth(:,:))
depth = 1000.0-depth
endif

if (dem_test_contact) then
lat(:,:)=lat(:,:)-0.45 !-25000.1
lon(:,:)=lon(:,:)-0.45
!assign land cells at the N and S portions of the domain.
!input.nml: set coastal drift parameter to prevent grounding
do j=jsd,jed; do i=isd,ied
if (lat(i,j)<=(-5.e3) .or. lat(i,j)>=(220.e3)) wet(i,j)=0.0
enddo;enddo

mid=50.e3

do j=jsd,jed; do i=isd,ied
if (lon(i,j)>65 .or. lat(i,j)==mid) then
vo(i,j)=0.0
else
if (lat(i,j)>mid) then
vo(i,j)=-ibvo
else
vo(i,j)=ibvo
endif
endif
enddo;enddo


endif

if (grounding_test) then
Expand Down Expand Up @@ -513,13 +541,15 @@ program icebergs_driver

if (saverestart) call icebergs_save_restart(bergs)

!Deallocate all memory and disassociated pointer
call icebergs_end(bergs)

call cpu_time(time_finish)

if (mpp_pe()==0) then
write(*,*) 'Simulation finished in ', (time_finish - time_begin)/60., 'minutes'
end if

call mpp_exit()

!Deallocate all memory and disassociated pointer
call icebergs_end(bergs)

end program icebergs_driver
76 changes: 42 additions & 34 deletions src/icebergs.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1176,6 +1176,7 @@ subroutine calculate_force_dem(bergs, berg, other_berg, current_bond, &

if (current_bond%nstress<0) then
!if fail under shear with normal compression, then should still feel contact force
damping_coef = bergs%dem_damping_coef*sqrt(bergs%dem_K_damp*M1*M2/(M1+M2))
Fd_x = Fd_x-damping_coef*ur; Fd_y = Fd_y-damping_coef*vr !linear damping force
F_x = F_x + Fn_x; F_y = F_y + Fn_y !linear forces
endif
Expand Down Expand Up @@ -2089,7 +2090,7 @@ subroutine accel_explicit_inner_mts(bergs, berg, i, j, xi, yj, lat, uvel, vvel,
do grdi = max(berg%ine-1,bergs%grd%isd+1),min(berg%ine+1,bergs%grd%ied)
other_berg=>bergs%list(grdi,grdj)%first
do while (associated(other_berg))
if (other_berg%id>0 .and. other_berg%conglom_id.eq.berg%conglom_id) then
if (other_berg%id>0 .and. other_berg%conglom_id.eq.berg%conglom_id .and. other_berg%n_bonds<bergs%max_bonds) then
if (bergs%dem) then
call calculate_unbonded_same_conglom_dem_force(bergs, berg, other_berg, &
IA_x, IA_y, IAd_x, IAd_y, uvel0, vvel0, uvel0, vvel0)
Expand Down Expand Up @@ -3392,12 +3393,13 @@ subroutine thermodynamics(bergs)
if (bergs%use_mixed_melting .or. bergs%allow_bergs_to_roll) then
N_bonds=0.
if (bergs%iceberg_bonds_on) then
N_bonds=this%n_bonds
! Determining number of bonds
current_bond=>this%first_bond
do while (associated(current_bond)) ! loop over all bonds
N_bonds=N_bonds+1.0
current_bond=>current_bond%next_bond
enddo
! current_bond=>this%first_bond
! do while (associated(current_bond)) ! loop over all bonds
! N_bonds=N_bonds+1.0
! current_bond=>current_bond%next_bond
! enddo
endif
if (this%static_berg .eq. 1) N_bonds=N_max !Static icebergs melt like ice shelves
endif
Expand Down Expand Up @@ -7303,35 +7305,40 @@ subroutine evolve_icebergs_mts(bergs)
do while (associated(berg)) ! loop over all bergs
!only evolve non-static bergs that overlap, or are part of a conglom that overlaps, the computational domain:
if (berg%static_berg .lt. 0.5 .and. (berg%conglom_id.ne.0 .or. bergs%force_convergence)) then
latn = berg%lat ; lonn = berg%lon
uvel1=berg%uvel ; vvel1=berg%vvel !V_n (previous cycle)
axn = 0.0 ; ayn = 0.0 !note these are redefined at the start of subroutine accel_mts
bxn = 0.0 ; byn = 0.0
i=berg%ine ; j=berg%jne
xi=berg%xi ; yj=berg%yj

call accel_mts(bergs, berg, i, j, xi, yj, latn, uvel1, vvel1, uvel1, vvel1, dt, rx, ry, &
ax1, ay1, axn, ayn, bxn, byn, save_bond_energy,Fec_x, Fec_y, Fdc_x, Fdc_y)

if (Fdc_x .ne. 0. .or. Fdc_y .ne. 0.) had_collision=.true.
if (monitor_energy .and. last_iter) &
call mts_energy_part_1(bergs,berg,ax1,ay1,axn,ayn,bxn,byn,Fec_x,Fec_y,Fdc_x,Fdc_y)

! Saving all the iceberg variables
berg%axn=axn; berg%ayn=ayn
berg%bxn=bxn; berg%byn=byn

if (bergs%force_convergence) then
berg%uvel_prev=berg%uvel+(dt*ax1); berg%vvel_prev=berg%vvel+(dt*ay1) !the new velocity
if (ii==1) usum=usum+berg%uvel_old**2 + berg%vvel_old**2
usum1=usum1+berg%uvel_prev**2+berg%vvel_prev**2
usum2=usum2+(berg%uvel_prev-berg%uvel_old)**2+(berg%vvel_prev-berg%vvel_old)**2
else
berg%uvel=berg%uvel+(dt*ax1); berg%vvel=berg%vvel+(dt*ay1)
if (ii==1 .or. berg%static_berg==0.1) then
latn = berg%lat ; lonn = berg%lon
uvel1=berg%uvel ; vvel1=berg%vvel !V_n (previous cycle)
axn = 0.0 ; ayn = 0.0 !note these are redefined at the start of subroutine accel_mts
bxn = 0.0 ; byn = 0.0
i=berg%ine ; j=berg%jne
xi=berg%xi ; yj=berg%yj

call accel_mts(bergs, berg, i, j, xi, yj, latn, uvel1, vvel1, uvel1, vvel1, dt, rx, ry, &
ax1, ay1, axn, ayn, bxn, byn, save_bond_energy,Fec_x, Fec_y, Fdc_x, Fdc_y)

if (Fdc_x .ne. 0. .or. Fdc_y .ne. 0. .and. bergs%force_convergence) then
had_collision=.true.
berg%static_berg=0.1
endif
if (monitor_energy .and. last_iter) &
call mts_energy_part_1(bergs,berg,ax1,ay1,axn,ayn,bxn,byn,Fec_x,Fec_y,Fdc_x,Fdc_y)

! Saving all the iceberg variables
berg%axn=axn; berg%ayn=ayn
berg%bxn=bxn; berg%byn=byn

if (bergs%force_convergence) then
berg%uvel_prev=berg%uvel+(dt*ax1); berg%vvel_prev=berg%vvel+(dt*ay1) !the new velocity
if (ii==1) usum=usum+berg%uvel_old**2 + berg%vvel_old**2
usum1=usum1+berg%uvel_prev**2+berg%vvel_prev**2
usum2=usum2+(berg%uvel_prev-berg%uvel_old)**2+(berg%vvel_prev-berg%vvel_old)**2
else
berg%uvel=berg%uvel+(dt*ax1); berg%vvel=berg%vvel+(dt*ay1)

!The final velocity from the previous cycle, which can be used
!during post-processing to calculate kinetic energy and momentum
berg%uvel_prev=berg%uvel ; berg%vvel_prev=berg%vvel
!The final velocity from the previous cycle, which can be used
!during post-processing to calculate kinetic energy and momentum
berg%uvel_prev=berg%uvel ; berg%vvel_prev=berg%vvel
endif
endif
endif
berg=>berg%next
Expand All @@ -7345,6 +7352,7 @@ subroutine evolve_icebergs_mts(bergs)
do while (associated(berg)) ! loop over all bergs
if (berg%static_berg .lt. 0.5) then
berg%uvel_old=berg%uvel_prev; berg%vvel_old=berg%vvel_prev
if (last_iter .and. berg%static_berg==0.1) berg%static_berg=0.
endif
berg=>berg%next
enddo
Expand Down
14 changes: 7 additions & 7 deletions tests/a68/long_run.nml
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
&icebergs_nml

Rearth=<re>
!6363827 !6.378e6 !6.371229e6 !6.36e6 !6.378e6

A68_test=.true.
A68_xdisp=<xd>
A68_ydisp=<yd>
Expand All @@ -40,9 +40,9 @@

constant_interaction_LW=.true.

traj_name=<name>_1p.nc
bond_traj_name=<name>_bonds_1p.nc
traj_name=<name>.nc
bond_traj_name=<name>_bonds.nc

break_bonds_on_sub_steps=.true.

power_ground=.false.
Expand Down Expand Up @@ -121,7 +121,7 @@
coastal_drift=0. !(0.)Vel added to currents to cause bergs to drift away from land cells
tidal_drift=0. !(0.)Amp of stoch. tidal vel added to currents to cause random berg drift
Runge_not_Verlet=.false. !(T) True=Runge Kutta, False=Verlet.
use_mixed_melting=.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #)
use_mixed_melting=.true. !.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #)
apply_thickness_cutoff_to_gridded_melt=.true.!(F)Prevnt melt for Ho < melt_cuttoff on grid
apply_thickness_cutoff_to_bergs_melt=.true. !(F)Prevnt melt for Ho < melt_cuttoff on bergs
melt_cutoff=10.0 ! (-1.0) Min ocean thickness for melting to occur (no apply for vals < 0)
Expand Down Expand Up @@ -153,7 +153,7 @@
verbose_hrs=1!24000 ! (24) Period between verbose messages
traj_sample_hrs=0!1 !(24) sampling period for trajectory storage
traj_write_hrs=0!6 !(480) Period for writing sampled trajectories to disk
save_bond_traj=.true. !.false. !(F) If T, saves bond trajectories in a separate file
save_bond_traj=.false. !(F) If T, saves bond trajectories in a separate file
budget=.false. !(T) Calculate budgets
debug=.false. !(F) Turn on debugging
really_debug=.false. !(F) Turn on debugging
Expand Down Expand Up @@ -211,7 +211,7 @@
/

&fms_nml
clock_grain='NONE'
clock_grain='MODULE'
domains_stack_size = 2000000
clock_flags='SYNC'
/
48 changes: 22 additions & 26 deletions tests/a68/input.nml → tests/a68/long_run_no_inc_contact_dist.nml
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,11 @@
data_dir='/lustre/f2/dev/gfdl/Alexander.Huth/MOM6-examples/src/icebergs/tests/a68/data/'
Old_wrong_Rearth_and_SSH=.false.
!mom_Rearth=.true.
Rearth=6363827
Rearth=<re>
!6363827 !6.378e6
no_wind=.false.
tau_is_velocity=.true.!.false.
ibdt=1800.0 !seconds
ibdt=<ibdt> !3600.0 !1800.0 !seconds
ibuo=0.0 !ocean x velocity
ibvo=0.0 !ocean y velocity
ibui=0.0 !ice x velocity
Expand All @@ -25,24 +25,24 @@

&icebergs_nml

Rearth=6363827
!6363827 !6.378e6 !6.371229e6 !6.36e6 !6.378e6
Rearth=<re>

A68_test=.true.
A68_xdisp=-37.51
A68_ydisp=-55.2166
A68_xdisp=<xd>
A68_ydisp=<yd>
short_step_mts_grounding=.true.

tau_is_velocity=.true.!.false.
ocean_drag_scale=17.8 !20. !1
tau_is_velocity=.true.
ocean_drag_scale=<od>

skip_first_outer_mts_step=.true.
add_curl_to_torque=.false.

constant_interaction_LW=.true.

traj_name=longrun_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166_1p.nc
bond_traj_name=longrun_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166_bonds_1p.nc
traj_name=<name>.nc
bond_traj_name=<name>_bonds.nc

break_bonds_on_sub_steps=.true.

power_ground=.false.
Expand All @@ -68,27 +68,23 @@
!defaults given in parentheses:
new_mts=.true.
mts=.true. !(F) for using the Multiple Time Stepping Velocity Verlet scheme
mts_sub_steps=90 !200 !75 !75 !150 !50 !(-1) # of mts sub-steps
mts_sub_steps=<ss> !200 !75 !75 !150 !50 !(-1) # of mts sub-steps
!mts_sub_steps=10
force_convergence=.true. !(F) experimental MTS convergence scheme to better conserve momentum during collisions
convergence_tolerance=1e-4 !8 !(1.e-8) tolerance for the MTS force_convergence scheme

contact_distance=4.e3 !(0) unbonded bergs contact at max(contact_dist,sum of their radii)
contact_distance=0.!4.e3 !(0) unbonded bergs contact at max(contact_dist,sum of their radii)
contact_spring_coef=1.e-7!(default=spring_coeff) Berg spring coeff for collisions
debug_write=.false. !(F)sets traj_sample_hours=traj_write_hours & includes halo bergs in write
cdrag_grounding=1.e4 !0.!1.e6!3.16e6!1.e6!1.e3!3.16e6 !for power-law-type friction
cdrag_grounding=<gc> !0.!1.e6!3.16e6!1.e6!1.e3!3.16e6 !for power-law-type friction
h_to_init_grounding=0!30!15!50. !15.
ewsame=.false. !(F) set T if periodic and 2 PEs along the x (zonal) direction
remove_unused_bergs=.true. !(T) remove unneeded bergs after PEs transfers

fracture_criterion='stress' !('none') 'strain_rate','strain',or 'none'
frac_thres_scaling = 1 !0000000
frac_thres_n = 18e3 !25000.0 !5000.0 !4250. !12000.0 !8000.0 !5575. !3000. !2.1e4 !12.
frac_thres_t = 100e3 !1000. !2000.0 !2500. !6000.0 !4000.0 !5575. !1750. !35000. !12.

!!fracture_criterion='strain'
!!frac_thres_n = 0.0162
!!frac_thres_t = 0.0457
frac_thres_scaling = 1
frac_thres_n = <ns>e3
frac_thres_t = 1000e3

halo=3 !(4) halo width
Lx=360 !28.875!45000. ! (360.) Length of X-domain. For periodicity (use a huge number for non-periodic)
Expand Down Expand Up @@ -125,7 +121,7 @@
coastal_drift=0. !(0.)Vel added to currents to cause bergs to drift away from land cells
tidal_drift=0. !(0.)Amp of stoch. tidal vel added to currents to cause random berg drift
Runge_not_Verlet=.false. !(T) True=Runge Kutta, False=Verlet.
use_mixed_melting=.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #)
use_mixed_melting=.true. !.false.!(F) 'T' partly determ. melt frm 3eq model(uses BergParams&Bond #)
apply_thickness_cutoff_to_gridded_melt=.true.!(F)Prevnt melt for Ho < melt_cuttoff on grid
apply_thickness_cutoff_to_bergs_melt=.true. !(F)Prevnt melt for Ho < melt_cuttoff on bergs
melt_cutoff=10.0 ! (-1.0) Min ocean thickness for melting to occur (no apply for vals < 0)
Expand Down Expand Up @@ -155,9 +151,9 @@

verbose=.false. !(F) Be verbose to stderr
verbose_hrs=1!24000 ! (24) Period between verbose messages
traj_sample_hrs=0!1 !(24) sampling period for trajectory storage
traj_write_hrs=0!6 !(480) Period for writing sampled trajectories to disk
save_bond_traj=.true. !.false. !(F) If T, saves bond trajectories in a separate file
traj_sample_hrs=1 !0!1 !(24) sampling period for trajectory storage
traj_write_hrs=1 !0!6 !(480) Period for writing sampled trajectories to disk
save_bond_traj=.false. !(F) If T, saves bond trajectories in a separate file
budget=.false. !(T) Calculate budgets
debug=.false. !(F) Turn on debugging
really_debug=.false. !(F) Turn on debugging
Expand Down Expand Up @@ -215,7 +211,7 @@
/

&fms_nml
clock_grain='NONE'
clock_grain='MODULE'
domains_stack_size = 2000000
clock_flags='SYNC'
/
40 changes: 40 additions & 0 deletions tests/a68/long_run_no_inc_contact_dist.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
#!/bin/bash

Rearth=6363827
odrag=(17.8)
ss=90
name=(lrun_nicd_od_17.8_tn_ng_ydn.1_ss90_ns18_rx-37.51_ry-55.2166)
ibdt=1800.0

gcoef=(1.e4)
nstress=(18)
A68dx=-37.51
A68dy=-55.2166

odrag=(10.35)
ss=180
name=(lrun_nicd_hrtest)
ibdt=3600.0
nstress=(24.) #(22.5)
A68dx=-37.5
A68dy=-55.1

for ((i = 0; i < ${#odrag[@]}; ++i)); do
drag=(${odrag[$i]})
fname=(${name[$i]})
gc=(${gcoef[$i]})
ns=(${nstress[$i]})

echo $drag
echo $fname
echo $Rearth
echo $ss
echo $A68dx
echo $A68dy
echo $ns

sed "s/<od>/$drag/; s/<re>/$Rearth/; s/<ss>/$ss/; s/<ns>/$ns/; s/<gc>/$gc/; s/<xd>/$A68dx/; s/<yd>/$A68dy/; s/<ibdt>/$ibdt/; s/<name>/$fname/g" \
long_run_no_inc_contact_dist.nml > input.nml
./RUN1p
# &>/dev/null &
done
Loading

0 comments on commit 2686245

Please sign in to comment.