SUBROUTINE conduction(ngrid,nlayer,ptimestep,
     $                   pplay,pt,pzlay,pzlev,
     $                   zdtconduc,tsurf)
      IMPLICIT NONE

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

c-----------------------------------------------------------------------
c   declarations:
c-----------------------------------------------------------------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "surfdat.h"

c   arguments:
c   ----------

      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)

c   local:
c   ------

      INTEGER ilayer,ig
      REAL alpha(nlayermx)
      REAL lambda(nlayermx),muvol(nlayermx)
      REAL C(nlayermx),D(nlayermx),den(nlayermx)
      REAL tsurf
c      REAL pdt(nlayermx),tsurf
      REAL pdt(nlayermx)

c   constants used locally
c    ---------------------
c     The atmospheric conductivity is a function of temperature T :
c      conductivity = Akk* T**skk
      real,parameter :: Akk=5.63E-05
      real,parameter :: skk=1.12
      
      logical firstcall
      save firstcall
      data firstcall /.true./
c-----------------------------------------------------------------------
c   calcul des coefficients alpha et lambda
c
c-----------------------------------------------------------------------

      IF (firstcall) THEN
        write (*,*) 'coeff to compute molecular conductivity Akk,skk'
        write(*,*) Akk,skk
        firstcall = .false.
      END IF
     
c       DO ilayer=1,nlayer
c           do ig=1,ngrid

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

c           enddo
c       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
c!!!!! alerte !!!!!c
c!!!!! zlev n'est pas declare a nlev !!!!!
c!!!!! ---->
          muvol(nlayer)=pplay(ig,nlayer)/(r*pt(ig,nlayer)) 
          alpha(nlayer)=cpp*(muvol(nlayer)/ptimestep)*
     $                      2*(pzlay(ig,nlayer)-pzlev(ig,nlayer))      
c!!!!! <----
c!!!!! alerte !!!!!c
c        write(*,*) lambda(1),muvol(1),tsurf,pt(1,1)   


c--------------------------------------------------------------------
c
c     calcul des coefficients C et D
c
c-------------------------------------------------------------------

      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      C(nlayer)=((C(nlayer-1)+pt(ig,nlayer-1)-pt(ig,nlayer))
c     $            *lambda(nlayer) + phitop) / den(nlayer) 
      C(nlayer)=C(nlayer-1)+pt(ig,nlayer-1)-pt(ig,nlayer) 
      C(nlayer)=(C(nlayer)*lambda(nlayer)+phitop(ig)) / den(nlayer) 
       	 	
c----------------------------------------------------------------------
c
c      calcul de la nouvelle temperature ptconduc
c
c----------------------------------------------------------------------

      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 
       
c      write(*,*) alpha(1),lambda(2),den(1),C(1),D(1),pdt(1)   

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

c-----------------------------------------------------------------------
c
c     calcul de la tendance zdtconduc
c
c-----------------------------------------------------------------------
    
      DO ilayer=1,nlayer
         zdtconduc(ig,ilayer) = pdt(ilayer) / ptimestep
c         write(*,*)'zdtconduc,c,d',ilayer,zdtconduc(ig,ilayer),C(ilayer)
c     $,D(ilayer)
c         write(*,*) 'pt',ilayer,pt(ig,ilayer)
      ENDDO

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