c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine o3_retrieval(p_lvl, t_lvl, sza, a0, a1, a2, b, p0, p2,
     1   w, nval_obs, q0, q_a, cov_a, cov_e, rad_m0, kern_m0, p_cld,
     2   p_terr, cldfrac, adjust_r, kern_mr, q, rmix, nval_s, nval,
     3   kern_s0, kern0, kern_s, kern, avkern, cov, sigma0, sigma,
     4   qtot0, qtot, res0, res, dr, tol_res, tol_q, n, nc, ns, nitr,
     5   itrc, ssonly, pcode)
c
c     -- ALERT:  nc must equal ncl; if nc other than 10 is desired, ncl
c        should also be changed to that value.
c
c     -- O3_retrieval/Ozone retrieval using Rodgers method/N. Nath
c        /May '01
c
c        Modified to conform with changes in o3_snglscatt (channel 
c        dependent sza & localized variables)/N. Nath/July '01
c
c        dyn: Dynamic allocation for rad_s, kern_sp, kern_x, cov_x, 
c        okern_x, pindex w.r.t. nc/August '01
c
c        dynmod0: Definitions of nval_m0 and kern_m0 are changed such 
c        that they are no longer additive with their single scattering 
c        counterparts (only the radiances are additive); also dynamical 
c        allocations for new local variables rad_m, rad, & f
c        /September '01
c
c        dynmod: Input changed from nval_m0 to rad_m0, & nval_m is not
c        computed & consequently not output; this is done to avoid
c        taking logarithm of very small nval_m0, which for computational
c        error may become negative/September '01
c
c        dynext: Scope extended to include (channel-dependent) cloud top
c        & terrain pressures (p_cld, p_terr), which serve as reflecting
c        surfaces and provide lower limits of atmosphere; single 
c        scattering radiance & kernel, thereby, take on a trailing
c        dimension corresponding to the reflecting level/September '01
c
c        Added mixing ratio to output list/March '02
c
c        Took out all dynamic allocations; also, added reflectivity 
c        adjustments due to re-determination of total ozone via 
c        determination of profile/April '02
c
c     -- Function:  This routine combines single- and multiple-
c        scattering contributions to n-values and weighting functions,
c        along with the Rodgers algorithm, to produce the ozone vertical
c        profile that best fits the SBUV ozone data.  The details of the 
c        various quantities are to be found in the lower-level routines.
c        The retrieval iteration is carried out until a desired 
c        convergence in the O3 profile and/or in the measured data is
c        achieved.  Following this, the averaging kernel and the 
c        covariance of the final solution are computed.  The ozone 
c        profile and all altitude-related quantities are arranged from 
c        bottom to top of the atmosphere, and the channels are arranged
c        from short to long wavelengths.
c
c     -- Notes:  (1) Ozone in the top layers is assumed to be
c        distributed according to the power law:  q proportional to 
c        p**(1/sigma), with sigma determined from nsig layers below the
c        topmost layer (n).  The amount in layer n (extending to 
c        infinity) is directly related with this sigma and the amount 
c        in layer n-1.  The kernel, in conformity with this idea, only
c        extends up to layer n-1.  (2) nval_m0 and kern_m0 are assumed 
c        to be zero for the first four channels, but places are kept 
c        for them in case something changes.  (3) Optical coefficients
c        are not computed in this routine as they only depend on the 
c        channel wavelengths and not on the individual profile.
c
c     -- Input & Output Variables:
c
            real        p_lvl(n), t_lvl(n), sza(nc),
                              !atmospheric pressure (atmos), temperature
                              !  (K), 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                  nval_obs(nc),
                              !observed n-value (-100*log_10(radiance/
                              !  flux)) (input)
     4                  q0(n), q_a(n), cov_a(n-1,n-1), cov_e(nc,nc),
                              !first guess & a-priori profile of layer
                              !  ozone amount (DU) (w/ nth value
                              !  containing the amount above level n), 
                              !  a-priori covariance matrix & 
                              !  measurement-error covariance matrix
                              !  (input)
     5                  rad_m0(nc), kern_m0(nc,n-1)
                              !multiple scattering radiance & kernel
                              !  computed at q0 (not needed if ssonly is
                              !  true) (input)

         real           p_cld, p_terr, cldfrac(nc)
                              !cloud & terrain pressure (atmos), cloud
                              !  fraction (input)

         real           adjust_r, kern_mr(nc)
                              !change in reflectivity per unit change
                              !  in total ozone, change in multiple
                              !  scattering n-value per unit change in 
                              !  reflectivity (input)

         real           q(n), rmix(n-1),
                              !retrieved layer ozone profile (DU) &
                              !  ozone volume mixing ratio (ppmv)
                              !  (output)
     1                  nval_s(nc), nval(nc), 
                              !retrieved single-scattering & total
                              !  n-values (output)
     2                  kern_s0(nc,n-1), kern0(nc,n-1),
                              !initial values of single-scattering &
                              !  total kernels
     3                  kern_s(nc,n-1), kern(nc,n-1),
                              !retrieved single-scattering & total
                              !  kernels (multiple-scattering kernel 
                              !  fixed at kern_m0) (output)
     4                  avkern(n-1,n-1), cov(n-1,n-1),
                              !averaging kernel & covariance of the
                              !  final solution (output)
     5                  sigma0, sigma, qtot0, qtot, res0(nc), res(nc)
                              !initial and final values of: ratio of 
                              !  ozone scale height to pressure scale
                              !  scale height, total ozone (DU), &
                              !  n-value residues (output)

            real        dr    !change in reflectivity due to change in 
                              !  total ozone (output)

            real        tol_res, tol_q
                              !error tolerances in residue & fractional 
                              !  change in q from previous iteration
                              !  (input)

            integer     n, nc, ns, nitr
                              !number of levels, number of channels,
                              !  number of subchannels in each channel
					!  (<-> slit), maximum number of
                              !  iterations (input)

            integer     itrc  !number of iterations required for 
                              !  convergence (output)

            logical     ssonly!if 'true', only single scattering (input)

            integer     pcode !print code: 0 - all, 1 - summary only
