      module solred_module
      implicit none
      private

      integer,parameter :: lmax=59

      real,parameter :: dtor=atan(1.)/45.
      real,parameter :: rtod=1./dtor
      real,parameter :: re=6371.e+5 
      real,parameter :: pi=180./rtod
      real,parameter :: boltz=1.38e-16
      real,parameter :: a6300=5.15e-3
      real,parameter :: a1d=6.81e-3
      real,parameter :: quench=2.3e-11
      real,parameter :: erg=1.602e-12
      real,parameter :: c1=23.5*dtor            

      real,save :: sigabs(9,lmax),sigion(13,lmax)

      public :: solred
      contains
!-----------------------------------------------------------------------
      subroutine solred(slt,zpht,xno,xno2,xnn2,tn,f107,f107a,
     |  day,glat,glon,kmx,sr63,qn2p,qop,glyb)
!
      implicit none
!
! Args:
      integer,intent(in) :: kmx
      real,intent(in) :: slt,f107,f107a,day,glat,glon
      real,dimension(kmx),intent(in) :: zpht,xno,xno2,xnn2,tn
      real,dimension(kmx),intent(out) :: sr63,qn2p,qop,glyb
!
! Local:
      integer :: k,l,iscale
      real,dimension(lmax) :: brn2,bro2,sig
      real,dimension(kmx) :: qtin,qsr,qo2p,txno,txno2,txnn2,gz,rho,sht
      real,dimension(kmx) :: cno,cno2,cnn2
      real,dimension(lmax) :: wave1,wave2,sflux
      real :: xmb,dle,dles,dlec,glatr,glonr,seczi,zi,alt,xp,y,tt,yerf,
     |  rg,zg,shg,bfd,chapo,chan2,chaps,taur,taup,taur1,taur2,taur3,
     |  ctaur,ctaur1,ctaur2,ctaur3,ctaup,rsp,rsq,cspi,cspj,taubt,chao2,
     |  tau,ctau,thng1,thng2,thng3,thng
      real :: xmo2,xmn2,xmo,tmm,cngo2,cngn2,cngo

      data brn2/
     |              36 * 0.00,
     |             0.01, 0.04, 0.04, 0.03, 0.05, 0.05,
     |             0.15, 0.20, 0.20, 0.25, 0.32, 0.34,
     |             11 * 0.36/
      data bro2/
     |              30 * 0.00,
     |             .025, .036, .045, .120, .155, .189,
     |             .230, 0.20, 0.20, 0.20, 0.23, 0.25,
     |             0.29, 16 * 0.33/

      sr63=0. ; qn2p=0. ; qop=0. ; glyb=0.
      iscale = 0
!
! Scalars wave1,wave2,sflux are output by ssflux:
!
      call ssflux(iscale,f107,f107a,wave1,wave2,sflux)

      do k=1,kmx              
        sr63(k) = 0.
        glyb(k) = 1.e-20 ! 5th component of e5577
        qop(k)=0.                                                  
        qo2p(k)=0.
        qn2p(k)=0.                   
        qtin(k)=1.E-20
        qsr(k)=1.E-20       
        gz(k)=980.665*(1.-3.14466e-4*zpht(k))
        rho(k)=(32.*xno2(k)+28.*xnn2(k)+16.*xno(k))*1.66e-24
        xmb=(28.*xnn2(k)+32.*xno2(k)+16.*xno(k))/(xnn2(k)+
     |    xno2(k)+xno(k))
        sht(k)=boltz*tn(k)/(xmb*gz(k)*1.66e-24)
      enddo
!
! cno,cno2,cnn2 (kmx) are output by column:
!
      call column(zpht,xno,xno2,tn,xnn2,sht,gz,cno,cno2,cnn2,kmx)
