#--Script performing bias-correction of RCP scenarios with the AMIP runs
#--over SST and SIC.
#--Th. Lurton, IPSL, 07/2018

#--working directory
ici=${GENCMIP6_CCCWORKDIR}/CMIP6_LIMIT/
if [ ! -d ${ici} ] ; then mkdir -p ${ici} ; fi

gridfile=/ccc/work/cont003/gencmip6/lurtont/CMIP6_LIMIT/GRID/grid_MR_CMIP5

#--AMIP input files
dirinAMIP=/ccc/work/cont003/gencmip6/lurtont/CMIP6_LIMIT/AMIP/

#--performing horizontal remapping for AMIP files between 2006 and 2017
fileintsoloceAMIP=${dirinAMIP}tos_input4MIPs_SSTsAndSeaIce_CMIP_PCMDI-AMIP-1-1-4_gn_187001-201712.nc
fileinpourcsicAMIP=${dirinAMIP}siconc_input4MIPs_SSTsAndSeaIce_CMIP_PCMDI-AMIP-1-1-4_gn_187001-201712.nc

fileinSSTyearsAMIP=${ici}SST_AMIP_2006_2017.nc
fileinSICyearsAMIP=${ici}SIC_AMIP_2006_2017.nc

cdo -seldate,2006-01-01,2017-12-31 ${fileintsoloceAMIP} ${fileinSSTyearsAMIP}
cdo -seldate,2006-01-01,2017-12-31 ${fileinpourcsicAMIP} ${fileinSICyearsAMIP}

#--AMIP files remapped to CMIP5-MR grid
fileremaptsoloceAMIP=${ici}amipbc_sst_remapped_2006_2017.nc
fileremappourcsicAMIP=${ici}amipbc_sic_remapped_2006_2017.nc

echo "Remapping AMIP SST file towards CMIP5 MR grid, for years 2006-2017..."
cdo remapcon,${gridfile} ${fileinSSTyearsAMIP} ${fileremaptsoloceAMIP}
echo "Remapping AMIP SIC file towards CMIP5 MR grid, for years 2006-2017..."
cdo remapcon,${gridfile} ${fileinSICyearsAMIP} ${fileremappourcsicAMIP}


#--loop on RCPs
for rcp in rcp26 rcp45 rcp60 rcp85
do

if [ ${rcp} == "rcp26" ] ; then dirinRCP=/ccc/store/cont003/dsmipsl/p86maf/dmf_import/IGCM_OUT/IPSLCM5A-MR/PROD/rcp26/v5.rcp26MR1/ATM/Analyse/TS_MO/ ; fi
if [ ${rcp} == "rcp45" ] ; then dirinRCP=/ccc/store/cont003/dsmipsl/p86maf/dmf_import/IGCM_OUT/IPSLCM5A-MR/PROD/rcp45/v5.rcp45MR1/ATM/Analyse/TS_MO/ ; fi
if [ ${rcp} == "rcp60" ] ; then dirinRCP=/ccc/store/cont003/dsm/p86denv/IGCM_OUT/IPSLCM5A-MR/PROD/rcp60/v5.rcp60MR1/ATM/Analyse/TS_MO/ ; fi
if [ ${rcp} == "rcp85" ] ; then dirinRCP=/ccc/store/cont003/dsmipsl/p86maf/dmf_import/IGCM_OUT/IPSLCM5A-MR/PROD/rcp85/v5.rcp85MR1/ATM/Analyse/TS_MO/ ; fi

rcpdir=${ici}Scenarios/${rcp}/
if [ ! -d ${rcpdir} ] ; then mkdir -p ${rcpdir} ; fi

rcpsplitdir=${rcpdir}split/
if [ ! -d ${rcpsplitdir} ] ; then mkdir -p ${rcpsplitdir} ; fi

#--RCP input files
fileintsoloceRCP=${dirinRCP}v5.${rcp}MR1_20060101_21001231_1M_tsol_oce.nc
fileinpourcsicRCP=${dirinRCP}v5.${rcp}MR1_20060101_21001231_1M_pourc_sic.nc
fileinpourclicRCP=${dirinRCP}v5.${rcp}MR1_20060101_21001231_1M_pourc_lic.nc

#--splitting years
echo ${rcp}": splitting years for SST..."
cdo splityear ${fileintsoloceRCP} ${rcpsplitdir}v5.${rcp}MR1_tsol_oce_
echo ${rcp}": splitting years for SIC..."
cdo splityear ${fileinpourcsicRCP} ${rcpsplitdir}v5.${rcp}MR1_pourc_sic_
echo ${rcp}": splitting years for LIC..."
cdo splityear ${fileinpourclicRCP} ${rcpsplitdir}v5.${rcp}MR1_pourc_lic_

