#!/bin/bash

ulimit -s unlimited 
#
EXPNO=100007
EXPNO1=100006 # EXPNO with rerun files
EXPERIMENTSYEAR=1980
EXPERIMENTEYEAR=2010

# Restart:
RESTARTYEAR=1980           # Set to ${EXPERIMENTSYEAR} for initial run 
RESTARTMONTH=01            # Month (01, 02 etc.);  Set to 01 for initial run
#
EXPERIMENTEYEAR1=`expr ${EXPERIMENTEYEAR} + 1`

#set -ex
set -e
export MPLARGS="verbose=1;check=0"
export VPP_MBX_SIZE=64000000
export VPP_STATS=10
#
# Job file to run echam model on euler cluster
# =============================================
#
F_RECLUNIT=BYTE ; export F_RECLUNIT
MPIPROGINF=detail ; export MPIPROGINF
F_SYSLEN=600 ; export F_SYSLEN
#
NCPUS=24
NPROCA=4
NPROCB=6
NPROMA=16
#
RES=42
LEV=39

#
# Pressure levels for output files <EXPNO>_chem_m_<YYYY><MM>.nc
# =============================================================
rm -f afterburner_levels.nml
cat > afterburner_levels.nml << EOL 
&AFTERBURNER_LEVELS
  LEVELS = 1, 5, 10, 30, 50, 100, 160, 250, 400, 630, 1000,  1585, 2512,  3981, 
  6310, 10000, 15000, 20000, 25000, 35000, 45000, 55000, 65000, 75000, 80000, 
  85000, 90000, 95000, 98000, 100000 /
EOL
#
LSOCOL=.TRUE.     # Switch for SOCOL model
LCHEM=.TRUE.      # Switch for chemistry module
#
ECHAM5DIR=${HOME}/SOCOL        	 #location of the model code
WORKDIR=${HOME}/BC     			 #location of the Boundary Conditions
RUNDIR=${HOME}/OUTPUT               #location of the output
DAT=$RUNDIR/run$EXPNO
DPATH=$DAT/

if [ ! -d $DAT ]
then
  mkdir $DAT
fi

#
################################################################################
#
cd $DPATH               # output and rerun files are written into $DPATH

cp $ECHAM5DIR/bin/echam5 $ECHAM5DIR/bin/echam5_${EXPNO}

MODEL=$ECHAM5DIR/bin/echam5  # <--- path to the executable
#
INI=$WORKDIR/T${RES}
#INIAMIP=$WORKDIR/datT${RES}/amip2
INIHADLEY=$WORKDIR/T${RES}/hadley
INISOCOL=$WORKDIR/SOCOL
#
################################################################################
#
# specification of file structures
#
# stop execution after the first run time error
F_ERRCNT=0
export F_ERRCNT
#
################################################################################
#
rm -f unit.?? sst* ice* rrtadata *.codes atmout chem_initial
#
ln -s  ${INI}/T${RES}L${LEV}_jan_spec.nc      unit.23
ln -s  ${INI}/T${RES}_jan_surf.nc             unit.24
#
################################################################################
#
#  Initial file for chemical species (MEZON):

#ln -s  ${INISOCOL}/T42/T42L39_2000_aer_ccmi_spinup.nc  chem_initial
ln -s  ${INISOCOL}/T42/T42L39_ksenia_spinup_m.nc  chem_initial
#
################################################################################
## for climatological sst and ice (LAMIP=F) use:
#ln -s  ${INIAMIP}/T${RES}_amip2sst_clim.nc    unit.20
#ln -s  ${INIAMIP}/T${RES}_amip2sic_clim.nc    unit.96
##
# for AMIP (variable) sst and ice (LAMIP=T) use:

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
    ln -s  ${INIHADLEY}/T${RES}_HadISST_sst_${YY}.nc sst${YY}
    ln -s  ${INIHADLEY}/T${RES}_HadISST_sic_${YY}.nc ice${YY}
    YY=`expr ${YY} + 1`
done
#
################################################################################
#
ln -s  ${INI}/T${RES}_O3clim2.nc    unit.21
ln -s  ${INISOCOL}/T42/misc/OH_ACCENT_T${RES}.nc unit.22
ln -s  ${INI}/T${RES}_VLTCLIM.nc    unit.90
ln -s  ${INI}/T${RES}_VGRATCLIM.nc  unit.91
ln -s  ${INI}/T${RES}_TSLCLIM2.nc   unit.92
ln -s  ${INI}/surrta_data           rrtadata
#
################################################################################
#
# Boundary conditions for SOCOL:

