Index: unk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/compute_tendencies_mod.F90	(revision 2892)
+++ 	(revision )
@@ -1,67 +1,0 @@
-!
-! $Id $
-!
-SUBROUTINE compute_tendencies(tendencies_h2o_ice,min_h2o_ice_Y1,&
-     min_h2o_ice_Y2,iim_input,jjm_input,ngrid,tendencies_h2o_ice_phys)
-
-      IMPLICIT NONE
-
-
-!=======================================================================
-!
-!  Compute the tendencies of the evolution of water ice over the years
-!
-!=======================================================================
-
-!   arguments:
-!   ----------
-
-!   INPUT
-
-     INTEGER, intent(in) :: iim_input,jjm_input,ngrid                             ! # of grid points along longitude/latitude/ total
-     REAL, intent(in) , dimension(iim_input+1,jjm_input+1):: min_h2o_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
-     REAL, intent(in) , dimension(iim_input+1,jjm_input+1):: min_h2o_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
-
-!   OUTPUT
-     REAL, intent(out) , dimension(iim_input+1,jjm_input+1) :: tendencies_h2o_ice ! LON x LAT field : difference between the minima = evolution of perenial ice
-     REAL, intent(out) , dimension(ngrid)   :: tendencies_h2o_ice_phys            ! physical point field : difference between the minima = evolution of perenial ice
-
-!   local:
-!   ------
-
-     INTEGER :: i,j,ig0                                                           ! loop variable
-
-!=======================================================================
-
-
-!  We compute the difference
-  tendencies_h2o_ice(:,:)=min_h2o_ice_Y2(:,:)-min_h2o_ice_Y1(:,:)
-
-!  If the difference is too small; there is no evolution
-  DO j=1,jjm_input+1
-    DO i = 1, iim_input
-       if(abs(tendencies_h2o_ice(i,j)).LT.1.0E-10) then
-          tendencies_h2o_ice(i,j)=0.
-       endif
-    ENDDO
-  ENDDO
-
-!  We reorganise the difference on the physical grid
-  tendencies_h2o_ice_phys(1)=tendencies_h2o_ice(1,1)
-
-  ig0 = 2
-  DO j=2,jjm_input
-    DO i = 1, iim_input
-       tendencies_h2o_ice_phys(ig0)  =tendencies_h2o_ice(i,j)
-       ig0= ig0 + 1
-    ENDDO
-  ENDDO
-
-  tendencies_h2o_ice_phys(ig0) = tendencies_h2o_ice(1,jjm_input+1)
-
-END SUBROUTINE compute_tendencies
-
-
-
-
-
Index: /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/comsoil_h_PEM.F90	(revision 2893)
@@ -42,5 +42,4 @@
     allocate(inertiedat_PEM(ngrid,nsoilmx_PEM)) 
     allocate(water_reservoir(ngrid))
-    allocate(water_reservoir(ngrid))
   end subroutine ini_comsoil_h_PEM
 
Index: /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/conf_pem.F90	(revision 2893)
@@ -14,10 +14,12 @@
 #endif
   
-  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, ice_criterion, ps_criterion, &
+  USE temps_mod_evol, ONLY: year_bp_ini, dt_pem, water_ice_criterion, co2_ice_criterion, ps_criterion, &
                 Max_iter_pem, evol_orbit_pem
   USE comsoil_h_pem, only: soil_pem,fluxgeo,water_reservoir_nom
   USE adsorption_mod,only: adsorption_pem
   CHARACTER(len=20),parameter :: modname ='conf_pem' 
+
 !PEM parameter
+
   year_bp_ini=0.
   CALL getin('year_bp_ini', year_bp_ini)
@@ -26,6 +28,9 @@
   CALL getin('dt_pem', dt_pem)
 
-  ice_criterion=0.2
-  CALL getin('ice_criterion', ice_criterion)
+  water_ice_criterion=0.2
+  CALL getin('water_ice_criterion', water_ice_criterion)
+
+  co2_ice_criterion=0.2
+  CALL getin('co2_ice_criterion', co2_ice_criterion)
 
   ps_criterion = 0.15
