SUBROUTINE conduction(ngrid,nlayer,ptimestep,    &
                          pplay,pt,pzlay,pzlev,   &
                          zdtconduc,tsurf)

      use comcstfi_mod, only: cpp, r
      use surfdat_h, only: phitop

      IMPLICIT NONE

!=======================================================================
!
!   Molecular thermal conduction
!
!   N. Descamp, F. Forget 05/1999
!
!=======================================================================

!-----------------------------------------------------------------------
!   declarations:
!-----------------------------------------------------------------------

#include "dimensions.h"

!   arguments:
!   ----------

      INTEGER ngrid,nlayer
      REAL ptimestep
      REAL pt(ngrid,nlayer),zdtconduc(ngrid,nlayer)
      REAL pzlay(ngrid,nlayer),pzlev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL zt(ngrid,nlayer)

!   local:
!   ------

      INTEGER ilayer,ig
      REAL alpha(nlayer)
      REAL lambda(nlayer),muvol(nlayer)
      REAL C(nlayer),D(nlayer),den(nlayer)
      REAL tsurf
!      REAL pdt(nlayer),tsurf
      REAL pdt(nlayer)

!   constants used locally
!    ---------------------
!     The atmospheric conductivity is a function of temperature T :
!      conductivity = Akk* T**skk
      real,parameter :: Akk=5.63E-05
      real,parameter :: skk=1.12

      logical, save :: firstcall=.true.
!-----------------------------------------------------------------------
!   calcul des coefficients alpha et lambda
!
!-----------------------------------------------------------------------

      IF (firstcall) THEN
        write (*,*) 'coeff to compute molecular conductivity Akk,skk'
        write(*,*) Akk,skk
        firstcall = .false.
      END IF

!       DO ilayer=1,nlayer
!           do ig=1,ngrid

!           zt(ig,ilayer)=pt(ig,ilayer)+pdt(ig,ilayer)*ptimestep

!           enddo
!       ENDDO


      DO ig=1,ngrid

        lambda(1) = Akk*tsurf**skk / pzlay(ig,1)
        DO ilayer = 2 , nlayer
          lambda(ilayer) = Akk * pt(ig,ilayer)**skk /   &
                                (pzlay(ig,ilayer)-pzlay(ig,ilayer-1))
        ENDDO

        DO ilayer=1,nlayer-1
          muvol(ilayer)=pplay(ig,ilayer)/(r*pt(ig,ilayer))
          alpha(ilayer)=cpp*(muvol(ilayer)/ptimestep)*  &
                                (pzlev(ig,ilayer+1)-pzlev(ig,ilayer))
        ENDDO
!!!!!! alerte !!!!!c
!!!!!! zlev n'est pas declare a nlev !!!!!
!!!!!! ---->
          muvol(nlayer)=pplay(ig,nlayer)/(r*pt(ig,nlayer))
          alpha(nlayer)=cpp*(muvol(nlayer)/ptimestep)*  &
                                2*(pzlay(ig,nlayer)-pzlev(ig,nlayer))
!!!!!! <----
!!!!!! alerte !!!!!c
!        write(*,*) lambda(1),muvol(1),tsurf,pt(1,1)


!--------------------------------------------------------------------
!
!     calcul des coefficients C et D
!
!-------------------------------------------------------------------

      den(1)=alpha(1)+lambda(2)+lambda(1)
      C(1)=lambda(1)*(tsurf-pt(ig,1))+lambda(2)*(pt(ig,2)-pt(ig,1))
      C(1)=C(1)/den(1)
      D(1)=lambda(2)/den(1)

      DO ilayer = 2,nlayer-1
         den(ilayer)=alpha(ilayer)+lambda(ilayer+1)
         den(ilayer)=den(ilayer)+lambda(ilayer)*(1-D(ilayer-1))

         C(ilayer) =lambda(ilayer+1)*(pt(ig,ilayer+1)-pt(ig,ilayer))    &
            +lambda(ilayer)*(pt(ig,ilayer-1)-pt(ig,ilayer)+C(ilayer-1))
         C(ilayer) =C(ilayer)/den(ilayer)

         D(ilayer) =lambda(ilayer+1) / den(ilayer)
      ENDDO

      den(nlayer)=alpha(nlayer) + lambda(nlayer) * (1-D(nlayer-1))
!      C(nlayer)=((C(nlayer-1)+pt(ig,nlayer-1)-pt(ig,nlayer))
!     $            *lambda(nlayer) + phitop) / den(nlayer)
      C(nlayer)=C(nlayer-1)+pt(ig,nlayer-1)-pt(ig,nlayer)
      C(nlayer)=(C(nlayer)*lambda(nlayer)+phitop) / den(nlayer)
    !   C(nlayer)=(C(nlayer)*lambda(nlayer)+phitop(ig)) / den(nlayer)

!----------------------------------------------------------------------
!
!      calcul de la nouvelle temperature ptconduc
!
!----------------------------------------------------------------------

      DO ilayer=1,nlayer
         pdt(ilayer)=0.
      ENDDO

         pdt(nlayer)=C(nlayer)
!         pt(nlayer)=ttop
!         write(*,*)'pt',pt(nlayer)

      DO ilayer=nlayer-1,1,-1
         pdt(ilayer)=C(ilayer)+D(ilayer)*pdt(ilayer+1)
      ENDDO

!      write(*,*) alpha(1),lambda(2),den(1),C(1),D(1),pdt(1)

!       write(*,*)'15: c et d',C(nlayer),D(nlayer),ptconduc(nlayer),
!     $    pt(ig,nlayer),ig
!       write(*,*)'14: c et d',C(14),D(14),ptconduc(14),
!     $    pt(ig,14),ig
!       write(*,*)'3: c et d',C(3),D(3),ptconduc(3),
!     $    pt(ig,3)
!       write(*,*)'2: c et d',C(2),D(2),ptconduc(2),
!     $    pt(ig,2)

!-----------------------------------------------------------------------
!
!     calcul de la tendance zdtconduc
!
!-----------------------------------------------------------------------

      DO ilayer=1,nlayer
         zdtconduc(ig,ilayer) = pdt(ilayer) / ptimestep
!         write(*,*)'zdtconduc,c,d',ilayer,zdtconduc(ig,ilayer),C(ilayer)
!     $,D(ilayer)
!         write(*,*) 'pt',ilayer,pt(ig,ilayer)
      ENDDO

      ENDDO             ! boucle sur ngrid
!       write(*,*)'zdtconduc',zdtconduc(ngrid,nlayer),ngrid
!        write(*,*)'*********************************'
      RETURN
      END