# Greenhouse gases and ODSs:
############################
rm -f co2_y ch4_y n2o_y odsmem_y odscll_y odscls_y odsbr_y

ln -s  ${INISOCOL}/GHG_ODS/co2_socol_1950-2100.txt         co2_y
ln -s  ${INISOCOL}/GHG_ODS/ch4_socol_1950-2100.txt         ch4_y
ln -s  ${INISOCOL}/GHG_ODS/n2o_socol_1950-2100.txt         n2o_y
ln -s  ${INISOCOL}/GHG_ODS/odsmem_socol_1950-2100.txt      odsmem_y
ln -s  ${INISOCOL}/GHG_ODS/odscls_socol_1950-2100.txt      odscls_y
ln -s  ${INISOCOL}/GHG_ODS/odscll_socol_1950-2100.txt      odscll_y
ln -s  ${INISOCOL}/GHG_ODS/odsbr_socol_1950-2100.txt       odsbr_y

# Solar irradiation, Schumann-Runge bands/Lyman alpha line heating 
# parameterisation, lookup-tables for photolysis rates:
#######################################################
rm -f sun_irrad_y sun_par_y photolysis???? photolysis_mean photolysis_delta_e_corr

ln -s ${INISOCOL}/CCMI_solar_data/irrad_echam5/sun_irrad_C1_1950_2011.txt       sun_irrad_y
ln -s ${INISOCOL}/CCMI_solar_data/extra_heating/sun_par_C1_1950-2011.txt           sun_par_y

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
  ln -s ${INISOCOL}/CCMI_solar_data/photolysis/phot_C1_${YY}_nrltest.nc photolysis${YY}
  YY=`expr ${YY} + 1`
done

ln -s ${INISOCOL}/CCMI_solar_data/photolysis/photolysis_socol_wose_mean_1977_1998.nc photolysis_mean # mean 1977-1998 (two solar cycles)

# Delta-E photolysis corrections:
ln -s ${INISOCOL}/sun/L${LEV}_photolysis_delta_e_corr_socol.nc photolysis_delta_e_corr

# Stratospheric aerosols:
#########################
rm -f strataer*

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
    ln -s  ${INISOCOL}/T42/strat_aerosols/T${RES}L${LEV}_sad_nd_ext_omega_g_socol_${YY}.nc strataer${YY}
    YY=`expr ${YY} + 1`
done

ln -s  ${INISOCOL}/T42/strat_aerosols/T${RES}L${LEV}_sad_nd_ext_omega_g_socol_mean_1995_2002.nc strataerbg

# Troposperic aerosol climatology:
##################################
rm -f tropoaer*

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
    ln -s  ${INISOCOL}/T42/tropo_aerosols/T${RES}L${LEV}_GISS_AEROSOLS_${YY}.nc tropoaer${YY}
    YY=`expr ${YY} + 1`
done

ln -s ${INISOCOL}/T42/tropo_aerosols/T${RES}L${LEV}_tropoaero_socol_clim.nc tropoaerclim

# CO and NOx emissions:
#######################
rm -f co_nox_emiss_surf_aircr???? nox_lightning lightn_fac

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
if [ $YY -le 2007 ]; then
    ln -s  ${INISOCOL}/T42/co_nox/T${RES}L${LEV}_co_nox_emiss_surf_aircr_socol_${YY}.nc co_nox_emiss_surf_aircr${YY}
else
    ln -s  ${INISOCOL}/T42/co_nox/T${RES}L${LEV}_co_nox_emiss_surf_aircr_socol_${YY}_ccmi.nc co_nox_emiss_surf_aircr${YY}
fi
    YY=`expr ${YY} + 1`
done
ln -s  ${INISOCOL}/T42/co_nox/T${RES}L${LEV}_nox_lightning_socol.nc nox_lightning

#CCMI
ln -s  ${INISOCOL}/T42/co_nox/T${RES}_fac_lightning_socol.nc lightn_fac
#####################################################################################################
#
#  Input files for tropospheric chemistry:
rm -f  *_emissions  photolysis_rates_messy c5h8_emissions???? ch2o_emissions???? ch3cooh_emissions????

YYY=`expr ${RESTARTYEAR} - 1`


