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

go default
ppl conset ,2

cancel region; cancel view; cancel data/all;can mode logo
define symbol print=no
define symbol filename=moc

use "../1-resolution/ORCA05-KAB109_mesh_mask.nc"
use INALT01-KJD308_1m_19480101_20071231_grid_V_amoc.nc
use 1_INALT01-KJD308_1m_19480101_20071231_grid_V_amoc.nc

!! Compute time-mean over 60 years
define axis/t="01-jan-1948":"31-dec-2007":1/units="months"/cal=noleap/edges tinalt
let mocb2 = zomsfatl[d=2,j=1:510,gt=tinalt@asn] 
!the latitude for jmax = 511 is left out, because it is south of j=510
let mocn2 = zomsfatl[d=3,gt=tinalt@asn]

!! Flip base model amoc variable vertically
! the cdfamoc as well as the cdfmoc function produce amocs for the base model and the nest that are flipped vertically with respect to each other

def ax/z/unit=meter dep=nav_lev[d=1]
let ksz      = `mocn2,r=ksize`
let z_rev    = SAMPLEK(z[gz=dep],ksz-K[gz=dep]+1)
let mocn2_rev = SAMPLEK(mocn2,ksz-K[gz=dep]+1)
let mocn2_revz = mocn2_rev[gz=dep@asn]

!! For Flipping base model amoc:
!let z_rev    = SAMPLEK(z[gz=dep],ksz-K[gz=dep]+1)
!def ax/z/from_data/units=meters zdepr=z_rev*-1
!let mocb2_revz = mocb2_rev[gz=zdepr@asn]

! Define grid for latitude
def ax/y/unit=degree_n latb=nav_lat[d=2,j=1:510]
def ax/y/unit=degree_n latn=nav_lat[d=3]
let mocb = mocb2[gz=dep@asn,gy=latb@asn]
let mocn = mocn2_revz[gy=latn@asn]


use/order=jk agu06_AOVERT_1980_2009_mm_4b_mod.nc
set axis/units=degrees_N lat

IF ($program_name"0|PyFerret>1|*>0") THEN 
   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


define axis/t="01-jan-1980":"31-dec-2009":1/units="months"/cal=noleap/edges tfesom
!def ax/z/from_data/units=meters/depth zdepf=Z[g=goy1]*-1
!let fesom=G_OVERT[gt=tfesom@asn,gz=zdepf@asn,d=4]
let fesom=G_OVERT[gt=tfesom@asn,d=4]

define view/text=0.7/y=0.65:1.0/x=0:0.7 upup
define view/text=0.7/y=0.4:0.915/x=0:0.7 uplo
define view/text=0.7/y=0.25:0.6/x=0:0.7 loup
define view/text=0.7/y=0.0:0.515/x=0:0.7 lolo
define view/text=0.7/y=0:1/x=0.6:1.0 le



define view/text=0.7/y=0.7:1.0/x=0:0.7 upup
define view/text=0.7/y=0.5:0.94/x=0:0.7 uplo
define view/text=0.7/y=0.4:0.7/x=0:0.7 loup
define view/text=0.7/y=0.2:0.64/x=0:0.7 lolo
define view/text=0.7/y=0:1/x=0.55:0.95 le
define view/text=0.7/y=0:0.35/x=0:0.7 lo


set win/asp=1
set region/t="01-jan-1988":"31-dec-2007"
let diff=ignore0(fesom[t=@ave])
set view upup
fill/lev=(-20,0,20)/hlim=-35:65/vlim=200:0/pal=lightgray/nokey/nol/set diff
ppl ylab 
ppl xlab
ppl axset 1 0 1 1
ppl axlabp 0,-1
ppl axlint,2,5
ppl fill
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=10/size=0.08/ov diff
label -35 -50 -1 0 0.18 (a) FESOM
set view uplo
fill/lev=(-20,0,20)/hlim=-35:65/vlim=5500:200/pal=lightgray/nokey/nol/set diff
ppl ylab Depth (m)
ppl xlab 
ppl axlabp 0,-1
ppl axset 0 1 1 1
ppl axlint,2,2
ppl fill
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=10/size=0.08/ov diff

