      module qrj_module
      use params_module,only: nlevp1,nlonp4
      implicit none
!     -------------------------------------------------------------
!     Titan TGCM version : latest mods (02/15/05) S. W. BOUGHER
!     Titan TGCM version : latest mods (03/24-29/05) S. W. BOUGHER
!     Titan TGCM version : latest mods (04/06-07/05) S. W. BOUGHER
!     (check units of qtotal: erg/gm/sec))
!     (working : 04/07/05): S. W. BOUGHER
!     -- lmax = 59 (7-SRC,7-SRB,1-LYA,37-EUV,7-XRAY). [1-1750A only]
!     -- lmax = 59 (Solar cycle variability included)
!     -- lmax = 16-59 (1050-1A   : EUV-XRAY)
!     -- lmax = 1-15  (1750-1050 : SRB-SRC)
!     -------------------------------------------------------------
!     -- X-Sections for N2(A), IN2(tot) from tables. 
!     -- X-Sections for CH4(A),  ICH4(tot) from tables. 
!     -- X-Sections for H2(A), IH2(tot) from tables. 
!     -- Branchings for Ich4-ch4p, Ich4-ch3p, and In2-np from tables. 
!     -------------------------------------------------------------
!
! Coefficients for wavelength scans in qrj.
! These are referenced in routines qrj, init_sflux, init_sigmas, and 
! set_qrj_coeff (all in qrj_module.F).
!
      integer,parameter :: lmax=59, l1=16
      real,parameter :: sfactor=0.01099
      real :: 
     |  euveff(nlevp1),  !
     |  srceff(nlevp1),  !
     |  sigeuv(3,lmax),  !
     |  rlmeuv(lmax),    !
     |  feuv(lmax),      !
     |  fsrc(l1-1),      !
     |  sigsrc(l1-1),    !
     |  rlmsrc(l1-1),    !
     |  sigmas(3,lmax),  !
     |  sflux(lmax),     !
     |  brnp(lmax),      !
     |  brch4p(lmax),    !
     |  brch3p(lmax)     !
!
! For euvac, init_sigmas
      integer,parameter :: neuv=37
      real,dimension(neuv) ::
     |  wleuv1,wleuv2,
     |  sigan2, sigin2,
     |  brn2np
