SUBROUTINE computeice_table(ngrid,nslope,nsoil_PEM,rhowatersurf_ave,rhowatersoil_ave,ice_table)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the ice table depth knowing the yearly average water
!!! density at the surface and at depth.
!!! Computations are made following the methods in Schorgofer et al., 2005
!!! 
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#ifndef CPP_STD
    USE comsoil_h_PEM, only: layer_PEM                             ! Depth of the vertical grid 
    implicit none 

    integer,intent(in) :: ngrid,nslope,nsoil_PEM                   ! Size of the physical grid, number of subslope, number of soil layer in the PEM 
    real,intent(in) :: rhowatersurf_ave(ngrid,nslope)              ! Water density at the surface, yearly averaged [kg/m^3]
    real,intent(in) :: rhowatersoil_ave(ngrid,nsoil_PEM,nslope)    ! Water density at depth, computed from clapeyron law's (Murchy and Koop 2005), yearly averaged  [kg/m^3]
    real,intent(inout) :: ice_table(ngrid,nslope)                  ! ice table depth [m]
    real :: z1,z2                                                  ! intermediate variables used when doing a linear interpolation between two depths to find the root 
    integer ig, islope,isoil                                       ! loop variables
    real :: diff_rho(nsoil_PEM)                                    ! difference of water vapor density between the surface and at depth [kg/m^3]


     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 (ig.eq.400) then 
         write(*,*) 'ig,islope,diff',ig,islope,diff_rho(:)
         write(*,*) 'rho surf',rhowatersurf_ave(ig,islope)
         write(*,*) 'rho soil',rhowatersoil_ave(ig,:,islope)
         endif
         if(diff_rho(1) > 0) then                                   ! ice is at the surface
           ice_table(ig,islope) = 0.
         else
           do isoil = 1,nsoil_PEM -1                                ! general case, we find the ice table depth by doing a linear approximation between the two depth, and then solve the first degree equation to find the root
             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