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

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

      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"

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(ngrid,nlayer)  ! wind
      REAL zdvelmolvis(ngrid,nlayer)
c   local:
c   ------

      INTEGER ilayer,ig
      REAL zvel(nlayer)
      REAL zt(nlayer)
      REAL alpha(nlayer)
      REAL lambda(nlayer),muvol(nlayer)
      REAL C(nlayer),D(nlayer),den(nlayer)
      REAL pdvelm(nlayer),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) / den(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