!script by bencat.xuan@gmail.com




define axis/t=0:`1-1/12`:`1/12` monthax







!!!!! load in data from 3 calibrations

use ../model_and_output_data/ironDyn1_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc	!1 FeDyn1
use ../model_and_output_data/ironMask_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc	!2 FeMAsk
use ../model_and_output_data/noIron_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc	!3 noFe
set mem/size=2048

use "../model_and_output_data/ironDyn1_spinup_fct_drift/data/O_fe_dissolved.nc" !4
let griddata=o_dfe[d=1,l=1]
let noFe_abio=0.6979944957

let FeMask_kfe=5.020096e-5
let FeMask_kfe_D=0.0004060349
let FeMask_abio=0.6

let FeDyn1_kfe=0.0005845
let FeDyn1_kfe_D=0.0003765222
let FeDyn1_abio=0.6

let FeDyn2_kfe=0.000352
let FeDyn2_kfe_D=0.0005815124
let FeDyn2_abio=0.6


let dummy = compressl(O_no3[d=1,l=1,k=1]*0.)
def sym oper1 = l=1,k=1:2@din,gxyz=dummy
def sym oper2 = l=300,k=1:2@din,gxyz=dummy
def sym oper3 = l=1,k=2@din,gxyz=dummy
def sym oper4 = l=300,k=2@din,gxyz=dummy
def sym oper5= k=1:2@din,x=@din,y=@din
def sym oper6= k=1:2@din,x=@din,y=@din,l=@sbx:5
def sym oper7=x=@din,y=@din,k=@din
def sym oper8=x=@din,y=@din,k=@din,l=@sbx:5
def sym s2y=3600*24*365 !*6.625*12
def sym molN2PgN=14.0067/1e15
def sym molN2TgN=14.0067/1e12
def sym molN2GtC=12*6.625/1e15
def sym pgN2pgC=12/14.0067*6.625
def sym GtC2pgN=14.0067/12/6.625

let k1n=0.7
let redptn=1/16
let k1p=k1n*redptn
let bbio = 1.066  ! b
let cbio = 1.0  ! c [1/deg_C]

! calculation



go ilandscape2x1_ext.jnl

ppl axset 1,1,1,0 ! supress right hand axis
ppl axnsig 2,2

go prepareFigure Fig1
set win 1
set view v1
let hmin=1800
let hmax=2105
let vmin=-6
let vmax=0.6
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set (o_phytnpp[d=1,($oper6)]-o_phytnpp[d=1,($oper5),l=1])*($s2y)*($molN2GtC)+(o_diaznpp[d=1,($oper6)]-o_diaznpp[d=1,($oper5),l=1])*($s2y)*($molN2GtC),(o_phytnpp[d=2,($oper6)]-o_phytnpp[d=2,($oper5),l=1])*($s2y)*($molN2GtC)+(o_diaznpp[d=2,($oper6)]-o_diaznpp[d=2,($oper5),l=1])*($s2y)*($molN2GtC),(o_phytnpp[d=3,($oper6)]-o_phytnpp[d=3,($oper5),l=1])*($s2y)*($molN2GtC)+(o_diaznpp[d=3,($oper6)]-o_diaznpp[d=3,($oper5),l=1])*($s2y)*($molN2GtC)
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl plot
go tilegline_nu.jnl 0.1 0.1 0.24 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.30 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.36 8 0.2 "NoFe"
label `tposi*1.20`,`vmin+(vmax-vmin)/2`,0,90,.15 (Pg N yr^-^1)
label `tposi*(-0.18)`,`vmin+(vmax-vmin)/2`,0,90,0.15  @SGD@SR total NPP (Gt C yr^-^1)
label `tposi*(-0.18)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (a) 
ppl %yaxis `int(vmin*($GtC2pgN)*100)/100`,`int(vmax*($GtC2pgN)*100)/100`,`int((vmax-vmin)/10*($GtC2pgN)*100)/100`,`tposi`,0,2,1,+1
ppl axnsig 2,1
set view v2
let hmin=1800
let hmax=2105
let vmin=-0.1
let vmax=1.1
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set (o_diaznpp[d=1,($oper6)]-o_diaznpp[d=1,($oper5),l=1])*($s2y)*($molN2GtC),(o_diaznpp[d=2,($oper6)]-o_diaznpp[d=2,($oper5),l=1])*($s2y)*($molN2GtC),(o_diaznpp[d=3,($oper6)]-o_diaznpp[d=3,($oper5),l=1])*($s2y)*($molN2GtC)
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl plot
go tilegline_nu.jnl 0.1 0.1 0.24 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.30 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.36 8 0.2 "NoFe"
label `tposi*1.20`,`vmin+(vmax-vmin)/2`,0,90,.15  (Pg N yr^-^1)
label `tposi*(-0.18)`,`vmin+(vmax-vmin)/2`,0,90,0.15 @SGD@SR diaz PP (Gt C yr^-^1)
label `tposi*(-0.18)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (b) 
ppl %yaxis `int(vmin*($GtC2pgN)*100)/100`,`int(vmax*($GtC2pgN)*100)/100`,`int((vmax-vmin)/10*($GtC2pgN)*100)/100`,`tposi`,0,2,1,+1

go finalizeFigure_psthicken

list o_diaznpp[d=1,($oper5)]/(o_diaznpp[d=1,($oper5)]+o_phytnpp[d=1,($oper5)]),o_diaznpp[d=2,($oper5)]/(o_diaznpp[d=2,($oper5)]+o_phytnpp[d=2,($oper5)]),o_diaznpp[d=3,($oper5)]/(o_diaznpp[d=3,($oper5)]+o_phytnpp[d=3,($oper5)])


