SUBROUTINE euvheat(pt,pdt,pplev,pplay,zzlay,dist_sol,
     $     mu0,ptimestep,ptime,zday,pq,pdq,pdteuv)

       IMPLICIT NONE
c=======================================================================
c   subject:
c   --------
c   Computing heating rate due to EUV absorption
c
c   author:  MAC 2002
c   ------
c
c   input:
c   ----- 
c   dist_sol		  sun-Mars distance (AU)
c   mu0(ngridmx)           
c   pplay(ngrid,nlayer)   pressure at middle of layers (Pa)
c
c   output:
c   -------
c
c   pdteuv(ngrid,nlayer)      Heating rate (K/s)
c
c=======================================================================
c
c    0.  Declarations :
c    ------------------
c
#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "comdiurn.h"
#include "param.h"
#include "param_v3.h"
#include "chimiedata.h"
#include "tracer.h"
#include "conc.h"
c-----------------------------------------------------------------------
c    Input/Output
c    ------------

      REAL pplay(ngridmx,nlayermx)
      REAL zzlay(ngridmx,nlayermx)
      real pplev(ngridmx,nlayermx+1)
      REAL pt(ngridmx,nlayermx)
      REAL pdt(ngridmx,nlayermx)
      real zday
      REAL dist_sol
      real mu0(ngridmx)
      real pq(ngridmx,nlayermx,nqmx)
      real pdq(ngridmx,nlayermx,nqmx)
      real ptimestep,ptime


      REAL pdteuv(ngridmx,nlayermx)
c
c    Local variables :
c    -----------------
      integer nespeuv
      parameter (nespeuv=6)

      INTEGER l,iq,ig,n,inif ! , ngrid, nlayer
      real rm(nlayermx,nespeuv)		! number density (cm-3)
      real zq(ngridmx,nlayermx,nqmx) ! local updated tracer quantity
      real zt(ngridmx,nlayermx) ! local updated atmospheric temperature
      real zlocal(nlayermx)
      real zenit
      real aux1(nlayermx)
      real aux2(nlayermx)
      real jtot(nlayermx)
      real dens						! amu/cm-3
      real tx(nlayermx)
      real tmean
      
! tracer indexes for the EUV heating:
      integer,parameter :: i_co2=1
      integer,parameter :: i_o2=2
      integer,parameter :: i_o=3
      integer,parameter :: i_h2=4
      integer,parameter :: i_h2o=5
      integer,parameter :: i_h2o2=6
      
! Tracer indexes in the GCM:
      integer,save :: g_co2=0
      integer,save :: g_o=0
      integer,save :: g_o2=0
      integer,save :: g_h2=0
      integer,save :: g_h2o2=0
      integer,save :: g_h2o=0


      logical,save :: firstcall=.true.

! Initializations and sanity checks:

      if (firstcall) then
        if ((nespeuv.gt.nqmx).or.(nespeuv.gt.ncomp)) then
          print*,'******* Dimension problem in EUVHEAT ********'
          STOP
        endif

        ! identify the indexes of the tracers we'll need
        g_co2=igcm_co2
        if (g_co2.eq.0) then
          write(*,*) "euvheat: Error; no CO2 tracer !!!"
          stop
        endif
        g_o=igcm_o
        if (g_o.eq.0) then
          write(*,*) "euvheat: Error; no O tracer !!!"
          stop
        endif
        g_o2=igcm_o2
        if (g_o2.eq.0) then
          write(*,*) "euvheat: Error; no O2 tracer !!!"
          stop
        endif
        g_h2=igcm_h2
        if (g_h2.eq.0) then
          write(*,*) "euvheat: Error; no H2 tracer !!!"
          stop
        endif
        g_h2o2=igcm_h2o2
        if (g_h2o2.eq.0) then
          write(*,*) "euvheat: Error; no H2O2 tracer !!!"
          stop
        endif
        g_h2o=igcm_h2o_vap
        if (g_h2o.eq.0) then
          write(*,*) "euvheat: Error; no water vapor tracer !!!"
          stop
        endif
        
        firstcall= .false.
!        print*,'EUV',nlayer,nlayermx,ngrid,ngridmx   
      endif ! of if (firstcall)

c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      ! build local updated values of tracers and temperature
      do l=1,nlayermx
        do ig=1,ngridmx
          ! chemical species
          zq(ig,l,g_co2)=pq(ig,l,g_co2)+pdq(ig,l,g_co2)*ptimestep
          zq(ig,l,g_o2)=pq(ig,l,g_o2)+pdq(ig,l,g_o2)*ptimestep
          zq(ig,l,g_o)=pq(ig,l,g_o)+pdq(ig,l,g_o)*ptimestep
          zq(ig,l,g_h2)=pq(ig,l,g_h2)+pdq(ig,l,g_h2)*ptimestep
          zq(ig,l,g_h2o2)=pq(ig,l,g_h2o2)+pdq(ig,l,g_h2o2)*ptimestep
          zq(ig,l,g_h2o)=pq(ig,l,g_h2o)+pdq(ig,l,g_h2o)*ptimestep
          ! atmospheric temperature
          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
        enddo
      enddo

      call flujo(solarcondate)
c     call flujo(solarcondate+zday/365.) ! version with EUV time evolution

      do ig=1,ngridmx

        zenit=acos(mu0(ig))*180./acos(-1.)
 
        do l=1,nlayermx
          dens=pplay(ig,l)/(rnew(ig,l)*zt(ig,l)) / 1.66e-21
          rm(l,i_co2)  = zq(ig,l,g_co2)  * dens / mmol(g_co2)
          rm(l,i_o2)   = zq(ig,l,g_o2)   * dens / mmol(g_o2)
          rm(l,i_o)    = zq(ig,l,g_o)    * dens / mmol(g_o)
          rm(l,i_h2)   = zq(ig,l,g_h2)   * dens / mmol(g_h2)
          rm(l,i_h2o)  = zq(ig,l,g_h2o)  * dens / mmol(g_h2o)
          rm(l,i_h2o2) = zq(ig,l,g_h2o2) * dens / mmol(g_h2o2)
        enddo

c        zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))
c     &            *Rnew(ig,1)*zt(ig,1)/g
        zlocal(1)=zzlay(ig,1)
        zlocal(1)=zlocal(1)/1000.
        tx(1)=zt(ig,1)

        do l=2,nlayermx
          tx(l)=zt(ig,l)
          tmean=tx(l)
          if(tx(l).ne.tx(l-1))
     &    tmean=(tx(l)-tx(l-1))/log(tx(l)/tx(l-1))
c          zlocal(l)= zlocal(l-1)
c     &   -log(pplay(ig,l)/pplay(ig,l-1))*Rnew(ig,l-1)*tmean/g/1000.
        zlocal(l)=zzlay(ig,l)/1000.
        enddo

        call hrtherm (rm(1,i_co2),rm(1,i_o2),rm(1,i_o),
     &                rm(1,i_h2),rm(1,i_h2o),rm(1,i_h2o2),
     &                aux1,aux2,tx,nlayermx,zlocal,
     &                solarcondate+zday/365.,zenit,jtot)

        do l=1,nlayermx
          pdteuv(ig,l)=0.16*jtot(l)/10.
     &            /(cpnew(ig,l)*pplay(ig,l)/(rnew(ig,l)*zt(ig,l)))
     &                  *(1.52/dist_sol)**2 

        enddo	

      enddo  ! of do ig=1,ngridmx

      return
      end