use "tavg.01981.01.01.nc"


let plastic = o_plastic[d=1] !particles/m2, including the missing plastic
! convert to number of aggregates that are full
let nagg_phyt1 = plastic/1 !assume 1 to be conservative !number of full organic aggregates/m2
! find total number of aggregates leaving surface ocean per day
let tot_nagg = (o_detrexp[d=1,k=1]+o_detrzexp[d=1,k=1])*G_dzt[k=1]*0.03*60*60*6.625*12.011/8.8e-6 !#aggs exported/m2/hour
!find ratio of full organic aggs to total organic aggs
let clean_out = nagg_phyt1/tot_nagg ! number of hours to clean out all the plastic

go "figure_scripts/prepareFigure.jnl" number_hours_to_removal
go fland
go log_key 1, 10 "clean_out[k=1,l=31]"
go fland

go "figure_scripts/finalizeFigure_psthicken.jnl"

let clean_out_rate= plastic/clean_out*24*365 !particles removed m2/y

let input = plastic[x=@din,y=@din,l=@ddc*60*60*24*365] + clean_out_rate[x=@din,y=@din]

!convert to tons of plastic
let erik = 355.4e8/485e10 ! this conversion factor produces much heavier plastic!


go "figure_scripts/prepareFigure.jnl" plastic_removal_2100
set view ul
go fland
fill/pal=blue_darkred/nolab/ov/key/lev=30h/set clean_out[k=1,l=31]
ppl shakey 1, 1, 0.12
ppl shade
go fland

set view ur !units of particles removed/m2/y
go fland
fill/pal=blue_darkred/nolab/ov/key/lev=30c/set clean_out_rate[k=1,l=120]-clean_out_rate[k=1,l=31]
ppl shakey 1, 1, 0.12, 5 
ppl shade
go fland

set view ll
pl/nolab/thick=2 clean_out_rate[x=@ave,y=@ave]

set view lr
pl/nolab/thick=2 input*erik*1e-6 !tons of MP per year 
go "figure_scripts/finalizeFigure_psthicken.jnl"

set win 1
sha clean_out_rate[k=1]*(o_velz[k=1,g=o_temp@mod]+o_gmvelz[k=1,g=o_temp@mod])*60*60*24*365