C     FLUX_N.FOR
C
      SUBROUTINE FSPEC ( NH, NL, NABS, WLT, CONC, AIR, SIG,
     *   HH, PATH, EXT, SSAM, GM, ALT, F8, SZA, ALB, E0, E, QF, EGR )
c------------------------------------------------------------------
c       Calculation of solar fluxes for the whole model spectrum
c
c  >>>  input:
c
c    NH   - number of levels
c    NL   - number of spectral gridpoints
c    NABS - number of the absorbers
c    WLT  - wavelength grid
c    CONC - concentrations of the absorbers
c    AIR  - air density
c    SIG  - abs. cross sections of the absorbers
c    HH   - height of homogeneous atm. for each absorber
c    YCOL - O2 and O3 opt. paths from each level to space
c    EXT  - extinction coef. of Mie particles
c    SSAM - single scat. albedo  --"--
c    GM   - asymm. factor  --"--
c    ALT  - model altitudes
c    F8   - TOA solar flux
c    SZA  - zenith angle
c    ALB  - surface albedo
c
c -------------------------------------------------------------
c  >>>  output:
c
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c
c ----------------------------------------------------------------
c  >>>  internal:
c
C    DZ     - thickness of each layer
c    TAU    - opt. depths of the layers
c    SSA    - single scat. alb.  ---"---
c    GE     - asymm. factor  ---"---
C    TAUD   - reduced opt. path from each level to space
C    QFF    - actinic flux top to bottom at given wlt
c ==================================================================
c
      include 'jpl.par'

      REAL  WLT(NL), CONC(NH,NABS), AIR(NH), SIG(NABS,NH,NL),
     *      EXT(NH,NL),
     *      SSAM(NH,NL), GM(NH,NL), F8(NL), ALT(NH), HH(NABS),
     *      QF(NL,NH), E(NH,NL), EGR(NH,NL), E0(NH,NL),
     *      DZ(Nlev), TAU(Nlev), SSA(Nlev), GE(Nlev), TAUD(Nlev),
     *      QFF(Nlev), PATH(NH,NH), CC(MNABS)
      AMU0 = COSD(SZA)
      AMU = SQRT(1224.*AMU0*AMU0 + 1.)/35.

c      NL1=11
c      IF(NL.LT.11) NL1=NL
      DO IH = 1, NH-1
         DZ(IH) = ALT(IH+1) - ALT(IH)
      ENDDO
C
C  >>             Geometric paths from each level to space:
C
c      DO IH = 1, NH
c         NAH = NH-IH+1
c         CALL  GPATH ( NAH, AMU0, ALT(IH), PATH(IH,IH) )
c      ENDDO
c      CALL  YCOLM ( NH, AMU0, ALT, HH, CONC, NABS, PATH, YCOL )

C  \
C ***   Calculation of fluxes with the O2 absorption O2 (up to 245 nm)
C  /
      DO 1 L = 1, NL
         SIGR=SR(WLT(L))
         E8 = F8(L)*AMU
         DO IH = 1, NH-1
            DO IA = 1, NABS
               CC(IA) = CONC(IH,IA)
            ENDDO
            CALL  OPTPAR ( DZ(IH), CC, SIG(1,IH,L), NABS, AIR(IH),
     *        SIGR, EXT(IH,L), SSAM(IH,L), GM(IH,L), TAU(IH), SSA(IH),
     *        GE(IH) )
         ENDDO
         CALL   TAUDIR ( NH, TAU, SSA, GE, AMU0, ALT, HH, SIG(1,1,L),
     *          CONC, NABS, PATH, TAUD) !'sigr' and 'air' NOT USED GM1/00
         CALL   FZMON ( NH, AMU, ALB, TAU, TAUD, SSA, GE, E8,
     *          E0(1,L), E(1,L), QFF, EGR(1,L) )
         DO J=1,NH
            QF(L,J)=QFF(J)
         ENDDO
 1     CONTINUE
c      IF ( NL.LT.12 ) GO TO 5

C
C ***  Calculation of fluxes without O2 abs.(245 nm and up)
C
c      DO 2 L = 12, NL
c         SIGR=SR(WLT(L))
c         E8 = F8(L)*AMU
c         DO IH = 1, NH-1
c            CALL  OPTPAR ( DZ(IH), CONC(IH,2), SIG(2,IH,L), 1, AIR(IH),
c     *        SIGR, EXT(IH,L), SSAM(IH,L), GM(IH,L), TAU(IH), SSA(IH),
c     *        GE(IH) )
c         ENDDO
c         do ih = 1, nh
c            sigO3(ih) = SIG(2,ih,L)
c         enddo
c         CALL   TAUDIR ( NH, TAU, SSA, GE, AMU0, ALT, HH(2), sigO3,
c     *          CONC(1,2), 1, PATH, TAUD, sigr, air )
c         CALL   FZMON ( NH, AMU, ALB, TAU, TAUD, SSA, GE, E8,
c     *          E0(1,L), E(1,L), QFF, EGR(1,L) )
c         DO J=1,NH
c            QF(L,J)=QFF(J)
c         ENDDO
c 2       CONTINUE
c 5    CONTINUE
      RETURN
      END
c ---------------------------------------------------------------
      SUBROUTINE FZMON ( NH, AMU, ALB, TAU, TAUD, SSA, GE, E8,
     *    E0, E, QF, EGR )
c ================================================================
c
C **  Calculation of monochromatic fluxes, top to bottom   ****
C     method delta-Eddington, after Zdunkowski et.al. 1980
c
c >> input parameters:
c
c    NH     - # of alt. levels
c    AMU  - cos SZA
c    ALB    - surface albedo
C    TAU    - NORMAL optical depths of the layers
c    TAUD  -  reduced (by factor (1-om*f)) slant opt.
c             paths from TOA to i-th level
c    SSA    - Single Scat. Albedo of the layers
c    GE     - asymmetry factor g
c    E8     - flux (with cosine) at the TOA
c ---------------------------------------------------------------
c >> output parameters:
c
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c
c ===================================================================
c >> internal:
c
c    T0,R0,TD,RD - direct and diffuse layer transm. and refl.
c
c ===================================================================

      include 'jpl.par'

      REAL  TAU(NH), TAUD(NH), SSA(NH), GE(NH), E0(NH), E(NH),
     *QF(NH), EGR(NH)
      REAL   T0(Nlev), R0(Nlev), TD(Nlev), RD(Nlev)
      DOUBLE PRECISION U1,TT,G,F,OM,U2,BETA,BETA0,AL1,AL2,AL3,AL4,EPS,
     *CM,CE,ZN1,ZN2
      U1 = 7.D0/4.D0
      E0(NH)=E8*EXP(-TAUD(NH))
      NH1=NH-1
      DO 1 II=1,NH1
         IH=NH-II
         TT = TAU(IH)
         G = GE(IH)
         F = G*G
         OM = SSA(IH)
         IF ( abs(1.D0-OM) .le. .00001D0 ) OM = 0.99999D0
