SUBROUTINE pemetat0(startpem_file,filename,ngrid,nsoil_GCM,nsoil_PEM,nslope,timelen,timestep,TI_PEM,tsoil_PEM_yr1,tsoil_PEM,ice_table, &
                      tsurf_ave_yr1,tsurf_ave,q_co2,q_h2o,ps_inst,tsurf_inst,tsoil_inst,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice, &
                      global_ave_pressure, m_co2_regolith_phys,deltam_co2_regolith_phys, watersurf_ave,watersoil_ave)
   
   use iostart_PEM, only:  open_startphy, close_startphy, get_field, get_var
   use comsoil_h_PEM, only: layer_PEM, mlayer_PEM,n_1km,fluxgeo,inertiedat_PEM
   use comsoil_h, only:  volcapa,inertiedat 
   use soil_evolution_mod, only: soil_pem,soil_pem_CN
   use adsorption_mod, only : regolith_co2adsorption
   USE temps_mod_evol, ONLY: year_PEM

  implicit none
 
  character(len=*), intent(in) :: filename
  LOGICAL,intent(in) :: startpem_file
  character*8 :: fichnom
  integer,intent(in) :: ngrid
  integer,intent(in) :: nsoil_GCM
  integer,intent(in) :: nsoil_PEM
  integer,intent(in) :: nslope
  integer,intent(in) :: timelen
  real, intent(in) :: tsurf_ave_yr1(ngrid,nslope)                          ! surface temperature at the first year of GCM call
  real,intent(in) :: tsurf_ave(ngrid,nslope)                  ! surface temperature at the current year
  real,intent(in) :: q_co2(ngrid,timelen)                     ! MMR tracer co2 [kg/kg]
  real,intent(in) :: q_h2o(ngrid,timelen)                     ! MMR tracer h2o [kg/kg]
  real,intent(in) :: ps_inst(ngrid,timelen)                        ! surface pressure [Pa]
  real,intent(in) :: tsurf_inst(ngrid,nslope,timelen) ! soil (mid-layer) temperature
  real,intent(in) :: timestep	    ! time step
  real,intent(in) :: tend_h2oglaciers(ngrid,nslope)
  real,intent(in) :: tend_co2glaciers(ngrid,nslope)
  real,intent(in) :: co2ice(ngrid,nslope)
  real,intent(in) :: waterice(ngrid,nslope) 
  real, intent(in) :: tsoil_PEM_yr1(ngrid,nsoil_PEM,nslope) 
  real, intent(in) :: watersurf_ave(ngrid,nslope)     ! surface water ice density, yearly averaged  (kg/m^3)
  real, intent(inout) :: watersoil_ave(ngrid,nsoil_PEM,nslope)    ! surface water ice density, yearly averaged (kg/m^3)
  real, intent(in) :: global_ave_pressure                       ! mean average pressure on the planet
! outputs

  real,intent(inout) :: TI_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) thermal inertia (SI)
  real,intent(inout) :: tsoil_PEM(ngrid,nsoil_PEM,nslope) !soil (mid-layer) temperature (K)
  real,intent(inout) :: ice_table(ngrid,nslope) ! (m)
  real,intent(inout) :: tsoil_inst(ngrid,nsoil_PEM,nslope,timelen) ! instantaneous soil (mid-layer) temperature (k)
  real,intent(inout) :: m_co2_regolith_phys(ngrid,nsoil_PEM,nslope)   ! mass of co2 adsorbed (kg/m^2)
  real,intent(out) :: deltam_co2_regolith_phys(ngrid)                 ! mass of co2 that is exchanged due to adsorption desorption 