c
c     -- Other Variables:
c
            integer     nl, ncl
            parameter   (nl = 81, ncl = 10)

            integer     ng, nerf
            parameter   (ng = 2, nerf = 3)

            integer     nsig 
            parameter   (nsig = 8)
                              !number of layers used in determining
                              !  sigma
c
            real        re, h
            parameter   (re = 6.378e8, h = 7.996e5)
                              !radius of Earth and pressure scale height
                              !  at the surface at standard temperature

            real        cons
            parameter   (cons = -43.429448)
                              !-100/log(10.0)
c
            real        z_phys(nl), p_o3(nl-1), fac(nl-1)
                              !altitude corresponding to p_lvl (cm), 
                              !  ozone layer pressure corresponding to q
                              !  (atm), & q to p_o3 ratio

            real        rad_e(ncl,nl), kern_ep(ncl,nl-1,nl)
                              !computed, extended single-scattering
                              !  radiance over flux (I/F), & single-
                              !  scattering kernel (d/dp_o3 I/F)
                              !  (output of 'o3_snglscatt')

            real        rad_sc(ncl), kern_spc(ncl,nl-1),
     1                  rad_st(ncl), kern_spt(ncl,nl-1),
     2                  rad_s(ncl), kern_sp(ncl,nl-1)
                              !cloud, terrain, & combined single-
                              !  scattering radiance & kernel

            real        z_cld(ncl), z_terr(ncl), fc, ft
                              !cloud & terrain height (cm), corres-
                              !  ponding fractions

            integer     lc, lt!location of cloud & terrain

            real        dnval_m(ncl), rad_m(ncl), rad(ncl), f(ncl)
                              !change in multiple-scattering nvalue,
                              !  multiple-scattering radiance, total 
                              !  radiance, & the ratio rad_s/rad

            real        q_in(nl-1)
                              !input q into the update program

            real        kern_x(ncl,nl-1 ), cov_x(ncl,ncl), 
     1                  okern_x(ncl,nl-1)
                              !kern*cov_a, LU-factored form of
                              !  o = kern*cov_a*kern_t + cov_e, &
                              !  o**-1*kern_x

            real        z(nl), delz
                              !pressure-scaled height (cm), layer 
                              !  thickness (cm) 

            real        logq(nsig)
                              !log of q0 & q for the top nsig layers  

            integer     pindex(ncl)
                              !row permutn index of kern*cov_a*kern_t

            real        slope, intcpt
                              !slope & intercept of linear fit to logq

            integer     itr, ic, j, k, l, isig

            real        resmax, resmean, diffqmax, temp
                              !residue maximum, residue st. devn., 
                              !  maximum fractional difference of q from 
                              !  previous q, & a temporary variable

            integer     icmax, jmax
                              !channel number at resmax, layer number at 
                              !  diffqmax         

            real        dbl   !double/single precision number
            real*4      sngl  !single precision number
