!Calculate rates of plastic deposition on seafloor

use tavg.01981.01.01.nc 

set memory/size=100

!total export in mol N/m3/s

let det_exp = o_detrexp
let detz_exp = o_detrzexp

!only count what hits the bottom
! Phytoplankton detritus first
let export2 = IF G_kmt EQ 2 THEN det_exp[k=2] else 0
let export3 = IF G_kmt EQ 3 THEN det_exp[k=3] else 0
let export4 = IF G_kmt EQ 4 THEN det_exp[k=4] else 0
let export5 = IF G_kmt EQ 5 THEN det_exp[k=5] else 0
let export6 = IF G_kmt EQ 6 THEN det_exp[k=6] else 0
let export7 = IF G_kmt EQ 7 THEN det_exp[k=7] else 0
let export8 = IF G_kmt EQ 8 THEN det_exp[k=8] else 0
let export9 = IF G_kmt EQ 9 THEN det_exp[k=9] else 0
let export10 = IF G_kmt EQ 10 THEN det_exp[k=10] else 0
let export11 = IF G_kmt EQ 11 THEN det_exp[k=11] else 0
let export12 = IF G_kmt EQ 12 THEN det_exp[k=12] else 0
let export13 = IF G_kmt EQ 13 THEN det_exp[k=13] else 0
let export14 = IF G_kmt EQ 14 THEN det_exp[k=14] else 0
let export15 = IF G_kmt EQ 15 THEN det_exp[k=15] else 0
let export16 = IF G_kmt EQ 16 THEN det_exp[k=16] else 0
let export17 = IF G_kmt EQ 17 THEN det_exp[k=17] else 0
let export18 = IF G_kmt EQ 18 THEN det_exp[k=18] else 0
let export19 = IF G_kmt EQ 19 THEN det_exp[k=19] else 0

let tot_exp = export2+export3+export4+export5+export6+export7+export8+export9+export10+export11+export12+export13+export14+export15+export16+export17+export18+export19

! Now zooplankton poo
let exportz2 = IF G_kmt EQ 2 THEN detz_exp[k=2] else 0
let exportz3 = IF G_kmt EQ 3 THEN detz_exp[k=3] else 0
let exportz4 = IF G_kmt EQ 4 THEN detz_exp[k=4] else 0
let exportz5 = IF G_kmt EQ 5 THEN detz_exp[k=5] else 0
let exportz6 = IF G_kmt EQ 6 THEN detz_exp[k=6] else 0
let exportz7 = IF G_kmt EQ 7 THEN detz_exp[k=7] else 0
let exportz8 = IF G_kmt EQ 8 THEN detz_exp[k=8] else 0
let exportz9 = IF G_kmt EQ 9 THEN detz_exp[k=9] else 0
let exportz10 = IF G_kmt EQ 10 THEN detz_exp[k=10] else 0
let exportz11 = IF G_kmt EQ 11 THEN detz_exp[k=11] else 0
let exportz12 = IF G_kmt EQ 12 THEN detz_exp[k=12] else 0
let exportz13 = IF G_kmt EQ 13 THEN detz_exp[k=13] else 0
let exportz14 = IF G_kmt EQ 14 THEN detz_exp[k=14] else 0
let exportz15 = IF G_kmt EQ 15 THEN detz_exp[k=15] else 0
let exportz16 = IF G_kmt EQ 16 THEN detz_exp[k=16] else 0
let exportz17 = IF G_kmt EQ 17 THEN detz_exp[k=17] else 0
let exportz18 = IF G_kmt EQ 18 THEN detz_exp[k=18] else 0
let exportz19 = IF G_kmt EQ 19 THEN detz_exp[k=19] else 0

let tot_expz = exportz2+exportz3+exportz4+exportz5+exportz6+exportz7+exportz8+exportz9+exportz10+exportz11+exportz12+exportz13+exportz14+exportz15+exportz16+exportz17+exportz18+exportz19

let plastic = o_plastic !particles/m3 in the surface box, not including the missing plastic
! Kooi et al. 2017 model oscillations between 0-80 m depth so I assume uniform concentration through the surface box (50 meter depth)