C ***  Calculation of T and R for each layer:  ********
c                          ...  coefficients:
         BETA0 = 3.D0/7.D0*(1.D0-G)
         BETA  = .5D0 - 3.D0*AMU*G/(4.D0*(1.D0+G))
 922     continue
         U2 = U1*( 1.D0- (1.D0-OM)/(7.D0*OM*BETA0 ))
c !!!! Correction 8/27/95 : if OM is small, U2 can be < 0 which
c      leads to negative diffuse fluxes; hence,

         if ( U2 .le. 1.D0 ) U2 = U1
         AL1 = U1*( 1.D0- OM*(1.D0-BETA0) )
         AL2 = U2*OM*BETA0
         AL3 = (1.D0-F)*OM*BETA
         AL4 = (1.D0-F)*OM*(1.D0-BETA)
         EPS = DSQRT( (AL1-AL2)*(AL1+AL2) )
         if ( om*f+eps*amu .eq. 1.D0 .or. om*f-eps*amu .eq. 1.D0 ) then
            om = om*.99D0
            goto 922
         endif
         CM  = AL2/(AL1+EPS)
         CE  = DEXP(-EPS*TT)
         ZN1=1.D0/( (1.D0-OM*F - EPS*AMU )*( 1.D0-OM*F + EPS*AMU) )
         GAM1 = ( (1.D0-OM*F)*AL3 - AMU*(AL1*AL3+AL2*AL4) )*ZN1
         GAM2 = ( -(1.D0-OM*F)*AL4 - AMU*(AL1*AL4+AL2*AL3) )*ZN1
         tau_scld = TAU(IH)*(1.-OM*F)/AMU
         if ( tau_scld .le. 20. ) then
             TDB  = EXP( - tau_scld )
         else
             TDB = 0.
         endif
c           ... now, we can calculate T, R and dir. flux:
         ZN2=1./((1.D0-CE*CM)*(1.D0+CE*CM))
         if ( AL2 .lt. 0. ) then
            RD(IH) = 0.
         else
            RD(IH) = CM*(1.D0-CE)*(1.D0+CE)*ZN2
         endif
         TD(IH) = CE*(1.D0-CM)*(1.D0+CM)*ZN2
         R0(IH) = -GAM2*RD(IH) - GAM1*TDB*TD(IH) + GAM1
         T0(IH) = -GAM2*TD(IH) - GAM1*TDB*RD(IH) + GAM2*TDB
C **   Calculation of direct fluxes
         if ( TAUD(IH) .le. 20. ) then
            E0(IH)=E8*EXP(-TAUD(IH))
         else
            E0(IH) = 0.
         endif
 1    CONTINUE
      call  flx_iter ( NH, AMU, ALB, R0, T0, RD, TD,
     *                        E8, E0, E, QF, EGR )
      END

      SUBROUTINE FSPEC_o ( NH, NL, Low_Lev, NABS, WLT, CONC, AIR,
     * SIG,HH, PATH, EXT,SSAM, GM,ALT, F8,SZA, ALB, E0, E, QF, EGR )
c------------------------------------------------------------------
c       Calculation of solar fluxes for the whole model spectrum;
c                   SZA > 90 deg.
c                   ============
c  >>>  input:
c
c    NH   - number of levels
c    NL   - number of spectral gridpoints
c    NABS - number of the absorbers
c    WLT  - wavelength grid
c    CONC - concentrations of the absorbers
c    AIR  - air density
c    SIG  - abs. cross sections of the absorbers
c    HH   - height of homogeneous atm. for each absorber
c    YCOL - O2 and O3 opt. paths from each level to space
c    EXT  - extinction coef. of Mie particles
c    SSAM - single scat. albedo  --"--
c    GM   - asymm. factor  --"--
c    ALT  - model altitudes
c    F8   - TOA solar flux
c    SZA  - zenith angle
c    ALB  - surface albedo
c
c -------------------------------------------------------------
c  >>>  output:
c
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c
c ----------------------------------------------------------------
c  >>>  internal:
c
C    DZ     - thickness of each layer
c    TAU    - opt. depths of the layers
c    SSA    - single scat. alb.  ---"---
c    GE     - asymm. factor  ---"---
C    TAUD   - reduced opt. path from each level to space
C    QFF    - actinic flux top to bottom at given wlt
c ==================================================================
c
       include 'jpl.par'

       REAL  WLT(NL), CONC(NH,NABS), AIR(NH), SIG(NABS,NH,NL),
     *       EXT(NH,NL),
     *       SSAM(NH,NL), GM(NH,NL), F8(NL), ALT(NH), HH(NABS),
     *       QF(NL,NH), E(NH,NL), EGR(NH,NL), E0(NH,NL),
     *       DZ(Nlev), TAU(Nlev), SSA(Nlev), GE(Nlev), TAUD(Nlev),
     *       QFF(Nlev), PATH(NH,NH), CC(MNABS)
      AMU0 = COSD(180.-SZA)
      AMU = SQRT(1224.*AMU0*AMU0 + 1.)/35.
c      Low_Lev = IMIN ( NH, ALT, AMU0 )
c      NL1=11
c      IF(NL.LT.11) NL1=NL
      DO IH = 1, NH-1
          DZ(IH) = ALT(IH+1) - ALT(IH)
      ENDDO

