!script by bencat.xuan@gmail.com
!!!!! load in data from 3 calibrations


use ../model_and_output_data/noIron_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc
use ../model_and_output_data/ironMask_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc
use ../model_and_output_data/ironDyn1_spinup_fct_drift/1800_2100/tavg.01801.01.01.nc





set mem/size=2048

def sym op1=x=@din,y=@din,k=@din
def sym op2=x=@din,y=@din,k=1:2@din
def sym op3=x=@din,y=@din,k=2@din
def sym op4=k=2@din
def sym s2yC=3600*24*365*6.625*12
def sym s2y=3600*24*365
def sym molN2PgN=14.0067/1e15
def sym molN2TgN=14.0067/1e12
def sym gN2gC=12/14.0067*6.625
def sym molN2gC=12*6.625

def sym iopind=x=75e:110E@din,y=0:20n@din
def sym ioppa=x=115w:60w@din,y=8s:20n@din
def sym iopatl=x=0:20E@din,y=20s:0@din

def sym opind=x=75e:110E,y=0:20n
def sym oppa=x=115w:60w,y=8s:20n
def sym opatl=x=0:20E,y=20s:0

let nofe_exp_1800=o_DETRexp[d=1,l=1,($op4)]
let nofe_exp_2100=o_DETRexp[d=1,l=300,($op4)]
let nofe_expdiff=nofe_exp_2100-nofe_exp_1800
let/bad=-999 nofe_exp_pa=if G_MSKHR eq 2 then o_DETRexp[d=1] else (-999)
let nofe_exp_odz=o_DETRexp[d=1,($op4),($iopind)]+ nofe_exp_pa[d=1,($op4),($ioppa)] +o_DETRexp[d=1,($op4),($iopatl)]


let femask_exp_1800=o_DETRexp[d=2,l=1,($op4)]
let femask_exp_2100=o_DETRexp[d=2,l=300,($op4)]
let femask_expdiff=femask_exp_2100-femask_exp_1800
let/bad=-999 femask_exp_pa=if G_MSKHR eq 2 then o_DETRexp[d=2] else (-999)
let femask_exp_odz=o_DETRexp[d=2,($op4),($iopind)]+femask_exp_pa[d=2,($op4),($ioppa)]+o_DETRexp[d=2,($op4),($iopatl)]



let fedyn1_exp_1800=o_DETRexp[d=3,l=1,($op4)]
let fedyn1_exp_2100=o_DETRexp[d=3,l=300,($op4)]
let fedyn1_expdiff=fedyn1_exp_2100-fedyn1_exp_1800
let/bad=-999 fedyn1_exp_pa=if G_MSKHR eq 2 then o_DETRexp[d=3] else (-999)
let fedyn1_exp_odz=o_DETRexp[d=3,($op4),($iopind)]+fedyn1_exp_pa[d=3,($op4),($ioppa)]+o_DETRexp[d=3,($op4),($iopatl)]


def sym iscale="(-inf)(-24,24,4)(inf)"
def sym iscalennpdiff="(-inf)(-24,24,8)(inf)"




go ilandscape3x1_ext.jnl


set win 1
go prepareFigure Fig4
go transparency
set vie v1
sha/nolab/set/lev=($iscale)/pal=blue_orange/nokey nofe_expdiff[X=30E:30E]*($s2y)*($molN2gC)
ppl ylab "diff Export 2100-1800"
ppl axlabp,1,-1
ppl xlab "NoFe"
ppl shade
go land
contour/overlay/size=0.1/pen=1/nolab/lev=($iscalennpdiff) nofe_expdiff[X=30E:30E]*($s2y)*($molN2gC)
go box_dashed 75, 110, 0, 22,14
go box_dashed 245, 300, -8, 20,14
go box_dashed 360, 380, -20, 0,14
label -8,111,-1,0.2 (a)

set vie v2
sha/nolab/set/lev=($iscale)/pal=blue_orange/nokey femask_expdiff[X=30E:30E]*($s2y)*($molN2gC)
ppl xlab "FeMask"
ppl axlabp,1,0
ppl shade
go land
contour/overlay/size=0.1/pen=1/nolab/lev=($iscalennpdiff) femask_expdiff[X=30E:30E]*($s2y)*($molN2gC)
go box_dashed 75, 110, 0, 22,14
go box_dashed 245, 300, -8, 20,14
go box_dashed 360, 380, -20, 0,14
label -8,111,-1,0.2 (b)

set vie v3
sha/nolab/set/lev=($iscale)/pal=blue_orange fedyn1_expdiff[X=30E:30E]*($s2y)*($molN2gC)
ppl xlab "FeDyn"
ppl axlabp,1,0
ppl shade
go land
contour/overlay/size=0.1/pen=1/nolab/lev=($iscalennpdiff) fedyn1_expdiff[X=30E:30E]*($s2y)*($molN2gC)
go box_dashed 75, 110, 0, 22,14
go box_dashed 245, 300, -8, 20,14
go box_dashed 360, 380, -20, 0,14
label -8,111,-1,0.2 (c)
go finalizeFigure_psthicken

!go prepareFigure nfct_exp_total
go transparency
set win 2
set vie v1
plot/nolab/hlim=1800:2105/set  (fedyn1_exp_odz[l=@sbx:5]-fedyn1_exp_odz[l=1])*($s2y)*($molN2gC)/1e12 ,(femask_exp_odz[l=@sbx:5]-femask_exp_odz[l=1])*($s2y)*($molN2gC)/1e12  ,(nofe_exp_odz[l=@sbx:5]-nofe_exp_odz[l=1])*($s2y)*($molN2gC)/1e12
ppl ylab "diff. export in all ODZ regions (Tg C yr^-^1)"
ppl pen 1, 10
ppl pen 2, 9
ppl pen 3, 8
ppl axlabp,1,-1
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"
!go finalizeFigure_psthicken


stat (fedyn1_exp_odz[l=300]-fedyn1_exp_odz[l=1])*($s2y)*($molN2gC)/1e12
stat (femask_exp_odz[l=300]-femask_exp_odz[l=1])*($s2y)*($molN2gC)/1e12
stat (nofe_exp_odz[l=300]-nofe_exp_odz[l=1])*($s2y)*($molN2gC)/1e12