!
! Ferret/PyFerret Script to plot Figure 9
!
! 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 the ACT section transport has been calculated beforehand (see extract scripts)

go default
ppl conset ,2
ppl axlint,2,2

Set mem/size=500

cancel region; cancel view; cancel data/all;cancel var/all; can mode logo

define view/x=0:0.725/y=0.42:1/text=0.72 up
define view/x=0:0.725/y=0:0.58/text=0.72 lo
define view/x=0.57:0.95/y=0:1/text=0.72 ri


define symbol print=no
define symbol filename=compare_ACT4trans

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

use agu06_5d_19880101_20071231_ACT4transport.nc
use 1_INALT01-KJD308_5d_19880101_20071231_ACT4transport.nc
use "../ACT/ACTtransport_insitu.nc"
use "../ACT/ACTtransport_proxy.nc"
use "../ACT/ACT-auc_2010-2013-5d_full_gridded_ts.nc"

let cjet=calcjet455/1e6
let cbox=calcbox260/1e6
let cauc=auc/1e6
let cauc1000=auc1000/1e6

define axis/t="01-jan-1988":"31-dec-2013":1/units="days"/edges/cal=noleap t1day
let ACTbox1d=ACTbox[gt=t1day]/1e6
let ACTjet1d=ACTjet[gt=t1day]/1e6

define axis/t="01-jan-1988":"31-dec-2013":5/units="days"/edges/cal=noleap t5day
let ACTbox5d=ACTbox1d[gt=t5day@ave]
let ACTjet5d=ACTjet1d[gt=t5day@ave]

define axis/t="01-jan-1988":"31-dec-2013":1/units="months"/edges/cal=noleap t1mon
let ACTbox1m=ACTbox1d[gt=t1mon@ave]
let ACTjet1m=ACTjet1d[gt=t1mon@ave]
let cbox1m=cbox[gt=t1mon@ave]
let cjet1m=cjet[gt=t1mon@ave]
let caucm=cauc[gt=t1mon@ave]
let cauc1000m=cauc1000[gt=t1mon@ave]

define axis/t="01-jan-1988":"31-dec-2013":1/units="years"/edges/cal=noleap t1yr
let ACTjet1y=ACTjet1d[gt=t1yr@ave]
let cbox1y=cbox[gt=t1yr@ave]
let cjet1y=cjet[gt=t1yr@ave]

set view up

plot/nol/d=1/thick=1/col=black/vlim=-120:20/hlim="01-jan-1988":"31-dec-2013"/set($op) cbox
ppl axlabp,-1,-1
ppl txlabp,-1,-1
ppl ylab Sv
ppl plot
plot/nol/d=1/thick=2/col=black/ov cbox1m[l=@shn:23]
let cbox=calcbox225/1e6
plot/nol/d=2/thick=1/col=red/ov($op) cbox
plot/nol/d=2/thick=2/col=red/ov cbox1m[l=@shn:23]
plot/nol/d=3/thick=1/col=blue/ov($op) ACTbox5d
plot/nol/d=4/thick=2/col=blue/ov ACTbox1m[l=@shn:23]

plot/nol/d=1/thick=1/col=black/ov($op) cauc
plot/nol/d=1/thick=2/col=black/ov caucm[l=@shn:23]
!plot/nol/d=1/thick=2/col=black/dashed/ov cauc1000m[l=@shn:23]
plot/nol/d=2/thick=1/col=red/ov($op) cauc
plot/nol/d=2/thick=2/col=red/ov caucm[l=@shn:23]
!plot/nol/d=2/thick=2/col=red/dashed/ov cauc1000m[l=@shn:23]
let aucm=auc[gt=t1mon@ave]
plot/nol/d=5/thick=1/col=blue/ov($op) aucm
label  -25000,30,-1,0,.18,(a) Tbox

set view lo

plot/nol/d=1/thick=1/col=black/vlim=-120:0/hlim="01-jan-1988":"31-dec-2013"/set($op) cjet
ppl axlabp,-1,-1
ppl txlabp,-1,-1
ppl ylab Sv
ppl plot
plot/nol/d=1/thick=2/col=black/ov cjet1m[l=@shn:23]
let cjet=calcjet400/1e6
plot/nol/d=2/thick=1/col=red/ov($op) cjet
plot/nol/d=2/thick=2/col=red/ov cjet1m[l=@shn:23]
plot/nol/d=3/thick=1/col=blue/ov($op) ACTjet5d
plot/nol/d=4/thick=2/col=blue/ov ACTjet1m[l=@shn:23]
label  -25000,10,-1,0,.18,(b) Tjet

!let cbox=calcbox250/1e6
!list/d=1 cjet[l=@ave],cjet[l=@var]^0.5,cbox[l=@ave],cbox[l=@var]^0.5,cauc[l=@ave],cauc[l=@var]^0.5,cauc1000[l=@ave],cauc1000[l=@var]^0.5
!let cbox=calcbox/1e6
!list/d=2 cjet[l=@ave],cjet[l=@var]^0.5,cbox[l=@ave],cbox[l=@var]^0.5,cauc[l=@ave],cauc[l=@var]^0.5,cauc1000[l=@ave],cauc1000[l=@var]^0.5

