Index: trunk/LMDZ.MARS/util/analyse_netcdf.py
===================================================================
--- trunk/LMDZ.MARS/util/analyse_netcdf.py	(revision 3459)
+++ trunk/LMDZ.MARS/util/analyse_netcdf.py	(revision 3459)
@@ -0,0 +1,84 @@
+############################################################
+### Python script to analyse a NetCDF file for debugging ###
+############################################################
+
+### This script gives useful information about a NetCDF file
+### to help for debugging. For each variable, it outputs the
+### dimensions, the min & max values, the average value and
+### warns the user in case of NaN or negative values.
+### The file name is asked to the user in the terminal.
+
+import os
+import readline
+import glob
+from netCDF4 import Dataset
+import numpy as np
+
+############################################################
+### Setup readline for file name autocompletion
+def complete(text,state):
+    line = readline.get_line_buffer().split()
+    # Use glob to find all matching files/directories for the current text
+    if '*' not in text:
+        text += '*'
+    matches = glob.glob(os.path.expanduser(text))
+    # Add '/' if the match is a directory
+    matches = [match + '/' if os.path.isdir(match) else match for match in matches]
+    
+    try:
+        return matches[state]
+    except IndexError:
+        return None
+
+### Function to analyze a variable in a NetCDF file
+def analyze_variable(variable):
+    # Get the data for the variable
+    data = variable[:]
+    
+    # Calculate min, max and mean
+    data_min = np.nanmin(data) # Min value ignoring NaN
+    data_max = np.nanmax(data) # Max value ignoring NaN
+    data_mean = np.nanmean(data) # Mean value ignoring NaN
+    
+    # Check if there are any NaN values
+    has_nan = np.isnan(data).any()
+
+    # Check for negative values
+    has_negative = (data < 0).any()
+    
+    # Print the results
+    print(f"\nAnalysis of variable: {variable.name}")
+    print(f"  Dimensions: {variable.dimensions}")
+    print(f"  Min value : {data_min:>12.6e}")
+    print(f"  Max value : {data_max:>12.6e}")
+    print(f"  Mean value: {data_mean:>12.6e}")
+    if has_nan:
+        print(f"  \033[91mContains NaN values!\033[0m")
+    if has_negative:
+        print(f"  \033[93mWarning: contains negative values!\033[0m")
+
+### Main function
+def analyze_netcdf():
+    # Ask for the file name
+    readline.set_completer(complete)
+    readline.parse_and_bind('tab: complete')
+    file = input("Enter the name of the NetCDF file: ")
+    
+    # Open the NetCDF file
+    try:
+        dataset = Dataset(file,mode='r')
+    except FileNotFoundError:
+        print(f"File '{file}' not found.")
+        return
+    
+    # Iterate through all variables in the dataset to analyze them
+    for variable_name in dataset.variables:
+        variable = dataset.variables[variable_name]
+        analyze_variable(variable)
+    
+    # Close the NetCDF file
+    dataset.close()
+
+### Call the main function
+analyze_netcdf()
+
Index: trunk/LMDZ.MARS/util/display_netcdf.py
===================================================================
--- trunk/LMDZ.MARS/util/display_netcdf.py	(revision 3459)
+++ trunk/LMDZ.MARS/util/display_netcdf.py	(revision 3459)
@@ -0,0 +1,119 @@
+##############################################################
+### Python script to visualize a variable in a NetCDF file ###
+##############################################################
+
+### This script can display any variable of a NetCDF file.
+### The file name, the variable to display and eventually the
+### dimension are asked to the user in the terminal.
+
+import os
+import readline
+import glob
+from netCDF4 import Dataset
+import matplotlib.pyplot as plt
+import numpy as np
+
+##############################################################
+### Setup readline for file name autocompletion
+def complete(text,state):
+    line = readline.get_line_buffer().split()
+    # Use glob to find all matching files/directories for the current text
+    if '*' not in text:
+        text += '*'
+    matches = glob.glob(os.path.expanduser(text))
+    # Add '/' if the match is a directory
+    matches = [match + '/' if os.path.isdir(match) else match for match in matches]
+    
+    try:
+        return matches[state]
+    except IndexError:
+        return None
+
+### Function to handle autocomplete for variable names
+def complete_variable_names(variable_names):
+    def completer(text, state):
+        options = [name for name in variable_names if name.startswith(text)]
+        if state < len(options):
+            return options[state]
+        else:
+            return None
+    return completer
+
+### Function to visualize a variable from a NetCDF file
+def visualize_variable():
+    # Ask for the NetCDF file name
+    readline.set_completer(complete)
+    readline.parse_and_bind('tab: complete')
+    file = input("Enter the name of the NetCDF file: ")
+    
+    # Open the NetCDF file
+    try:
+        dataset = Dataset(file,mode='r')
+    except FileNotFoundError:
+        print(f"File '{file}' not found.")
+        return
+
+    # Display available variables
+    variable_names = list(dataset.variables.keys())
+    print("Available variables:\n",variable_names)
+    
+    # Ask for the variable to display
+    readline.set_completer(complete_variable_names(variable_names))
+    variable_name = input("\nEnter the name of the variable you want to visualize: ")
+    
+    # Check if the variable exists
+    if variable_name not in dataset.variables:
+        print(f"Variable '{variable_name}' not found in the dataset.")
+        dataset.close()
+        return
+    
+    # Extract the selected variable
+    variable = dataset.variables[variable_name][:]
+    
+    # Extract latitude, longitude and altitude
+    latitude = dataset.variables['latitude'][:]
+    longitude = dataset.variables['longitude'][:]
+    
+    # Check if the variable has altitude and time dimensions
+    dimensions = dataset.variables[variable_name].dimensions
+    print(f"\nDimensions of '{variable_name}': {dimensions}")
+    
+    # If the variable has a time dimension, ask for the time index
+    if 'Time' in dimensions:
+        time_index = int(input(f"Enter the time index (0 to {variable.shape[0] - 1}): "))
+    else:
+        time_index = None
+    
+    # If the variable has an altitude dimension, ask for the altitude index
+    if 'altitude' in dimensions:
+        altitude = dataset.variables['altitude'][:]
+        altitude_index = int(input(f"Enter the altitude index (0 to {altitude.shape[0] - 1}): "))
+    else:
+        altitude_index = None
+    
+    # Prepare the 2D slice for plotting
+    if time_index is not None and altitude_index is not None:
+        data_slice = variable[time_index,altitude_index,:,:]
+    elif time_index is not None:
+        data_slice = variable[time_index,:,:]
+    elif altitude_index is not None:
+        data_slice = variable[altitude_index,:,:]
+    else:
+        data_slice = variable[:,:]
+    
+    # Plot the selected variable
+    plt.figure(figsize=(10,6))
+    plt.contourf(longitude,latitude,data_slice,cmap='viridis')
+    plt.colorbar(label=f"{variable_name.capitalize()} (units)") # Adjust units based on your data
+    plt.xlabel('Longitude (degrees)')
+    plt.ylabel('Latitude (degrees)')
+    plt.title(f"{variable_name.capitalize()} visualization")
+    
+    # Show the plot
+    plt.show()
+    
+    # Close the NetCDF file
+    dataset.close()
+
+### Call the main function
+visualize_variable()
Index: trunk/LMDZ.MARS/util/launch_1Dchained.sh
===================================================================
--- trunk/LMDZ.MARS/util/launch_1Dchained.sh	(revision 3459)
+++ trunk/LMDZ.MARS/util/launch_1Dchained.sh	(revision 3459)
@@ -0,0 +1,106 @@
+#!/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 to be simulated ------------
+## set the number of Martian years:
+n_myears=1000
+
+#------------------ Modify here the name of PCM exe --------------------
+## fill in the name of executable for PCM:
+exePCM="testphys1d_29_phymars_seq.e"
+################################################################
+
+
+#------ 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 [ ! -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
+    # 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 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 ended."
Index: trunk/LMDZ.MARS/util/launch_orb_1Dchained.sh
===================================================================
--- trunk/LMDZ.MARS/util/launch_orb_1Dchained.sh	(revision 3459)
+++ trunk/LMDZ.MARS/util/launch_orb_1Dchained.sh	(revision 3459)
@@ -0,0 +1,191 @@
+#!/bin/bash
+################################################################
+### Launching script for a chained simulation of 1D PCM runs ###
+###                following orbital changes                 ###
+################################################################
+
+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)/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!"
+    # 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 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 ended."
