source: trunk/LMDZ.COMMON/libf/evolution/update_soil.F90 @ 2889

   SUBROUTINE update_soil(ngrid,nslope,nsoil_GCM,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
 USE vertical_layers_mod, ONLY: ap,bp
 USE comsoil_h_PEM, only: n_1km
 implicit none
! Input: 
 INTEGER,INTENT(IN) ::  ngrid, nslope,nsoil_GCM, 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 ::  inertie_averaged = 250 ! Mellon et al. 2000
 REAL ::  ice_inertia = 1200  ! Inertia of ice
 REAL ::  P610 = 610.0       ! current average pressure on Mars [Pa]
 REAL :: TI_breccia = 750.   ! THermal inertia of Breccia following Wood 2009 [SI]
 REAL :: TI_bedrock = 2300.  ! Thermal inertia of Bedrock following Wood 2009 [SI]

! 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) = inertie_averaged
       endif
        write(*,*) 'ig,islope',ig,islope,inertie_thresold,TI_PEM(ig,1,islope)
       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
   enddo
 enddo


! 3. Build new Thermal Inertia 

do  islope=1,nslope
   do ig = 1,ngrid
     do iloop = 1,nsoil_GCM
         TI_PEM(ig,iloop,islope) = regolith_inertia(ig,islope)
     enddo
     if(regolith_inertia(ig,islope).lt.TI_breccia_new) 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_new**2))))
            do iloop=nsoil_GCM+2,n_1km
              TI_PEM(ig,iloop,islope) = TI_breccia_new
            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 ! 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,n_1km
       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)=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,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