Index: trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3602)
@@ -557,2 +557,7 @@
 == 23/01/2025 == JBC
 Few small optimizations
+
+== 27/01/2025 == JBC
+- Complete rework of "read_data_PCM_mod.F90" to reduce as much as possible memory consumption which caused problems and the model failure.
+- In case of very low surface shifting, checking if the index is inside the bounds of 'tsoil' and 'mlayer' and replace them by the value at surface if necessary.
+- Few cleanings.
Index: trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90	(revision 3602)
@@ -243,5 +243,5 @@
 ! ------
 integer                             :: ig, isoil, islope, ishift, ilag, iz
-real                                :: a, z, zshift_surfloc
+real                                :: a, z, zshift_surfloc, tsoil_minus, mlayer_minus
 real, dimension(ngrid,nsoil,nslope) :: tsoil_old
 
@@ -265,5 +265,4 @@
                 ! The "new soil" temperature is set to tsurf
                 tsoil(ig,:ishift,islope) = tsurf(ig,islope)
-
                 do isoil = ishift + 1,nsoil
                     ! Position in the old discretization of the depth
@@ -274,7 +273,14 @@
                         iz = iz + 1
                     enddo
+                    if (iz == 0) then
+                        tsoil_minus = tsurf(ig,islope)
+                        mlayer_minus = 0.
+                    else
+                        tsoil_minus = tsoil_old(ig,iz,islope)
+                        mlayer_minus = mlayer_PEM(iz - 1)
+                    endif
                     ! Interpolation of the temperature profile from the old discretization
-                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
-                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus)
+                    tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus
                 enddo
             endif
@@ -299,7 +305,14 @@
                         iz = iz + 1
                     enddo
+                    if (iz == 0) then
+                        tsoil_minus = tsurf(ig,islope)
+                        mlayer_minus = 0.
+                    else
+                        tsoil_minus = tsoil_old(ig,iz,islope)
+                        mlayer_minus = mlayer_PEM(iz - 1)
+                    endif
                     ! The "new lag" temperature is set to the ice temperature just below
-                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
-                    tsoil(ig,:ilag,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus)
+                    tsoil(ig,:ilag,islope) = a*(z - mlayer_minus) + tsoil_minus
                 endif
 
@@ -318,7 +331,14 @@
                         iz = iz + 1
                     enddo
+                    if (iz == 0) then
+                        tsoil_minus = tsurf(ig,islope)
+                        mlayer_minus = 0.
+                    else
+                        tsoil_minus = tsoil_old(ig,iz,islope)
+                        mlayer_minus = mlayer_PEM(iz - 1)
+                    endif
                     ! Interpolation of the temperature profile from the old discretization
-                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1))
-                    tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope)
+                    a = (tsoil_old(ig,iz + 1,islope) - tsoil_minus)/(mlayer_PEM(iz) - mlayer_minus)
+                    tsoil(ig,isoil,islope) = a*(z - mlayer_minus) + tsoil_minus
                 enddo
 
Index: trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3602)
@@ -46,5 +46,5 @@
 
 !------------------------------------------------------------------------------------------------------
-SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_beg,ice_table_thickness)
+SUBROUTINE computeice_table_equilibrium(ngrid,nslope,nsoil,watercaptag,rhowatersurf_ave,rhowatersoil_ave,regolith_inertia,ice_table_depth,ice_table_thickness)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -72,5 +72,5 @@
 ! Ouputs
 ! ------
-real, dimension(ngrid,nslope), intent(out) :: ice_table_beg       ! ice table depth [m]
+real, dimension(ngrid,nslope), intent(out) :: ice_table_depth     ! ice table depth [m]
 real, dimension(ngrid,nslope), intent(out) :: ice_table_thickness ! ice table thickness [m]
 ! Locals
@@ -84,22 +84,22 @@
 ! Code
 ! ----
-previous_icetable_depth = ice_table_beg
+previous_icetable_depth = ice_table_depth
 do ig = 1,ngrid
     if (watercaptag(ig)) then
-        ice_table_beg(ig,:) = 0.
+        ice_table_depth(ig,:) = 0.
         ice_table_thickness(ig,:) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
     else
         do islope = 1,nslope
-            ice_table_beg(ig,islope) = -1.
+            ice_table_depth(ig,islope) = -1.
             ice_table_thickness(ig,islope) = 0.
             do isoil = 1,nsoil
                 diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)
             enddo
