!
! Ferret/PyFerret Script to plot Figure 11
!
! Biastoch, Sein, Durgadoo, Wang, Danilov
! Simulating the Agulhas system in global ocean models – nesting vs. multi-resolution 
! unstructured meshes, Ocean Modelling, doi: 10.1016/J.OCEMOD-D-17-00174
!
! Note that Agulhas leakage needs to be calculated through a Lagrangian estimate
go default
cancel region

define symbol print=no
define symbol filename=compare_FESOM-NEMO_shift
IF ($program_name"0|PyFerret>1|*>0") THEN 
   cancel text
   SET TEXT/FONT=arial hlab 
   SET TEXT/FONT=arial vlab
   SET TEXT/FONT=arial movable
ELSE
   if ($print) then set mode metafile ($filename).plt endif
ENDIf

use Ariane_FESOM_clim.nc
use Ariane_FESOM_shift.nc
use AL_INALT01_309_340_341.nc
use INALT01_SHW40p.nc
use agu08_sbarotrob_1981_2020_glob_mm_INALT01.nc
use agu09_sbarotrob_1981_2020_glob_mm_INALT01.nc

define view/x=0.00:0.55/y=0.425:1 ul
define view/x=0.45:1.00/y=0.425:1 ur
define view/x=0.00:0.55/y=0:0.575 ll
define view/x=0.45:1.00/y=0:0.575 lr

define axis/t="01-jan-0001":"31-dec-0062":1/units="years"/edges/cal=noleap tFclim
define axis/t="01-jan-0031":"31-dec-0061":1/units="years"/edges/cal=noleap tFshft
define axis/t="01-jan-0031":"31-dec-0060":1/units="years"/edges/cal=noleap tIclim
define axis/t="01-jan-0031":"31-dec-0060":1/units="years"/edges/cal=noleap tIshft

define axis/t="01-jan-0031":"31-dec-0060":1/units="months"/edges/cal=noleap tFmon
let ps=psi[gt=tFmon@asn]

define axis/t="01-jan-0031":"31-dec-0060":1/units="years"/edges/cal=noleap tann

let AL_FESOM_clim=ignore0(WEST[d=1,gt=tFclim@asn]+NW[d=1,gt=tFclim@asn])
let AL_FESOM_shft=ignore0(WEST[d=2,gt=tFshft@asn]+NW[d=2,gt=tFshft@asn])
let AL_INALT_clim=KJD309[d=3,gt=tIclim@asn]
let AL_INALT_shft=KJD340[d=3,gt=tIshft@asn]

let ACC_FESOM_clim=ps[d=5,i=615:634@max,j=1:250@max,gt=tann@ave]
let ACC_FESOM_shft=ps[d=6,i=615:634@max,j=1:250@max,gt=tann@ave]
let ACC_INALT_clim=acc309[d=4,gt=tIclim@asn]
let ACC_INALT_shft=acc340[d=4,gt=tIshft@asn]


set region/t="01-jan-0030":"31-dec-0060"

set win/asp=1

set view ul
plot/vlim=8:28/hlim="01-jan-0030":"31-dec-0060"/thick=3/nol/dashed/set AL_FESOM_clim
ppl ylab Sv
ppl axlabp,1,-1
ppl plot
plot/ov/nol/thick=3/col=black AL_FESOM_shft
LET datelabelpos=TT[T="01-JAN-0032"@ITP]; label ($PPL$XFIRST1),26,-1,0,.14,(a) FESOM AL

set view ur
plot/vlim=8:28/hlim="01-jan-0030":"31-dec-0060"/thick=3/nol/col=red/dashed/set AL_INALT_clim
ppl ylab Sv
ppl axlabp,1,1
ppl plot
plot/ov/nol/thick=3/col=red   AL_INALT_shft
LET datelabelpos=TT[T="01-JAN-0032"@ITP]; label ($PPL$XFIRST1),26,-1,0,.14,(b) NEMO AL