!
!     write(6,"('solred after column: cno=',/,(6e12.4))") cno
!     write(6,"('solred after column: cno2=',/,(6e12.4))") cno2
!     write(6,"('solred after column: cnn2=',/,(6e12.4))") cnn2

      dle=atan(tan(c1)*sin(2.*pi*(day-80.)/365.))
      dles=sin(dle)
      dlec=cos(dle)
      glatr=glat*dtor
      glonr=glon*dtor
      seczi=1./(dles*sin(glatr)+cos(glatr)*cos(pi*(slt-12.)/12.)*dlec)
      zi=acos(1./seczi)*rtod
!
! 7/11/95: Return if solar zenith angle > 110 deg:
!
      if (zi > 110.) return

      do k=1,kmx
        alt=zpht(k)
        xp=(re+alt*1.e+5)/sht(k)
        y=0.5*xp*(cos(zi*dtor))**2
        tt=sqrt(y)
        if(tt > 8.) then
          yerf=0.56498823/(0.06651874+tt)
        else
          yerf=(1.0606963+0.55643831*tt)/(1.0619896+1.7245609*tt+tt*tt)
        endif

        if (zi > 90.) then
          rg=(re+alt*1.e+5)*sin(zi*dtor)
          zg=rg-re
          if (zg <= 10.e+5) zg=10.e+5
!
! Scalars xmo2,xmn2,xmo,tmm,cngo2,cngn2,cngo are output by tancomp:
!
          call tancomp(glat,glon,zg,xmo2,xmn2,xmo,tmm,cngo2,cngn2,cngo,
     |      zpht,xno,xno2,tn,xnn2,sht,cno,cno2,cnn2,kmx)

          shg=boltz*tmm/(28.*1.66e-24*gz(k))
          bfd=0.5*pi*rg/shg
          chapo=sqrt(bfd*16.)
          chao2=sqrt(bfd*32.)
          chan2=sqrt(bfd*28.)
          txno2(k)=chao2*(2.*cngo2-xno2(k)*shg*yerf/32.)
          txnn2(k)=chan2*(2.*cngn2-xnn2(k)*shg*yerf/28.)
          txno(k)=chapo*(2.*cngo-xno(k)*shg*yerf/16.)
        else ! zi <= 90
          chaps=sqrt(0.5*pi*xp)*yerf
          txno2(k)=cno2(k)*chaps
          txnn2(k)=cnn2(k)*chaps
          txno(k)=cno(k)*chaps
        endif ! zi

        taur=sigabs(1,49)*txno(k)+sigabs(3,49)*txnn2(k)+
     |    sigabs(2,49)*txno2(k)
        taup=sigabs(1,20)*txno(k)+sigabs(3,20)*txnn2(k)+
     |    sigabs(2,20)*txno2(k)
        taur1=1.3*taur
        taur2=2.0*taur
        taur3=2.5*taur
        if (taur  > 9.) taur=9.
        if (taup  > 9.) taup=9.
        if (taur1 > 9.) taur1=9.
        if (taur2 > 9.) taur2=9.
        if (taur3 > 9.) taur3=9.
        ctaur=exp(-taur)
        ctaur1=exp(-taur1)
        ctaur2=exp(-taur2)
        ctaur3=exp(-taur3)
        ctaup=exp(-taup)
        rsp=2.4*ctaur/(ctaur+2.*(ctaur1+ctaur2+ctaur3))
        rsq=1.5*ctaur/(ctaur+2.*(ctaur1+ctaur2+ctaur3)+taup/taur*ctaup)
        cspi=1.0+rsp
        cspj=1.0+rsq
c
c glyb(kmx) = photo dissociation of o2 by solar lyman-beta, returned
c             for use as 5th source of greenline e5577 emission:
c
        taubt = sigabs(1,17)*txno(k)+sigabs(2,17)*txno2(k)+
     |    sigabs(3,17)*txnn2(k)
        glyb(k) = 0.1*sflux(17)*xno2(k)*sigabs(2,17)*exp(-taubt)
