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,ptime,zday,tsurf,zzlev,zzlay,
     &     pt,pq,pu,pv,pdt,pdq,
     $     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
      implicit none

      include "callkeys.h"

      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)
      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) :: zday
      REAL,INTENT(in) :: dist_sol
      REAL,INTENT(in) :: mu0(ngrid)
      REAL,INTENT(in) :: pq(ngrid,nlayer,nq)
      REAL,INTENT(in) :: ptimestep
      REAL,INTENT(in) :: ptime
      REAL,INTENT(in) :: tsurf(ngrid)
      REAL,INTENT(in) :: pu(ngrid,nlayer),pv(ngrid,nlayer)
      REAL,INTENT(in) :: pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)

      REAL,INTENT(out) :: zdteuv(ngrid,nlayer)
      REAL,INTENT(out) :: zdtconduc(ngrid,nlayer)
      REAL,INTENT(out) :: zdumolvis(ngrid,nlayer)
      REAL,INTENT(out) :: zdvmolvis(ngrid,nlayer)
      REAL,INTENT(out) :: zdqmoldiff(ngrid,nlayer,nq)
      REAL*8,INTENT(out) :: PhiEscH,PhiEscH2,PhiEscD

      INTEGER :: l,ig
      logical,save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
      
      if (firstcall) then
        ! default scheme for molecular diffusion: legacy
        moldiff_scheme=1
        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
      ! (tendencies are added later on, even if parametrization is not called)
      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,ptime,zday,pq,pdq,zdteuv)
      endif

      if (callconduct) THEN
        call conduction(ngrid,nlayer,ptimestep,pplay,pplev,pt,zdteuv,
     $                   tsurf,zzlev,zzlay,zdtconduc)
      endif

      if (callmolvis) THEN
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                zdteuv,zdtconduc,pu,
     $                   tsurf,zzlev,zzlay,zdumolvis)
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                zdteuv,zdtconduc,pv,
     $                   tsurf,zzlev,zzlay,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,
     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff,
     &                   PhiEscH,PhiEscH2,PhiEscD)
       else
        ! new MPF (modified pass flow) scheme
        call moldiff_MPF(ngrid,nlayer,nq,
     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff,
     &                   PhiEscH,PhiEscH2,PhiEscD)
       endif ! of if (moldiff_scheme==1)
      endif

      end subroutine thermosphere

      END MODULE thermosphere_mod