dirout=${rcpdir}/OUT_anocorr/
if [ ! -d ${dirout} ] ; then mkdir -p ${dirout} ; fi

tmpavgsstAMIP=${ici}tmpavgsstAMIP.nc
tmpavgsicAMIP=${ici}tmpavgsicAMIP.nc
tmpavgsstRCP=${ici}tmpavgsstRCP.nc
tmpavgsicRCP=${ici}tmpavgsicRCP.nc
tmpavglicRCP=${ici}tmpavglicRCP.nc
cdo ymonmean -seldate,2006-01-01,2017-12-31 ${fileremaptsoloceAMIP} ${tmpavgsstAMIP}
cdo ymonmean -seldate,2006-01-01,2017-12-31 ${fileremappourcsicAMIP} ${tmpavgsicAMIP}
cdo ymonmean -seldate,2006-01-01,2017-12-31 ${fileintsoloceRCP} ${ici}sstRCP.nc
cdo ymonmean -seldate,2006-01-01,2017-12-31 ${fileinpourcsicRCP} ${ici}sicRCP.nc
cdo ymonmean -seldate,2006-01-01,2017-12-31 ${fileinpourclicRCP} ${tmpavglicRCP}
ncrename -v siconc,pourc_sic ${tmpavgsicAMIP}
ncrename -v tos,tsol_oce ${tmpavgsstAMIP} -O ${ici}tmp.nc
cdo expr,'tsol_oce=tsol_oce+273.15' ${ici}tmp.nc ${tmpavgsstAMIP}
rm -f ${ici}tmp.nc
ncrename -v time,time_counter -d time,time_counter ${tmpavgsstAMIP}
ncrename -v time,time_counter -d time,time_counter ${tmpavgsicAMIP}
ncrename -v pourc_lic,pourc_sic ${tmpavglicRCP} -O ${ici}tmp.nc
ncks -C -O -x -v aire ${ici}tmp.nc ${tmpavglicRCP}
rm -f ${ici}tmp.nc
ncks -C -O -x -v aire ${ici}sstRCP.nc ${tmpavgsstRCP}
ncks -C -O -x -v aire ${ici}sicRCP.nc ${tmpavgsicRCP}
rm -f ${ici}sstRCP.nc ${ici}sicRCP.nc

for year in {2015..2100}
do

echo ${rcp}": Correcting anomaly for year "${year}"..."

fileyeartsoloceRCP=${rcpsplitdir}v5.${rcp}MR1_tsol_oce_${year}.nc
fileouttsoloce=${dirout}v5.${rcp}MR1_tsol_oce_${year}_anocorr.nc

fileyearpourcsicRCP=${rcpsplitdir}v5.${rcp}MR1_pourc_sic_${year}.nc
fileoutpourcsic=${dirout}v5.${rcp}MR1_pourc_sic_${year}_anocorr.nc

ncks -C -O -x -v aire ${fileyeartsoloceRCP} ${ici}tmpyearSST.nc
ncks -C -O -x -v aire ${fileyearpourcsicRCP} ${ici}tmpyearSIC.nc

#--performing anomaly corrections
echo "... for SST..."
cdo add -sub ${ici}tmpyearSST.nc ${tmpavgsstRCP} ${tmpavgsstAMIP} ${fileouttsoloce}
echo "... for SIC, as SIC = SIC(RCP) + <SICLIC(AMIP)> - <SIC(RCP)> - <LIC(RCP)>..."
cdo add -sub ${ici}tmpyearSIC.nc ${tmpavgsicRCP} ${tmpavgsicAMIP} ${ici}tmp.nc
cdo sub ${ici}tmp.nc ${tmpavglicRCP} ${fileoutpourcsic}
rm -f ${ici}tmp.nc
#--clipping of pourc_sic values between 0 and 100
ncap2 -s 'where(pourc_sic>100.) pourc_sic=100;' ${fileoutpourcsic} -O ${ici}tmp.nc
ncap2 -s 'where(pourc_sic<0.) pourc_sic=0;' ${ici}tmp.nc -O ${fileoutpourcsic}
rm -f ${ici}tmp.nc ${ici}tmpyearSST.nc ${ici}tmpyearSIC.nc

#--end loop on years
done

echo "Cleaning up..."
rm -f ${tmpavgsstAMIP} ${tmpavgsicAMIP} ${tmpavgsstRCP} ${tmpavgsicRCP} ${tmpavglicRCP}
rm -rf ${rcpsplitdir}

#--end loop on RCPs
done

#--cleaning up remapped AMIPs
echo "Cleaning up..."
rm -f ${ici}amipbc_sst_2006_2017_remapped.nc
rm -f ${ici}amipbc_sic_2006_2017_remapped.nc
echo "...OK!"