SUBROUTINE molvis(ngrid,nlayer,ptimestep,
     $                   pplay,pt,pzlay,pzlev,
     $                   pdtconduc,pvel,tsurf,zdvelmolvis)
      IMPLICIT NONE

c=======================================================================
c
c   Molecular Viscosity Diffusion
c   
c   Based on conduction.F  by N. Descamp, F. Forget 05/1999
c
c   modified by M. Angelats i Coll
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),pdtconduc(ngrid,nlayer)
      REAL pzlay(ngrid,nlayer),pzlev(ngrid,nlayer+1),pplay(ngrid,nlayer)
      REAL pvel(ngridmx,nlayermx)  ! wind
      REAL zdvelmolvis(ngridmx,nlayermx)
c   local:
c   ------

      INTEGER ilayer,ig
      REAL zvel(nlayermx)
      REAL zt(nlayermx)
      REAL alpha(nlayermx)
      REAL lambda(nlayermx),muvol(nlayermx)
      REAL C(nlayermx),D(nlayermx),den(nlayermx)
      REAL pdvelm(nlayermx),tsurf      
      REAL fac, m

c   constants used locally
c    ---------------------
c     The atmospheric conductivity is a function of temperature T :
c      conductivity = Akk* T**skk
c     Molecular viscosity is related to thermal conductivity by:
c       conduc = 0.25*(9*gamma - 5)* Cv * molvis
c     where gamma = Cp/Cv.  For dry air.

      real,parameter :: Akk=5.63E-05  
      real,parameter :: skk=1.12
!      real,parameter :: Akk=2.6E-05  
!      real,parameter :: skk=1.3
      real,parameter :: velsurf=0.0

      logical firstcall
      save firstcall
      data firstcall /.true./
c-----------------------------------------------------------------------
c   calcul des coefficients alpha et lambda
c
c-----------------------------------------------------------------------

      IF (firstcall) THEN
        write (*,*) 'molvis : coeff of molecular viscosity skk'
        write(*,*) skk
        firstcall = .false.
      END IF

      DO ig=1,ngrid  

        zt(1)=pt(ig,1)+pdtconduc(ig,1)*ptimestep
        zvel(1)=pvel(ig,1)  
  
        DO ilayer = 2 , nlayer
           
          zt(ilayer)=pt(ig,ilayer)+pdtconduc(ig,ilayer)*ptimestep
          zvel(ilayer)=pvel(ig,ilayer)
        ENDDO
    
        fac=0.25*(9.*cpp-5.*(cpp-r)) 
        lambda(1) = Akk*tsurf**skk / pzlay(ig,1)/fac  

         DO ilayer = 2 , nlayer
           
          fac=(9.*cpp-5.*(cpp-r))/4.
          lambda(ilayer) = Akk/fac * zt(ilayer)**skk /
     $                      (pzlay(ig,ilayer)-pzlay(ig,ilayer-1)) 
        ENDDO

        DO ilayer=1,nlayer-1
          muvol(ilayer)=pplay(ig,ilayer)/(r*zt(ilayer)) 
          alpha(ilayer)=(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*zt(nlayer)) 
!OLD
!          alpha(nlayer)=(muvol(nlayer)/ptimestep)*
!     $                      2*(pzlay(ig,nlayer)-pzlev(ig,nlayer))  
c!!!!! <----
c!!!!! alerte !!!!!c
c        write(*,*) lambda(1),muvol(1),tsurf,pt(1,1)   

!NEW TB16:
          
          alpha(nlayer)=(muvol(nlayer)/ptimestep)*
     $                      (pzlev(ig,nlayer)+10000. 
     &                               -pzlev(ig,nlayer))


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

      den(1)=alpha(1)+lambda(2)+lambda(1)
      C(1)=lambda(1)*(velsurf-zvel(1))+lambda(2)*(zvel(2)-zvel(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)*(zvel(ilayer+1)-zvel(ilayer)) 
     $   +lambda(ilayer)*(zvel(ilayer-1)-zvel(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)+zvel(nlayer-1)-zvel(nlayer) 
      C(nlayer)=(C(nlayer)*lambda(nlayer)+phitop(ig)) / den(nlayer) 
       	 	
c----------------------------------------------------------------------
c
c      calcul de la nouvelle temperature pdvelm
c
c----------------------------------------------------------------------

      DO ilayer=1,nlayer
         pdvelm(ilayer)=0. 
      ENDDO
    
         pdvelm(nlayer)=C(nlayer)
!         pt(nlayer)=ttop
!         write(*,*)'pt',pt(nlayer)
         
      DO ilayer=nlayer-1,1,-1
         pdvelm(ilayer)=C(ilayer)+D(ilayer)*pdvelm(ilayer+1)     
      ENDDO        

c-----------------------------------------------------------------------
c
c     calcul de la tendance zdvelmolvis
c
c-----------------------------------------------------------------------
    
      DO ilayer=1,nlayer
         zdvelmolvis(ig,ilayer) = pdvelm(ilayer) / ptimestep
      ENDDO
      ENDDO             ! boucle sur ngrid

      RETURN
      END