Index: trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3492)
+++ trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3493)
@@ -54,6 +54,6 @@
 use orbit_param_criterion_mod,  only: orbit_param_criterion
 use recomp_orb_param_mod,       only: recomp_orb_param
-use ice_table_mod,              only: porefillingice_depth, porefillingice_thickness, end_ice_table_porefilling, &
-                                      ini_ice_table_porefilling, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium,compute_massh2o_exchange_ssi
+use ice_table_mod,              only: icetable_depth, icetable_thickness, end_ice_table, ice_porefilling, &
+                                      ini_ice_table, icetable_equilibrium, icetable_dynamic, computeice_table_equilibrium, compute_massh2o_exchange_ssi
 use soil_thermalproperties_mod, only: update_soil_thermalproperties
 use time_phylmdz_mod,           only: daysec, dtphys
@@ -204,30 +204,28 @@
 
 ! Variables for surface and soil
-real, dimension(:,:),     allocatable :: tsurf_avg                          ! Physic x SLOPE field: Averaged Surface Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsoil_avg                          ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
-real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries               ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries          ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
-real, dimension(:,:,:,:), allocatable :: tsoil_anom                         ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
-real, dimension(:,:),     allocatable :: tsurf_avg_yr1                      ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
-real, dimension(:,:),     allocatable :: Tsoil_locslope                     ! Physic x Soil: Intermediate when computing Tsoil [K]
-real, dimension(:),       allocatable :: Tsurf_locslope                     ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries       ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
-real, dimension(:,:),     allocatable :: watersurf_density_avg              ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
-real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries   ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
-real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg          ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
-real, dimension(:,:),     allocatable :: Tsurfavg_before_saved              ! Surface temperature saved from previous time step [K]
-real, dimension(:),       allocatable :: delta_co2_adsorbded                ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
-real, dimension(:),       allocatable :: delta_h2o_adsorbded                ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
-real                                  :: totmassco2_adsorbded               ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
-real                                  :: totmassh2o_adsorbded               ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
-logical                               :: bool_sublim                        ! logical to check if there is sublimation or not
-logical, dimension(:,:),  allocatable :: co2ice_disappeared                 ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
-real, dimension(:,:),     allocatable :: porefillingice_thickness_prev_iter ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
-real, dimension(:),       allocatable :: delta_h2o_icetablesublim           ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
-
-! Dynamic ice table
-real, dimension(:,:),     allocatable :: ss_ice_pf                          ! the amout of porefilling in each layer in each grid [m^3/m^3]
-real, dimension(:,:),     allocatable :: ss_ice_pf_new                      ! the amout of porefilling in each layer in each grid after the ss module[m^3/m^3]
-real, dimension(:,:),     allocatable :: porefillingice_depth_new           ! new pure ice table depth
+real, dimension(:,:),     allocatable :: tsurf_avg                        ! Physic x SLOPE field: Averaged Surface Temperature [K]
+real, dimension(:,:,:),   allocatable :: tsoil_avg                        ! Physic x SOIL x SLOPE field: Averaged Soil Temperature [K]
+real, dimension(:,:,:),   allocatable :: tsurf_PCM_timeseries             ! ngrid x SLOPE x TIMES field: Surface Temperature in timeseries [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_phys_PEM_timeseries        ! IG x SOIL x SLOPE x TIMES field: Non averaged Soil Temperature [K]
+real, dimension(:,:,:,:), allocatable :: tsoil_anom                       ! IG x SOIL x SLOPE x TIMES field: Amplitude between instataneous and yearly average soil temperature at the beginning [K]
+real, dimension(:,:),     allocatable :: tsurf_avg_yr1                    ! Physic x SLOPE field: Averaged Surface Temperature of first call of the PCM [K]
+real, dimension(:,:),     allocatable :: Tsoil_locslope                   ! Physic x Soil: Intermediate when computing Tsoil [K]
+real, dimension(:),       allocatable :: Tsurf_locslope                   ! Physic x Soil: Intermediate surface temperature to compute Tsoil [K]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_timeseries     ! Physic x Soil x Slope x Times water soil density, time series [kg /m^3]
+real, dimension(:,:),     allocatable :: watersurf_density_avg            ! Physic x Slope, water surface density, yearly averaged [kg/m^3]
+real, dimension(:,:,:,:), allocatable :: watersoil_density_PEM_timeseries ! Physic x Soil x Slope x Times, water soil density, time series [kg/m^3]
+real, dimension(:,:,:),   allocatable :: watersoil_density_PEM_avg        ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3]
+real, dimension(:,:),     allocatable :: Tsurfavg_before_saved            ! Surface temperature saved from previous time step [K]
+real, dimension(:),       allocatable :: delta_co2_adsorbded              ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]
+real, dimension(:),       allocatable :: delta_h2o_adsorbded              ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]
+real                                  :: totmassco2_adsorbded             ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]
+real                                  :: totmassh2o_adsorbded             ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]
+logical                               :: bool_sublim                      ! logical to check if there is sublimation or not
+logical, dimension(:,:),  allocatable :: co2ice_disappeared               ! logical to check if a co2 ice reservoir already disappeared at a previous timestep
+real, dimension(:,:),     allocatable :: icetable_thickness_prev_iter     ! ngrid x nslope: Thickness of the ice table at the previous iteration [m]
+real, dimension(:),       allocatable :: delta_h2o_icetablesublim         ! ngrid x Total mass of the H2O that has sublimated / condenses from the ice table [kg]
+real, dimension(:),       allocatable :: porefill                         ! Pore filling (output) to compute the dynamic ice table
+real                                  :: ssi_depth                        ! Ice table depth (output) to compute the dynamic ice table
+
 ! Some variables for the PEM run
 real, parameter :: year_step = 1   ! Timestep for the pem
