Index: trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3019)
+++ trunk/LMDZ.COMMON/libf/evolution/pemetat0.F90	(revision 3039)
@@ -1,31 +1,27 @@
-   SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
+SUBROUTINE pemetat0(filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,ice_table, ice_table_thickness, &
                       tsurf_ave_yr1,tsurf_ave_yr2,q_co2,q_h2o,ps_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
                       global_ave_pressure,watersurf_ave,watersoil_ave, m_co2_regolith_phys,deltam_co2_regolith_phys,  &
                       m_h2o_regolith_phys,deltam_h2o_regolith_phys, water_reservoir)
    
-   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
-   use comsoil_h_PEM, only: soil_pem,layer_PEM, mlayer_PEM,fluxgeo,inertiedat_PEM,water_reservoir_nom,depth_breccia,depth_bedrock,index_breccia,index_bedrock
-   use comsoil_h, only:  volcapa,inertiedat 
-   use adsorption_mod, only : regolith_adsorption,adsorption_pem
-   USE temps_mod_evol, ONLY: year_PEM
-   USE ice_table_mod, only: computeice_table_equilibrium
-   USE constants_marspem_mod,only: alpha_clap_h2o,beta_clap_h2o, &
-                                   TI_breccia,TI_bedrock
-   use soil_thermalproperties_mod, only: update_soil_thermalproperties
-   use tracer_mod, only: mmol,igcm_h2o_vap ! tracer names and molar masses
+  use iostart_PEM,                only: open_startphy, close_startphy, get_field, get_var
+  use comsoil_h_PEM,              only: soil_pem, layer_PEM, mlayer_PEM, fluxgeo, inertiedat_PEM, water_reservoir_nom, depth_breccia, depth_bedrock, index_breccia, index_bedrock
+  use comsoil_h,                  only: volcapa,inertiedat 
+  use adsorption_mod,             only: regolith_adsorption, adsorption_pem
+  use ice_table_mod,              only: computeice_table_equilibrium
+  use constants_marspem_mod,      only: alpha_clap_h2o, beta_clap_h2o, TI_breccia, TI_bedrock
+  use soil_thermalproperties_mod, only: update_soil_thermalproperties
+  use tracer_mod,                 only: mmol, igcm_h2o_vap ! tracer names and molar masses
 
 #ifndef CPP_STD
-      USE comcstfi_h, only: r, mugaz
+    use comcstfi_h,   only: r, mugaz
+    use surfdat_h,    only: watercaptag
 #else
-      USE comcstfi_mod, only: r, mugaz
+    use comcstfi_mod, only: r, mugaz
 #endif
 
-
-#ifndef CPP_STD   
-   use surfdat_h, only: watercaptag
-#endif
-
   implicit none
- 
+
+  include "callkeys.h"
+
   character(len=*), intent(in) :: filename                    ! name of the startfi_PEM.nc
   integer,intent(in) :: ngrid                                 ! # of physical grid points
@@ -73,5 +69,4 @@
    real :: alph_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computation []
    real :: beta_tmp(ngrid,nsoil_PEM-1)                                ! Intermediate for tsoil computatio []                               
-   real :: year_PEM_read                                              ! Year of the PEM previous run
    LOGICAL :: startpem_file                                           ! boolean to check if we read the startfile or not
 #ifdef CPP_STD   
@@ -97,5 +92,4 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-
 !0. Check if the start_PEM exist.
 
@@ -110,13 +104,5 @@
    call open_startphy(filename)
 
-   call get_var("Time",year_PEM_read,found)
-   year_PEM=INT(year_PEM_read)
-   if(.not.found) then
-      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
-   else
-      write(*,*)'year_PEM of startpem=', year_PEM
-   endif
-
-    if(soil_pem) then
+    if (soil_pem) then
  
 !1. Thermal Inertia
@@ -161,5 +147,5 @@
 ENDDO ! islope
 
-print *,'PEMETAT0: THERMAL INERTIA DONE'
+write(*,*) 'PEMETAT0: THERMAL INERTIA DONE'
 
 ! b. Special case for inertiedat, inertiedat_PEM
@@ -204,7 +190,5 @@
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
 !2. Soil Temperature
-
 DO islope=1,nslope
   write(num,fmt='(i2.2)') islope
@@ -260,15 +244,10 @@
         enddo
       enddo
-
-    
 ENDDO ! islope
 