c
      sngl(dbl) = dbl
c
      if (nc .ne. ncl) then
         write (*, *) "Error in o3_retrieval:  ncl must equal nc"
         write (*, *) "nc, ncl: ", nc, ncl
      end if
c
c           Initialize before entering the iterative loop:  define z;
c           compute sigma0 & qtot0; set q = q0 & sigma = sigma0
c
      delz = h*log(p_lvl(1)/p_lvl(n))/(n - 1)
      do j = 1,n
         z(j) = (j - 1)*delz
      end do

      do isig = 1,nsig
         logq(isig) = log(q0(n-nsig-1+isig))
      end do
      call linear_fit(z(n-nsig), logq, slope, intcpt, nsig)
      sigma0 = -1.0/(h*slope)

      qtot0 = q0(n)           !assumed q0(n) contains total amount of
                              !  ozone above level n
      do j = n-1,1,-1
         qtot0 = qtot0 + q0(j)
      end do

      do j = 1,n-1
         q(j) = q0(j)
      end do
      sigma = sigma0
c
c           Define z_phys & fac
c
      call altitude(p_lvl, t_lvl, z_phys, n)

      do j = 1,n-1
         fac(j) = ((1 + z_phys(j)/re)*(1 + z_phys(j+1)/re)*delz)*1e3
                              !ratio of q over p_o3: q in DU, p_o3 in 
                              !  atm
      end do
c
c           Solve the retrieval problem by iteration
c
      do itr = 1,nitr

         itrc = itr           !number of iterations required for 
                              !  convergence
c
c              o  Determine single-scattering n-value and kernel
c
c                 -  Convert q into p_o3
c
         do j = 1,n-1
            p_o3(j) = q(j)/fac(j)
         end do
c
c                 - Call single-scattering forward model
c
         call o3_snglscatt(p_lvl, t_lvl, z_phys, sza, a0, a1, a2, b,
     1     p0, p2, w, p_o3, sigma, rad_e, kern_ep, n, nc, ns, ng, nerf,
     2     itr)
c
c                 - Compute radiance & kernel appropriate to cloud
c                   & terrain pressure
c
         do ic = 1,nc
            z_cld(ic) = h*log(p_lvl(1)/p_cld)
            lc = int(z_cld(ic)/delz) + 1
                              !cloud lies between levels lc & lc+1
            fc = (z_cld(ic) - z(lc))/delz
            rad_sc(ic) = fc*rad_e(ic,lc+1) + (1.0 - fc)*rad_e(ic,lc)
            do j = 1,n-1
               kern_spc(ic,j) = fc*kern_ep(ic,j,lc+1) + 
     1                            (1.0 - fc)*kern_ep(ic,j,lc)
                              !weighted average of lc & lc+1
            end do
            z_terr(ic) = h*log(p_lvl(1)/p_terr)
            lt = int(z_terr(ic)/delz) + 1
                              !terrain lies between levels lt & lt+1
            ft = (z_terr(ic) - z(lt))/delz
            rad_st(ic) = ft*rad_e(ic,lt+1) + (1.0 - ft)*rad_e(ic,lt)
            do j = 1,n-1
               kern_spt(ic,j) = ft*kern_ep(ic,j,lt+1) + 
     1                            (1.0 - ft)*kern_ep(ic,j,lt)
            end do
         end do
c
c                 - Combine radiance & kernel from cloud & terrain 
c                   levels
c
         do ic = 1,nc
            rad_s(ic) = cldfrac(ic)*rad_sc(ic) + 
     1                    (1.0 - cldfrac(ic))*rad_st(ic)
            do j = 1,n-1
               kern_sp(ic,j) = cldfrac(ic)*kern_spc(ic,j) + 
     1                           (1.0 - cldfrac(ic))*kern_spt(ic,j)
            end do
         end do
