! This PyFerret script plots FIG_S3 of the Paper "Artificial Upwelling - A New Narrative" 
! and was written by M. Juerchott (mjuerchott@geomar.de) on the 12th September 2022.
! In order to run this script without error you have to change the path to the model output files to your personal folder structure.
 
 
! ***  Only run it in -nodisplay -nojnl mode ***
 
 
! Save current state
 
	define REGION/DEFAULT "save"
	set data/save
	set grid/save
 
! Clean everything
 
	cancel data/all
	cancel variable/all
	cancel region/all
 
 
! Script function

! Select model output files

	use /Users/mjuerchott/runs/2021_02_Pipes_noLand_COU85/tavg.02025.01.01.nc		![d=1]		ArtUp_8.5
	use /Users/mjuerchott/runs/2021_11_Pipes_noLand_COU60/tavg.02025.01.01.nc		![d=2]		ArtUp_6.0
	use /Users/mjuerchott/runs/2021_11_Pipes_noLand_COU45/tavg.02025.01.01.nc		![d=3]		ArtUp_4.5
	use /Users/mjuerchott/runs/2021_11_Pipes_noLand_COU26/tavg.02025.01.01.nc		![d=4]		ArtUp_2.6
	use /Users/mjuerchott/runs/2021_02_Pipes_noLand_CTRL/tavg.02025.01.01.nc		![d=5]		ArtUp_0.0

	use /Users/mjuerchott/runs/2021_02_noPipes_noLand_COU85_cont/tavg.02025.01.01.nc	![d=6]		REF_8.5
	use /Users/mjuerchott/runs/2021_11_noPipes_noLand_COU60_cont/tavg.02025.01.01.nc	![d=7]		REF_6.0
	use /Users/mjuerchott/runs/2021_11_noPipes_noLand_COU45_cont/tavg.02025.01.01.nc	![d=8]		REF_4.5
	use /Users/mjuerchott/runs/2021_11_noPipes_noLand_COU26_cont/tavg.02025.01.01.nc	![d=9]		REF_2.6
	use /Users/mjuerchott/runs/2021_02_noPipes_noLand_CTRL_cont/tavg.02025.01.01.nc		![d=10]		REF_0.0



! calculation ArtUp - REF files and convert unit into Pg carbon

	let xDIC_85 =	(O_DIC[d=1]  - O_DIC[d=6])	  *(12 *10^(-15))
	let xDIC_60 =	(O_DIC[d=2]  - O_DIC[d=7])        *(12 *10^(-15))
	let xDIC_45 =	(O_DIC[d=3]  - O_DIC[d=8])        *(12 *10^(-15))
	let xDIC_26 =	(O_DIC[d=4]  - O_DIC[d=9])        *(12 *10^(-15))
	let xDIC_00 =	(O_DIC[d=5]  - O_DIC[d=10])       *(12 *10^(-15))

	let/title="RPC 8.5"/units="<#916>C [Pg]" DIC_85  = xDIC_85[I=@DIN,J=@DIN,K=@DIN]
	let/title="RCP 6.0"/units="<#916>C [Pg]" DIC_60  = xDIC_60[I=@DIN,J=@DIN,K=@DIN]
	let/title="RPC 4.5"/units="<#916>C [Pg]" DIC_45  = xDIC_45[I=@DIN,J=@DIN,K=@DIN]
	let/title="RPC 2.6"/units="<#916>C [Pg]" DIC_26  = xDIC_26[I=@DIN,J=@DIN,K=@DIN]
	let/title="RPC 0.0"/units="<#916>C [Pg]" DIC_00  = xDIC_00[I=@DIN,J=@DIN,K=@DIN]