c  >>             Geometric paths from each level to space:
c      DO IH = Low_Lev, NH
c         CALL  GPATH_o ( IH, NH, AMU0, ALT, PATH(1,IH) )
c      ENDDO
c  >>           O2 and O3 optical paths from each level to space:
c      CALL  YCOLM_o ( NH, Low_Lev, AMU0, ALT, HH, CONC, NABS,
c     *                PATH, YCOL )
c  \
C ***   Calculation of fluxes with the O2 absorption O2 (up to 245 nm)
c  /
      DO 1 L = 1, NL
         SIGR=SR(WLT(L))
         E8 = F8(L)*AMU
         DO IH = 1, NH-1
            DO IA = 1, NABS
               CC(IA) = CONC(IH,IA)
            ENDDO
            CALL  OPTPAR ( DZ(IH), CC, SIG(1,IH,L), NABS, AIR(IH),
     *        SIGR, EXT(IH,L), SSAM(IH,L), GM(IH,L), TAU(IH), SSA(IH),
     *        GE(IH) )
         ENDDO
         CALL   TAUDIR_o ( NH, Low_Lev, TAU, SSA, GE, AMU0, ALT,
     *          HH, SIG(1,1,L), CONC, NABS, PATH, TAUD) 
C                                    'sigr' and 'air' NOT USED GM1/00
         CALL   FZMON_o ( NH, Low_Lev, AMU, ALB, TAU, TAUD, SSA,
     *          GE, E8, E0(1,L), E(1,L), QFF, EGR(1,L) )
          DO J=1,NH
            QF(L,J)=QFF(J)
          ENDDO
 1     CONTINUE
c      IF ( NL.LT.12 ) GO TO 5
C ***  Calculation of fluxes without O2 abs.(245 nm and up)
c      DO 2 L = 12, NL
c         SIGR=SR(WLT(L))
c         E8 = F8(L)*AMU
c         DO IH = 1, NH-1
c            CALL  OPTPAR ( DZ(IH), CONC(IH,2), SIG(2,IH,L), 1, AIR(IH),
c     *        SIGR, EXT(IH,L), SSAM(IH,L), GM(IH,L), TAU(IH), SSA(IH),
c     *        GE(IH) )
c         ENDDO
c         do ih = 1, nh
c            sigO3(ih) = SIG(2,ih,L)
c         enddo
c         CALL   TAUDIR_o ( NH, Low_Lev, TAU, SSA, GE, AMU0, ALT,
c     *  HH(2), sigO3, CONC(1,2), 1, PATH, TAUD, sigr, air )
c         CALL   FZMON_o ( NH, Low_Lev, AMU, ALB, TAU, TAUD, SSA,
c     *          GE, E8, E0(1,L), E(1,L), QFF, EGR(1,L) )
c          DO J=1,NH
c            QF(L,J)=QFF(J)
c          ENDDO
c 2     CONTINUE
c 5     CONTINUE
       RETURN
       END
c ============================================================
      SUBROUTINE FZMON_o ( NH, Low_Lev, AMU, ALB, TAU, TAUD,
     *    SSA, GE, E8, E0, E, QF, EGR )
c ================================================================
c
C **  Calculation of monochromatic fluxes, top to bottom   ****
C     method delta-Eddington, after Zdunkowski et.al. 1980
c                     SZA > 90 deg
c >> input parameters:
c
c    NH     - # of alt. levels
c    AMU  - cos SZA
c    ALB    - surface albedo
C    TAU    - NORMAL optical depths of the layers
c    TAUD  -  reduced (by factor (1-om*f)) slant opt.
c             paths from TOA to i-th level
c    SSA    - Single Scat. Albedo of the layers
c    GE     - asymmetry factor g
c    E8     - flux (with cosine) at the TOA
c ---------------------------------------------------------------
c >> output parameters:
c
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c
c ===================================================================
c >> internal:
c
c    T0,R0,TD,RD - direct and diffuse layer transm. and refl.
c
c ===================================================================

      include 'jpl.par'

      REAL  TAU(NH), TAUD(NH), SSA(NH), GE(NH), E0(NH), E(NH),
     *QF(NH), EGR(NH)
      REAL   T0(Nlev), R0(Nlev), TD(Nlev), RD(Nlev)
      U1 = 7./4.
      if ( TAUD(NH) .le. 20. ) then
          E0(NH)=E8*EXP(-TAUD(NH))
      else
          E0(NH) = 0.
      endif
      NH1=NH-1
      DO 1 II=1,NH1
         IH=NH-II
         TT = TAU(IH)
         G = GE(IH)
         F = G*G
         OM = SSA(IH)
         IF ( abs(1.-OM) .le. .00001 ) OM = 0.99999
C ***  Calculation of T and R for each layer:  ********
c                          ...  coefficients:
         BETA0 = 3./7.*(1.-G)
         BETA  = .5 - 3.*AMU*G/4./(1.+G)
 922     continue
         U2 = 7./4.*( 1.- (1.-OM)/7./OM/BETA0 )
c !!!! Correction 8/27/95 : if OM is small, U2 can be < 0 which
c      leads to negative diffuse fluxes; hence,

         if ( U2 .le. 1. ) U2 = U1

         AL1 = U1*( 1.- OM*(1.-BETA0) )
         AL2 = U2*OM*BETA0
         AL3 = (1.-F)*OM*BETA
         AL4 = (1.-F)*OM*(1.-BETA)
         EPS = SQRT( (AL1-AL2)*(AL1+AL2) )
         if ( om*f+eps*amu .eq. 1. .or. om*f-eps*amu .eq. 1. ) then
            om = om*.99
            goto 922
         endif
         CM  = AL2/(AL1+EPS)
         CE  = EXP(-EPS*TT)
         GAM1 = ( (1.-OM*F)*AL3 - AMU*(AL1*AL3+AL2*AL4) )/
     *        ( (1.-OM*F - EPS*AMU )*( 1.-OM*F + EPS*AMU) )
         GAM2 = ( -(1.-OM*F)*AL4 - AMU*(AL1*AL4+AL2*AL3) )/
     *        ( (1.-OM*F - EPS*AMU )*( 1.-OM*F + EPS*AMU) )
         tau_scld = TAU(IH)*(1.-OM*F)/AMU
         if ( tau_scld .le. 20. ) then
             TDB  = EXP( - tau_scld )
         else
             TDB = 0.
         endif