let ncp_fedyn=(o_phytnpp[d=1]+o_diaznpp[d=1]-o_phytrecy[d=1]-o_excret[d=1]-o_detrremi[d=1])*($molN2GtC)*($s2y)
let ncp_femask=(o_phytnpp[d=2]+o_diaznpp[d=2]-o_phytrecy[d=2]-o_excret[d=2]-o_detrremi[d=2])*($molN2GtC)*($s2y)
let ncp_nofe=(o_phytnpp[d=3]+o_diaznpp[d=3]-o_phytrecy[d=3]-o_excret[d=3]-o_detrremi[d=3])*($molN2GtC)*($s2y)


go isquare2x2_ext.jnl
ppl axset 1,1,1,1


go prepareFigure Fig2
set win 3
set view v1
let hmin=1800
let hmax=2105
let vmin=30
let vmax=58
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set o_no3fix[d=1,($oper6)]*($s2y)*($molN2TgN),o_no3fix[d=2,($oper6)]*($s2y)*($molN2TgN),o_no3fix[d=3,($oper6)]*($s2y)*($molN2TgN)
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl ylab "abs N_2 fix (Tg N yr^-^1)" 
ppl axnsig 2,2
ppl plot
label `tposi*(-0.22)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (a) 
go tilegline_nu.jnl 0.1 0.1 0.24 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.30 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.36 8 0.2 "NoFe"


set view v2
let hmin=1800
let hmax=2105
let vmin=-2
let vmax=22
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set (o_no3fix[d=1,($oper6)]-o_no3fix[d=1,($oper5),l=1])*($s2y)*($molN2TgN),(o_no3fix[d=2,($oper6)]-o_no3fix[d=2,($oper5),l=1])*($s2y)*($molN2TgN),(o_no3fix[d=3,($oper6)]-o_no3fix[d=3,($oper5),l=1])*($s2y)*($molN2TgN)
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl axnsig 2,2
ppl ylab "@SGD@SR N_2 fix (Tg N yr^-^1)" 
ppl plot
label `tposi*(-0.22)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (b) 
go tilegline_nu.jnl 0.1 0.1 0.24 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.30 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.36 8 0.2 "NoFe"



let FeDyn1=if o_o2[d=1] lt 0.005 then 1 else o_o2[d=1]*0
let FeMask=if o_o2[d=2] lt 0.005 then 1 else o_o2[d=2]*0
let NoFe=if o_o2[d=3] lt 0.005 then 1 else o_o2[d=3]*0
let ocean= o_o2[l=1,d=1]*0+1
let oceanvol=ocean[($oper7)]



set vie v3
let hmin=1800
let hmax=2105
let vmin=-15.0
let vmax=23.5
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set (o_DENITR[d=1,($oper8)]-o_DENITR[d=1,($oper7),l=1])*($s2y)*($molN2TgN),(o_DENITR[d=2,($oper8)]-o_DENITR[d=2,($oper7),l=1])*($s2y)*($molN2TgN),(o_DENITR[d=3,($oper8)]-o_DENITR[d=3,($oper7),l=1])*($s2y)*($molN2TgN)
ppl axnsig 2,3
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl ylab "@SGD@SR denitr (Tg N yr^-^1)" 
ppl plot
label `tposi*(-0.22)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (c) 
go tilegline_nu.jnl 0.1 0.1 0.64 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.70 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.76 8 0.2 "NoFe"


set vie v4
let hmin=1800
let hmax=2105
let vmin=-6e-2
let vmax=1.5e-2
let tposi=2654000
plot/nolab/hlimits=`hmin`:`hmax`/vlim=`vmin`:`vmax`/set (FeDyn1[($oper8)]-FeDyn1[($oper7),l=1])/oceanvol*100,(FeMask[($oper8)]-FeMask[($oper7),l=1])/oceanvol*100,(NoFe[($oper8)]-NoFe[($oper7),l=1])/oceanvol*100
ppl axnsig 2,1
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl ylab "@SGD@SR ODZ volume in ocean" 
ppl plot
label `tposi*(-0.22)`,`vmax+(vmax-vmin)*0.05`,0,0,0.15 (d) 
go tilegline_nu.jnl 0.1 0.1 0.24 10 0.2 "FeDyn"
go tilegline_nu.jnl 0.1 0.1 0.30 9 0.2 "FeMask"
go tilegline_nu.jnl 0.1 0.1 0.36 8 0.2 "NoFe"
go finalizeFigure_psthicken


list o_no3fix[d=1,($oper5),l=300]/o_no3fix[d=1,($oper5),l=1],o_no3fix[d=2,($oper5),l=300]/o_no3fix[d=2,($oper5),l=1],o_no3fix[d=3,($oper5),l=300]/o_no3fix[d=3,($oper5),l=1]


list (FeDyn1[($oper7),l=300]-FeDyn1[($oper7),l=1])/oceanvol*100,(FeMask[($oper7),l=300]-FeMask[($oper7),l=1])/oceanvol*100,(NoFe[($oper7),l=300]-NoFe[($oper7),l=1])/oceanvol*100



list (FeDyn1[($oper7),l=300]-FeDyn1[($oper7),l=1])/oceanvol*100,(FeMask[($oper7),l=300]-FeMask[($oper7),l=1])/oceanvol*100,(NoFe[($oper7),l=300]-NoFe[($oper7),l=1])/oceanvol*100