Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3093)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3096)
@@ -111,2 +111,10 @@
 == 18/10/2023 == JBC
 The optional file to define the wanted outputs in "diagpem.nc" is now "diagpem.def" (instead of "diagfi.def") + Some updates in the files of deftank.
+
+== 23/10/2023 == JBC
+The management of files during the chained simulation of PCM/PEM runs has been simplified:
+    - "tmp_PEMyears.txt" and "info_run_PEM.txt" have been merged into one file called "info_PEM.txt";
+    - "reshape_XIOS_output.F90" now creates directly the "data_PCM_Y*.nc" files needed by the PEM;
+    - where it is relevant, 'GCM' has been replaced by 'PCM' in the files naming;
+    - the files in deftank have been updated consequently.
+Following r3095, 'iniorbit' is now a subroutine of "planete_h.F90".
Index: trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3093)
+++ trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 3096)
@@ -42,7 +42,7 @@
 
 ! #---------- Martian years parameters from launching script
-open(100,file = 'tmp_PEMyears.txt',status = 'old',form = 'formatted',iostat = ierr)
+open(100,file = 'info_PEM.txt',status = 'old',form = 'formatted',iostat = ierr)
 if (ierr /= 0) then
-    write(*,*) 'Cannot find required file "tmp_PEMyears.txt"!'
+    write(*,*) 'Cannot find required file "info_PEM.txt"!'
     write(*,*) 'It should be created by the launching script...'
     stop