c           ... now, we can calculate T, R and dir. flux:
         if ( AL2 .lt. 0. ) then
            RD(IH) = 0.
         else
            RD(IH) = CM*(1.-CE)*(1.+CE)/(1.-CE*CM)/(1.+CE*CM)
         endif
         TD(IH) = CE*(1.-CM)*(1.+CM)/(1.-CE*CM)/(1.+CE*CM)
         R0(IH) = -GAM2*RD(IH) - GAM1*TDB*TD(IH) + GAM1
         T0(IH) = -GAM2*TD(IH) - GAM1*TDB*RD(IH) + GAM2*TDB
C **   Calculation of direct fluxes
         if ( TAUD(IH) .le. 20. .and. ih .ge. Low_Lev ) then
            E0(IH)=E8*EXP(-TAUD(IH))
         else
            E0(IH) = 0.
         endif
C
 1    CONTINUE
C ------------------------------------------------

      call  flx_iter_o ( NH, Low_Lev, AMU, ALB, R0, T0, RD, TD,
     *                        E8, E0, E, QF, EGR )

      END

c ===========================================

        Subroutine  flx_iter ( NH, AMU, ALB, R0, T0, RD, TD,
     *                         E8, E0, E, QF, EGR )
c
c    Calculation of fluxes by iterations
c
c >> input parameters:
c
c    NH     - # of alt. levels
c    AMU    - cos SZA
c    ALB    - surface albedo
c    T0,R0,TD,RD - direct and diffuse layer transm. and refl.
c ===
c >> output parameters:
c
c    E8     - flux (with cosine) at the TOA
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c ===
c >> internal parameters:
c
c    FUPold, FDold - upward and downward diffuse fluxes (i-1 iter.)
c    FUPnew, FDnew - upward and downward diffuse fluxes (i iter.)
c
c ====================================================================
c
         parameter ( MMNH = 100 )
         real   R0(NH), T0(NH), RD(NH), TD(NH), E0(NH),
     *          E(NH), QF(NH), EGR(NH), FUPold(MMNH), FUPnew(MMNH),
     *          FDold(MMNH), FDnew(MMNH)

        cutoff = 1.e-9      !  total transmittion cutoff
        eps = 0.01          !  iteration accuracy

        do i = 1, NH
           FDold(i) = 0.
           FUPold(i) = 0.
           FDnew(i) = 0.
           FUPnew(i) = 0.
        enddo

        NH1 = NH-1
        FDnew(NH) = 0.
c >> Begin iterations
        ITER = 0
 100    continue
c >> calculate FD going down:
        do ih = 1, NH1
           i = NH - ih
           FDnew(i) = FDnew(i+1)*TD(i) + FUPold(i)*RD(i) +
     +        E0(i+1)*T0(i)

c ... terminate calculations if (FDnew+FUPold++E0)/E8 < cutoff
           if ( (FDnew(i)+FUPold(i)+E0(i))/E8 .le.cutoff) then
             iterm = i
             goto 1
           endif
        enddo
        iterm = 1
  1     continue
 9     format (2x, i3, f8.5)

c >> Now, go up from the cutoff level:

       if ( iterm .eq. 1 ) then
         FUPnew(1) = ALB*( E0(1)+FDnew(1) )
       else
         FUPnew(iterm) = E0(iterm)*R0(iterm-1) +
     +                    FDnew(iterm)*RD(iterm-1)
       endif
       do i = iterm, NH1
          FUPnew(i+1) = FUPnew(i)*TD(i) + FDnew(i+1)*RD(i) +
     +            E0(i+1)*R0(i)
       enddo
c ==================================================================
c >> Check if the desired accuracy is reached:
       do i = iterm, NH
          ftot =  FDnew(i) + FUPnew(i)
          if ( ftot .gt. 100. ) then
             Ferr = ( ftot - FDold(i) - FUPold(i) )/ftot
          else
             Ferr = 0.
          endif

          if ( Ferr .gt. eps ) then
             do j = 1, NH
               FDold(j) = FDnew(j)
               FUPold(j) =  FUPnew(j)
             enddo
             ITER = ITER+1
             goto 100
          endif
        enddo
Cc ... calculate fluxes:

       do ih = 1, NH-iterm+1
          i = NH-ih+1
          EDIF =  FDnew(i) + FUPnew(i)
          E(i) = E0(i) + EDIF
          EGR(i) = E0(i)+FDnew(i)-FUPnew(i)
          QF(i) = 2.*EDIF + E0(i)/AMU
       enddo

       if ( iterm .gt. 1 ) then
          do i = 1, iterm
             E(i)   = 0.
             EGR(i) = 0.
             QF(i)  = 0.
          enddo
       endif

       end

c ============================================================

        Subroutine  flx_iter_o ( NH, Low_Lev, AMU, ALB, R0, T0,
     *                         RD, TD, E8, E0, E, QF, EGR )
c
c    Calculation of fluxes by iterations for SZA > 90 deg.
c
c >> input parameters:
c
c    NH     - # of alt. levels
c    AMU    - cos SZA
c    ALB    - surface albedo
c    T0,R0,TD,RD - direct and diffuse layer transm. and refl.
c ===
c >> output parameters:
c
c    E8     - flux (with cosine) at the TOA
c    E0     - direct solar flux (with cosine) at each level
c    E      - flux (with cosine) at each level
c    QF     - actinic flux (without cos)
c    EGR    - flux gradient
c ===
c >> internal parameters:
c
c    FUPold, FDold - upward and downward diffuse fluxes (i-1 iter.)
c    FUPnew, FDnew - upward and downward diffuse fluxes (i iter.)
c
c ====================================================================
c
         parameter ( MMNH = 100 )
         real   R0(NH), T0(NH), RD(NH), TD(NH), E0(NH),
     *          E(NH), QF(NH), EGR(NH), FUPold(MMNH), FUPnew(MMNH),
     *          FDold(MMNH), FDnew(MMNH)

        cutoff = 1.e-9      !  total transmittion cutoff
        eps = 0.01          !  iteration accuracy

        do i = 1, NH
           FDold(i) = 0.
           FUPold(i) = 0.
           FDnew(i) = 0.
           FUPnew(i) = 0.
        enddo

        NH1 = NH-1
c >> Begin iterations
        ITER = 0
        iterm = Low_Lev+1
 100    continue