c
        do l=1,lmax
          tau=sigabs(1,l)*txno(k)+sigabs(2,l)*txno2(k)+sigabs(3,l)*
     |      txnn2(k)
          ctau=exp(-tau)
          thng1=xno2(k)*sigion(2,l)*ctau*sflux(l)
          thng2=xnn2(k)*sigion(3,l)*ctau*sflux(l)
          thng3=xno(k)*sigion(1,l)*ctau*sflux(l)
          qo2p(k)=qo2p(k)+thng1*(1.-bro2(l)+rsq)
          qn2p(k)=qn2p(k)+thng2*(1.-brn2(l)+rsp)
          qop(k)=qop(k)+thng3*(1.+rsp)+thng1*bro2(l)
          if (l <= 11) then
            tau=sigabs(2,l)*txno2(k)
            thng=sflux(l)*sigabs(2,l)*exp(-tau)
            sr63(k)=sr63(k)+xno2(k)*thng*a6300/(a1d*(1.+quench*xnn2(k)/
     |        a1d))
          endif
        enddo ! l=1,lmax
        qtin(k)=qo2p(k)+qn2p(k)+qop(k)
      enddo ! do k=1,kmx

!     write(6,"('solred: glat=',f8.2,' glon=',f8.2,' kmx=',i4,' sr63=',
!    |  /,(6e12.4))") glat,glon,kmx,sr63
!     call shutdown('debug sr63')

      end subroutine solred
!-----------------------------------------------------------------------
      subroutine ssflux(iscale,f107,f107a,wave1,wave2,sflux)
      implicit none