@@ -51,13 +56,11 @@
        print*,'Adsorption must be used when soil_pem = T'
        call abort_physic(modname,"Adsorption must be used when soil_pem = T",1)
-        endif
+  endif
   
-
   if ((not(soil_pem)).and.fluxgeo.gt.0.) then
        print*,'Soil is not activated but Flux Geo > 0.'
        call abort_physic(modname,"Soil is not activated but Flux Geo > 0.'",1)
-        endif
+  endif
   
-
    water_reservoir_nom = 1e4
   CALL getin('water_reservoir_nom',water_reservoir_nom)
Index: /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/criterion_pem_stop_mod.F90	(revision 2893)
@@ -15,5 +15,5 @@
 SUBROUTINE criterion_waterice_stop(cell_area,ini_surf,qsurf,STOPPING,ngrid,initial_h2o_ice)
 
-  USE temps_mod_evol, ONLY: ice_criterion
+  USE temps_mod_evol, ONLY: water_ice_criterion
   use comslope_mod, ONLY: subslope_dist,nslope
 
@@ -32,5 +32,5 @@
   INTEGER, intent(in) :: ngrid                  ! # of grid physical grid points 
   REAL,    intent(in) :: cell_area(ngrid)       ! physical point field : Area of the cells
-  REAL,    intent(in) ::  qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
+  REAL,    intent(in) :: qsurf(ngrid,nslope)          ! physical point field : Actual density of water ice
   REAL,    intent(in) :: ini_surf
   REAL,    intent(in) :: initial_h2o_ice(ngrid,nslope)
@@ -60,12 +60,12 @@
   
 !   check of the criterion
-  if(present_surf.LT.ini_surf*(1-ice_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+ice_criterion)) then
+  if(present_surf.LT.ini_surf*(1-water_ice_criterion) .OR. &
+     present_surf.GT.ini_surf*(1+water_ice_criterion)) then
     STOPPING=.TRUE. 
     print *, "Reason of stopping : The surface of water ice sublimating reach the threshold:"
     print *, "Current surface of water ice sublimating=", present_surf
     print *, "Initial surface of water ice sublimating=", ini_surf
-    print *, "Percentage of change accepted=", ice_criterion*100
-    print *, "present_surf<ini_surf*(1-ice_criterion)", (present_surf.LT.ini_surf*(1-ice_criterion))
+    print *, "Percentage of change accepted=", water_ice_criterion*100
+    print *, "present_surf<ini_surf*(1-water_ice_criterion)", (present_surf.LT.ini_surf*(1-water_ice_criterion))
   endif
 
@@ -80,5 +80,5 @@
 SUBROUTINE criterion_co2_stop(cell_area,ini_surf,qsurf,STOPPING_ice,STOPPING_ps,ngrid,initial_h2o_ice,global_ave_press_GCM,global_ave_press_new,nslope)
 
-  USE temps_mod_evol, ONLY: ice_criterion,ps_criterion
+  USE temps_mod_evol, ONLY: co2_ice_criterion,ps_criterion
   use comslope_mod, ONLY: subslope_dist
 
@@ -127,12 +127,12 @@
   
 !   check of the criterion
-  if(present_surf.LT.ini_surf*(1-ice_criterion) .OR. &
-     present_surf.GT.ini_surf*(1+ice_criterion)) then
+  if(present_surf.LT.ini_surf*(1-co2_ice_criterion) .OR. &
+     present_surf.GT.ini_surf*(1+co2_ice_criterion)) then
     STOPPING_ice=.TRUE. 
     print *, "Reason of stopping : The surface of co2 ice sublimating reach the threshold:"
     print *, "Current surface of co2 ice sublimating=", present_surf
     print *, "Initial surface of co2 ice sublimating=", ini_surf