c >> calculate FUP going from iterm to top:

       if ( iterm .eq. 1 ) then
         FUPnew(1) = ALB*FDold(1)
       else
         FUPnew(iterm) = E0(iterm-1)*T0(iterm-1) +
     +                    FDold(iterm)*RD(iterm-1)
       endif

       do i = iterm, NH1
          FUPnew(i+1) = FUPnew(i)*TD(i) + FDold(i+1)*RD(i) +
     +            E0(i)*T0(i)
       enddo

c >> now, calculate FD going from top to cutoff level:

        do ih = 1, NH1
           i = NH - ih
           FDnew(i) = FDnew(i+1)*TD(i) + FUPnew(i)*RD(i) +
     +        E0(i)*R0(i)

c ... terminate calculations if (FDnew+FUPold+E0)/E8 < cutoff
           if ( (FDnew(i)+FUPnew(i)+E0(i))/E8 .le. cutoff)then
             iterm = i
             goto 1
           endif
        enddo
        iterm = 1
  1     continue
c ==================================================================
c >> Check if the desired accuracy is reached:

       do i = iterm, NH
          ftot = FDnew(i)+FUPnew(i)
          if ( ftot .gt. 100. ) then
             Ferr = ( ftot - FDold(i) - FUPold(i) )/ftot
          else
             Ferr = 0.
          endif

          if ( Ferr .gt. eps ) then
             do j = 1, NH
               FDold(j) = FDnew(j)
               FUPold(j) =  FUPnew(j)
             enddo
             ITER = ITER+1
             goto 100
          endif
        enddo

c ... calculate fluxes:

       do ih = 1, NH-iterm+1
          i = NH-ih+1
          EDIF =  FDnew(i) + FUPnew(i)
          E(i) = E0(i) + EDIF
          EGR(i) = -E0(i)+FDnew(i)-FUPnew(i)
          QF(i) = 2.*EDIF + E0(i)/AMU
       enddo

       if ( iterm .gt. 1 ) then
          do i = 1, iterm
             E(i)   = 0.
             EGR(i) = 0.
             QF(i)  = 0.
          enddo
       endif

       end
        SUBROUTINE HI_LO_INT ( Nz, Nzz, z, zz, A, AA )
c
c -------------------------------------------------------
c
c        Linear interpolation from high resolution to
c        low resolution grid
c  !!!   arguments should be in increasing order  !!!
c  !!!   DIMENSION OF z is Nz+1 !!!!
c -------------------------------------------------------
c
        INTEGER Nz, Nzz
        REAL z(Nz+1), zz(Nzz), A(Nz), AA(Nzz)
c
        iz = 1
        do while ( z(iz) .le. zz(1) )
           iz = iz+1
        enddo

        do izz = 1, Nzz-1
           iz = iz-1
           zb = zz(izz)
           ze = zz(izz+1)
           As = 0.
           zs = z(iz)
           do while ( z(iz) .le. ze .and. iz .le. Nz )
              As = As + A(iz)*( z(iz+1)-z(iz) )
              iz = iz+1
            enddo
            AA(izz) = As/( z(iz)-zs )
         enddo

        END
        SUBROUTINE LINTERP ( Nz, Nzz, z, zz, A, AA )
c
c -------------------------------------------------------
c
c        Linear interpolation and extrapolation
c  !!!   arguments should be in increasing order  !!!
c
c -------------------------------------------------------
c
        INTEGER Nz, Nzz
        REAL z(Nz), zz(Nzz), A(Nz), AA(Nzz)
c
        INTEGER i,j
        REAL zz1, z1, z2
c
        DO  i = 1, Nzz
            IF (z(1).GT.zz(i)) THEN
              AA(i) = A(1) + (A(2)-A(1))*
     *    (zz(i)-z(1))/(z(2)-z(1))
              i0 = i
            ELSE
              i0 = 1
              GOTO 11
            END IF
        END DO
   11   CONTINUE
c
        zNz = z(Nz)
        DO i = i0, Nzz
         zz1 = zz(i)
         IF (zz1.LE.zNz) THEN
         DO j=1,Nz-1
          z1 = z(j)
          z2 = z(j+1)
          IF ( (zz1 .GE. z1) .AND. (zz1 .LE. z2) ) THEN
           AA(i) = A(j) + (A(j+1)-A(j))*(zz1-z1)/(z2-z1)
          ENDIF
         ENDDO
         ELSE
           AA(i) = A(Nz) + (A(Nz)-A(Nz-1))*
     *   (zz(i)-zNz)/(zNz-z(Nz-1))
         END IF
        ENDDO
        RETURN
        END
        SUBROUTINE LINTERP_LOG ( Nz, Nzz, z, zz, A, AA )
c
c -------------------------------------------------------
c
c        Linear interpolation and extrapolation
C        FOR EXPONENT-LIKE FUNCTIONS A(z)
c  !!!   arguments should be in increasing order  !!!
c
c -------------------------------------------------------
c
        INTEGER Nz, Nzz
        REAL z(Nz), zz(Nzz), A(Nz), AA(Nzz)
        INTEGER i
        DO  i = 1, Nz
            A(i) = ALOG(A(i))
        enddo
        CALL  LINTERP ( Nz, Nzz, z, zz, A, AA )
        DO  i = 1, Nzz
            AA(i) = EXP(AA(i))
        enddo
        RETURN
        END
      SUBROUTINE  OPTPAR ( DZ, CONC, SIG, NABS, AIR, SIGR, EXT,
     *                     SSAM, GM, TAU, SSA, GE )
C ------------------------------------------------------------
c     Program to calculate opt. parameters (tau, ssa and g)
c     of each layer
c ------------------------------------------------------------
c >>> input:
c
c   DZ    - layer geom. thickness
C   CONC  - concentrations of the absorbers
c   SIG   - cross-sections      --"--
c   NABS  - number of the absorbers
c   AIR   - air density
c   SIGR  - Rayleigh scat. cross-section
c   EXT   - extinction opt. depth of Mie scat.
c   SSAM  - SSA of Mie scat.
c   GM    - assym. factor of Mie scat.
c---------------------------------------------------------------
c >>> output:
c
c   TAU   - total opt. depth of the layer
c   SSA   - effect. SSA of the layer
c   GE    - effect. assym. fact. of the layer
c
c ----------------------------------------------------------------
c
       REAL  CONC(NABS), SIG(NABS)
       TAUA = 0.
       DO IA = 1, NABS
          TAUA = TAUA+CONC(IA)*SIG(IA)
       ENDDO
       TAUA = ( TAUA + EXT*(1.-SSAM) )*DZ
       TAUAIR = SIGR*AIR*DZ
       TAUMIE = EXT*SSAM*DZ
       TAUR   = TAUAIR + TAUMIE
       TAU = TAUA + TAUR
       SSA = TAUR/TAU
       GE = GM*TAUMIE/TAUR

       END
