# Script générant les fichiers ozone hybride CNRM pour les ssp du Tier 2 (en l'absence de l'ozone officiel UReading)
# -----------------------------------
# Th. Lurton - janvier 2019
# -----------------------------------
# Ce script est à exécuter sur Ciclad
# inclus dans un job.
# -----------------------------------

# Chemins à personnaliser
ici=/data/tlurton/CMIP6/O3_Tier2/
repout=${ici}Tier2_OUT/

# Chemin des données fournies par le CNRM
repCNRM=/data/tlurton/CMIP6/OzoneCNRM/

# Grande boucle sur les années
# On a 2100 du CNRM, mais pas 2100 d'UReading, ce qui nous oblige à ne pas traiter la dernière année. Au final on dupliquera 2100 = 2099.
for year in {2015..2099}
do

echo "========================================================================="
echo "================================= "${year}" ================================="
echo "========================================================================="

# PHASE 1
# INTERPOLATION VERTICALE NIVEAUX CNRM --> NIVEAUX UREADING

# --> ici boucle sur tous les ssp
# On n'aura pas besoin de ssp585

for ssp in "ssp119" "ssp126" "ssp245" "ssp370" "ssp434" "ssp460" "ssp534-over"
do

echo "Processing year "${year}" for "${ssp}"..."

fileCNRM=${repCNRM}o3_HOMEPmon_CNRM-ESM2-1_${ssp}_r1i1p1f2_gr_${year}01-${year}12.nc

filetmp1=${ici}${ssp}_${year}_CNRM_hPa.nc
filetmp2=${ici}${ssp}_${year}_CNRM_hPa_interpURea.nc

# Attention !
# Les niveaux plev des fichiers CNRM sont en Pa
# Les niveaux plev des fichiers UReading sont en hPa

# Niveaux CNRM entre 1000 hPa 3 Pa :
# [Pa] 100000,92500,85000,70000,60000,50000,40000,30000,25000,20000,17000,15000,13000,11500,10000,9000,8000,7000,5000,3000,2000,1500,1000,700,500,300,200,150,100,70,50,40,30,20,15,10,7,5,3
# [hPa] 1000,925,850,700,600,500,400,300,250,200,170,150,130,115,100,90,80,70,50,30,20,15,10,7,5,3,2,1.5,1,.7,.5,.4,.3,.2,.15,.1,.07,.05,.03

# Niveaux Reading, liste totale
# [hPa] 1000,925,850,800,780,750,700,650,600,500,450,400,350,300,285,250,200,170,150,130,115,100,90,80,70,60,50,40,35,30,25,20,15,10,7,5,4,3,2,1.5,1,0.7,0.5,0.4,0.3,0.2,0.15,0.1,0.07,0.05,0.04,0.03,0.02,0.015,0.01,0.007,0.005,0.004,0.003,0.002,0.0015,0.001,0.0008,0.0005,0.0003,0.0001

# Niveaux Reading entre 1000 hPa et 3 Pa
# [hPa] 1000,925,850,800,780,750,700,650,600,500,450,400,350,300,285,250,200,170,150,130,115,100,90,80,70,60,50,40,35,30,25,20,15,10,7,5,4,3,2,1.5,1,0.7,0.5,0.4,0.3,0.2,0.15,0.1,0.07,0.05,0.04,0.03

# Conversion des plev CNRM vers hPa, en utilisant le fichier "CNRM_axis_en_hPa" : overwrite de la dimension verticale
echo "... Converting CNRM pressure levels towards hPa..."
cdo setzaxis,CNRM_axis_en_hPa ${fileCNRM} ${filetmp1}

# Projection sur les niveaux Reading
echo "... Interpolating CNRM file onto UReading pressure levels..."
cdo intlevel,1000,925,850,800,780,750,700,650,600,500,450,400,350,300,285,250,200,170,150,130,115,100,90,80,70,60,50,40,35,30,25,20,15,10,7,5,4,3,2,1.5,1,0.7,0.5,0.4,0.3,0.2,0.15,0.1,0.07,0.05,0.04,0.03 ${filetmp1} ${filetmp2}

rm -f ${filetmp1}

# --> Fin des deux boucles sur les ssp
done


# PHASE 2
# "RÈGLE DE TROIS"

