MODULE molvis_mod
      
IMPLICIT NONE
      
CONTAINS

      SUBROUTINE molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,pdt, &
                        pvel,zzlev,zzlay,zdvelmolvis)
      
      use comcstfi_mod, only: r, cpp, mugaz
      use callkeys_mod, only: phitop_molvis,zztop
      use conc_mod,     only: lambda
      use gases_h

!=======================================================================
!
!   Molecular Viscosity Diffusion
!  
!   Based on the molvis.F module from the Mars PCM
!   Adapted by Maxime Maurice (2025)  
!   The kiunematic viscosity mu (kg/m/s) is calculated from
!   the thermal conductivity lambda (m/s^2) using the formula:
!   lambda / mu = (9*Cp + 5*Cv) / 4
!               = (9*Cp + 5*(Cp - R)) / 4
!
!=======================================================================

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

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

      integer,intent(in) :: ngrid ! number of atmospheric columns
      integer,intent(in) :: nlayer ! number of atmospheric layers
      REAL,intent(in)    :: ptimestep
      REAL,intent(in)    :: pplay(ngrid,nlayer)
      REAL,intent(in)    :: pplev(ngrid,nlayer+1)
      REAL,intent(in)    :: zzlay(ngrid,nlayer)
      REAL,intent(in)    :: zzlev(ngrid,nlayer+1)
      real,intent(in)    :: pt(ngrid,nlayer)
      REAL,intent(in)    :: pvel(ngrid,nlayer)
      real,intent(in)    :: pdt(ngrid,nlayer)

      real,intent(out)   :: zdvelmolvis(ngrid,nlayer)

!   local:
!   ------

      INTEGER l,ig, nz     ! layer index, grid index, and number of layers
      real fac             ! conversion factor between thermal conductivity and molecular viscosity
      REAL zvel(nlayer)    ! wind
      real zt(nlayer)      ! temperature (K)
      REAL alpha(nlayer)   ! physical coefficient
      REAL mu(nlayer)      ! dynamical viscosity (kg/m/s)
      real muvol(nlayer)   ! molecular density (m^-3)
      REAL C(nlayer)       ! Finite-difference coefficient 1
      real D(nlayer)       ! Finite-difference coefficient 2
      real den(nlayer)     ! density
      REAL pdvelm(nlayer)  ! wind tendency
      REAL zlay(nlayer)
      real zlev(nlayer+1)
      REAL,PARAMETER :: velsurf =0.0

!-----------------------------------------------------------------------
!   calcul des coefficients alpha et mu
!-----------------------------------------------------------------------

      nz = nlayer

      fac = (9 * cpp - 5 * (cpp - r)) / 4

      do ig=1,ngrid

        zt(1)   = pt(ig,1) + pdt(ig,1) * ptimestep
        zvel(1) = pvel(ig,1)
        zlay(1) = zzlay(ig,1)
        zlev(1) = zzlev(ig,1)

        do l=2,nz
          zt(l)   = pt(ig,l) + (pdt(ig,l) + pdt(ig,l)) * ptimestep
          zvel(l) = pvel(ig,l)
          zlay(l) = zzlay(ig,l)
          zlev(l) = zzlev(ig,l)
        enddo

        zlev(nz+1) = zlev(nz) + zztop
        
        mu(1) = lambda(ig,1) / fac

        DO l=2,nz
          mu(l)=lambda(ig,l)/fac
        ENDDO
    
        DO l=1,nz-1
          muvol(l) = pplay(ig,l) / (r*zt(l)) 
          alpha(l) = (muvol(l) / ptimestep) * (zlev(l+1) - zlev(l))
        ENDDO

        muvol(nz) = pplay(ig,nz) / (r*zt(nz)) 
        alpha(nz) = (muvol(nz) / ptimestep) * (zlev(nz+1) - zlev(nz))

!--------------------------------------------------------------------
!
!     C and D coefficients
!
!-------------------------------------------------------------------

        den(1) = alpha(1) + mu(2) + mu(1)
        C(1)   = mu(1) * (velsurf - zvel(1)) + mu(2) * (zvel(2) - zvel(1))
        C(1)   = C(1) / den(1)	     
        D(1)   = mu(2) / den(1)           
   
        DO l = 2, nz-1
           den(l) = alpha(l) + mu(l+1)
           den(l) = den(l) + mu(l) * (1 - D(l-1))
         
           C(l) = mu(l+1) * (zvel(l+1) - zvel(l)) + mu(l) * (zvel(l-1) - zvel(l) + C(l-1))    
           C(l) = C(l) / den(l)           

           D(l) = mu(l+1) / den(l)
        ENDDO 

        den(nz) = alpha(nz) + mu(nz) * (1 - D(nz-1))
        C(nz)   = C(nz-1) + zvel(nz-1) - zvel(nz) 
        C(nz)   = (C(nz) * mu(nz) + phitop_molvis) / den(nz) 
       	 	
!----------------------------------------------------------------------
!
!      Wind tendency calculation
!
!-----------------------------------------------------------------------
    
        DO l=1,nz
          pdvelm(l) = 0.
        ENDDO
          pdvelm(nz) = C(nz)
        DO l=nz-1,1,-1
          pdvelm(l) = C(l) + D(l) * pdvelm(l+1)
        ENDDO 

        DO l=1,nz
          zdvelmolvis(ig,l) = pdvelm(l) / ptimestep
        ENDDO

      ENDDO             ! boucle sur ngrid

      END SUBROUTINE molvis

      END MODULE molvis_mod