! local
   real :: tsoil_startPEM(ngrid,nsoil_PEM,nslope)
   real :: TI_startPEM(ngrid,nsoil_PEM,nslope)
   LOGICAL :: found
   integer :: iloop,ig,islope,it,isoil
   REAL :: TI_breccia = 750.
   REAL :: TI_bedrock = 2300.
   real :: kcond
   real :: delta
   CHARACTER*2 :: num
   real :: tsoil_saved(ngrid,nsoil_PEM)
   real :: tsoil_tmp_yr1(ngrid,nsoil_PEM,nslope)
   real :: tsoil_tmp_yr2(ngrid,nsoil_PEM,nslope)
   real :: tsoil_tmp(ngrid,nsoil_PEM,nslope)
   real :: alph_tmp(ngrid,nsoil_PEM-1)
   real :: beta_tmp(ngrid,nsoil_PEM-1)
   real :: co2_ads_prev(ngrid)
   real :: year_PEM_read
   real :: alpha_clap_h2o = -6143.7
   real :: beta_clap_h2o = 28.9074

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: read start_pem. Need a specific iostart_PEM
!!!
!!! Order: 1. Thermal Inertia 
!!!        2. Soil Temperature
!!!        3. Ice table
!!!        4. Mass of CO2 adsorbed
!!!
!!! /!\ This order must be respected !
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!1. Run

if (startpem_file) then
   ! open pem initial state file:
   call open_startphy(filename)

    if(soil_pem) then

   call get_var("Time",year_PEM_read,found)
   year_PEM=INT(year_PEM_read)
   if(.not.found) then
      write(*,*)'PEMetat0: failed loading year_PEM; take default=0'
   else
      write(*,*)'year_PEM of startpem=', year_PEM
   endif
 
!1. Thermal Inertia
! a. General case

DO islope=1,nslope
  write(num,fmt='(i2.2)') islope
  call get_field("TI_PEM_slope"//num,TI_startPEM(:,:,islope),found)
   if(.not.found) then
      write(*,*)'PEM settings: failed loading <TI_PEM_slope'//num//'>'
      write(*,*)'will reconstruct the values of TI_PEM'

      do ig = 1,ngrid
          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
!!! transition
             delta = 50.
             TI_PEM(ig,nsoil_GCM+1,islope) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))
            do iloop=nsoil_GCM+2,n_1km
              TI_PEM(ig,iloop,islope) = TI_breccia
            enddo      
          else ! we keep the high ti values
            do iloop=nsoil_GCM+1,n_1km
              TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
            enddo 
          endif ! TI PEM and breccia comparison
!! transition
          delta = 1000.
          TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
               (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
               ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
          do iloop=n_1km+2,nsoil_PEM
            TI_PEM(ig,iloop,islope) = TI_bedrock
          enddo   
      enddo
   else
     do iloop = nsoil_GCM+1,nsoil_PEM
       TI_PEM(:,iloop,islope) = TI_startPEM(:,iloop,islope)  ! ! 1st layers can change because of the presence of ice at the surface, so we don't change it here.
     enddo
   endif ! not found
ENDDO ! islope

print *,'PEMETAT0: THERMAL INERTIA DONE'

! b. Special case for inertiedat, inertiedat_PEM
call get_field("inertiedat_PEM",inertiedat_PEM,found)
if(.not.found) then
 do iloop = 1,nsoil_GCM
   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
 enddo 
!!! zone de transition
delta = 50.
do ig = 1,ngrid
if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))

 do iloop = nsoil_GCM+2,n_1km 
       inertiedat_PEM(ig,iloop) = TI_breccia
  enddo

else
   do iloop=nsoil_GCM+1,n_1km
      inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
   enddo 
endif ! comparison ti breccia
enddo!ig

!!! zone de transition
delta = 1000.
do ig = 1,ngrid
inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
enddo

 do iloop = n_1km+2, nsoil_PEM
   do ig = 1,ngrid
      inertiedat_PEM(ig,iloop) = TI_bedrock
   enddo
 enddo
endif ! not found

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!2. Soil Temperature

DO islope=1,nslope
  write(num,fmt='(i2.2)') islope
   call get_field("tsoil_PEM_slope"//num,tsoil_startPEM(:,:,islope),found)
  if(.not.found) then
      write(*,*)'PEM settings: failed loading <tsoil_PEM_slope'//num//'>'
      write(*,*)'will reconstruct the values of Tsoil'
      do ig = 1,ngrid
        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))   

       do iloop=nsoil_GCM+2,n_1km
            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
      enddo
        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))  

       do iloop=n_1km+2,nsoil_PEM
            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
      enddo
    enddo

   else