# --> Boucle sur les ssp du Tier 2 uniquement : ssp119, ssp434, ssp460 et ssp534-over

for sspT2 in "ssp119" "ssp434" "ssp460" "ssp534-over"
do

# choix des ssp-Tier1 associés
if [ ${sspT2} = "ssp119" ]
then
sspT1a=ssp126
sspT1b=ssp126
fi
if [ ${sspT2} = "ssp434" ]
then
sspT1a=ssp126
sspT1b=ssp245
fi
if [ ${sspT2} = "ssp460" ]
then
sspT1a=ssp245
sspT1b=ssp370
fi
if [ ${sspT2} = "ssp534-over" ]
then
sspT1a=ssp126
sspT1b=ssp245
fi

fileCNRMT1a=${ici}${sspT1a}_${year}_CNRM_hPa_interpURea.nc
fileCNRMT1b=${ici}${sspT1b}_${year}_CNRM_hPa_interpURea.nc
fileCNRMT2=${ici}${sspT2}_${year}_CNRM_hPa_interpURea.nc

if [ ${year} -le 2049 ]
then
fileURead1aori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1a}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1a}-1-0_gn_201501-204912.nc
fileURead1bori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1b}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1b}-1-0_gn_201501-204912.nc
else
fileURead1aori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1a}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1a}-1-0_gn_205001-209912.nc
fileURead1bori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1b}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1b}-1-0_gn_205001-209912.nc
fi

fileUReadayear=${ici}o3_URead_${sspT1a}_${year}.nc
fileUReadbyear=${ici}o3_URead_${sspT1b}_${year}.nc

fileUReadatronq=${ici}o3_URead_${sspT1a}_${year}_tronq.nc
fileUReadbtronq=${ici}o3_URead_${sspT1b}_${year}_tronq.nc

# Fichier contenant les latitudes UReading exactes
latURead=${ici}lat_URead.nc

filevsa=${ici}o3_${sspT2}_${year}_vs_${sspT1a}.nc
filevsb=${ici}o3_${sspT2}_${year}_vs_${sspT1b}.nc

echo "... Selecting year "${year}" in Tier-1 UReading ssp "${sspT1a}"..."
# Extraction de l'année d'intérêt du ssp-Tier1 de UReading
cdo selyear,${year} ${fileURead1aori} ${fileUReadayear}
# Pas besoin de la seconde borne Tier1 si on est dans le cas ssp119 :
if [ ${sspT1a} != ${sspT1b} ]
then
echo "... Selecting year "${year}" in Tier-1 UReading ssp "${sspT1b}"..."
cdo selyear,${year} ${fileURead1bori} ${fileUReadbyear}
fi

rm -f ${latURead}
ncks -s 'lat=double(lat)' ${fileUReadayear} ${latURead}

echo "... Chopping up pressure levels, retaining only > 3 Pa..."
# Troncature des niveaux verticaux
cdo sellevel,1000,925,850,800,780,750,700,650,600,500,450,400,350,300,285,250,200,170,150,130,115,100,90,80,70,60,50,40,35,30,25,20,15,10,7,5,4,3,2,1.5,1,0.7,0.5,0.4,0.3,0.2,0.15,0.1,0.07,0.05,0.04,0.03 ${fileUReadayear} ${fileUReadatronq}
# On vient d'enlever les 14 niveaux supérieurs du fichier UReading. 66-14 = 52 niveaux restent.
# On fait la même chose pour la seconde borne Tier1, sauf dans le cas ssp119, où on n'en a pas besoin.
if [ ${sspT1a} != ${sspT1b} ]
then
cdo sellevel,1000,925,850,800,780,750,700,650,600,500,450,400,350,300,285,250,200,170,150,130,115,100,90,80,70,60,50,40,35,30,25,20,15,10,7,5,4,3,2,1.5,1,0.7,0.5,0.4,0.3,0.2,0.15,0.1,0.07,0.05,0.04,0.03 ${fileUReadbyear} ${fileUReadbtronq}
fi

# Remappage horizontal des fichiers CNRM
fileCNRMT1arem=${ici}o3_CNRM_${sspT1a}_${year}_remap.nc
fileCNRMT1brem=${ici}o3_CNRM_${sspT1b}_${year}_remap.nc
fileCNRMT2rem=${ici}o3_CNRM_${sspT2}_${year}_remap.nc