-            if (diff_rho(1) > 0) then ! ice is at the surface
-                ice_table_beg(ig,islope) = 0.
+            if (diff_rho(1) > 0.) then ! ice is at the surface
+                ice_table_depth(ig,islope) = 0.
                 do isoilend = 2,nsoil - 1
                     if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
-                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
-                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
+                        call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end)
+                        ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                         exit
                     endif
@@ -107,12 +107,12 @@
             else
                 do isoil = 1,nsoil - 1 ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
-                    if (diff_rho(isoil) < 0 .and. diff_rho(isoil + 1) > 0.) then
-                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_beg(ig,islope))
+                    if (diff_rho(isoil) < 0. .and. diff_rho(isoil + 1) > 0.) then
+                        call findroot(diff_rho(isoil),diff_rho(isoil + 1),mlayer_PEM(isoil),mlayer_PEM(isoil + 1),ice_table_depth(ig,islope))
                         ! Now let's find the end of the ice table
                         ice_table_thickness(ig,islope) = layer_PEM(nsoil) ! Let's assume an infinite ice table (true when geothermal flux is set to 0.)
                         do isoilend = isoil + 1,nsoil - 1
-                            if (diff_rho(isoilend) > 0 .and. diff_rho(isoilend + 1) < 0.) then
-                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend),mlayer_PEM(isoilend + 1),ice_table_end)
-                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_beg(ig,islope)
+                            if (diff_rho(isoilend) > 0. .and. diff_rho(isoilend + 1) < 0.) then
+                                call findroot(diff_rho(isoilend),diff_rho(isoilend + 1),mlayer_PEM(isoilend - 1),mlayer_PEM(isoilend),ice_table_end)
+                                ice_table_thickness(ig,islope) = ice_table_end - ice_table_depth(ig,islope)
                                 exit
                             endif
@@ -128,10 +128,10 @@
 do islope = 1,nslope
     do ig = 1,ngrid
-        if (ice_table_beg(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
+        if (ice_table_depth(ig,islope) > previous_icetable_depth(ig,islope) .and. previous_icetable_depth(ig,islope) >= 0) then
             call ice_thermal_properties(.false.,1.,regolith_inertia(ig,islope),wice_inertia)
             stretch = (regolith_inertia(ig,islope)/wice_inertia)**2
-            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_beg(ig,islope) - &
-                                             previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
-            ice_table_beg(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_beg(ig,islope) - previous_icetable_depth(ig,islope))/stretch
+            ice_table_thickness(ig,islope) = ice_table_thickness(ig,islope) + (ice_table_depth(ig,islope) - &
+                                             previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch)
+            ice_table_depth(ig,islope) = previous_icetable_depth(ig,islope) + (ice_table_depth(ig,islope) - previous_icetable_depth(ig,islope))/stretch
         endif
     enddo
Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3602)
@@ -615,8 +615,7 @@
 
 delta_h2o_icetablesublim = 0.
-call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth, &
-              icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,        &
-              ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,watersurf_density_avg,         &
-              watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
+call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
+              tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,ps_avg_global_ini,d_h2oice,d_co2ice,co2_ice,h2o_ice,            &
+              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys,delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif)
 deallocate(tsurf_avg_yr1)
 
Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3602)
@@ -7,8 +7,7 @@
 !=======================================================================
 
-SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness, &
-                    ice_porefilling,tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,                           &
-                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,                                      &
-                    d_h2oice,d_co2ice,co2_ice,h2o_ice,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
+SUBROUTINE pemetat0(filename,ngrid,nsoil_PCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,icetable_depth,icetable_thickness,ice_porefilling, &
+                    tsurf_avg_yr1,tsurf_avg_yr2,q_co2,q_h2o,ps_timeseries,ps_avg_global,d_h2oice,d_co2ice,co2_ice,h2o_ice,                         &
+                    watersurf_avg,watersoil_avg,m_co2_regolith_phys,deltam_co2_regolith_phys,m_h2o_regolith_phys,deltam_h2o_regolith_phys,stratif)
 
 use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var, inquire_dimension, inquire_dimension_length
@@ -210,5 +209,5 @@
         enddo ! islope
 