for ((YY=${YYY}; YY<=${EXPERIMENTEYEAR1}; YY++)); do
  if [ $YY -le 2000 ]; then
    ln -s ${INISOCOL}/T42/tropchem/ACCMIP_emissions_Isoprene_T42_${YY}.nc  c5h8_emissions${YY}
    ln -s ${INISOCOL}/T42/tropchem/ACCMIP_emissions_Formaldehyde_T42_${YY}.nc  ch2o_emissions${YY}
  else
    ln -s ${INISOCOL}/T42/tropchem/RCP60/ACCMIP_emissions_RCP60_Isoprene_T42_${YY}.nc  c5h8_emissions${YY}
    ln -s ${INISOCOL}/T42/tropchem/RCP60/ACCMIP_emissions_RCP60_Formaldehyde_T42_${YY}.nc  ch2o_emissions${YY}
  fi
done

for ((YY=${YYY}; YY<=${EXPERIMENTEYEAR1}; YY++)); do
  if [ $YY -le 2008 ]; then
    ln -s ${INISOCOL}/T42/tropchem/ACCMIP_emissions_CH3COOH_T42_${YY}.nc  ch3cooh_emissions${YY}
  else
    ln -s ${INISOCOL}/T42/tropchem/ACCMIP_emissions_CH3COOH_T42_2008.nc  ch3cooh_emissions${YY}
  fi
done
#
ln -s ${INISOCOL}/T42/tropchem/Photolysis_tropchem_T${RES}.nc photolysis_rates_messy
#

#
# SO2, OCS surface emission, DMS sea concentration
###################################################
rm -f so2_surf_emiss* ocs_surf_emiss* dms_sea_con* cs2_surf_emiss* band_echam* so2_cvol_emiss* so2_vol1_emiss*

ln -s ${INISOCOL}/T42/sulfur/band_echam5_pw_socol_condense.dat band_echam5_dryrad

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
    ln -s  ${INISOCOL}/T42/sulfur/T42_SO2_surf_emiss_${YY}_CMIP6.nc  so2_surf_emiss${YY}
#    ln -s  ${INISOCOL}/T42/sulfur/T${RES}L39_OCS_surf_emiss_2000.nc ocs_surf_emiss${YY}       #no file                          

    for MM in 01 02 03 04 05 06 07 08 09 10 11 12; do
    ln -s  ${INISOCOL}/T42/sulfur/T42_so2_emiss_volc_2000${MM}.nc so2_vol1_emiss${YY}${MM}    
done

    YY=`expr ${YY} + 1`
done
ln -s  ${INISOCOL}/T42/sulfur/T${RES}L39_DMS_sea_concen.nc dms_sea_con
ln -s  ${INISOCOL}/T42/sulfur/T${RES}L39_CS2_surf_emiss.nc cs2_surf_emiss
ln -s  ${INISOCOL}/T42/sulfur/T${RES}_so2_cvol_emiss.nc    so2_cvol_emiss

# NMVOC emissions:
########################
rm -f nmvoc_emiss_biogenic???? nmvoc_emiss_bb???? nmvoc_emiss_anthrop????
rm -f nmvoc_emiss_biogenic nmvoc_emiss_bb nmvoc_emiss_anthrop

YY=`expr ${RESTARTYEAR} - 1`

for ((YY=${YY}; YY<=${EXPERIMENTEYEAR1}; YY++)); do
if [ $YY -le 2000 ]; then
    ln -s ${INISOCOL}/T42/nmhc/ACCMIP_emissions_biogenic_NMHC_T42_${YY}.nc  nmvoc_emiss_biogenic${YY}
    ln -s ${INISOCOL}/T42/nmhc/ACCMIP_emissions_anthropogenic_NMHC_T42_${YY}.nc  nmvoc_emiss_anthrop${YY}
    ln -s ${INISOCOL}/T42/nmhc/ACCMIP_emissions_biomassburning_NMHC_T42_${YY}.nc  nmvoc_emiss_bb${YY}
else
    ln -s ${INISOCOL}/T42/nmhc/ACCMIP_emissions_biogenic_NMHC_T42_2000.nc  nmvoc_emiss_biogenic${YY}
    ln -s ${INISOCOL}/T42/nmhc/RCP60/ACCMIP_emissions_anthropogenic_RCP60_NMHC_T42_${YY}.nc  nmvoc_emiss_anthrop${YY}
    ln -s ${INISOCOL}/T42/nmhc/RCP60/ACCMIP_emissions_biomassburning_RCP60_NMHC_T42_${YY}.nc  nmvoc_emiss_bb${YY}
