MODULE thermosphere_mod
      
      IMPLICIT NONE
      
      integer,save :: moldiff_scheme ! 1:legacy scheme, 2: MPF scheme
!$OMP THREADPRIVATE(moldiff_scheme)

      CONTAINS
      
      subroutine thermosphere(ngrid,nlayer,nq,
     &     pplev,pplay,dist_sol,
     $     mu0,ptimestep,zday,tsurf,zzlev,zzlay,
     &     pt,pq,pu,pv,pdt,pdq,pdu,pdv,
     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
     $     PhiEscH,PhiEscH2,PhiEscD)

      use ioipsl_getin_p_mod, only: getin_p
      use moldiff_red_mod, only: moldiff_red ! old molecular diffusion scheme
      use moldiff_MPF_mod, only: moldiff_MPF ! new molecular diffusion scheme
      use euvheat_mod, only: euvheat
      use conduction_mod, only: conduction
      use molvis_mod, only: molvis
      use conc_mod, only: rnew, cpnew
      USE comcstfi_h, only: r, cpp
      use mod_phys_lmdz_para, only : is_master
      use callkeys_mod, only: calleuv, callconduct, callmoldiff
      use callkeys_mod, only: callmolvis
      implicit none

      integer,intent(in) :: ngrid ! number of atmospheric columns
      integer,intent(in) :: nlayer ! number of atmospheric layers
      integer,intent(in) :: nq ! number of advected tracers
      REAL,INTENT(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
      REAL,INTENT(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
      REAL,INTENT(in) :: zzlay(ngrid,nlayer) ! altitude of mid-layer (m)
      REAL,INTENT(in) :: zzlev(ngrid,nlayer+1) ! altitude of layer boundaries (m)
      REAL,INTENT(in) :: pt(ngrid,nlayer) ! input temperature (K)
      REAL,INTENT(in) :: zday ! Mars date (sols since Ls=0)
      REAL,INTENT(in) :: dist_sol ! Sun-Mars distance (AU)
      REAL,INTENT(in) :: mu0(ngrid) ! cosine of solar zenith angle
      REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! input tracer mixing ratio (kg/kg_air)
      REAL,INTENT(in) :: ptimestep ! physics time step (s)
      REAL,INTENT(in) :: tsurf(ngrid) ! surface temperature (K)
      REAL,INTENT(in) :: pu(ngrid,nlayer) ! input zonal wind (m/s)
      REAL,INTENT(in) :: pv(ngrid,nlayer) ! input meridional wind (m/s)

      REAL,INTENT(out) :: zdteuv(ngrid,nlayer) ! tendency due to EUV (K/s)
      REAL,INTENT(out) :: zdtconduc(ngrid,nlayer) ! tendency due to conduction (K/s)
      REAL,INTENT(out) :: zdumolvis(ngrid,nlayer) ! tendency due to viscosity (m/s/s)
      REAL,INTENT(out) :: zdvmolvis(ngrid,nlayer) ! tendency due to viscosity (m/s/s)
      REAL,INTENT(out) :: zdqmoldiff(ngrid,nlayer,nq) ! tendency due to diffusion (kg/kg_air/s)
      REAL*8,INTENT(out) :: PhiEscH ! total escape flux of H (s-1)
      REAL*8,INTENT(out) :: PhiEscH2 ! total escape flux of H2 (s-1)
      REAL*8,INTENT(out) :: PhiEscD ! total escape flux of D (s-1)

      REAL,INTENT(inout) :: pdt(ngrid,nlayer) ! tendency on temperature (K/s)
      REAL,INTENT(inout) :: pdq(ngrid,nlayer,nq) ! tendency on tracers (kg/kg_air/s)
      REAL,INTENT(inout) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s/s)
      REAL,INTENT(inout) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s/s)

      INTEGER :: l,ig
      logical,save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
      
      if (firstcall) then
        ! default scheme for molecular diffusion: improved MPF scheme
        moldiff_scheme=2
        call getin_p("moldiff_scheme",moldiff_scheme)
        if (is_master) then
          write(*,*)"thermosphere: moldiff_scheme=",moldiff_scheme
        endif
        firstcall=.false.
      endif ! of if (firstcall)

      ! initialize tendencies to zero in all cases
      zdteuv(1:ngrid,1:nlayer)=0
      zdtconduc(1:ngrid,1:nlayer)=0
      zdumolvis(1:ngrid,1:nlayer)=0
      zdvmolvis(1:ngrid,1:nlayer)=0
      zdqmoldiff(1:ngrid,1:nlayer,1:nq)=0
      
      if (calleuv) then
        call euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay,
     $               mu0,ptimestep,zday,pq,pdq,zdteuv)
        ! update tendency on temperature
        pdt(:,:)=pdt(:,:)+zdteuv(:,:)
      endif

      if (callconduct) THEN
        call conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,pdt,
     $                   tsurf,zzlev,zzlay,zdtconduc)
        ! update tendency on temperature
        pdt(:,:)=pdt(:,:)+zdtconduc(:,:)
      endif

      if (callmolvis) THEN
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                pdt,pu,tsurf,zzlev,zzlay,zdumolvis)
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                pdt,pv,tsurf,zzlev,zzlay,zdvmolvis)
        ! update tendencies on winds
        pdu(:,:)=pdu(:,:)+zdumolvis(:,:)
        pdv(:,:)=pdv(:,:)+zdvmolvis(:,:)
      endif

      if (callmoldiff) THEN
       if (moldiff_scheme==1) then
        ! old "legacy" scheme
        call moldiff_red(ngrid,nlayer,nq,
     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     &                   zdqmoldiff,PhiEscH,PhiEscH2,PhiEscD)
       else
        ! new MPF (modified pass flow) scheme
        call moldiff_MPF(ngrid,nlayer,nq,
     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     &                   zdqmoldiff,PhiEscH,PhiEscH2,PhiEscD)
       endif ! of if (moldiff_scheme==1)

       ! update tendencies on tracers
       pdq(:,:,:)=pdq(:,:,:)+zdqmoldiff(:,:,:)
       
      endif ! of if (callmoldiff)

      end subroutine thermosphere

      END MODULE thermosphere_mod