! plots atmospheric CO2, freshwater forcing, net primary production, Atlantic meridional overturning circulation and limitation factors for PO4 and Fe
! revised script: added HOSING scenario for pure freshwater forcing (additionally to CONTROL, CO2X2, CO2X4, CO2X2_HOSING and CO2X4_HOSING)

cancel data/all
cancel window/all
set memory/size=200

set window/quality=high/thicken=1.2/aspect=1.25:axis 1


!----------------SIMULATIONSDATEN LADEN----------------------------
use "/p/projects/climber3/hofmann/RL_NEW/CO2X4/snap_ocean_1000.nc" 1 
use "/p/projects/climber3/hofmann/RL_NEW/CO2X4/history.nc" 2
use "/p/projects/climber3/hofmann/RL_NEW/CO2X2/snap_ocean.nc" 3
use "/p/projects/climber3/hofmann/RL_NEW/CO2X2/history.nc" 4
use "/p/projects/climber3/hofmann/RL_NEW/CO2X4_HOSING/snap_ocean_hosing_1000.nc" 5 
use "/p/projects/climber3/hofmann/RL_NEW/CO2X4_HOSING/history.nc" 6
use "/p/projects/climber3/hofmann/RL_NEW/CONTROL/snap_ocean_1000.nc" 7
use "/p/projects/climber3/hofmann/RL_NEW/CONTROL/history.nc" 8
use "/p/projects/climber3/hofmann/RL_NEW/CO2X2_HOSING/snap_ocean_hosing.nc" 9 
use "/p/projects/climber3/hofmann/RL_NEW/CO2X2_HOSING/history.nc" 10

use "/p/projects/climber3/hofmann/RL_1000/snap_ocean_hosing_ink_1000.nc" 11
use "/p/projects/climber3/hofmann/RL_1000/history_1000.nc" 12

!----------------GRAFIKEINSTELLUNGEN-------------------------------
! for 6 plots on a page (2 columns)
define viewport/xlimits=0,0.55/ylim=0.05,0.35 leftlow
define viewport/xlimits=0.45,0.99/ylim=0.05,0.35 rightlow
define viewport/xlimits=0,0.55/ylim=0.36,0.66 leftmid
define viewport/xlimits=0.45,0.99/ylim=0.36,0.66 rightmid
define viewport/xlimits=0,0.55/ylim=0.67,0.97 leftup
define viewport/xlimits=0.45,0.99/ylim=0.67,0.97 rightup

! for 8 plots on a page (2 columns)
define viewport/xlimits=0,0.55/ylim=0.05,0.27 pos14
define viewport/xlimits=0.45,0.99/ylim=0.05,0.27 pos24
define viewport/xlimits=0,0.55/ylim=0.28,0.50 pos13
define viewport/xlimits=0.45,0.99/ylim=0.28,0.50 pos23
define viewport/xlimits=0,0.55/ylim=0.51,0.73 pos12
define viewport/xlimits=0.45,0.99/ylim=0.51,0.73 pos22
define viewport/xlimits=0,0.55/ylim=0.74,0.96 pos11
define viewport/xlimits=0.45,0.99/ylim=0.74,0.96 pos21
!----------------FORMELN---------------------------------
define axis/T= -50:600:1  ngt! -50:780:1  ngt
define grid/T=ngt  nngt
let ll=L[gt=nngt@asn]
!
let argT1 =  32.0 - temp[d=1]
let flimit1 = 1.e-3 * (max(0.08,argT1/(argT1 + 20.0)))
let fvalue1 = tracer_13[d=1]/(tracer_13[d=1] + flimit1)
!
let plimit1 = tracer_06[d=1]/((1./122.)+tracer_06[d=1])
!
let argT3 =  32.0 - temp[d=3]
let flimit3 = 1.e-3 * (max(0.08,argT3/(argT3 + 20.0)))
let fvalue3 = tracer_13[d=3]/(tracer_13[d=3] + flimit3)
!
let plimit3 = tracer_06[d=3]/((1./122.)+tracer_06[d=3])
!
let argT5 =  32.0 - temp[d=5]
let flimit5 = 1.e-3 * (max(0.08,argT5/(argT5 + 20.0)))
let fvalue5 = tracer_13[d=5]/(tracer_13[d=5] + flimit5)
!
let plimit5 = tracer_06[d=5]/((1./122.)+tracer_06[d=5])
!
let argT7 =  32.0 - temp[d=7]
let flimit7 = 1.e-3 * (max(0.08,argT7/(argT7 + 20.0)))
let fvalue7 = tracer_13[d=7]/(tracer_13[d=7] + flimit7)
!
let plimit7 = tracer_06[d=7]/((1./122.)+tracer_06[d=7])
!
let argT9 =  32.0 - temp[d=9]
let flimit9 = 1.e-3 * (max(0.08,argT9/(argT9 + 20.0)))
let fvalue9 = tracer_13[d=9]/(tracer_13[d=9] + flimit9)
!
let plimit9 = tracer_06[d=9]/((1./122.)+tracer_06[d=9])
!
let argT11 =  32.0 - temp[d=11]
let flimit11 = 1.e-3 * (max(0.08,argT11/(argT11 + 20.0)))
let fvalue11 = tracer_13[d=11]/(tracer_13[d=11] + flimit11)
!
let plimit11 = tracer_06[d=11]/((1./122.)+tracer_06[d=11])
!
!----------------ZEITACHSE----------------------------------------
!
define axis/units=years/t=-50:1000:1 nt1
define grid/T=nt1 ntime1
define axis/units=years/t=0:1000:1 nt2
define grid/T=nt2 ntime2
!
let tt=T[gt=ntime1@asn]
let t0=tt[T="01-jan-1901"]
let t1=tt[T="01-jan-2040"]
let t2=tt[T="01-jan-2180"]
let t3=tt[T="01-jan-1970"]
!
!----------------GRAFISCHE AUSGABE---------------------------------