fi
done

# QBO:
######
rm -f qbo????

YY=`expr ${RESTARTYEAR} - 1`
while [ ${YY} -le ${EXPERIMENTEYEAR1} ]; do
    ln -s  ${INISOCOL}/QBO/QBO_socol_${YY}.nc qbo${YY}
    YY=`expr ${YY} + 1`
done

#  Input files for Drydep
rm -f soilpH surf_para
ln -s  ${INISOCOL}/T42/misc/soilpH_T${RES}.nc  soilpH
ln -s  ${INISOCOL}/T42/misc/surf_para_T${RES}.nc  surf_para
#
#File for Berylium production

cp $WORKDIR/10Be/ber_lut_GCR.txt .
cp $WORKDIR/10Be/ber7_lut_GCR.txt .
cp $WORKDIR/10Be/be10_spe69.txt . 
cp $WORKDIR/10Be/be7_spe69.txt . 
cp $WORKDIR/10Be/Pc_1997.txt .
cp $WORKDIR/10Be/phi_1937-2019.txt .

#################################################################################
# Prepare rerun files (if necessary):

if [ ${RESTARTYEAR} -gt ${EXPERIMENTSYEAR} -o ${RESTARTMONTH} -gt 01 ]; then
    if [ ${RESTARTMONTH} -eq 01 ]; then
	MM1=12
	YY1=`expr ${RESTARTYEAR} - 1`
    else
	MM1=`expr ${RESTARTMONTH} - 1`
	YY1=${RESTARTYEAR}
    fi
    if [ ${MM1} -le 9 ]; then
	MM1=0${MM1}
    fi

    #here you can specify the path for the reference restart files
    if [ ${EXPNO1} -ne ${EXPNO} ]; then
       cp -f $HOME/OUTPUT/run${EXPNO1}/${EXPNO1}_rerun_${YY1}${MM1}_echam.nc.gz $DPATH/${EXPNO}_rerun_${YY1}${MM1}_echam.nc.gz
       cp -f $HOME/OUTPUT/run${EXPNO1}/${EXPNO1}_rerun_${YY1}${MM1}_chem1.nc.gz $DPATH/${EXPNO}_rerun_${YY1}${MM1}_chem1.nc.gz
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_aero2 $DPATH/rerun_${EXPNO}_aero2
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_chem2 $DPATH/rerun_${EXPNO}_chem2
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_tracer $DPATH/rerun_${EXPNO}_tracer
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_chem_m $DPATH/rerun_${EXPNO}_chem_m
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_echam $DPATH/rerun_${EXPNO}_echam
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_photo $DPATH/rerun_${EXPNO}_photo 
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_drydep $DPATH/rerun_${EXPNO}_drydep 
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_drydep2 $DPATH/rerun_${EXPNO}_drydep2
       cp -f $HOME/OUTPUT/run${EXPNO1}/rerun_${EXPNO1}_wetdep $DPATH/rerun_${EXPNO}_wetdep
    fi
    
#     a) ECHAM5:
    if [ ! -s ${EXPNO}_rerun_${YY1}${MM1}_echam.nc.gz ]; then
     	echo ${EXPNO}_rerun_${YY1}${MM1}_echam.nc.gz does not exist!
    	exit
    fi
    cp -f ${EXPNO}_rerun_${YY1}${MM1}_echam.nc.gz rerun_${EXPNO}_echam.gz
    gunzip -f rerun_${EXPNO}_echam.gz
 
    if [ ${LSOCOL} = .TRUE. -a ${LCHEM} = .TRUE. ]; then 
#     b) MEZON:
    	if [ ! -s ${EXPNO}_rerun_${YY1}${MM1}_chem1.nc.gz ]; then
    	    echo ${EXPNO}_rerun_${YY1}${MM1}_chem1.nc.gz does not exist!
    	    exit
    	fi
    	cp -f ${EXPNO}_rerun_${YY1}${MM1}_chem1.nc.gz rerun_${EXPNO}_tracer.gz
    	gunzip -f rerun_${EXPNO}_tracer.gz

    fi

    RERUN=.TRUE.     # Rerun switch; .false. for initial run, .true. for rerun
else
    RERUN=.FALSE.
fi

#    RERUN=.FALSE.

rm -f fortran_error_messages

#
#######################################################################################################################################################
#
# Set namelist and call model:  

YY=${RESTARTYEAR}

while [ ${YY} -le ${EXPERIMENTEYEAR} ]; do