c
c                 - Convert radiance & kernel into n-value & 
c                   d(n-value)/dq
c
         do ic = 1,nc
            nval_s(ic) = cons*log(rad_s(ic))
                              !cons = -100/log(10)
            do j = 1,n-1
               kern_s(ic,j) = (cons/rad_s(ic))*kern_sp(ic,j)/fac(j)
                              !fac(j) is needed because kern_sp is
                              !  w.r.t. p
               if (itr .eq. 1) kern_s0(ic,j) = kern_s(ic,j)
            end do
         end do
c
c              o  Add multiple-scattering contributions to n-value and 
c                 kernel (note:  both nval_m0 and kern_m0 are zero for
c                 first 4 channels, wavelength < 290 nm)
c
         do ic = 1,4
            nval(ic) = nval_s(ic)
         end do
         do ic = 5,nc
            if (ssonly) then
                  nval(ic) = nval_s(ic)
               else
                  dnval_m(ic) = 0.0
                  do j = 1,n-1
                     dnval_m(ic) = dnval_m(ic) + kern_m0(ic,j)*
     1                 (q(j) - q0(j))
                              !change in nval_m; assumed kern_m is fixed
                              !  at kern_m0, computed at q0 
                  end do
                  if (itr .gt. 1) then
                     dr = adjust_r*(qtot - qtot0)
                              !total change in reflectivity,
                              !  such that the change in mult. scatt. 
                              !  n-value in 331 nm channel due to
                              !  the change in qtot is zero
                     dnval_m(ic) = dnval_m(ic) + kern_mr(ic)*dr
                              !resulting change in dnval_m
                  end if
                  rad_m(ic) = rad_m0(ic)*exp(dnval_m(ic)/cons)
                  rad(ic) = rad_s(ic) + rad_m(ic)
                  nval(ic) = cons*log(rad(ic))
                  f(ic) = rad_s(ic)/rad(ic)
            endif
         end do
c
         do ic = 1,4
            do j = 1,n-1
               kern(ic,j) = kern_s(ic,j)
               if (itr .eq. 1) kern0(ic,j) = kern(ic,j)
            end do
         end do
         do ic = 5,nc
            do j = 1,n-1
               if (ssonly) then
                     kern(ic,j) = kern_s(ic,j)
                  else
                     kern(ic,j) = f(ic)*kern_s(ic,j) + 
     1                              (1.0 - f(ic))*kern_m0(ic,j)
               end if
               if (itr .eq. 1) kern0(ic,j) = kern(ic,j)
            end do
         end do
c
c              o  Compute residue & exit loop if desired limit is 
c                 attained
c
         do ic = 1,nc
            res(ic) = nval_obs(ic) - nval(ic)
            if (itr .eq. 1) res0(ic) = res(ic)
                              !save for output
         end do
         resmax = 0.0
         do ic = 1,nc
            temp = abs(res(ic))
            if (temp .gt. resmax) then
              resmax = temp
              icmax = ic
            end if
         end do
         resmean = 0.0
         do ic = 1,nc
            resmean = resmean + res(ic)**2
         end do
         resmean = sqrt(resmean/nc)
c <<
         if (pcode .eq. 0) then
c            write (7, *) "nval_fit: ",(sngl(nval(ic)), ic=1,nc)
c            write (7, *) "res: ", (sngl(res(ic)), ic=1,nc)
            write (7, *) "resmax, icmax: ", sngl(resmax), icmax
            write (7, *) "resmean: ", sngl(resmean)
         end if
c >>
         if ((resmax .lt. tol_res) .and. (itr .gt. 1)) go to 100
c
c              o  Update q (Rodgers method)
c
         do j = 1,n-1
            q_in(j) = q(j)
         end do
         call update(q_in, nval, kern, nval_obs, q_a, cov_a, cov_e, q,
     1     kern_x, cov_x, n-1, nc, pindex, .true.)