! convert to number of aggregates that are full
let nagg_phyt1 = plastic/20 !1 !assume 1 to be conservative !number of full organic aggregates/m3

! find total number of aggregates leaving the ocean per day
! count export reaching the bottom only
! Assume 3% aggregation rate, 100% egestion is poo pellets
! aggregate and pellet to gC conversion factors: 8.8, 1.27e-6 and 1.7e-5 
let tot_nagg = (tot_exp[k=1])*0.03*60*60*24*6.625*12.011/8.8e-6 !diatom #aggs removed/m3/day
let tot_naggz = (tot_expz[k=1])*60*60*24*6.626*12.011/.000017 !#aggs removed/m3/day

! Or count all export to look at water column or surface
!let tot_nagg = (o_detrexp)*0.03*60*60*24*6.625*12.011/8.8e-6 !diatom #aggs removed/m3/day
!let tot_naggz = (o_detrzexp)*60*60*24*6.626*12.011/.000017 !#aggs removed/m3/day

list tot_nagg[x=@din,y=@din,k=1]*365 !really from the surface, not the bottom of k=1
list tot_naggz[x=@din,y=@din,k=1]*365

!find ratio of full organic aggs to total organic aggs
let clean_out = nagg_phyt1/(tot_nagg+tot_naggz) ! number of days to clean out all the plastic
!let clean_out = nagg_phyt1/(tot_nagg)
!let clean_out = nagg_phyt1/(tot_naggz)
sha clean_out[k=1,l=31,y=40S:40N]

go "figure_scripts/prepareFigure.jnl" number_days_to_removal_2010
go fland
go log_key 1, 1000 "clean_out[k=1,l=31]"
go fland
go "figure_scripts/finalizeFigure_psthicken.jnl"

let clean_out_rate= plastic/clean_out[k=1]/365 !particles removed m3/y

!convert to tons of plastic
let erik = 355.4e8/485e10 ! g/particles this conversion factor produces much heavier plastic!
let sebil = 236000/51.2e12 ! tons/particles

list  clean_out_rate[x=@din,y=@din,l=1:31@sum]*erik*1e-6 !plastic removed between 1980-2010
list  clean_out_rate[x=@din,y=@din,l=1:31@sum]*sebil

go "figure_scripts/prepareFigure.jnl" mp_removed_tons_per_year
pl/nolab/thick=2 clean_out_rate[x=@din,y=@din,k=1]*erik*1e-6 !tons of MP per year 
go "figure_scripts/finalizeFigure_psthicken.jnl"
list  clean_out_rate[x=@din,y=@din]*sebil !*erik*1e-6 
list  clean_out_rate[x=@din,y=@din]*erik*1e-6 

go "figure_scripts/prepareFigure.jnl" mp_removed_to_sediments_map_2010
sha/nolab/pal=yellow_green_blue clean_out_rate[k=1,l=1:31@sum] !particles MP seafloor 1980-2010/m2
go fland
go "figure_scripts/finalizeFigure_psthicken.jnl"

go "figure_scripts/prepareFigure.jnl" del_number_years_to_removal
sha/nolab/pal=blue_darkred/lev=(-inf)(-10,10,1)(inf) (clean_out[k=1,l=121]-clean_out[k=1,l=31])/365
go fland
go "figure_scripts/finalizeFigure_psthicken.jnl"


let tot_nagg = (o_detrexp)*0.03*60*60*24*6.625*12.011/8.8e-6 !diatom #aggs removed/m3/day
let tot_naggz = (o_detrzexp)*60*60*24*6.626*12.011/.000017 !#aggs removed/m3/day

let export = (tot_nagg+tot_naggz) !particles exported/day at each level
let min_exp = 0.25*export[k=1]
let zero_at = export - min_exp
 
go "figure_scripts/prepareFigure.jnl" export_depth_2010
sha/nolab/pal=yellow_green_blue/lev=(0,1000,50)(inf) zero_at[l=31,z=@LOC:0]
con/ov/nolab/lev=(10000,10000) plastic[l=31]*1e6
go fland
go "figure_scripts/finalizeFigure_psthicken.jnl"