echo $EXPERIMENTEYEAR
    echo 'rabotaet'

    for ((MM=${RESTARTMONTH}; MM<=12; MM++)); do
       MM2=${MM}
       if [ ${MM} -le 9 ]; then
	   MM2=0${MM}
fi

       echo ${MM}!!!!!!!!!!!! ${EXPNO}_${YY}${MM2}.01_chem_m.nc

	if [ ${YY} -gt ${RESTARTYEAR} -o ${MM} -ge ${RESTARTMONTH} ]; then

	    if [ -s namelist.echam ]; then
		rm -f namelist.echam
	fi
		
#  namelist control variables and output control for grid variables
#  spectral variables are written out by default except liq. water
#  for production runs set LABORT=.FALSE.
#
	    cat > namelist.echam << EOF
&RUNCTL
  LSOCOL=${LSOCOL}
  LRESUME=$RERUN,
  OUT_DATAPATH = "$DPATH"
  OUT_EXPNAME  = "$EXPNO"
  OUT_FILETYPE = 2
  TRAC_FILETYPE = 2
  DT_START  = ${EXPERIMENTSYEAR},01,01,12,0,0
  DT_STOP   = ${EXPERIMENTEYEAR},01,01,12,0,0
  NO_CYCLES = 1
  PUTDATA   = 6, 'hours', 'first', 0
  PUTRERUN  = 1, 'months', 'last', 0
  DELTA_TIME = 900
  LAMIP=.TRUE.
  LMIDATM=.TRUE.
  LABORT=.FALSE.
  NPROCA=${NPROCA}
  NPROCB=${NPROCB}
  NPROMA=${NPROMA}
/
&SOCOLCTL
  LCHEM=${LCHEM}
  lradcoup=.true.
  lo3_coupl=.true.
  lch4_coupl=.true.
  ln2o_coupl=.true.
  lodscl_coupl=.true.
  lh2o_coupl=.true.
  lsrb_lya_heating=.true.
  nastep=20
  CYEAR=2000
  lqbonudg=.true.
  lhetchem=.true.
  lvolcano=.false.
  linteractive_drydep=.true.
  linteractive_wetdep=.true.
  lscav_gas=.true.
  lscav_cv=.true.
  lscav_ls=.true.
  lscav_aer=.true.
  iscav_easy=3
  coeff_para=2
  interactivelnox=.true.
  nlev_fbc=1
  nlev_mbc=1
  nlev_drydep=1
  ln2o5hydro = .true.  
  lgn2o5_const = .false.
  ln2o5_dry = .true.
/
&PHYSCTL
  LCOND= .true.
  LCONV= .true.
/
&DYNCTL
  VCHECK=235.
  SPDRAG=0.926E-4
/
&RADCTL
  IAERO=10
  LGADSRH=.TRUE.
/
EOF


	    # Call model:
	    date
	    #$MODEL
	    mpirun $MODEL


	    # Save restart files:
            # a) ECHAM5:
	    cp -f rerun_${EXPNO}_echam ${EXPNO}_rerun_${YY}${MM2}_echam.nc
	    gzip -f ${EXPNO}_rerun_${YY}${MM2}_echam.nc

            # b) MEZON:
	    if [ ${LSOCOL} = .TRUE. -a ${LCHEM} = .TRUE. ]; then
		cp -f rerun_${EXPNO}_tracer ${EXPNO}_rerun_${YY}${MM2}_chem1.nc
		gzip -f ${EXPNO}_rerun_${YY}${MM2}_chem1.nc
	    fi

	    RERUN=.true.

           #conversion to monthly means    
           cdo monmean ${EXPNO}_${YY}${MM2}.01_aero2.nc ${EXPNO}_${YY}${MM2}.01_aero2_m.nc
           rm ${EXPNO}_${YY}${MM2}.01_aero2.nc
           cdo monmean ${EXPNO}_${YY}${MM2}.01.nc ${EXPNO}_${YY}${MM2}.01_m.nc
           rm ${EXPNO}_${YY}${MM2}.01.nc
           cdo monmean ${EXPNO}_${YY}${MM2}.01_tracer.nc ${EXPNO}_${YY}${MM2}.01_tracer_m.nc
           rm ${EXPNO}_${YY}${MM2}.01_tracer.nc
           echo aero2 is converted to monthly

	fi

    done

    YY=`expr ${YY} + 1`
    RESTARTMONTH=1

done

ls -lt

