SUBROUTINE update_soil(ngrid,nslope,nsoil_PEM,tendencies_waterice,waterice,p_avg_new,ice_depth,TI_PEM)
#ifndef CPP_STD
 USE comsoil_h, only:  inertiedat, volcapa
 USE comsoil_h_PEM, only: layer_PEM,n_1km,inertiedat_PEM,depth_breccia,depth_bedrock,index_breccia,index_bedrock,reg_thprop_dependp
 USE vertical_layers_mod, ONLY: ap,bp
 USE comsoil_h_PEM, only: n_1km
 USE constants_marspem_mod,only: TI_breccia,TI_bedrock, TI_regolith_avg
 implicit none
! Input: 
 INTEGER,INTENT(IN) ::  ngrid, nslope, nsoil_PEM
 REAL,INTENT(IN) :: p_avg_new
 REAL,INTENT(IN) :: tendencies_waterice(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 ::  inertie_thresold = 800. ! look for ice
 REAL ::  ice_inertia = 1200  ! Inertia of ice
 REAL ::  P610 = 610.0       ! current average pressure on Mars [Pa]

! Local variables:

 INTEGER :: ig,islope,iloop,iref,k
 REAL :: regolith_inertia(ngrid,nslope) ! TI of the regolith 
 REAL :: delta
 REAL :: TI_breccia_new
!  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 islope = 1,nslope
   regolith_inertia(:,islope) = inertiedat_PEM(:,1)
   do ig = 1,ngrid
      if((tendencies_waterice(ig,islope).lt.-1e-5).and.(waterice(ig,islope).eq.0)) then
              regolith_inertia(ig,islope) = TI_regolith_avg
       endif
      if(reg_thprop_dependp) then
          if(TI_PEM(ig,1,islope).lt.inertie_thresold) then
            regolith_inertia(ig,islope) = regolith_inertia(ig,islope)*(p_avg_new/P610)**0.3
          endif 
         TI_breccia_new = TI_breccia*(p_avg_new/P610)**0.3
      else 
         TI_breccia_new = TI_breccia
      endif
   enddo
 enddo


! 3. Build new Thermal Inertia 

do  islope=1,nslope
   do ig = 1,ngrid
     do iloop = 1,index_breccia
         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
     enddo
     if(regolith_inertia(ig,islope).lt.TI_breccia_new) then 
!!! transition
             delta = depth_breccia
             TI_PEM(ig,index_breccia+1,islope) = sqrt((layer_PEM(index_breccia+1)-layer_PEM(index_breccia))/ &
            (((delta-layer_PEM(index_breccia))/(TI_PEM(ig,index_breccia,islope)**2))+ &
                      ((layer_PEM(index_breccia+1)-delta)/(TI_breccia_new**2))))
            do iloop=index_breccia+2,index_bedrock
              TI_PEM(ig,iloop,islope) = TI_breccia_new
            enddo      
      else ! we keep the high ti values
            do iloop=index_breccia+1,index_bedrock
              TI_PEM(ig,iloop,islope) = TI_PEM(ig,index_breccia,islope)
            enddo 
       endif ! TI PEM and breccia comparison
!! transition
       delta = depth_bedrock
       TI_PEM(ig,index_bedrock+1,islope) = sqrt((layer_PEM(index_bedrock+1)-layer_PEM(index_bedrock))/ &
            (((delta-layer_PEM(index_bedrock))/(TI_PEM(ig,index_bedrock,islope)**2))+ &
            ((layer_PEM(index_bedrock+1)-delta)/(TI_bedrock**2))))
       do iloop=index_bedrock+2,nsoil_PEM
            TI_PEM(ig,iloop,islope) = TI_bedrock
       enddo   
      enddo ! ig
ENDDO ! islope

!  4. Build new TI in case of ice table
! a) For the regolith
  do ig=1,ngrid
    do islope=1,nslope
!      do iloop = 1,index_bedrock
!       TI_PEM(ig,iloop,islope) = TI_PEM(ig,1,islope)
!      enddo
  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)=ice_inertia
      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,1,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
#endif
      END