C --------------------------------------------------------------
       SUBROUTINE  TAUDIR ( NH, TAU, SSA, GE, AMU, ALT, HH, SIG,
     *                      CONC, NABS, PATH, TAUD) 
C                                  'sigr' and 'air' NOT USED GM 1/00
C --------------------------------------------------------------
C    Calculation of the reduced optical path for direct beam
c    from each alt. to space
c
C  >>>  input:
c
c    NH   - number of levels
c    TAU  - normal optical depths of the layers
c    SSA  - single scattering albedo of the layers
c    GE   - asymm. param. g  --"--
c    AMU  = cos(SZA)
c    ALT  - model altitudes
c    HH   - height of the homogeneous atm. at model top for each absorber
c    SIG  - abs. cross-section
c    CONC - concentrations of the absorbers
c    NABS - number of absorbers
c    PATH - geometric paths through each layer
c --------------------------------------------------------------
c  >>>  output:
c
c    TAUD  - opt. path for direct beam
c ===============================================================
       include 'jpl.par'

       REAL   TAU(NH), SSA(NH), GE(NH), ALT(NH), TAUD(NH), PATH(NH,NH),
     *        CONC(NH,NABS), F(MNH), DZ(MNH),
     *        HH(NABS), SIG(NABS,NH) ! 'air(nh)' NOT USED GM 1/00
       DATA   REAR/6.371E8/
        amu35 = sqrt(1224.*amu*amu+1.)/35.
       DO I = 1, NH-1
          F(I) = GE(I)*GE(I)
          DZ(I) = ALT(I+1) - ALT(I)
       ENDDO
       DO I = 1, NH-1
          TTT = 0.
          qq = 0.
          DO J = I, NH-1
             TTT = TTT + TAU(J)*(1.-SSA(J)*F(J))/DZ(J)*PATH(J,I)
          ENDDO
c     Calculate opt. depth above model top by Chapman's function:

c                ... first, SZA at model top:
         RL=REAR+ALT(I)
c  from sin theorem:
         RG=RL*SQRT(1.-AMU*AMU)
         RTOP = REAR + ALT(NH)
         VV=RG/RTOP
         AMU1=SQRT(1.-VV*VV)

         TTOP = 0.
         DO IA = 1, NABS
           YC8 = CONC(NH-1,IA)*HH(IA)*CHAP(RTOP,AMU1,HH(IA))*SIG(IA,NH)
           TTOP = TTOP + YC8
         ENDDO

         TAUD(I) = TTT + TTOP
      ENDDO

      TTOP = 0.
      DO IA = 1, NABS
        YC8 = CONC(NH-1,IA)*HH(IA)*CHAP(RTOP,AMU,HH(IA))*SIG(IA,NH)
        TTOP = TTOP + YC8
      ENDDO

      TAUD(NH) = TTOP

      END
C --------------------------------------------------------------
       SUBROUTINE  TAUDIR_o ( NH, Low_Lev, TAU, SSA, GE, AMU,
     *   ALT, HH, SIG, CONC, NABS, PATH, TAUD) 
C                                 'sigr' and 'air' NOT USED GM1/00
C --------------------------------------------------------------
C    Calculation of the reduced optical path for direct beam
c    from each alt. to space ,   SZA > 90 deg.
c
C  >>>  input:
c
c    NH   - number of levels
c    TAU  - normal optical depths of the layers
c    SSA  - single scattering albedo of the layers
c    GE   - asymm. param. g  --"--
c    AMU  = cos(SZA)
c    ALT  - model altitudes
c    HH   - height of the homogeneous atm. at model top for each absorber
c    SIG  - abs. cross-section
c    CONC - concentrations of the absorbers
c    NABS - number of absorbers
c    PATH - geometric paths through each layer
c --------------------------------------------------------------
c  >>>  output:
c
c    TAUD  - opt. path for direct beam
c ===============================================================
       include 'jpl.par'

       REAL   TAU(NH), SSA(NH), GE(NH), ALT(NH), TAUD(NH), PATH(NH,NH),
     *        CONC(NH,NABS), F(MNH), DZ(MNH),
     *        HH(NABS), SIG(NABS,NH) ! 'air(nh)' NOT USED GM 1/00
       DATA   REAR/6.371E8/
        amu35 = sqrt(1224.*amu*amu+1.)/35.
 101   format(2x, i5, 3f8.2, 4e10.1)
       DO I = 1, NH-1
          F(I) = GE(I)*GE(I)
          DZ(I) = ALT(I+1) - ALT(I)
       ENDDO

       DO I = Low_Lev, NH
          TTT = 0.
          qq = 0.
          DO J = 1, NH-1
             TTT = TTT + TAU(J)*(1.-SSA(J)*F(J))/DZ(J)*PATH(J,I)
          ENDDO
c     Calculate opt. depth above model top by Chapman's function:
C
c                ... first, SZA at model top:
         RL=REAR+ALT(I)
c  from sin theorem:
         RG=RL*SQRT(1.-AMU*AMU)
         RTOP = REAR + ALT(NH)
         VV=RG/RTOP
         AMU1=SQRT(1.-VV*VV)

         TTOP = 0.
         DO IA = 1, NABS
           YC8 = CONC(NH-1,IA)*HH(IA)*CHAP(RTOP,AMU1,HH(IA))*SIG(IA,NH)
           TTOP = TTOP + YC8
         ENDDO

         TAUD(I) = TTT + TTOP
      ENDDO

      END
C --------------------------------------------------------------
      SUBROUTINE  YCOLM ( NH, AMU, ALT, HH, CONC, NABS, PATH, YCOL )
