! Aeolus output: large, synoptic and meso scale energy transport
!
! 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 Energytransport-data.nc
!spawn ncatted -a axis\,level\,c\,c\,Z Energytransport-data.nc
!spawn cdo yseasmean Energytransport-data.nc Energytransport-data_sm.nc
!spawn cdo settaxis,2013-02-01,00:00,3months Energytransport-data_sm.nc Energytransport-data_sm_taxis.nc 
!spawn ncks -A -v levelb $3 Energytransport-data_sm_taxis.nc
!spawn ncatted -a edges\,level\,c\,c\,levelb Energytransport-data_sm_taxis.nc
!spawn rm -f Energytransport-data.nc Energytransport-data_sm.nc
!use "Energytransport-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: energytransport.plt
endif

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

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

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

set win 1
! meridional energy 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 cp*rho*exp(-1*z[k=1:5]/H0)*LARGE_TRANSPORT_TY[d=1,l=$1,i=@ave,k=1:5]
shade/nolabels/levels=20 cp*rho_z*LARGE_TRANSPORT_TY[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 energy transport (J/m^2/s)"
PPL SHASET RESET

set viewport UR
shade/nolabels/levels=20 cp*rho_z*SYNOP_TRANSPORT_TY[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 energy transport (J/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 cp*rho_z*MESO_TRANSPORT_TY[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 energy transport (J/m^2/s)"
PPL SHASET RESET

set viewport LR
shade/nolabels/levels=20 cp*rho_z*LARGE_TRANSPORT_TY[d=1,l=$1,i=@ave,k=1:`kmax-1`]+cp*rho_z*SYNOP_TRANSPORT_TY[d=1,l=$1,i=@ave,k=1:`kmax-1`]+cp*rho_z*MESO_TRANSPORT_TY[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 energy transport (J/m^2/s)"
PPL SHASET RESET

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

!!

set win 2
! zonal energy transport
set viewport UL
shade/nolabels/levels=20 cp*rho_z*LARGE_TRANSPORT_TX[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 energy transport (J/m^2/s)"
PPL SHASET RESET

set viewport UR
shade/nolabels/levels=20 cp*rho_z*SYNOP_TRANSPORT_TX[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 energy transport (J/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 cp*rho_z*MESO_TRANSPORT_TX[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 energy transport (J/m^2/s)"
PPL SHASET RESET

set viewport LR
shade/nolabels/levels=20 cp*rho_z*LARGE_TRANSPORT_TX[d=1,l=$1,i=@ave,k=1:`kmax-1`]+cp*rho_z*SYNOP_TRANSPORT_TX[d=1,l=$1,i=@ave,k=1:`kmax-1`]+cp*rho_z*MESO_TRANSPORT_TX[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 energy transport (J/m^2/s)"
PPL SHASET RESET

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


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