! predictor corrector: restart from year 1 of the GCM and build the evolution of
! tsoil at depth

    tsoil_tmp_yr1(:,:,islope) = tsoil_startPEM(:,:,islope)
    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
    call soil_pem_routine(ngrid,nsoil_PEM,.false.,TI_PEM(:,:,islope),timestep,tsurf_ave_yr1(:,islope),tsoil_tmp_yr1(:,:,islope),alph_tmp,beta_tmp)
    tsoil_tmp_yr2(:,:,islope) = tsoil_tmp_yr1(:,:,islope)     
    call soil_pem_routine(ngrid,nsoil_PEM,.true.,TI_PEM(:,:,islope),timestep,tsurf_ave(:,islope),tsoil_tmp_yr2(:,:,islope),alph_tmp,beta_tmp)


     do iloop = nsoil_GCM+1,nsoil_PEM
       tsoil_PEM(:,iloop,islope) = tsoil_tmp_yr2(:,iloop,islope)
     enddo
   endif

    do it = 1,timelen
        do isoil = nsoil_GCM+1,nsoil_PEM
        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
        enddo
     enddo
      do isoil = nsoil_GCM+1,nsoil_PEM
        do ig = 1,ngrid
        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
        enddo
      enddo

    
ENDDO

print *,'PEMETAT0: SOIL TEMP  DONE'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!3. Ice Table



     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)


print *,'PEMETAT0: ICE TABLE  DONE'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!4. CO2 Adsorption
DO islope=1,nslope
  write(num,fmt='(i2.2)') islope
   call get_field("mco2_reg_ads_slope"//num,m_co2_regolith_phys(:,:,islope),found)
   if(.not.found) then
      write(*,*)'PEM settings: failed loading <m_co2_regolith_phys>'
      write(*,*)'will reconstruct the values of co2 adsorbded'
   m_co2_regolith_phys(:,:,:) = 0.
   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)

   deltam_co2_regolith_phys(:) = 0.
   exit
   else
   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)

   endif
ENDDO


print *,'PEMETAT0: CO2 adsorption done  '

endif ! soil_pem

call close_startphy

else !No startfi, let's build all by hand

    year_PEM=0

    if(soil_pem) then

!a) Thermal inertia
   do islope = 1,nslope
      do ig = 1,ngrid

          if(TI_PEM(ig,nsoil_GCM,islope).lt.TI_breccia) then 
!!! transition
             delta = 50.
             TI_PEM(ig,nsoil_GCM+1,islope) =sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
            (((delta-layer_PEM(nsoil_GCM))/(TI_PEM(ig,nsoil_GCM,islope)**2))+ &
                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))

          do iloop=nsoil_GCM+2,n_1km
            TI_PEM(ig,iloop,islope) = TI_breccia
         enddo      
         else ! we keep the high ti values
           do iloop=nsoil_GCM+1,n_1km
                  TI_PEM(ig,iloop,islope) = TI_PEM(ig,nsoil_GCM,islope)
           enddo 
         endif

!! transition
             delta = 1000.
             TI_PEM(ig,n_1km+1,islope) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
            (((delta-layer_PEM(n_1km))/(TI_PEM(ig,n_1km,islope)**2))+ &
                      ((layer_PEM(n_1km+1)-delta)/(TI_breccia**2))))
          do iloop=n_1km+2,nsoil_PEM
            TI_PEM(ig,iloop,islope) = TI_bedrock
         enddo   
      enddo
enddo

 do iloop = 1,nsoil_GCM
   inertiedat_PEM(:,iloop) = inertiedat(:,iloop)
 enddo
!!! zone de transition
delta = 50.
do ig = 1,ngrid
      if(inertiedat_PEM(ig,nsoil_GCM).lt.TI_breccia) then 
inertiedat_PEM(ig,nsoil_GCM+1) = sqrt((layer_PEM(nsoil_GCM+1)-layer_PEM(nsoil_GCM))/ &
            (((delta-layer_PEM(nsoil_GCM))/(inertiedat(ig,nsoil_GCM)**2))+ &
                      ((layer_PEM(nsoil_GCM+1)-delta)/(TI_breccia**2))))


 do iloop = nsoil_GCM+2,n_1km 

       inertiedat_PEM(ig,iloop) = TI_breccia
   
 enddo
else 
   do iloop = nsoil_GCM+1,n_1km 
       inertiedat_PEM(ig,iloop) = inertiedat_PEM(ig,nsoil_GCM)
    enddo

endif
enddo