C --------------------------------------------------------------
C    Calculation of  the O2 and O3 opt. paths
c
C  >>>  input:
c
c    NH   - number of levels
c    AMU  = cos(SZA)
c    ALT  - model altitudes
c    HH   - height of the homogeneous atm. at model top for each absorber
c    CONC - concentrations of the absorbers
c    NABS - number of absorbers
c    PATH - geometric paths through each layer
c --------------------------------------------------------------
c  >>>  output:
c
c    YCOL  - O2 and O3 opt. path
c ===============================================================
      include 'jpl.par'
      REAL ALT(NH), PATH(NH,NH),
     *     CONC(NH,NABS), YCOL(NH,NABS), DZ(MNH),
     *     HH(NABS), YC(MNABS)
      DATA   REAR/6.371E8/
      DO I = 1, NH-1
         DZ(I) = ALT(I+1) - ALT(I)
      ENDDO
      RTOP = REAR + ALT(NH)
      DO I = 1, NH-1
         DO IA = 1, NABS
            YC(IA) = 0.
         ENDDO
         DO J = I, NH-1
            DO IA = 1, NABS
               YC(IA) = YC(IA) + CONC(J,IA)*PATH(J,I)
            ENDDO
         ENDDO
C
C     Calculate opt. depth above model top by Chapman's function:
C
C                ... first, SZA at model top:
         RL=REAR+ALT(I)
C        from sin theorem:
C
         RG=RL*SQRT(1.-AMU*AMU)
         VV=RG/RTOP
         AMU1=SQRT(1.-VV*VV)
         DO IA = 1, NABS
            YC8 = CONC(NH,IA)*HH(IA)*CHAP(RTOP,AMU1,HH(IA))
            YCOL(I,IA) = YC(IA) + YC8
         ENDDO
      ENDDO
      DO IA = 1, NABS
         YC8 = CONC(NH,IA)*HH(IA)*CHAP(RTOP,AMU,HH(IA))
         YCOL(NH,IA) = YC8
      ENDDO
      RETURN
      END
C --------------------------------------------------------------
       SUBROUTINE  YCOLM_o ( NH, Low_Lev, AMU, ALT, HH, CONC,
     &        NABS, PATH, YCOL )
C --------------------------------------------------------------
C    Calculation of the O2 and O3 opt. paths for SZA > 90 deg.
c
C  >>>  input:
c
c    NH   - number of levels
c    AMU  = cos(SZA)
c    ALT  - model altitudes
c    HH   - height of the homogeneous atm. at model top for each absorber
c    CONC - concentrations of the absorbers
c    NABS - number of absorbers
c    PATH - geometric paths through each layer
c --------------------------------------------------------------
c  >>>  output:
c
c    YCOL  - O2 and O3 opt. path
c ===============================================================
       include 'jpl.par'

       REAL   ALT(NH), PATH(NH,NH),
     *        CONC(NH,NABS), YCOL(NH,NABS), DZ(MNH),
     *        HH(NABS), YC(MNABS)
       DATA   REAR/6.371E8/
       DO I = 1, NH-1
          DZ(I) = ALT(I+1) - ALT(I)
       ENDDO

       DO I = Low_Lev, NH
          DO IA = 1, NABS
             YC(IA) = 0.
          ENDDO
          DO J = 1, NH-1
             DO IA = 1, NABS
                YC(IA) = YC(IA) + CONC(J,IA)*PATH(J,I)
             ENDDO
          ENDDO
c     Calculate opt. depth above model top by Chapman's function:

c                ... first, SZA at model top:
         RL=REAR+ALT(I)
c  from sin theorem:
         RG=RL*SQRT(1.-AMU*AMU)
         RTOP = REAR + ALT(NH)
         VV=RG/RTOP
         AMU1=SQRT(1.-VV*VV)

         DO IA = 1, NABS
            YC8 = CONC(NH-1,IA)*HH(IA)*CHAP(RTOP,AMU1,HH(IA))
            YCOL(I,IA) = YC(IA) + YC8
         ENDDO

      ENDDO
      END
      SUBROUTINE SYSTRB(X,A,B,C,D,N)
C *** Solution of the system of linear alg. eqs  ***
C ***         "obratnaja progonka"  method
C     A(1)*X(2)+B(1)*X(1)+D(1)=0
C     ------------------
C     A(I)*X(I+1)+B(I)*X(I)+C(I)*X(I-1)+D(I)=0    I=2,..,N-1
C     ------------------
C     B(N)*X(N)+C(N)*X(N-1)+D(N)=0
C     ****
C     X     - output (solution), DIM=N
C     A,B,C - coeff.,     DIM=N, A(N)=C(1)=0
C     D     - free coef., DIM=N
C     G     - int. array, DIM=100
C     N     - system order
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIMENSION X(N),A(N),B(N),C(N),D(N),G(100)
      G1=B(N)
      X(N)=-D(N)/G1
      G(N)=-C(N)/G1
      N1=N-1
      DO 1 II=1,N1
         I=N-II
         G1=B(I)+A(I)*G(I+1)
         X(I)=-(A(I)*X(I+1)+D(I))/G1
         G(I)=-C(I)/G1
  1   CONTINUE
      DO 2 I=2,N
         X(I)=X(I-1)*G(I)+X(I)
  2   CONTINUE
      RETURN
      END

C-------------------------------------------------------------
      SUBROUTINE PHRT ( DW, LIR, Q, SIG, S)
C **  Photorate calc. in specified spectr. int. and alt. :
      DIMENSION DW(LIR),Q(LIR),SIG(LIR)
      S=0.
      DO 1 L=1,LIR
         S=S + Q(L)*SIG(L)*DW(L)
    1 CONTINUE
      RETURN
      END
C--------------------------------------------------------------
 
C-----------------------------------------------------------
      SUBROUTINE NIC1S ( NH, NHU, NABS, YCOL, SF1,SF2,SF3 )
C**** Calculation of the       !  O2  in SRB (SF1)****
C     phot. rates using        !  NO  in the (0-0) and (1-0)(SF2)
C     Nicolet's par-zations:   !  H2O in SRB (SF3)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL  YCOL(NH,NABS), SF1(NHU), SF2(NHU), SF3(NHU)
      DO 3 IIH=1,NHU
         IH  = NHU-IIH+1
         IHY = NH-IIH+1
         O2N = YCOL(IHY,1)
         O3N = YCOL(IHY,2)
