#--this is a script that does it all

dirout="OUTPUT"
if [ ! -d "$dirout" ]; then
mkdir $dirout
fi

#--first clean up things
rm -f solarforcing_picontrol_fx_3.1_SSI_regridded_2_TSI_anomaly.nc
rm -f solarforcing_picontrol_fx_3.1_SSI_regridded_6_TSI_anomaly.nc
rm -f solarforcing_ref_day_3.1_SSI_regridded_2_TSI_anomaly.nc
rm -f solarforcing_ref_day_3.1_SSI_regridded_6_TSI_anomaly.nc
rm -f solarforcing_ext_day_3.1_SSI_regridded_2_TSI_anomaly.nc
rm -f solarforcing_ext_day_3.1_SSI_regridded_6_TSI_anomaly.nc

#--run matlab script for piControl
matlab -r "run  CMIP6_piControl_2bands.m, quit"
matlab -r "run  CMIP6_piControl_6bands.m, quit"

#--run matlab script for Historical and Future
matlab -r "run  CMIP6_historical_2bands.m, quit"
matlab -r "run  CMIP6_historical_6bands.m, quit"

#--run matlab script for Future Sensitivity
matlab -r "run  CMIP6_sensitivity_2bands.m, quit"
matlab -r "run  CMIP6_sensitivity_6bands.m, quit"

output=${dirout}"/Control"
if [ ! -d "$output" ]; then
mkdir $output
fi

#--as computed by CMIP6 for piControl
filein2bands=./solarforcing_picontrol_fx_3.1_SSI_regridded_2_TSI_anomaly.nc
filein6bands=./solarforcing_picontrol_fx_3.1_SSI_regridded_6_TSI_anomaly.nc

cdo duplicate,360 $filein2bands ${output}/solarforcing_v3.1_piControl_ave_2bands_360days.nc
cdo duplicate,360 $filein6bands ${output}/solarforcing_v3.1_piControl_ave_6bands_360days.nc
cdo duplicate,365 $filein2bands ${output}/solarforcing_v3.1_piControl_ave_2bands_365days.nc
cdo duplicate,365 $filein6bands ${output}/solarforcing_v3.1_piControl_ave_6bands_365days.nc
cdo duplicate,366 $filein2bands ${output}/solarforcing_v3.1_piControl_ave_2bands_366days.nc
cdo duplicate,366 $filein6bands ${output}/solarforcing_v3.1_piControl_ave_6bands_366days.nc

#--now the same duplication but to match Gregorian calendar
outputG=${dirout}"/ControlGregorian"
if [ ! -d "$outputG" ]; then
mkdir $outputG
fi

for year in {1850..2500}
do 
nbdays=`python -c "import calendar ; print 366 if calendar.isleap(${year}) else 365"`
rm -f ${outputG}/solarforcing_v3.1_piControl_ave_2bands_${year}.nc
rm -f ${outputG}/solarforcing_v3.1_piControl_ave_6bands_${year}.nc
ln -s ../Control/solarforcing_v3.1_piControl_ave_2bands_${nbdays}days.nc ${outputG}/solarforcing_v3.1_piControl_ave_2bands_${year}.nc
ln -s ../Control/solarforcing_v3.1_piControl_ave_6bands_${nbdays}days.nc ${outputG}/solarforcing_v3.1_piControl_ave_6bands_${year}.nc
done

rm -f ${filein2bands} ${filein6bands}

#--present-day average for unofficial CMIP6 pdControl_ave

filein2bands=./solarforcing_ref_day_3.1_SSI_regridded_2_TSI_anomaly.nc
filein6bands=./solarforcing_ref_day_3.1_SSI_regridded_6_TSI_anomaly.nc

#year1=1995
#year2=2004
year1=1991
year2=2010
start_time=`python -c "from datetime import date; print (date($year1,1,1)-date(1850,1,1)).days+1"`
end_time=`python -c "from datetime import date; print (date($year2,12,31)-date(1850,1,1)).days"`
echo $year1 $year2 $start_time $end_time

fileout2bands=solarforcing_ref_day_3.1_SSI_regridded_2_TSI_anomaly_${year1}-${year2}.nc
fileout6bands=solarforcing_ref_day_3.1_SSI_regridded_6_TSI_anomaly_${year1}-${year2}.nc

#--somewhat dirty average as in principle ssi_frac should be weighted by tsi !
#--note that seltimestep starts with 1
cdo timmean -seltimestep,$start_time,$end_time $filein2bands $fileout2bands
cdo timmean -seltimestep,$start_time,$end_time $filein6bands $fileout6bands

