C       ZDUN_97.FOR
C
        SUBROUTINE PHOTS ( NR, NL, Len_t, NH, NABS, NO2, NH2O, GAS,
     *  ALT, wb, we, wlt, SZA, ALB, EXT, SSAM, GM, SIG, CONC, AIR,  !'T' not used  GM 1/00
     *  F8, HH, YCOL, E0, E, QF, EGR, PJ, xs_all, Num_beg,
     *  Len_X )
c >>>  Fluxes obtained by iterations <<<
c ---------------------------------------------------------------------
c         Calculation of photodissociation rates
c                    spherical geometry
c                 SZA > 90 deg. allowed !!
c
c  >>>  input:
c
C     NR     -  number of the photodissociation processes
C     NL     -  number of spectal grid points to the right of 200 nm
C
C     NH     -  number of altitude levels
C
C     NHU    -  number of levels where the contribution of
c               wavelengths < 200 nm (SRB) is considered
c     NABS   -  number of the absorbers
c     NO2    -  number of the reaction of the O2 photolysis
c     NH2O   -  number of the reaction of the H2O photolysis
c     ALT    -  altitude grid
c     WLT    -  spectral grid (wavelength) array
c     SZA    -  solar zenith angle
C     ALB    -  surface albedo
c     EXT    -  extinction coef. of Mie particles
c     SSAM   -  single scat. albedo  --"--
c     GM     -  asymm. factor  --"--
C     SIG    -  absorption cross-sections of the absorbers (O2 and O3)
C
c     SECH   -  effective cross-sections of photodissociating gases
c               ( x quantum yield )
c     CONC   -  concentrations of the absorbers
c     AIR    -  air density
C     F8     -  insolation at the upper boundary (photons/(cm^2s))
c     HH     -  height of homogeneous atm. for each absorber
C     LBEG   -  numbers of the lower photodiss. boundaries of each comp.
C     NIR    -  length of the photolysis spectrum for each comp.
C
c     YCOL   -  O2 and O3 opt. paths from each level to space
c ----------------------------------------------------------------------
c
c  >>>  output:
c
C     QF     -  actinic flux in the i-th spectral interval, photon/cm^2/s
C
c     E0     -  direct solar flux (with cosine) at each level
c     E      -  flux (with cosine) at each level
c     EGR    -  flux gradient for calc. of heating rates
C     PJ     -  photolysis rates
C
c ========================================================================
c
      include 'jpl.par'

      REAL  ALT(NH), WLT(Len_t), EXT(NH,NL), SSAM(NH,NL), GM(NH,NL),
     *      SIG(NABS,NH,Len_t), XS_all(Len_t,NH,NR), CONC(NH,NABS),
     *      HH(NABS), E0(NH,NL), E(NH,NL), QF(NL,NH), EGR(NH,NL),
     *      PJ(NR,NH), F8(Len_t), AIR(NH),
     *      SF1(Nlev), SF2(Nlev), SF3(Nlev), ! 'CFD(Nlev)' NOT USED GM 1/00
     *      YCOL(NH,NABS), 
     *      DW(Len_tm), wb(Len_t), we(Len_t),   ! 'T(NH)' NOT USED GM 1/00
     $      path(Nlev,Nlev)

      REAL dto2(Nlev), xso2(Nlev)

      integer   Num_beg(NR), Len_X(NR)
      character*32  GAS(NR)

      if ( SZA .le. 90. ) then
         if ( SZA .gt. 89.9 ) SZA = 89.9

c ++++  8/18/97  ++++++++++++++++++++++++++++++++++++++++++++++++++

         AMU0 = COSD(SZA)

         Low_lev=1         
         NHU=NH-Low_lev+1

C        Geometric paths from each level to space:
         DO IH = 1, NH
            NAH = NH-IH+1
            CALL  GPATH ( NAH, AMU0, ALT(IH), PATH(IH,IH) )
         ENDDO
         CALL  YCOLM ( NH, AMU0, ALT, HH, CONC, NABS, PATH, YCOL )

         do iw=1,NSRB
            call Dschuk ( NHU, YCOL(Low_lev,1), iw, dto2(Low_lev),
     $                  xso2(Low_lev) )
            do ih=Low_lev, NH
               sig(1,ih,iw)=dto2(ih)
               xs_all(iw,ih,NO2)=xso2(ih)
            enddo
         enddo

c ++++  8/18/97   +++++++++++++++++++++++++++++++++++++

         CALL  FSPEC (NH,NL,NABS,WLT,CONC,AIR,SIG,
     *        HH,PATH,EXT,SSAM,GM,ALT,F8,SZA,ALB,E0,E,
     *        QF,EGR )

      else if ( SZA .gt. 90. ) then
         if ( SZA .lt. 90.1 ) SZA = 90.1

