#!/bin/bash
################################################################
### Launching script for a chained simulation of 1D PCM runs ###
###                following orbital changes                 ###
################################################################
set -e
trap 'echo -e "\033[31mError: an issue occurred in the script on line $LINENO! Please review the command and try again.\033[0m"' ERR


################################################################
# Modify here the parameters for the simulation
###############################################
# Path to the arch.env to source:
source ../trunk/LMDZ.COMMON/arch.env

# Set the number of Earth years before present to start:
year_earth_bp_ini=-100

# Set the number of years to be simulated, either Martian or Earth years (> 0):
n_mars_years=100.
#n_earth_years=300

# Name of executable for the PCM:
exePCM="testphys1d_32_phymars_seq.e"

# Name of the orbital data file:
orb_data="obl_ecc_lsp.asc"
################################################################


echo "The launching script is starting!"
echo "The output file is \"launch.log\"."
exec > launch.log 2>&1

# Save the current value of LC_NUMERIC and set it to a locale that uses a dot as the decimal separator
OLD_LC_NUMERIC=$LC_NUMERIC
LC_NUMERIC=en_US.UTF-8

# To end the launching script with error
errlaunch() {
    # Restore the previous value of LC_NUMERIC
    LC_NUMERIC=$OLD_LC_NUMERIC

    date
    echo "Error: an issue occured in the launching script!"
    exit 1
}

# To convert Earth years into Mars years
convertyears() {
    myear=686.9725      # Number of Earth days in Martian year
    eyear=365.256363004 # Number of days in Earth year
    convert_years=$(echo "$myear/$eyear" | bc -l)
    convert_years=$(printf "%.4f" $convert_years) # Rounding to the 4th decimal to respect the precision of Martian year
    if [ -v n_mars_years ]; then
        n_myear=$n_mars_years
        echo "Number of years to be simulated: $n_myear Martian years."
    elif [ -v n_earth_years ]; then
        n_myear=$(echo "$n_earth_years/$convert_years" | bc -l)
        echo "Number of years to be simulated: $n_earth_years Earth years = $n_myear Martian years."
    fi
}

# Check if the set-up is correct
dir=`pwd`
machine=`hostname`
address=`whoami`
if [ -v n_mars_years ] && [ ! -z "$n_mars_years" ]; then
    if [ $(echo "$n_mars_years <= 0." | bc -l) -eq 1 ]; then
        echo "Error: 'n_mars_years' must be > 0!"
        errlaunch
    fi
elif [ -v n_earth_years ] && [ ! -z "$n_earth_years" ]; then
    if [ $(echo "$n_earth_years <= 0." | bc -l) -eq 1 ]; then
        echo "Error: 'n_earth_years' must be > 0!"
        errlaunch
    fi
else
    echo "Error: the number of years to be simulated is not set!"
    errlaunch
fi
if [ ! -f "$exePCM" ]; then
    echo "Error: file \"$exePCM\" does not exist in $dir!"
    errlaunch
fi
if [ ! -f "$orb_data" ]; then
    echo "Error: file \"$orb_data\" does not exist in $dir!"
    errlaunch
fi
if [ ! -f "startfi.nc" ]; then
    echo "Error: file \"startfi.nc\" does not exist in $dir!"
    errlaunch
fi
if [ ! -d "logs" ]; then
    mkdir logs
fi
if [ ! -d "starts" ]; then
    mkdir starts
fi
if [ ! -d "diags" ]; then
    mkdir diags
fi

# Initialization
dir=`pwd`
machine=`hostname`
address=`whoami`
echo "This is a chained simulation for PCM 1D runs in $dir on $machine."
convertyears
i_myear=0
iPCM=1
if [ -f "startfi.nc" ]; then
    cp startfi.nc starts/
fi
if [ -f "start1D.txt" ]; then
    cp start1D.txt starts/
fi