-    print *, "Percentage of change accepted=", ice_criterion*100
-    print *, "present_surf<ini_surf*(1-ice_criterion)", (present_surf.LT.ini_surf*(1-ice_criterion))
+    print *, "Percentage of change accepted=", co2_ice_criterion*100
+    print *, "present_surf<ini_surf*(1-co2_ice_criterion)", (present_surf.LT.ini_surf*(1-co2_ice_criterion))
   endif
 
@@ -155,8 +155,2 @@
 
 END MODULE
-
-
-
-
-
-
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 2893)
@@ -143,4 +143,5 @@
 
       character (len = *), parameter :: FILE_NAME = "startfi_evol.nc" !Name of the file used for initialsing the PEM
+      character*2 str2
       integer :: ncid, varid,status                                !Variable for handling opening of files
       integer :: phydimid, subdimid, nlayerdimid, nqdimid          !Variable ID for Netcdf files
@@ -158,23 +159,6 @@
 
 ! Variable for h2o_ice evolution
-
-      REAL , dimension(:,:), allocatable :: tendencies_h2o_ice     ! LON x LAT field : Tendency of evolution of perenial ice over a year
-      REAL, dimension(:),allocatable  :: tendencies_h2o_ice_phys   ! physical point field : Tendency of evolution of perenial ice over a year
-
-      REAL , dimension(:,:), allocatable :: tendencies_co2_ice     ! LON x LAT field : Tendency of evolution of perenial co2 ice over a year
-      REAL, dimension(:),allocatable  :: tendencies_co2_ice_phys   ! physical point field : Tendency of evolution of perenial co2 ice over a year
-
-      REAL :: ini_surf                                             ! Initial surface of sublimating water ice
-      REAL :: ini_surf_h2o                                             ! Initial surface of sublimating water ice
-      REAL, dimension(:),allocatable  :: initial_h2o_ice           ! physical point field : Logical array indicating sublimating point
-
+      REAL :: ini_surf_h2o                                         ! Initial surface of sublimating water ice
       REAL :: ini_surf_co2                                         ! Initial surface of sublimating co2 ice
-      REAL, dimension(:),allocatable  :: initial_co2_ice           ! physical point field : Logical array indicating sublimating point of co2 ice
-
-      REAL , dimension(:,:), allocatable :: min_h2o_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
-      REAL , dimension(:,:), allocatable :: min_h2o_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
-
-      REAL , dimension(:,:), allocatable :: min_co2_ice_s_1        ! LON x LAT field : minimum of water ice at each point for the first year
-      REAL , dimension(:,:), allocatable :: min_co2_ice_s_2        ! LON x LAT field : minimum of water ice at each point for the second year
 
       REAL :: global_ave_press_GCM           ! constant: global average pressure retrieved in the GCM [Pa]
@@ -193,5 +177,4 @@
       LOGICAL :: STOPPING_pressure
       INTEGER :: criterion_stop  ! which criterion is reached ? 1= h2o ice surf, 2 = co2 ice surf, 3 = ps, 4 = orb param
-
 
       REAL, dimension(:,:,:),allocatable  :: q_co2_GCM ! Lon x Lat x Time : mass mixing ratio of co2 in the first layer [kg/kg]
@@ -430,5 +413,4 @@
    
      call surfini(ngrid,qsurf)
-     call surfini(ngrid,qsurf)
 
 #else
@@ -452,8 +434,8 @@
 
 if (nslope.eq.1) then
-def_slope(1) = 0
-def_slope(2) = 0
-def_slope_mean=0
-subslope_dist(:,1) = 1.
+ def_slope(1) = 0
+ def_slope(2) = 0
+ def_slope_mean=0
+ subslope_dist(:,1) = 1.
 endif
 
@@ -500,6 +482,6 @@
      allocate(flag_co2flow_mesh(ngrid))
 