C        Calculation of SF1
         IF(O2N .LE. 1.E19) then
           SF1(IH)=1.1E-7*EXP(-1.97E-10*O2N**.522)
         else if (O2N .le. 1.e30) then
           SF1(IH)=1.45E8/O2N**.83
         else
           SF1(ih) = 0.
         endif
C          Calculation of SF2
           IF(O2N .le. 2.E19) then
             SF2(IH)=4.75E-6*EXP(-1.5E-20*O2N)
           else
             TAU=1.E-8*O2N**.38+5.E-19*O3N
             if ( TAU .le. 20.) then
               SF2(IH)=4.5E-6*EXP(-TAU)
             else
               SF2(IH) = 0.
             endif
           endif
C          Calculation of SF3
           TAU = 1.E-7*O2N**.35
           if ( TAU .le. 20.) then
              SF3(IH) = 1.3E-6*EXP(-TAU)
           else
              SF3(IH) = 0.
           endif
    3 CONTINUE
      RETURN
      END
c =============================================================
      FUNCTION SR(WLT)
C *** Rayleigh scattering cross-section:
      WN=1000./WLT
      POK=1.+1.E-8*(6432.8+2949810./(146.-WN**2)+25540./(41.-WN**2))
      SR=1.349E-21*(POK**2-1.)**2*WN**4
      RETURN
      END
c  *** Programs for calculating geometric path from i-th level ***
c  ***        to space in full spherical geometry            ***
C vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
      FUNCTION CHAP(Z,AMU,H)
C ------------------------------------------------------
C *** Chapman's  function  for   calculation of OPTICAL
C *** path for certain gas from alt. Z to space
C *** Approximation by [Smith and Smith, ????]
C ------------------------------------------------------
C     H   - the homogeneous atmosphere height
C     AMU - solar zenith angle
C ------------------------------------------------------
      DATA  REAR/6.371E8/ , PI/3.141592653/
      DATA A/1.0607/,B/0.55644/,C/1.06199/,D/1.72456/,
     *     F/0.56499/,G/0.06652/
      XP=(REAR+Z)/H
      Y=SQRT(0.5*XP)*AMU
      IF (Y.GE.8.) THEN
         CHAP=SQRT(0.5*PI*XP)*F/(G+Y)
      ELSE
         CHAP=SQRT(0.5*PI*XP)*(A+B*Y)/(C+D*Y+Y*Y)
      ENDIF
      RETURN
      END
C ---------------------------------------------------------------
      SUBROUTINE GPATH_o ( NBH, NH, AMU, ALT, PATH )
C  ---------------------------------------------------------------
C ***    Calculation of the set of geometric paths through
c ***    each layer from Z(ILEV) down to ZMIN
C ***    at angles > 90 degrees
C  ---------------------------------------------------------------
c  >>>  input:
c
C   NBH  = ILEV  - number of levels below current level ILEV
c   ALT  - model altitude grid
c   AMU  = cos(180-SZA)
c -------------------------------------------------
c  >>>  output:
c
c   PATH - geom. paths for each layer
c
c --------------------------------------------------
c
      REAL   ALT(NH), PATH(NH)
C ----------------------------------------------------------------
C
      DATA  REAR/6.371E8/

      do i = 1, NBH
         PATH(i) = 0.
      enddo


      RL = REAR + ALT(NBH)
      RG=RL*SQRT(1.-AMU*AMU)

      IF(RL-RG.GT.10.)  THEN
          S=SQRT(RL*RL-RG*RG)
      ELSE
          S=0.
      ENDIF

      RLL=RL
      DO 1 IH=NBH-1, 1, -1
         DZ = ALT(IH+1) - ALT(IH)
         RLL=RLL-DZ
         IF(RLL.LE.RG) GOTO 2
         SS=SQRT(RLL*RLL-RG*RG)
         PATH(IH) = 2.*(S-SS)
         S = SS
  1   CONTINUE
      IH=IH+1
  2   CONTINUE
      PATH(IH) = 2.*S

      if ( NBH .lt. NH ) then
         NAH = NH - NBH + 1
         call  GPATH ( NAH, AMU, ALT(NBH), PATH(NBH) )
      endif

      END
C ---------------------------------------------------------------
      SUBROUTINE GPATH(NAH,AMU,ALT,PATH )
C  ---------------------------------------------------------------
C ***    Calculation of the set of geometric paths through
c ***    each layer from Z(ILEV) to TOA
C ***    at angles less or equal 90 degrees
C  ---------------------------------------------------------------
c  >>>  input:
c
C   NAH  = NH-ILEV  - number of levels above current level ILEV
c   ALT  - model altitude grid
c   AMU  = cos(SZA)
c -------------------------------------------------
c  >>>  output:
c
c   PATH - geom. paths for each layer
c
c --------------------------------------------------
c
      REAL ALT(NAH),PATH(NAH)
C ----------------------------------------------------------------
C
      DATA  REAR/6.371E8/
      RL = REAR + ALT(1)
      RG=RL*SQRT(1.-AMU*AMU)
      IF(RL-RG.GT.10.)  THEN
         S=SQRT(RL*RL-RG*RG)
      ELSE
         S=0.
      ENDIF
      RLL=RL
      DO 1 IH=1, NAH-1
         DZ = ALT(IH+1) - ALT(IH)
         RLL=RLL+DZ
         SS=SQRT(RLL*RLL-RG*RG)
         PATH(IH) = SS-S
         IF(PATH(IH).EQ.0.)THEN
            PRINT *,' SS=',SS,' S=',S
            PAUSE
         ENDIF
         S = SS
  1   CONTINUE
      RETURN
      END


      FUNCTION IMIN ( NH, ALT, AMU )
C -------------------------------------------------
C    Number of the lowest directly illuminated level
c --------------------------------------------------
c
c >>> input:
c
c     NH   - total number of levels
c     ALT  - level altitudes
c     AMU  - cos(180-SZA)
c ----------------------------------------------------
c
      REAL  ALT(NH)
      DATA  REAR/6.371E8/
      SINSZA = SQRT(1.-AMU*AMU)
      ZM = REAR*(1./SINSZA-1.)
      DO I = 1, NH-1
         IF ( ALT(I) .GE. ZM ) GOTO 1
      ENDDO
 1    CONTINUE
      IMIN = I
      END