!
! Definitions (S.C. Solomon, 12/88):
! ISCALE   =0 for contrast ratio method, =1 for linear interpolation
! F107     daily 10.7 cm flux
! F107A    81-day centered average 10.7 cm flux
! WAVE1    longwave bound of spectral intervals
! WAVE2    shortwave bound of spectral intervals (= WAVE1 for indiv. lin
! SFLUX    scaled solar flux returned by subroutine
! LMAX     dimension of WAVE1, WAVE2, and SFLUX arrays, must be <= LM
! WAVEL    = WAVE1
! WAVES    = WAVE2
! RFLUX    low solar activity reference flux
! XFLUX    high solar activity flux
! SCALE1   scaling factors for H LyB-keyed chromospheric emissions
! SCALE2   scaling factors for FeXVI-keyed coronal emissions
! LM       dimension of above arrays, currently = 59
! SRA      'A' value for S-R continuum scaling formula
! SRB      'B' value for S-R continuum scaling formula
! B1       fit coefficients for H LyB
! B2       fit coefficients for FeXVI
! R1       enhancement ratio for H LyB
! R2       enhancement ratio for FeXVI
! SFNORM   normalization factor for scaling flux shortwards of 250A
!
! Args:
      integer,intent(in) :: iscale
      real,intent(in) :: f107,f107a
      real,dimension(lmax),intent(out) :: wave1,wave2,sflux
!
! Local:
      integer :: l,m
      real :: b1(3),b2(3), sra(8),srb(8), r1,r2,frat,sfnorm
      real,dimension(lmax) :: wavel,waves,rflux,xflux,scale1,scale2
      real,dimension(lmax) :: sigao,sigao2,sigan2,sigio,sigin2,
     |  sigio2,sigaco,sigico,sigaco2,sigico2,sigio4s,sigio2d,sigio2p,
     |  sigio4p,sigio2q,sigihe,sigin,sigih     
      logical :: first=.true.

      data b1/1.0, 0.0138, 0.005/, b2/1.0, 0.59425, 0.3811/

      DATA 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/          
      DATA 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/
      DATA RFLUX/  370.45,  203.69,   96.00,   69.71,   50.70,   26.67, 
     >              17.21,    8.26,   12.86,    4.10,    5.20,  333.80,
     >               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,
     >             0.0009,  0.0003,   1.E-6,   3.E-9,   1.E-11/
      DATA XFLUX/  464.20,  241.50,  131.50,  101.90,   81.32,   48.71,
     >              37.16,   21.14,   30.70,   11.20,   12.00,  438.80,
     >               6.50,    1.60,    6.40,    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,
     >              0.025,    0.03,   1.E-3,   3.E-5,   1.E-6/
      DATA SCALE1/35347.5, 33095.6, 18040.6, 13733.0, 12564.2, 7121.38,
     >            6608.74, 5779.89, 8009.80, 3186.34, 3033.78,  47555.,
     >            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/
      DATA 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.30,    0.01,   3.E-4,   1.E-5/
      DATA SRA/     0.536,   0.216,   0.203,   0.184,   0.175,   0.126,
     >              0.114,   0.073/
      DATA SRB/     334.0,   189.0,    82.2,    57.2,    38.8,    18.1,
     >               9.46,    3.30/
      DATA SIGAO /  18 * 0.00,                                          
     >             0.00, 0.00, 2.12, 4.18, 4.38, 4.23,                  
     >             4.28, 4.18, 4.18, 8.00,11.35,10.04,                  
     >            12.21,12.22,12.23,11.90,12.17,12.13,                  
     >            11.91,11.64,11.25,11.21, 9.64, 9.95,                  
     >             8.67, 7.70, 7.68, 6.61, 7.13, 6.05,                  
     >             5.30, 2.90, 1.60, 0.59, 0.16, 0.05,                  
     >             0.51, 0.07, .012, .002, .0002/                       
!                                                                       
      DATA SIGAO2/ 0.50, 1.50, 3.40, 6.00,10.00,13.00,                  
     >            15.00,12.00, 2.20, 0.30, 3.00, 0.01,                  
! 1/30/04: changed siga02(17) from 1.0 to 2.2, as per roble,
!          for lyman-beta source of E5577 (ie5577(5) > 0)
     >             0.30, 0.10, 1.00, 1.10, 2.20, 1.60,                  
     >            16.53, 4.00,15.54, 9.85,20.87,27.09,                  
     >            26.66,25.18,21.96,29.05,25.00,26.27,                  
     >            26.02,25.80,26.10,25.04,22.00,25.59,                  
     >            24.06,21.59,20.40,19.39,18.17,18.40,                  
     >            17.19,16.80,16.80,15.10,15.70,13.20,                  
     >            10.60, 7.10, 4.00, 1.18, 0.32, 0.10,                  
     >             1.02, 0.14, .024, .004, .0004/                       
!                                                                       
      DATA SIGAN2/  18 * 0.00,                                          
     >            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/                       
!                                                                       
      DATA SIGIO /  18 * 0.00,                                          
     >             0.00, 0.00, 2.12, 4.18, 4.38, 4.23,                  
     >             4.28, 4.18, 4.18, 8.00,11.35,10.04,                  
     >            12.21,12.22,12.23,11.90,12.17,12.13,                  
     >            11.91,11.64,11.25,11.21, 9.64, 9.95,                  
     >             8.67, 7.70, 7.68, 6.61, 7.13, 6.05,                  
     >             5.30, 2.90, 1.60, 0.59, 0.16, 0.05,                  
     >             0.51, 0.07, .012, .002, .0002/                       
!                                                                       
      DATA SIGIO2/ 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.27, 0.00, 1.00,                  
     >            12.22, 2.50, 9.34, 4.69, 6.12, 9.39,                  
     >            11.05, 9.69, 8.59,23.81,23.00,22.05,                  
     >            25.94,25.80,26.10,25.04,22.00,25.59,                  
     >            24.06,21.59,20.40,19.39,18.17,18.40,                  
     >            17.19,16.80,16.80,15.10,15.70,13.20,                  
     >            10.60, 7.10, 4.00, 1.18, 0.32, 0.10,                  
     >             1.02, 0.14, .024, .004, .0004/                       
!                                                                       
      DATA SIGIN2/  18 * 0.00,                                          
     >             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/                       
!                                                                       
!  CO                                                                   
      DATA SIGACO/6*0.,1.92,3.53,5.48,8.02,10.02,11.70,11.01,           
     | 12.52,12.47,13.61,15.43,15.69,18.01,19.92,20.09,21.61,22.28,     
     | 22.52,22.41,18.42,18.60,19.78,25.59,24.45,25.98,26.28,15.26,     
     |  33.22,21.35,22.59,37.64,49.44,28.50,52.9,3*0.0,16*0./           
      DATA SIGICO/6*0.,1.92,3.53,5.48,8.02,10.02,11.70,11.01,           
     | 12.52,12.47,13.61,15.43,15.69,18.01,19.92,20.09,21.44,22.31,     
     | 21.38,21.62,16.93,16.75,17.01,17.04,16.70,17.02,12.17,9.20,      
     | 15.44,11.38,17.13,11.70,6*0.0,16*0./                             
!  CO2                                                                  
      DATA SIGACO2/6*0.,4.42,7.51,11.03,14.98,17.88,21.21,20.00,        
     | 23.44,23.44,23.88,25.70,25.81,27.52,28.48,29.27,31.61,33.20,     
     | 34.21,34.00,25.31,25.86,25.88,25.96,21.76,22.48,53.96,26.48,     
     | 21.79,31.83,12.84,49.06,70.89,29.91,34.41,3*0.0,16*0./           
      DATA SIGICO2/6*0.,4.42,7.51,11.03,14.98,17.88,21.21,20.00,        
     | 23.44,23.44,23.88,25.70,25.81,27.52,28.48,29.27,31.61,33.20,     
     | 34.21,34.00,20.16,21.27,21.14,21.72,17.71,17.02,50.39,20.00,     
     | 17.07,21.53,10.67,19.66,6*0.0,16*0./                             
!  O+(4S)                                                               
      DATA SIGIO4S/6*0.,.32,1.03,1.62,1.95,2.15,2.33,2.23,2.23,         
     | 2.45,2.61,2.81,2.77,2.99,3.15,3.28,3.39,3.50,3.58,3.46,3.67,     
     | 3.74,3.73,4.04,4.91,4.20,4.18,4.18,4.28,4.23,4.38,4.18,2.12,     
     | 5*.00,16*0./                                                     
!  O+(2D)                                                               
      DATA SIGIO2D/6*0.,.34,1.14,2.00,2.62,3.02,3.39,3.18,3.62,         
     | 3.63,3.98,4.37,4.31,4.75,5.04,5.23,5.36,5.47,5.49,5.30,5.51,     
     | 5.50,5.50,5.52,6.44,3.80,12*.00,16*0./                           
!  O+(2P)                                                               
      DATA SIGIO2P/6*0.,.22,.75,1.30,1.70,1.95,2.17,2.04,2.32,          
     | 2.32,2.52,2.74,2.70,2.93,3.06,3.13,3.15,3.16,3.10,3.02,3.05,     
     | 2.98,2.97,.47,14*.00,16*0./                                      
!  O+(4P)                                                               
      DATA SIGIO4P/6*0.,.10,.34,.58,.73,.82,.89,.85,.91,.91,.93,        
     | .92,.92,.55,24*.00,16*0./                                        
!  O+(2P*)                                                              
      DATA SIGIO2Q/6*0.,.03,.27,.46,.54,.56,.49,.52,.41,.41,            
     | 28*.00,16*0./                                                    
!  HE+                                                                  
      DATA SIGIHE/6*0.,.21,.53,1.02,1.71,2.16,2.67,2.38,                
     | 3.05,3.05,3.65,4.35,4.25,5.51,6.53,7.09,.72,21*.00,16*0./        
!  N+  BANKS AND KOCKARTS X1.E-17                                       
      DATA SIGIN/6*0.,0.1,0.2,0.25,0.35,0.4,0.5,0.5,0.6,0.6,0.65,       
     |0.8,0.7,1.,1.,1.,1.1,1.15,1.2,1.1,1.2,1.2,1.2,1.2,1.2,1.1,1.1,1.1,
     |1.0,1.0,1.0,1.0,1.0,5*0.,16*0./                                   
!  H+  BANKS AND KOCKARTS X1.E-18                                       
      DATA SIGIH/6*0.,0.05,0.02,0.05,0.12,0.16,0.20,0.23,0.27,          
     |0.27,0.36,0.44,0.53,0.8,0.9,1.0,1.4,1.6,1.8,1.8,2.,2.2,2.3,2.8,   
     |3.1,3.5,3.8,4.,4.,4.1,4.8,5.8,6.12,5*0.,16*0./                    
!                                                                       
      sflux=0.
      if (iscale == 0) then
        r1 =  b1(1) + b1(2)*(f107a-71.5) + b1(3)*(f107-f107a+3.9)
        r2 =  b2(1) + b2(2)*(f107a-71.5) + b2(3)*(f107-f107a+3.9)
        do l=1,lmax
          if (l < 9) then
            sflux(l) = sra(l) * f107 + srb(l)
          else
            if (l .eq. 12) then
              sflux(l) = 332. + 0.6 * (F107-65.)
            else
              sflux(l) = (rflux(l) + ((r1-1.)*scale1(l)
     |                              + (r2-1.)*scale2(l)) / 1000.)
              if (sflux(l) < 0.0) sflux(l) = 0.0
            endif
          endif
        enddo
      else
        frat = (f107-68.) / (243.-68.)
        do l=1,lmax
          sflux(l) = rflux(l) + (xflux(l)-rflux(l)) * frat
        enddo
      endif

      sfnorm = 2. - (f107-68.) / (243.-68.)
      if (sfnorm < 1.0) sfnorm = 1.0
!
      do l=1,lmax
        wave1(l) = wavel(l)
        wave2(l) = waves(l)
        if (wave1(l) < 251. .and. wave2(l) > 49.)
     |     sflux(l) = sflux(l) * sfnorm
        sflux(l) = sflux(l) * 1.e9
      enddo

      if (first) then
        first = .false.
        do l=1,lmax
          m=lmax-l+1
          sigabs(1,l) = sigao(l)  * 1.e-18
          sigabs(2,l) = sigao2(l) * 1.e-18
          sigabs(3,l) = sigan2(l) * 1.e-18
          sigion(1,l) = sigio(l)  * 1.e-18
          sigion(2,l) = sigio2(l) * 1.e-18
          sigion(3,l) = sigin2(l) * 1.e-18
          sigabs(4,l) = sigaco(m) * 1.e-18
          sigion(4,l) = sigico(m) * 1.e-18
          sigabs(5,l) = sigaco2(m)* 1.e-18
          sigion(5,l) = sigico2(m)* 1.e-18
          sigion(6,l) = sigio4s(m)* 1.e-18
          sigion(7,l) = sigio2d(m)* 1.e-18
          sigion(8,l) = sigio2p(m)* 1.e-18
          sigion(9,l) = sigio4p(m)* 1.e-18
          sigion(10,l)= sigio2q(m)* 1.e-18
          sigion(11,l)= sigihe(m) * 1.e-18
          sigion(12,l)= sigin(m)  * 1.e-17
          sigion(13,l)= sigih(m)  * 1.e-18
        enddo
      endif
      end subroutine ssflux
!-----------------------------------------------------------------------
      subroutine column(zpht,xno,xno2,tn,xnn2,sht,gz,cno,cno2,cnn2,kmx)
!
! Args:
      integer,intent(in) :: kmx
      real,dimension(kmx),intent(in) :: zpht,xno,xno2,tn,xnn2,sht,gz
      real,dimension(kmx),intent(out) :: cno,cno2,cnn2
!
! Local:
      integer :: k,k1
      real :: zqht(kmx),shtcp,alp1,alp2,alp3

      do k=1,kmx
        zqht(k) = zpht(k)*1.e+5
      enddo
      shtcp=1.38e-16*tn(kmx)/(1.66e-24*gz(kmx))
      cno2(kmx)=xno2(kmx)*shtcp/32.
      cno(kmx)=xno(kmx)*shtcp/16.
      cnn2(kmx)=xnn2(kmx)*shtcp/28.

      do k=kmx-1,1,-1
        k1=k+1
        alp1=alog(xno2(k1)/xno2(k))/(zqht(k1)-zqht(k))
        alp2=alog(xno(k1)/xno(k))/(zqht(k1)-zqht(k))
        alp3=alog(xnn2(k1)/xnn2(k))/(zqht(k1)-zqht(k))
        cno2(k)=cno2(k1)+xno2(k)*(exp(alp1*(zqht(k1)-zqht(k)))-1.)/alp1
        cnn2(k)=cnn2(k1)+xnn2(k)*(exp(alp3*(zqht(k1)-zqht(k)))-1.)/alp3
        if (abs(alp2) < 1.e-10) then
          cno(k)=cno(k1)+xno(k)*(zqht(k1)-zqht(k))
        else
          cno(k)=cno(k1)+xno(k)*(exp(alp2*(zqht(k1)-zqht(k)))-1.)/alp2
        endif
      enddo ! k=kmx-1,1,-1

!     write(6,"('column returning: cno=',/,(6e12.4))") cno
!     write(6,"('column returning: cno2=',/,(6e12.4))") cno2
!     write(6,"('column returning: cnn2=',/,(6e12.4))") cnn2

      end subroutine column
!-----------------------------------------------------------------------
      subroutine tancomp(glat,glon,zg,xmo2,xmn2,xmo,tmm,cngo2,cngn2,cngo
     |  ,zpht,xno,xno2,tn,xnn2,sht,cno,cno2,cnn2,kmx)
!
! Args:
      integer,intent(in) :: kmx
      real,intent(in) :: zg
      real,intent(in) :: glat,glon ! for debug only
      real,dimension(kmx),intent(in) :: zpht,xno,xno2,tn,xnn2,sht,
     |  cno,cno2,cnn2
      real,intent(out) :: xmo2,xmn2,xmo,tmm,cngo2,cngn2,cngo
!
! Local:
      integer :: k,kk
      real :: zu,abd

      zu=zg*1.e-5
      if (zu < zpht(1)) then
        xmo2=xno2(1)*exp((-zg+zpht(1)*1.e+5)/sht(1))
        xmn2=xnn2(1)*exp((-zg+zpht(1)*1.e+5)/sht(1))
        xmo=xno(1)
        tmm=tn(1)
        cngo2=xmo2*sht(1)
        cngn2=xmn2*sht(1)
        cngo=cno(1)
      else ! zu >= zpht(1)
        kk = 0
        do k=1,kmx-1
          if (zu >= zpht(k).and.zu <= zpht(k+1)) then
            kk = k
            exit
          endif
        enddo
!
! If zu > zpht(kmx), then the above conditional did not go true,
!   and k==kmx on exit from the above loop (kk==0), then the abd 
!   statement below will be out-of-bounds referencing zpht(k+1). 
! In this case, I am setting kk=kmx-1 to avoid this problem.
! Open up the print statement below if you want to see when this happens.
!
        if (kk == 0) then ! k==kmx
!         write(6,"('>>> WARNING tancomp: zu > zpht(kmx): zu=',e15.7,
!    |      ' zpht(kmx)=',e15.7,' glat=',f8.2,' glon=',f8.2)") 
!    |      zu,zpht(kmx),glat,glon
!         write(6,"('Am stopping to avoid referencing zpht(kmx+1)')")
!         call shutdown('tancomp (solred module)')
          kk = kmx-1
        endif

        k = kk
        abd=(zu-zpht(k))/(zpht(k+1)-zpht(k))
        xmo2=xno2(k)*exp(alog(xno2(k+1)/xno2(k))*abd)
        xmn2=xnn2(k)*exp(alog(xnn2(k+1)/xnn2(k))*abd)
        xmo=xno(k)*exp(alog(xno(k+1)/xno(k))*abd)
        tmm=tn(k)+(tn(k+1)-tn(k))*(zu-zpht(k))/(zpht(k+1)-zpht(k))
        cngo2=cno2(k)*exp(alog(cno2(k+1)/cno2(k))*abd)
        cngn2=cnn2(k)*exp(alog(cnn2(k+1)/cnn2(k))*abd)
        cngo=cno(k)*exp(alog(cno(k+1)/cno(k))*abd)
      endif ! zu

      end subroutine tancomp
!-----------------------------------------------------------------------
      end module solred_module
