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)

      use conc_mod, only: rnew, cpnew
      USE comcstfi_h
      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 pplay(ngrid,nlayer)
      real pplev(ngrid,nlayer+1)
      REAL zzlay(ngrid,nlayer)
      real zzlev(ngrid,nlayer+1)
      REAL pt(ngrid,nlayer)
      real zday
      REAL dist_sol
      real mu0(ngrid)
      real pq(ngrid,nlayer,nq)
      real ptimestep
      real ptime
      real tsurf(ngrid)
      REAL pu(ngrid,nlayer),pv(ngrid,nlayer)
      REAL pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)

      REAL zdteuv(ngrid,nlayer)
      REAL zdtconduc(ngrid,nlayer)
      REAL zdumolvis(ngrid,nlayer)
      REAL zdvmolvis(ngrid,nlayer)
      real zdqmoldiff(ngrid,nlayer,nq)

      INTEGER l,ig
      logical,save :: firstcall=.true.

      if (firstcall) then
        if (.not. tracer) then
          do l=1,nlayer
            do ig=1,ngrid
              rnew(ig,l)=r
              cpnew(ig,l)=cpp
            enddo
          enddo
        endif
        firstcall= .false.
      endif

      if (calleuv) then
        zdteuv(1:ngrid,1:nlayer)=0
        call euvheat(ngrid,nlayer,nq,pt,pdt,pplev,pplay,zzlay,
     $               mu0,ptimestep,ptime,zday,pq,pdq,zdteuv)
      endif

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

      if (callmolvis) THEN
        zdumolvis(1:ngrid,1:nlayer)=0
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                zdteuv,zdtconduc,pu,
     $                   tsurf,zzlev,zzlay,zdumolvis)
        zdvmolvis(1:ngrid,1:nlayer)=0
        call molvis(ngrid,nlayer,ptimestep,pplay,pplev,pt,
     &                zdteuv,zdtconduc,pv,
     $                   tsurf,zzlev,zzlay,zdvmolvis)
      endif

      if (callmoldiff) THEN
        zdqmoldiff(1:ngrid,1:nlayer,1:nq)=0
        call moldiff_red(ngrid,nlayer,nq,
     &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff)
      endif

      return
      end