-        write(*,*) 'PEMETAT0: THERMAL INERTIA done'
+        write(*,*) 'PEMETAT0: TI done'
 
 ! b. Special case for inertiedat, inertiedat_PEM
@@ -277,5 +276,5 @@
             enddo
         enddo ! islope
-        write(*,*) 'PEMETAT0: SOIL TEMP done'
+        write(*,*) 'PEMETAT0: TSOIL done'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -367,12 +366,12 @@
                     !!! transition
                     delta = depth_breccia
-                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
-                                                               (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
-                                                               ((layer_PEM(index_breccia+1)-delta)/(TI_breccia**2))))
-                    do iloop=index_breccia + 2,index_bedrock
+                    TI_PEM(ig,index_breccia + 1,islope) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                   &
+                                                               (((delta - layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2)) + &
+                                                               ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
+                    do iloop = index_breccia + 2,index_bedrock
                         TI_PEM(ig,iloop,islope) = TI_breccia
                     enddo
                 else ! we keep the high ti values
-                    do iloop=index_breccia + 1,index_bedrock
+                    do iloop = index_breccia + 1,index_bedrock
                         TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
                     enddo
@@ -380,7 +379,7 @@
                 !! transition
                 delta = depth_bedrock
-                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
-                                                           (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
-                                                           ((layer_PEM(index_bedrock+1)-delta)/(TI_breccia**2))))
+                TI_PEM(ig,index_bedrock + 1,islope) = sqrt((layer_PEM(index_bedrock + 1) - layer_PEM(index_bedrock))/                   &
+                                                           (((delta - layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2)) + &
+                                                           ((layer_PEM(index_bedrock + 1) - delta)/(TI_breccia**2))))
                 do iloop = index_bedrock + 2,nsoil_PEM
                     TI_PEM(ig,iloop,islope) = TI_bedrock
@@ -398,5 +397,5 @@
                 inertiedat_PEM(ig,index_breccia + 1) = sqrt((layer_PEM(index_breccia + 1) - layer_PEM(index_breccia))/                &
                                                             (((delta - layer_PEM(index_breccia))/(inertiedat(ig,index_breccia)**2)) + &
-                                                            ((layer_PEM(index_breccia + 1)-delta)/(TI_breccia**2))))
+                                                            ((layer_PEM(index_breccia + 1) - delta)/(TI_breccia**2))))
                 do iloop = index_breccia + 2,index_bedrock
                     inertiedat_PEM(ig,iloop) = TI_breccia
@@ -448,5 +447,5 @@
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             enddo
-            write(*,*) 'PEMETAT0: Ice table done'
+            write(*,*) 'PEMETAT0: ICE TABLE done'
         else if (icetable_dynamic) then
             call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,d_h2oice,h2o_ice,ps_avg_global,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM)
@@ -454,5 +453,5 @@
                 call ini_tsoil_pem(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_avg_yr2(:,islope),tsoil_PEM(:,:,islope))
             enddo
-            write(*,*) 'PEMETAT0: Ice table done'
+            write(*,*) 'PEMETAT0: ICE TABLE done'
         endif
 
@@ -468,5 +467,5 @@
         endif
 
-        write(*,*) 'PEMETAT0: CO2 adsorption done'
+        write(*,*) 'PEMETAT0: CO2 & H2O adsorption done'
     endif !soil_pem
 endif ! of if (startphy_file)
Index: trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3601)
+++ trunk/LMDZ.COMMON/libf/evolution/read_data_PCM_mod.F90	(revision 3602)
@@ -27,5 +27,5 @@
 ! Purpose: Read initial confitions file from the PCM
 !
-! Authors: RV & LL
+! Authors: JBC
 !=======================================================================
 
@@ -52,28 +52,14 @@
 real, dimension(ngrid,nsoilmx,nslope,timelen), intent(out) :: watersoil_density_timeseries ! Water density timeseries in the soil layer [kg/m^3]
 ! Local variables
