TENSTREAM=$HOME/tenstream/build

THETA=40
Nz=60
LOGDIR=logs_Nz$Nz/
mkdir -p $LOGDIR

cat $HOME/libRadtran/data/atmmod/afglus.dat | tail -n 21 > afglus.dat
OPT="-atm_filename afglus.dat -skip_load_LUT -theta0 $THETA -Nz $Nz "

make -C $TENSTREAM ex_pprts_rrtm_iterations -j || exit
[ ! -e $LOGDIR/log.solar.twostr      ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -twostr_only -thermal no                                     |& tee $LOGDIR/log.solar.twostr
[ ! -e $LOGDIR/log.solar.rrtmg       ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -rrtmg_only  -thermal no                                     |& tee $LOGDIR/log.solar.rrtmg
[ ! -e $LOGDIR/log.solar.disort      ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -disort_only -thermal no -disort_streams 32                  |& tee $LOGDIR/log.solar.disort
[ ! -e $LOGDIR/log.solar.disort6     ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -disort_only -thermal no -disort_streams 6                   |& tee $LOGDIR/log.solar.disort6

[ ! -e $LOGDIR/log.thermal.twostr    ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -twostr_only -solar no                                       |& tee $LOGDIR/log.thermal.twostr
[ ! -e $LOGDIR/log.thermal.schwarz5  ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -twostr_only -solar no -schwarzschild -schwarzschild_Nmu 5   |& tee $LOGDIR/log.thermal.schwarz5
[ ! -e $LOGDIR/log.thermal.schwarz50 ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -twostr_only -solar no -schwarzschild -schwarzschild_Nmu 50  |& tee $LOGDIR/log.thermal.schwarz50
[ ! -e $LOGDIR/log.thermal.rrtmg     ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -rrtmg_only  -solar no                                       |& tee $LOGDIR/log.thermal.rrtmg
[ ! -e $LOGDIR/log.thermal.disort    ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -disort_only -solar no -disort_streams 32                    |& tee $LOGDIR/log.thermal.disort
[ ! -e $LOGDIR/log.thermal.disort6   ] && $TENSTREAM/bin/ex_pprts_rrtm_iterations $OPT -disort_only -solar no -disort_streams 6                     |& tee $LOGDIR/log.thermal.disort6

rm $LOGDIR/*{atm,clean}*
for f in $LOGDIR/log.*.*
do
  grep 'edn' $f > ${f}.clean
  grep "zt.*plev.*Tlev" $f > ${f}.atm${Nz}
done

cat > plot.py << EOF
#!/bin/env python3
from pylab import *
import glob
figure(figsize=(12,12))

benchmark = glob.glob('${LOGDIR}/log.solar.disort.clean')[0]
k, edir0, edn0, eup0, abso0 = np.genfromtxt(benchmark, usecols=(0,2,4,6,7), unpack=True )

atmfile = glob.glob('${LOGDIR}/log.*.atm${Nz}')[0]
zt, plev, Tlev = np.genfromtxt(atmfile, usecols=(2,4,6), unpack=True )

subplot(421); cla(); xlabel('$\ E_{dn}^{solar}$')
subplot(422); cla(); xlabel('$\ E_{up}^{solar}$')
subplot(423); cla(); xlabel('$\Delta\ E_{dn}^{solar}$')
subplot(424); cla(); xlabel('$\Delta\ E_{up}^{solar}$')

for f in sorted(glob.glob('${LOGDIR}/log.solar.*.clean')):
  ident = f.split('.')[2]
  k, edir, edn, eup, abso = np.genfromtxt(f, usecols=(0,2,4,6,7), unpack=True )
  subplot(421)
  plot(edir+edn, plev, label=ident)
  subplot(422)
  plot(eup, plev, label=ident)
  subplot(423)
  plot(edir+edn-edir0-edn0, plev, label=ident)
  subplot(424)
  plot(eup-eup0, plev, label=ident)

benchmark = glob.glob('${LOGDIR}/log.thermal.disort.clean')[0]
k, edn0, eup0, abso0 = np.genfromtxt(benchmark, usecols=(0,2,4,5), unpack=True )

subplot(425); cla(); xlabel('$\ E_{dn}^{thermal}$')
subplot(426); cla(); xlabel('$\ E_{up}^{thermal}$')
subplot(427); cla(); xlabel('$\Delta\ E_{dn}^{thermal}$')
subplot(428); cla(); xlabel('$\Delta\ E_{up}^{thermal}$')

for f in sorted(glob.glob('${LOGDIR}/log.thermal.*.clean')):
  ident = f.split('.')[2]
  k, edn, eup, abso = np.genfromtxt(f, usecols=(0,2,4,5), unpack=True )
  subplot(425)
  plot(edn, plev, label=ident)
  subplot(426)
  plot(eup, plev, label=ident)
  subplot(427)
  plot(edn-edn0, plev, label=ident)
  subplot(428)
  plot(eup-eup0, plev, label=ident)

[ ax.legend(loc='best', framealpha=.5, fontsize='small') for ax in gcf().get_axes() ]
[ ax.set_ylabel('pressure [hPa]') for ax in gcf().get_axes() ]
[ ax.set_ylim(reversed(ax.get_ylim())) for ax in gcf().get_axes() ]
[ ax.axhline(700, ls='--', color='.5') for ax in gcf().get_axes() ]
[ ax.axhline(800, ls='--', color='.5') for ax in gcf().get_axes() ]
tight_layout()
savefig('diff_plot_Nz_${Nz}.pdf', bbox_inches='tight')
EOF

ipython3 plot.py

[ -e article.md ] && touch article.md