-       flag_co2flow(:,:) = 0.     
-       flag_co2flow_mesh(:) = 0.
+     flag_co2flow(:,:) = 0.     
+     flag_co2flow_mesh(:) = 0.
 
 
@@ -518,6 +500,4 @@
      call nb_time_step_GCM("data_GCM_Y1.nc",timelen)
 
-     allocate(min_h2o_ice_s_1(iim+1,jjm+1))
-     allocate(min_co2_ice_s_1(iim+1,jjm+1))
      allocate(vmr_co2_gcm(iim+1,jjm+1,timelen))
      allocate(q_h2o_GCM(iim+1,jjm+1,timelen))
@@ -552,5 +532,5 @@
      print *, "Downloading data Y1..."
 
-     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm, min_h2o_ice_s_1,min_co2_ice_s_1,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
+     call read_data_GCM("data_GCM_Y1.nc",timelen, iim,jjm,vmr_co2_gcm,ps_GCM_yr1,min_co2_ice_slope_1,min_h2o_ice_slope_1,&   
                        nslope,tsurf_ave_yr1,tsoil_ave_yr1, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
                        watersurf_density_timeseries,watersoil_density_timeseries)
@@ -560,6 +540,4 @@
      print *, "Downloading data Y1 done"
 
-     allocate(min_h2o_ice_s_2(iim+1,jjm+1))
-     allocate(min_co2_ice_s_2(iim+1,jjm+1))
      allocate(min_co2_ice_slope_2(iim+1,jjm+1,nslope))
      allocate(min_h2o_ice_slope_2(iim+1,jjm+1,nslope))
@@ -567,5 +545,5 @@
      print *, "Downloading data Y2"
 
-     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm ,min_h2o_ice_s_2,min_co2_ice_s_2,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
+     call read_data_GCM("data_GCM_Y2.nc",timelen,iim,jjm,vmr_co2_gcm,ps_GCM,min_co2_ice_slope_2,min_h2o_ice_slope_2, &
                   nslope,tsurf_ave,tsoil_ave, tsurf_GCM_timeseries,tsoil_GCM_timeseries,TI_GCM,q_co2_GCM,q_h2o_GCM,co2_ice_GCM_slope, &      
                   watersurf_density_timeseries,watersoil_density_timeseries)
@@ -639,5 +617,4 @@
       deallocate(TI_GCM)
       deallocate(tsurf_GCM_timeseries)
-
 
    if(soil_pem) then
@@ -654,5 +631,4 @@
            CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,watersoil_density_timeseries(:,:,l,islope,t),watersoil_density_phys_PEM_timeseries(:,l,islope,t))
            ENDDO
-
         ENDDO
         DO l=nsoilmx+1,nsoilmx_PEM
@@ -664,5 +640,5 @@
         ENDDO
       ENDDO
-       watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
+      watersoil_density_phys_PEM_ave(:,:,:) = SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
       watersurf_density_phys_ave(:,:) = SUM(watersurf_density_phys_timeseries(:,:,:),3)/timelen
       deallocate(watersurf_density_timeseries)
@@ -685,9 +661,4 @@
 
 !----- Compute tendencies from the PCM run
-
-     allocate(tendencies_h2o_ice(iim+1,jjm+1))
-     allocate(tendencies_h2o_ice_phys(ngrid))
-     allocate(tendencies_co2_ice(iim+1,jjm+1))
-     allocate(tendencies_co2_ice_phys(ngrid))
      allocate(tendencies_co2_ice_slope(iim+1,jjm+1,nslope))
      allocate(tendencies_co2_ice_phys_slope(ngrid,nslope))
@@ -697,10 +668,4 @@
 
 !  Compute the tendencies of the evolution of ice over the years