@@ -557,5 +555,5 @@
 allocate(delta_co2_adsorbded(ngrid))
 allocate(co2ice_disappeared(ngrid,nslope))
-allocate(porefillingice_thickness_prev_iter(ngrid,nslope))
+allocate(icetable_thickness_prev_iter(ngrid,nslope))
 allocate(delta_h2o_icetablesublim(ngrid))
 allocate(delta_h2o_adsorbded(ngrid))
@@ -583,6 +581,6 @@
 call end_adsorption_h_PEM
 call ini_adsorption_h_PEM(ngrid,nslope,nsoilmx_PEM)
-call end_ice_table_porefilling
-call ini_ice_table_porefilling(ngrid,nslope)
+call end_ice_table
+call ini_ice_table(ngrid,nslope,nsoilmx_PEM)
 
 if (soil_pem) then
@@ -646,9 +644,9 @@
 endif
 
-call pemetat0("startpem.nc",ngrid,nsoilmx,nsoilmx_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM,porefillingice_depth, &
-              porefillingice_thickness,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys,ps_timeseries,          &
-              tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,global_avg_press_PCM,              &
-              watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,delta_co2_adsorbded,                &
-              h2o_adsorbded_phys,delta_h2o_adsorbded,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,tsoil_phys_PEM_timeseries,tend_h2o_ice,tend_co2_ice,co2_ice,h2o_ice,               &
+              global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbded_phys,          &
+              delta_co2_adsorbded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)
 
 ! We save the places where h2o ice is sublimating
@@ -928,18 +926,24 @@
 ! II_d.3 Update the ice table
         if (icetable_equilibrium) then
-            write(*,*) "Compute ice table"
-            porefillingice_thickness_prev_iter = porefillingice_thickness
-            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            write(*,*) "Compute ice table (equilibrium method)"
+            icetable_thickness_prev_iter = icetable_thickness
+            call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_avg,watersoil_density_PEM_avg,TI_PEM(:,1,:),icetable_depth,icetable_thickness)
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         else if (icetable_dynamic) then
-            write(*,*) "Compute ice table"
-            porefillingice_thickness_prev_iter = porefillingice_thickness
-            call dyn_ss_ice_m_wrapper(ngrid,nsoilmx,TI_PEM,ps,mmol(igcm_h2o_vap),tsoil_PEM,porefillingice_depth,ss_ice_pf,ss_ice_pf_new,porefillingice_depth_new)
-            
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            write(*,*) "Compute ice table (dynamic method)"
+            allocate(porefill(nsoilmx_PEM))
+            do ig = 1,ngrid
+                do islope = 1,nslope
+                    call dyn_ss_ice_m(icetable_depth(ig,islope),tsurf_avg(ig,islope),tsoil_PEM(ig,:,islope),nsoilmx_PEM,TI_PEM(ig,1,nslope),ps(ig),sum(q_h2o_PEM_phys(ig,:))/size(q_h2o_PEM_phys,2),ice_porefilling(ig,:,islope),porefill,ssi_depth)
+                    icetable_depth(ig,islope) = ssi_depth
+                    ice_porefilling(ig,:,islope) = porefill
+                enddo
+            enddo
+            deallocate(porefill)
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,icetable_thickness_prev_iter,icetable_thickness,icetable_depth,tsurf_avg, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
 
 ! II_d.4 Update the soil thermal properties
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tend_h2o_ice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,TI_PEM)
 
 ! II_d.5 Update the mass of the regolith adsorbed
@@ -952,7 +956,7 @@
             totmassh2o_adsorbded = 0.
             do ig = 1,ngrid
-                do islope =1, nslope
+                do islope = 1,nslope
                     do l = 1,nsoilmx_PEM 
-                        if (l==1) then
+                        if (l == 1) then
                             totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
@@ -960,7 +964,7 @@
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         else
-                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
+                            totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
-                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &
+                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &
                                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         endif
@@ -987,9 +991,9 @@
         call writediagpem(ngrid,'tsurf_slope'//str2,'tsurf','K',2,tsurf(:,islope))
         if (icetable_equilibrium) then
-            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
-            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
+            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
+            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
         else if (icetable_dynamic) then
-            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,porefillingice_depth(:,islope))
-            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,porefillingice_thickness(:,islope))
+            call writediagpem(ngrid,'ssi_depth_slope'//str2,'ice table depth','m',2,icetable_depth(:,islope))
+            call writediagpem(ngrid,'ssi_thick_slope'//str2,'ice table depth','m',2,icetable_thickness(:,islope))
         endif
 
@@ -1205,5 +1209,5 @@
 call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist)
 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM, &
-             TI_PEM,porefillingice_depth,porefillingice_thickness,       &
+             TI_PEM,icetable_depth,icetable_thickness,ice_porefilling,   &
              co2_adsorbded_phys,h2o_adsorbded_phys,h2o_ice,stratif)
 write(*,*) "restartpem.nc has been written"
@@ -1227,5 +1231,5 @@
 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved)
 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg)
-deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,porefillingice_thickness_prev_iter)
+deallocate(delta_co2_adsorbded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim,icetable_thickness_prev_iter)
 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif)
 !----------------------------- END OUTPUT ------------------------------