echo "... Remapping CNRM file for "${sspT2}" and year "${year}"..."
cdo remapcon,grid_URead_144x96 ${fileCNRMT2} ${fileCNRMT2rem}
echo "... Remapping CNRM file for "${sspT1a}" and year "${year}"..."
cdo remapcon,grid_URead_144x96 ${fileCNRMT1a} ${fileCNRMT1arem}
rm -f ${fileCNRMT1a}
if [ ${sspT1a} != ${sspT1b} ]
then
echo "... Remapping CNRM file for "${sspT1b}" and year "${year}"..."
cdo remapcon,grid_URead_144x96 ${fileCNRMT1b} ${fileCNRMT1brem}
fi

# Fix des latitudes, qui ne tombent pas juste
echo "... Fixing latitudes..."
ncks -3 -A -v lat ${latURead} ${fileCNRMT2rem}
ncks -3 -A -v lat ${latURead} ${fileCNRMT1arem}
if [ ${sspT1a} != ${sspT1b} ]
then
ncks -3 -A -v lat ${latURead} ${fileCNRMT1brem}
fi

# La règle de trois en soi : l'ozone s'appelle "o3" dans les fichiers CNRM, et "vmro3" dans les fichiers UReading.
echo "... Performing scaling : "${sspT2}" created by scaling against "${sspT1a}", for year "${year}"..."
cdo -div -mul -selname,o3 ${fileCNRMT2rem} -selname,vmro3 ${fileUReadatronq} -selname,o3 ${fileCNRMT1arem} ${filevsa}
if [ ${sspT1a} != ${sspT1b} ]
then
echo "... Performing scaling : "${sspT2}" created by scaling against "${sspT1b}", for year "${year}"..."
cdo -div -mul -selname,o3 ${fileCNRMT2rem} -selname,vmro3 ${fileUReadbtronq} -selname,o3 ${fileCNRMT1brem} ${filevsb}
fi

# Fin des boucles ssp-Tier2
done


# PHASE 3
# "CHAPEAUTAGE"

if [ ${year} -le 2049 ]
then
fileURead1ori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1a}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1a}-1-0_gn_201501-204912.nc
fileURead2ori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1b}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1b}-1-0_gn_201501-204912.nc
else
fileURead1ori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1a}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1a}-1-0_gn_205001-209912.nc
fileURead2ori=/prodigfs/project/input4MIPs/CMIP6/ScenarioMIP/UReading/UReading-CCMI-${sspT1b}-1-0/atmos/mon/vmro3/gn/v20181101/vmro3_input4MIPs_ozone_ScenarioMIP_UReading-CCMI-${sspT1b}-1-0_gn_205001-209912.nc
fi

fileURead1year=${ici}o3_URead_${sspT1a}_${year}.nc
fileURead2year=${ici}o3_URead_${sspT1b}_${year}.nc

echo "... Preparing hat for year "${year}"..."
if [ ${sspT1a} != ${sspT1b} ]
then
cdo selyear,${year} ${fileURead1ori} ${fileURead1year}
cdo selyear,${year} ${fileURead2ori} ${fileURead2year}
else
cdo selyear,${year} ${fileURead1ori} ${fileURead1year}
fi

fileURead1chapo=${ici}o3_URead_${sspT1a}_${year}_chapo.nc
fileURead2chapo=${ici}o3_URead_${sspT1b}_${year}_chapo.nc

if [ ${sspT1a} != ${sspT1b} ]
then
cdo sellevel,0.02,0.015,0.01,0.007,0.005,0.004,0.003,0.002,0.0015,0.001,0.0008,0.0005,0.0003,0.0001 ${fileURead1year} ${fileURead1chapo}
cdo sellevel,0.02,0.015,0.01,0.007,0.005,0.004,0.003,0.002,0.0015,0.001,0.0008,0.0005,0.0003,0.0001 ${fileURead2year} ${fileURead2chapo}
else
cdo sellevel,0.02,0.015,0.01,0.007,0.005,0.004,0.003,0.002,0.0015,0.001,0.0008,0.0005,0.0003,0.0001 ${fileURead1year} ${fileURead1chapo}
fi

