! 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,l=24,d=1]+o_detrexp_B[k=11,l=24,d=1])*mult1*ntc*1000*350 ! mmol C/y
let model_pic=o_caco3exp[k=11,l=24,d=1]*mult1*1000*350 ! mmol C/y
let model_opl=o_oplexp[k=11,l=24,d=1]*mult1*1000*350 ! mmol Si/y

! Keller model data
use "$WORK/Keller_transient/tavg.01775.01.01.nc"
let model_poc2=(o_detrexp[k=11,l=24,d=2])*mult1*ntc*1000*350 ! mmol C/y
let model_pic2=o_caco3exp[k=11,l=24,d=2]*mult1*1000*350 ! mmol C/y

go prepareFigure fluxes_comp
go "$HOME/opt_postprocessing_scripts/go/ilandscape3x2.jnl"
ppl axlsze,0.15,0.15
ppl axlint 4 2
ppl labset,0.15,0.15,0.15,0.15
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set view v1
shade/nokey/set/pal=purple_red/levels=(0,300,50)(inf)/nolab model_poc
!ppl shakey 1,1,0.2, , , , , , ,`($ppl$ylen)*(0.15)`
ppl shade
go fland;go land black
go "$HOME/Scripts/David_scripts/readhonjotraps.jnl"
go polymark polygon/pal=purple_red/fill/ov/nolab/levels=(0,300,50)(inf)/line hlon hlat hfco circle 1.0
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "POC (mmol C m^-^2 y^-^1)"

set view v4
let sammodel1 = samplexy(model_poc,hlon,hlat)
plot/set/nolab/vs/hlog/vlog/hlim=1:1e+3/vlim=1:1e+3/sym=2/col=red/thick=1 hfco,sammodel1
ppl xlab "Observations (mmol C m^-^2 y^-^1)"
ppl ylab "Model (mmol C m^-^2 y^-^1)"
ppl plot
def ax/x=1:10000:10 dux
def g/x=dux dug
plot/ov/nolab/hlog/vlog/line=7 x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 2*x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 0.5*x[g=dug]

let sdiff = (sammodel1-hfco)^2
let sdiffsum = sdiff[x=@sum]
let sdiffngd = sdiff[x=@ngd]
let sdiffave = (sdiffsum/sdiffngd)^0.5
label 0.1,3.7,-1,0,.2 rms @AS `(INT(sdiffave*100))/100`
label 0.1,3.4,-1,0,.2 n @AS `sdiffngd`

let p = log(hfco)
let q = log(sammodel1)

go regressx

let sint = INT(`intercep`*100)/100
let sslo = INT(`slope`*100)/100
let srsq = INT(`rsquare`*100)/100

let sammodel2 = samplexy(model_poc2,hlon,hlat)
plot/ov/nolab/vs/hlog/vlog/sym=1/col=blue/thick=3 hfco,sammodel2

!label 2,4.5,0,0,.1 All: log(m) = `sint`+`sslo`*log(s), R^2=`srsq`
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set view v2
shade/nokey/set/pal=purple_red/levels=(0,300,50)(inf)/nolab model_pic[d=1]
ppl shakey 1,0,0.15, , , , , ,`($ppl$ylen)*(0.15)`,`($ppl$ylen)*(0.25)`
ppl shade
go fland;go land black
go polymark polygon/pal=purple_red/fill/ov/nolab/levels=(0,300,50)(inf)/line hlon[d=2] hlat[d=2] hfci[d=2] circle 1.0
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "PIC (mmol C m^-^2 y^-^1)"

set view v5
let sammodel2 = samplexy(model_pic,hlon,hlat)
plot/set/nolab/vs/hlog/vlog/hlim=1:1e+3/vlim=1:1e+3/sym=2/col=red/thick=1 hfci,sammodel2
ppl xlab "Observations (mmol C m^-^2 y^-^1)"
ppl ylab "Model (mmol C m^-^2 y^-^1)"
ppl plot
def ax/x=1:10000:10 dux
def g/x=dux dug
plot/ov/nolab/hlog/vlog/line=7 x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 2*x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 0.5*x[g=dug]

let sdiff = (sammodel2-hfci)^2
let sdiffsum = sdiff[x=@sum]
let sdiffngd = sdiff[x=@ngd]
let sdiffave = (sdiffsum/sdiffngd)^0.5
label 0.1,3.7,-1,0,.2 rms @AS `(INT(sdiffave*100))/100`
label 0.1,3.4,-1,0,.2 n @AS `sdiffngd`

let p = log(hfci)
let q = log(sammodel2)

go regressx

let sint = INT(`intercep`*100)/100
let sslo = INT(`slope`*100)/100
let srsq = INT(`rsquare`*100)/100
let sammodel3 = samplexy(model_pic2,hlon,hlat)
plot/ov/nolab/vs/hlog/vlog/sym=1/col=blue/thick=3 hfci,sammodel3

!label 2,4.5,0,0,.1 All: log(m) = `sint`+`sslo`*log(s), R^2=`srsq`
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
set view v3
shade/set/pal=purple_red/levels=(0,300,50)(inf)/nolab model_opl[d=1]
ppl shakey 0,0,0.15, , , ,`($ppl$xlen)*1.2`,`($ppl$xlen)*1.3`,`($ppl$ylen)*.05` ,`($ppl$ylen)*.95`
ppl shade
go fland;go land black
go "$HOME/Scripts/David_scripts/readhonjotraps.jnl"
go polymark polygon/pal=purple_red/fill/ov/nolab/levels=(0,300,50)(inf)/line hlon[d=2] hlat[d=2] hfsi[d=2] circle 1.0
label/nouser `($ppl$xlen)*.5` `($ppl$ylen)*(1.1)`  0,0,0.2 "Opal (mmol Si m^-^2 y^-^1)"


set view v6
let sammodel3 = samplexy(model_poc,hlon,hlat)
plot/set/nolab/vs/hlog/vlog/hlim=1:1e+3/vlim=1:1e+3/sym=2/col=red/thick=1 hfsi,sammodel3
ppl xlab "Observations (mmol Si m^-^2 y^-^1)"
ppl ylab "Model (mmol Si m^-^2 y^-^1)"
ppl plot
def ax/x=1:10000:10 dux
def g/x=dux dug
plot/ov/nolab/hlog/vlog/line=7 x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 2*x[g=dug]
plot/ov/nolab/hlog/vlog/line=7/dash 0.5*x[g=dug]

let sdiff = (sammodel3-hfsi)^2
let sdiffsum = sdiff[x=@sum]
let sdiffngd = sdiff[x=@ngd]
let sdiffave = (sdiffsum/sdiffngd)^0.5
label 0.1,3.7,-1,0,.2 rms @AS `(INT(sdiffave*100))/100`
label 0.1,3.4,-1,0,.2 n @AS `sdiffngd`

let p = log(hfsi)
let q = log(sammodel3)

go regressx

let sint = INT(`intercep`*100)/100
let sslo = INT(`slope`*100)/100
let srsq = INT(`rsquare`*100)/100

!label 2,4.5,0,0,.1 All: log(m) = `sint`+`sslo`*log(s), R^2=`srsq`
go finalizeFigure_psthicken