! Script to plot fluxes and overlay trap data on it

set mem/size=100

! ! ! ! !  CONVERSION FACTORS  ! ! ! ! ! ! ! 
! convert from per second to per day
let mult = 60*60*24

! convert from per day to per year
let mult1 = mult*365

! convert nitrogen to carbon
let ntc = 6.625

! convert mol C to grams C
let mtg = 12.011
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! model data
use tavg.01775.01.01.nc

let model_poc=(o_detrexp[k=11]+o_detrexp_B[k=11])*mult1*ntc*1000*350 ! mmol C/y
let model_pic=o_caco3exp[k=11]*mult1*1000*350 ! mmol C/y
let model_opl=o_oplexp[k=11]*mult1*1000*350 ! mmol Si/y

list (model_poc[x=@din,y=@din,l=53]-model_poc[x=@din,y=@din,l=25])/model_poc[x=@din,y=@din,l=25]
list (model_pic[x=@din,y=@din,l=53]-model_pic[x=@din,y=@din,l=25])/model_pic[x=@din,y=@din,l=25]
list (model_opl[x=@din,y=@din,l=53]-model_opl[x=@din,y=@din,l=25])/model_opl[x=@din,y=@din,l=25]

go prepareFigure fluxes_2300
go "$HOME/opt_postprocessing_scripts/go/ilandscape3x1.jnl"
ppl axlsze,0.2,0.2
ppl axlint 4 2
set view v1
shade/set/pal=blue_darkred/levels=(-inf)(-100,100,25)(inf)/nolab (model_poc[d=1,l=53]-model_poc[d=1,l=25])
ppl shakey 1,0,0.2, , , , , ,`($ppl$ylen)*(0.1)`,`($ppl$ylen)*(0.15)`
ppl shade
go fland;go land black
con/ov/levels=(-inf)(-100,100,50)(inf)/nolab (model_poc[d=1,l=53]-model_poc[d=1,l=25])
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "POC (mmol C m^-^2 y^-^1)"

set view v2
shade/set/pal=blue_darkred/levels=(-inf)(-100,100,25)(inf)/nolab (model_pic[d=1,l=53]-model_pic[d=1,l=25])
ppl shakey 1,0,0.2, , , , , ,`($ppl$ylen)*(0.1)`,`($ppl$ylen)*(0.15)`
ppl shade
go fland;go land black
con/ov/levels=(-inf)(-100,100,50)(inf)/nolab (model_pic[d=1,l=53]-model_pic[d=1,l=25])
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "PIC (mmol C m^-^2 y^-^1)"

set view v3
shade/set/pal=blue_darkred/levels=(-inf)(-300,300,100)(inf)/nolab (model_opl[d=1,l=53]-model_opl[d=1,l=25])
ppl shakey 1,0,0.2, , , , , ,`($ppl$ylen)*(0.1)`,`($ppl$ylen)*(0.15)`
ppl shade
go fland;go land black
con/ov/levels=(-inf)(-300,300,150)(inf)/nolab (model_opl[d=1,l=53]-model_opl[d=1,l=25])
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "Opal (mmol Si m^-^2 y^-^1)"

go finalizeFigure_psthicken