#!/bin/bash
################################################################
### Launching script for a chained simulation of 1D PCM runs ###
################################################################

echo "The launching script is starting!"
echo "The output file is \"loglaunch.txt\"."
if [ "$1" = "bg" ]; then
    date
else
    nohup "$0" bg > loglaunch.txt 2>&1 &
    exit 1
fi

# A few parameters that might need be changed depending on your setup:
# Path to the arch.env to source:
source ../trunk/LMDZ.COMMON/arch.env

# 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

#--------- Modify here the number of years for initialization ----------
## set the number of Earth years before present to start:
eyears_bp_ini=-1000

#---------- Modify here the number of years to be simulated ------------
## set the number of Martian years:
n_myears=500

#------------------ Modify here the name of PCM exe --------------------
## fill in the name of executable for PCM:
exePCM="testphys1d_29_phymars_para.e"

#----------- Modify here the name of the orbital data file -------------
## fill in the name of executable for PCM:
orb_data="obl_ecc_lsp.asc"

#------ Check if files/directories necessary for the script exist ------
dir=`pwd`
machine=`hostname`
address=`whoami`
if [ ! -f "$exePCM" ]; then
    echo "Error: file \"$exePCM\" does not exist in $dir!"
    exit 1
fi
if [ ! -f "$orb_data" ]; then
    echo "Error: file \"$orb_data\" does not exist in $dir!"
    exit 1
fi
if [ ! -f "startfi.nc" ]; then
    echo "Error: file \"startfi.nc\" does not exist in $dir!"
    exit 1
fi
if [ ! -d "out_PCM" ]; then
    mkdir out_PCM
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 runs in $dir on $machine."
echo "Number of years to be simulated: $n_myears Martian years."
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
i_myear=0
iPCM=1
cp startfi.nc starts/
if [ -f "star1D.nc" ]; then
    cp star1D.txt starts/
fi

#---------------- Main loop to to run PCM year by year -----------------
while [ $i_myear -lt $n_myears ]; do
    # Get the new values for the orbital parameters
    yearlask=$(echo "$eyears_bp_ini + $i_myear*$convert_years" | 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!"
    # Run the PCM
    echo "Run PCM $iPCM: year $iPCM/$n_myears..."
    ./$exePCM > out_runPCM${iPCM} 2>&1
    if [ ! -f "restartfi.nc" ]; then # Check if run ended abnormally
        echo "Error: the run PCM $iPCM has crashed!"
        exit 1
    fi
    # Copy data files and prepare the next run
    mv out_runPCM${iPCM} out_PCM/run${iPCM}
    mv diagfi.nc diags/diagfi${iPCM}.nc
    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 "The launching script has terminated."