SET TEXT/FONT=symbol
set window/xpix=1680/ypix=1920 1
!
set view pos11
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
ppl axset 0,1,1,0
plot/nolabels/thickness=2/color=red/vlimits=200:1250/hlimits=-50:1000:100 pco2diag[gt=ntime1@asn,d=2]
plot/over/nolabels/thickness=2/color=(100,55,0) pco2diag[gt=ntime1@asn,d=4]
plot/over/nolabels/thickness=2/color=black pco2diag[gt=ntime1@asn,d=8]
let tt=T[gt=ntime1@asn]
let t0=tt[T="01-jan-1901"]
let t1=tt[T="01-jan-2040"]
let t2=tt[T="01-jan-2180"]
let t3=tt[T="01-jan-1970"]
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser 4.7,2.2,0,0,.2 "a)"
label/nouser -0.9,0.7,0,90,.2 "pCO<sub>2</sub> (ppmv)"
!
set view pos21
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
ppl axset 0,1,1,0
let FWT4 = (0.05/4.) * (pco2diag[gt=ntime1@asn,d=6]/280.)
let val4 = if (ll GT 50 AND ll LT 330) then FWT4 else 0
let FWT2 = (0.05/4.) * (pco2diag[gt=ntime1@asn,d=10]/280.)
let val2 = if (ll GT 50 AND ll LT 190) then FWT2 else 0
plot/nolabels/thickness=2/color=blue/hlimits=-50:1000:100 val4 ! [l=1:500]
plot/over/nolabels/thickness=2/color=LightBlue val2 ! [l=1:500]
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser 4.7,2.2,0,0,.2 "b)"
label/nouser -0.9,1.2,0,90,.15 "NA freshw. add. F<sub>FWT</sub> (Sv)"
!
set view pos12
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
ppl axset 0,1,1,0
plot/nolabels/thickness=2/color=red/hlimits=-50:1000:100/vlimits=0:24 atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=2]
plot/over/nolabels/thickness=2/color=(100,55,0) atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=4]
plot/over/nolabels/thickness=2/color=blue atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=6]
plot/over/nolabels/thickness=2/color=black atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=8]
plot/over/nolabels/thickness=2/color=lightblue atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=10]
plot/over/dash=(0.1,0.1,0.1,0.1)/nolabels/thickness=2/color=green atl_over[gt=ntime1@asn,y=30n:60n@max,z=500:1500@max,d=12]
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser -0.9,1.3,0,90,.175 "Atlantic Overturning (Sv)"
label/nouser 4.7,2.2,0,0,.2 "c)"
!
set view pos22
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
ppl axset 0,1,1,0
plot/nolabels/thickness=2/color=red/hlimits=-50:1000:100 glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=2]*1.e-15
plot/over/nolabels/thickness=2/color=(100,55,0) glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=4]*1.e-15
plot/over/nolabels/thickness=2/color=blue glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=6]*1.e-15
plot/over/nolabels/thickness=2/color=black glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=8]*1.e-15
plot/over/nolabels/thickness=2/color=lightblue glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=10]*1.e-15
plot/over/dash=(0.4,0.4,0.4,0.4)/nolabels/thickness=2/color=green glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=12]*1.e-15
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser -0.9,1.3,0,90,.2 "NPP (GtC/a)"
label/nouser 4.7,2.2,0,0,.2 "d)"
!
set view pos13
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
plot/nolabels/thickness=2/color=red/hlimits=-50:1000:100/vlimits=4.8:8.5 glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=2]*1.e-15
plot/over/nolabels/thickness=2/color=(100,55,0) glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=4]*1.e-15
plot/over/nolabels/thickness=2/color=blue glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=6]*1.e-15
plot/over/nolabels/thickness=2/color=black glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=8]*1.e-15
plot/over/nolabels/thickness=2/color=lightblue glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=10]*1.e-15
plot/over/dash=(0.1,0.1,0.1,0.1)/nolabels/thickness=2/color=green glb_xpp[gt=ntime1@asn,x=130w:80w@din,y=20s:20n@din,d=12]*1.e-15
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser 4.7,2.2,0,0,.2 "@CRe)"
label/nouser -0.9,1.3,0,90,.2 "@CRNPP (GtC/a) at EEP"
!
set view pos23
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
plot/nolabels/thickness=2/color=red/hlimits=-50:1000:100/y=40s:80s/vlimits=13:18 glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=2]*1.e-15
plot/over/nolabels/thickness=2/color=(100,55,0)/y=40s:80s glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=4]*1.e-15
plot/over/nolabels/thickness=2/color=blue/y=40s:80s glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=6]*1.e-15
plot/over/nolabels/thickness=2/color=black/y=40s:80s glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=8]*1.e-15
plot/over/nolabels/thickness=2/color=lightblue/y=40s:80s glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=10]*1.e-15
plot/over/dash=(0.1,0.1,0.1,0.1)/nolabels/thickness=2/color=green/y=40s:80s glb_xpp[gt=ntime1@asn,i=@din,j=@din,d=12]*1.e-15
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser -0.9,1.3,0,90,.2 "@CRNPP (GtC/a) at SO"
label/nouser 4.7,2.2,0,0,.2 "@CRf)"
!
set view pos14
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
!plot plimit1[x=@ave,y=@ave,d=1,z=100]
plot/nolabels/thickness=2/color=red/vlimits=0.88:0.98:0.01/hlimits=-50:1000:100 plimit1[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=(100,55,0) plimit3[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=blue plimit5[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=black plimit7[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=lightblue plimit9[gt=ntime2@asn,i=@ave,j=@ave,z=100]
plot/over/dash=(0.1,0.1,0.1,0.1)/nolabels/thickness=2/color=green plimit11[gt=ntime2@asn,i=@ave,j=@ave,z=100]
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser 4.7,2.2,0,0,.2 "g)"
label/nouser -0.9,1.3,0,90,.2 "P availability L<sub>PO4</sub>"
label/nouser 2.2,-0.8,0,0,.2 "@CR Year of simulation"
!
set view pos24
GO margins  0.5 1.0 2 1.0
ppl labset, 0.12,0.12,0.12,0.4
ppl axlsze 0.14,0.14
!plot fvalue1[x=@ave,y=@ave,d=1,z=100]
plot/nolabels/thickness=2/color=red/vlimits=0.47:0.62:0.01/hlimits=-50:1000:100 fvalue1[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=(100,55,0) fvalue3[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=blue fvalue5[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=black fvalue7[gt=ntime1@asn,i=@ave,j=@ave,z=100]
plot/over/nolabels/thickness=2/color=lightblue fvalue9[gt=ntime2@asn,i=@ave,j=@ave,z=100]
plot/over/dash=(0.1,0.1,0.1,0.1)/nolabels/thickness=2/color=green fvalue11[gt=ntime2@asn,i=@ave,j=@ave,z=100]
plot/over/vs/line/color=black/nolabel {`T0`,`T0`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T1`,`T1`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T2`,`T2`},{($YAXIS_MIN),($YAXIS_MAX)}
plot/over/vs/line/color=black/nolabel {`T3`,`T3`},{($YAXIS_MIN),($YAXIS_MAX)}
label/nouser 4.7,2.2,0,0,.2 "@CRh)"
label/nouser -0.9,1.3,0,90,.2 "@CR Fe availability L<sub>Fe</sub>"
label/nouser 2.8,-0.8,0,0,.2 "@CR Year of simulation"

!
define view/x=0.30:0.40/y=0.82:0.92/axes legend
set view legend
plot/vs/sym=dot/nolab/noax/vlim=0:1/hlim=0:1/color=white {0,1}, {0,1}
plot/over/vs/nolab/line=1/thickness=3 {0.1,0.8}, {1.00,1.00}
annotate/norm/xpos=.1/ypos=1.10/halign=-1/valign=0/size=0.2 "CONTROL"
plot/over/vs/nolab/line/color=(100,55,0) {0.1,0.8}, {0.80,0.80}
annotate/norm/xpos=.1/ypos=0.90/halign=-1/valign=0/size=0.2 "2xCO2"
plot/over/vs/nolab/line=2/thickness=3 {0.1,0.8}, {0.60,0.60}
annotate/norm/xpos=.1/ypos=0.70/halign=-1/valign=0/size=0.2 "4xCO2"
plot/over/vs/nolab/line=5/thickness=3 {0.1,0.8}, {0.40,0.40}
annotate/norm/xpos=.1/ypos=0.50/halign=-1/valign=0/size=0.2 "2xCO2HOSING"
plot/over/vs/nolab/line=4/thickness=3 {0.1,0.8}, {0.20,0.20}
annotate/norm/xpos=.1/ypos=0.30/halign=-1/valign=0/size=0.2 "4xCO2HOSING"
plot/over/vs/nolab/line=3/thickness=3 {0.1,0.25}, {0.00,0.00}
plot/over/vs/nolab/line=3/thickness=3 {0.375,0.525}, {0.00,0.00}
plot/over/vs/nolab/line=3/thickness=3 {0.65,0.8}, {0.00,0.00}
annotate/norm/xpos=.1/ypos=0.10/halign=-1/valign=0/size=0.2 "HOSING"
!
frame/file="figures/COMFORT/paper/COMFORT_means_forcings-limits-NPP_rev.pdf"