! Aeolus water transport: large, synoptic and meso scale
!
! Arguments
! $1 = Season index: 1=DJF, 2=MAM, 3=JJA, 4=SON, ANNUAL=1:4@ave
! $2 = Season: DJF, MAM, JJA, SON, ANNUAL
! $3 = dataset

cancel data/all

!spawn cdo -seldate\,0052-01-01T00:00:00\,0104-12-31T23:59:59 $3 Watert-data.nc
!spawn ncatted -a axis\,level\,c\,c\,Z Watert-data.nc
!spawn cdo yseasmean Watert-data.nc Watert-data_sm.nc
!spawn cdo settaxis,2013-02-01,00:00,3months Watert-data_sm.nc Watert-data_sm_taxis.nc 
!spawn ncks -A -v levelb $3 Watert-data_sm_taxis.nc
!spawn ncatted -a edges\,level\,c\,c\,levelb Watert-data_sm_taxis.nc
!spawn rm -f Watert-data.nc Watert-data_sm.nc
!use "Watert-data_sm_taxis.nc" !d=1 AEOLUS DATA
use $3 !d=1 AEOLUS DATA

set memory/size=10000

if `"($PROGRAM_NAME)" eq "Ferret"` then
set mode metafile: watertransport.plt
endif


!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!! Make the plots! !!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!

let cp = 1.0e3
let rho = 1.3
let H0 = 8000

let kk=k[d=1,gz=LARGE_TRANSPORT_QY]
let kmax=kk[k=@max]
! `kmax-1` is top-most troposphere level
let zz=z[d=1,gz=LARGE_TRANSPORT_QY]
let p_s=1000 ! surface pressure
let rho_z=rho*zz[k=1:`kmax-1`]/p_s ! air density

set win 1
! meridional water transport
set viewport UL
! Sigh. This tries to compute air density based on meter heights.
! But we already converted the data set z axis to pressure heights.
!shade/nolabels/levels=20 rho*exp(-1*z[k=1:5]/H0)*LARGE_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:5]
shade/nolabels/levels=20 rho_z*LARGE_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Large scale meridional water transport (kg/m^2/s)"
PPL SHASET RESET

set viewport UR
shade/nolabels/levels=20 rho_z*SYNOP_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Synoptic scale meridionnal water transport (kg/m^2/s)"
!---timestemp--------------------------------------------------
let xpos=`($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.1 "@ASPIK/L.Caesar<nl>($SESSION_DATE)"
PPL SHASET RESET

set viewport LL
shade/nolabels/levels=20 rho_z*MESO_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Meso scale meridional water transport (kg/m^2/s)"
PPL SHASET RESET

set viewport LR
!shade/nolabels/levels=20 rho*exp(-1*z[k=1:5]/H0)*LARGE_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:5]+rho*exp(-1*z[k=1:5]/H0)*SYNOP_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:5]+rho*exp(-1*z[k=1:5]/H0)*MESO_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:5]
shade/nolabels/levels=20 rho_z*LARGE_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]+rho_z*SYNOP_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]+rho_z*MESO_TRANSPORT_QY[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Total meridional water transport (kg/m^2/s)"
PPL SHASET RESET

if `"($PROGRAM_NAME)" eq "PyFerret"` then
!frame/file="watertransport_1.ps" Sigh. Postscript export is broken with Qt5
frame/file="watertransport_1.pdf"
spawn pdf2ps watertransport_1.pdf watertransport_1.ps
spawn rm watertransport_1.pdf
endif
!!

set win 2
! zonal water transport
set viewport UL
shade/nolabels/levels=20 rho_z*LARGE_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Large scale zonal water transport (kg/m^2/s)"
PPL SHASET RESET

set viewport UR
shade/nolabels/levels=20 rho_z*SYNOP_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Synoptic scale zonal water transport (kg/m^2/s)"
!---timestemp--------------------------------------------------
let xpos=`($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.1 "@ASPIK/L.Caesar<nl>($SESSION_DATE)"
PPL SHASET RESET

set viewport LL
shade/nolabels/levels=20 rho_z*MESO_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Meso scale zonal water transport (kg/m^2/s)"
PPL SHASET RESET

set viewport LR
shade/nolabels/levels=20 rho_z*LARGE_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]+rho_z*SYNOP_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]+rho_z*MESO_TRANSPORT_QX[d=1,l=$1,i=@ave,k=1:`kmax-1`]
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Latitude"
!---ylabel-------------------------------------------------
let xpos=`-0.7*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Height (mbar)"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Aeolus $2 zonal mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Total zonal water transport (kg/m^2/s)"
PPL SHASET RESET

if `"($PROGRAM_NAME)" eq "PyFerret"` then
!frame/file="watertransport_2.ps" Sigh. Postscript export is broken with Qt5
frame/file="watertransport_2.pdf"
spawn pdf2ps watertransport_2.pdf watertransport_2.ps
spawn rm watertransport_2.pdf
endif


!!!!! The end !!!!!
if `"($PROGRAM_NAME)" eq "Ferret"` then
  CANCEL MODE METAFILE
  spawn Fprint -l cps -o watertransport_1.ps -R watertransport.plt
  spawn Fprint -l cps -o watertransport.ps -R watertransport.plt.~*~
  spawn cat watertransport_1.ps >> watertransport.ps
else
  spawn cat watertransport_?.ps > watertransport.ps
endif
spawn rm watertransport_?.ps
spawn ps2pdf watertransport.ps watertransport_seas_$2.pdf
spawn rm watertransport.plt
spawn rm watertransport.plt*
spawn rm ferret.jnl
spawn rm watertransport.ps
!cancel data/all
!spawn rm -f Watert-data_sm_taxis.nc