SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,p_avg_new,&
                          ice_depth,TI_PEM)

 USE comsoil_h, only:  inertiedat, volcapa
 USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM
 USE vertical_layers_mod, ONLY: ap,bp
 implicit none
! Input: 
 INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
 REAL,INTENT(IN) :: tend_h2oglaciers(ngrid,nslope),tend_co2glaciers(ngrid,nslope)
 REAL,INTENT(IN) :: p_avg_new
 REAL,INTENT(IN) :: co2ice(ngrid,nslope)
 REAL,INTENT(IN) :: waterice(ngrid,nslope)
 REAL,INTENT(in) :: ice_depth(ngrid,nslope)

! Outputs:

 REAL,INTENT(INOUT) :: TI_PEM(ngrid,nsoil_PEM,nslope)
  
! Constants:

 REAL ::  alpha = 0.2
 REAL ::  beta = 1.08e7
 REAL ::  To = 273.15
 REAL ::  R =  8.314
 REAL ::  L =  51058.
 REAL ::  inertie_thresold = 800. ! look for ice
 REAL ::  inertie_co2glaciers = 2120 ! Mellon et al. 2000
 REAL ::  inertie_averaged = 250 ! Mellon et al. 2000
 REAL ::  ice_inertia = 2120  ! Inertia of ice
 REAL ::  surfaceice_inertia = 800  ! Inertia of ice
 REAL ::  P610 = 610.0
 REAL ::  m_h2o = 18.01528E-3
 REAL ::  m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
 REAL ::  m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
 REAL ::  A,B,mmean             

! Local variables:

 INTEGER :: ig,islope,iloop,iref,k
 REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith 
 REAL :: d(ngrid,nsoil_PEM,nslope)
 REAL :: Total_surface
 INTEGER :: ispermanent_co2glaciers(ngrid,nslope)
 INTEGER :: ispermanent_h2oglaciers(ngrid,nslope)

! 0. Initialisation

  do ig = 1,ngrid
    do islope = 1,nslope
     if((abs(tend_h2oglaciers(ig,islope)).lt.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)).lt.1e-5).and.(abs(co2ice(ig,islope)).gt.0)) then
        ispermanent_co2glaciers(ig,islope) = 1
     else
        ispermanent_co2glaciers(ig,islope) = 0
     endif
    enddo

  enddo

 ispermanent_co2glaciers(:,:) = 0
 ispermanent_h2oglaciers(:,:) = 0

!  1.Ice TI feedback

    do islope = 1,nslope
     call soil_TIfeedback_PEM(ngrid,nsoil_PEM,waterice(:,islope),   TI_PEM(:,:,islope))
    enddo

! 2. Modification of the regolith thermal inertia.

  do ig=1,ngrid
    do islope=1,nslope
     do iloop =1,n_1km
       d(ig,iloop,islope) = ((inertiedat_PEM(ig,iloop)*inertiedat_PEM(ig,iloop))/(volcapa*alpha*P610**0.6))**(-1/(0.11*log10(P610/beta)))
   if(TI_PEM(ig,iloop,islope).lt.inertie_thresold) then !           we are modifying the regolith properties, not ice
          TI_PEM(ig,iloop,islope) = sqrt(volcapa*alpha*(p_avg_new**0.6)* d(ig,iloop,islope)**(-0.11*log10(p_avg_new/beta)))

   endif 
     enddo
   enddo
enddo

!  3. Build new TI for the PEM
! a) For the regolith
  do ig=1,ngrid
    do islope=1,nslope

  if (ice_depth(ig,islope).gt. -1.e-10) then
! 3.0 FIrst if permanent ice
   if (ice_depth(ig,islope).lt. 1e-10) then
      do iloop = 1,nsoil_PEM
       TI_PEM(ig,iloop,islope)=max(ice_inertia,inertiedat_PEM(ig,iloop))
      enddo
     else 
      ! 4.1 find the index of the mixed layer
      iref=0 ! initialize iref
      do k=1,nsoil_PEM ! loop on layers
        if (ice_depth(ig,islope).ge.layer_PEM(k)) then
          iref=k ! pure regolith layer up to here
        else
         ! correct iref was obtained in previous cycle
         exit
        endif
      enddo

      ! 4.2 Build the new ti
      do iloop=1,iref
         TI_PEM(ig,iloop,islope) =TI_PEM(ig,iloop,islope)
      enddo
      if (iref.lt.nsoil_PEM) then
        if (iref.ne.0) then
          ! mixed layer
           TI_PEM(ig,iref+1,islope)=sqrt((layer_PEM(iref+1)-layer_PEM(iref))/ &
            (((ice_depth(ig,islope)-layer_PEM(iref))/(TI_PEM(ig,iref,islope)**2))+ &
                      ((layer_PEM(iref+1)-ice_depth(ig,islope))/(ice_inertia**2))))
        else ! first layer is already a mixed layer
          ! (ie: take layer(iref=0)=0)
          TI_PEM(ig,1,islope)=sqrt((layer_PEM(1))/ &
                          (((ice_depth(ig,islope))/(TI_PEM(ig,1,islope)**2))+ &
                           ((layer_PEM(1)-ice_depth(ig,islope))/(ice_inertia**2))))
        endif ! of if (iref.ne.0)        
        ! lower layers of pure ice
        do iloop=iref+2,nsoil_PEM
          TI_PEM(ig,iloop,islope)=ice_inertia
        enddo
      endif ! of if (iref.lt.(nsoilmx))
      endif ! permanent glaciers
      endif ! depth > 0
     enddo !islope
    enddo !ig

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