Index: /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3263)
+++ /trunk/LMDZ.COMMON/libf/evolution/adsorption_mod.F90	(revision 3264)
@@ -218,6 +218,11 @@
     do iloop = 1,index_breccia
        if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-             deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
-                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop))
+             if(iloop==1) then
+                 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
+                                                   *(layer_PEM(iloop))
+             else
+                 deltam_reg_complete(ig,iloop,islope) = (dm_h2o_regolith_slope(ig,iloop,islope) - m_h2o_completesoil(ig,iloop,islope)) &
+                                                   *(layer_PEM(iloop) - layer_PEM(iloop-1))
+             endif
        else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
              deltam_reg_complete(ig,iloop,islope) = 0.
@@ -387,6 +392,11 @@
     do iloop = 1,index_breccia
        if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
-             deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
-                                                   *(layer_PEM(iloop+1) - layer_PEM(iloop))
+             if(iloop == 1) then 
+                 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
+                                                   *(layer_PEM(iloop))               
+             else
+                 deltam_reg_complete(ig,iloop,islope) = (dm_co2_regolith_slope(ig,iloop,islope) - m_co2_completesoil(ig,iloop,islope)) &
+                                                   *(layer_PEM(iloop) - layer_PEM(iloop-1))
+             endif
        else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
              deltam_reg_complete(ig,iloop,islope) = 0.
Index: /trunk/LMDZ.COMMON/libf/evolution/changelog.txt
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3263)
+++ /trunk/LMDZ.COMMON/libf/evolution/changelog.txt	(revision 3264)
@@ -247,2 +247,5 @@
 == 07/03/2024 == JBC
 Adaptation of threshold values for ice management (in particular 'inf_h2oice_threshold') to more realistic values.
+
+== 13/03/2024 == LL
+Correction of small mistake when computing the total mass of CO2/H2O stored in the soil; the wrong layer index was used
Index: /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3263)
+++ /trunk/LMDZ.COMMON/libf/evolution/ice_table_mod.F90	(revision 3264)
@@ -137,5 +137,5 @@
 
 
-   SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsoil,delta_m_h2o)
+   SUBROUTINE compute_massh2o_exchange_ssi(ngrid,nslope,nsoil_PEM,former_ice_table_thickness,new_ice_table_thickness,ice_table_depth,tsurf,tsoil,delta_m_h2o)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!
@@ -147,5 +147,6 @@
     use comsoil_h_PEM, only: mlayer_PEM
     use comslope_mod, only: subslope_dist,def_slope_mean
-    use constants_marspem_mod,only:porosity
+    use constants_marspem_mod, only: porosity
+    use vdifc_mod, only: compute_Tice
 #ifndef CPP_STD
     use comcstfi_h,   only: pi
@@ -160,5 +161,6 @@
     real,intent(in) :: former_ice_table_thickness(ngrid,nslope)    ! ice table thickness at the former iteration [m]
     real,intent(in) :: new_ice_table_thickness(ngrid,nslope)       ! ice table thickness at the current iteration [m]
-    real,intent(out) :: ice_table_depth(ngrid,nslope)              ! ice table depth [m]
+    real,intent(in) :: ice_table_depth(ngrid,nslope)              ! ice table depth [m]
+    real,intent(in) :: tsurf(ngrid,nslope)                         ! Surface temperature [K]
     real,intent(in) :: tsoil(ngrid,nsoil_PEM,nslope)               ! Soil temperature [K]
 !   outputs
@@ -177,26 +179,6 @@
    do ig = 1,ngrid
       do islope = 1,nslope
-         if(ice_table_depth(ig,islope).gt.0.) then
-            if (ice_table_depth(ig,islope).lt. 1e-10) then
-               Tice(ig,islope) = tsoil(ig,1,islope)
-            else
-               iref=0 ! initialize iref
-               do ilay=1,nsoil_PEM ! loop on layers to find the beginning of the ice table
-                  if (ice_table_depth(ig,islope).ge.mlayer_PEM(ilay)) then
-                     iref=ilay ! pure regolith layer up to here
-                  else
-                  ! correct iref was obtained in previous cycle
-                     exit
-                  endif
-               enddo
-               if(iref.eq.nsoil_PEM) then
-                  Tice(ig,islope) = tsoil(ig,nsoil_PEM,islope)
-               else
-                  Tice(ig,islope) = (tsoil(ig,iref+1,islope) - tsoil(ig,iref,islope))/(mlayer_PEM(iref+1) - mlayer_PEM(iref))* &
-                                    (ice_table_depth(ig,islope) - mlayer_PEM(iref)) +  tsoil(ig,iref,islope)           
-               endif         
-               rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
-            endif
-         endif
+         call compute_Tice(nsoil_PEM,tsoil(ig,:,islope),tsurf(ig,islope), ice_table_depth(ig,islope), Tice(ig,islope))
+         rho(ig,islope) = -3.5353e-4*Tice(ig,islope)**2+   0.0351* Tice(ig,islope) +  933.5030 ! Rottgers, 2012
       enddo
    enddo      
Index: /trunk/LMDZ.COMMON/libf/evolution/pem.F90
===================================================================
--- /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3263)
+++ /trunk/LMDZ.COMMON/libf/evolution/pem.F90	(revision 3264)
@@ -599,8 +599,15 @@
         do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
-                totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
+                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)
-                totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
+                   totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
                                        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))* &
+                                       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))* &
+                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                endif
             enddo
         enddo
@@ -840,5 +847,5 @@
             porefillingice_thickness_prev_iter = porefillingice_thickness
             call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
-            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
+            call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsurf_ave, tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
         endif
 ! II_d.4 Update the soil thermal properties
@@ -855,9 +862,16 @@
             do ig = 1,ngrid
                 do islope =1, nslope
-                    do l = 1,nsoilmx_PEM - 1
-                        totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l + 1) - layer_PEM(l))* &
-                                               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 + 1) - layer_PEM(l))* &
-                                               subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                    do l = 1,nsoilmx_PEM 
+                        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)
+                            totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &
+                                       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))* &
+                                       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))* &
+                                       subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
+                        endif
                     enddo
                 enddo