-
-      call compute_tendencies(tendencies_h2o_ice,min_h2o_ice_s_1,&
-             min_h2o_ice_s_2,iim,jjm,ngrid,tendencies_h2o_ice_phys)
-
-      call compute_tendencies(tendencies_co2_ice,min_co2_ice_s_1,&
-             min_co2_ice_s_2,iim,jjm,ngrid,tendencies_co2_ice_phys)
 
       call compute_tendencies_slope(tendencies_co2_ice_slope,min_co2_ice_slope_1,&
@@ -725,6 +690,4 @@
 !---------------------------- Save initial PCM situation --------------------- 
 
-     allocate(initial_h2o_ice(ngrid))
-     allocate(initial_co2_ice(ngrid))
      allocate(initial_co2_ice_sublim_slope(ngrid,nslope))
      allocate(initial_co2_ice_slope(ngrid,nslope))
@@ -733,5 +696,4 @@
 ! We save the places where water ice is sublimating
 ! We compute the surface of water ice sublimating
-  ini_surf=0.
   ini_surf_co2=0.
   ini_surf_h2o=0.
@@ -739,10 +701,4 @@
   do i=1,ngrid
     Total_surface=Total_surface+cell_area(i)
-    if (tendencies_h2o_ice_phys(i).LT.0) then
-         initial_h2o_ice(i)=1.
-         ini_surf=ini_surf+cell_area(i)
-    else
-         initial_h2o_ice(i)=0.         
-    endif
     do islope=1,nslope
       if (tendencies_co2_ice_phys_slope(i,islope).LT.0) then
@@ -767,5 +723,4 @@
 
      print *, "Total initial surface of co2ice sublimating (slope)=", ini_surf_co2
-     print *, "Total initial surface of h2o ice sublimating=", ini_surf
      print *, "Total initial surface of h2o ice sublimating (slope)=", ini_surf_h2o
      print *, "Total surface of the planet=", Total_surface
@@ -813,6 +768,4 @@
       enddo
     enddo
-
-
 
     if(adsorption_pem) then
@@ -889,8 +842,10 @@
       do i=1,ngrid
         global_ave_press_new = global_ave_press_new -g*cell_area(i)*delta_co2_adsorbded(i)/Total_surface
-       enddo 
+      enddo 
+      print *, 'Global average pressure old time step',global_ave_press_old
+      print *, 'Global average pressure new time step',global_ave_press_new
      endif
-       print *, 'Global average pressure old time step',global_ave_press_old
-       print *, 'Global average pressure new time step',global_ave_press_new
+
+     call WRITEDIAGFI(ngrid,'ps_ave','Global average pressure','Pa',0,global_ave_press_new)
 
 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
@@ -961,4 +916,9 @@
      call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
 
+     DO islope=1, nslope
+       write(str2(1:2),'(i2.2)') islope
+       call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,'H2O ice','kg.m-2',2,qsurf_slope(:,igcm_h2o_ice,islope))
+       call WRITEDIAGFI(ngrid,'tendencies_h2o_ice_slope'//str2,'H2O ice tend','kg.m-2.year-1',2,tendencies_h2o_ice_phys_slope(:,islope))
+     ENDDO
 
       print *, "Evolution of co2 ice"
@@ -976,12 +936,12 @@
       print *, "Co2 glacier flows"
 
-
-
     call co2glaciers_evol(timelen,ngrid,nslope,iflat,subslope_dist,def_slope_mean,vmr_co2_pem_phys,ps_phys_timeseries,&
                          global_ave_press_GCM,global_ave_press_new,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
 
-
-
-
+     DO islope=1, nslope
+       write(str2(1:2),'(i2.2)') islope
+       call WRITEDIAGFI(ngrid,'co2ice_slope'//str2,'CO2 ice','kg.m-2',2,co2ice_slope(:,islope))
+       call WRITEDIAGFI(ngrid,'tendencies_co2_ice_slope'//str2,'CO2 ice tend','kg.m-2.year-1',2,tendencies_co2_ice_phys_slope(:,islope))
+     ENDDO
 
 !------------------------
@@ -1063,4 +1023,8 @@
       enddo
 
+     DO islope=1, nslope
+       write(str2(1:2),'(i2.2)') islope
+       call WRITEDIAGFI(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf_slope(:,islope))
+     ENDDO
 
     if(soil_pem) then
@@ -1090,11 +1054,8 @@
          enddo
        enddo
-
-
      enddo
   enddo
      tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
      watersoil_density_phys_PEM_ave(:,:,:)= SUM(watersoil_density_phys_PEM_timeseries(:,:,:,:),4)/timelen
-
 
   print *, "Update of soil temperature done"
@@ -1134,5 +1095,5 @@
 
       print *, "Adaptation of the new co2 tendencies to the current pressure"
-      call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
+      call recomp_tend_co2_slope(tendencies_co2_ice_phys_slope,tendencies_co2_ice_phys_slope_ini,co2ice_slope,vmr_co2_gcm_phys,vmr_co2_pem_phys,ps_phys_timeseries,&     
                                global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
 
@@ -1154,6 +1115,4 @@
 
       year_iter=year_iter+dt_pem
-
-      
 
       print *, "Checking all the stopping criterion."
@@ -1191,6 +1150,4 @@
       endif
 
-
-
       global_ave_press_old=global_ave_press_new
 
@@ -1223,12 +1180,7 @@
       ENDDO ! of DO ig=1,ngrid
 
-
 ! H2O ice
-
-
-
      DO ig=1,ngrid
        if(watercaptag(ig)) then
-        
          watercap_sum=0.
          DO islope=1,nslope
@@ -1247,5 +1199,4 @@
      enddo
 
-
       DO ig = 1,ngrid
         qsurf(ig,igcm_h2o_ice) = 0.
@@ -1258,13 +1209,11 @@
 
      DO ig=1,ngrid
-!       DO islope=1,nslope
-         if(qsurf(ig,igcm_h2o_ice).GT.500) then
+       DO islope=1,nslope
+         if(qsurf_slope(ig,igcm_h2o_ice,islope).GT.500) then
             watercaptag(ig)=.true.
-       DO islope=1,nslope
             qsurf_slope(ig,igcm_h2o_ice,islope)=qsurf_slope(ig,igcm_h2o_ice,islope)-250
             water_reservoir(ig)=water_reservoir(ig)+250
-       ENDDO
          endif
-!       enddo
+       enddo
      enddo
 
@@ -1277,5 +1226,4 @@
        ENDDO
          endif
-
      enddo
 
@@ -1289,11 +1237,4 @@
       ENDDO ! of DO ig=1,ngrid
 
-
-
-
-
-      
-
-
 ! III_a.2  Tsoil update (for startfi)
 
@@ -1304,5 +1245,4 @@
       TI_GCM_phys(:,:,:)=TI_GCM_start(:,:,:)
    endif !soil_pem
-
 
 #ifndef CPP_STD
@@ -1339,5 +1279,4 @@
       ENDDO
     ENDDO
-
 
       DO ig = 1,ngrid
@@ -1410,5 +1349,4 @@
      DO ig=1,ngrid
        if(watercaptag(ig)) then
-         print *, "INNN"
          WC_sum=0.
          DO islope=1,nslope
Index: /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 2893)
@@ -91,5 +91,5 @@
 write(*,*)'Is start PEM?',startpem_file
 
-startpem_file = .false.
+startpem_file = .true.
 !1. Run
 
@@ -317,7 +317,7 @@
       do ig=1,ngrid
         if(watercaptag(ig)) then
-           water_reservoir=water_reservoir_nom
+           water_reservoir(ig)=water_reservoir_nom
         else
-           water_reservoir=0.
+           water_reservoir(ig)=0.
         endif
       enddo
@@ -483,10 +483,4 @@
 endif !soil_pem
 
-endif ! of if (startphy_file)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-
-
 !. e) water reservoir
 
@@ -502,4 +496,8 @@
 
 #endif
+
+endif ! of if (startphy_file)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
 if(soil_pem) then
Index: /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/read_data_GCM.F90	(revision 2893)
@@ -2,5 +2,5 @@
 ! $Id $
 !
-SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,min_h2o_ice_s,min_co2_ice_s,vmr_co2_gcm,ps_GCM, & 
+SUBROUTINE read_data_GCM(fichnom,timelen, iim_input,jjm_input,vmr_co2_gcm,ps_GCM, & 
              min_co2_ice_slope,min_h2o_ice_slope,nslope,tsurf_ave,tsoil_ave,tsurf_gcm,tsoil_gcm,TI_ave,q_co2_GCM,q_h2o_GCM,co2_ice_slope, &
              watersurf_density,watersoil_density)
@@ -31,6 +31,4 @@
 
 ! Ouputs
-  REAL, INTENT(OUT) ::  min_h2o_ice_s(iim_input+1,jjm_input+1) ! Minimum of h2o ice, mesh averaged of the year  [kg/m^2]
-  REAL, INTENT(OUT) ::  min_co2_ice_s(iim_input+1,jjm_input+1) ! Minimum of co2 ice, mesh averaged  of the year [kg/m^2]
   REAL, INTENT(OUT) ::  min_co2_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of co2 ice  per slope of the year [kg/m^2]
   REAL, INTENT(OUT) ::  min_h2o_ice_slope(iim_input+1,jjm_input+1,nslope) ! Minimum of h2o ice per slope of the year [kg/m^2]
@@ -78,8 +76,5 @@
       B=1/m_noco2
 
-      allocate(co2_ice_s(iim+1,jjm+1,timelen))
       allocate(h2o_ice_s_slope(iim+1,jjm+1,nslope,timelen))
-      allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
-      allocate(h2o_ice_s(iim+1,jjm+1,timelen))
       allocate(watercap_slope(iim_input+1,jjm_input+1,nslope,timelen))
 
@@ -90,15 +85,4 @@
   CALL err(NF90_OPEN(var,NF90_NOWRITE,fID),"open",var)
 
-     print *, "Downloading data for h2oice ..."
-
-! Get h2o_ice_s of the concatenated file
-  CALL get_var3("h2o_ice_s"   ,h2o_ice_s)
-
-     print *, "Downloading data for h2oice done"
-     print *, "Downloading data for co2ice ..."
-
-  CALL get_var3("co2ice"   ,co2_ice_s)
-
-     print *, "Downloading data for co2ice done"
      print *, "Downloading data for vmr co2..."
 
@@ -202,9 +186,4 @@
 
 ! Compute the minimum over the year for each point
-  print *, "Computing the min of h2o_ice"
-  min_h2o_ice_s(:,:)=minval(h2o_ice_s,3)
-  print *, "Computing the min of co2_ice"
-  min_co2_ice_s(:,:)=minval(co2_ice_s,3)
-
   print *, "Computing the min of h2o_ice_slope"
   min_h2o_ice_slope(:,:,:)=minval(h2o_ice_s_slope+watercap_slope,4)
@@ -227,8 +206,4 @@
   DO i=1,iim+1
     DO j = 1, jjm+1
-       if (min_co2_ice_s(i,j).LT.0) then
-          min_h2o_ice_s(i,j)  = 0.
-          min_co2_ice_s(i,j)  = 0.
-       endif
        DO islope=1,nslope
           if (min_co2_ice_slope(i,j,islope).LT.0) then
@@ -261,9 +236,5 @@
   ENDDO
 
-
-      deallocate(co2_ice_s)
       deallocate(h2o_ice_s_slope)
-      deallocate(h2o_ice_s)
-
 
   CONTAINS
Index: /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope.F90	(revision 2893)
@@ -2,5 +2,5 @@
 ! $Id $
 !
-SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
+SUBROUTINE recomp_tend_co2_slope(tendencies_co2_ice_phys,tendencies_co2_ice_phys_ini,co2ice_slope,vmr_co2_gcm,vmr_co2_pem,ps_GCM_2,global_ave_press_GCM,global_ave_press_new,timelen,ngrid,nslope)
 
       IMPLICIT NONE
@@ -23,4 +23,5 @@
   REAL, intent(in) :: global_ave_press_new
   REAL, intent(in) ::  tendencies_co2_ice_phys_ini(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
+  REAL, intent(in) :: co2ice_slope(ngrid,nslope)
 
 !   OUTPUT
@@ -39,5 +40,5 @@
   alpha=23.3494
 
-  coef=669*24*3600*eps*sigma/L
+  coef=669*88875*eps*sigma/L
 
 ! Evolution of the water ice for each physical point
@@ -45,8 +46,9 @@
     do islope=1,nslope
       ave=0.
-      if(abs(tendencies_co2_ice_phys(i,islope)).gt.1e-4) then
+!      if(abs(tendencies_co2_ice_phys(i,islope)).gt.1e-4) then
+      if(co2ice_slope(i,islope).gt.1e-4 .and. abs(tendencies_co2_ice_phys(i,islope)).gt.1e-5) then
         do t=1,timelen
            ave=ave+(beta/(alpha-log(vmr_co2_gcm(i,t)*ps_GCM_2(i,t)/100.)))**4  &
-              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*global_ave_press_GCM/global_ave_press_new/100.)))**4
+              -(beta/(alpha-log(vmr_co2_pem(i,t)*ps_GCM_2(i,t)*(global_ave_press_new/global_ave_press_GCM)/100.)))**4
          enddo
       endif
Index: /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2892)
+++ /trunk/LMDZ.COMMON/libf/evolution/temps_mod_evol.F90	(revision 2893)
@@ -3,11 +3,12 @@
 IMPLICIT NONE  
 