# To modify orbital parameters of a file "startfi.nc"
modify_startfi_orb() {
    # Get the new values for the orbital parameters
    yearlask=$(echo "($year_earth_bp_ini + $i_myear*$convert_years)/1000." | bc -l)
    found=false
    read -r y1 obl1 ecc1 lsp1 < "$orb_data"
    while read -r y2 obl2 ecc2 lsp2; do
        if [ "$first_line" = true ]; then
            first_line=false
            continue # Skip the first line as it's already read outside the loop
        fi
        if [ "$(echo "$y1 >= $yearlask && $yearlask > $y2" | bc)" -eq 1 ]; then
            found=true
            # Calculate linear interpolations for each orbital parameter
            new_obl=$(echo "($obl2 - $obl1)/($y2 - $y1)*($yearlask - $y1) + $obl1" | bc -l)
            new_ecc=$(echo "($ecc2 - $ecc1)/($y2 - $y1)*($yearlask - $y1) + $ecc1" | bc -l)
            if [ "$(echo "$lsp2 - $lsp1 > 300. || $lsp2 - $lsp1 < -300." | bc)" -eq 1 ]; then # If modulo is "crossed" through the interval
                if [ "$(echo "$lsp1 < $lsp2" | bc -l)" -eq 1 ]; then # Lsp should be decreasing
                    lsp1=$(echo "$lsp1 + 360" | bc -l)
                else # Lsp should be increasing
                    lsp2=$(echo "$lsp2 + 360" | bc -l)
                fi
            fi
            new_Lsp=$(echo "($lsp2 - $lsp1)/($y2 - $y1)*($yearlask - $y1) + $lsp1" | bc -l)
            echo "New obliquity    = $new_obl"
            echo "New eccentricity = $new_ecc"
            echo "New Lsp          = $new_Lsp"
            echo "Done!"
            break
        fi
        y1=$y2
        obl1=$obl2
        ecc1=$ecc2
        lsp1=$lsp2
    done < "$orb_data"
    if [ "$found" = false ]; then
        echo "Error: the current year could not be found in the file \"$orb_data\"!"
        exit 1
    fi
    # Compute the modified values in the "starfi.nc"
    year_day=669.
    halfaxe=227.94
    pi=$(echo "4.*a(1.)" | bc -l)
    degrad=$(echo "180./$pi" | bc -l)
    periheli=$(echo "$halfaxe*(1. - $new_ecc)" | bc -l)
    aphelie=$(echo "$halfaxe*(1. + $new_ecc)" | bc -l)
    tan=$(echo "s(0.5*$new_Lsp/$degrad)/c(0.5*$new_Lsp/$degrad)" | bc -l)
    zx0=$(echo "-2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
    if [ $(echo "$zx0 <= 0." | bc -l) -eq 1 ]; then
        zx0=$(echo "$zx0 + 2.*$pi" | bc -l)
    fi
    peri_day=$(echo "$year_day*(1. - ($zx0 - $new_ecc*s($zx0))/(2.*$pi))" | bc -l)
    ztheta=$(echo "($new_iniLs - $new_Lsp)/$degrad" | bc -l)
    tan=$(echo "s(0.5*$ztheta)/c(0.5*$ztheta)" | bc -l)
    zx0=$(echo "2.*a($tan*sqrt((1. - $new_ecc)/(1. + $new_ecc)))" | bc -l)
    xref=$(echo "$zx0 - $new_ecc*s($zx0)" | bc -l)
    # Modify the orbital parameters
      # controle(15) = periheli ! min. Sun-Mars distance (Mkm)  ̃206.66
      # controle(16) = aphelie  ! max. Sun-Mars distance (Mkm)  ̃249.22
      # controle(17) = peri_day ! date of perihelion (sols since N. spring)
      # controle(18) = obliquit ! Obliquity of the planet (deg)  ̃23.98
    ncap2 -O -s "controle(17)=$new_obl" \
             -s "controle(14)=$periheli" \
             -s "controle(15)=$aphelie" \
             -s "controle(16)=$peri_day" \
             "startfi.nc" "startfi.nc"
    echo "New periheli = $periheli"
    echo "New aphelie  = $aphelie"
    echo "New peri_day = $peri_day"
    echo "Done!"
}


# Main loop to to run PCM year by year
while [ $i_myear -lt $n_myear ]; do
    # Run the PCM
    echo "Run $iPCM: year $iPCM/$n_mars_years..."
    ./$exePCM > run.log 2>&1
    if [ ! -f "restartfi.nc" ] || ! (tail -n 100 run.log | grep -iq "everything is cool!"); then # Check if it ended abnormally
        echo "Error: the run $iPCM crashed!"
        errlaunch
    fi
    # Copy data files and prepare the next run
    mv run.log logs/run${iPCM}.log
    if [ -f "diagfi.nc" ]; then
        mv diagfi.nc diags/diagfi${iPCM}.nc
    fi
    if [ -f "diagsoil.nc" ]; then
        mv diagsoil.nc diags/diagsoil${iPCM}.nc
    fi
    if [ -f "stats.nc" ]; then
        mv stats.nc diags/stats${iPCM}.nc
    fi
    if [ -f "Xdiurnalave.nc" ]; then
        mv Xdiurnalave.nc diags/Xdiurnalave${iPCM}.nc
    fi
    cp restartfi.nc starts/startfi${iPCM}.nc
    mv restartfi.nc startfi.nc
    if [ -f "restart1D.txt" ]; then
        cp restart1D.txt starts/restart1D${iPCM}.txt
        mv restart1D.txt start1D.txt
    fi
    ((iPCM++))
    ((i_myear++))
    echo "Done!"
done

# Restore the previous value of LC_NUMERIC
LC_NUMERIC=$OLD_LC_NUMERIC

date
echo "Success: the launching script completed normally!"
exit 0