SUBROUTINE regolith_co2adsorption_old(ngrid,nslope,nsoil_PEM,ps,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,m_co2_completesoil,delta_mreg)
      use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km
      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
      use comslope_mod, only : subslope_dist,def_slope_mean
      IMPLICIT NONE
! Input: 
 INTEGER,INTENT(IN) :: ngrid, nslope, nsoil_PEM
 REAL,INTENT(IN) :: ps(ngrid)                                      ! Average surface pressure [Pa]
 REAL,INTENT(IN) :: tsoil_PEM(ngrid,nsoil_PEM,nslope)               ! Mean Soil Temperature [K]
 REAL,INTENT(IN) :: TI_PEM(ngrid,nsoil_PEM,nslope)                  ! Mean Thermal Inertia [USI]
 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope) !tendancies on the glaciers
 REAL,INTENT(IN) :: waterice(ngrid,nslope)                          ! water ice at the surface [KG/m^2]
 REAL,INTENT(IN) :: co2ice(ngrid,nslope)                            ! co2 ice at the surface [Kg/m^2] 

! Outputs:
 REAL,INTENT(INOUT) :: m_co2_completesoil(ngrid,nsoil_PEM,nslope)
 REAL,INTENT(INOUT) :: delta_mreg(ngrid)


  
! Constants:

 REAL ::  as = 1.9e4      ! Bulher & Piqueux 2020, assuming exchange only between 0 and 868m, last layer of the PEM
 REAL ::  delta = 5.749e1 ! Zent & Quinn, 1995
 REAL ::  beta =  -4.0711 ! Zent & Quinn, 1995
 REAL ::  gama =  0.2788  ! Zent & Quinn, 1995
 REAL ::  inertie_thresold = 2000. ! look for ice
 REAL ::  rho_regolith = 2000. ! density of the reoglith, buhler and piqueux 2021
! Input: Local 
 INTEGER :: ig,islope,iloop,iref,k
 REAL :: dm_co2_regolith_slope(ngrid,nsoil_PEM,nslope)   ! elementary mass adsorded per mesh per slope
 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)
 REAL :: deltam_reg_complete(ngrid,n_1km,nslope)
 REAL :: deltam_reg_slope(ngrid,nslope)

dm_co2_regolith_slope(:,:,:) = 0
delta_mreg(:) = 0.
!0. Initialization
  do ig = 1,ngrid
    do islope = 1,nslope
     if((abs(tend_h2oglaciers(ig,islope)).gt.1e-5).and.(abs(waterice(ig,islope)).gt.0)) then
        ispermanent_h2oglaciers(ig,islope) = 1
     else
        ispermanent_h2oglaciers(ig,islope) = 0
     endif

     if((abs(tend_co2glaciers(ig,islope)).gt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
        ispermanent_co2glaciers(ig,islope) = 1
     else
        ispermanent_co2glaciers(ig,islope) = 0
     endif
    enddo
   enddo
! 1. we compute the mass of co2 adsorded in each layer of the meshes  

 do ig = 1,ngrid
  do islope = 1,nslope
    do iloop = 1,n_1km
     if((TI_PEM(ig,iloop,islope).lt.inertie_thresold).and.(ispermanent_h2oglaciers(ig,islope).eq.0).and.(ispermanent_co2glaciers(ig,islope).eq.0)) then
     dm_co2_regolith_slope(ig,iloop,islope) = rho_regolith*as*delta*(ps(ig)**gama)*(tsoil_PEM(ig,iloop,islope))**beta 
     else
        if(abs(m_co2_completesoil(ig,iloop,islope)).lt.1-10) then !!! we are at first call
          dm_co2_regolith_slope(ig,iloop,islope) = rho_regolith*as*delta*(ps(ig)**gama)*(tsoil_PEM(ig,iloop,islope))**beta 
        else ! no change: permanent ice stick the atoms of CO2
          dm_co2_regolith_slope(ig,iloop,islope) = m_co2_completesoil(ig,iloop,islope)
        endif
     endif
  enddo
 enddo
enddo


! 2. Check the exchange between the atmosphere and the regolith

  do ig = 1,ngrid
   delta_mreg(ig) = 0.
   do islope = 1,nslope
    deltam_reg_slope(ig,islope) = 0.
    do iloop = 1,n_1km
       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)) 
       else ! NO EXCHANGE AS ICE BLOCK THE DYNAMIC!
             deltam_reg_complete(ig,iloop,islope) = 0.
       endif
       deltam_reg_slope(ig,islope) = deltam_reg_slope(ig,islope) + deltam_reg_complete(ig,iloop,islope)
    enddo
   delta_mreg(ig) = delta_mreg(ig) + deltam_reg_slope(ig,islope)*subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
   enddo
  enddo
   m_co2_completesoil(:,:,:) = dm_co2_regolith_slope(:,:,:)


!=======================================================================
      RETURN
      END