set region/t="01-jan-0035":"31-dec-0040"
list AL_FESOM_clim[l=@ave],AL_FESOM_shft[l=@ave],(AL_FESOM_shft[l=@ave]-AL_FESOM_clim[l=@ave])/AL_FESOM_clim[l=@ave]*100
list AL_INALT_clim[l=@ave],AL_INALT_shft[l=@ave],(AL_INALT_shft[l=@ave]-AL_INALT_clim[l=@ave])/AL_INALT_clim[l=@ave]*100
list ACC_FESOM_clim[l=@ave],ACC_FESOM_shft[l=@ave],(ACC_FESOM_shft[l=@ave]-ACC_FESOM_clim[l=@ave])/ACC_FESOM_clim[l=@ave]*100
list ACC_INALT_clim[l=@ave],ACC_INALT_shft[l=@ave],(ACC_INALT_shft[l=@ave]-ACC_INALT_clim[l=@ave])/ACC_INALT_clim[l=@ave]*100

set region/t="01-jan-0045":"31-dec-0050"
list AL_FESOM_clim[l=@ave],AL_FESOM_shft[l=@ave],(AL_FESOM_shft[l=@ave]-AL_FESOM_clim[l=@ave])/AL_FESOM_clim[l=@ave]*100
list AL_INALT_clim[l=@ave],AL_INALT_shft[l=@ave],(AL_INALT_shft[l=@ave]-AL_INALT_clim[l=@ave])/AL_INALT_clim[l=@ave]*100
list ACC_FESOM_clim[l=@ave],ACC_FESOM_shft[l=@ave],(ACC_FESOM_shft[l=@ave]-ACC_FESOM_clim[l=@ave])/ACC_FESOM_clim[l=@ave]*100
list ACC_INALT_clim[l=@ave],ACC_INALT_shft[l=@ave],(ACC_INALT_shft[l=@ave]-ACC_INALT_clim[l=@ave])/ACC_INALT_clim[l=@ave]*100
cancel region

set region/t="01-jan-0030":"31-dec-0060"

set view ll
plot/vlim=170:240/hlim="01-jan-0030":"31-dec-0060"/thick=3/nol/dashed/set ACC_FESOM_clim
ppl ylab Sv
ppl axlabp,-1,-1
ppl plot
plot/ov/nol/thick=3/col=black ACC_FESOM_shft
LET datelabelpos=TT[T="01-JAN-0032"@ITP]; label ($PPL$XFIRST1),232,-1,0,.14,(c) FESOM ACC

set view lr
plot/vlim=170:240/hlim="01-jan-0030":"31-dec-0060"/thick=3/dashed/col=red/nol/set ACC_INALT_clim
ppl ylab Sv
ppl axlabp,1,1
ppl plot
plot/ov/nol/thick=3/col=red ACC_INALT_shft
LET datelabelpos=TT[T="01-JAN-0032"@ITP]; label ($PPL$XFIRST1),232,-1,0,.14,(d) NEMO ACC

!let acc=acc309[d=4]
!go my_right_axis_plot acc*0 "170,300,10" "ACC (Sv)" 0.14 " " "thick=2" " 0 0 100" 
!go my_right_axis_plot acc "170,300,10" "ACC (Sv)" 0.14 " " "thick=2" " 0 0 100" 
!let acc=acc340[d=4]
!go my_right_axis_plot acc "170,300,10" "ACC (Sv)" 0.14 " " "thick=3" " 0 0 100" 
!cancel region

!let acc=ps[d=5,i=615:634@max,j=1:250@max,gt=tann@ave]
!go my_right_axis_plot acc "170,300,10" "ACC (Sv)" 0.14 " " "thick=2" " 0 100 0" 

!let acc=ps[d=6,i=615:634@max,j=1:250@max,gt=tann@ave]
!go my_right_axis_plot acc "170,300,10" "ACC (Sv)" 0.14 " " "thick=3" " 0 100 0" 

IF ($print) THEN
IF ($program_name"0|PyFerret>1|*>0") THEN 
   FRAME/FILE=($filename).pdf
ELSE
   CANCEL MODE META
   sp Fprint -o ($filename).ps -l cps -p portrait ($filename).plt
ENDIF
ENDIF