! pipe area over time

	let pipe_mask_85 = if O_ZPIPE[d=1] GT 50 then 1 else 0
	let pipe_mask_60 = if O_ZPIPE[d=2] GT 50 then 1 else 0
	let pipe_mask_45 = if O_ZPIPE[d=3] GT 50 then 1 else 0
	let pipe_mask_26 = if O_ZPIPE[d=4] GT 50 then 1 else 0
	let pipe_mask_00 = if O_ZPIPE[d=5] GT 50 then 1 else 0

	let pipe_area_m2_85 = (G_AREAT[d=1] * pipe_mask_85)
	let pipe_area_m2_60 = (G_AREAT[d=2] * pipe_mask_60)
	let pipe_area_m2_45 = (G_AREAT[d=3] * pipe_mask_45)
	let pipe_area_m2_26 = (G_AREAT[d=4] * pipe_mask_26)
	let pipe_area_m2_00 = (G_AREAT[d=5] * pipe_mask_00)

	let/unit="[km<&sup2>]" pipe_area_km2_85 = pipe_area_m2_85[i=@SUM,j=@SUM] / 1000000
	let/unit="[km<&sup2>]" pipe_area_km2_60 = pipe_area_m2_60[i=@SUM,j=@SUM] / 1000000
	let/unit="[km<&sup2>]" pipe_area_km2_45 = pipe_area_m2_45[i=@SUM,j=@SUM] / 1000000
	let/unit="[km<&sup2>]" pipe_area_km2_26 = pipe_area_m2_26[i=@SUM,j=@SUM] / 1000000
	let/unit="[km<&sup2>]" pipe_area_km2_00 = pipe_area_m2_00[i=@SUM,j=@SUM] / 1000000


! intensity calculation C-uptake / pipe area

	let intensity_85 = (DIC_85 / pipe_area_km2_85) * 10^7
	let intensity_60 = (DIC_60 / pipe_area_km2_60) * 10^7
	let intensity_45 = (DIC_45 / pipe_area_km2_45) * 10^7
	let intensity_26 = (DIC_26 / pipe_area_km2_26) * 10^7
	let intensity_00 = (DIC_00 / pipe_area_km2_00) * 10^7



! plot routine

	cancel windows/all
	set windows/new
	set text/font='Times New Roman'
	set text/isiz=4


	go landscape2x2.jnl
	define view/xlimits=0.2,0.7/ylimits=0.05,0.55 ll4
	define view/xlimits=0.2,0.7/ylimits=0.5,1 ul4

	set view ul4
	plot/nolabel/color=black/hlimit=2020:2101:10/vlimit=0:3:0.5/thick=2 intensity_85
	plot/nolabel/overlay/color=blue/thick=2 intensity_60
	plot/nolabel/overlay/color=red/thick=2 intensity_45
	plot/nolabel/overlay/color=green/thick=2 intensity_26
	plot/nolabel/overlay/color=(100,68,0)/thick=2 intensity_00

	label/nouser -1,1.7,0,90,0.13 "(<#916>C [Pg] / pipe covered area [km<&sup2>]) *10<sup>7</sup>"
	label/nouser 2.5,-0.8,0,0,0.13 "[years]"



! plot legend

	set view ll4

	set text/font='Times New Roman'/color=black
	label/nouser 0.2,8.5,-1,0,0.11 - RCP 8.5

	set text/font='Times New Roman'/color=blue
	label/nouser 0.2,8.2,-1,0,0.11 - RCP 6.0

	set text/font='Times New Roman'/color=red
	label/nouser 0.2,7.9,-1,0,0.11 - RCP 4.5

	set text/font='Times New Roman'/color=green
	label/nouser 0.2,7.6,-1,0,0.11 - RCP 2.6

	set text/font='Times New Roman'/color=(100,68,0)
	label/nouser 0.2,7.3,-1,0,0.11 - no emission

	set text/font='Times New Roman'/color=black




frame/file = "FIG_S3.png"



! Cancel everything
 
	cancel data/all
	cancel variable/all
	cancel region/all
 
! Restore current state
 
	set data/restore
	set grid/restore