list/d=3 ACTjet5d[l=@ave],ACTjet5d[l=@var]^0.5,ACTbox5d[l=@ave],ACTbox5d[l=@var]^0.5

let p=cbox1m[d=1,l=@shn:23]
let q=cbox1m[d=2,l=@shn:23]
go regresst
list rsquare

let p=cjet1m[d=1,l=@shn:23]
let q=cjet1m[d=2,l=@shn:23]
list rsquare


cancel data/all

use "../4-5-psi/agu06_sbarotrob_1951_2007_glob_mm_INALT01.nc"
use "../4-5-psi/INALT01-KJD308_1m_19480101_20071231_grid_U_psi.nc"
use "../1-resolution/ORCA05-KAB109_mesh_mask.nc"
use "../1-resolution/1_mesh_mask_B8N.nc"

define axis/t="01-jan-1948":"31-dec-2007":1/units="months"/edges/cal=noleap tmon
let/units="Sv" psi_inalt=psi[d=2,gt=tmon@asn]

define axis/t="01-jan-1951":"31-dec-2007":1/units="months"/edges/cal=noleap tmon2
let/units="Sv" psi_fesom=psi[d=1,gt=tmon2@asn]

let/units=degrees_E lamt=glamt[d=3]
let/units=degrees_N phit=gphit[d=3]

cancel region
!set region/i=612:655/j=170:198
let mske=0*tmask[d=3,k=1,i=612:655,j=170:198]+1-tmask[d=3,k=20,i=612:655@shf:-7,j=170:198]
let msks=0*tmask[d=3,k=1,i=612:655,j=170:198]+1-tmask[d=3,k=20,i=612:655,j=170:198@shf:7]
let msk= if mske eq 1 or msks eq 1 then 1 else 0
!shade msk
!shade/nol/pal=black/lev=(0,0.1,0.1)/d=3/ov tmask[k=1]

let psitmp_inalt=ignore0(psi_inalt[i=612:655,j=170:198]*msk[l=1])
let psimax_inalt=psitmp_inalt[i=@min]
let psitmp_fesom=ignore0(psi_fesom[i=612:655,j=170:198]*msk[l=1])
let psimax_fesom=psitmp_fesom[i=@min]
let ylatdat=GPHIT[d=3,j=170:198,i=615,l=1]
!fill/line/key psitmp_inalt[l=@ave]
save/file=tmp_psimax.nc/clobber psimax_inalt,psimax_fesom,ylatdat
cancel var psimax_inalt psimax_fesom ylatdat
use tmp_psimax.nc

define axis/y/units=degrees_N ylat=ylatdat
define axis/t="01-jan-1988":"31-dec-2007":1/units="years"/edges/cal=noleap tyear
let fesom=psimax_fesom[gy=ylat@asn,gt=tyear@ave]
let inalt=psimax_inalt[gy=ylat@asn,gt=tyear@ave]
let fstd=fesom[l=@var]^0.5
let fmean=fesom[t="01-jan-1988":"31-dec-2007"@ave]
let istd=inalt[l=@var]^0.5
let imean=inalt[t="01-jan-1988":"31-dec-2007"@ave]

set view ri
plot/trans/vlim=37S:28S/hlim=-30:-100/thick=3/nol/set fmean
ppl ylab Latitude
ppl xlab Sv
ppl axlabp -1,1
ppl axlint,4,1
ppl plot
plot/trans/nol/thick=1/col=black/dashed/ov fmean-fstd
plot/trans/nol/thick=1/col=black/dashed/ov fmean+fstd
!go fill_between poly/over/palette=lightgray/trans/nol fmean-fstd ylat fmean+fstd
plot/trans/nol/thick=3/col=red/ov imean
plot/trans/nol/thick=1/col=red/dashed/ov imean-istd
plot/trans/nol/thick=1/col=red/dashed/ov imean+istd

!plot/trans/nol/thick=2/col=black/dashed/ov fmean
!plot/trans/nol/thick=2/col=red/dashed/ov imean

!label  -85.6,-35.7,0,0,.18,@C004X
label  -37,-36.7,0,0,.18,(c)

IF ($program_name"0|PyFerret>1|*>0") THEN 
   label -70,-31,-1,,0.14 FESOM
   SET TEXT/FONT=arial/col=red
   label -70,-31.5,-1,,0.14 NEMO
   CANCEL TEXT/FONT=arial/col=red
ENDIf



list/y=35S fmean,imean

IF ($program_name"0|PyFerret>1|*>0") THEN 
   SET TEXT/FONT=arial/col=white
   CANCEL TEXT/FONT=arial/col=white
ENDIf


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