Index: trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3096)
+++ trunk/LMDZ.COMMON/libf/evolution/info_PEM_mod.F90	(revision 3096)
@@ -0,0 +1,53 @@
+MODULE info_PEM_mod
+
+implicit none
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE info_PEM(year_iter,criterion_stop,i_myear,n_myear)
+
+!=======================================================================
+!
+! Purpose: Update the first line of "info_PEM.txt" to count the number of simulated Martian years
+!          Write in "info_PEM.txt" the reason why the PEM stopped and the number of simulated years
+!
+! Author: RV, JBC
+!=======================================================================
+
+use time_evol_mod, only: convert_years
+
+implicit none
+
+!----- Arguments
+integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
+integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
+
+!----- Local variables
+logical       :: ok
+integer       :: cstat
+character(10) :: ich1, ich2, fch
+
+!----- Code
+inquire(file = 'info_PEM.txt', exist = ok)
+if (ok) then
+    write(ich1,'(i0)') i_myear
+    write(ich2,'(i0)') n_myear
+    write(fch,'(f0.4)') convert_years ! 4 digits afetr to the right of the decimal point to respect the precision of Martian year in "launch_pem.sh"
+    call execute_command_line('sed -i "1s/.*/'//trim(ich1)//' '//trim(ich2)//' '//trim(fch)//'/" info_PEM.txt',cmdstat = cstat)
+    if (cstat > 0) then
+        error stop 'info_PEM: command exection failed!'
+    else if (cstat < 0) then
+        error stop 'info_PEM: command execution not supported!'
+    endif
+    open(1,file = 'info_PEM.txt',status = "old",position = "append",action = "write")
+    write(1,*) year_iter, i_myear, criterion_stop
+    close(1)
+else
+    error stop 'The file ''info_PEM.txt'' does not exist and cannot be updated!'
+endif
+
+END SUBROUTINE info_PEM
+
+END MODULE info_PEM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/info_run_PEM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/info_run_PEM_mod.F90	(revision 3093)
+++ 	(revision )
@@ -1,46 +1,0 @@
-MODULE info_run_PEM_mod
-
-implicit none
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
-
-!=======================================================================
-!
-! Purpose: Write in a file called info_run_PEM.txt the reason why the PEM did stop and the number of extrapolation year done
-!          Update the file tmp_PEMyears.txt to count the number of simulated Martian years
-!
-! Author: RV, JBC
-!=======================================================================
-
-use time_evol_mod, only: convert_years
-
-implicit none
-
-!----- Arguments
-integer, intent(in) :: year_iter, criterion_stop ! # of year and reason to stop
-integer, intent(in) :: i_myear, n_myear          ! Current simulated Martian year and maximum number of Martian years to be simulated
-
-!----- Local variables
-logical :: ok
-
-!----- Code
-inquire(file = 'info_run_PEM.txt', exist = ok)
-if (ok) then
-    open(12,file = 'info_run_PEM.txt',status = "old",position = "append",action = "write")
-else
-    open(12,file = 'info_run_PEM.txt',status = "new",action = "write")
-endif
-write(12,*) year_iter, i_myear, criterion_stop
-close(12)
-
-open(100,file = 'tmp_PEMyears.txt',status = 'replace')
-write(100,*) i_myear, n_myear, convert_years
-close(100)
-
-END SUBROUTINE info_run_PEM
-
-END MODULE info_run_PEM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/nb_time_step_GCM_mod.F90	(revision 3093)
+++ 	(revision )
@@ -1,100 +1,0 @@
-MODULE nb_time_step_GCM_mod
-
-use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
-                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
-                  nf90_inquire_dimension, nf90_close
-
-implicit none
-
-character(256) :: msg, var, modname
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE nb_time_step_GCM(fichnom,timelen)
-
-implicit none
-
-!=======================================================================
-!
-! Purpose: Read in the data_GCM_Yr*.nc the number of time step
-!
-! Author: RV
-!=======================================================================
-
-include "dimensions.h"
-
-!=======================================================================
-! Arguments:
-character(len = *), intent(in) :: fichnom !--- FILE NAME
-!=======================================================================
-!   Local Variables
-integer        :: iq, fID, vID, idecal, ierr
-integer        :: timelen ! number of times stored in the file
-!-----------------------------------------------------------------------
-modname = "nb_time_step_GCM"
-
-!  Open initial state NetCDF file
-var = fichnom
-call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-
-ierr = nf90_inq_varid (fID, "temps", vID)
-if (ierr /= nf90_noerr) then
-    write(*,*)"read_data_GCM: Le champ <temps> est absent"
-    write(*,*)"read_data_GCM: J essaie <time_counter>"
-    ierr = nf90_inq_varid (fID, "time_counter", vID)
-    if (ierr /= nf90_noerr) then
-        write(*,*)"read_data_GCM: Le champ <time_counter> est absent"
-        write(*,*)"read_data_GCM: J essaie <Time>"
-        ierr = nf90_inq_varid (fID, "Time", vID)
-        if (ierr /= nf90_noerr) then
-            write(*,*)"read_data_GCM: Le champ <Time> est absent"
-            write(*,*)trim(nf90_strerror(ierr))
-            call abort_gcm("nb_time_step_GCM", "", 1)
-        endif
-        ! Get the length of the "Time" dimension
-        ierr = nf90_inq_dimid(fID,"Time",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
-    else
-        ! Get the length of the "time_counter" dimension
-        ierr = nf90_inq_dimid(fID,"time_counter",vID)
-        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
-    endif
-else
-    ! Get the length of the "temps" dimension
-    ierr = nf90_inq_dimid(fID,"temps",vID)
-    ierr = nf90_inquire_dimension(fID,vID,len = timelen)
-endif
-
-call error_msg(NF90_CLOSE(fID),"close",fichnom)
-
-write(*,*) "The number of timestep of the PCM run data=", timelen
-
-END SUBROUTINE nb_time_step_GCM
-
-!=======================================================================
-
-SUBROUTINE error_msg(ierr,typ,nam)
-
-implicit none
-
-integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
-character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
-character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
-
-if (ierr == nf90_noerr) return
-select case(typ)
-    case('inq');   msg = "Field <"//trim(nam)//"> is missing"
-    case('get');   msg = "Reading failed for <"//trim(nam)//">"
-    case('open');  msg = "File opening failed for <"//trim(nam)//">"
-    case('close'); msg = "File closing failed for <"//trim(nam)//">"
-    case default
-        write(*,*) 'There is no message for this error.'
-        error stop
-end select
-call abort_gcm(trim(modname),trim(msg),ierr)
-
-END SUBROUTINE error_msg
-
-END MODULE nb_time_step_GCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90	(revision 3096)
+++ trunk/LMDZ.COMMON/libf/evolution/nb_time_step_PCM_mod.F90	(revision 3096)
@@ -0,0 +1,100 @@
+MODULE nb_time_step_PCM_mod
+
+use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
+                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
+                  nf90_inquire_dimension, nf90_close
+
+implicit none
+
+character(256) :: msg, var, modname
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE nb_time_step_PCM(fichnom,timelen)
+
+implicit none
+
+!=======================================================================
+!
+! Purpose: Read in the data_PCM_Yr*.nc the number of time step
+!
+! Author: RV
+!=======================================================================
+
+include "dimensions.h"
+
+!=======================================================================
+! Arguments:
+character(len = *), intent(in) :: fichnom !--- FILE NAME
+!=======================================================================
+!   Local Variables
+integer        :: iq, fID, vID, idecal, ierr
+integer        :: timelen ! number of times stored in the file
+!-----------------------------------------------------------------------
+modname = "nb_time_step_PCM"
+
+!  Open initial state NetCDF file
+var = fichnom
+call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+ierr = nf90_inq_varid (fID, "temps", vID)
+if (ierr /= nf90_noerr) then
+    write(*,*)"read_data_PCM: Le champ <temps> est absent"
+    write(*,*)"read_data_PCM: J essaie <time_counter>"
+    ierr = nf90_inq_varid (fID, "time_counter", vID)
+    if (ierr /= nf90_noerr) then
+        write(*,*)"read_data_PCM: Le champ <time_counter> est absent"
+        write(*,*)"read_data_PCM: J essaie <Time>"
+        ierr = nf90_inq_varid (fID, "Time", vID)
+        if (ierr /= nf90_noerr) then
+            write(*,*)"read_data_PCM: Le champ <Time> est absent"
+            write(*,*)trim(nf90_strerror(ierr))
+            call abort_gcm("nb_time_step_PCM", "", 1)
+        endif
+        ! Get the length of the "Time" dimension
+        ierr = nf90_inq_dimid(fID,"Time",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
+    else
+        ! Get the length of the "time_counter" dimension
+        ierr = nf90_inq_dimid(fID,"time_counter",vID)
+        ierr = nf90_inquire_dimension(fID,vID,len = timelen)
+    endif
+else
+    ! Get the length of the "temps" dimension
+    ierr = nf90_inq_dimid(fID,"temps",vID)
+    ierr = nf90_inquire_dimension(fID,vID,len = timelen)
+endif
+
+call error_msg(NF90_CLOSE(fID),"close",fichnom)
+
+write(*,*) "The number of timestep of the PCM run data=", timelen
+
+END SUBROUTINE nb_time_step_PCM
+
+!=======================================================================
+
+SUBROUTINE error_msg(ierr,typ,nam)
+
+implicit none
+
+integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
+character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
+character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
+
+if (ierr == nf90_noerr) return
+select case(typ)
+    case('inq');   msg = "Field <"//trim(nam)//"> is missing"
+    case('get');   msg = "Reading failed for <"//trim(nam)//">"
+    case('open');  msg = "File opening failed for <"//trim(nam)//">"
+    case('close'); msg = "File closing failed for <"//trim(nam)//">"
+    case default
+        write(*,*) 'There is no message for this error.'
+        error stop
+end select
+call abort_gcm(trim(modname),trim(msg),ierr)
+
+END SUBROUTINE error_msg
+
+END MODULE nb_time_step_PCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3093)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3096)
@@ -4,5 +4,5 @@
 !    I_b READ of start_evol.nc and starfi_evol.nc
 !    I_c Subslope parametrisation
-!    I_d READ GCM data and convert to the physical grid
+!    I_d READ PCM data and convert to the physical grid
 !    I_e Initialization of the PEM variable and soil
 !    I_f Compute tendencies & Save initial situation
@@ -48,5 +48,5 @@
                                         TI_PEM, inertiedat_PEM,           & ! soil thermal inertia
                                         tsoil_PEM, mlayer_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths
-                                        fluxgeo,                          & ! Geothermal flux for the PEM and GCM
+                                        fluxgeo,                          & ! Geothermal flux for the PEM and PCM
                                         water_reservoir                     ! Water ressources
 use adsorption_mod,               only: regolith_adsorption, adsorption_pem,        & ! Bool to check if adsorption, main subroutine
@@ -63,9 +63,9 @@
 use soil_settings_PEM_mod,        only: soil_settings_PEM
 use compute_tendencies_slope_mod, only: compute_tendencies_slope
-use info_run_PEM_mod,             only: info_run_PEM
+use info_PEM_mod,                 only: info_PEM
 use interpolate_TIPEM_TIGCM_mod,  only: interpolate_TIPEM_TIGCM
-use nb_time_step_GCM_mod,         only: nb_time_step_GCM
+use nb_time_step_PCM_mod,         only: nb_time_step_PCM
 use pemetat0_mod,                 only: pemetat0
-use read_data_GCM_mod,            only: read_data_GCM
+use read_data_PCM_mod,            only: read_data_PCM
 use recomp_tend_co2_slope_mod,    only: recomp_tend_co2_slope
 use soil_pem_compute_mod,         only: soil_pem_compute
@@ -83,5 +83,5 @@
     use tracer_mod,         only: noms,igcm_h2o_ice, igcm_co2, mmol, igcm_h2o_vap ! Tracer names and molar masses
     use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
-    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit
+    use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit, iniorbit
     use paleoclimate_mod,   only: albedo_perenialco2
     use comcstfi_h,         only: pi, rad, g, cpp, mugaz, r
@@ -120,5 +120,5 @@
 parameter(ngridmx = 2 + (jjm - 1)*iim - 1/jjm)
 
-! Same variable names as in the GCM
+! Same variable names as in the PCM
 integer, parameter :: nlayer = llm ! Number of vertical layer
 integer            :: ngrid        ! Number of physical grid points
@@ -126,7 +126,7 @@
 integer            :: day_ini      ! First day of the simulation
 real               :: pday         ! Physical day
-real               :: time_phys    ! Same as GCM
-real               :: ptimestep    ! Same as GCM
-real               :: ztime_fin    ! Same as GCM
+real               :: time_phys    ! Same as PCM
+real               :: ptimestep    ! Same as PCM
+real               :: ztime_fin    ! Same as PCM
 
 ! Variables to read start.nc
@@ -166,5 +166,5 @@
 real                                :: global_ave_press_new ! constant: Global average pressure of current time step
 real, dimension(:,:),   allocatable :: zplev_new            ! Physical x Atmospheric field : mass of the atmospheric  layers in the pem at current time step [kg/m^2]
-real, dimension(:,:),   allocatable :: zplev_gcm            ! same but retrieved from the gcm [kg/m^2]
+real, dimension(:,:),   allocatable :: zplev_gcm            ! same but retrieved from the PCM [kg/m^2]
 real, dimension(:,:,:), allocatable :: zplev_new_timeseries ! Physical x Atmospheric x Time: same as zplev_new, but in times series [kg/m ^2]
 real, dimension(:,:,:), allocatable :: zplev_old_timeseries ! same but with the time series, for oldest time step
@@ -175,8 +175,8 @@
 integer                             :: criterion_stop       ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
 real, save                          :: A, B, mmean          ! Molar mass: intermediate A, B for computations of the  mean molar mass of the layer [mol/kg]
-real, dimension(:,:),   allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
+real, dimension(:,:),   allocatable :: vmr_co2_gcm          ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
 real, dimension(:,:),   allocatable :: vmr_co2_pem_phys     ! Physics x Times  co2 volume mixing ratio used in the PEM
-real, dimension(:,:),   allocatable :: q_co2_PEM_phys       ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from GCM [kg/kg]
-real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from GCM [kg/kg]
+real, dimension(:,:),   allocatable :: q_co2_PEM_phys       ! Physics x Times co2 mass mixing ratio in the first layer computed in the PEM, first value comes from PCM [kg/kg]
+real, dimension(:,:),   allocatable :: q_h2o_PEM_phys       ! Physics x Times: h2o mass mixing ratio computed in the PEM, first value comes from PCM [kg/kg]
 integer                             :: timelen              ! # time samples
 real                                :: ave                  ! intermediate varibale to compute average
@@ -189,10 +189,10 @@
 real, dimension(:,:),   allocatable :: min_h2o_ice_1          ! ngrid field: minimum of water ice at each point for the first year [kg/m^2]
 real, dimension(:,:),   allocatable :: min_h2o_ice_2          ! ngrid field: minimum of water ice at each point for the second year [kg/m^2]
-real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field: co2 ice given by the GCM  [kg/m^2]
+real, dimension(:,:,:), allocatable :: co2_ice_GCM            ! Physics x NSLOPE x Times field: co2 ice given by the PCM  [kg/m^2]
 real, dimension(:,:),   allocatable :: initial_co2_ice_sublim ! physical point field: Logical array indicating sublimating point of co2 ice
 real, dimension(:,:),   allocatable :: initial_h2o_ice        ! physical point field: Logical array indicating if there is water ice at initial state
 real, dimension(:,:),   allocatable :: initial_co2_ice        ! physical point field: Logical array indicating if there is co2 ice at initial state
 real, dimension(:,:),   allocatable :: tendencies_co2_ice     ! physical point x slope field: Tendency of evolution of perenial co2 ice over a year
-real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the GCM
+real, dimension(:,:),   allocatable :: tendencies_co2_ice_ini ! physical point x slope field x nslope: Tendency of evolution of perenial co2 ice over a year in the PCM
 real, dimension(:,:),   allocatable :: tendencies_h2o_ice     ! physical point x slope  field: Tendency of evolution of perenial h2o ice
 real, dimension(:,:),   allocatable :: flag_co2flow           ! (ngrid,nslope): Flag where there is a CO2 glacier flow
@@ -207,5 +207,5 @@
 real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
 real, dimension(:,:,:,:), allocatable :: tsoil_GCM_timeseries               ! IG x SLOPE XTULES field : NOn averaged Soil Temperature [K]
-real, dimension(:,:),     allocatable :: tsurf_ave_yr1                      ! Physic x SLOPE field : Averaged Surface Temperature of first call of the gcm [K]
+real, dimension(:,:),     allocatable :: tsurf_ave_yr1                      ! Physic x SLOPE field : Averaged Surface Temperature of first call of the PCM [K]
 real, dimension(:,:),     allocatable :: TI_locslope                        ! Physic x Soil: Intermediate thermal inertia  to compute Tsoil [SI]
 real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: intermediate when computing Tsoil [K]
@@ -346,5 +346,5 @@
 #endif
 
-! In the gcm, these values are given to the physic by the dynamic.
+! In the PCM, these values are given to the physic by the dynamic.
 ! Here we simply read them in the startfi_evol.nc file
 status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
@@ -459,8 +459,8 @@
 !------------------------
 ! I   Initialization
-!    I_d READ GCM data and convert to the physical grid
-!------------------------
-! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the GCM run, saving only the minimum value
-call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
+!    I_d READ PCM data and convert to the physical grid
+!------------------------
+! First we read the evolution of water and co2 ice (and the mass mixing ratio) over the first year of the PCM run, saving only the minimum value
+call nb_time_step_PCM("data_PCM_Y1.nc",timelen)
 
 allocate(tsoil_ave(ngrid,nsoilmx,nslope))
@@ -491,12 +491,12 @@
 
 write(*,*) "Downloading data Y1..."
-call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
+call read_data_PCM("data_PCM_Y1.nc",timelen, iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_1,min_h2o_ice_1, &
                    tsurf_ave_yr1,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,           &
                    co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
 write(*,*) "Downloading data Y1 done"
 
-! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the GCM run, saving only the minimum value
+! Then we read the evolution of water and co2 ice (and the mass mixing ratio) over the second year of the PCM run, saving only the minimum value
 write(*,*) "Downloading data Y2"
-call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
+call read_data_PCM("data_PCM_Y2.nc",timelen,iim,jjm_value,ngrid,nslope,vmr_co2_gcm,ps_timeseries,min_co2_ice_2,min_h2o_ice_2, &
                    tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys,              &
                    co2_ice_GCM,watersurf_density_ave,watersoil_density_timeseries)
@@ -1068,5 +1068,5 @@
 enddo
 
-! Conserving the tracers mass for GCM start files
+! Conserving the tracers mass for PCM start files
 do nnq = 1,nqtot
     do ig = 1,ngrid
@@ -1139,5 +1139,5 @@
 write(*,*) "restartpem.nc has been written"
 
-call info_run_PEM(year_iter,criterion_stop,i_myear,n_myear)
+call info_PEM(year_iter,criterion_stop,i_myear,n_myear)
 
 write(*,*) "The PEM has run for", year_iter, "Martian years."
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_GCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_GCM_mod.F90	(revision 3093)
+++ 	(revision )
@@ -1,368 +1,0 @@
-MODULE read_data_GCM_mod
-
-use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
-                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
-                  nf90_inquire_dimension, nf90_close
-
-implicit none
-
-character(256) :: msg, var, modname ! for reading
-integer        :: fID, vID          ! for reading
-
-!=======================================================================
-contains
-!=======================================================================
-
-SUBROUTINE read_data_GCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
-                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
-                         watersurf_density_ave,watersoil_density)
-use comsoil_h,             only: nsoilmx
-use comsoil_h_PEM,         only: soil_pem
-use constants_marspem_mod, only: m_co2, m_noco2
-
-implicit none
-
-!=======================================================================
-!
-! Purpose: Read initial confitions file from the GCM
-!
-! Authors: RV & LL
-!=======================================================================
-
-include "dimensions.h"
-
-!=======================================================================
-! Arguments:
-
-character(len = *), intent(in) :: fichnom                             !--- FILE NAME
-integer,            intent(in) :: timelen                             ! number of times stored in the file
-integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
-! Ouputs
-real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
-real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
-real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the gcm [m^3/m^3]
-real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
-real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
-real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
-real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
-!SOIL
-real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
-real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
-real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
-real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
-real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
-!===============================================================================
-!   Local Variables
-character(12)                   :: start_file_type = "earth" ! default start file type
-real, dimension(:), allocatable :: time      ! times stored in start
-integer                         :: indextime ! index of selected time
-integer                         :: edges(4), corner(4)
-integer                         :: i, j, l, t          ! loop variables
-real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
-integer                         :: islope              ! loop for variables
-character(2)                    :: num                 ! for reading sloped variables
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
-real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
-real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
-
-!-----------------------------------------------------------------------
-modname="read_data_gcm"
-
-A = (1/m_co2 - 1/m_noco2)
-B = 1/m_noco2
-
-write(*,*) "Opening ", fichnom, "..."
-
-!  Open initial state NetCDF file
-var = fichnom
-call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
-
-write(*,*) "Downloading data for vmr co2..."
-
-call get_var3("co2_layer1"   ,q_co2_dyn)
-
-write(*,*) "Downloading data for vmr co2 done"
-write(*,*) "Downloading data for vmr h20..."
-
-call get_var3("h2o_layer1"   ,q_h2o_dyn)
-
-write(*,*) "Downloading data for vmr h2o done"
-write(*,*) "Downloading data for surface pressure ..."
-
-call get_var3("ps"   ,ps_GCM)
-
-write(*,*) "Downloading data for surface pressure done"
-write(*,*) "nslope=", nslope
-
-if (nslope > 1) then
-    write(*,*) "Downloading data for co2ice_slope ..."
-
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
-    enddo
-    
-    write(*,*) "Downloading data for co2ice_slope done"
-    write(*,*) "Downloading data for h2o_ice_s_slope ..."
-    
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
-    enddo
-
-    write(*,*) "Downloading data for h2o_ice_s_slope done"
-
-#ifndef CPP_STD
-    write(*,*) "Downloading data for watercap_slope ..."
-    do islope = 1,nslope
-        write(num,'(i2.2)') islope
-        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
-!        watercap_slope(:,:,:,:) = 0.
-    enddo
-    write(*,*) "Downloading data for watercap_slope done"
-#endif
-
-    write(*,*) "Downloading data for tsurf_slope ..."
-
-    do islope=1,nslope
-        write(num,'(i2.2)') islope
-        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
-    enddo
-
-    write(*,*) "Downloading data for tsurf_slope done"
-
-#ifndef CPP_STD
-    if (soil_pem) then
-        write(*,*) "Downloading data for soiltemp_slope ..."
-
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
-        enddo
-
-        write(*,*) "Downloading data for soiltemp_slope done"
-        write(*,*) "Downloading data for watersoil_density ..."
-
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
-        enddo
-
-        write(*,*) "Downloading data for  watersoil_density  done"
-        write(*,*) "Downloading data for  watersurf_density  ..."
-
-        do islope = 1,nslope
-            write(num,'(i2.2)') islope
-            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
-        enddo
-
-        write(*,*) "Downloading data for  watersurf_density  done"
-
-    endif !soil_pem
-#endif
-else !nslope = 1 no slope, we copy all the values
-    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
-    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
-    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
-    
-#ifndef CPP_STD
-!    call get_var3("watercap", watercap_slope(:,:,1,:))
-    watercap_slope(:,:,1,:) = 0.
-    watersurf_density_dyn(:,:,:,:) = 0.
-    watersoil_density_dyn(:,:,:,:,:) = 0.
-#endif
-
-    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
-endif !nslope=1
-
-! Compute the minimum over the year for each point
-write(*,*) "Computing the min of h2o_ice_slope"
-min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
-write(*,*) "Computing the min of co2_ice_slope"
-min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
-
-!Compute averages
-write(*,*) "Computing average of tsurf"
-tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
-
-#ifndef CPP_1D
-    do islope = 1,nslope
-        do t = 1,timelen
-            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
-        enddo
-    enddo
-#endif
-
-if (soil_pem) then
-    write(*,*) "Computing average of tsoil"
-    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
-    write(*,*) "Computing average of watersurf_density"
-    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
-endif
-
-! By definition, a density is positive, we get rid of the negative values
-where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
-where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
-
-do i = 1,iim + 1
-    do j = 1,jjm_input + 1
-        do t = 1, timelen
-            if (q_co2_dyn(i,j,t) < 0) then
-                q_co2_dyn(i,j,t) = 1.e-10
-            else if (q_co2_dyn(i,j,t) > 1) then
-                q_co2_dyn(i,j,t) = 1.
-            endif
-            if (q_h2o_dyn(i,j,t) < 0) then
-                q_h2o_dyn(i,j,t) = 1.e-30
-            else if (q_h2o_dyn(i,j,t) > 1) then
-                q_h2o_dyn(i,j,t) = 1.
-            endif
-            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
-            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
-        enddo
-    enddo
-enddo
-
-#ifndef CPP_1D
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
-
-    do islope = 1,nslope
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
-        if (soil_pem) then
-            do l = 1,nsoilmx
-                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
-                do t=1,timelen
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
-                enddo
-            enddo
-        endif !soil_pem
-        do t = 1,timelen
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
-        enddo
-    enddo
-
-    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
-#else
-    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
-    ps_timeseries(1,:) = ps_GCM(1,1,:)
-    q_co2(1,:) = q_co2_dyn(1,1,:)
-    q_h2o(1,:) = q_h2o_dyn(1,1,:)
-    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
-    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
-    if (soil_pem) then
-        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
-        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
-        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
-    endif !soil_pem
-    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
-    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
-    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
-#endif
-
-END SUBROUTINE read_data_gcm
-
-SUBROUTINE check_dim(n1,n2,str1,str2)
-
-implicit none
-
-integer,            intent(in) :: n1, n2
-character(len = *), intent(in) :: str1, str2
-
-character(256) :: s1, s2
-
-if (n1 /= n2) then
-    s1 = 'value of '//trim(str1)//' ='
-    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
-    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
-    call abort_gcm(trim(modname),trim(msg),1)
-endif
-
-END SUBROUTINE check_dim
-
-!=======================================================================
-
-SUBROUTINE get_var1(var,v)
-
-implicit none
-
-character(len = *), intent(in)  :: var
-real, dimension(:), intent(out) :: v
-
-call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
-
-END SUBROUTINE get_var1
-
-!=======================================================================
-
-SUBROUTINE get_var3(var,v) ! on U grid
-
-implicit none
-
-character(len = *),     intent(in)  :: var
-real, dimension(:,:,:), intent(out) :: v
-
-call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
-
-END SUBROUTINE get_var3
-
-!=======================================================================
-
-SUBROUTINE get_var4(var,v)
-
-implicit none
-
-character(len = *),       intent(in)  :: var
-real, dimension(:,:,:,:), intent(out) :: v
-
-call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
-call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
-
-END SUBROUTINE get_var4
-
-!=======================================================================
-
-SUBROUTINE error_msg(ierr,typ,nam)
-
-implicit none
-
-integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
-character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
-character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
-
-if (ierr == nf90_noerr) return
-select case(typ)
-    case('inq');   msg="Field <"//trim(nam)//"> is missing"
-    case('get');   msg="Reading failed for <"//trim(nam)//">"
-    case('open');  msg="File opening failed for <"//trim(nam)//">"
-    case('close'); msg="File closing failed for <"//trim(nam)//">"
-    case default
-        write(*,*) 'There is no message for this error.'
-        error stop
-end select
-call abort_gcm(trim(modname),trim(msg),ierr)
-
-END SUBROUTINE error_msg
-
-END MODULE read_data_GCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3096)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3096)
@@ -0,0 +1,368 @@
+MODULE read_data_PCM_mod
+
+use netcdf, only: nf90_open, NF90_NOWRITE, nf90_noerr, nf90_strerror, &
+                  nf90_get_var, nf90_inq_varid, nf90_inq_dimid,       &
+                  nf90_inquire_dimension, nf90_close
+
+implicit none
+
+character(256) :: msg, var, modname ! for reading
+integer        :: fID, vID          ! for reading
+
+!=======================================================================
+contains
+!=======================================================================
+
+SUBROUTINE read_data_PCM(fichnom,timelen,iim_input,jjm_input,ngrid,nslope,vmr_co2_gcm_phys,ps_timeseries,           &
+                         min_co2_ice,min_h2o_ice,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,q_co2,q_h2o,co2_ice_slope, &
+                         watersurf_density_ave,watersoil_density)
+use comsoil_h,             only: nsoilmx
+use comsoil_h_PEM,         only: soil_pem
+use constants_marspem_mod, only: m_co2, m_noco2
+
+implicit none
+
+!=======================================================================
+!
+! Purpose: Read initial confitions file from the PCM
+!
+! Authors: RV & LL
+!=======================================================================
+
+include "dimensions.h"
+
+!=======================================================================
+! Arguments:
+
+character(len = *), intent(in) :: fichnom                             !--- FILE NAME
+integer,            intent(in) :: timelen                             ! number of times stored in the file
+integer                        :: iim_input, jjm_input, ngrid, nslope ! number of points in the lat x lon dynamical grid, number of subgrid slopes
+! Ouputs
+real, dimension(ngrid,nslope),         intent(out) :: min_co2_ice      ! Minimum of co2 ice  per slope of the year [kg/m^2]
+real, dimension(ngrid,nslope),         intent(out) :: min_h2o_ice      ! Minimum of h2o ice per slope of the year [kg/m^2]
+real, dimension(ngrid,timelen),        intent(out) :: vmr_co2_gcm_phys ! Physics x Times  co2 volume mixing ratio retrieve from the PCM [m^3/m^3]
+real, dimension(ngrid,timelen),        intent(out) :: ps_timeseries    ! Surface Pressure [Pa]
+real, dimension(ngrid,timelen),        intent(out) :: q_co2            ! CO2 mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,timelen),        intent(out) :: q_h2o            ! H2O mass mixing ratio in the first layer [kg/m^3]
+real, dimension(ngrid,nslope,timelen), intent(out) :: co2_ice_slope    ! co2 ice amount per  slope of the year [kg/m^2]
+!SOIL
+real, dimension(ngrid,nslope),                 intent(out) :: tsurf_ave             ! Average surface temperature of the concatenated file [K]
+real, dimension(ngrid,nsoilmx,nslope),         intent(out) :: tsoil_ave             ! Average soil temperature of the concatenated file [K]
+real, dimension(ngrid,nslope,timelen),         intent(out) :: tsurf_gcm             ! Surface temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: tsoil_gcm             ! Soil temperature of the concatenated file, time series [K]
+real, dimension(ngrid,nslope),                 intent(out) :: watersurf_density_ave ! Water density at the surface [kg/m^3]
+real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density     ! Water density in the soil layer, time series [kg/m^3]
+!===============================================================================
+!   Local Variables
+character(12)                   :: start_file_type = "earth" ! default start file type
+real, dimension(:), allocatable :: time      ! times stored in start
+integer                         :: indextime ! index of selected time
+integer                         :: edges(4), corner(4)
+integer                         :: i, j, l, t          ! loop variables
+real                            ::  A , B, mmean       ! Molar Mass of co2 and no co2, A;B intermediate variables to compute the mean molar mass of the layer
+integer                         :: islope              ! loop for variables
+character(2)                    :: num                 ! for reading sloped variables
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: h2o_ice_s_dyn         ! h2o ice per slope of the concatenated file [kg/m^2]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watercap_slope
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_gcm           ! CO2 volume mixing ratio in the first layer  [mol/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_GCM                ! Surface Pressure [Pa]
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_co2_ice_dyn
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: min_h2o_ice_dyn
+real, dimension(iim_input + 1,jjm_input + 1,nslope)                 :: tsurf_ave_dyn         ! Average surface temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_ave_dyn         ! Average soil temperature of the concatenated file [K]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_gcm_dyn         ! Surface temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_gcm_dyn         ! Soil temperature of the concatenated file, time series [K]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_co2_dyn             ! CO2 mass mixing ratio in the first layer [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: q_h2o_dyn             ! H2O mass mixing ratio in the first layer [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: co2_ice_slope_dyn     ! co2 ice amount per  slope of the year [kg/m^2]
+real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: watersurf_density_dyn ! Water density at the surface, time series [kg/m^3]
+real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: watersoil_density_dyn ! Water density in the soil layer, time series [kg/m^3]
+real, dimension(ngrid,nslope,timelen)                               :: watersurf_density     ! Water density at the surface, time series [kg/m^3]
+
+!-----------------------------------------------------------------------
+modname="read_data_PCM"
+
+A = (1/m_co2 - 1/m_noco2)
+B = 1/m_noco2
+
+write(*,*) "Opening ", fichnom, "..."
+
+!  Open initial state NetCDF file
+var = fichnom
+call error_msg(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
+
+write(*,*) "Downloading data for vmr co2..."
+
+call get_var3("co2_layer1"   ,q_co2_dyn)
+
+write(*,*) "Downloading data for vmr co2 done"
+write(*,*) "Downloading data for vmr h20..."
+
+call get_var3("h2o_layer1"   ,q_h2o_dyn)
+
+write(*,*) "Downloading data for vmr h2o done"
+write(*,*) "Downloading data for surface pressure ..."
+
+call get_var3("ps"   ,ps_GCM)
+
+write(*,*) "Downloading data for surface pressure done"
+write(*,*) "nslope=", nslope
+
+if (nslope > 1) then
+    write(*,*) "Downloading data for co2ice_slope ..."
+
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
+    enddo
+    
+    write(*,*) "Downloading data for co2ice_slope done"
+    write(*,*) "Downloading data for h2o_ice_s_slope ..."
+    
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
+    enddo
+
+    write(*,*) "Downloading data for h2o_ice_s_slope done"
+
+#ifndef CPP_STD
+    write(*,*) "Downloading data for watercap_slope ..."
+    do islope = 1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("watercap_slope"//num,watercap_slope(:,:,islope,:))
+!        watercap_slope(:,:,:,:) = 0.
+    enddo
+    write(*,*) "Downloading data for watercap_slope done"
+#endif
+
+    write(*,*) "Downloading data for tsurf_slope ..."
+
+    do islope=1,nslope
+        write(num,'(i2.2)') islope
+        call get_var3("tsurf_slope"//num,tsurf_gcm_dyn(:,:,islope,:))
+    enddo
+
+    write(*,*) "Downloading data for tsurf_slope done"
+
+#ifndef CPP_STD
+    if (soil_pem) then
+        write(*,*) "Downloading data for soiltemp_slope ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("soiltemp_slope"//num,tsoil_gcm_dyn(:,:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for soiltemp_slope done"
+        write(*,*) "Downloading data for watersoil_density ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var4("Waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for  watersoil_density  done"
+        write(*,*) "Downloading data for  watersurf_density  ..."
+
+        do islope = 1,nslope
+            write(num,'(i2.2)') islope
+            call get_var3("Waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
+        enddo
+
+        write(*,*) "Downloading data for  watersurf_density  done"
+
+    endif !soil_pem
+#endif
+else !nslope = 1 no slope, we copy all the values
+    call get_var3("h2o_ice_s", h2o_ice_s_dyn(:,:,1,:))
+    call get_var3("co2ice", co2_ice_slope_dyn(:,:,1,:))
+    call get_var3("tsurf", tsurf_gcm_dyn(:,:,1,:))
+    
+#ifndef CPP_STD
+!    call get_var3("watercap", watercap_slope(:,:,1,:))
+    watercap_slope(:,:,1,:) = 0.
+    watersurf_density_dyn(:,:,:,:) = 0.
+    watersoil_density_dyn(:,:,:,:,:) = 0.
+#endif
+
+    if (soil_pem) call get_var4("soiltemp",tsoil_gcm_dyn(:,:,:,1,:))
+endif !nslope=1
+
+! Compute the minimum over the year for each point
+write(*,*) "Computing the min of h2o_ice_slope"
+min_h2o_ice_dyn(:,:,:) = minval(h2o_ice_s_dyn + watercap_slope,4)
+write(*,*) "Computing the min of co2_ice_slope"
+min_co2_ice_dyn(:,:,:) = minval(co2_ice_slope_dyn,4)
+
+!Compute averages
+write(*,*) "Computing average of tsurf"
+tsurf_ave_dyn(:,:,:) = sum(tsurf_gcm_dyn(:,:,:,:),4)/timelen
+
+#ifndef CPP_1D
+    do islope = 1,nslope
+        do t = 1,timelen
+            call gr_dyn_fi(1,iim_input+1,jjm_input+1,ngrid,watersurf_density_dyn(:,:,islope,t),watersurf_density(:,islope,t))
+        enddo
+    enddo
+#endif
+
+if (soil_pem) then
+    write(*,*) "Computing average of tsoil"
+    tsoil_ave_dyn(:,:,:,:) = sum(tsoil_gcm_dyn(:,:,:,:,:),5)/timelen
+    write(*,*) "Computing average of watersurf_density"
+    watersurf_density_ave(:,:) = sum(watersurf_density(:,:,:),3)/timelen
+endif
+
+! By definition, a density is positive, we get rid of the negative values
+where (min_co2_ice_dyn < 0.) min_co2_ice_dyn = 0.
+where (min_h2o_ice_dyn < 0.) min_h2o_ice_dyn = 0.
+
+do i = 1,iim + 1
+    do j = 1,jjm_input + 1
+        do t = 1, timelen
+            if (q_co2_dyn(i,j,t) < 0) then
+                q_co2_dyn(i,j,t) = 1.e-10
+            else if (q_co2_dyn(i,j,t) > 1) then
+                q_co2_dyn(i,j,t) = 1.
+            endif
+            if (q_h2o_dyn(i,j,t) < 0) then
+                q_h2o_dyn(i,j,t) = 1.e-30
+            else if (q_h2o_dyn(i,j,t) > 1) then
+                q_h2o_dyn(i,j,t) = 1.
+            endif
+            mmean = 1/(A*q_co2_dyn(i,j,t) + B)
+            vmr_co2_gcm(i,j,t) = q_co2_dyn(i,j,t)*mmean/m_co2
+        enddo
+    enddo
+enddo
+
+#ifndef CPP_1D
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,vmr_co2_gcm,vmr_co2_gcm_phys)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_GCM,ps_timeseries)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_co2_dyn,q_co2)
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,q_h2o_dyn,q_h2o)
+
+    do islope = 1,nslope
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_co2_ice_dyn(:,:,islope),min_co2_ice(:,islope))
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope))
+        if (soil_pem) then
+            do l = 1,nsoilmx
+                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_ave_dyn(:,:,l,islope),tsoil_ave(:,l,islope))
+                do t=1,timelen
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_gcm_dyn(:,:,l,islope,t),tsoil_gcm(:,l,islope,t))
+                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density(:,l,islope,t))
+                enddo
+            enddo
+        endif !soil_pem
+        do t = 1,timelen
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_GCM_dyn(:,:,islope,t),tsurf_GCM(:,islope,t))
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2_ice_slope_dyn(:,:,islope,t),co2_ice_slope(:,islope,t))
+        enddo
+    enddo
+
+    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_ave_dyn,tsurf_ave)
+#else
+    vmr_co2_gcm_phys(1,:) = vmr_co2_gcm(1,1,:)
+    ps_timeseries(1,:) = ps_GCM(1,1,:)
+    q_co2(1,:) = q_co2_dyn(1,1,:)
+    q_h2o(1,:) = q_h2o_dyn(1,1,:)
+    min_co2_ice(1,:) = min_co2_ice_dyn(1,1,:)
+    min_h2o_ice(1,:) = min_h2o_ice_dyn(1,1,:)
+    if (soil_pem) then
+        tsoil_ave(1,:,:) = tsoil_ave_dyn(1,1,:,:)
+        tsoil_gcm(1,:,:,:) = tsoil_gcm_dyn(1,1,:,:,:)
+        watersoil_density(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
+    endif !soil_pem
+    tsurf_GCM(1,:,:) = tsurf_GCM_dyn(1,1,:,:)
+    co2_ice_slope(1,:,:) = co2_ice_slope_dyn(1,1,:,:)
+    tsurf_ave(1,:) = tsurf_ave_dyn(1,1,:)
+#endif
+
+END SUBROUTINE read_data_PCM
+
+SUBROUTINE check_dim(n1,n2,str1,str2)
+
+implicit none
+
+integer,            intent(in) :: n1, n2
+character(len = *), intent(in) :: str1, str2
+
+character(256) :: s1, s2
+
+if (n1 /= n2) then
+    s1 = 'value of '//trim(str1)//' ='
+    s2 = ' read in starting file differs from parametrized '//trim(str2)//' ='
+    write(msg,'(10x,a,i4,2x,a,i4)')trim(s1),n1,trim(s2),n2
+    call abort_gcm(trim(modname),trim(msg),1)
+endif
+
+END SUBROUTINE check_dim
+
+!=======================================================================
+
+SUBROUTINE get_var1(var,v)
+
+implicit none
+
+character(len = *), intent(in)  :: var
+real, dimension(:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var1
+
+!=======================================================================
+
+SUBROUTINE get_var3(var,v) ! on U grid
+
+implicit none
+
+character(len = *),     intent(in)  :: var
+real, dimension(:,:,:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var3
+
+!=======================================================================
+
+SUBROUTINE get_var4(var,v)
+
+implicit none
+
+character(len = *),       intent(in)  :: var
+real, dimension(:,:,:,:), intent(out) :: v
+
+call error_msg(NF90_INQ_VARID(fID,var,vID),"inq",var)
+call error_msg(NF90_GET_VAR(fID,vID,v),"get",var)
+
+END SUBROUTINE get_var4
+
+!=======================================================================
+
+SUBROUTINE error_msg(ierr,typ,nam)
+
+implicit none
+
+integer,            intent(in) :: ierr !--- NetCDF ERROR CODE
+character(len = *), intent(in) :: typ  !--- TYPE OF OPERATION
+character(len = *), intent(in) :: nam  !--- FIELD/FILE NAME
+
+if (ierr == nf90_noerr) return
+select case(typ)
+    case('inq');   msg="Field <"//trim(nam)//"> is missing"
+    case('get');   msg="Reading failed for <"//trim(nam)//">"
+    case('open');  msg="File opening failed for <"//trim(nam)//">"
+    case('close'); msg="File closing failed for <"//trim(nam)//">"
+    case default
+        write(*,*) 'There is no message for this error.'
+        error stop
+end select
+call abort_gcm(trim(modname),trim(msg),ierr)
+
+END SUBROUTINE error_msg
+
+END MODULE read_data_PCM_mod
Index: trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 3093)
+++ trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90	(revision 3096)
@@ -3,8 +3,8 @@
 !=======================================================================
 !
-! Purpose: Read XIOS files, and convert them into the correct GCM grid
+! Purpose: Read XIOS files, and convert them into the correct PCM grid
 !          XIOS  longitudes start at -180 but stop before -180 (not duplicated)
 !          We basically add the last point, and complete the XIOS file. Looped
-!          over the two GCM runs
+!          over the two PCM runs
 !
 ! Authors: RV & LL
@@ -35,5 +35,5 @@
     if (state /= nf90_noerr) call handle_err(state)
 
-    state = nf90_create(path = "datareshaped"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
+    state = nf90_create(path = "data_PCM_Y"//str2//".nc", cmode=or(nf90_noclobber,nf90_64bit_offset), ncid = ncid2)
     if (state /= nf90_noerr) call handle_err(state)