! Comparing Era Interim Data to Coarse Bling and Aeolus output: Albedo direct, Albedo diffusive, sea ice cover
!
! Arguments
! $1 = Season index: 1=DJF, 2=MAM, 3=JJA, 4=SON, 1:4@ave=ANNUAL
! $2 = Season: DJF, MAM, JJA, SON, ANNUAL
! $3 = dataset aeolus model
! $4 = dataset atmos model down
! $5 = dataset flux

cancel data/all

!spawn cdo -seldate\,0052-01-01T00:00:00\,0104-12-31T23:59:59 $3 Albedo-data.nc
!spawn ncatted -a axis\,level\,c\,c\,Z Albedo-data.nc
!spawn cdo yseasmean Albedo-data.nc Albedo-data_sm.nc
!! set taxis, such that Aeolus and Era have same time axis for the climatologies and can be compared.
!spawn cdo settaxis\,2013-02-01\,00:00\,3months Albedo-data_sm.nc Albedo-data_sm_taxis.nc
!! cdo has omitted level bounds, extract them from the original input and append to private data file.
!spawn ncks -A -v levelb $3 Albedo-data_sm_taxis.nc
!spawn ncatted -a edges\,level\,c\,c\,levelb Albedo-data_sm_taxis.nc
!spawn rm -f Albedo-data.nc Albedo-data_sm.nc
!use "Albedo-data_sm_taxis.nc" !d=1 AEOLUS DATA
use $3 !d=1 AEOLUS DATA
use "/p/projects/climber3/caesar/CM2M/multi_year_seasonal_means/00530101.atmos-dly_sm.nc" !d=2 CM2M DATA
use "/p/projects/climber3/caesar/CM2M/multi_year_seasonal_means/00530101.flux-dly_sm.nc" !d=3 CM2M DATA
use "/p/projects/climber3/caesar/Era_interim_data/noleap_data/multi_year_season_mean/albedo_sm.nc" !d=4
use "/p/projects/climber3/caesar/CM2M/multi_year_seasonal_means/00530101.atmos_model_down-dly_sm.nc" !d=5 CM2M DATA

!spawn cdo -seldate\,0052-01-01T00:00:00\,0104-12-31T23:59:59 $4 Albedo-data_down.nc
!spawn ncatted -a axis\,level\,c\,c\,Z Albedo-data_down.nc
!spawn cdo yseasmean Albedo-data_down.nc Albedo-data_down_sm.nc
!spawn cdo settaxis\,2013-02-01\,12:00\,3months Albedo-data_down_sm.nc Albedo-data_down_sm_taxis.nc 
!spawn ncks -A -v levelb $4 Albedo-data_down_sm_taxis.nc
!spawn ncatted -a edges\,level\,c\,c\,levelb Albedo-data_down_sm_taxis.nc
!spawn rm -f Albedo-data_down.nc Albedo-data_down_sm.nc
!use "Albedo-data_down_sm_taxis.nc" !d=6 ATMOS DOWN DATA
use $4 !d=6 ATMOS DOWN DATA

use "/p/projects/climber3/caesar/Era_interim_data/noleap_data/multi_year_season_mean/sea_ice_cover_daily_mean_monthly_means_sm.nc" !d=7

!spawn cdo -seldate\,0052-01-01T00:00:00\,0104-12-31T23:59:59 $5 Albedo-data_flux.nc
!spawn ncatted -a axis\,level\,c\,c\,Z Albedo-data_flux.nc
!spawn cdo yseasmean Albedo-data_flux.nc Albedo-data_flux_sm.nc
!spawn cdo settaxis,2013-02-01,12:00,3months Albedo-data_flux_sm.nc Albedo-data_flux_sm_taxis.nc 
!spawn ncks -A -v levelb $5 Albedo-data_flux_sm_taxis.nc
!spawn ncatted -a edges\,level\,c\,c\,levelb Albedo-data_flux_sm_taxis.nc
!spawn rm -f Albedo-data_flux.nc Albedo-data_flux_sm.nc
!use "Albedo-data_flux_sm_taxis.nc" !d=8 FMS FLUX DATA
use $5 !d=8 FMS FLUX DATA

set memory/size=10000

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

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