-integer                                                             :: i, j, l, t, islope        ! 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
-character(2)                                                        :: num                       ! Hor 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
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: perennial_co2ice
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: vmr_co2_dyn               ! CO2 volume mixing ratio in the first layer [mol/m^3]
-real, dimension(iim_input + 1,jjm_input + 1,timelen)                :: ps_dyn                   ! Surface pressure [Pa]
-real, dimension(iim_input + 1,jjm_input + 1)                        :: ps_avg_dyn                ! Averaged surface pressure of the concatenated file [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_avg_dyn             ! Averaged surface temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope)         :: tsoil_avg_dyn             ! Averaged soil temperature of the concatenated file [K]
-real, dimension(iim_input + 1,jjm_input + 1,nslope,timelen)         :: tsurf_dyn                ! Surface temperature of the concatenated file, time series [K]
-real, dimension(iim_input + 1,jjm_input + 1,nsoilmx,nslope,timelen) :: tsoil_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,nslope)                 :: watersurf_density_dyn_avg ! Water density at the surface, dynamic grid, yearly averaged [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]
+integer                                                             :: i, j, l, islope ! Loop variables
+real                                                                :: A, B            ! Intermediate variables to compute the mean molar mass of the layer
+character(2)                                                        :: num             ! For reading sloped variables
+
+real, dimension(:,:),     allocatable :: var2_read                             ! Variables for reading (2 dimensions)
+real, dimension(:,:,:),   allocatable :: var3_read_1, var3_read_2, var3_read_3 ! Variables for reading (3 dimensions)
+real, dimension(:,:,:,:), allocatable :: var4_read                             ! Variables for reading (4 dimensions)
 !-----------------------------------------------------------------------
 
-!------------------ Year 1 ------------------
+!------------------------------- Year 1 --------------------------------
 ! Open the NetCDF file
 write(*,*) "Opening "//filename1//"..."
@@ -81,70 +67,121 @@
 
 ! Download the data from the file
+allocate(var2_read(iim_input + 1,jjm_input + 1),var3_read_1(iim_input + 1,jjm_input + 1,timelen),var3_read_2(iim_input + 1,jjm_input + 1,timelen))
+
 if (nslope == 1) then ! There is no slope
-    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
+    ! CO2 ice
+    !--------
+    call get_var3("co2ice",var3_read_1)
+    where (var3_read_1 < 0.) var3_read_1 = 0.
     write(*,*) "Data for co2_ice downloaded."
-
-    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
+#ifndef CPP_STD
+    call get_var3("perennial_co2ice",var3_read_2)
+    write(*,*) "Data for perennial_co2ice downloaded."
+#endif
+
+    ! Compute the minimum over the year for each point
+    var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1))
+#else
+    min_co2_ice(1,1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of co2_ice computed."
+
+    ! H2O ice
+    !--------
+    call get_var3("h2o_ice_s",var3_read_1)
+    where (var3_read_1 < 0.) var3_read_1 = 0.
     write(*,*) "Data for h2o_ice_s downloaded."
-
-    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
+#ifndef CPP_STD
+    call get_var3("watercap",var3_read_2)
+    write(*,*) "Data for watercap downloaded."
+#endif
+
+    ! Compute the minimum over the year for each point
+    var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1))
+#else
+    min_h2o_ice(1,1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of h2o_ice computed."
+
+    ! Tsurf
+    !------
+    call get_var3("tsurf",var3_read_1)
     write(*,*) "Data for tsurf downloaded."
 
-#ifndef CPP_STD
-    call get_var3("watercap",watercap(:,:,1,:))
-    write(*,*) "Data for watercap downloaded."
-
-    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
-    write(*,*) "Data for perennial_co2ice downloaded."
-#endif
+    ! Compute average over the year for each point
+    var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1))
+#else
+    tsurf_avg_yr1(1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Average of tsurf computed."
 else ! We use slopes
     do islope = 1,nslope
         write(num,'(i2.2)') islope
-        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
+
+        ! CO2 ice
+        !--------
+        call get_var3("co2ice_slope"//num,var3_read_1)
+        where (var3_read_1 < 0.) var3_read_1 = 0.
         write(*,*) "Data for co2_ice_slope"//num//" downloaded."
-        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
+#ifndef CPP_STD
+        call get_var3("perennial_co2ice_slope"//num,var3_read_2)
+        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
+#endif
+
+        ! Compute the minimum over the year for each point
+        var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,1))
+#else
+        min_co2_ice(1,islope,1) = var2_read(1,1)
+#endif
+        write(*,*) "Min of co2_ice computed for slope "//num//"."
+
+        ! H2O ice
+        !--------
+        call get_var3("h2o_ice_s_slope"//num,var3_read_1)
+        where (var3_read_1 < 0.) var3_read_1 = 0.
         write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
-        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
+#ifndef CPP_STD
+        call get_var3("watercap_slope"//num,var3_read_2)
+        write(*,*) "Data for watercap_slope"//num//" downloaded."
+#endif
+
+        ! Compute the minimum over the year for each point
+        var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,1))
+#else
+        min_h2o_ice(1,islope,1) = var2_read(1,1)
+#endif
+        write(*,*) "Min of h2o_ice computed for slope "//num//"."
+
+        ! Tsurf
+        !------
+        call get_var3("tsurf_slope"//num,var3_read_1)
         write(*,*) "Data for tsurf_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
-        write(*,*) "Data for watercap_slope"//num//" downloaded."
-        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
-        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
-#endif
+
+        ! Compute average over the year for each point
+        var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,islope))
+#else
+        tsurf_avg_yr1(1,islope) = var2_read(1,1)
+#endif
+        write(*,*) "Average of tsurf computed for slope "//num//"."
     enddo
 endif
 
 ! Close the NetCDF file