!
! Heating and ionization terms set by qrj, and used by other routines.
! These are allocated for task subdomains by alloc_q (called from allocdata)
!
      real,dimension(:,:,:),allocatable ::
     |  rj,     ! solar photodissociation of ch4  (#/cm3/sec)
     |  qtef,   ! solar photodissociation of n2
     |  qtotal, ! total heating (erg/gm/sec)
     |  qch4p,  ! ch4+ (#/cm3/sec)
     |  qch3p,  ! ch3+
     |  qn2p,   ! n2+
     |  qnp     ! n+
!
      contains
!-----------------------------------------------------------------------
      subroutine qrj(scch4,sch2,scn2,tn,ch4,h2,xnmbari,fcp,
     |  lev0,lev1,lon0,lon1,lat)
!
!    13-arguments to qrj call-- titan
!
! Calculate heating and dissociation rates.
!
      use input_module,only: f107, f107a
      use init_module,only: sfeps,idn
      use cons_module,only: t0,expz,avo,
     |  rmassinv_ch4,rmassinv_h2,rmassinv_n2,expzmid,gask,
     |  expzmid_inv,p0,check_exp
      use bndry_module,only: fb,b
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) ::
     |  scch4,sch2,scn2,  ! chapman integrals
     |  tn,ch4,h2,        ! tn and species mass mixing ratios
     |  xnmbari,          ! p0*e(-z)*barm/kT at interfaces
     |  fcp               ! cp from cpktkm at interfaces
!
! VT vampir tracing:
!
#ifdef VT
#include "VT.inc"
#endif
!
! Local:
      integer :: k,i,l,nlevs,ier
!     New SRB Absorption by CH4 for titans:
      real,parameter ::  hc = 1.9845E-16
      real :: rlmeuvinv(lmax),rlmsrcinv(l1-1),factor,fmin,fmax
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  r,      ! s14
     |  rp,     ! s15
     |  taur,   ! s1
     |  tauq,   ! s2
     |  tau1,   ! s3
     |  tau2,   ! s4
     |  tau3,   ! s5
     |  etaur,  ! s6
     |  etauq   ! s7
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  ch4i,   ! ch4 at interfaces (s1)
     |  h2i,    ! h2  at interfaces (s2)
     |  n2i,    ! n2  at interfaces (s3)
     |  tni     ! tn  at interfaces (s6)
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  zch4,      ! ch4 total ionization (s8)
     |  zn2,       ! n2  total ionization (s10)
     |  zch4_brch4p,  ! ch4  ionization * brch4p (s11)
     |  zch4_brch3p,  ! ch4  ionization * brch3p (s12)
     |  zn2_brn2p,  ! n2  ionization * (1-brnp) (s13)
     |  zn2_brnp,  !  n2  ionization * brnp  (s13)
     |  sigchap,   ! (s9)
     |  fexp,      ! 
     |  qeuv,      ! QEUV (erg/cm3/sec)
     |  qsrc,      ! QSRC (erg/cm3/sec)
     |  qtot       ! QTOT (erg/cm3/sec)
      real,dimension(lev0:lev1,lon0:lon1) ::
     |  sum1,   ! sum(ch4,h2,n2)(Fluxes)   
     |  sum2,   ! sum(ch4,h2,n2)((sigma*psi/rmass)
     |  sum3,   ! sum(ch4)(sigma*Fluxes)   
     |  sum4    ! sum(ch4)(sigma*Fluxes*energy) 
      logical,parameter :: 
     |  debug=.true.     ! insert print statements
!
! expo() (util.F) is used only if check_exp is true. This will avoid
! NaNS fpe, but will degrade performance. Check_exp is in cons.F.
!
      real,external :: expo ! used only when check_exp is set (util.F)
!
!     if (debug) write(6,"('Enter qrj: lat=',i3,' lon0,1=',2i3)") 
!    |  lat,lon0,lon1
!
#ifdef VT
      code = 118 ; state = 'qrj' ; activity='ModelCode'
      call vtbegin(118,ier)
#endif
!
! Exec:
      nlevs = lev1-lev0+1

!     write(6,"('qrj: lat=',i2,' lon0,1=',2i3,' idn(lon0:lon1)=',/,
!    |  (15i3))") lat,lon0,lon1,idn(lon0:lon1)
!
! --------------------------------------------------------------------
!     call addfsech('SCCH4',' ',' ',scch4(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('SCH2',' ',' ',sch2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('SCN2',' ',' ',scn2(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('XNMBARI',' ',' ',xnmbari(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
! --------------------------------------------------------------------
!
! Will multiply by inverses:
      do i=1,lmax
        rlmeuvinv(i) = 1./rlmeuv(i)
      enddo
      do i=1,l1-1 ! 1,15
        rlmsrcinv(i) = 1./rlmsrc(i)
      enddo
!
! COMPUTE S14 = RSQ, S15 = RSP
! S1 = TAU(R),  S2 = TAU(Q)
!
!     write(6,"('qrj: lat=',i2,' lon0,1=',2i3,' idn(lon0:lon1)=',/,
!    |  (12i3))") lat,lon0,lon1,idn(lon0:lon1)
!
! ---------------------------------------------------------------
! Optical Depth Calculations: Is this a PE enhancement scheme?
! Assess Carefully Before Using.
! ---------------------------------------------------------------
      do i=lon0,lon1
        do k=lev0,lev1
        taur(k,i) = sigeuv(1,49)*scch4(k,i)+sigeuv(2,49)*sch2(k,i)+
     |                sigeuv(3,49)*scn2(k,i)
        tauq(k,i) = sigeuv(1,20)*scch4(k,i)+sigeuv(2,20)*sch2(k,i)+
     |                sigeuv(3,20)*scn2(k,i)
          tau1(k,i) = 1.3*taur(k,i)
          tau2(k,i) = 2.0*taur(k,i)
          tau3(k,i) = 2.5*taur(k,i)
          if (taur(k,i) > 9.) taur(k,i) = 9.
          if (tauq(k,i) > 9.) tauq(k,i) = 9.
          if (tau1(k,i) > 9.) tau1(k,i) = 9.
          if (tau2(k,i) > 9.) tau2(k,i) = 9.
          if (tau3(k,i) > 9.) tau3(k,i) = 9.
!
! TAU(N) = EXP(-TAU(N)) FOR N = 1,3,1
!
          tau1(k,i) = exp(-tau1(k,i))
          tau2(k,i) = exp(-tau2(k,i))
          tau3(k,i) = exp(-tau3(k,i))
!
! taur = EXP(-TAU(R)), tauq = EXP(-TAU(Q))
!
          etaur(k,i) = exp(-taur(k,i))
          etauq(k,i) = exp(-tauq(k,i))
!
! r = r(n2p) :  Are these enchancements for PE yields?
!
          r(k,i) = etaur(k,i)+2.*(tau1(k,i)+tau2(k,i)+tau3(k,i))
          r(k,i) = 2.4*etaur(k,i)/r(k,i)
        enddo
      enddo
!     if (debug) write(6,"(/,'qrj after tau: lat=',i3)") lat
!
!     call addfsech('TAUQ',' ',' ',tauq,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('TAUR',' ',' ',taur,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('R'   ,' ',' ',r   ,lon0,lon1,nlevs,nlevs,lat)
! ---------------------------------------------------------------
!
! ch4,h2 at interface levels: mass mixing ratios
      do i=lon0,lon1
        do k=lev0,lev1-1
          ch4i (k+1,i) = 0.5*(ch4(k,i)+ch4(k+1,i))
          h2i (k+1,i) = 0.5*(h2(k,i)+h2(k+1,i))
        enddo
      enddo
!
! Bottom boundary: mass mixing ratios
      do i=lon0,lon1
        ch4i(1,i) = .5*((b(i,1,1)+1.)*ch4(1,i)+b(i,1,2)*h2(1,i)+fb(i,1))
        h2i(1,i) = .5*(b(i,2,1)*ch4(1,i)+(b(i,2,2)+1.)*h2(1,i)+fb(i,2))
      enddo
!
! N2 at interfaces:  mass mixing ratios
      do k=lev0,lev1
        n2i(k,:) = 1.-ch4i(k,:)-h2i(k,:)
      enddo 
!     call addfsech('CH4I'  ,' ',' ',ch4i,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('H2I'  ,' ',' ',h2i,lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('N2I'  ,' ',' ',n2i,lon0,lon1,nlevs,nlevs,lat)

!     if (debug) write(6,"(/,'qrj after interfaces: lat=',i3)") lat
!
! ---------------------------------------------------------------------
! Initialize heating, I, J terms on current processor domain: 7-fields
! (global module data for use by other routines)
      rj    (lev0:lev1,lon0:lon1,lat) = 0.
      qtef  (lev0:lev1,lon0:lon1,lat) = 0.  
      qtotal(lev0:lev1,lon0:lon1,lat) = 0.
      qch4p (lev0:lev1,lon0:lon1,lat) = 0.
      qch3p (lev0:lev1,lon0:lon1,lat) = 0.
      qn2p  (lev0:lev1,lon0:lon1,lat) = 0.
      qnp   (lev0:lev1,lon0:lon1,lat) = 0.
!
! Initialize ionization rates: 6-fields
! (array ops on local arrays dimensioned (lon0:lon1,lev0:lev1)):
      zch4 (lev0:lev1,lon0:lon1) = 0.          ! ch4 total ionization 
      zn2 (lev0:lev1,lon0:lon1) = 0.           ! n2  total ionization 
      zch4_brch4p (lev0:lev1,lon0:lon1) = 0.   ! ch4  ionization * brch4p 
      zch4_brch3p (lev0:lev1,lon0:lon1) = 0.   ! ch4  ionization * brch3p
      zn2_brnp (lev0:lev1,lon0:lon1) = 0.      ! n2  ionization * brnp 
      sigchap (lev0:lev1,lon0:lon1) = 0.       ! SRC intermediate product
      fexp (lev0:lev1,lon0:lon1) = 0.          ! SRC intermediate product

!     if (debug) write(6,"(/,'qrj after initializations: lat=',i3)") lat
!
! ---------------------------------------------------------------------
!
! Wavelength loop --------------------------------------------------
! Summation over wavelength: soft-XRAY and EUV bins (44)
! Multiply Q by (1/r**2) mean heliocentric distance of Saturn 
!
      do l=l1,lmax     ! 16,59 (44 wavelengths) [1050A-1A]
        sum1(:,:) = 0. ! sum(ch4,h2,n2)(Fluxes). FLUXES.
        sum2(:,:) = 0. ! sum(ch4,h2,n2)(sigma*psi/rmass). DENSITIES.
        do i=lon0,lon1
          do k=lev0,lev1
            sum1(k,i) = sum1(k,i)+sigeuv(1,l)*scch4(k,i)+
     |                            sigeuv(2,l)*sch2(k,i)+
     |                            sigeuv(3,l)*scn2(k,i)
!
!           if (.not.check_exp) then
!             sum1(k,i) = feuv(l)*exp(-sum1(k,i))
!           else
!             sum1(k,i) = feuv(l)*expo(-sum1(k,i),0)
!           endif
              sum1(k,i) = feuv(l)*exp(-sum1(k,i))*sfactor

        sum2(k,i) = sum2(k,i)+sigeuv(1,l)*ch4i(k,i)*rmassinv_ch4+
     |                          sigeuv(2,l)*h2i(k,i)*rmassinv_h2+
     |                          sigeuv(3,l)*n2i(k,i)*rmassinv_n2
          enddo
        enddo
! -------------------------------------------------------------------
! Longitude and column domain of current process:
! Longitude-Column Loop
! -------------------------------------------------------------------
        do i=lon0,lon1
          do k=lev0,lev1
            qtotal(k,i,lat) = qtotal(k,i,lat)+hc*rlmeuvinv(l)*
     |        sum1(k,i)*sum2(k,i)
!
! Add ionization for ch4 and n2: total values!
            zch4(k,i) = zch4(k,i)+sigmas(1,l)*sum1(k,i)*ch4i(k,i)*
     |        rmassinv_ch4
            zn2(k,i) = zn2(k,i)+sigmas(3,l)*sum1(k,i)*n2i(k,i)*
     |        rmassinv_n2
!
! Photo dissoc of n2:
            qtef(k,i,lat) = qtef(k,i,lat)+(sigeuv(3,l)-sigmas(3,l))*
     |        sum1(k,i)*n2i(k,i)*rmassinv_n2
!
! Photo dissoc of ch4: Must add SRC component also!
            rj(k,i,lat) = rj(k,i,lat)+(sigeuv(1,l)-sigmas(1,l))*
     |        sum1(k,i)*ch4i(k,i)*rmassinv_ch4
!
! Partitioning Ionization for  ch4(2) and n2(2):
!  ** CH4 paths for ionization ------------------------------
            zch4_brch4p(k,i) = zch4_brch4p(k,i)+
     |          sigmas(1,l)*sum1(k,i)*ch4i(k,i)*rmassinv_ch4*
     |          brch4p(l)
            zch4_brch3p(k,i) = zch4_brch3p(k,i)+
     |          sigmas(1,l)*sum1(k,i)*ch4i(k,i)*rmassinv_ch4*
     |          brch3p(l)
!  ** IN2 path for N+ ionization ----------------------------
            zn2_brnp(k,i) = zn2_brnp(k,i)+
     |          sigmas(3,l)*sum1(k,i)*n2i(k,i)*rmassinv_n2*
     |          brnp(l)
          enddo
        enddo
! ------------------------------------------------------------
!     Wavelength loop1 ----------------------------------------
      enddo ! l=l1,lmax
!     if (debug) write(6,"(/,'qrj after wavelength loop1: lat=',i3)") 
!    | lat
! -------------------------------------------------------------------
!
! Multiply Q by efficiency factor: Titan value from Fox (1991)
! -- Redo entire efficiency calculation (De La Haye and Bell, 2005)
! -- Check Against 1-D and 3-D Model Heating (qeuv) : erg/cm3/sec
!
      do i=lon0,lon1
        do k=lev0,lev1
          qtotal(k,i,lat) = qtotal(k,i,lat)*euveff(k)*avo
          qeuv(k,i)=qtotal(k,i,lat)*xnmbari(k,i)/avo
        enddo
      enddo
!
!    1-1050A Neutral Heating at Titan (erg/cm3/sec): Scaled Incl. Efficiency
!     call addfsech('Q1',' ',' ',qeuv(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!
!  --------------------------------------------------------------
! Final Contributions to ICH4+, ICH3+, IN2+, IN+, JN2:
! How to add PE contributions to N2+ and N+?
!  --------------------------------------------------------------
      do i=lon0,lon1
        do k=lev0,lev1
!  JN2 Dissociation Rate 
          qtef(k,i,lat) = qtef(k,i,lat)*xnmbari(k,i)
!  IN2 to N2+ Ionization rate
!         qn2p(k,i,lat) = (zn2(k,i)*(1.+r(k,i))-
!    |      zn2_brnp(k,i))*xnmbari(k,i) 
          qn2p(k,i,lat) = (zn2(k,i)-zn2_brnp(k,i))*
     |      xnmbari(k,i) 
!  IN2 to N+ Ionization rate
          qnp(k,i,lat)  = zn2_brnp(k,i)*xnmbari(k,i)
!  ICH4 to CH4+ Ionization rate
          qch4p(k,i,lat) = zch4_brch4p(k,i)*xnmbari(k,i)
!  ICH4 to CH3+ Ionization rate
          qch3p(k,i,lat) = zch4_brch3p(k,i)*xnmbari(k,i)
        enddo
      enddo
!  --------------------------------------------------------------
! 
! tn at interfaces:
      do i=lon0,lon1
        tni(1,i) = tn(lev1,i) ! tn bottom boundary in top slot
        do k=lev0+1,lev1-1
          tni(k,i) = .5*(tn(k-1,i)+tn(k,i))
        enddo
        tni(lev1,i) = tn(lev1-1,i) ! nlevp1 <- nlev
      enddo
!
!---------------------------------------------------------------
! Wavelength loop 
! Summation over wave length: 1750-1050A (SRB, LYA, SRC Start)
!---------------------------------------------------------------
! Fluxes at Titan are scaled by mean Saturn heliocentric distance
! Multiply Q by (1/r**2) mean heliocentric distance of Saturn 
! Double Check this Scheme before proceeding to 1-D validation
!
!  Initializations outside lambda loop for full SRC region summation
!
        sum3(:,:) = 0. ! sum(ch4)(sigma*F). SIGMA*FLUXES.
        sum4(:,:) = 0. ! sum(ch4)(sigma*F*hv). SIGMA*FLUXES*ENERGY
!
      do l=1,l1-1 ! 1->15   1750-1050A
        do i=lon0,lon1
          do k=lev0,lev1
            sigchap(k,i) = sigsrc(l)*scch4(k,i)
            fexp(k,i) = fsrc(l)*exp(-sigchap(k,i))*sfactor
!
            sum3(k,i) = sum3(k,i)+ sigsrc(l)*fexp(k,i)
            sum4(k,i) = sum4(k,i)+ hc*rlmsrcinv(l)*fexp(k,i)*
     |                  sigsrc(l)
          enddo
        enddo
      enddo
!     if (debug) write(6,"(/,'qrj after wavelength loop2: lat=',i3)")
!    | lat
!---------------------------------------------------------------
!  JCH4 Dissociation Rate modification : xnmbari(k,i) factor
!---------------------------------------------------------------
! Longitude and column domain of current process:
! Longitude-Column Loop
! -------------------------------------------------------------------
!
! Update q : From CH4 SRB-SRC Absorption Longward of 1050A:
! -- Multiply Q by efficiency factor: Titan value from Fox (1991)
! -- convert to erg/gm/sec units (internal)
! -- convert to erg/cm3/sec units (comparisons)
! Update rj: From CH4 SRB-SRC Dissociation Longward of 1050A:
!
      do i=lon0,lon1
        do k=lev0,lev1
          qsrc(k,i)=sum4(k,i)*ch4i(k,i)*rmassinv_ch4*srceff(k)
     |       *xnmbari(k,i)
          qtotal(k,i,lat) = qtotal(k,i,lat)+ sum4(k,i)*avo
     |       *ch4i(k,i)*rmassinv_ch4*srceff(k)
          qtot(k,i)=  qeuv(k,i)+qsrc(k,i)
          rj(k,i,lat) = (rj(k,i,lat)+sum3(k,i)*ch4i(k,i)*
     |       rmassinv_ch4)*xnmbari(k,i)
        enddo
      enddo
!
!    1050-1750A Neutral Heating at Titan (erg/cm3/sec): Scaled Incl. Efficiency
!     call addfsech('Q2',' ',' ',qsrc(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!    TOTAL-Neutral Heating at Titan (erg/cm3/sec): Scaled Incl. Efficiency
!     call addfsech('Q3',' ',' ',qtot(:,lon0:lon1),
!    |  lon0,lon1,nlevs,nlevs,lat)
!
!---------------------------------------------------------------
! Contributions from schumann runge bands: Earth only!
! Strobel Formula used.
!---------------------------------------------------------------
!
!    FINAL QRJ SUBROUTINE CHECKS -------------------------------
!    -- units:  Q(erg/gm/sec), Productions (#/cm3/sec)
!     call addfsech('Q'     ,' ',' ',qtotal(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('JCH4'    ,' ',' ',rj(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('JN2'   ,' ',' ',qtef(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('QN2P'  ,' ',' ',qn2p(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('QNP'   ,' ',' ',qnp(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('QCH4P' ,' ',' ',qch4p(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     call addfsech('QCH3P' ,' ',' ',qch3p(lev0:lev1,lon0:lon1,lat),
!    |  lon0,lon1,nlevs,nlevs,lat)
!     if (debug) write(6,"(/,'qrj at end of qrj sub : lat=',i3)") lat
!
!---------------------------------------------------------------
#ifdef VT
       code = 118 ; state = 'qrj' ; activity='ModelCode'
       call vtend(118,ier)
#endif
      end subroutine qrj
!---------------------------------------------------------------
      subroutine init_sflux
      use input_module,only: f107,f107a
      use init_module,only: sfeps
!
! Flux initialization once per time step, called from advance.
!
! Local:
      real :: euvflx(neuv),wave1(lmax),wave2(lmax)
      real :: flya,hlybr,fexvir,hlya,heiew,xuvfac
      integer :: iscale,n,nn
!
! Use f107,f107a from input_mod (these are either provided by
! the user, or obtained from GPI database, see input_module.F and
! gpi_module.F).
!
      hlybr = 0.
      fexvir = 0.
      hlya = 3.e+11+0.4E+10*(f107-70.)
      heiew = 0.
      xuvfac = 4.0 - (f107-68.0) / (243.0-68.0)
      if (xuvfac .lt. 1.0) xuvfac = 1.0
!
! ISCALE   =0 for Hinteregger contrast ratio method
!          =1 for Hinteregger linear interpolation
!          =2 for Tobiska EUV91 model
!          =3 for Woods & Rottman 10 Nov. 1988 measurement
!          =4 for Woods & Rottman 20 Jun. 1989 measurement
! F107     daily 10.7 cm flux (1.E-22 W m-2 Hz-1)
! F107A    81-day centered average 10.7 cm flux
! HLYBR    ratio of H Ly-b 1026A flux to solar minimum value (optional)
! FEXVIR   ratio of Fe XVI 335A flux to solar minimum value (optional)
! HLYA     H Lyman-alpha flux (photons cm-2 s-1) (optional)
! HEIEW    He I 10830A equivalent width, (milliAngstroms) (optional)
! XUVFAC   factor for scaling flux 16-250A (optional)
! WAVE1    longwave bound of spectral intervals (Angstroms)
! WAVE2    shortwave bound of intervals (= WAVE1 for indiv. lines)
! SFLUX    scaled solar flux returned by subroutine (photons cm-2 s-1)
!
      iscale = 1
      call ssflux(iscale,f107,f107a,hlybr,fexvir,hlya,
     |            heiew, xuvfac, wave1, wave2, sflux)
      call euvac(f107,f107a,euvflx)
!
! Transfer values of sflux to appropriate slots
!
      do n = l1,lmax
        feuv(n) = sflux(n)*sfeps
      enddo
      do n = 1,l1-1
        fsrc(n)   = sflux(n)*sfeps
      enddo
      do n=1,neuv
        nn = n+15
        feuv(nn) = euvflx(n)*sfeps
      enddo
!
!     write(6,"('  n   wave1   wave2       sflux      ',
!    |  'euvflx        fsrc        feuv')")
!     do n=1,lmax                   ! 1->59
!       if (n < l1) then            ! 1->15
!         write(6,"(i3,2f8.2,e12.4,12x,e12.4,e12.4)") 
!    |      n,wave1(n),wave2(n),sflux(n),fsrc(n),feuv(n)
!       elseif (n < l1+neuv-1) then ! 16->52
!         write(6,"(i3,2f8.2,e12.4,e12.4,12x,e12.4)") 
!    |      n,wave1(n),wave2(n),sflux(n),euvflx(n-l1+1),feuv(n)
!       else                        ! 53->59
!         write(6,"(i3,2f8.2,e12.4,12x,12x,e12.4)") 
!    |      n,wave1(n),wave2(n),sflux(n),feuv(n)
!       endif 
!     enddo
      end subroutine init_sflux
!-----------------------------------------------------------------
      subroutine init_qrj
!
! Called once per run, from init_module.
!
! TOTAL ABSORPTION COEFFICIENTS
!   sigeuv(1,N) = SIGA(CH4)
!   sigeuv(2,N) = SIGA(H2)
!   sigeuv(3,N) = SIGA(N2)
!
! TOTAL IONIZATION COEFFICIENTS
!   SIGMAS(1,N) = SIGI(CH4)
!   SIGMAS(2,N) = SIGI(H2)
!   SIGMAS(3,N) = SIGI(N2)
!
      RLMEUV = (/0.17250E-04, 0.16750E-04, 0.16250E-04, 0.15750E-04,
     |           0.15250E-04, 0.14750E-04, 0.14250E-04, 0.13750E-04,
     |           0.13250E-04, 0.12750E-04, 0.12250E-04, 0.12157E-04,
     |           0.11750E-04, 0.11250E-04, 0.10750E-04, 0.10250E-04,
     |           0.10319E-04, 0.10257E-04, 0.97500E-05, 0.97702E-05,
     |           0.92500E-05, 0.87500E-05, 0.82500E-05, 0.77500E-05,
     |           0.78936E-05, 0.77041E-05, 0.76515E-05, 0.72500E-05,
     |           0.70331E-05, 0.67500E-05, 0.62500E-05, 0.62973E-05,
     |           0.60976E-05, 0.57500E-05, 0.58433E-05, 0.55437E-05,
     |           0.52500E-05, 0.47500E-05, 0.46522E-05, 0.42500E-05,
     |           0.37500E-05, 0.36807E-05, 0.32500E-05, 0.30378E-05,
     |           0.30331E-05, 0.27500E-05, 0.28415E-05, 0.25630E-05,
     |           0.22500E-05, 0.17500E-05, 0.12500E-05, 0.75000E-06,
     |           0.41000E-06, 0.27500E-06, 0.19500E-06, 0.12000E-06,
     |           0.60000E-07, 0.30000E-07, 0.15000E-07/)
!
!    CH4 Total Absorption (59 bins)------------------------
!    Updated (1200-1400 A) from Chen and Wu (2004)
!    (where 121.6 : 17.86E-18 is new)
!    |             0.13, 3.29,12.83,19.07,18.32,17.86,
!    (where 121.6 : 31.0E-18 is old)
!    |             0.13, 3.29,12.83,19.07,18.32,31.00,
!
      sigeuv(1,:) = (/
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.0008,
     |             0.13, 3.29,12.83,19.07,18.32,31.00,
     |            19.20,18.00,28.40,29.11,30.12,32.70,
     |            38.19,41.15,48.00,48.97,48.33,45.92,
     |            46.47,45.72,45.46,42.93,41.07,39.28,
     |            34.99,35.28,33.18,30.70,31.05,27.92,
     |            24.64,21.45,18.77,14.70, 9.78,10.17,
     |             7.83, 6.06, 6.03, 4.36, 5.05, 3.86,
     |             2.79, 1.50, 0.59, 0.20, 0.77, 0.38,
     |             0.16, 0.00, 0.00, 0.00, 0.00/)
!
!    H2 Total Absorption (59-bins)------------------------
!
      sigeuv(2,:) = (/
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,
     |             0.05, 0.00, 0.19, 0.48, 8.25, 8.75,
     |             7.34, 8.64,10.76,11.37,10.73, 9.90,
     |             8.45, 8.46, 7.81, 6.86, 7.02, 6.17,
     |             5.36, 4.50, 3.87, 3.03, 1.90, 1.94,
     |             1.42, 1.02, 1.02, 0.73, 0.84, 0.60,
     |             0.43, 0.21, 0.08, 0.01, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00/)
!
!    N2 Total Absorption (59-bins)------------------------
!
      sigeuv(3,:) = (/
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,
     |            36.16, 0.70,16.99,46.63,15.05,30.71,
     |            19.26,26.88,35.46,30.94,26.30,29.75,
     |            23.22,23.20,23.10,22.38,23.20,24.69,
     |            24.53,21.85,21.80,21.07,17.51,18.00,
     |            13.00,11.60,11.60,10.30,10.60, 9.70,
     |             8.00, 4.40, 1.90, 0.60, 0.24, 1.16,
     |             0.48, 0.09, .015, .003, .0003/)
!
!    CH4 Total Ionization (59-bins)----------------------
!
      SIGMAS(1,:) = (/
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.48, 0.14,13.86,25.53,40.28,44.69,
     |            44.61,44.80,44.80,42.93,41.07,39.28,
     |            34.99,35.28,33.18,30.70,31.05,27.92,
     |            24.64,21.45,18.77,14.70, 9.78,10.17,
     |             7.83, 6.06, 6.03, 4.36, 5.05, 3.86,
     |             2.79, 1.50, 0.59, 0.20, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00/)
!
!    H2 to H2+ Ionization (59-bins)----------------------
!    Nearly 100% Yield.
!
      SIGMAS(2,:) = (/
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,
     |             0.00, 0.00, 0.00, 0.00, 0.09, 6.63,
     |             5.07, 7.07, 8.62, 9.76,10.73, 9.70,
     |             8.29, 8.30, 7.67, 6.74, 6.90, 6.05,
     |             5.25, 4.05, 3.78, 2.89, 1.74, 1.78,
     |             1.30, 0.90, 0.90, 0.65, 0.75, 0.55,
     |             0.40, 0.20, 0.08, 0.01, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00/)
!
!
!    N2 Total Ionization (59-bins)----------------------
!
      SIGMAS(3,:) = (/
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,
     |             0.00, 0.00, 0.00, 0.00, 0.00,16.75,
     |            10.18,18.39,23.77,23.20,23.00,25.06,
     |            23.22,23.20,23.10,22.38,23.20,24.69,
     |            24.53,21.85,21.80,21.07,17.51,18.00,
     |            13.00,11.60,11.60,10.30,10.60, 9.70,
     |             8.00, 4.40, 1.90, 0.60, 0.24, 1.16,
     |             0.48, 0.09, .015, .003, .0003/)
!
!    N2 Ionization Branching ratio yielding N+ (59-bins)
!
      BRNP = (/
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
     |             0.01, 0.04, 0.04, 0.03, 0.05, 0.05,
     |             0.15, 0.20, 0.20, 0.25, 0.32, 0.34,
     |             0.36, 0.36, 0.36, 0.36, 0.36, 0.36,
     |             0.36, 0.36, 0.36, 0.36, 0.36/)
!
!    CH4 Ionization Branching ratio yielding CH4+ (59-bins)
!
      BRCH4P = (/
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             1.00, 1.00, 1.00, 0.95, 0.64, 0.51,
     |             0.52, 0.51, 0.50, 0.50, 0.49, 0.49,
     |             0.48, 0.48, 0.47, 0.48, 0.48, 0.48,
     |             0.46, 0.46, 0.45, 0.43, 0.39, 0.39,
     |             0.38, 0.35, 0.35, 0.32, 0.33, 0.31,
     |             0.30, 0.26, 0.25, 0.25, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00/)
!
!    CH4 Ionization Branching ratio yielding CH3+ (59-bins)
!    
      BRCH3P = (/
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.05, 0.36, 0.47,
     |             0.47, 0.47, 0.48, 0.49, 0.49, 0.49,
     |             0.50, 0.50, 0.50, 0.50, 0.50, 0.50,
     |             0.49, 0.48, 0.47, 0.45, 0.42, 0.42,
     |             0.40, 0.38, 0.38, 0.35, 0.36, 0.33,
     |             0.32, 0.27, 0.26, 0.25, 0.00, 0.00,
     |             0.00, 0.00, 0.00, 0.00, 0.00/)

!
!    Titan heating efficiency (check wavelength dep.): 27%
!    Redo detailed ion-neutral chemical calculations here!
!     Ingo 3-D values -----------
!     euveff(:) = 0.50
!     srceff(:) = 0.50
!     Virginee 1-D values -----------
      euveff(:) = 0.27
      srceff(:) = 0.27
      end subroutine init_qrj
!---------------------------------------------------------------
      subroutine init_sigmas
!
! Called once per run from init_module.
!
      integer :: m,n,nn
      do m=1,3
        do n=1,lmax
          sigeuv(m,n)=sigeuv(m,n)*1.E-18
          sigmas(m,n)=sigmas(m,n)*1.E-18
        enddo
      enddo
!  ------------------------------------------------------------
!    Methane is main (only) absorber longward of 1050A for Titan.
!  ------------------------------------------------------------
      do n = 1,l1-1
        rlmsrc(n) = rlmeuv(n)
        sigsrc(n) = sigeuv(1,n)
      enddo
!     Overlay EUVAC Cross-sections on sigeuv(:,:) and sigmas(:,:)
!     N2-Stuff on init-euvac only!
      do n = 1,neuv
        nn = n+15    ! 16:52
!     Total N2 absorption
        sigeuv(3,nn) = sigan2(n)
!     Total N2 ionization
        sigmas(3,nn) = sigin2(n)
!     Branching for N2 ionization yielding N+
        brnp(nn) = brn2np(n)
      enddo
      end subroutine init_sigmas
!-----------------------------------------------------------------------
      subroutine init_euvac
!
! Called once per run from init_module.
!
! lambdas:
      wleuv1 = 
     +    (/1000.00, 1031.91, 1025.72,  950.00,  977.02,  900.00,
     +       850.00,  800.00,  750.00,  789.36,  770.41,  765.15,
     +       700.00,  703.36,  650.00,  600.00,  629.73,  609.76,
     +       550.00,  584.33,  554.31,  500.00,  450.00,  465.22,
     +       400.00,  350.00,  368.07,  300.00,  303.78,  303.31,
     +       250.00,  284.15,  256.30,  200.00,  150.00,  100.00,
     +        50.00/)
      wleuv2 =
     +    (/1050.00, 1031.91, 1025.72, 1000.00,  977.02,  950.00,
     +       900.00,  850.00,  800.00,  789.36,  770.41,  765.15,
     +       750.00,  703.36,  700.00,  650.00,  629.73,  609.76,
     +       600.00,  584.33,  554.31,  550.00,  500.00,  465.22,
     +       450.00,  400.00,  368.07,  350.00,  303.78,  303.31,
     +       300.00,  284.15,  256.30,  250.00,  200.00,  150.00,
     +       100.00/)
c
c sigan2 : Updated Total N2-absorption cross sections (50-1050A)
      sigan2 =
     +    (/0.00e+00, 0.00e+00, 0.00e+00, 5.10e-17, 2.24e-18, 9.68e-18,
     +      2.02e-17, 1.70e-17, 3.36e-17, 1.65e-17, 1.42e-17, 1.20e-16,
     +      2.47e-17, 2.65e-17, 3.18e-17, 2.33e-17, 2.34e-17, 2.28e-17,
     +      2.28e-17, 2.24e-17, 2.41e-17, 2.45e-17, 2.35e-17, 2.32e-17,
     +      2.17e-17, 1.64e-17, 1.69e-17, 1.39e-17, 1.17e-17, 1.17e-17,
     +      1.05e-17, 1.09e-17, 1.02e-17, 8.39e-18, 4.96e-18, 2.26e-18,
     +      7.20e-19/)
c
c sigin2 : Updated Total N2-ionization cross sections (50-1050A)
      sigin2 =
     +    (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00,
     +      0.00e+00, 0.00e+00, 1.43e-17, 8.86e-18, 8.50e-18, 6.58e-17,
     +      1.51e-17, 2.55e-17, 2.92e-17, 2.33e-17, 2.34e-17, 2.28e-17,
     +      2.28e-17, 2.24e-17, 2.41e-17, 2.45e-17, 2.35e-17, 2.32e-17,
     +      2.17e-17, 1.64e-17, 1.69e-17, 1.39e-17, 1.17e-17, 1.17e-17,
     +      1.05e-17, 1.09e-17, 1.02e-17, 8.39e-18, 4.96e-18, 2.26e-18,
     +      7.20e-19/)
c
c N2 -> N+ : N2-ionization branching ratios yielding N+ (50-1050A)
      brn2np =
     +    (/0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00,
     +      0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00,
     +      0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00, 0.00e+00,
     +      0.00e+00, 0.00e+00, 0.00e+00, 1.30e-03, 3.00e-02, 4.60e-02,
     +      4.50e-02, 1.05e-01, 9.00e-02, 1.63e-01, 2.13e-01, 2.13e-01,
     +      3.00e-01, 2.57e-01, 3.35e-01, 3.77e-01, 3.64e-01, 3.46e-01,
     +      3.85e-01/)
      end subroutine init_euvac
!-----------------------------------------------------------------------
      subroutine euvac(F107,F107A,EUVFLX)
!
! This EUV flux model uses the F74113 solar reference spectrum and
! ratios determined from Hinteregger's SERF1 model. It uses the daily
! F10.7 flux (F107) and the 81 day mean (F107A) as a proxy for scaling
! The fluxes are returned in EUVFLX and correspond to the 37 wavelength
! bins of Torr et al. [1979] Geophys. Res. Lett. p771.
! See Richards et al. [1994] J. Geophys. Res. p8981 for details.
!
! Called from init_sflux, once per timestep.
!
! F107   = input daily 10.7 cm flux index. (e.g. 74)
! F107A  = input 81 day average of daily F10.7 centered on current day
! EUVFLX = output array for EUV flux in units of photons/cm2/sec.
!
! Reduced argument list for Titan:
!     real :: wleuv1,wleuv2,sigan2,sigin2,brn2np
!
! Args:
      real,intent(in) :: f107,f107a
      real,intent(out) :: euvflx(neuv)
!
! Local:
      integer :: i
      real :: AFAC(neuv),F74113(neuv),flxfac
      real :: afac_tmp(neuv),f74113_tmp(neuv)
!
! F74113 reference spectrum (doubled below 150-250 A, tripled <150)
! Will be multiplied by 1.0E9 later
!
      F74113=
     | (/1.200,0.450,4.800,3.100,0.460,0.210,1.679,0.800,6.900,
     |   0.965,0.650,0.314,0.383,0.290,0.285,0.452,0.720,1.270,
     |   0.357,0.530,1.590,0.342,0.230,0.360,0.141,0.170,0.260,
     |   0.702,0.758,1.625,3.537,3.000,4.400,1.475,3.500,2.100,
     |   2.467/)
!
! Scaling factors(Ai) for the EUV flux
      AFAC=
     | (/1.0017E-02,7.1250E-03,1.3375E-02,1.9450E-02,2.7750E-03,
     |   1.3768E-01,2.6467E-02,2.5000E-02,3.3333E-03,2.2450E-02,
     |   6.5917E-03,3.6542E-02,7.4083E-03,7.4917E-03,2.0225E-02,
     |   8.7583E-03,3.2667E-03,5.1583E-03,3.6583E-03,1.6175E-02,
     |   3.3250E-03,1.1800E-02,4.2667E-03,3.0417E-03,4.7500E-03,
     |   3.8500E-03,1.2808E-02,3.2750E-03,4.7667E-03,4.8167E-03,
     |   5.6750E-03,4.9833E-03,3.9417E-03,4.4167E-03,5.1833E-03,
     |   5.2833E-03,4.3750E-03/)
!
! 6/20/03 bf: Reverse order of F74113 and AFAC to correspond to descending
! wavelength bins (as in wleuv1,2 and WAVEL,WAVES):
      do i=1,neuv
        afac_tmp(i) = afac(neuv-i+1)  
        f74113_tmp(i) = f74113(neuv-i+1)  
      enddo
      afac = afac_tmp      ! whole array
      f74113 = f74113_tmp
!
! Loop through the wavelengths calculating the scaling factors and
!   the resulting solar flux.
! The scaling factors are restricted to be greater than 0.8
!
      DO 50 I=1,neuv
        FLXFAC=(1.0 + AFAC(I) * (0.5*(F107+F107A) - 80.0))
        IF(FLXFAC.LT.0.8) FLXFAC=0.8
        EUVFLX(I)=F74113(I) * FLXFAC * 1.0E9
 50   CONTINUE
      end subroutine euvac
!-----------------------------------------------------------------------
      subroutine ssflux (iscale, f107, f107a, hlybr, fexvir, hlya,
     |                   heiew, xuvfac, wave1, wave2, sflux)
      integer,parameter :: lmax=59
!
! Args:
      integer,intent(in) :: iscale
      real,intent(in) :: f107,f107a,hlybr,fexvir,hlya,heiew,xuvfac
      real,intent(out) :: wave1(lmax),wave2(lmax),sflux(lmax) 
!
! Local:
      real ::
     |          wavel(lmax), waves(lmax), rflux(lmax), xflux(lmax),
     |          scale1(lmax), scale2(lmax),
     |          tchr0(lmax), tchr1(lmax), tchr2(lmax),
     |          tcor0(lmax), tcor1(lmax), tcor2(lmax),
     |          war1(lmax), war2(lmax),
     |          b1(3), b2(3)
      real :: frat,r1,r2,hlymod,heimod,xuvf
      integer :: l
!
! regression coefficients which reduce to solar min. spectrum:
      b1 = (/1.0, 0.0138, 0.005/) 
      b2 = (/1.0, 0.59425, 0.3811/)
!
! 'best fit' regression coefficients, commented out, for reference:
!     b1 = (/1.31, 0.01106, 0.00492/, B2/-6.618, 0.66159, 0.38319/)
!
      WAVEL =   (/1750.00, 1700.00, 1650.00, 1600.00, 1550.00, 1500.00,
     |            1450.00, 1400.00, 1350.00, 1300.00, 1250.00, 1215.67,
     |            1200.00, 1150.00, 1100.00, 1050.00, 1031.91, 1025.72,
     |            1000.00,  977.02,  950.00,  900.00,  850.00,  800.00,
     |             789.36,  770.41,  765.15,  750.00,  703.31,  700.00,
     |             650.00,  629.73,  609.76,  600.00,  584.33,  554.37,
     |             550.00,  500.00,  465.22,  450.00,  400.00,  368.07,
     |             350.00,  303.78,  303.31,  300.00,  284.15,  256.30,
     |             250.00,  200.00,  150.00,  100.00,   50.00,   32.00,
     |              23.00,   16.00,    8.00,    4.00,    2.00/)
      WAVES =   (/1700.00, 1650.00, 1600.00, 1550.00, 1500.00, 1450.00,
     |            1400.00, 1350.00, 1300.00, 1250.00, 1200.00, 1215.67,
     |            1150.00, 1100.00, 1050.00, 1000.00, 1031.91, 1025.72,
     |             950.00,  977.02,  900.00,  850.00,  800.00,  750.00,
     |             789.36,  770.41,  765.15,  700.00,  703.31,  650.00,
     |             600.00,  629.73,  609.76,  550.00,  584.33,  554.37,
     |             500.00,  450.00,  465.22,  400.00,  350.00,  368.07,
     |             300.00,  303.78,  303.31,  250.00,  284.15,  256.30,
     |             200.00,  150.00,  100.00,   50.00,   32.00,   23.00,
     |              16.00,    8.00,    4.00,    2.00,    1.00/)
      RFLUX =    (/322.00,  168.00,   95.00,   62.00,   44.00,   25.00,
     |              16.90,   11.80,   19.50,    4.10,   11.10,  249.00,
     |               2.78,    0.70,    3.07,    3.64,    3.18,    4.38,
     |               1.78,    5.96,    4.22,    4.43,    1.93,    0.87,
     |               0.79,    0.24,    0.20,    0.17,    0.39,    0.22,
     |               0.17,    1.50,    0.45,    0.48,    1.58,    0.80,
     |               0.51,    0.31,    0.18,    0.39,    0.21,    0.74,
     |               0.87,    6.00,    0.24,    0.84,    0.10,    0.27,
     |               0.92,    1.84,    0.13,    0.38,  0.0215,  0.0067,
     |              1.E-3,   2.E-3,   1.E-5,   5.E-8,   1.E-10/)
      XFLUX =    (/354.00,  191.00,  110.00,   76.00,   55.00,   28.00,
     |              19.60,   14.30,   25.30,    5.00,   17.20,  401.00,
     |               6.26,    1.51,    6.11,    8.66,    9.04,   13.12,
     |               4.42,   13.18,   12.03,   13.29,    5.01,    2.18,
     |               1.59,    0.67,    0.43,    0.43,    0.72,    0.46,
     |               0.48,    3.02,    1.46,    1.02,    4.86,    1.59,
     |               1.57,    1.67,    0.36,    0.99,    2.20,    1.39,
     |               5.63,   11.28,    2.50,    4.14,    3.16,    0.59,
     |               3.70,    4.85,    0.34,    1.15,    0.18,    0.08,
     |             2.5E-2,   5.E-2,   8.E-4,   3.E-5,   5.E-7/)
      SCALE1 = (/     0.0,     0.0,     0.0,     0.0,     0.0,     0.0,
     |                0.0,     0.0,     0.0,     0.0,     0.0,     0.0,
     |            1692.09,  405.95, 1516.20, 2731.70, 3314.57, 4375.00,
     |            1316.91, 3621.91, 3908.56, 4432.54, 1541.21,  531.73,
     |             364.83,    0.00,  116.00,  129.41,  162.48,   94.07,
     |              41.29,  709.50,    0.00,  268.47, 1561.05,  367.64,
     |             290.06,  184.36,    0.00,   86.15,    7.50,    0.00,
     |               0.00, 2220.00,    0.00,   61.00,    0.00,   86.95,
     |             206.00,  135.89,   60.35,  157.12,    7.06,    0.75,
     |               0.00,    0.00,    0.00,    0.00,    0.00/)
      SCALE2 =  (/   0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    5.34,    0.00,    0.00,    0.00,    0.54,
     |               3.30,    0.00,   12.60,    0.00,    0.00,    0.00,
     |               5.34,   11.63,    2.28,    5.56,   24.93,    8.16,
     |              60.69,    0.00,   28.20,   45.90,   40.80,    1.27,
     |              35.47,   42.80,    1.12,    6.19,    1.26,    0.69,
     |               0.23,    0.46,  7.6E-3,  2.9E-4,  4.8E-6/)
      TCHR0 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,-4.290E+00,-5.709E+00,-8.493E+00,
     |-1.161E+00,-3.429E+00,-5.464E+00,-6.502E+00,-1.912E+00,-4.034E-01,
     |-1.448E-01, 0.000E+00,-9.702E-02,-6.591E-02,-2.338E-02,-1.273E-01,
     |-2.406E-01,-3.351E-01, 0.000E+00,-1.465E+00,-2.405E+00,-7.975E-02,
     |-4.197E-01,-1.971E-01, 0.000E+00,-5.895E-02,-5.815E-03, 0.000E+00,
     | 0.000E+00, 2.138E-01, 0.000E+00,-7.713E-02, 0.000E+00,-3.035E-02,
     |-2.039E-01,-1.749E-01,-1.041E-01,-2.638E-01,-1.094E-02, 0.000E+00,
     |       0.0,       0.0,       0.0,       0.0,       0.0/)
      TCHR1 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,-3.023E-13,-3.745E-13,-5.385E-13,
     |-1.211E-13,-3.868E-13,-3.646E-13,-4.125E-13,-1.527E-13,-4.753E-14,
     |-3.411E-14, 0.000E+00,-1.190E-14,-1.034E-14,-1.343E-14,-1.539E-14,
     |-5.174E-14,-6.934E-14, 0.000E+00,-1.215E-13,-1.537E-13,-2.024E-14,
     |-4.596E-14,-1.562E-14, 0.000E+00,-1.221E-14,-1.123E-15, 0.000E+00,
     | 0.000E+00,-2.263E-13, 0.000E+00,-1.508E-14, 0.000E+00,-1.744E-14,
     |-2.100E-14,-1.805E-14,-8.224E-15,-1.919E-14,-7.944E-16, 0.000E+00,
     |       0.0,       0.0,       0.0,       0.0,       0.0/)
      TCHR2 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0, 3.275E-11, 4.057E-11, 6.160E-11,
     | 1.312E-11, 4.189E-11, 4.167E-11, 4.716E-11, 1.654E-11, 5.150E-12,
     | 3.901E-12, 0.000E+00, 1.289E-12, 1.120E-12, 1.455E-12, 1.667E-12,
     | 5.604E-12, 7.931E-12, 0.000E+00, 1.317E-11, 1.757E-11, 2.194E-12,
     | 4.978E-12, 1.693E-12, 0.000E+00, 1.324E-12, 1.285E-13, 0.000E+00,
     | 0.000E+00, 2.586E-11, 0.000E+00, 1.724E-12, 0.000E+00, 1.889E-12,
     | 2.400E-12, 2.063E-12, 8.911E-13, 2.193E-12, 9.090E-14, 0.000E+00,
     |       0.0,       0.0,       0.0,       0.0,       0.0/)
      TCOR0 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0, 0.000E+00, 0.000E+00, 0.000E+00,
     | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00,-6.060E-02,
     | 0.000E+00,-3.399E-02, 0.000E+00, 0.000E+00, 0.000E+00, 4.866E-02,
     |-1.762E-01, 0.000E+00,-2.412E-01, 0.000E+00, 0.000E+00, 0.000E+00,
     |-4.743E-01,-9.713E-01, 5.891E-02,-1.263E-01,-1.246E+00, 2.870E-01,
     |-4.659E+00, 0.000E+00,-1.058E+00,-3.821E+00,-1.874E+00, 0.000E+00,
     |-1.896E+00,-8.505E-01,-2.101E-04,-2.012E-01,-6.097E-02,-2.925E-02,
     |-4.875E-03,       0.0,       0.0,       0.0,       0.0/)
      TCOR1 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0, 0.000E+00, 0.000E+00, 0.000E+00,
     | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 2.877E-03,
     | 0.000E+00, 1.760E-03, 0.000E+00, 0.000E+00, 0.000E+00, 3.313E-04,
     | 3.643E-03, 0.000E+00, 5.225E-03, 0.000E+00, 0.000E+00, 0.000E+00,
     | 4.085E-03, 1.088E-02, 8.447E-04, 3.237E-03, 1.907E-02, 2.796E-03,
     | 4.460E-02, 0.000E+00, 1.007E-02, 3.481E-02, 1.604E-02, 0.000E+00,
     | 2.029E-02, 2.160E-02, 6.342E-04, 3.594E-03, 5.503E-04, 2.687E-04,
     | 4.479E-05,       0.0,       0.0,       0.0,       0.0/)
      TCOR2 = (/
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0,       0.0,       0.0,       0.0,
     |       0.0,       0.0,       0.0, 0.000E+00, 0.000E+00, 0.000E+00,
     | 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.846E-03,
     | 0.000E+00, 1.127E-03, 0.000E+00, 0.000E+00, 0.000E+00, 1.891E-04,
     | 2.326E-03, 0.000E+00, 2.801E-03, 0.000E+00, 0.000E+00, 0.000E+00,
     | 2.446E-03, 7.121E-03, 5.204E-04, 1.983E-03, 1.204E-02, 1.721E-03,
     | 2.911E-02, 0.000E+00, 7.177E-03, 2.272E-02, 9.436E-03, 0.000E+00,
     | 1.316E-02, 1.398E-02, 4.098E-04, 2.328E-03, 3.574E-04, 1.745E-04,
     | 2.909E-05,       0.0,       0.0,       0.0,       0.0/)
      WAR1 =  (/     0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    3.80,    6.25,    4.93,    6.06,
     |               2.70,    7.07,    8.62,    9.60,    4.54,    2.37,
     |               0.82,    0.33,    0.24,    0.67,    0.28,    0.55,
     |               1.56,    1.11,    0.77,    1.32,    1.71,    0.44,
     |               1.11,    0.95,    0.39,    0.81,    2.00,    1.49,
     |               6.81,    5.07,    1.63,    5.62,    2.08,    0.59,
     |               3.89,    5.19,    0.35,    1.18,    0.099,   0.04,
     |               0.007,   0.00,    0.00,    0.00,    0.00/)
      WAR2 =  (/     0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,    0.00,    0.00,    0.00,    0.00,
     |               0.00,    0.00,   20.80,   17.90,    9.30,   14.30,
     |               6.90,   12.00,   15.60,   18.60,   10.10,    4.30,
     |              12.40,    8.00,    3.60,    1.80,    0.50,    1.40,
     |               3.90,    2.60,    1.60,    3.40,    4.10,    0.70,
     |               4.30,    4.30,    3.80,    2.60,    6.08,    1.35,
     |              12.60,    9.78,    2.96,   10.20,    4.11,    6.68,
     |               6.62,    8.07,    0.47,    1.73,    0.17,    0.075,
     |               0.012,   0.00,    0.00,    0.00,    0.00/)
C
C Linear Interpolation between SC#21REFW and F79050:
C
      FRAT = (F107-68.) / (243.-68.)
      DO 200 L=1,LMAX
        SFLUX(L) = RFLUX(L) + (XFLUX(L)-RFLUX(L)) * FRAT
  200 CONTINUE
C
C Hinteregger contrast ratio method:
C
      IF (ISCALE .EQ. 0) THEN
        IF (HLYBR .GT. 0.001) THEN
          R1 = HLYBR
        ELSE
          R1 =  B1(1) + B1(2)*(F107A-71.5) + B1(3)*(F107-F107A+3.9)
        ENDIF
        IF (FEXVIR .GT. 0.001) THEN
          R2 = FEXVIR
        ELSE
          R2 =  B2(1) + B2(2)*(F107A-71.5) + B2(3)*(F107-F107A+3.9)
        ENDIF
        DO 100 L=13,LMAX
          SFLUX(L) = (RFLUX(L) + ((R1-1.)*SCALE1(L)
     |                         +  (R2-1.)*SCALE2(L)) / 1000.)
  100   CONTINUE
      ENDIF
C
C Tobiska EUV91 Method:
C
      IF (ISCALE .EQ. 2) THEN
        IF (HLYA .GT. 0.001) THEN
          HLYMOD = HLYA
        ELSE
          IF (HEIEW .GT. 0.001) THEN
             HLYMOD = HEIEW * 3.77847E9 + 8.40317E10
          ELSE
             HLYMOD = 8.70E8 * F107 + 1.90E11
             HLYMOD = 8.70E8 * F107 + 1.90E11
          ENDIF
        ENDIF
        IF (HEIEW .GT. 0.001) THEN
          HEIMOD = HEIEW * 3.77847E9 + 8.40317E10
        ELSE
          HEIMOD = HLYMOD
        ENDIF
        DO 500 L=16,55
          SFLUX(L) = TCHR0(L) + TCHR1(L)*HLYMOD + TCHR2(L)*HEIMOD
     |             + TCOR0(L) + TCOR1(L)*F107 + TCOR2(L)*F107A
 500    CONTINUE
      ENDIF
C
C Woods and Rottman (10 Nov. 1988) spectrum:
C
      IF (ISCALE .EQ. 3) THEN
        DO 550 L=15,55
        SFLUX(L) = WAR1(L)
  550   CONTINUE
      ENDIF
C
C Woods and Rottman (20 June 1989) spectrum:
C
      IF (ISCALE .EQ. 4) THEN
        DO 560 L=15,55
        SFLUX(L) = WAR2(L)
  560   CONTINUE
      ENDIF
C
C Substitute in H Lyman-alpha and XUVFAC if provided:
C -- Update of Lyman-alpha fluxes over solar cycle
C
      IF (HLYA .GT. 0.001) SFLUX(12) = HLYA / 1.E9
      IF (XUVFAC .GT. 0.001) THEN
        XUVF = XUVFAC
      ELSE
        XUVF = 1.0
      ENDIF
C
C Convert from gigaphotons to photons, etc.:
C
      DO 600 L=1,LMAX
        WAVE1(L) = WAVEL(L)
        WAVE2(L) = WAVES(L)
        IF (SFLUX(L) .LT. 0.0) SFLUX(L) = 0.0
        IF (WAVEL(L).LT.251.0 .AND. WAVES(L).GT.15.0)
     |     SFLUX(L)=SFLUX(L)*XUVF
        SFLUX(L) = SFLUX(L) * 1.E9
  600 CONTINUE
C
      end subroutine ssflux
!-----------------------------------------------------------------------
      subroutine alloc_q(lon0,lon1,lat0,lat1)
!
! Args:
      integer,intent(in) :: lon0,lon1,lat0,lat1
!
! Local:
      integer :: istat
!
! Allocate (7-fields): rj, qtef, qtotal, qch4p, qch3p, qn2p, qnp 
!
      allocate(rj(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' rj: stat=',i3)") istat
      allocate(qtef(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qtef: stat=',i3)") istat
      allocate(qtotal(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qtotal: stat=',i3)") istat
      allocate(qch4p(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qch4p: stat=',i3)") istat
      allocate(qch3p(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qch3p: stat=',i3)") istat
      allocate(qn2p(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qn2p: stat=',i3)") istat
      allocate(qnp(nlevp1,lon0:lon1,lat0:lat1),stat=istat)
      if (istat /= 0) write(6,"('>>> alloc_q: error allocating',
     |  ' qnp: stat=',i3)") istat
      end subroutine alloc_q
!-----------------------------------------------------------------------
      end module qrj_module