fileUReadchapo=${ici}o3_URead_${year}_chapo.nc

# Moyenne des chapeaux des deux ssp
if [ ${sspT1a} != ${sspT1b} ]
then
echo "... Averaging hats for Tier 1 UReading "${sspT1a}" and "${sspT1b}", for "${year}"..."
cdo mulc,.5 -add -selname,vmro3 ${fileURead1chapo} -selname,vmro3 ${fileURead2chapo} ${fileUReadchapo}
else
echo "... Using solely hat from Tier 1 UReading "${sspT1a}", for year "${year}"..."
cp ${fileURead1chapo} ${fileUReadchapo}
fi

echo "... Cleaning up..."
if [ ${sspT1a} != ${sspT1b} ]
then
rm -f ${fileURead1year} ${fileURead2year} ${fileURead1chapo} ${fileURead2chapo}
else
rm -f ${fileURead1year} ${fileURead1chapo}
fi

# Ajout du chapeau aux ssp Tier2 créés à l'étape 2

for sspT2 in "ssp119" "ssp434" "ssp460" "ssp534-over"
do

# choix des ssp-Tier1 associés
if [ ${sspT2} = "ssp119" ]
then
sspT1a=ssp126
sspT1b=ssp126
fi
if [ ${sspT2} = "ssp434" ]
then
sspT1a=ssp126
sspT1b=ssp245
fi
if [ ${sspT2} = "ssp460" ]
then
sspT1a=ssp245
sspT1b=ssp370
fi
if [ ${sspT2} = "ssp534-over" ]
then
sspT1a=ssp126
sspT1b=ssp245
fi

filevsa=${ici}o3_${sspT2}_${year}_vs_${sspT1a}.nc
filevsb=${ici}o3_${sspT2}_${year}_vs_${sspT1b}.nc

fileUReadfinala=${repout}${sspT2}/o3_${sspT2}_${year}_vs_${sspT1a}.nc
fileUReadfinalb=${repout}${sspT2}/o3_${sspT2}_${year}_vs_${sspT1b}.nc

echo "... Renaming O3 variable and fixing NetCDF-3 inconsistency..."
ncrename -v o3,vmro3 ${filevsa}
ncks -3 ${filevsa} ${ici}tmp.nc
mv ${ici}tmp.nc ${filevsa}
if [ ${sspT1a} != ${sspT1b} ]
then
ncrename -v o3,vmro3 ${filevsb}
ncks -3 ${filevsb} ${ici}tmp.nc
mv ${ici}tmp.nc ${filevsb}
fi

echo "... Topping newly created "${sspT2}" (against "${sspT1a}") for year "${year}" with hat..."
cdo merge ${filevsa} ${fileUReadchapo} ${fileUReadfinala}
if [ ${sspT1a} != ${sspT1b} ]
then
echo "... Topping newly created "${sspT2}" (against "${sspT1b}") for year "${year}" with hat..."
cdo merge ${filevsb} ${fileUReadchapo} ${fileUReadfinalb}
fi

# Je moyenne les deux versions
echo "... Merging "${sspT2}" against "${sspT1a}" and "${sspT2}" against "${sspT1b}" into a single file..."
filemerge=${repout}${sspT2}/o3_${sspT2}_${year}.nc
if [ ${sspT1a} != ${sspT1b} ]
then
cdo add -mulc,.5 -selname,vmro3 ${fileUReadfinala} -mulc,.5 -selname,vmro3 ${fileUReadfinalb} ${filemerge}
else
cp ${fileUReadfinala} ${filemerge}
fi

# Petite entourloupe finale pour avoir précisément les niveaux de pression d'origine (entre-temps, on est passés en double avec petite erreur d'arrondi...)
echo "... Fixing z-axis for year "${year}"..."
cdo setzaxis,URead_axis_en_hPa ${filemerge} ${repout}tmp.nc
mv ${repout}tmp.nc ${filemerge}

# Fix des latitudes
ncks -3 -A -v lat ${latURead} ${filemerge}

# Fin de la boucle sur les sspT2
done

rm -f ${fileUReadchapo}

# Ménage
rm -f ${ici}o3*${year}*.nc ${ici}ssp*${year}*.nc

# Fin de la grande boucle sur les années
done