!
! Ferret/PyFerret Script to plot Figure 7
!
! 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 is extracted beforehand (see extract scripts)

go default
ppl conset ,2
Set mem/size=500

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


define view/x=0.00:0.43/text=0.8 xeins
define view/x=0.28:0.71/text=0.8 xzwei
define view/x=0.56:0.99/text=0.8 xdrei

!define view/x=0.00:0.44/y=0:0.75/text=0.8 xeins
!define view/x=0.275:0.715/y=0:0.75/text=0.8 xzwei
!define view/x=0.55:0.99/y=0:0.75/text=0.8 xdrei

define view/x=0:0.99/y=0.55:1/text=0.8 up


define symbol print=no
define symbol filename=compare_ACT4_1988-2007

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
ELSE
   if ($print) then set mode metafile ($filename).plt endif
ENDIf

define axis/t="01-jan-1988":"31-dec-2007":1/edges/unit="months"/cal=noleap tmon


define region/t="01-jan-1988":"31-dec-2007" tim
set region tim

!------ FESOM data -----------------

use agu06_5d_19880101_20071231_ACT4.nc

set win/asp=0.7

set view xeins
fill/pal=lightgray/lev=(0,400,400)/nokey/nol/hlim=0:280/vlim=5000:0/set vcross[l=@ave]
ppl axlabp,-1,-1
ppl xlab km
ppl ylab depth (m)
ppl fill
contour/ov/nol/col=black/lev="(-20,0,5)(0,20,5) line(-20,-5,5)" vcross[l=@ave]
contour/ov/nol/col=blue/lev="(-200,-30,10) line(-200,-30,10)" vcross[l=@ave]
contour/ov/nol/col=red/lev="(1,200,1)" vcross[x=0:150,l=@ave]
cancel region
go fill_between poly/over/palette=black/nol bat[l=1] xax bat[l=1]*0+5000
set region tim

!------ INALT01 data -----------------

use 1_INALT01-KJD308_5d_19880101_20071231_ACT4.nc


set view xzwei
fill/pal=lightgray/lev=(0,400,400)/nokey/nol/hlim=0:280/vlim=5000:0/set vcross[l=@ave]
ppl axlabp,-1,0
ppl xlab km
ppl ylab depth (m)
ppl fill
contour/ov/nol/col=black/lev="(-20,0,5)(0,20,5) line(-20,-5,5)" vcross[l=@ave]
contour/ov/nol/col=blue/lev="(-200,-30,10) line(-200,-30,10)" vcross[l=@ave]
contour/ov/nol/col=red/lev="(1,200,1)" vcross[x=0:150,l=@ave]
cancel region
go fill_between poly/over/palette=black/nol bat[l=1] xax bat[l=1]*0+5000
set region tim



!----- calculate and plot transports -----

! calculate transports

let velN = if vcrosstrans gt 0 then vcrosstrans
let velS = if vcrosstrans lt 0 then vcrosstrans

let AC= velS[x=@din,     k=@sum,gt=tmon@ave]/1e6
let AUC=velN[x=0:200@din,k=@sum,gt=tmon@ave]/1e6

!let/title="AC ACT Barotropic Transport (Sv)"           AC_baro_full = vcrosstrans[x=@din,k=@sum]/1e6
!let/title="AC ACT Barotropic Northward Transport (Sv)" AC_baro_Northward = locN[x=@din,k=@sum]/1e6
!let/title="AC ACT Barotropic Southward Transport (Sv)" AC_baro_Southward = locS

!let/title="AC ACT 0-250Km Northward Transport (Sv)" AC_Northward = 
!let/title="AC ACT 0-250Km Southward Transport (Sv)" AC_Southward = locS[x=0:280@din,k=@sum]/1e6


!set view up

!plot/nol/d=1/thick=2/col=black/vlim=-120:40/hlim="01-jan-1988":"31-dec-2014"/set AC,AUC
!ppl axlabp,-1,-1
!ppl txlabp,-1,-1
!ppl ylab Sv
!ppl plot
!plot/nol/d=2/thick=2/col=red/ov AC,AUC

!list/d=1 AC[l=@ave], AUC[l=@ave]
!list/d=2 AC[l=@ave], AUC[l=@ave]

!list/d=1 AC_baro_full[l=@ave],AC_baro_Northward[l=@ave],AC_baro_Southward[l=@ave],AC_Northward[l=@ave],AC_Southward[l=@ave]
cancel region

!------ ACT data -----------------

set view xdrei
!set mem/size=1000
!use/order=zxt "../ACT/ACT_2010-2013_full_gridded_ts.nc"
!define axis/x/units="km" dist=X[gx=distance]/1000
!let vcross=v[gx=dist@asn]*100

!let vave=vcross[z=0:50@ave]
!let vmin=vave[x=@min]
!let vmin0=vave-vmin
!let vcore=vmin0[x=@loc:0]

!Let vcoremean=vcore[t=@ave]
!Let vcoresdev=vcore[t=@var]^0.5

!let vcross_normal=if vcore lt vcoremean+vcoresdev then vcross
!let vcross_extrem=if vcore gt vcoremean+2*vcoresdev then vcross

!let act_normal=vcross_normal[t=@ave]
!let act_extrem=vcross_extrem[t=@ave]
let act_mean=vcross[t=@ave]

!save/file="../ACT/ACT_2010-2013-mean_full_gridded_ts.nc"/clobber act_mean,act_normal,act_extrem
!cancel var act_normal,act_extrem; cancel data/all
!define axis/t="01-may-2010":"31-jan-2013":1/unit="months" tmon
!let vcross=v[gx=dist@asn,gt=tmon@ave]*100
!save/file="../ACT/ACT_2010-2013-1m_full_gridded_ts.nc" vcross
use "../ACT/ACT_2010-2013-1m_full_gridded_ts.nc"
let act_mean=vcross[l=@ave]

fill/pal=lightgray/lev=(0,400,400)/nokey/nol/hlim=0:280/vlim=0:5000/set act_mean
ppl axlabp,-1,1
ppl xlab km
ppl ylab depth (m)
ppl fill
contour/ov/nol/col=black/lev="(-20,0,5)(0,20,5) line(-20,-10,5)" act_mean
contour/ov/nol/col=blue/lev="(-200,-30,10) line(-200,-30,10)" act_mean
contour/ov/nol/col=red/lev="(1,200,1)" act_mean[x=0:150]
contour/ov/nol/col=green/lev="(-5)" act_mean

set view xeins
contour/ov/nol/col=green/lev="(-5)" act_mean
set view xzwei
contour/ov/nol/col=green/lev="(-5)" act_mean
set view xdrei


use etopo5
let bathy=rose*(-1)
let seclim="27.5,-33.3,28.9,-35.7,100"
go vertical_section_bat  black rose*(-1) `seclim`
!label 20,4800,-1,0,.18,@C006ACT (2010-2013)


IF ($program_name"0|PyFerret>1|*>0") THEN 
   SET TEXT/FONT=arial/col=white
   label 20,4800,-1,0,.14, (c) ACT !(2010-2013)
   Set view xzwei
   label 20,4800,-1,0,.14, (b) NEMO !(1990-1999)
   Set view xeins
   label 20,4800,-1,0,.14, (a) FESOM !(1990-1999)
   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