!!! zone de transition
delta = 1000.
do ig = 1,ngrid
inertiedat_PEM(ig,n_1km+1) = sqrt((layer_PEM(n_1km+1)-layer_PEM(n_1km))/ &
            (((delta-layer_PEM(n_1km))/(inertiedat_PEM(ig,n_1km)**2))+ &
                      ((layer_PEM(n_1km+1)-delta)/(TI_bedrock**2))))
enddo



 do iloop = n_1km+2, nsoil_PEM
   do ig = 1,ngrid
      inertiedat_PEM(ig,iloop) = TI_bedrock
   enddo
 enddo

print *,'PEMETAT0: TI  DONE'

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!b) Soil temperature 
    do islope = 1,nslope
     do ig = 1,ngrid
        kcond = (TI_PEM(ig,nsoil_GCM+1,islope)*TI_PEM(ig,nsoil_GCM+1,islope))/volcapa
        tsoil_PEM(ig,nsoil_GCM+1,islope) = tsoil_PEM(ig,nsoil_GCM,islope) + fluxgeo/kcond*(mlayer_PEM(nsoil_GCM)-mlayer_PEM(nsoil_GCM-1))

       do iloop=nsoil_GCM+2,n_1km
            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,nsoil_GCM+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(nsoil_GCM))
      enddo
        kcond = (TI_PEM(ig,n_1km+1,islope)*TI_PEM(ig,n_1km+1,islope))/volcapa
        tsoil_PEM(ig,n_1km+1,islope) = tsoil_PEM(ig,n_1km,islope) + fluxgeo/kcond*(mlayer_PEM(n_1km)-mlayer_PEM(n_1km-1))

       do iloop=n_1km+2,nsoil_PEM
            kcond = (TI_PEM(ig,iloop,islope)*TI_PEM(ig,iloop,islope))/volcapa
            tsoil_PEM(ig,iloop,islope) = tsoil_PEM(ig,n_1km+1,islope) + fluxgeo/kcond*(mlayer_PEM(iloop-1)-mlayer_PEM(n_1km))
      enddo
     enddo
    
    do it = 1,timelen
        do isoil = nsoil_GCM+1,nsoil_PEM
        tsoil_inst(:,isoil,islope,it) = tsoil_PEM(:,isoil,islope)
        enddo
     enddo
      do isoil = nsoil_GCM+1,nsoil_PEM
        do ig = 1,ngrid
        watersoil_ave(ig,isoil,islope) = exp(alpha_clap_h2o/tsoil_PEM(ig,isoil,islope) + beta_clap_h2o)/tsoil_PEM(ig,isoil,islope)
        enddo
      enddo
enddo
print *,'PEMETAT0: TSOIL DONE  '
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!c) Ice table    
     call computeice_table(ngrid,nslope,nsoil_PEM,watersurf_ave,watersoil_ave, ice_table)
     call update_soil(ngrid,nslope,nsoil_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,global_ave_pressure,ice_table,TI_PEM)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

print *,'PEMETAT0: Ice table DONE  '



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!d) Regolith adsorbed

   m_co2_regolith_phys(:,:,:) = 0.
   call regolith_co2adsorption(ngrid,nslope,nsoil_PEM,timelen,ps_inst,tsoil_PEM,TI_PEM,tend_h2oglaciers,tend_co2glaciers,co2ice,waterice,q_co2,q_h2o,m_co2_regolith_phys, deltam_co2_regolith_phys)
   deltam_co2_regolith_phys(:) = 0.

print *,'PEMETAT0: CO2 adsorption done  '

endif !soil_pem

endif ! of if (startphy_file)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
     if(soil_pem) then

    DO ig = 1,ngrid
      DO islope = 1,nslope

     write(*,*) 'ig,islope ,ice table=',ig,islope,ice_table(ig,islope)

      ENDDO
    ENDDO

!! small test
    DO ig = 1,ngrid
      DO islope = 1,nslope
        DO iloop = 1,nsoil_PEM
           if(isnan(tsoil_PEM(ig,iloop,islope))) then 
         write(*,*) "failed nan construction", ig, iloop, islope
           stop
           endif
        ENDDO
      ENDDO
     ENDDO

     endif!soil_pem

      write(*,*) "construction ok, no nan" 
     

   END SUBROUTINE