SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table)
#ifndef CPP_STD
    USE comsoil_h, only:  inertiedat, volcapa
    USE vertical_layers_mod, ONLY: ap,bp
    USE comsoil_h_PEM, only: layer_PEM 
    implicit none 

    integer,intent(in) :: ngrid,nslope,nsoil_PEM
    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)    ! surface temperature, timeseries [K]
    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! soil water density, yearly average [kg/m^3]
    real,intent(inout) :: ice_table(ngrid,nslope)                 ! ice table [m]
    real :: z1,z2
    integer ig, islope,isoil
    real :: diff_rho(nsoil_PEM)                    ! difference of vapor content


     do ig = 1,ngrid
      do islope = 1,nslope
           ice_table(ig,islope) = -10.

         do isoil = 1,nsoil_PEM
           diff_rho(isoil) = rhowatersurf_ave(ig,islope) - rhowatersoil_ave(ig,isoil,islope)

         enddo

         if(diff_rho(1) > 0) then
           ice_table(ig,islope) = 0.
         else
           do isoil = 1,nsoil_PEM -1
             if((diff_rho(isoil).lt.0).and.(diff_rho(isoil+1).gt.0.)) then
               z1 = (diff_rho(isoil) - diff_rho(isoil+1))/(layer_PEM(isoil) - layer_PEM(isoil+1))
               z2 = -layer_PEM(isoil+1)*z1 +  diff_rho(isoil+1)
               ice_table(ig,islope) = -z2/z1
               exit
             endif
           enddo
          endif
        enddo
      enddo

    

!=======================================================================
      RETURN

#endif
      END