c ++++  8/18/97  ++++++++++++++++++++++++++++++++++++++++++++++++++

         AMU0 = COSD(180.-SZA)

C        Geometric paths from each level to space:

         Low_Lev = IMIN ( NH, ALT, AMU0 )
         NHU=NH-Low_lev+1

         DO IH = Low_Lev, NH
            CALL  GPATH_o ( IH, NH, AMU0, ALT, PATH(1,IH) )
         ENDDO
c >>   O2 and O3 optical paths from each level to space:

         CALL  YCOLM_o ( NH, Low_Lev, AMU0, ALT, HH, CONC, NABS,
     *                PATH, YCOL )

         do iw=1,NSRB
            call Dschuk ( NHU, YCOL(Low_lev,1), iw, dto2(Low_lev),
     $                  xso2(Low_lev) )
            do ih=Low_lev, NH
               sig(1,ih,iw)=dto2(ih)
               xs_all(iw,ih,NO2)=xso2(ih)
            enddo
         enddo

c ++++  8/18/97   +++++++++++++++++++++++++++++++++++++

         CALL  FSPEC_o ( NH, NL, Low_Lev, NABS, WLT,
     *        CONC, AIR, SIG, HH, PATH, EXT, SSAM, GM, ALT,
     *        F8, SZA, ALB, E0, E, QF, EGR )
      endif
C
C     Photorates for O2, NO and H2O from Low_Lev and up:
C
      NUP1 = NH - Low_Lev + 1
      CALL  NIC1S ( NH, NUP1, NABS, YCOL,  SF1,SF2,SF3 )

c      nachos=25
c      write(*,777) ( i, wlt(i), f8(i), qf(i,nachos),
c     $          xs_all(i,nachos,NO2), i=1,29 )

c 777  format(i4,0pf8.2,1pe10.2,1pe10.2,1pe10.2)

c -----------------------------

      DO il = 1, Len_t
         dw(il) = we(il) - wb(il)
      ENDDO

      NNO = 0
      NH2O = 0
      NSO2 = 0
      NH2SO4 = 0
      NSO3=0

      DO IR=1,NR
         IF ( GAS(IR) .EQ. 'HCL     -->H       +CL' ) NHCL=IR
         IF ( GAS(IR) .EQ. 'H2O     -->H       +OH' ) THEN
               NH2O = IR
               GOTO 70
         ELSE IF ( GAS(IR) .EQ. 'NO      -->N       +O' ) THEN
c         IF ( GAS(IR) .EQ. 'NO      -->N       +O' ) THEN
               NNO = IR
               GOTO 70
         ELSE IF ( GAS(IR) .EQ. 'CH4     -->CH3     +H' ) THEN
               NCH4 = IR
               GOTO 70
C         ELSE IF ( GAS(IR) .EQ. 'SO2     -->SO      +O' ) THEN
C               NSO2 = IR
C               GOTO 70
         ELSE IF ( GAS(IR) .EQ. 'H2SO4   -->SO2     +H2O2' ) THEN
               NH2SO4 = IR
               GOTO 70
C         ELSE IF ( GAS(IR) .EQ. 'SO3     -->SO2     +O' ) THEN
C               NSO3 = IR
C               GOTO 70
         ENDIF
         DO IH=1,NH
            CALL PHRT(DW(Num_beg(ir)),len_X(ir),qf(Num_beg(ir),IH),
     $            XS_all(Num_beg(ir),IH,ir),PJ(ir,IH) )
         enddo

  70     CONTINUE
      ENDDO

c ----------------------
      NHU = NH-Low_Lev+1
      NHD = Low_Lev-1
    
      IF ( NH2O .NE. 0 ) then
         DO ih=1,NHD
            PJ(NH2O,ih)=0.
         enddo
         do ih = NHD+1, NH
             PJ(NH2O,IH)=SF3(IH-NHD)
         enddo
      endif
      IF ( NNO .NE. 0 ) then
         DO ih=1,NHD
            PJ(NNO,ih)=0.
         enddo
         do ih = NHD+1, NH
            PJ(NNO,IH)=SF2(IH-NHD)
         enddo
      endif
C      IF ( NSO2 .NE. 0 ) then
C         DO ih=1,NH
C            PJ(NSO2,ih)=0.
C         enddo
C      endif
C      IF ( NSO3 .NE. 0 ) then
C         DO ih=1,NH
C            PJ(NSO3,ih)=0.
C         enddo
C      endif
      IF ( NH2SO4 .NE. 0 ) then
         DO ih=1,NH
            PJ(NH2SO4,ih)=PJ(NHCL,ih)*0.3
         enddo
      endif
C
C    For 0-60km range, PJ(CH4) is practically 0, hence:
C
      IF ( NCH4 .NE. 0 ) then
         DO ih=1,NH
            PJ(NCH4,ih)=0.
         enddo
      endif
      RETURN
      END