set win 1
! Albedo visible/infrared radiation clear sky/direct
set viewport UR
shade/nolabels/levels=20 ALBEDO_VIS_DIR[d=6,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---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 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible direct"
!---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 UL
shade/nolabels/levels=(-inf)(`($lev_min)+($lev_del)`,`($lev_max)-($lev_del)`,($lev_del))(inf) AL[d=4,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Era Interim $2 climatology (1979-2013)"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Surface albedo"
PPL SHASET RESET

set viewport LL
shade/nolabels/levels=($lev_text) ALBEDO_VIS_DIR[d=5,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "CM2M $2 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible direct"
PPL SHASET RESET

set viewport LR
plot/nolabels AL[d=4,l=$1,i=@ave], ALBEDO_VIS_DIR[d=6,l=$1,i=@ave],ALBEDO_VIS_DIR[d=5,i=@ave,l=$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.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Albedo"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "$2 Zonal Means"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible direct"
!---key--------------------------------------------------
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.9*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P1Era Interim"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.85*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P2Aeolus"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.8*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P3CM2M"
PPL SHASET RESET

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

set win 2
! Albedo visible/infrared radiation cloudy conditions/diffusive
set viewport UR
shade/nolabels/levels=20 ALBEDO_VIS_DIF[d=6,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---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 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible diffusive"
!---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 UL
shade/nolabels/levels=(-inf)(`($lev_min)+($lev_del)`,`($lev_max)-($lev_del)`,($lev_del))(inf) AL[d=4,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Era Interim $2 climatology (1979-2013)"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Surface albedo"
PPL SHASET RESET

set viewport LL
shade/nolabels/levels=($lev_text) ALBEDO_VIS_DIF[d=5,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "CM2M $2 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible diffusive"
PPL SHASET RESET

set viewport LR
plot/nolabels AL[d=4,l=$1,i=@ave], ALBEDO_VIS_DIF[d=6,l=$1,i=@ave],ALBEDO_VIS_DIF[d=5,i=@ave,l=$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.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Albedo"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "$2 Zonal Means"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Albedo visible diffusive"
!---key--------------------------------------------------
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.9*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P1Era Interim"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.85*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P2Aeolus"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.8*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P3CM2M"
PPL SHASET RESET

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

!!

set win 3
! Sea ice cover
set viewport UR
shade/nolabels/levels=20 ICE_MASK[d=8,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---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 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Sea-ice cover (0-1)"
!---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 UL
shade/nolabels/levels=(-inf)(`($lev_min)+($lev_del)`,`($lev_max)-($lev_del)`,($lev_del))(inf) CI[d=7,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Era Interim $2 climatology (1979-2013)"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Sea-ice cover (0-1)"
PPL SHASET RESET

set viewport LL
shade/nolabels/levels=($lev_text) ICE_MASK[d=3,l=$1]; go land
!---xlabel-------------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.5*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Longitude"
!---ylabel-------------------------------------------------
let xpos=`-0.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Latitude"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "CM2M $2 mean"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Sea-ice cover (0-1)"
PPL SHASET RESET

set viewport LR
plot/nolabels CI[d=7,l=$1,i=@ave], ICE_MASK[d=8,l=$1,i=@ave],ICE_MASK[d=3,i=@ave,l=$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.5*($PPL$XORG)`
let ypos=`0.5*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,0,90,.15 "Sea-ice cover"
!---titel 1 - above--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`($PPL$YLEN)+0.2*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "$2 Zonal Means"
!---titel 2 - below--------------------------------------------
let xpos=`0.5*($PPL$XLEN)`
let ypos=`-0.8*($PPL$YORG)`
label/nouser `xpos`,`ypos`,0,0,.15 "Sea-ice cover (0-1)"
!---key--------------------------------------------------
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.9*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P1Era Interim"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.85*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P2Aeolus"
let xpos=`0.05*($PPL$XLEN)`
let ypos=`0.8*($PPL$YLEN)`
label/nouser `xpos`,`ypos`,-1,0,.1 "@P3CM2M"
PPL SHASET RESET

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


!!

!!!!! The end !!!!!
if `"($PROGRAM_NAME)" eq "Ferret"` then
  CANCEL MODE METAFILE
  spawn Fprint -l cps -o albedo_1.ps -R albedo.plt
  spawn Fprint -l cps -o albedo.ps -R albedo.plt.~*~
  spawn cat albedo_1.ps >> albedo.ps
  spawn rm albedo.plt*
else
  spawn cat albedo_?.ps > albedo.ps
endif
spawn ps2pdf albedo.ps albedo_seas_$2.pdf
spawn rm albedo_?.ps albedo.ps 
spawn rm ferret.jnl
cancel data/all
!spawn rm Albedo-data_sm_taxis.nc Albedo-data_down_sm_taxis.nc Albedo-data_flux_sm_taxis.nc