!
! Ferret/PyFerret Script to plot Figure 10
!
! Biastoch, Sein, Durgadoo, Wang, Danilov
! Modelling the eddying Agulhas system – 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

cancel region
use Ariane_INALT01-KJD308_5d_1948-2005.nc
use Ariane_FESOM_9y.nc
show data

define symbol print=no
define symbol filename=compare_FESOM-NEMO

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

set win/asp=1
let/title="NEMO"   NEMO=ignore0(AL[d=1])
let/title="FESOM" FESOM=ignore0(WEST[d=2]+NW[d=2])
plot/vlim=0:75/thick=2/nol/set FESOM,NEMO
ppl ylab Sv
ppl plot
plot/ov/nol/col=1/thick=2 EAST[d=2]
plot/ov/nol/col=2/thick=2 EAST[d=1]
!plot/ov/nol/col=2/thick=1 LOST[d=1]
!plot/ov/nol/col=1/thick=1 LOST[d=2]
plot/ov/nol/col=2/thick=2 EXCEPTMNDS[d=1]
plot/ov/nol/col=1/thick=2 MNDS[d=2]

set region/t="01-jan-1965":"31-dec-2000"
let p=T[gt=NEMO]
let q=NEMO
go regresst
list slope*10,q[l=@ave]
plot/ov/nol/dash/thick=1/col=1 qhat
label/user 200000,22.5,-1,0,.18,Agulhas Leakage
let AL_INALT01=q-qhat
save/file=tmp.nc/clobber AL_INALT01

let q=EAST[d=1]
list slope*10
plot/ov/nol/dash/thick=1/col=1 qhat
label/user 200000,51,-1,0,.18,Agulhas Return Current


let q= EXCEPTMNDS[d=1]
list slope*10
plot/ov/nol/dash/thick=1/col=1 qhat
label/user 200000,70,-1,0,.18,Agulhas Current
let AC_INALT01=q-qhat
save/file=tmp.nc/app AC_INALT01

let p=T[gt=FESOM]
let q=FESOM
list slope*10,q[l=@ave]
plot/ov/nol/dash/thick=1/col=2 qhat
let AL_FESOM=q-qhat
save/file=tmp.nc/app AL_FESOM


let q=EAST[d=2]
list slope*10
plot/ov/nol/dash/thick=1/col=2 qhat

let q= MNDS[d=2]
list slope*10
plot/ov/nol/dash/thick=1/col=2 qhat
let AC_FESOM=q-qhat
save/file=tmp.nc/append AC_FESOM

!label/user 190000,5,-1,0,.14,Lost


!cancel region/t
!use Ariane_FESOM_9y.nc

!let/title="FESOM" FESOM=ignore0(WEST[d=3]+NW[d=3])
!plot/ov/thick=1/col=4/nol FESOM
!plot/ov/nol/col=4/thick=1 EAST[d=3]
!plot/ov/nol/col=4/thick=1 LOST[d=3]
!plot/ov/nol/col=4/thick=1 MNDS[d=3]

!set region/t="01-jan-1965":"31-dec-2000"
!let p=T[gt=FESOM]
!let q=FESOM
!list slope*10,q[l=@ave]
!plot/ov/nol/dash/thick=1/col=4 qhat

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