c
c              o  Compute sigma and qtot
c
c ..              -  Quit iteration if q in any top layer is negative

         do isig = 1,nsig
            if (q(n-nsig-1+isig) .le. 0.0) then
               q(n) = q_in(n) !do not determine sigma, use q_in(n) for
                              !  q(n)
               qtot = q(n)    !determine qtot
               do j = n-1,1,-1
                  qtot = qtot + q(j)
               end do
               write (*, *)
               write (*, *)
     1           "negative ozone at top layers: itr stopped at ", itr
               write (*, *) "q: ", (sngl(q(j)), j=1,n)
               go to 100
            end if
         end do
c ..
         do isig = 1,nsig
            logq(isig) = log(q(n-nsig-1+isig))
         end do
         call linear_fit(z(n-nsig), logq, slope, intcpt, nsig)
         sigma = -1.0/(h*slope)

         do isig = nsig/2+1,nsig
            q(n-nsig-1+isig) = exp(slope*z(n-nsig-1+isig) + intcpt)
                              !replace q(n-nsig/2) to q(n-1) by fitted
                              !  values; this is done for smooth 
                              !  transition of the profile, and it is 
                              !  expected that the change is very small
         end do
         q(n) = q(n-1)*(q(n-1)/q(n-2))/(1.0 - q(n-1)/q(n-2))
                              !geometric series with ratio r = q(n-1)/
                              !  q(n-2); the same value is obtained by
                              !  using the formula (d/dz qtot at 
                              !  level n)*h*sigma (in fact this equality 
                              !  prescribes the value of the indicated 
                              !  derivative)
         qtot = q(n)
         do j = n-1,1,-1
            qtot = qtot + q(j)
         end do
c
c              o  Exit loop if desired accuracy in q is attained
c
         diffqmax = 0.0
         do j = 1,n-1
            temp = abs(1.0 - q(j)/q_in(j))
            if (temp .gt. diffqmax) then
               diffqmax = temp
               jmax = j
            end if
         end do
c <<
         if (pcode .eq. 0) then
            if (itr .eq. 1) then
              write (7, *)
              write (7, *) "begin iteration"
            end if
            write (7, *)
            write (7, *) "iteration: ", itr
            write (7, *)
c            write (7, *) "q: ", (sngl(q(j)), j=1,n)
            write (7, *) "sigma, qtot: ", sngl(sigma), sngl(qtot)
            write (7, *) "diffqmax, jmax: ", sngl(diffqmax), jmax
         end if
c >>
      
         if (diffqmax .lt. tol_q) go to 100

      end do

  100 continue

      if (pcode .eq. 1) then
         write (7, *)
         write (7, *) "itrc: ", itrc
         write (7, *)
         write (7, *) "nval_fit: ",(sngl(nval(ic)), ic=1,nc)
         write (7, *) "res: ", (sngl(res(ic)), ic=1,nc)
         write (7, *) "resmax, icmax: ", sngl(resmax), icmax
         write (7, *) "sigma, qtot: ", sngl(sigma), qtot
         write (7, *) "diffqmax, jmax: ", sngl(diffqmax), jmax
      end if
c
c           Compute ozone mixing ratio (in ppmv)
c
      do j = 1,n-1
         p_o3(j) = q(j)/fac(j)
         rmix(j) = ( p_o3(j)/sqrt(p_lvl(j)*p_lvl(j+1)) )*1e6
      end do
c
c           Compute averaging kernel & covariance of the solution
c
c              o  First compute (kern*cov_a*kern_t + cov_e)**-1*kern_x
c
      do j = 1,n-1
         do k = 1,nc
            okern_x(k,j) = kern_x(k,j)
         end do
      end do
      call inverse(cov_x, okern_x, pindex, nc, n-1, .true.)
                              !remember okern_x is overwritten with 
                              !  the desired result
c
c              o  Left-multiply kern with transpose of okern_x 
c
      do l = 1,n-1
         do j = 1,n-1         !averaging kernel is not symmetric
            avkern(j,l) = 0.0
            do k = 1,nc
               avkern(j,l) = avkern(j,l) + okern_x(k,j)*kern(k,l)
            end do
         end do
      end do

      do l = 1,n-1
         do j = 1,l           !solution covariance is symmetric
            cov(j,l) = cov_a(j,l)
            do k = 1,nc
               cov(j,l) = cov(j,l) - okern_x(k,j)*kern_x(k,l)
            end do
            cov(l,j) = cov(j,l)
         end do
      end do
c
c
      return
      end