-print *,'PEMETAT0: SOIL TEMP  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+write(*,*) 'PEMETAT0: SOIL TEMP  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !3. Ice Table
-
-
   call get_field("ice_table",ice_table,found)
    if(.not.found) then
@@ -282,10 +261,8 @@
    endif
 
-print *,'PEMETAT0: ICE TABLE  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+write(*,*) 'PEMETAT0: ICE TABLE  DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !4. CO2 & H2O Adsorption
-
  if(adsorption_pem) then
   DO islope=1,nslope
@@ -318,13 +295,10 @@
        deltam_h2o_regolith_phys(:) = 0.  
     endif
-print *,'PEMETAT0: CO2 & H2O adsorption done  '
+write(*,*) 'PEMETAT0: CO2 & H2O adsorption DONE'
     endif
 endif ! soil_pem
 
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. 5 water reservoir
-
 #ifndef CPP_STD   
    call get_field("water_reservoir",water_reservoir,found)
@@ -342,12 +316,9 @@
 #endif
 
-
 call close_startphy
 
 else !No startfi, let's build all by hand
 
-    year_PEM=0
-
-    if(soil_pem) then
+    if (soil_pem) then
 
 !a) Thermal inertia
@@ -407,5 +378,4 @@
 enddo
 
-
 !!! zone de transition
 delta = depth_bedrock
@@ -424,9 +394,7 @@
  enddo
 
-print *,'PEMETAT0: TI  DONE'
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
+write(*,*) 'PEMETAT0: TI DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !b) Soil temperature 
     do islope = 1,nslope
@@ -458,31 +426,24 @@
         enddo
      enddo
- 
-      do isoil = nsoil_GCM+1,nsoil_PEM
-        do ig = 1,ngrid
-        watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
+
+        do isoil = nsoil_GCM+1,nsoil_PEM
+          do ig = 1,ngrid
+            watersoil_ave(ig,isoil,islope) = exp(beta_clap_h2o/tsoil_PEM(ig,isoil,islope) + alpha_clap_h2o)/tsoil_PEM(ig,isoil,islope)*mmol(igcm_h2o_vap)/(mugaz*r)
+          enddo
         enddo
-      enddo
 enddo !islope
-print *,'PEMETAT0: TSOIL DONE  '
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!c) Ice table    
-     call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
-     call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
-     do islope = 1,nslope
-     call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
-     enddo
-    
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-print *,'PEMETAT0: Ice table DONE  '
-
-
+write(*,*) 'PEMETAT0: TSOIL DONE'
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!c) Ice table
+       call computeice_table_equilibrium(ngrid,nslope,nsoil_PEM,watercaptag,watersurf_ave,watersoil_ave,TI_PEM(:,1,:),ice_table,ice_table_thickness)
+       call update_soil_thermalproperties(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,waterice,global_ave_pressure,ice_table,ice_table_thickness,TI_PEM)
+       do islope = 1,nslope
+         call soil_pem_ini(ngrid,nsoil_PEM,TI_PEM(:,:,islope),tsurf_ave_yr2(:,islope),tsoil_PEM(:,:,islope))
+       enddo
+       write(*,*) 'PEMETAT0: Ice table DONE'
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !d) Regolith adsorbed
-
  if(adsorption_pem) then
     m_co2_regolith_phys(:,:,:) = 0.
@@ -496,18 +457,10 @@
  endif
 
-print *,'PEMETAT0: CO2 adsorption done  '
-
+write(*,*) 'PEMETAT0: CO2 adsorption DONE'
 endif !soil_pem
 
-
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-    
-
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !. e) water reservoir
-
 #ifndef CPP_STD   
-
       do ig=1,ngrid
         if(watercaptag(ig)) then
@@ -517,5 +470,4 @@
         endif
       enddo
-
 #endif
 
@@ -527,7 +479,5 @@
       DO islope = 1,nslope
         DO iloop = 1,nsoil_PEM
-           if(isnan(tsoil_PEM(ig,iloop,islope))) then 
-         call abort_pem("PEM - pemetat0","NAN detected in Tsoil",1) 
-           endif
+          if (isnan(tsoil_PEM(ig,iloop,islope))) call abort_pem("PEM - pemetat0","NaN detected in Tsoil",1) 
         ENDDO
       ENDDO
@@ -535,5 +485,3 @@
 endif!soil_pem
 
-
-
-   END SUBROUTINE
+END SUBROUTINE
