Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2201)
+++ /trunk/LMDZ.MARS/README	(revision 2202)
@@ -2803,2 +2803,6 @@
 == 18/12/2019 == MV
 Update of the rocketduststorm parametrization: clean-up of the code + back to the sub-timesteps method for the Van Leer transport, homogeneous with the one used in topmons_mod.F.
+
+== 20/12/2019 == AB
+Add the observer simulator program simu_MCS_temp.F90 to the LMDZ.MARS/util/ directory, as well as a simu_MCS_temp.def. This program reads a MCS file (binned by Luca Montabone), extracts the daytime and nighttime temperatures from a GCM simulation file, and outputs a netcdf file in the same format as the MCS file.
+More details in simu_MCS_temp.def and in the preface of simu_MCS_temp.F90.
Index: /trunk/LMDZ.MARS/util/simu_MCS_temp.F90
===================================================================
--- /trunk/LMDZ.MARS/util/simu_MCS_temp.F90	(revision 2202)
+++ /trunk/LMDZ.MARS/util/simu_MCS_temp.F90	(revision 2202)
@@ -0,0 +1,1663 @@
+program simu_MCS_temp
+! program for a MCS observer simulator of the temperature data
+! author : Antoine Bierjon, June-December 2019
+! PLEASE GIVE ME YOUR FEEDBACKS ON THIS PROGRAM ! :D
+!
+!===================================================================================================
+!     PREFACE
+!===================================================================================================
+! This program loads Observer's data file (MCS-type, binned by Luca Montabone), serving as a 
+! reference, then extract the dayside and nightside temperature data from either :
+!       - a GCM simulation (diagfi.nc) that covers the observation period,
+!       - a GCM stats file (stats.nc) after extending it to cover the observation period.
+!
+! Since the MCS data is binned by 5°-intervals of Ls, the GCM simulation's data is interpolated at 
+! the Observer's spatial coordinates, at every GCM sol value in each 5°-interval of Ls, before 
+! being averaged to create output bins for each interval. The binning in sols also takes into
+! account the variability of Local Time in the bin and try to represent it too.
+!
+! Eventually, the program outputs a netcdf file, filled with the output bin's data and edited with
+! the same format than the MCS file (and with the same dimensions' declaration order).
+!
+! Minimal requirements :
+!     - the MCS file must contain the variables dtemp,dtimeave,dtimemax,dtimemin,
+!                                               ntemp,ntimeave,ntimemax,ntimemin
+!     - the GCM file must : contain the variable temp (atmospheric temperature) ;
+!                           have the same altitude type as the MCS file ;
+!                           cover completely the NON-BINNED observation period
+!     - if the GCM file is a stats file, please name it with "stats" in it (to be changed in the future)
+!
+! See also NOTA BENE in section 1.3
+! More details can also be found in my internship report (in French) :
+!           /planeto/abierjon/rapportstage2019/Rapport_de_PFE_Master_Bierjon_2019.pdf
+! or directly in my office !
+!===================================================================================================
+
+
+use netcdf
+
+!===================================================================================================
+! 0. Variable declarations
+!===================================================================================================
+
+implicit none ! for no implicitly typed variables
+
+! Files:
+character(len=256) :: gcmfile ! GCM simulation file
+logical :: is_stats = .false. ! to check if the GCM file is a stats.nc or a diagfi.nc file
+character(len=256) :: obsfile ! observation file = MCS/Observer data file
+character(len=256) :: outfile ! output file
+
+! NetCDF stuff
+integer :: status ! NetCDF routines return code
+character (len=256) :: error_text ! to store the error text
+integer :: gcmfid ! NetCDF gcm file ID
+integer :: obsfid ! NetCDF observation file ID
+integer :: outfid ! NetCDF output file ID
+integer :: GCMvarid, OBSvarid, LT_id, outvarid ! to store the ID of a variable
+integer :: lat_dimid_obs,lon_dimid_obs,alt_dimid_obs,time_dimid_obs ! dimensions' ID in OBS and output files
+integer :: lat_dimid_gcm,lon_dimid_gcm,alt_dimid_gcm,time_dimid_gcm ! dimensions' ID in GCM file
+integer :: GCMvarshape(4), OBSvarshape(4), LTshape(3) ! to store a variable's coordinates order
+integer :: nb
+
+real,dimension(:),allocatable :: GCMlon, OBSlon ! longitude in the GCM & the Observer files
+integer GCMlonlen, OBSlonlen ! # of grid points along GCMlon & OBSlon
+
+real,dimension(:),allocatable :: GCMlat, OBSlat ! latitude in the GCM & the Observer files
+integer GCMlatlen, OBSlatlen ! # of grid points along GCMlat & OBSlat
+
+real,dimension(:),allocatable :: GCMalt, OBSalt ! can be geometric heights or pressure
+integer GCMaltlen, OBSaltlen ! # of grid point along GCMalt & OBSalt
+character (len=256) :: units ! to store the units attribute of altitude
+character(len=1) :: GCMalttype, OBSalttype ! altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
+
+real,dimension(:),allocatable :: GCMtime, OBSLs ! time in the GCM diagfi (sols) & the Observer files (Ls)
+real,dimension(:),allocatable :: GCMstatstime ! time in the GCM stats file (LT at lon 0°)
+integer GCMtimelen, GCMstatstimelen, OBSLslen ! # of points along GCMtime, GCMstatstime, OBSLs
+real :: starttimeoffset=0 ! offset (in sols) wrt Ls=0 of sol 0 in GCM file
+
+character (len=64) :: GCMvarname, OBSvarname ! to store the name of the variable to retrieve
+real,dimension(:,:,:,:),allocatable :: GCM_var,OBS_var ! the 4D variable extracted from GCM & OBS files
+character (len=64) :: OBSLTave_name, OBSLTmax_name, OBSLTmin_name ! to store the name of the average, max and min of LT (ex : dtimeave, ntimemax...)
+real,dimension(:,:,:),allocatable :: OBSLT, OBSLTmax, OBSLTmin ! 3D variables extracted from obsfile (average, max and min of LT in a bin)
+real :: GCMmiss_val, OBSmiss_val, LTmiss_val ! value to denote non-existant data
+real :: extr_value ! to store the result of the extraction subroutine
+
+real :: OBSdeltaLs ! difference of Ls between each observation bin
+real :: sol, maxsol, minsol ! sol date corresponding to Observed Ls and LT
+integer :: m_maxsol, m_minsol ! indexes of the maximum and minimum GCM sol taken in a bin for interpolation
+
+real,dimension(:,:,:),allocatable :: LTmod_arraymax, LTmod_arraymin ! declaration of the type of the function LTmod_array
+real :: LTmod ! declaration of the type of the function LTmod
+real :: maxLTvar ! maximum variability of LT within one MCS bin
+integer :: LTcount ! nb of LT samples on which the interpolation is performed, for every LT interval
+
+integer :: solcount ! number of GCM sol integer values in one Ls interval
+real,dimension(:),allocatable :: int_sol_list ! list of the integer values of GCM sol
+real,dimension(:),allocatable :: sol_list ! list of the sol values used for the interpolation
+integer :: solerrcount ! nb of GCM missing values during interpolation, which are removed for the binning 
+integer :: errcount = 0 ! total number of GCM missing values
+real :: solbinned_value ! to store the extracted value averaged on sol samples, which is finally put in the output bin
+real :: GCMdeltat ! timestep of the GCM simulation (in Martian hours)
+
+real :: lon_val, lat_val, alt_val, Ls_val, LT_val ! where and when the output file is written at
+
+integer :: i,j,k,l ! loop iteration indexes
+integer :: m,m_int,n ! sol binning loops iteration index
+
+real, dimension(:,:,:,:), allocatable :: outvar ! outvar(,,,): 4D array to store the output variable's data
+
+!===================================================================================================
+! 1. Read the observation file (Observer data)
+!===================================================================================================
+!===============================================================================
+! 1.1 Open the Observer data file obsfile
+!===============================================================================
+! Ask the user to give a netcdf observation file
+write(*,*) "-> Enter observation file name :"
+read(*,*) obsfile
+
+status=nf90_open(obsfile,nf90_nowrite,obsfid) ! nowrite mode=the program can only read the opened file
+error_text="Error: could not open file "//trim(obsfile)
+call status_check(status,error_text)
+
+!===============================================================================
+! 1.2 Get grids for dimensions lon,lat,alt,time from the observation file
+!===============================================================================
+! OBS Latitude
+status=nf90_inq_dimid(obsfid,"latitude",lat_dimid_obs)
+error_text="Failed to find Observer latitude dimension"
+call status_check(status,error_text)
+
+status=nf90_inquire_dimension(obsfid,lat_dimid_obs,len=OBSlatlen)
+error_text="Failed to find Observer latitude length"
+call status_check(status,error_text)
+allocate(OBSlat(OBSlatlen))
+
+status=nf90_inq_varid(obsfid,"latitude",OBSvarid)
+error_text="Failed to find Observer latitude ID"
+call status_check(status,error_text)
+
+! Read OBSlat
+status=nf90_get_var(obsfid,OBSvarid,OBSlat)
+error_text="Failed to load OBSlat"
+call status_check(status,error_text)
+
+
+! OBS Longitude
+status=nf90_inq_dimid(obsfid,"longitude",lon_dimid_obs)
+error_text="Failed to find Observer longitude dimension"
+call status_check(status,error_text)
+
+status=nf90_inquire_dimension(obsfid,lon_dimid_obs,len=OBSlonlen)
+error_text="Failed to find Observer longitude length"
+call status_check(status,error_text)
+allocate(OBSlon(OBSlonlen))
+
+status=nf90_inq_varid(obsfid,"longitude",OBSvarid)
+error_text="Failed to find Observer longitude ID"
+call status_check(status,error_text)
+
+! Read GCMlon
+status=nf90_get_var(obsfid,OBSvarid,OBSlon)
+error_text="Failed to load OBSlon"
+call status_check(status,error_text)
+
+
+! OBS Time (Ls)
+status=nf90_inq_dimid(obsfid,"time",time_dimid_obs)
+error_text="Failed to find Observer time (Ls) dimension"
+call status_check(status,error_text)
+
+status=nf90_inquire_dimension(obsfid,time_dimid_obs,len=OBSLslen)
+error_text="Failed to find Observer time (Ls) length"
+call status_check(status,error_text)
+allocate(OBSLs(OBSLslen))
+
+status=nf90_inq_varid(obsfid,"time",OBSvarid)
+error_text="Failed to find Observer time (Ls) ID"
+call status_check(status,error_text)
+
+! Read OBSLs
+status=nf90_get_var(obsfid,OBSvarid,OBSLs)
+error_text="Failed to load OBSLs"
+call status_check(status,error_text)
+
+! Get the observation timestep between bins
+OBSdeltaLs=OBSLs(2)-OBSLs(1)
+
+! OBS Altitude
+status=nf90_inq_dimid(obsfid,"altitude",alt_dimid_obs)
+error_text="Failed to find Observer altitude dimension"
+call status_check(status,error_text)
+
+status=nf90_inquire_dimension(obsfid,alt_dimid_obs,len=OBSaltlen)
+error_text="Failed to find Observer altitude length"
+call status_check(status,error_text)
+allocate(OBSalt(OBSaltlen))
+
+status=nf90_inq_varid(obsfid,"altitude",OBSvarid)
+error_text="Failed to find Observer altitude ID"
+call status_check(status,error_text)
+
+! Read OBSalt
+status=nf90_get_var(obsfid,OBSvarid,OBSalt)
+error_text="Failed to load OBSalt"
+call status_check(status,error_text)
+
+! Check altitude attribute "units" to find out altitude type and compare with the GCM file
+status=nf90_get_att(obsfid,OBSvarid,"units",units)
+error_text="Failed to load Observer altitude units attribute"
+call status_check(status,error_text)
+! an unknown and invisible character is placed just after the unit's 
+! characters in the Observer file so we only take the first characters
+! corresponding to the sought unit
+if (trim(units(1:2)).eq."Pa") then
+  units="Pa"
+  OBSalttype='p' ! pressure coordinate
+else if (trim(units(1:2)).eq."m") then
+  units="m"
+  OBSalttype='z' ! altitude coordinate
+else
+  write(*,*)" I do not understand this unit ",trim(units)," for Observer altitude!"
+  stop
+endif
+
+!===============================================================================
+! 1.3 Extraction of the wanted variables from the observation file
+!===============================================================================
+
+! Observer variables to extract :
+!******************** NOTA BENE (cf sections 1.3 and 4)*************************
+! We execute the program a first time with the daytime values, and then a second
+! time with the nighttime values. However, a lot of instructions below are
+! independent of this distinction, hence we execute them only in the first loop,
+! when OBSvarname="dtemp". It also implies a lot of IF (when it's only executed in 
+! the first loop) or CASE (when it's executed in both loops but depends on wether 
+! it is daytime or nighttime).
+! This should be changed one day with a cleaner version of the code...
+!*******************************************************************************
+OBSvarname="dtemp" ! we begin with daytime temperature
+write(*,*)"" ; write(*,*) "Beginning the 1st loop, on daytime values"; write(*,*)""
+DAY_OR_NIGHT: DO ! (the end of the loop is in section 4.)
+
+  ! Check that the observation file contains that variable
+  status=nf90_inq_varid(obsfid,trim(OBSvarname),OBSvarid)
+  error_text="Failed to find variable "//trim(OBSvarname)//" in "//trim(obsfile)
+  call status_check(status,error_text)
+
+  ! Sanity checks on the variable
+  status=nf90_inquire_variable(obsfid,OBSvarid,ndims=nb,dimids=OBSvarshape)
+  error_text="Failed to obtain information on variable "//trim(OBSvarname)
+  call status_check(status,error_text)
+
+  ! Check that it is a 4D variable
+  if (nb.ne.4) then
+    write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for ", trim(OBSvarname)
+    stop
+  endif
+  ! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
+  if (OBSvarshape(1).ne.lon_dimid_obs) then
+    write(*,*) "Error, expected first dimension of ", trim(OBSvarname), " to be longitude!"
+    stop
+  endif
+  if (OBSvarshape(2).ne.lat_dimid_obs) then
+    write(*,*) "Error, expected second dimension of ", trim(OBSvarname), " to be latitude!"
+    stop
+  endif
+  if (OBSvarshape(3).ne.alt_dimid_obs) then
+    write(*,*) "Error, expected third dimension of ", trim(OBSvarname), " to be altitude!"
+    stop
+  endif
+  if (OBSvarshape(4).ne.time_dimid_obs) then
+    write(*,*) "Error, expected fourth dimension of ", trim(OBSvarname), " to be time!"
+    stop
+  endif
+
+  !write(*,*)"Dimensions ID in the obsfile are :"  
+  !write(*,*)"lon_dimid=",lon_dimid_obs
+  !write(*,*)"lat_dimid=",lat_dimid_obs
+  !write(*,*)"alt_dimid=",alt_dimid_obs
+  !write(*,*)"time_dimid=",time_dimid_obs
+
+  ! Length allocation for each dimension of the 4D variable
+  allocate(OBS_var(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+
+  ! Load datasets
+  status=nf90_get_var(obsfid,OBSvarid,OBS_var)
+  error_text="Failed to load "//trim(OBSvarname)//" from the obsfile"
+  call status_check(status,error_text)
+  write(*,*) "Loaded ",trim(OBSvarname)," from the obsfile"
+
+  ! Get OBS_var missing_value attribute
+  status=nf90_get_att(obsfid,OBSvarid,"_FillValue",OBSmiss_val)
+  error_text="Failed to load missing_value attribute"
+  call status_check(status,error_text)
+  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
+
+!===============================================================================
+! 1.4 Extraction of the OBS Local Time (LT)
+!===============================================================================
+  SELECT CASE (OBSvarname)
+  CASE ("dtemp")
+    OBSLTave_name = "dtimeave"
+    OBSLTmax_name = "dtimemax"
+    OBSLTmin_name = "dtimemin"
+  CASE ("ntemp")
+    OBSLTave_name = "ntimeave"
+    OBSLTmax_name = "ntimemax"
+    OBSLTmin_name = "ntimemin"
+  END SELECT
+  
+  status=nf90_inq_varid(obsfid,trim(OBSLTave_name),LT_id)
+  error_text="Failed to find Observer local time ("//trim(OBSLTave_name)//") dimension"
+  call status_check(status,error_text)
+  ! Sanity checks on the variable
+  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
+  error_text="Failed to obtain information on variable "//trim(OBSLTave_name)
+  call status_check(status,error_text)
+  ! Check that it is a 3D variable
+  if (nb.ne.3) then
+    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of "//trim(OBSLTave_name)//"!"
+    stop
+  endif
+  ! Check that its dimensions are indeed lon,lat,time (in the right order)
+  if (LTshape(1).ne.lon_dimid_obs) then
+    write(*,*) "Error, expected first dimension of ", trim(OBSLTave_name), " to be longitude!"
+    stop
+  endif
+  if (LTshape(2).ne.lat_dimid_obs) then
+    write(*,*) "Error, expected second dimension of ", trim(OBSLTave_name), " to be latitude!"
+    stop
+  endif
+  if (LTshape(3).ne.time_dimid_obs) then
+    write(*,*) "Error, expected third dimension of ", trim(OBSLTave_name), " to be time!"
+    stop
+  endif
+  ! Length allocation for each dimension of the 3D variable
+  allocate(OBSLT(OBSlonlen,OBSlatlen,OBSLslen))
+  ! Load datasets
+  status=nf90_get_var(obsfid,LT_id,OBSLT)
+  error_text="Failed to load "//trim(OBSLTave_name)//" from the obsfile"
+  call status_check(status,error_text)
+  write(*,*) "Loaded ", trim(OBSLTave_name), " from the obsfile"
+  ! Get LT missing_value attribute
+  status=nf90_get_att(obsfid,LT_id,"_FillValue",LTmiss_val)
+  error_text="Failed to load missing_value attribute"
+  call status_check(status,error_text)
+  write(*,'("  with missing_value attribute : ",1pe12.5)')LTmiss_val
+
+
+  ! Get the maximum values of LT in the OBS bins
+  status=nf90_inq_varid(obsfid,trim(OBSLTmax_name),LT_id)
+  error_text="Failed to find Observer max local time ("//trim(OBSLTmax_name)//") dimension"
+  call status_check(status,error_text)
+  ! Sanity checks on the variable
+  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
+  error_text="Failed to obtain information on variable "//trim(OBSLTmax_name)
+  call status_check(status,error_text)
+  ! Check that it is a 3D variable
+  if (nb.ne.3) then
+    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of ", trim(OBSLTmax_name), "!"
+    stop
+  endif
+  ! Check that its dimensions are indeed lon,lat,time (in the right order)
+  if (LTshape(1).ne.lon_dimid_obs) then
+    write(*,*) "Error, expected first dimension of ", trim(OBSLTmax_name), " to be longitude!"
+    stop
+  endif
+  if (LTshape(2).ne.lat_dimid_obs) then
+    write(*,*) "Error, expected second dimension of ", trim(OBSLTmax_name), " to be latitude!"
+    stop
+  endif
+  if (LTshape(3).ne.time_dimid_obs) then
+    write(*,*) "Error, expected third dimension of ", trim(OBSLTmax_name), " to be time!"
+    stop
+  endif
+  ! Length allocation for each dimension of the 3D variable
+  allocate(OBSLTmax(OBSlonlen,OBSlatlen,OBSLslen))
+  ! Load datasets
+  status=nf90_get_var(obsfid,LT_id,OBSLTmax)
+  error_text="Failed to load "//trim(OBSLTmax_name)//" from the obsfile"
+  call status_check(status,error_text)
+  write(*,*) "Loaded ", trim(OBSLTmax_name), " from the obsfile"
+
+
+  ! Get the minimum values of LT in the OBS bins
+  status=nf90_inq_varid(obsfid,trim(OBSLTmin_name),LT_id)
+  error_text="Failed to find Observer min local time ("//trim(OBSLTmin_name)//") dimension"
+  call status_check(status,error_text)
+  ! Sanity checks on the variable
+  status=nf90_inquire_variable(obsfid,LT_id,ndims=nb,dimids=LTshape)
+  error_text="Failed to obtain information on variable "//trim(OBSLTmin_name)
+  call status_check(status,error_text)
+  ! Check that it is a 3D variable
+  if (nb.ne.3) then
+    write(*,*) "Error, expected a 3D (lon-lat-time) variable for of ", trim(OBSLTmin_name), "!"
+    stop
+  endif
+  ! Check that its dimensions are indeed lon,lat,time (in the right order)
+  if (LTshape(1).ne.lon_dimid_obs) then
+    write(*,*) "Error, expected first dimension of ", trim(OBSLTmin_name), " to be longitude!"
+    stop
+  endif
+  if (LTshape(2).ne.lat_dimid_obs) then
+    write(*,*) "Error, expected second dimension of ", trim(OBSLTmin_name), " to be latitude!"
+    stop
+  endif
+  if (LTshape(3).ne.time_dimid_obs) then
+    write(*,*) "Error, expected third dimension of ", trim(OBSLTmin_name), " to be time!"
+    stop
+  endif
+  ! Length allocation for each dimension of the 3D variable
+  allocate(OBSLTmin(OBSlonlen,OBSlatlen,OBSLslen))
+  ! Load datasets
+  status=nf90_get_var(obsfid,LT_id,OBSLTmin)
+  error_text="Failed to load "//trim(OBSLTmin_name)//" from the obsfile"
+  call status_check(status,error_text)
+  write(*,*) "Loaded ", trim(OBSLTmin_name), " from the obsfile"
+
+  
+  SELECT CASE (OBSvarname)
+  CASE ("dtemp")
+    maxLTvar=maxval((OBSLTmax(:,:,:)-OBSLTmin(:,:,:)))
+    write(*,*)"The max daytime LT variability is ",maxLTvar,"hours (within a same bin)"
+  CASE ("ntemp")
+    ! in case of night local times, we replace hours from a [0;24[ interval to a 
+    ! [-12;12[ interval in which midnight = 0, with the subroutine LTmod
+    allocate(LTmod_arraymax(OBSlonlen,OBSlatlen,OBSLslen))
+    allocate(LTmod_arraymin(OBSlonlen,OBSlatlen,OBSLslen))
+    do i=1,OBSlonlen
+      do j=1,OBSlatlen
+        do l=1,OBSLslen
+          LTmod_arraymax(i,j,l) = LTmod(OBSLTmax(i,j,l))
+          LTmod_arraymin(i,j,l) = LTmod(OBSLTmin(i,j,l))
+        enddo
+      enddo
+    enddo
+    maxLTvar=maxval(LTmod_arraymax(:,:,:) - LTmod_arraymin(:,:,:))
+    write(*,*)"The max nighttime LT variability is ",maxLTvar,"hours (within a same bin)"
+  END SELECT
+  
+  LTcount=floor(maxLTvar) ! integer nb of hours fitting in the broadest interval of LT
+
+!===================================================================================================
+! 2. Read the input file (GCM simulation)
+!===================================================================================================
+  IF (OBSvarname.EQ."dtemp") THEN ! reading the GCM file is done only in the first loop
+  
+!===============================================================================
+! 2.1 Open the GCM simulation file gcmfile
+!===============================================================================
+    ! Ask the user to give a netcdf input file
+    write(*,*)"";write(*,*) "-> Enter input file name (GCM simulation) :"
+    read(*,*) gcmfile
+    if (index(gcmfile,"stats").ne.0) then ! if the GCM file name contains "stats", we assume it is a stats file
+    ! => A REMPLACER PAR UN TEST SUR LA DIMENSION TIME ? REGARDER DANS AUTRES UTILITAIRES 
+      is_stats = .true.
+      write(*,*)"The GCM file is recognized as a stats file."
+    else
+      write(*,*)"The GCM file is recognized as a diagfi/concatnc file."
+    endif
+
+    ! Open GCM file
+    status=nf90_open(gcmfile,nf90_nowrite,gcmfid)
+    ! nowrite mode=the program can only read the opened file
+    error_text="Failed to open datafile "//trim(gcmfile)
+    call status_check(status,error_text)
+
+!===============================================================================
+! 2.2 Get grids for dimensions lon,lat,alt,time from the GCM file
+!===============================================================================
+    ! GCM Latitude
+    status=nf90_inq_dimid(gcmfid,"latitude",lat_dimid_gcm)
+    error_text="Failed to find GCM latitude dimension"
+    call status_check(status,error_text)
+
+    status=nf90_inquire_dimension(gcmfid,lat_dimid_gcm,len=GCMlatlen)
+    error_text="Failed to find GCM latitude length"
+    call status_check(status,error_text)
+    allocate(GCMlat(GCMlatlen))
+
+    status=nf90_inq_varid(gcmfid,"latitude",GCMvarid)
+    error_text="Failed to find GCM latitude ID"
+    call status_check(status,error_text)
+
+    ! Read GCMlat
+    status=nf90_get_var(gcmfid,GCMvarid,GCMlat)
+    error_text="Failed to load GCMlat"
+    call status_check(status,error_text)
+
+
+    ! GCM Longitude
+    status=nf90_inq_dimid(gcmfid,"longitude",lon_dimid_gcm)
+    error_text="Failed to find GCM longitude dimension"
+    call status_check(status,error_text)
+
+    status=nf90_inquire_dimension(gcmfid,lon_dimid_gcm,len=GCMlonlen)
+    error_text="Failed to find GCM longitude length"
+    call status_check(status,error_text)
+    allocate(GCMlon(GCMlonlen))
+
+    status=nf90_inq_varid(gcmfid,"longitude",GCMvarid)
+    error_text="Failed to find GCM longitude ID"
+    call status_check(status,error_text)
+
+    ! Read GCMlon
+    status=nf90_get_var(gcmfid,GCMvarid,GCMlon)
+    error_text="Failed to load GCMlon"
+    call status_check(status,error_text)
+
+
+    ! GCM Time
+    status=nf90_inq_dimid(gcmfid,"Time",time_dimid_gcm)
+    error_text="Failed to find GCM time (sols) dimension"
+    call status_check(status,error_text)
+
+    if (.not.is_stats) then 
+      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMtimelen)
+      error_text="Failed to find GCM time length"
+      call status_check(status,error_text)
+      allocate(GCMtime(GCMtimelen))
+    else 
+      status=nf90_inquire_dimension(gcmfid,time_dimid_gcm,len=GCMstatstimelen)
+      error_text="Failed to find GCM stats time length"
+      call status_check(status,error_text)
+      allocate(GCMstatstime(GCMstatstimelen))
+    endif
+
+    status=nf90_inq_varid(gcmfid,"Time",GCMvarid)
+    error_text="Failed to find GCM time ID"
+    call status_check(status,error_text)
+
+    ! Read GCMtime
+    if (.not.is_stats) then ! if GCM file is a diagfi, time is in sols at longitude 0°
+      status=nf90_get_var(gcmfid,GCMvarid,GCMtime)
+      error_text="Failed to load GCMtime (sols)"
+      call status_check(status,error_text)
+    else ! if GCM file is a stats, time is in LT at longitude 0°
+      status=nf90_get_var(gcmfid,GCMvarid,GCMstatstime)
+      error_text="Failed to load GCMstatstime (LT at lon 0)"
+      call status_check(status,error_text)
+    endif
+
+    ! Simulation time offset management
+    if (.not.is_stats) then
+      write(*,*) "Beginning date of the simulation file?"
+      write(*,*) "(i.e. number of sols since Ls=0 at the Time=0.0 in the GCM file)"
+      read(*,*) starttimeoffset
+      ! Add the offset to GCMtime(:)
+      GCMtime(:)=GCMtime(:)+starttimeoffset
+    endif 
+
+    ! Get the timestep of the simulation (in Martian hours)
+    if (.not.is_stats) then
+      GCMdeltat = (GCMtime(2)-GCMtime(1))*24.
+    else
+      GCMdeltat = (GCMstatstime(2)-GCMstatstime(1))
+    endif
+    write(*,*)"GCM simulation's timestep is ", GCMdeltat, " hours."
+
+    ! Check of temporal coherence between gcmfile & obsfile
+    call ls2sol(OBSLs(OBSLslen),maxsol) ! maximum date considered
+    call ls2sol(OBSLs(1),minsol) ! minimum date considered
+    if (.not.is_stats) then ! if it is a diagfi, we check the time coherence between the 2 files
+      if ((maxsol.gt.maxval(GCMtime)).or.(minsol.lt.minval(GCMtime))) then 
+        write(*,*)"Error : obsfile temporal bounds exceed the GCM simulation bounds."
+        write(*,*)"Please use a GCM file whose time interval contains the observation period."
+        stop
+      else
+        write(*,*)"Both files are temporally coherent. The program continues..."
+      endif
+
+    else ! if it is a stats, we create the array GCMtime array (in sols) covering the observation period
+         ! and filled with the mean GCM day stored in stats.nc
+
+      GCMtimelen = ((ceiling(maxsol)-floor(minsol)+1)+2) ! we add 2 days in the beginning and the end
+                                                         ! to be sure we cover the observation period
+      allocate(GCMtime(GCMstatstimelen * GCMtimelen))
+      do l=1,GCMtimelen
+        do m=1,GCMstatstimelen
+          GCMtime(m+(l-1)*GCMstatstimelen) = (floor(minsol)-1) + (l-1) + GCMstatstime(m)/24.
+        enddo
+      enddo
+      GCMtimelen = GCMstatstimelen * GCMtimelen
+      write(*,*)"GCMtime has been created from the stats.nc time and the observation period. The program continues..."
+    endif
+
+
+    ! GCM Altitude
+    status=nf90_inq_dimid(gcmfid,"altitude",alt_dimid_gcm)
+    error_text="Failed to find GCM altitude dimension"
+    call status_check(status,error_text)
+
+    status=nf90_inquire_dimension(gcmfid,alt_dimid_gcm,len=GCMaltlen)
+    error_text="Failed to find GCM altitude length"
+    call status_check(status,error_text)
+    allocate(GCMalt(GCMaltlen))
+
+    status=nf90_inq_varid(gcmfid,"altitude",GCMvarid)
+    error_text="Failed to find GCM altitude ID"
+    call status_check(status,error_text)
+
+    ! Read GCMalt
+    status=nf90_get_var(gcmfid,GCMvarid,GCMalt)
+    error_text="Failed to load GCMalt"
+    call status_check(status,error_text)
+
+    ! Check altitude attribute "units" to find out altitude type
+    status=nf90_get_att(gcmfid,GCMvarid,"units",units)
+    error_text="Failed to load GCM altitude units attribute"
+    call status_check(status,error_text)
+    if (trim(units).eq."Pa") then
+      GCMalttype='p' ! pressure coordinate
+    else if (trim(units).eq."m") then
+      GCMalttype='z' ! altitude coordinate
+    else
+      write(*,*)"I do not understand this unit ",trim(units)," for GCM altitude!"
+      if (OBSalttype.eq.'p') then
+        write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (pressure in Pa)"
+      else if (OBSalttype.eq.'z') then
+        write(*,*)"Please use zrecast to put the altitude in the same type as the MCS file (altitude in m)"
+      endif
+      stop
+    endif
+    if(OBSalttype.ne.GCMalttype) then
+      write(*,*)"Observer altitude type (", OBSalttype,") and ", &
+                "GCM altitude type (",GCMalttype,") don't match!"
+      stop
+    endif
+!===============================================================================
+! 2.3 Extraction of the wanted variables (the ones to compare with observer data)
+!===============================================================================
+    ! GCM variable to extract for both MCS dtemp & ntemp is : temp
+    GCMvarname="temp" ! temperature
+
+    ! Check that the GCM file contains that variable
+    status=nf90_inq_varid(gcmfid,trim(GCMvarname),GCMvarid)
+    error_text="Failed to find variable "//trim(GCMvarname)//" in "//trim(gcmfile)
+    call status_check(status,error_text)
+
+    ! Sanity checks on the variable
+    status=nf90_inquire_variable(gcmfid,GCMvarid,ndims=nb,dimids=GCMvarshape)
+    error_text="Failed to obtain information on variable "//trim(GCMvarname)
+    call status_check(status,error_text)
+
+    ! Check that it is a 4D variable
+    if (nb.ne.4) then
+      write(*,*) "Error, expected a 4D (lon-lat-alt-time) variable for ", trim(GCMvarname)
+      stop
+    endif
+    ! Check that its dimensions are indeed lon,lat,alt,time (in the right order)
+    if (GCMvarshape(1).ne.lon_dimid_gcm) then
+      write(*,*) "Error, expected first dimension of ", trim(GCMvarname), " to be longitude!"
+      stop
+    endif
+    if (GCMvarshape(2).ne.lat_dimid_gcm) then
+      write(*,*) "Error, expected second dimension of ", trim(GCMvarname), " to be latitude!"
+      stop
+    endif
+    if (GCMvarshape(3).ne.alt_dimid_gcm) then
+      write(*,*) "Error, expected third dimension of ", trim(GCMvarname), " to be altitude!"
+      stop
+    endif
+    if (GCMvarshape(4).ne.time_dimid_gcm) then
+      write(*,*) "Error, expected fourth dimension of ", trim(GCMvarname), " to be time!"
+      stop
+    endif
+
+    ! Length allocation for each dimension of the 4D variable
+    allocate(GCM_var(GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen))
+
+    ! Load datasets
+    if (.not.is_stats) then
+      status=nf90_get_var(gcmfid,GCMvarid,GCM_var)
+      error_text="Failed to load "//trim(GCMvarname)
+      call status_check(status,error_text)
+    else
+      ! if it is a stats file, we load only the first sol, and then copy it to all the other sols
+      status=nf90_get_var(gcmfid,GCMvarid,GCM_var(:,:,:,1:GCMstatstimelen))
+      error_text="Failed to load "//trim(GCMvarname)
+      call status_check(status,error_text)
+    !  write(*,*)"GCMstatstimelen = ", GCMstatstimelen
+    !  write(*,*)"GCMtimelen = ", GCMtimelen
+      do l=(GCMstatstimelen+1),GCMtimelen
+        if (modulo(l,GCMstatstimelen).ne.0) then
+          GCM_var(:,:,:,l) = GCM_var(:,:,:,modulo(l,GCMstatstimelen))
+        else ! if l is a multiple of GCMstatstimelen, since index modulo(l,GCMstatstimelen)=0 doesn't exist,
+             ! we make a special case
+          GCM_var(:,:,:,l) = GCM_var(:,:,:,GCMstatstimelen)
+        endif
+      enddo
+    endif
+    write(*,*) "Loaded ",trim(GCMvarname), " from the GCM file"
+
+    ! Get dataset's missing_value attribute
+    status=nf90_get_att(gcmfid,GCMvarid,"missing_value",GCMmiss_val)
+    error_text="Failed to load missing_value attribute"
+    call status_check(status,error_text)
+    write(*,'("  with missing_value attribute : ",1pe12.5)')GCMmiss_val
+
+  ENDIF ! end of IF (OBSvarname.EQ."dtemp")
+  
+!===================================================================================================
+! 3. Output file
+!===================================================================================================
+!===============================================================================
+! 3.1 Create the output file & initialize the coordinates
+!===============================================================================
+  IF (OBSvarname.EQ."dtemp") THEN ! creating the output file is done only in the first loop
+    ! Name of the outfile
+    if (.not.is_stats) then
+      outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMdiagfi.nc"
+    else
+      outfile=obsfile(1:index(obsfile, ".nc")-1)//"_GCMstats.nc"
+    endif
+
+    ! Creation of the outfile
+    status=nf90_create(outfile,nf90_clobber,outfid) ! NB: clobber mode=overwrite any dataset with the same file name !
+    error_text="Error: could not create file "//trim(outfile)
+    call status_check(status,error_text)
+    write(*,*)"";write(*,*)"-> Output file is: ",trim(outfile)
+    
+    ! Creation of the dimensions
+    call inidim(outfid,OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen,OBSlon,OBSlat,OBSalt,OBSLs,OBSalttype,&
+                lon_dimid_obs,lat_dimid_obs,alt_dimid_obs,time_dimid_obs)
+
+    !write(*,*)"Dimensions ID in the outfile are :"           
+    !write(*,*)"lon_dimid=",lon_dimid_obs
+    !write(*,*)"lat_dimid=",lat_dimid_obs
+    !write(*,*)"alt_dimid=",alt_dimid_obs
+    !write(*,*)"time_dimid=",time_dimid_obs
+    
+  ENDIF ! end of IF (OBSvarname.EQ."dtemp")
+  
+!===============================================================================
+! 3.2 Create the outfile variables (other variables than the coordinates)
+!===============================================================================
+  ! Switch to netcdf define mode
+  status=nf90_redef(outfid)
+  error_text="Error: could not switch to define mode in the outfile"
+  call status_check(status,error_text)
+
+  ! Insert the definition of the variables
+  status=nf90_def_var(outfid,trim(OBSvarname),nf90_float,OBSvarshape,outvarid)
+  error_text="Error: could not define the variable "//trim(OBSvarname)//" in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_def_var(outfid,trim(OBSLTave_name),nf90_float,LTshape,LT_id)
+  error_text="Error: could not define the variable "//trim(OBSLTave_name)//" in the outfile"
+  call status_check(status,error_text)
+
+  ! Write the attributes
+  !! Temperature
+  SELECT CASE (OBSvarname)
+  CASE ("dtemp")
+    status=nf90_put_att(outfid,outvarid,"long_name","Temperature - day side (interpolated from GCM)")
+  CASE ("ntemp")
+    status=nf90_put_att(outfid,outvarid,"long_name","Temperature - night side (interpolated from GCM)")
+  END SELECT
+  error_text="Error: could not put the long_name attribute for the variable "//trim(OBSvarname)//" in the outfile"
+  call status_check(status,error_text)
+  
+  status=nf90_put_att(outfid,outvarid,"units","K")
+  error_text="Error: could not put the units attribute for the variable "//trim(OBSvarname)//" in the outfile"
+  call status_check(status,error_text)
+  
+  status=nf90_put_att(outfid,outvarid,"_FillValue",OBSmiss_val)
+  error_text="Error: could not put the _FillValue attribute for the variable "//trim(OBSvarname)//" in the outfile"
+  call status_check(status,error_text)
+
+  write(*,*)trim(OBSvarname)," has been created in the outfile"
+  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
+
+  !! Local Time (average in bin)
+  SELECT CASE (OBSvarname)
+  CASE ("dtemp")
+    status=nf90_put_att(outfid,LT_id,"long_name","Local Time, evaluated as average local time in bin - day side [6h, 18h]")
+  CASE ("ntemp")
+    status=nf90_put_att(outfid,LT_id,"long_name","Local Time, evaluated as average local time in bin - night side [18h, 6h]")
+  END SELECT
+  error_text="Error: could not put the long_name attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
+  call status_check(status,error_text)
+  
+  status=nf90_put_att(outfid,LT_id,"units","hours")
+  error_text="Error: could not put the units attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
+  call status_check(status,error_text)
+  
+  status=nf90_put_att(outfid,LT_id,"_FillValue",LTmiss_val)
+  error_text="Error: could not put the _FillValue attribute for the variable "//trim(OBSLTave_name)//" in the outfile"
+  call status_check(status,error_text)
+
+  write(*,*)"Local Time (", trim(OBSLTave_name), ") has been created in the outfile"
+  write(*,'("  with missing_value attribute : ",1pe12.5)')OBSmiss_val
+
+  ! End the netcdf define mode (and thus enter the "data writing" mode)
+  status=nf90_enddef(outfid)
+  error_text="Error: could not close the define mode of the outfile"
+  call status_check(status,error_text)
+
+!===============================================================================
+! 3.3 Fill the output file
+!===============================================================================
+  allocate(outvar(OBSlonlen,OBSlatlen,OBSaltlen,OBSLslen))
+
+!!===============================================================
+!! 3.3.1 Do some checks
+!!===============================================================
+  Ls: do l=1,OBSLslen ! loop on all observed Solar longitudes
+    Ls_val=OBSLs(l) ! Ls_val=center of the output bin
+    if ((Ls_val.lt.0.).or.(Ls_val.gt.360.)) then
+      write(*,*) "Unexpected value for OBSLs: ",Ls_val
+      stop
+    endif
+
+    ! Convert the Ls bin into a sol interval on which the binning is done :
+    !-> get the index of the maximum value of GCM sol (m_maxsol) that is lower than Ls bin's superior bound (maxsol)
+    call ls2sol(Ls_val+OBSdeltaLs/2.,maxsol)
+    m_maxsol=1
+    do while (((GCMtime(m_maxsol)+0.5).lt.maxsol).AND.(m_maxsol.le.(GCMtimelen-1)))
+  !! The +0.5 is there to take into account the whole planet (lon=[-180°;180°]) and not just the lon=0° from the GCM
+      m_maxsol=m_maxsol+1
+    enddo
+    !-> get the index of the minimum value of GCM sol (m_minsol) that is greater than Ls bin's inferior bound (minsol)
+    call ls2sol(Ls_val-OBSdeltaLs/2.,minsol)
+    m_minsol=1
+    do while (((GCMtime(m_minsol)-0.5).le.minsol).AND.(m_minsol.le.(GCMtimelen-1)))
+  !! Same comment for the -0.5
+      m_minsol=m_minsol+1
+    enddo
+    if (m_minsol.gt.m_maxsol) then
+      write(*,*) "No value in gcmfile between sol=",minsol," and sol=",maxsol," (Ls=",Ls_val,"°)"
+      ! Write a missing_value to output
+      outvar(:,:,:,l)=OBSmiss_val
+      CYCLE Ls ! go directly to the next Ls
+    endif  
+    ! Get all the integer values of GCM sols that fit in this interval
+    solcount=floor(GCMtime(m_maxsol))-ceiling(GCMtime(m_minsol))+1 ! sols that are not fully in the interval are not counted
+    allocate(int_sol_list(solcount))
+    m_int=1
+  !  write(*,*) "GCMminsol=", GCMtime(m_minsol)
+  !  write(*,*)"GCMmaxsol=", GCMtime(m_maxsol)
+    do m=m_minsol,m_maxsol ! loop on all GCM sols of the bin
+      sol=GCMtime(m)
+      if (sol.eq.floor(sol)) then
+  !      write(*,*)"sol=",sol
+        int_sol_list(m_int)=sol
+        m_int=m_int+1
+      endif
+    enddo
+  !  write(*,*)"int_sol_list=",int_sol_list
+
+    latitude: do j=1,OBSlatlen ! loop on all observed latitudes
+      lat_val=OBSlat(j)
+      if ((lat_val.lt.-90.).or.(lat_val.gt.90.)) then
+        write(*,*) "Unexpected value for OBSlat: ",lat_val
+        stop
+      endif
+
+      longitude: do i=1,OBSlonlen ! loop on all observed longitudes
+        lon_val=OBSlon(i)
+        if ((lon_val.lt.-360.).or.(lon_val.gt.360.)) then
+          write(*,*) "Unexpected value for lon_val: ",lon_val
+          stop
+        endif
+        ! We want lon_val in [-180:180] for the subroutine extraction
+        if (lon_val.lt.-180.) lon_val=lon_val+360.
+        if (lon_val.gt.180.) lon_val=lon_val-360.
+
+        LT_val=OBSLT(i,j,l) ! find the Observer corresponding LT value at lon_val, lat_val, Ls_val
+        
+        if ((LT_val.lt.0.).or.(LT_val.gt.24.)) then
+          if (LT_val.eq.LTmiss_val) then
+  !            write(*,*) "Missing value in obsfile for LT_val"
+            ! Write a missing_value to output
+            outvar(i,j,:,l)=OBSmiss_val
+            CYCLE longitude ! go directly to the next longitude
+          else
+            write(*,*) "Unexpected value for LT_val: ",LT_val
+            stop
+          endif
+        endif
+        ! Generate the sol list for the interpolation
+        allocate(sol_list(solcount*LTcount))
+        call gen_sol_list(solcount,int_sol_list,LTcount,LT_val,OBSLTmax(i,j,l),OBSLTmin(i,j,l),lon_val,LTmod,OBSvarname,&
+                          sol_list)
+
+        altitude: do k=1,OBSaltlen ! loop on all observed altitudes
+          alt_val=OBSalt(k)
+          if (OBS_var(i,j,k,l).eq.OBSmiss_val) then
+  !          write(*,*) "Missing value in obsfile for ",OBSvarname
+            ! Write a missing_value to output
+            outvar(i,j,k,l)=OBSmiss_val
+            CYCLE altitude ! go directly to the next altitude
+          endif
+
+!!===============================================================
+!! 3.3.2 Compute GCM sol date corresponding to Observer Ls 
+!!       (via m_(min/max)sol) and LT (via OBSLT(min/max))
+!!===============================================================       
+          solerrcount=0
+          solbinned_value=0
+          sol_bin: do n=1,solcount*LTcount ! loop on all GCM sols of the bin
+            sol=sol_list(n)
+  !          write(*,*)"sol=",sol
+!!===============================================================  
+!! 3.3.3 Do the interpolation and binning for the given location
+!!===============================================================
+            call extraction(lon_val,lat_val,alt_val,sol,&
+                            GCMlonlen,GCMlatlen,GCMaltlen,GCMtimelen,&
+                            GCMlon,GCMlat,GCMalt,GCMtime,&
+                            GCM_var,GCMmiss_val,GCMalttype,GCMvarname,extr_value)
+
+            if (extr_value.eq.GCMmiss_val) then
+  !            write(*,*) "Missing value in gcmfile at lon=",lon_val,"; lat=",lat_val,"; alt=",alt_val,"; sol=",sol 
+              solerrcount=solerrcount+1  
+              CYCLE sol_bin ! go directly to the next GCM sol of the bin
+            endif
+            solbinned_value=solbinned_value+extr_value
+
+          enddo sol_bin ! end loop on all GCM sols of the bin
+          if ((solcount*LTcount-solerrcount).ne.0) then
+            solbinned_value=solbinned_value/(solcount*LTcount-solerrcount)
+          else
+  !          write(*,*)"No GCM value in this sol bin"
+            solbinned_value=OBSmiss_val
+          endif
+          ! Write value to output
+          outvar(i,j,k,l)=solbinned_value
+
+          errcount=errcount+solerrcount
+
+        enddo altitude ! end loop on observed altitudes
+
+        deallocate(sol_list)
+
+      enddo longitude ! end loop on observed longitudes
+    enddo latitude ! end loop on observed latitudes
+
+    deallocate(int_sol_list)
+
+  enddo Ls ! end loop on observed Solar longitudes
+
+  !write(*,*)"Nb of GCM missing values :",errcount
+
+!!===============================================================  
+!! 3.3.4 Write the data in the netcdf output file
+!!===============================================================
+  status = nf90_put_var(outfid, outvarid, outvar)
+  error_text="Error: could not write "//trim(OBSvarname)//" data in the outfile"
+  call status_check(status,error_text)
+
+  status = nf90_put_var(outfid, LT_id, OBSLT) ! write the MCS d/ntimeave as the output Local Time
+  error_text="Error: could not write Local Time data in the outfile"
+  call status_check(status,error_text)
+
+!!==============================================================================  
+!! 4. End of the Day/Night loop
+!!==============================================================================
+  IF (OBSvarname.EQ."dtemp") THEN
+    ! this is the end of the first loop (on daytime values)
+    ! and we still have to loop on nighttime values
+    OBSvarname="ntemp"
+    deallocate(OBS_var)
+    deallocate(OBSLT)
+    deallocate(OBSLTmax)
+    deallocate(OBSLTmin)
+    deallocate(outvar)
+    write(*,*)"";write(*,*)""; write(*,*) "Beginning the 2nd loop, on nighttime values" ; write(*,*)""
+    CYCLE DAY_OR_NIGHT
+  ELSE ! i.e. OBSvarname="ntemp"
+    ! this is the end of the second loop (on nighttime values) 
+    ! and thus the end of the program
+    EXIT DAY_OR_NIGHT
+  ENDIF
+ENDDO DAY_OR_NIGHT ! end of the day/night loop that begins in section 1.3
+
+!=============================================================================== 
+! 5. Close output file
+!===============================================================================
+status = nf90_close(outfid)
+error_text="Error: could not close file "//trim(outfile)
+call status_check(status,error_text)
+
+write(*,*)"";write(*,*)"-> Program completed!"
+
+end program simu_MCS_temp
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!===================================================================================================
+! Subroutines
+!===================================================================================================
+
+subroutine extraction(lon,lat,alt,sol,&
+                  lonlen,latlen,altlen,timelen,&
+                  longitude,latitude,altitude,time,&
+                  field,missing_value,alttype,varname,value)
+
+implicit none
+!================================================================
+! Arguments:
+!================================================================
+real,intent(in) :: lon  ! sought longitude (deg, in [-180:180])
+real,intent(in) :: lat  ! sought latitude (deg, in [-90:90])
+real,intent(in) :: alt  ! sought altitude (m or Pa)
+real,intent(in) :: sol  ! sought date (sols)
+integer,intent(in) :: lonlen
+integer,intent(in) :: latlen
+integer,intent(in) :: altlen
+integer,intent(in) :: timelen
+real,intent(in) :: longitude(lonlen)
+real,intent(in) :: latitude(latlen)
+real,intent(in) :: altitude(altlen)
+real,intent(in) :: time(timelen)
+real,intent(in) :: field(lonlen,latlen,altlen,timelen)
+real,intent(in) :: missing_value ! default value for "no data"
+character,intent(in) :: alttype ! altitude coord. type:'z' (altitude, m) 
+                                !                      'p' (pressure, Pa)
+character(len=*),intent(in) :: varname ! variable name (in GCM file)
+real,intent(out) :: value
+
+!================================================================
+! Local variables:
+!================================================================
+real,save :: prev_lon=-99 ! previous value of 'lon' routine was called with
+real,save :: prev_lat=-99 ! previous value of 'lat' routine was called with
+real,save :: prev_alt=-9.e20 ! ! previous value of 'alt'
+real,save :: prev_sol=-99 ! previous value of 'sol' routine was called with
+
+! encompasing indexes:
+integer,save :: ilon_inf=-1,ilon_sup=-1 ! along longitude
+integer,save :: ilat_inf=-1,ilat_sup=-1 ! along latitude
+integer,save :: ialt_inf=-1,ialt_sup=-1 ! along altitude
+integer,save :: itim_inf=-1,itim_sup=-1 ! along time
+
+! intermediate interpolated values
+real :: t_interp(2,2,2) ! after time interpolation
+real :: zt_interp(2,2) ! after altitude interpolation
+real :: yzt_interp(2) ! after latitude interpolation
+real :: coeff ! interpolation coefficient
+
+integer :: i
+
+! By default, set value to missing_value
+value=missing_value
+
+!================================================================
+! 1. Find encompassing indexes
+!================================================================
+if (lon.ne.prev_lon) then
+  do i=1,lonlen-1
+    if (longitude(i).le.lon) then
+      ilon_inf=i
+    else
+      exit
+    endif
+  enddo
+  ilon_sup=ilon_inf+1
+endif ! of if (lon.ne.prev_lon)
+!write(*,*) 'ilon_inf=',ilon_inf," longitude(ilon_inf)=",longitude(ilon_inf)
+
+if (lat.ne.prev_lat) then
+  ! recall that latitudes start from north pole
+  do i=1,latlen-1
+    if (latitude(i).ge.lat) then
+      ilat_inf=i
+    else
+      exit
+    endif
+  enddo
+  ilat_sup=ilat_inf+1
+endif ! of if (lat.ne.prev_lat)
+!write(*,*) 'ilat_inf=',ilat_inf," latitude(ilat_inf)=",latitude(ilat_inf)
+
+if (alt.ne.prev_alt) then
+  if (alttype.eq.'p') then ! pressures are ordered from max to min
+    !handle special case for alt not in the altitude(1:altlen) interval
+    if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen))) then
+      ialt_inf=-1
+      ialt_sup=-1
+      ! return to main program (with value=missing_value)
+!      write(*,*)"Problem in extraction : GCM alt p"
+      return
+    else ! general case
+      do i=1,altlen-1
+        if (altitude(i).ge.alt) then
+          ialt_inf=i
+        else
+          exit
+        endif
+      enddo
+      ialt_sup=ialt_inf+1
+    endif ! of if ((alt.gt.altitude(1)).or.(alt.lt.altitude(altlen)))
+  else ! altitudes (m) are ordered from min to max
+    !handle special case for alt not in the altitude(1:altlen) interval
+    if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen))) then
+      ialt_inf=-1
+      ialt_sup=-1
+      ! return to main program (with value=missing_value)
+!      write(*,*)"Problem in extraction : GCM alt z"
+      return
+    else ! general case
+      do i=1,altlen-1
+        if (altitude(i).le.alt) then
+          ialt_inf=i
+        else
+          exit
+        endif
+      enddo
+      ialt_sup=ialt_inf+1
+    endif ! of if ((alt.lt.altitude(1)).or.(alt.gt.altitude(altlen)))
+  endif ! of if (alttype.eq.'p') 
+endif ! of if (alt.ne.prev_alt)
+!write(*,*) 'ialt_inf=',ialt_inf," altitude(ialt_inf)=",altitude(ialt_inf)
+
+if (sol.ne.prev_sol) then
+  !handle special case for sol not in the time(1):time(timenlen) interval
+  if ((sol.lt.time(1)).or.(sol.gt.time(timelen))) then
+    itim_inf=-1
+    itim_sup=-1
+    ! return to main program (with value=missing_value)
+!    write(*,*)"Problem in extraction : GCM sol"
+    return
+  else ! general case
+    do i=1,timelen-1
+      if (time(i).le.sol) then
+        itim_inf=i
+      else
+        exit
+      endif
+    enddo
+    itim_sup=itim_inf+1
+  endif ! of if ((sol.lt.time(1)).or.(sol.gt.time(timenlen)))
+endif ! of if (sol.ne.prev_sol)
+!write(*,*) 'itim_inf=',itim_inf," time(itim_inf)=",time(itim_inf)
+!write(*,*) 'itim_sup=',itim_sup," time(itim_sup)=",time(itim_sup)
+
+!================================================================
+! 2. Interpolate
+!================================================================
+! check that there are no "missing_value" in the field() elements we need
+! otherwise return to main program (with value=missing_value)
+if (field(ilon_inf,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 1 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 2 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 3 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 4 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 5 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 6 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 7 in extraction"
+  return
+ endif
+if (field(ilon_inf,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 8 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_inf,ialt_inf,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 9 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_inf,ialt_inf,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 10 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_inf,ialt_sup,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 11 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_inf,ialt_sup,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 12 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_sup,ialt_inf,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 13 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_sup,ialt_inf,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 14 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_sup,ialt_sup,itim_inf).eq.missing_value) then
+!  write(*,*)"Error 15 in extraction"
+  return
+ endif
+if (field(ilon_sup,ilat_sup,ialt_sup,itim_sup).eq.missing_value) then
+!  write(*,*)"Error 16 in extraction"
+  return
+ endif
+
+!================================================================
+! 2.1 Linear interpolation in time
+!================================================================
+coeff=(sol-time(itim_inf))/(time(itim_sup)-time(itim_inf))
+t_interp(1,1,1)=field(ilon_inf,ilat_inf,ialt_inf,itim_inf)+ &
+                coeff*(field(ilon_inf,ilat_inf,ialt_inf,itim_sup)- &
+                       field(ilon_inf,ilat_inf,ialt_inf,itim_inf))
+t_interp(1,1,2)=field(ilon_inf,ilat_inf,ialt_sup,itim_inf)+ &
+                coeff*(field(ilon_inf,ilat_inf,ialt_sup,itim_sup)- &
+                       field(ilon_inf,ilat_inf,ialt_sup,itim_inf))
+t_interp(1,2,1)=field(ilon_inf,ilat_sup,ialt_inf,itim_inf)+ &
+                coeff*(field(ilon_inf,ilat_sup,ialt_inf,itim_sup)- &
+                       field(ilon_inf,ilat_sup,ialt_inf,itim_inf))
+t_interp(1,2,2)=field(ilon_inf,ilat_sup,ialt_sup,itim_inf)+ &
+                coeff*(field(ilon_inf,ilat_sup,ialt_sup,itim_sup)- &
+                       field(ilon_inf,ilat_sup,ialt_sup,itim_inf))
+t_interp(2,1,1)=field(ilon_sup,ilat_inf,ialt_inf,itim_inf)+ &
+                coeff*(field(ilon_sup,ilat_inf,ialt_inf,itim_sup)- &
+                       field(ilon_sup,ilat_inf,ialt_inf,itim_inf))
+t_interp(2,1,2)=field(ilon_sup,ilat_inf,ialt_sup,itim_inf)+ &
+                coeff*(field(ilon_sup,ilat_inf,ialt_sup,itim_sup)- &
+                       field(ilon_sup,ilat_inf,ialt_sup,itim_inf))
+t_interp(2,2,1)=field(ilon_sup,ilat_sup,ialt_inf,itim_inf)+ &
+                coeff*(field(ilon_sup,ilat_sup,ialt_inf,itim_sup)- &
+                       field(ilon_sup,ilat_sup,ialt_inf,itim_inf))
+t_interp(2,2,2)=field(ilon_sup,ilat_sup,ialt_sup,itim_inf)+ &
+                coeff*(field(ilon_sup,ilat_sup,ialt_sup,itim_sup)- &
+                       field(ilon_sup,ilat_sup,ialt_sup,itim_inf))
+
+!================================================================
+! 2.2 Vertical interpolation
+!================================================================
+if (((varname=='rho').or.(varname=='pressure')).and.(alttype=='z')) then
+  ! do the interpolation on the log of the quantity
+  coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
+  zt_interp(1,1)=log(t_interp(1,1,1))+coeff* &
+                             (log(t_interp(1,1,2))-log(t_interp(1,1,1)))
+  zt_interp(1,2)=log(t_interp(1,2,1))+coeff* &
+                             (log(t_interp(1,2,2))-log(t_interp(1,2,1)))
+  zt_interp(2,1)=log(t_interp(2,1,1))+coeff* &
+                             (log(t_interp(2,1,2))-log(t_interp(2,1,1)))
+  zt_interp(2,2)=log(t_interp(2,2,1))+coeff* &
+                             (log(t_interp(2,2,2))-log(t_interp(2,2,1)))
+  zt_interp(1:2,1:2)=exp(zt_interp(1:2,1:2))
+else ! general case
+  coeff=(alt-altitude(ialt_inf))/(altitude(ialt_sup)-altitude(ialt_inf))
+  zt_interp(1,1)=t_interp(1,1,1)+coeff*(t_interp(1,1,2)-t_interp(1,1,1))
+  zt_interp(1,2)=t_interp(1,2,1)+coeff*(t_interp(1,2,2)-t_interp(1,2,1))
+  zt_interp(2,1)=t_interp(2,1,1)+coeff*(t_interp(2,1,2)-t_interp(2,1,1))
+  zt_interp(2,2)=t_interp(2,2,1)+coeff*(t_interp(2,2,2)-t_interp(2,2,1))
+endif
+
+!================================================================
+! 2.3 Latitudinal interpolation
+!================================================================
+coeff=(lat-latitude(ilat_inf))/(latitude(ilat_sup)-latitude(ilat_inf))
+yzt_interp(1)=zt_interp(1,1)+coeff*(zt_interp(1,2)-zt_interp(1,1))
+yzt_interp(2)=zt_interp(2,1)+coeff*(zt_interp(2,2)-zt_interp(2,1))
+
+!================================================================
+! 2.4 longitudinal interpolation
+!================================================================
+coeff=(lon-longitude(ilon_inf))/(longitude(ilon_sup)-longitude(ilon_inf))
+value=yzt_interp(1)+coeff*(yzt_interp(2)-yzt_interp(1))
+
+end subroutine extraction
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine inidim(outfid,lonlen,latlen,altlen,timelen,lon,lat,alt,time,units_alt,&
+                     londimid,latdimid,altdimid,timedimid)
+!================================================================
+! Purpose:
+! Initialize a data file (NetCDF format)
+!================================================================
+! Remarks:
+! The NetCDF file remains open
+!================================================================
+use netcdf ! NetCDF definitions
+implicit none
+!================================================================
+! Arguments:
+!================================================================
+integer, intent(in):: outfid
+! outfid: [netcdf] file ID
+integer, intent(in):: lonlen
+! lonlen: longitude length (# of longitude values)
+integer, intent(in):: latlen
+! latlen: latitude length (# of latitude values)
+integer, intent(in):: altlen
+! altlen: altitude length (# of altitude values)
+integer, intent(in):: timelen
+! timelen: time length (# of time values)
+real, intent(in):: lon(lonlen)
+! lon(): longitude
+real, intent(in):: lat(latlen)
+! lat(): latitude
+real, intent(in):: alt(altlen)
+! alt(): altitude
+real, intent(in):: time(timelen)
+! time(): time (Ls)
+character(len=1), intent(in) :: units_alt
+! units_alt: altitude coord. type:'z' (altitude, m) 'p' (pressure, Pa)
+integer,intent(inout):: londimid
+! londimid: [netcdf] lon() (i.e.: longitude) ID in MCS and output files (they are the same)
+integer,intent(inout) :: latdimid
+! latdimid: [netcdf] lat() ID in MCS and output files (they are the same)
+integer,intent(inout):: altdimid
+! altdimid: [netcdf] alt() ID in MCS and output files (they are the same)
+integer,intent(inout):: timedimid
+! timedimid: [netcdf] time() ID in MCS and output files (they are the same)
+
+!================================================================
+! Local variables:
+!================================================================
+integer :: lonvarid,latvarid,altvarid,timevarid,status
+! *varid: [netcdf] ID of a variable
+! status: [netcdf]  return error code (from called subroutines)
+character(len=256) :: error_text
+integer :: d ! loop index on dimensions ID
+
+!===============================================================
+! 1. Define/write "dimensions" in the same order than MCS file
+!    and get their IDs
+!================================================================
+do d=1,4
+  if (altdimid.eq.d) then
+    status=nf90_def_dim(outfid, "altitude", altlen, altdimid)
+    error_text="Error: could not define the altitude dimension in the outfile"
+    call status_check(status,error_text)
+  else if (timedimid.eq.d) then
+    status=nf90_def_dim(outfid, "time", timelen, timedimid)
+    error_text="Error: could not define the time dimension in the outfile"
+    call status_check(status,error_text)
+  else if (latdimid.eq.d) then
+    status=nf90_def_dim(outfid, "latitude", latlen, latdimid)
+    error_text="Error: could not define the latitude dimension in the outfile"
+    call status_check(status,error_text)
+  else if (londimid.eq.d) then
+    status=nf90_def_dim(outfid, "longitude", lonlen, londimid)
+    error_text="Error: could not define the longitude dimension in the outfile"
+    call status_check(status,error_text)
+  endif
+enddo
+
+!================================================================
+! 2. Define "variables" and their attributes
+!================================================================
+!================================================================
+! 2.1 Write "longitude" (data and attributes)
+!================================================================
+
+! Insert the definition of the variable
+status=nf90_def_var(outfid,"longitude",nf90_float,(/londimid/),lonvarid)
+error_text="Error: could not define the longitude variable in the outfile"
+call status_check(status,error_text)
+
+! Write the attributes
+status=nf90_put_att(outfid,lonvarid,"long_name","longitude")
+error_text="Error: could not put the long_name attribute for the longitude variable in the outfile"
+call status_check(status,error_text)
+
+status=nf90_put_att(outfid,lonvarid,"units","degrees_east")
+error_text="Error: could not put the units attribute for the longitude variable in the outfile"
+call status_check(status,error_text)
+
+!================================================================
+! 2.2 "latitude"
+!================================================================
+
+! Insert the definition of the variable
+status=nf90_def_var(outfid,"latitude",nf90_float,(/latdimid/),latvarid)
+error_text="Error: could not define the latitude variable in the outfile"
+call status_check(status,error_text)
+
+! Write the attributes
+status=nf90_put_att(outfid,latvarid,"long_name","latitude")
+error_text="Error: could not put the long_name attribute for the latitude variable in the outfile"
+call status_check(status,error_text)
+
+status=nf90_put_att(outfid,latvarid,"units","degrees_north")
+error_text="Error: could not put the units attribute for the latitude variable in the outfile"
+call status_check(status,error_text)
+
+!================================================================
+! 2.3 Write "altitude" (data and attributes)
+!================================================================
+
+! Insert the definition of the variable
+status=nf90_def_var(outfid,"altitude",nf90_float,(/altdimid/),altvarid)
+error_text="Error: could not define the altitude variable in the outfile"
+call status_check(status,error_text)
+
+! Write the attributes
+if (units_alt.eq.'p') then ! pressure coordinate
+  status=nf90_put_att(outfid,altvarid,"long_name","pressure")
+  error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_put_att(outfid,altvarid,'units',"Pa")
+  error_text="Error: could not put the units attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_put_att(outfid,altvarid,'positive',"down")
+  error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_put_att(outfid,altvarid,'comment',&
+  "The vertical variable is in fact pressure, not altitude. We just call it altitude for easy reading in Ferret and Grads.")
+  error_text="Error: could not put the comment attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+else if (units_alt.eq.'z') then ! altitude coordinate
+  status=nf90_put_att(outfid,altvarid,"long_name","altitude")
+  error_text="Error: could not put the long_name attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_put_att(outfid,altvarid,'units',"m")
+  error_text="Error: could not put the units attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+  status=nf90_put_att(outfid,altvarid,'positive',"up")
+  error_text="Error: could not put the positive attribute for the altitude variable in the outfile"
+  call status_check(status,error_text)
+
+else
+  write(*,*)"I do not understand this unit type ",units_alt," for outfile altitude!"
+  stop
+end if
+
+!================================================================
+! 2.4 "time"
+!================================================================
+
+! Insert the definition of the variable
+status=nf90_def_var(outfid,"time",nf90_float,(/timedimid/),timevarid)
+error_text="Error: could not define the time variable in the outfile"
+call status_check(status,error_text)
+
+! Write the attributes
+status=nf90_put_att(outfid,timevarid,"long_name","Solar longitude")
+error_text="Error: could not put the long_name attribute for the time variable in the outfile"
+call status_check(status,error_text)
+
+status=nf90_put_att(outfid,timevarid,"units","days since 0000-01-1 00:00:00")
+error_text="Error: could not put the units attribute for the time variable in the outfile"
+call status_check(status,error_text)
+
+status=nf90_put_att(outfid,timevarid,"comment",&
+"Units is in fact degrees, but set to a dummy value of days for compatibility with Ferret and Grads.")
+error_text="Error: could not put the comment attribute for the time variable in the outfile"
+call status_check(status,error_text)
+
+!================================================================
+! 2.5 End netcdf define mode
+!================================================================
+status=nf90_enddef(outfid)
+error_text="Error: could not end the define mode of the outfile"
+call status_check(status,error_text)
+
+!================================================================
+! 3. Write "variables" (data of the dimension variables)
+!================================================================
+! "time"
+status=nf90_put_var(outfid,timevarid,time)
+error_text="Error: could not write the time variable's data in the outfile"
+call status_check(status,error_text)
+
+! "latitude"
+status=nf90_put_var(outfid,latvarid,lat)
+error_text="Error: could not write the latitude variable's data in the outfile"
+call status_check(status,error_text)
+
+! "longitude"
+status=nf90_put_var(outfid,lonvarid,lon)
+error_text="Error: could not write the longitude variable's data in the outfile"
+call status_check(status,error_text)
+
+! "altitude"
+status=nf90_put_var(outfid,altvarid,alt)
+error_text="Error: could not write the altitude variable's data in the outfile"
+call status_check(status,error_text)
+
+end subroutine inidim
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine ls2sol(ls,sol)
+
+implicit none
+!================================================================
+!  Arguments:
+!================================================================
+real,intent(in) :: ls
+real,intent(out) :: sol
+
+!================================================================
+!  Local:
+!================================================================
+double precision xref,zx0,zteta,zz
+!xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
+double precision year_day 
+double precision peri_day,timeperi,e_elips
+double precision pi,degrad 
+parameter (year_day=668.6d0) ! number of sols in a martian year
+parameter (peri_day=485.35d0) ! date (in sols) of perihelion
+!timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
+parameter (timeperi=1.90258341759902d0)
+parameter (e_elips=0.0934d0)  ! eccentricity of orbit
+parameter (pi=3.14159265358979d0)
+parameter (degrad=57.2957795130823d0)
+
+      if (abs(ls).lt.1.0e-5) then
+         if (ls.ge.0.0) then
+            sol = 0.0
+         else
+            sol = real(year_day)
+         end if
+         return
+      end if
+
+      zteta = ls/degrad + timeperi
+      zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
+      xref = zx0-e_elips*dsin(zx0)
+      zz = xref/(2.*pi)
+      sol = real(zz*year_day + peri_day)
+      if (sol.lt.0.0) sol = sol + real(year_day)
+      if (sol.ge.year_day) sol = sol - real(year_day)
+
+
+end subroutine ls2sol
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine gen_sol_list(solcount,int_sol_list,LTcount,LTave,LTmax,LTmin,lon,f_LT,LTname,&
+                        sol_list)
+!================================================================
+! Purpose:
+! Generate a list that is the combination of two lists :
+! - int_sol_list, which is a list of integer values of sol
+! - LT_list, which contains a number LTcount of LT values in the
+!   interval [LTmin;LTmax] which are evenly distributed around
+!   LTave (so that their average is LTave)
+!================================================================
+
+implicit none
+!================================================================
+! Arguments:
+!================================================================
+integer, intent(in) :: solcount ! nb of integer values of sol
+real, intent(in) :: int_sol_list(solcount) ! list of the integer values of sol
+integer, intent(in) :: LTcount ! nb of LT per sol to be interpolated
+real, intent(in) :: LTave ! average of LT
+real, intent(in) :: LTmax, LTmin ! bounds of the LT interval for the interpolation
+real, intent(in) :: lon ! longitude value for the transformation into a sol value at lon=0°
+external f_LT ! function that changes LT interval for night LT
+real :: f_LT
+character (len=64), intent(in) :: LTname ! to know if we have day or night values
+
+real, intent(out) :: sol_list(solcount*LTcount) ! all the sol values at lon=0° to interpolate
+
+!================================================================
+! Local variables:
+!================================================================
+integer :: N
+integer :: a,b,c ! loop iteration indexes
+real :: LT_list(LTcount)
+
+N = floor(LTcount/2.)
+
+!================================================================
+! 1. Filling of LT_list
+!================================================================
+SELECT CASE (LTname)
+CASE ("dtemp")
+  if (abs(LTave-LTmax).le.abs(LTave-LTmin)) then ! LTave is closer to LTmax than to LTmin
+  do a=1,N
+    LT_list(a)=LTave+a*abs(LTave-LTmax)/(N+1)
+    LT_list(N+a)=LTave-a*abs(LTave-LTmax)/(N+1)
+  enddo
+else ! LTave is closer to LTmin than to LTmax
+  do a=1,N
+    LT_list(a)=LTave+a*abs(LTave-LTmin)/(N+1)
+    LT_list(N+a)=LTave-a*abs(LTave-LTmin)/(N+1)
+  enddo
+endif
+CASE ("ntemp")
+  if (abs(f_LT(LTave)-f_LT(LTmax)).le.abs(f_LT(LTave)-f_LT(LTmin))) then ! LTave is closer to LTmax than to LTmin
+    do a=1,N
+      LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
+      LT_list(a)=mod(LT_list(a),24.) ! reput the LT in a [0;24[ interval if needed
+      LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmax))/(N+1)
+      LT_list(N+a)=mod(LT_list(N+a),24.)
+    enddo
+  else ! LTave is closer to LTmin than to LTmax
+    do a=1,N
+      LT_list(a)=LTave+a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
+      LT_list(a)=mod(LT_list(a),24.)
+      LT_list(N+a)=LTave-a*abs(f_LT(LTave)-f_LT(LTmin))/(N+1)
+      LT_list(N+a)=mod(LT_list(N+a),24.)
+    enddo
+  endif
+END SELECT
+
+if (mod(LTcount,2).eq.1) then ! if LTcount is an odd number
+  LT_list(LTcount)=LTave ! add LTave to the list
+endif
+
+!================================================================
+! 2. Combination of int_sol_list & LT_list into sol_list
+!================================================================
+c=1
+do a=1,solcount
+  do b=1,LTcount
+    sol_list(c)=int_sol_list(a)+(LT_list(b)-lon/15.)/24.
+    c=c+1
+  enddo
+enddo
+
+end subroutine gen_sol_list
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine status_check(status,error_text)
+
+use netcdf
+implicit none
+!================================================================
+!  Arguments:
+!================================================================
+integer,intent(in) :: status
+character(len=256),intent(in) :: error_text
+
+if (status.ne.nf90_noerr) then
+  write(*,*)trim(error_text)
+  write(*,*)trim(nf90_strerror(status))
+  stop
+endif
+
+end subroutine status_check
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+real function LTmod(LT)
+!================================================================
+! Purpose:
+! For night local times management, replace hours  
+! from a [0;24[ interval to a [-12;12[ interval in which
+! midnight = 0 (in order to ensure continuity when comparing
+! night local times)
+!================================================================
+real :: LT
+
+LTmod = 2*mod(LT,12.)-LT
+return
+end function LTmod
Index: /trunk/LMDZ.MARS/util/simu_MCS_temp.def
===================================================================
--- /trunk/LMDZ.MARS/util/simu_MCS_temp.def	(revision 2202)
+++ /trunk/LMDZ.MARS/util/simu_MCS_temp.def	(revision 2202)
@@ -0,0 +1,8 @@
+MCS_MY29_tempandLT_mo1.nc
+diagfi1_P.nc
+0
+
+##----------------------------------------------------------------
+# 1st line is the Observer file (containing dtemp,dtimeave,dtimemax,dtimemin,ntemp,ntimeave,ntimemax,ntimemin)
+# 2nd line is the GCM file : diagfi,concatnc or stats (containing temp and with the same altitude type as MCS)
+# 3rd line is the starttimeoffset of the GCM file (not needed for stats)