-  INTEGER   year_bp_ini     !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
-  INTEGER   dt_pem          !     dt_pem  : in years, the time step used by the PEM    
-  REAL      ice_criterion   !     ice_criterion : percentage of change of ice before stopping the PEM 
-  REAL      ps_criterion    !     ice_criterion : percentage of change of ice before stopping the PEM 
-  INTEGER   year_PEM        !     year written in startfiPEM.nc
-  INTEGER   Max_iter_pem    !     Maximal number of iteration when converging to a steady state, read in evol.def
-  LOGICAL   evol_orbit_pem  !     True if we want to follow the orbital parameters of ob_ex_lsp.asc, read in evol.def
+  INTEGER   year_bp_ini         !     year_bp_ini : Initial year of the simulation of the PEM (in evol.def)
+  INTEGER   dt_pem              !     dt_pem  : in Planetary years, the time step used by the PEM    
+  REAL      water_ice_criterion !     Percentage of change of the surface of water ice sublimating before stopping the PEM 
+  REAL      co2_ice_criterion   !     Percentage of change of the surface of co2 ice sublimating before stopping the PEM 
+  REAL      ps_criterion        !     Percentage of change of averaged surface pressure before stopping the PEM 
+  INTEGER   year_PEM            !     year written in startfiPEM.nc
+  INTEGER   Max_iter_pem        !     Maximal number of iteration when converging to a steady state, read in evol.def
+  LOGICAL   evol_orbit_pem      !     True if we want to follow the orbital parameters of ob_ex_lsp.asc, read in evol.def
 
 END MODULE temps_mod_evol
Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2892)
+++ /trunk/LMDZ.MARS/README	(revision 2893)
@@ -3860,2 +3860,10 @@
 Change back the modification of r2883 to its original form + Retrocompatibility for icelocation_mode.ne.5 in surfini
 
+== 10/02/2023 == RV
+Correction of recomp_tend_co2 formula and make it dependent on the presence of co2_ice.
+Add different stopping criterion for water_ice and co2_ice called : water_ice_criterion and co2_ice_criterion.
+Add the possibility to output a diagfi.nc each year for the amount of ice, the tendencies, tsurf and ps.
+Remove useless variables (not sloped) 
+Remove useless file
+Some cleaning
+