c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine o3_snglscatt(p_lvl, t_lvl, z_phys, sza, a0, a1, a2, b,
     1   p0, p2, w, p_o3, sigma, rad_av, kern_av, n, nc, ns, ng, nerf,
     2   itr)
c
c     -- O3_snglscatt/Ozone single scattering approximation/N. Nath
c        /March '01
c
c        Modified:  sza made dependent on channel, & several parameters
c        in calling list (that are computed only in first iteration & 
c        used in successive iterations) are moved into save area/N. Nath
c        /July '01
c
c        Extended:  Scope is extended so rad_av & kern_av have a 
c        trailing dimension corresponding to reflecting surface/bottom
c        of atmosphere/N. Nath/September '01
c
c        Extended2:  Calculations leading to kern_av are optimized
c        /N. Nath/November '01
c
c     -- Function:  This routine computes the back-scattered solar 
c        ultraviolet radiation appropriate to the conditions of the SBUV 
c        experiment.  It incorporates ozone absorption and Rayleigh 
c        scattering in the single scattering approximation.  It also 
c        computes the kernel for the determination of the ozone profile.  
c        Spherical geometry and gravitational correction are fully taken
c        into account.
c
c        The computations are based on the explicit, double-integral 
c        representation of the radiative transfer solution:
c
c              I = cons * Intgrl(t(z)*dtau_sc(z); z=0:inf),
c
c              t(z) = exp(- Intgrl(S(z",z)*dtau(z"); z"=z:inf),
c
c              cons = F * P(theta)/(4*pi),
c              dtau_sc(z) = beta(z)*(1/H)*p(z)*dz,
c
c              dtau(z") = dtau_o3(z") + dtau_sc(z")
c                 = alpha(z")*p_o3(z")*dz" + beta(z")*(1/H)*p(z")*dz"
c
c        where I is the intensity of radiation, t is the transmittance, 
c        F is the solar flux, P is the scattering phase function, 
c        dtau_o3 and dtau_sc are the ozone absorption and Rayleigh 
c        scattering differential optical depth components, alpha and 
c        beta are the ozone absorption coefficient and Rayleigh 
c        scattering coefficient in (atmos-cm)**-1 and atmos**-1, 
c        respectively, p_o3 and p are the ozone partial pressure and 
c        atmospheric pressure, in atmos, respectively, H is the 
c        atmospheric scale height in cm at standard temperature and 
c        zero height, and z and z" are the pressure-scaled heights 
c        corresponding to d(logp) = -dz/H and d(logp") = -dz"/H, 
c        respectively.  S(z",z) (when multiplied by dz") is the total,
c        incoming (slant) and outgoing (vertical) optical path through 
c        layer thickness dz", for radiation starting and ending at z" 
c        and scattered from the level z (<z"); it equals 1 + sec(sza"),
c        where sza", the solar zenith angle at z" is related to sza, the
c        solar zenith angle at z, by sin(sza") = ((R + z_phys)/(R + 
c        z"_phys))*sin(sza).  The coefficients alpha(z) and beta(z) 
c        include the gravity correction factor (1 + z_phys/R)**2, where
c        z_phys is the physical height corresponding to z and R is the 
c        radius of the Earth.
c
c        The differential change in I corresponding to changes in the
c        p_o3 profile is given by:
c
c              del[I] = Intgrl(K(z')*del[p_o3](z')*dz'; z'=0:inf),
c
c              K(z') = - cons*alpha(z')*
c                          Intgrl(S(z',z)*t(z)*dtau_sc(z); z=0:z')
c
c        where K is the kernel.
c
c     -- Implementation:  (1) The atmosphere is discretized into n 
c        levels, in equal increments of dz, with level 1 at the surface
c        and level n at the top.  There are n-1 layers between these 
c        levels, layer j lying between levels j and j+1, and an 
c        additional layer from level n to z=inf.  (2) The contributions 
c        to t(z), and subsequently to I, from layers 1 through n-1, are 
c        approximated by layer-averaged, mid-point physical quantities, 
c        in which temperature dependence of the absorption coefficients
c        and the change in gravity correction across the layer are 
c        properly taken into account.  Instead, for (relatively smaller) 
c        contributions from the topmost layer, extending to infinity, 
c        first it is assumed that p_o3(z) is proportional to p(z)**
c        (1/sigma), with a given sigma (i.e., p_o3 has the scale height
c        of sigma*H), and then changes in temperature and gravity 
c        correction across the layer are ignored.  Consequently, these 
c        later contributions are of analytic form and are easily 
c        integrable.  (3) Because t(z) involves z as an integration 
c        limit and the quantities in this integrand for layers 1 through
c        n-1 would be known only at midpoints, it is first computed at
c        the levels and then averaged over layers for computation of I.  
c        Accordingly, S(z",z), which occurs in the definition of t(z), 
c        is computed at levels z.  For similar reasons, S(z',z), which
c        occurs in K(z'), is computed at levels z'.  (Conceptually, this 
c        approach is simpler than other alternatives.)  (4) Notation-
c        wise, z_phys will be denoted at levels only; all other 
c        atmospheric quantities, except those identified with suffix 
c        'lvl', will generally refer to layers.  (5) sec(sza"), it can 
c        be shown, is exactly equal to sec(sza)/sqrt(1 + xi*(2 - xi)* 
c        (tan(sza))**2), where xi = ((z"_phys-z_phys)/R)/(1 + 
c        z"_phys/R).  Its evaluation in various situations is tied with 
c        manipulations of the basic relation (1 + z_phys/R)**-1 = 1 - 
c        Intgrl((T(z')/Ts)*dz'/R; z'=0:z), where T and Ts are the 
c        temperatures at z' and the standard temperature.  The term 
c        (2 - xi) is ignored for the topmost layer; it may or may not be
c        kept for the other layers.  (6) The kernel K is computed only 
c        for layers 1 through n-1; i.e., it is assumed that del[p_o3] 
c        will be defined only for these layers in the profile 
c        determination process, and a new value of sigma will be 
c        determined, after each iteration, with the values of p_o3 at 
c        uppermost layers, which will then be used in the evaluation of 
c        next I.
c
c     -- Assumptions:  (1) H = k*Ts/(m*g) is computed with mean 
c        molecular weight m = 28.964 amu.  For atmospheres with 
c        significantly different m, this value may need to be changed.  
c        (2) In general, terms of the order (z/R)**2, in relation to 1,
c        are ignored, while terms of the order z/R are carefully 
c        maintained everywhere.
c
c     -- Input & Output Variables:
c
            real        p_lvl(n), t_lvl(n), z_phys(n), sza(nc),
                              !atmospheric pressure (atmos), temperature
                              !  (K), & physical height (cm), and 
                              !  solar zenith angle (degree) (input)
     1                  a0(nc,ns), a1(nc,ns), a2(nc,ns), b(nc,ns),
     2                  p0(nc,ns), p2(nc,ns), w(nc,ns),
                              !o3 absorption coefficients (a0 in 
                              !  (atm-cm)**-1), Rayleigh scattering
                              !  coefficients (atm**-1), phase function
                              !  coefficients, & weights for subchannel
                              !  averaging (input)
     3                  p_o3(n-1), sigma
                              !o3 layer pressure (atm) and sigma (input)

            real        rad_av(nc,n), kern_av(nc,n-1,n)
                              !slit-averaged radiance and kernel
                              !  (trailing dimension: reflecting surface
                              !  /bottom of atmosphere) (output)

            integer     n, nc, ns, ng, nerf, itr
                              !number of levels, number of channels,
                              !  number of subchannels in each channel
					!  (<-> slit), dimension of (xg, wg) (4
                              !  for high accuracy, 2 otherwise), number
                              !  of terms in erfc computation (5 for
                              !  high accuracy, 3 otherwise), and
                              !  iteration number (note: slant factors
                              !  & scattering quantities are computed
                              !  in itr 1 only) (input)
c
c     -- Local Variables:
c
            integer     nl, ncl, nsl, ngl
            parameter   (nl = 81, ncl = 10, nsl = 8, ngl = 4)
c
            real        r, h, ts
            parameter   (r = 6.378e8, h = 7.996e5, ts = 273.16)
                              !radius of Earth, scale height at (z=0 & 
                              !  ts), standard temperature

            integer     nchannel  !nc*ns

            real        sec0(ncl), tan20(ncl), delz, phase(ncl,nsl)
                              !constants

            real        alpha(ncl*nsl,nl), beta(ncl*nsl,nl),
                              !effective absorption & scattering coeffs
                              !  incorporating gravity correction
                              !  factors
     1                  st_lvl(ncl,nl,nl), s_lvl(ncl,nl,nl-1),
     2                  arg(ncl,nl), shift(ncl,nl-1),
                              !slant factors (st_lvl & s_lvl), & related
                              !  quantities
     3                  deltau_sc(ncl*nsl,nl)
                              !dtau_sc

            real        xg(ngl), wg(ngl)
                              !Gaussian positions & weights (0,1) for
                              !  contribution to radiance from the 
                              !  topmost layer

            save        nchannel, sec0, tan20, delz, phase, alpha, beta,
     1                  st_lvl, s_lvl, arg, shift, deltau_sc, xg, wg
                              !save these quantities for itr > 1

            real        fac(nl), p(nl-1), t(nl-1), pt,
                              !gravity correction factor, layer-averaged 
                              !  atmospheric pressures & temperatures,
                              !  pressure at top 
     1                  xi0(nl,nl), xit_lvl(nl,nl-1), xi_lvl(nl,nl-1), 
     2                  xit_lvln,
                              !various xi's used in defining s's
     3                  arg_o3(ncl,nl), sto3_lvl(ncl,nl), pt_o3,
                              !quantities similar to arg, slant factor,
                              !  & top pressure, for o3
     4                  deltau_o3(ncl*nsl,nl), deltau(ncl*nsl,nl), 
     5                  tr_lvl(ncl*nsl,nl), tr(ncl*nsl,nl-1),
                              !dtau_o3, dtau, level and layer 
                              !  transmittances
     6                  rad(ncl*nsl,nl), kern(ncl*nsl,nl-1,nl)
                              !radiances & kernel (pre-average)

            real        delzr, delzh, delt, temp, sqrtsigma, 
     1                  del, eps(nl), c1, c2, sum(nl,nl)
                              !various work variables

            real        pi, sqrtpi
            parameter   (pi = 3.141592654, sqrtpi = 1.772453851)

            integer     i, j, k, l, ic, is, ig, i0
c
            real        xg4(4), wg4(4)
            data        xg4	/0.06943184, 0.33000948, 0.66999052, 
     1                         0.93056816/
            data        wg4   /0.17392742, 0.32607258, 0.32607258,
     1                         0.17392742/
                              !4-point Gaussian
            real        xg2(2), wg2(2)
            data        xg2   /0.21132487, 0.78867513/
            data        wg2   /0.5, 0.5/
                              !2-point Gaussian
c
            real        f, q3, f3, q5, f5, x
            parameter   (q3 = 0.47047, q5 = 0.3275911)
            f3(x) = ((0.7478556*x - 0.0958798)*x + 0.3480242)*x
                              !Chapman function, 3-term: (f3(x), q3)
            f5(x) = ((((1.061405429*x - 1.453152027)*x + 1.421413741)*x
     1         - 0.284496736)*x + 0.254829592)*x
                              !5-term: (f5(x), q5)
c
c           Bypass calculations if itr > 1
c
      if (itr .gt. 1) go to 100
c
c           Set up invariant quantities
c
c           o  Define constants
c
      nchannel = nc*ns        !total number of subchannels
      do ic = 1,nc
         sec0(ic) = 1.0/cos(min(max(sza(ic), 1e-3), 89.999)*pi/180.0)
                              !ensures proper behavior of algorithm at 
                              !  sza very close to 0 and 90 degrees
         tan20(ic) = sec0(ic)**2 - 1.0
      end do
      delz = h*log(p_lvl(1)/p_lvl(n))/(n - 1)
                              !layer thickness
      do ic = 1,nc
         do is = 1,ns
            phase(ic,is) = w(ic,is)*(p0(ic,is) + p2(ic,is)/sec0(ic)**2)
                              !weight times phase function
         end do
      end do
c
c           o  Define auxiliary quantities
c
      do k = 1,n
         fac(k) = 1.0 + z_phys(k)/r		
      end do
      do k = 1,n-1
         p(k) = sqrt(p_lvl(k)*p_lvl(k+1))
         t(k) = (t_lvl(k) + t_lvl(k+1))/2.0
      end do
      delzr = delz/r
      delzh = delz/h
      pt = p_lvl(n)           !p top
c
c           o  Define Gaussian positions & weights
c
      if (ng .eq. 4) then
            do ig = 1,4
               xg(ig) = xg4(ig)
               wg(ig) = wg4(ig)
            end do
         else if (ng .eq. 2) then
            do ig = 1,2
               xg(ig) = xg2(ig)
               wg(ig) = wg2(ig)
            end do
         else
            write (*, *) "Input error:  ng must equal 4 or 2"
      end if        
c
c           o  Define alpha and beta for all channels and components
c
      do ic = 1,nc            !channels with distinct wavelengths
         do is = 1,ns         !subchannels within each channel
            i = (ic-1)*ns +is !channels & subchannels in a linear chain
            do k = 1,n-1
               delt = t(k) - ts
               alpha(i,k) = ((a2(ic,is)*delt + a1(ic,is))*delt + 
     1           a0(ic,is))* fac(k)*fac(k+1)
               beta(i,k) = b(ic,is)*fac(k)*fac(k+1)
            end do
            delt = t_lvl(n) - ts
            alpha(i,n) = ((a2(ic,is)*delt + a1(ic,is))*delt + 
     1        a0(ic,is))*fac(n)**2
            beta(i,n) = b(ic,is)*fac(n)**2
         end do
      end do
c
c           o  Define xit_lvl ('t' for transpose) & xi_lvl, and compute
c              st_lvl & s_lvl, for levels 1 through n
c                 !note: the various xi's are only defined in parts of
c                 !  index spaces, where needed
c
      do k = 1,n              !both k & j in xi0 are levels
         xi0(k,k) = 0.0
         temp = 0.0
         do j = k+1,n
            temp = temp + t(j-1)/ts
            xi0(j,k) = temp*fac(k)*delzr
         end do
      end do
c
      do l = 1,n-1            !k in xit is a level, l a layer;
                              !  l <-> z" integration in tr
         do k = 1,l
            xit_lvl(k,l) = (xi0(l,k) + xi0(l+1,k))/2.0
         end do
      end do

      do k = 1,n-1            !j in xi is a level, k a layer;
                              !  k <-> z integration in kern
        do j = k+1,n
            xi_lvl(j,k) = (xi0(j,k) + xi0(j,k+1))/2.0
        end do
      end do
c
      do ic = 1,nc
         do l = 1,n-1         !k is a level, l is a layer
            do k = 1,l
               st_lvl(ic,k,l) = 1.0 + sec0(ic)/
     1           sqrt(1.0 + xit_lvl(k,l)*(2.0 - xit_lvl(k,l))*tan20(ic))
            end do
         end do
      end do

      do ic = 1,nc
         do k = 1,n-1         !j is a level, k is a layer
            do j = k+1,n
               s_lvl(ic,j,k) = 1.0 + sec0(ic)/
     1           sqrt(1.0 + xi_lvl(j,k)*(2.0 - xi_lvl(j,k))*tan20(ic))
            end do
         end do
      end do
c
c           o  Define equivalents of st_lvl for the topmost layer
c              (l = n)
c                 !note: these are equivalent, point quantities that go
c                 !  with deltau_sc for layer n (as given later), and 
c                 !  are results of integration over this layer, with
c                 !  xit_lvl for each level expressed in terms of 
c                 !  xit_lvl at level n; the integral involved is:
c                 !  Intgrl((exp(-x)/sqrt(1 + 2*a*x))*dx; x=0:inf), 
c                 !  which equals sqrt(pi)*arg*exp(arg**2)*erfc(arg),
c                 !  where arg = 1/sqrt(2*a) and the last two terms are 
c                 !  approximated by f(1.0/(1.0 + q*arg))
c
      xit_lvln = fac(n)*(t_lvl(n)/ts)*h/r
      do ic = 1,nc
         arg(ic,n) = 1.0/sqrt(2.0*xit_lvln*tan20(ic))
         if (nerf .eq. 5) then
               f = f5(1.0/(1.0 + q5*arg(ic,n)))
            else
               f = f3(1.0/(1.0 + q3*arg(ic,n)))
         end if
         st_lvl(ic,n,n) = 1.0 + sec0(ic)*sqrtpi*arg(ic,n)*f
      end do

      eps(n) = 0.0
      do k = n-1,1,-1
         del = fac(k)*(t(k)/ts)*delzr
         eps(k) = eps(k+1) + del
      end do
      do ic = 1,nc
         do k = n-1,1,-1
            shift(ic,k) = 1.0 + 2.0*eps(k)*tan20(ic)
            arg(ic,k) = 1.0/sqrt(2.0*((1.0 - eps(k))/shift(ic,k))*
     1        xit_lvln*tan20(ic))
                  !note: xit_lvlk = xit_lvln*(1.0 - eps(k)) + eps(k)
            if (nerf .eq. 5) then
                  f = f5(1.0/(1.0 + q5*arg(ic,k)))
               else
                  f = f3(1.0/(1.0 + q3*arg(ic,k)))
            end if
            st_lvl(ic,k,n) = 1.0 + sec0(ic)*(1.0/sqrt(shift(ic,k)))*
     1        sqrtpi*arg(ic,k)*f
         end do
      end do
c
c           o  Define deltau_sc and save for next iterations
c
      do i = 1,nchannel
         do l = 1,n-1
            deltau_sc(i,l) = beta(i,l)*p(l)*delzh
         end do
         deltau_sc(i,n) = beta(i,n)*pt
                              !layer n (equivalent thickness ~ h)
      end do
c
  100 continue
c
c           Compute tr, first at levels starting from the top, and then
c           at mid-points of layers
c
c           o  Define equivalents of st_lvl for o3 for the topmost layer
c              (l = n); other layers share the same st with 'sc'
c                             !index l = n is implied here
      sqrtsigma = sqrt(sigma)

      do ic = 1,nc
         arg_o3(ic,n) = arg(ic,n)/sqrtsigma
         if (nerf .eq. 5) then
               f = f5(1.0/(1.0 + q5*arg_o3(ic,n)))
            else
               f = f3(1.0/(1.0 + q3*arg_o3(ic,n)))
         end if
         sto3_lvl(ic,n) = 1.0 + sec0(ic)*sqrtpi*arg_o3(ic,n)*f
      end do

      do ic = 1,nc
         do k = n-1,1,-1
            arg_o3(ic,k) = arg(ic,k)/sqrtsigma
            if (nerf .eq. 5) then
                  f = f5(1.0/(1.0 + q5*arg_o3(ic,k)))
               else
                  f = f3(1.0/(1.0 + q3*arg_o3(ic,k)))
            end if
            sto3_lvl(ic,k) = 1.0 + sec0(ic)*(1.0/sqrt(shift(ic,k)))*
     1        sqrtpi*arg_o3(ic,k)*f
         end do
      end do
c
c           o  Define deltau_o3; add to deltau_sc to form deltau
c
      pt_o3 = p_o3(n-1)*exp(-delz/(h*sigma)/2.0)
                              !p_o3 top
      do i = 1,nchannel
         do l = 1,n-1
            deltau_o3(i,l) = alpha(i,l)*p_o3(l)*delz
            deltau(i,l) = deltau_o3(i,l) + deltau_sc(i,l)
         end do
         deltau_o3(i,n) = alpha(i,n)*pt_o3*sigma*h
                              !layer n (equivalent thickness ~ sigma*h)
         deltau(i,n) = deltau_o3(i,n) + deltau_sc(i,n)
      end do
c
c           o  Compute tr
c
      do i = 1,nchannel
         ic = (i-1)/ns + 1
         do k = n,1,-1        !this k is a level
            temp = sto3_lvl(ic,k)*deltau_o3(i,n) + 
     1               st_lvl(ic,k,n)*deltau_sc(i,n)
            do l = n-1,k,-1
               temp = temp + st_lvl(ic,k,l)*deltau(i,l)
            end do
            tr_lvl(i,k) = exp(-temp)
         end do
      end do

      do i = 1,nchannel
         do k = n-1,1,-1      !this k is a layer
            tr(i,k) = sqrt(tr_lvl(i,k)* tr_lvl(i,k+1))
         end do
      end do
c
c           Compute rad, the intensities for the various channels
c
      do i = 1,nchannel
c
c           o  Compute the contribution from the topmost layer
c                 !note: tr for layer n, as a function of z, is given by
c                 !  exp(-(c1*(p(z)/pt)**(1.0/sigma) + c2*(p(z)/pt))),
c                 !  and the integral runs over (p(z)/pt=0:1)
c
         ic = (i-1)/ns + 1
         c1 = sto3_lvl(ic,n)*deltau_o3(i,n)
         c2 = st_lvl(ic,n,n)*deltau_sc(i,n)
         temp = 0.0
         do ig = 1,ng
            temp = temp + wg(ig)*exp(-(c1*xg(ig)**(1.0/sigma) + 
     1        c2*xg(ig)))
         end do
         rad(i,n) = temp*deltau_sc(i,n)
c
c           o  Add contributions from all other layers; multiply the 
c              result by flux & phase function
c
         do k = n-1,1,-1
            rad(i,k) = rad(i,k+1) + tr(i,k)*deltau_sc(i,k)
         end do
                              !flux & phase incorporated while averaging

      end do
c           
c           Average subchannel rad's
c
      do ic = 1,nc
         i0 = (ic-1)*ns
         do k = 1,n
            temp = 0.0
            do is = 1,ns
               temp = temp + phase(ic,is)*rad(i0+is,k)
            end do
            rad_av(ic,k) = temp/(4.0*pi)
                              !flux assumed to be 1
         end do
      end do
c
c           Compute kern, the kernel
c
      do i = 1,nchannel

         ic = (i-1)/ns + 1
         do j = 1,n           !j is a level
            sum(j,j) = 0.0    !starting point of defining kernel: 
                              !  sum(j,k) = 0 for k = j; k, a level,
                              !  stands for reflecting surface/bottom of
                              !  atmosphere; sum(j,k) for k > j is also
                              !  0, but not needed
            do k = j-1,1,-1
               sum(j,k) = sum(j,k+1) + 
     1                      s_lvl(ic,j,k)*tr(i,k)*deltau_sc(i,k)
            end do
         end do

         do k = 1,n
            do j = k,n-1      !j here is a layer
               kern(i,j,k) = - alpha(i,j)*((sum(j,k) + sum(j+1,k))/2.0)
     1                           *delz
                              !this defines the non-zero elements of
                              !  kern; the additional delz arises from 
                              !  discretization
            end do
         end do
                              !flux & phase incorporated while averaging
      end do
c
c           Average subchannel kern's
c
      do ic = 1,nc
         i0 = (ic-1)*ns
         do k = 1,n
            do j = 1,k-1
               kern_av(ic,j,k) = 0.0
                              !no contribution to the kernel from layers
                              !  below the reflecting surface 
            end do
            do j = k,n-1
               temp = 0.0
               do is = 1,ns
                  temp = temp + phase(ic,is)*kern(i0+is,j,k)
               end do
               kern_av(ic,j,k) = temp/(4.0*pi)
            end do
         end do
      end do
c
      return
      end