#--duplicate 360 365 or 366 values for annual file 
cdo duplicate,360 $fileout2bands ${output}/solarforcing_v3.1_pdControl_ave_2bands_360days.nc
cdo duplicate,360 $fileout6bands ${output}/solarforcing_v3.1_pdControl_ave_6bands_360days.nc
cdo duplicate,365 $fileout2bands ${output}/solarforcing_v3.1_pdControl_ave_2bands_365days.nc
cdo duplicate,365 $fileout6bands ${output}/solarforcing_v3.1_pdControl_ave_6bands_365days.nc
cdo duplicate,366 $fileout2bands ${output}/solarforcing_v3.1_pdControl_ave_2bands_366days.nc
cdo duplicate,366 $fileout6bands ${output}/solarforcing_v3.1_pdControl_ave_6bands_366days.nc

#--now the same duplication but to match Gregorian calendar
for year in {1850..2500}
do 
nbdays=`python -c "import calendar ; print 366 if calendar.isleap(${year}) else 365"`
rm -f ${outputG}/solarforcing_v3.1_pdControl_ave_2bands_${year}.nc
rm -f ${outputG}/solarforcing_v3.1_pdControl_ave_6bands_${year}.nc
ln -s ../Control/solarforcing_v3.1_pdControl_ave_2bands_${nbdays}days.nc ${outputG}/solarforcing_v3.1_pdControl_ave_2bands_${year}.nc
ln -s ../Control/solarforcing_v3.1_pdControl_ave_6bands_${nbdays}days.nc ${outputG}/solarforcing_v3.1_pdControl_ave_6bands_${year}.nc
done

rm -f $fileout2bands $fileout6bands

#--split files for historical
filein2bands=./solarforcing_ref_day_3.1_SSI_regridded_2_TSI_anomaly.nc
filein6bands=./solarforcing_ref_day_3.1_SSI_regridded_6_TSI_anomaly.nc

output=${dirout}"/Historical"
if [ ! -d "$output" ]; then
mkdir $output
fi

#--note that the time axis starts from 0 with ncks.... so I put -1 here
end_time=-1

for year in {1850..2014}
do 
yearp1=$((year+1))
nbday=`python -c "from datetime import date; print (date($yearp1,1,1)-date($year,1,1)).days"`
start_time=$((end_time+1))
end_time=$((start_time+nbday-1))
echo $year $yearp1 $start_time $end_time
rm -f ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
rm -f ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
ncks -d time,$start_time,$end_time $filein2bands ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
ncks -d time,$start_time,$end_time $filein6bands ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
done 

#--split files for Future
output=${dirout}"/Future"
if [ ! -d "$output" ]; then
mkdir $output
fi

for year in {2015..2299}
do 
yearp1=$((year+1))
nbday=`python -c "from datetime import date; print (date($yearp1,1,1)-date($year,1,1)).days"`
start_time=$((end_time+1))
end_time=$((start_time+nbday-1))
echo $year $yearp1 $start_time $end_time
rm -f ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
rm -f ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
ncks -d time,$start_time,$end_time $filein2bands ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
ncks -d time,$start_time,$end_time $filein6bands ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
done 

rm -f $filein2bands $filein6bands

#--split files for Future Sensitivity
filein2bands=./solarforcing_ext_day_3.1_SSI_regridded_2_TSI_anomaly.nc
filein6bands=./solarforcing_ext_day_3.1_SSI_regridded_6_TSI_anomaly.nc

output=${dirout}"/FutureSensitivity"
if [ ! -d "$output" ]; then
mkdir $output
fi

#--reinitialising end_time to 31/12/2014
end_time=`python -c "from datetime import date; print (date(2014,12,31)-date(1850,1,1)).days"`

for year in {2015..2299}
do 
yearp1=$((year+1))
nbday=`python -c "from datetime import date; print (date($yearp1,1,1)-date($year,1,1)).days"`
start_time=$((end_time+1))
end_time=$((start_time+nbday-1))
echo $year $yearp1 $start_time $end_time
rm -f ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
rm -f ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
ncks -d time,$start_time,$end_time $filein2bands ${output}/solarforcing_v3.1_daily_2bands_${year}.nc
ncks -d time,$start_time,$end_time $filein6bands ${output}/solarforcing_v3.1_daily_6bands_${year}.nc
done 

rm -f $filein2bands $filein6bands