-write(*,*) "Closing "//filename1//"..."
 call error_msg(nf90_close(fID),"close",filename1)
-
-! Compute the minimum over the year for each point
-write(*,*) "Computing the min of h2o_ice..."
-where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
-min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
-write(*,*) "Computing the min of co2_ice..."
-where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
-min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
-
-! Compute averages over the year for each point
-write(*,*) "Computing the average of tsurf..."
-tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
-
-#ifndef CPP_1D
-    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
-    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,1))
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,1))
-    enddo
-#else
-    tsurf_avg_yr1(1,:) = tsurf_avg_dyn(1,1,:)
-    min_co2_ice(1,:,1) = min_co2_ice_dyn(1,1,:)
-    min_h2o_ice(1,:,1) = min_h2o_ice_dyn(1,1,:)
-#endif
-write(*,*) "Downloading data Y1 done!"
-
-!------------------ Year 2 ------------------
-
+write(*,*) '"'//filename1//'" processed!'
+
+!------------------------------- Year 2 --------------------------------
 ! Open the NetCDF file
 write(*,*) "Opening "//filename2//"..."
@@ -152,153 +189,261 @@
 
 ! Download the data from the file
-call get_var3("ps",ps_dyn)
+allocate(var3_read_3(iim_input + 1,jjm_input + 1,nsoilmx),var4_read(iim_input + 1,jjm_input + 1,nsoilmx,timelen))
+
+! Surface pressure
+!-----------------
+call get_var3("ps",var3_read_1)
 write(*,*) "Data for surface pressure downloaded."
 
-call get_var3("co2_layer1",q_co2_dyn)
-write(*,*) "Data for vmr co2 downloaded."
-
-call get_var3("h2o_layer1",q_h2o_dyn)
-write(*,*) "Data for vmr h2o downloaded."
+! Compute average over the year for each point
+var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,ps_timeseries)
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,ps_avg)
+#else
+    ps_timeseries(1,:) = var3_read_1(1,:)
+    ps_avg(1) = var2_read(1,1)
+#endif
+write(*,*) "Average of surface pressure computed."
+
+! CO2 vmr
+!--------
+call get_var3("co2_layer1",var3_read_1)
+where (var3_read_1 < 0.)
+    var3_read_1 = 1.e-10
+else where (var3_read_1 > 1.)
+    var3_read_1 = 1.
+end where
+#ifndef CPP_1D
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_co2)
+#else
+    q_co2(1,:) = var3_read_1(1,1,:)
+#endif
+A = (1./m_co2 - 1./m_noco2)
+B = 1./m_noco2
+vmr_co2_PCM = q_co2/(A*q_co2 + B)/m_co2
+write(*,*) "Data for CO2 vmr downloaded."
+
+! H2O vmr
+!--------
+call get_var3("h2o_layer1",var3_read_1)
+where (var3_read_1 < 0.)
+    var3_read_1 = 1.e-10
+else where (var3_read_1 > 1.)
+    var3_read_1 = 1.
+end where
+#ifndef CPP_1D
+    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var3_read_1,q_h2o)
+#else
+    q_h2o(1,:) = var3_read_1(1,1,:)
+#endif
+write(*,*) "Data for H2O vmr downloaded."
 
 if (nslope == 1) then ! There is no slope