let diffn=ignore0(mocn[j=3:657,t=@ave])
let diffb=ignore0(mocb[j=268:510,t=@ave])
set view loup
fill/lev=(-20,0,20)/hlim=-35:65/vlim=200:0/pal=lightgray/nokey/nol/set diffb
ppl ylab 
ppl xlab 
ppl axset 1 0 1 1
ppl axlabp 0,-1
ppl axlint,2,5
ppl fill
fill/lev=(-20,0,20)/pal=lightgray/nokey/nol/ov diffn
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=1/size=0.08/col=black/ov diffb
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=1/size=0.08/col=black/ov diffn
label -35 -50 -1 0 0.18 (b) NEMO
set view lolo
fill/lev=(-20,0,20)/hlim=-35:65/vlim=5500:200/pal=lightgray/nokey/nol/set diffb
ppl ylab Depth (m)
ppl xlab Latitude
ppl axset 0 1 1 1
ppl axlabp -1,-1
ppl axlint,2,2
ppl fill
fill/lev=(-20,0,20)/pal=lightgray/nokey/nol/ov diffn
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=1/size=0.08/col=black/ov diffb
cont/nol/lev=(-5,-1,1)(0,30,2)/spacing=1/size=0.08/col=black/ov diffn
shade/pal=black/ov/nol/lev=(0) missing(diffb[l=1,k=@shf:1],0)
shade/pal=black/ov/nol/lev=(0) missing(diffn[l=1,k=@shf:1],0)
ppl color 6 99 99 99
!label 60 5000 0 0 0.24 @C006(b)


Set view uplo
shade/pal=black/ov/nol/lev=(0) missing(diffb[l=1,k=@shf:1],0)
shade/pal=black/ov/nol/lev=(0) missing(diffn[l=1,k=@shf:1],0)
!label 60 5000 0 0 0.24 @C006(a)

set view le
!plot/vlim=5800:0/thick=2/nol/hlim=-4:20/set fesom[t=@ave,y=26.5],mocb[t=@ave,y=26.5]
plot/vlim=5800:0/thick=2/nol/hlim=-14:20/set fesom[t=@ave,y=26.5,z=@ddc]*1e3,mocb[t=@ave,y=26.5,z=@ddc]*1e3
ppl axlabp -1,1
ppl axlint,2,2
ppl axset 1 1 1 1
ppl ylab Depth (m)
ppl xlab mSv/m
ppl plot
plot/ov/dashed/col=black/thick=1/nol fesom[t=@ave,y=26.5]*0

IF ($program_name"0|PyFerret>1|*>0") THEN 
   label 2,2500,-1,,0.14 FESOM
   SET TEXT/FONT=arial/col=red
   label 2,2750,-1,,0.14 NEMO
   SET TEXT/FONT=arial/col=blue
   label 2,3000,-1,,0.14 RAPID
   CANCEL TEXT/FONT=arial/col=blue
ENDIf


cancel region
! RAPID data
use moc_vertical.nc
!plot/ov/col=blue/thick=2/nol STREAM_FUNCTION_MAR[t=@ave]
plot/ov/col=blue/thick=2/nol STREAM_FUNCTION_MAR[t=@ave,z=@ddc]*1e3
label 15.5 5600 0 0 0.24 @C001(c)

define axis/t="01-jan-2004":"31-dec-2015":1/units="months"/cal=noleap/edges trapid
let fesom26N=fesom[y=26.5,z=500:3000@max]
let inalt26N=mocb[y=26.5,z=500:3000@max]
let rapid26N=STREAM_FUNCTION_MAR[gt=trapid,z=500:3000@max]
save/file=AMOC_26N_series.nc/clobber fesom26N,inalt26N,rapid26N
cancel var fesom26N,inalt26N,rapid26N;use AMOC_26N_series.nc

set view lo
plot/nol/thick=2/hlim="01-jan-1980":"31-dec-2007"/set fesom26N[l=@shn:23]-fesom26N[l=@ave],inalt26N[l=@shn:23]-inalt26N[l=@ave]
ppl axlabp -1,-1
ppl axlint,5,4
ppl ylab Sv
ppl plot
!plot/ov/nol RAPID26N[l=@shn:23]-RAPID26N[l=@ave]
label 295000,2.2,-1,,0.18 (d)


! calculate trends and statistics

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

!let p=T[gt=fesom26N]
!let q=fesom26N
!let fesom26Nd=fesom26N[l=@shn:23]-qhat
!go regresst
!save/file=tmp.nc/clobber fesom26Nd
!let p=T[gt=inalt26N]
!let q=inalt26N
!let inalt26Nd=inalt26N[l=@shn:23]-qhat
!save/file=tmp.nc/ap inalt26Nd
!cancel var fesom26Nd, inalt26Nd; use tmp
!let p=inalt26Nd
!let q=fesom26Nd
!go regresst
!list rsquare^0.5

define axis/t="01-jan-1980":"31-dec-2007":1/units="years"/edges/cal=noleap tyear
!let fesoma=fesom26N[gt=tyear@ave]
let fesoma=fesom26N[l=@shn:23]
let q=fesoma
let p=T[gt=q]
go regresst
let fesomd=q-qhat
save/file=tmp_moc.nc/clobber fesoma,fesomd
!let inalta=inalt26N[gt=tyear@ave]
let inalta=inalt26N[l=@shn:23]
let q=inalta
let p=T[gt=q]
let inaltd=q-qhat
save/file=tmp_moc.nc/app inalta,inaltd
cancel data/all; cancel var fesoma,fesomd,inalta,inaltd

use tmp_moc
let p=fesoma
let q=inalta
go regresst
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