!
! Ferret/PyFerret Script to plot Figure 5
!
! 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
!

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

go default
ppl pen 3 10 10 10

define symbol print=no
define symbol filename=psi_series

IF ($program_name"0|PyFerret>1|*>0") THEN 
   SET TEXT/FONT=arial
ELSE
   if ($print) then set mode metafile ($filename).plt endif
ENDIf

! transports are read from streamfunction files 
! and defined as transport differences between two locations
!
use agu06_sbarotrob_1951_2007_glob_mm_INALT01.nc
use INALT01-KJD308_1m_19480101_20071231_grid_U_psi.nc
define axis/t="01-jan-1948":"31-dec-2007":1/units="months"/edges/cal=noleap tmon
define axis/t="01-jan-1951":"31-dec-2007":1/units="months"/edges/cal=noleap tmon2
define axis/t="01-jan-1951":"31-dec-2007":1/units="years"/edges/cal=noleap tyear

! Transports are defined as streamfunction differences
let acc=psi[i=444,j=71]-psi[i=439,j=117]
let moz=psi[i=664,j=215]-psi[i=656,j=220]
let itf=psi[i=118,j=226]-psi[i=64,j=237]

set region/t="01-jan-1958":"31-dec-2007"

set win/asp=1.2


set view yeins
let fesom=itf[d=1,gt=tmon2@asn]
let inalt=itf[d=2,gt=tmon@asn]
plot/nol/set/thick=2 fesom[l=@shn:23],inalt[l=@shn:23]
ppl ylab Sv
ppl plot
label 20000,-7.5,-1,,0.18 (a) Indonesian Throughflow

! calculate statistics and detrend
!
list fesom[l=@ave],fesom[l=@var]^0.5
list inalt[l=@ave],inalt[l=@var]^0.5
!list rsquare^0.5
!let fesoma=fesom[gt=tyear@ave]
let fesoma=fesom[l=@shn:23]
let q=fesoma
let p=T[gt=q]
go regresst
let fesomd=q-qhat
save/file=tmp_itf.nc/clobber fesoma,fesomd
!let inalta=inalt[gt=tyear@ave]
let inalta=inalt[l=@shn:23]
let q=inalta
let p=T[gt=q]
let inaltd=q-qhat
save/file=tmp_itf.nc/app inalta,inaltd


let p=fesom[gt=tyear@ave]
let q=inalt[gt=tyear@ave]

set view yzwei
let fesom=moz[d=1,gt=tmon2@asn]
let inalt=moz[d=2,gt=tmon@asn]
plot/nol/thick=2/set fesom[l=@shn:23],inalt[l=@shn:23]
ppl ylab Sv
ppl plot
label 20000,-3.5,-1,,0.18 (b) Mozambique Channel transport

list fesom[l=@ave],fesom[l=@var]^0.5
list inalt[l=@ave],inalt[l=@var]^0.5
!list rsquare^0.5
!let fesoma=fesom[gt=tyear@ave]
let fesoma=fesom[l=@shn:23]
let q=fesoma
let p=T[gt=q]
go regresst
let fesomd=q-qhat
save/file=tmp_moz.nc/clobber fesoma,fesomd
!let inalta=inalt[gt=tyear@ave]
let inalta=inalt[l=@shn:23]
let q=inalta
let p=T[gt=q]
let inaltd=q-qhat
save/file=tmp_moz.nc/app inalta,inaltd



set view ydrei
let fesom=acc[d=1,gt=tmon2@asn]
let inalt=acc[d=2,gt=tmon@asn]
plot/nol/thick=2/set fesom[l=@shn:23],inalt[l=@shn:23]
ppl ylab Sv
ppl plot
label 20000,160,-1,,0.18 (c) Drake Passage Transport

IF ($program_name"0|PyFerret>1|*>0") THEN 
   label 350000,142,-1,,0.18 FESOM
   SET TEXT/FONT=arial/col=red
   label 350000,125,-1,,0.18 NEMO
   CANCEL TEXT/FONT=arial/col=red
ENDIf


list fesom[l=@ave],fesom[l=@var]^0.5
list inalt[l=@ave],inalt[l=@var]^0.5
!list rsquare^0.5
!let fesoma=fesom[gt=tyear@ave]
let fesoma=fesom[l=@shn:23]
let q=fesoma
let p=T[gt=q]
go regresst
let fesomd=q-qhat
save/file=tmp_acc.nc/clobber fesoma,fesomd
!let inalta=inalt[gt=tyear@ave]
let inalta=inalt[l=@shn:23]
let q=inalta
let p=T[gt=q]
let inaltd=q-qhat
save/file=tmp_acc.nc/app inalta,inaltd


! calculate temporal correlations between timeseries
!
cancel data/all; cancel var fesoma,fesomd,inalta,inaltd
use tmp_itf
let p=fesoma
let q=inalta
go regresst
list rsquare^0.5
let p=fesomd
let q=inaltd
list rsquare^0.5

use tmp_moz
let p=fesoma
let q=inalta
list rsquare^0.5
let p=fesomd
let q=inaltd
list rsquare^0.5

use tmp_acc
let p=fesoma
let q=inalta
list rsquare^0.5
let p=fesomd
let q=inaltd
list rsquare^0.5

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