-    call get_var3("co2ice",co2_ice_slope_dyn(:,:,1,:))
+    ! CO2 ice
+    !--------
+    call get_var3("co2ice",var3_read_1)
+    where (var3_read_1 < 0.) var3_read_1 = 0.
     write(*,*) "Data for co2_ice downloaded."
-
-    call get_var3("h2o_ice_s",h2o_ice_s_dyn(:,:,1,:))
+#ifndef CPP_STD
+    call get_var3("perennial_co2ice",var3_read_2)
+    write(*,*) "Data for perennial_co2ice downloaded."
+#endif
+
+    ! Compute the minimum over the year for each point
+    var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,1,1))
+#else
+    min_co2_ice(1,1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of co2_ice computed."
+
+    ! H2O ice
+    !--------
+    call get_var3("h2o_ice_s",var3_read_1)
+    where (var3_read_1 < 0.) var3_read_1 = 0.
     write(*,*) "Data for h2o_ice_s downloaded."
-
-    call get_var3("tsurf",tsurf_dyn(:,:,1,:))
+#ifndef CPP_STD
+    call get_var3("watercap",var3_read_2)
+    write(*,*) "Data for watercap downloaded."
+#endif
+
+    ! Compute the minimum over the year for each point
+    var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,1,1))
+#else
+    min_h2o_ice(1,1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Min of h2o_ice computed."
+
+    ! Tsurf
+    !------
+    call get_var3("tsurf",var3_read_1)
     write(*,*) "Data for tsurf downloaded."
 
-#ifndef CPP_STD
-    call get_var3("watercap",watercap(:,:,1,:))
-    write(*,*) "Data for watercap downloaded."
-
-    call get_var3("perennial_co2ice",perennial_co2ice(:,:,1,:))
-    write(*,*) "Data for perennial_co2ice downloaded."
+    ! Compute average over the year for each point
+    var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg_yr1(:,1))
+#else
+    tsurf_avg_yr1(1,1) = var2_read(1,1)
+#endif
+    write(*,*) "Average of tsurf computed."
 
     if (soil_pem) then
-        call get_var4("soiltemp",tsoil_dyn(:,:,:,1,:))
+        ! Tsoil
+        !------
+        call get_var4("soiltemp",var4_read)
         write(*,*) "Data for soiltemp downloaded."
 
-        call get_var4("waterdensity_soil",watersoil_density_dyn(:,:,:,1,:))
+        ! Compute average over the year for each point
+        var3_read_3 = sum(var4_read,4)/timelen
+#ifndef CPP_1D
+        do l = 1,nsoilmx
+            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,1,:))
+        enddo
+        call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,1))
+#else
+        tsoil_timeseries(1,:,1,:) = var4_read(1,1,:,:)
+        tsoil_avg(1,:,1) = var3_read_3(1,:,1)
+#endif
+        write(*,*) "Average of tsoil computed."
+
+        ! Soil water density
+        !-------------------
+        call get_var4("waterdensity_soil",var4_read)
+#ifndef CPP_1D
+        do l = 1,nsoilmx
+            call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,1,:))
+        enddo
+#else
+        watersoil_density_timeseries(1,:,1,:) = var4_read(1,1,:,:)
+#endif
         write(*,*) "Data for waterdensity_soil downloaded."
 
-        call get_var3("waterdensity_surface",watersurf_density_dyn(:,:,1,:))
+        ! Surface water density
+        !----------------------
+        call get_var3("waterdensity_surface",var3_read_1)
         write(*,*) "Data for waterdensity_surface downloaded."
-    endif !soil_pem
-#endif
+
+        ! Compute average over the year for each point
+        var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,1))
+#else
+        watersurf_density_avg(1,1) = var2_read(1,1)
+#endif
+        write(*,*) "Average of waterdensity_surface computed."
+    endif
 else ! We use slopes
     do islope = 1,nslope
         write(num,'(i2.2)') islope
-        call get_var3("co2ice_slope"//num,co2_ice_slope_dyn(:,:,islope,:))
+
+        ! CO2 ice
+        !--------
+        call get_var3("co2ice_slope"//num,var3_read_1)
+        where (var3_read_1 < 0.) var3_read_1 = 0.
         write(*,*) "Data for co2_ice_slope"//num//" downloaded."
-        call get_var3("h2o_ice_s_slope"//num,h2o_ice_s_dyn(:,:,islope,:))
+#ifndef CPP_STD
+        call get_var3("perennial_co2ice_slope"//num,var3_read_2)
+        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
+#endif
+
+        ! Compute the minimum over the year for each point
+        var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_co2_ice(:,islope,2))
+#else
+        min_co2_ice(1,islope,2) = var2_read(1,1)
+#endif
+        write(*,*) "Min of co2_ice computed for slope "//num//"."
+
+        ! H2O ice
+        !--------
+        call get_var3("h2o_ice_s_slope"//num,var3_read_1)
+        where (var3_read_1 < 0.) var3_read_1 = 0.
         write(*,*) "Data for h2o_ice_s_slope"//num//" downloaded."
-        call get_var3("tsurf_slope"//num,tsurf_dyn(:,:,islope,:))
+#ifndef CPP_STD
+        call get_var3("watercap_slope"//num,var3_read_2)
+        write(*,*) "Data for watercap_slope"//num//" downloaded."
+#endif
+
+        ! Compute the minimum over the year for each point
+        var2_read = minval(var3_read_1 + var3_read_2,3)
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,min_h2o_ice(:,islope,2))
+#else
+        min_h2o_ice(1,islope,2) = var2_read(1,1)
+#endif
+        write(*,*) "Min of h2o_ice computed for slope "//num//"."
+
+        ! Tsurf
+        !------
+        call get_var3("tsurf_slope"//num,var3_read_1)
         write(*,*) "Data for tsurf_slope"//num//" downloaded."
-#ifndef CPP_STD
-        call get_var3("watercap_slope"//num,watercap(:,:,islope,:))
-        write(*,*) "Data for watercap_slope"//num//" downloaded."
-        call get_var3("perennial_co2ice_slope"//num,perennial_co2ice(:,:,islope,:))
-        write(*,*) "Data for perennial_co2ice_slope"//num//" downloaded."
+
+        ! Compute average over the year for each point
+        var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,tsurf_avg(:,islope))
+#else
+        tsurf_avg(1,islope) = var2_read(1,1)
+#endif
+        write(*,*) "Average of tsurf computed for slope "//num//"."
+
         if (soil_pem) then
-            call get_var4("soiltemp_slope"//num,tsoil_dyn(:,:,:,islope,:))
+            ! Tsoil
+            !------
+            call get_var4("soiltemp_slope"//num,var4_read)
             write(*,*) "Data for soiltemp_slope"//num//" downloaded."
-            call get_var4("waterdensity_soil_slope"//num,watersoil_density_dyn(:,:,:,islope,:))
+
+            ! Compute average over the year for each point
+            var3_read_3 = sum(var4_read,4)/timelen
+#ifndef CPP_1D
+            do l = 1,nsoilmx
+                call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),tsoil_timeseries(:,l,islope,:))
+            enddo
+            call gr_dyn_fi(nsoilmx,iim_input + 1,jjm_input + 1,ngrid,var3_read_3,tsoil_avg(:,:,islope))
+#else
+            tsoil_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
+            tsoil_avg(1,:,islope) = var3_read_3(1,:,1)
+#endif
+            write(*,*) "Average of tsoil computed for slope "//num//"."
+
+            ! Soil water density
+            !-------------------
+            call get_var4("waterdensity_soil_slope"//num,var4_read)
+#ifndef CPP_1D
+            do l = 1,nsoilmx
+                call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,var4_read(:,:,l,:),watersoil_density_timeseries(:,l,islope,:))
+            enddo
+#else
+            watersoil_density_timeseries(1,:,islope,:) = var4_read(1,1,:,:)
+#endif
             write(*,*) "Data for waterdensity_soil_slope"//num//" downloaded."
-            call get_var3("waterdensity_surface"//num,watersurf_density_dyn(:,:,islope,:))
-            write(*,*) "Data for waterdensity_surface_slope"//num//" downloaded."
+
+            ! Surface water density
+            !----------------------
+            call get_var3("waterdensity_surface"//num,var3_read_1)
+            write(*,*) "Data for waterdensity_surface"//num//" downloaded."
+
+            ! Compute average over the year for each point
+            var2_read = sum(var3_read_1,3)/timelen
+#ifndef CPP_1D
+            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,var2_read,watersurf_density_avg(:,islope))
+#else
+            watersurf_density_avg(1,islope) = var2_read(1,1)
+#endif
+            write(*,*) "Average of waterdensity_surface computed for slope "//num//"."
         endif
-#endif
     enddo
 endif
+deallocate(var2_read,var3_read_1,var3_read_2,var3_read_3,var4_read)
 
 ! Close the NetCDF file
-write(*,*) "Closing "//filename2//"..."
 call error_msg(nf90_close(fID),"close",filename2)
-
-! Compute the minimum over the year for each point
-write(*,*) "Computing the min of h2o_ice..."
-where (h2o_ice_s_dyn < 0.) h2o_ice_s_dyn = 0.
-min_h2o_ice_dyn = minval(h2o_ice_s_dyn + watercap,4)
-write(*,*) "Computing the min of co2_ice..."
-where (co2_ice_slope_dyn < 0.) co2_ice_slope_dyn = 0.
-min_co2_ice_dyn = minval(co2_ice_slope_dyn + perennial_co2ice,4)
-
-! Compute averages over the year for each point
-write(*,*) "Computing the average of tsurf..."
-tsurf_avg_dyn = sum(tsurf_dyn,4)/timelen
-
-write(*,*) "Computing the average of Ps..."
-ps_avg_dyn = sum(ps_dyn,3)/timelen
-
-if (soil_pem) then
-    write(*,*) "Computing the average of tsoil..."
-    tsoil_avg_dyn = sum(tsoil_dyn,5)/timelen
-    write(*,*) "Computing the average of waterdensity_surface..."
-    watersurf_density_dyn_avg = sum(watersurf_density_dyn,4)/timelen
-endif
-
-! By definition, we get rid of the negative values
-A = (1./m_co2 - 1./m_noco2)
-B = 1./m_noco2
-do i = 1,iim_input + 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-10
-            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_dyn(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_dyn,vmr_co2_PCM)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_dyn,ps_timeseries)
-    call gr_dyn_fi(timelen,iim_input + 1,jjm_input + 1,ngrid,ps_avg_dyn,ps_avg)
-    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,2))
-        call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,min_h2o_ice_dyn(:,:,islope),min_h2o_ice(:,islope,2))
-        if (soil_pem) then
-            do l = 1,nsoilmx
-                call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_avg_dyn(:,:,l,islope),tsoil_avg(:,l,islope))
-                do t = 1,timelen
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsoil_dyn(:,:,l,islope,t),tsoil_timeseries(:,l,islope,t))
-                    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersoil_density_dyn(:,:,l,islope,t),watersoil_density_timeseries(:,l,islope,t))
-                enddo
-            enddo
-            call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,watersurf_density_dyn_avg(:,:,islope),watersurf_density_avg(:,islope))
-        endif ! soil_pem
-    enddo
-    call gr_dyn_fi(nslope,iim_input + 1,jjm_input + 1,ngrid,tsurf_avg_dyn,tsurf_avg)
-#else
-    vmr_co2_PCM(1,:) = vmr_co2_dyn(1,1,:)
-    ps_timeseries(1,:) = ps_dyn(1,1,:)
-    ps_avg(1) = ps_avg_dyn(1,1)
-    q_co2(1,:) = q_co2_dyn(1,1,:)
-    q_h2o(1,:) = q_h2o_dyn(1,1,:)
-    min_co2_ice(1,:,2) = min_co2_ice_dyn(1,1,:)
-    min_h2o_ice(1,:,2) = min_h2o_ice_dyn(1,1,:)
-    if (soil_pem) then
-        tsoil_timeseries(1,:,:,:) = tsoil_dyn(1,1,:,:,:)
-        tsoil_avg(1,:,:) = tsoil_avg_dyn(1,1,:,:)
-        watersoil_density_timeseries(1,:,:,:) = watersoil_density_dyn(1,1,:,:,:)
-        watersurf_density_avg(1,:) = watersurf_density_dyn_avg(1,1,:)
-    endif ! soil_pem
-    tsurf_avg(1,:) = tsurf_avg_dyn(1,1,:)
-#endif
-write(*,*) "Downloading data Y2 done!"
+write(*,*) '"'//filename2//'" processed!'
 
 END SUBROUTINE read_data_PCM
@@ -392,2 +537,3 @@
 
 END MODULE read_data_PCM_mod
+
