c2--567--1---------2---------3---------4---------5---------6---------712
      subroutine profile
c
c     -- Profile/Sub-driver for ozone retrieval/N. Nath & C. Wellemeyer
c        /'01-'02
c
            implicit none
c
c     -- Input & Output Variables:
c
            integer     n, nc, ns, nitr, nref
c                             !number of levels, number of channels,
c                             !  number of subchannels in each channel
c	       		!  (<-> slit), maximum number of
c                             !  iterations
            parameter   (n = 81, nc = 10, ns = 8, nitr = 8)

            real        tol_res
                              !residue limit
            parameter   (tol_res = 0.01)

            logical     ssonly, normal

            integer     pcode !print code: 0 - all, 1 - less, 
                              !  2 - no iteration, 3 - summary
c
            real        q_a(n), t_lvl(n), sza,
                              !a-priori profile of layer ozone amount
                              !  atmospheric pressure (atmos),
                              !  atmospheric temperature (K), and
                              !  solar zenith angle (degree)
c    1                  a0(nc,ns), a1(nc,ns), a2(nc,ns), b(nc,ns),
c    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
     3                  cov_a(n-1,n-1)
                              !  a-priori covariance matrix & 

            real        q(n), !retrieved layer ozone profile (DU)
     1                  nval_s(nc), nval(nc), 
                              !retrieved single-scattering, multiple-
                              !  scattering, & total n-values
     2                  kern_s(nc,n-1), kern(nc,n-1),
                              !retrieved single-scattering & total
                              !  kernels (multiple-scattering kernel 
                              !  fixed at kern_m0)
     3                  kern_s0(nc,n-1), kern0(nc,n-1),
                              !retrieved single-scattering & total
                              !  kernels for first iteration 
     4                  avkern(n-1,n-1), cov(n-1,n-1),
                              !averaging kernel & covariance of the
                              !  final solution
     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
     6                  kern_m0(nc,n-1)
                              !multiple scattering kernel
                              !  computed at q0 (not needed if ssonly is
                              !  true)

            integer     itrc  !number of iterations required for 
                              !  convergence
c
            integer     nu    !number of Umkehr layers in a-priori
                              !  ozone & temperature profiles
            parameter   (nu = 13)

            integer     nc1, nuk
                              !beginning index of msr wavelengths
                              !no. of umkehr layers in msr kernel
            parameter   (nc1 = 5, nuk = 11)

            integer     n_fine !number of fine layers in multiple
                               !  scattered kernel
            parameter   (n_fine = 80)

            real        krngr_fine(nc,n_fine), krncl_fine(nc,n_fine)
                              !interpolated msr kernel in fine layers
                              !for terrain pressure and cloud pressure

            real        end1, endn
                              !end parameters for spline interpolation
            parameter   (end1 = 0.0, endn = 0.0)

            real        slope, sigmaa
                              !  slope in 'interpol_qo3'
                              !  slope expressed in sigma form

            real        lam1, lamn
                              !end parameters for spline fits for both
                              !  q_o3 & t interpolns
            parameter   (lam1 = 0.0, lamn = 0.0)

            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

            real        dr    !total change in reflectivity
c
c ..     adding 20-layer profile-related, & associated, variables
c
            real        rmix(80) !layer ozone mixing ratio (ppmv)

            real        q21(21), qa21(21), q021(21),
                              !21-layer retrieved, apriori, &
                              !  first guess profile
     1                  err_q20(20), qtot21, err_qtot,
                              !20-layer estimated measurement error
                              !  in q (percent), total q, & estimated
                              !  error in total q (percent) (same as
                              !  that in r
     2                  err_q(80),
                              !80-layer estimated error in q (percent)
     3                  kern20(nc,20), avkern20(20,20)
                              !10 ch x 20 layer transformed kernel, &
                              !  20-layer averaging kernel

            real        rmix_pres(15), err_pres(15)
                              !ozone mixing ratio (ppmcv) and its 
                              !  estimated measurement error (percent)
                              !  in prescribed levels
c
c     -- Other Variables:
c
            real        p_top, p_bottom, h
                              !pressure at top & bottom of atmosphere
                              !  (atmos), constant scale height (cm)
            parameter   (p_top = 1.0e-4, p_bottom = 1.0, h = 7.996e5)
c
            integer     i, k, j, l, ic
c
      include 'prfprm.com'
c
c           o  Interpolate q_a & t_a from Umkehr data
c
      normal = .false.
      ssonly = .false.
      call interpol_qo3(qu_a, p_top, q_a, slope, nu, n, lam1, lamn,
     1  normal)
      sigmaa = -(1.0/(h*slope))
c
      call interpol_t(tu_a, p_top, t_lvl, nu, n, lam1, lamn, normal)
c
c           o  Interpolate multiple scattered kernel to fine layers
c              for clear and cloud seperately, then mix
c
c     print*,'calling interpol_kern 1',p_top, p_terr
c     write(6,'(6(e12.4))') ((krngr_mu(i,j),i=5,10),j=1,11)
      call interpol_kern(krngr_mu, p_top, p_terr, krngr_fine, nc,
     1  nc1, nuk, n_fine, end1, endn)
c
c     print*,'msr terrain kernel'
c     write(6,'(6(e12.4))') ((krngr_fine(i,j),i=5,10),j=1,80)
c     stop 1
c
c     print*,'calling interpol_kern 2',p_top, p_cld
c     write(6,'(6(e12.4))') ((krncl_mu(i,j),i=5,10),j=1,11)
      call interpol_kern(krncl_mu, p_top, p_cld, krncl_fine, nc,
     1  nc1, nuk, n_fine, end1, endn)
c
c     print*,'msr cloud kernel'
c     write(6,'(6(e12.4))') ((krncl_fine(i,j),i=5,10),j=1,80)
c     stop 1
c
      do j = 1,nc
         do k = 1,n-1
            kern_m0(j,k) = cldfrac(j) * krncl_fine(j,k) +
     1         (1.0 - cldfrac(j)) * krngr_fine(j,k)
         enddo
      enddo
c
c     print*,'sampled msr kernel'
c     write(6,'(6(e12.4))') ((kern_m0(i,j),i=5,10),j=1,80)
c     stop 1
c
c           o  Define a priori covariance matrix
c
      do j = 1,n-1
         cov_a(j,j) = (sigmaq**2)*q_a(j)**2
                              !sigmaq is the fractional error in q
         do l = 1,j-1
            cov_a(j,l) = (sigmaq**2)*q_a(j)*q_a(l)*
     1        exp(-abs(j-l)/corrlen)
            cov_a(l,j) = cov_a(j,l)
         end do
      end do
c
c           o  Retrieve ozone amount in layers 1 through n
c
      adjust_r = - tozsns/rflsns
                              !the required change in reflectivity, such
                              !  that the change in mult. scatt. n-value
                              !  in 331 nm channel due to unit change in
                              !  qtot is zero
      do ic = 1,4
         kern_mr(ic) = 0.0
      end do
      do ic = 5,10
         kern_mr(ic) = rsmsr(ic-4)
                              !note: rsmsr(1) corresponds to channel 5
      end do

      call o3_retrieval(p_lvl, t_lvl, sza10, a0,a1 ,a2 ,b, p0, p2, w,
     1  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  write detail from profile retrieval to unit 15
c
c     p_lvl(81), t_lvl(81), nval_obs(10), q0(81), q_a(81),
c     rad_m0(10), kern_m0(10,80), kern_mu(10,11), sigmae, sigmaq,
c     corrlen, tol_q, q(81), nval_s(10), nval(10), kern_s0(10,80),
c     kern0(10,80), kern_s(10,80), kern(10,80), avkern(80,80),
c     cov(80,80), qu_a(13), tu_a(13), wvl_c(10), sigma0, sigmaa,
c     sigma, qtot0, qtot, res0(10), res(10), itrc
c
c     print*,'q0'
c     write(6,'(10f7.3)') q0
c     print*,'q_a'
c     write(6,'(10f7.3)') q_a
c     print*,'q'
c     write(6,'(10f7.3)') q
c
      call dtailp(p_lvl, t_lvl, nval_obs, q0, q_a, rad_m0, kern_m0,
     1  krngr_mu, krncl_mu, sigmae, sigmaq, corrlen, tol_q, q, nval_s,
     2  nval, kern_s0, kern0, kern_s, kern, avkern, cov, wvl_c, sigma0,
     3  sigmaa, sigma, qtot0, qtot, res0, res, itrc)
c
c     print*,'finished calling dtailp'
c     close(15)
c     stop 1
c
c ..        o  Convert 80-layer quantities into 20-layer quantitities;
c              determine mixing ratio & its error at prescribed levels;
c              and load output buffer with 20-layer profile data
c
      call convert20(q, q_a, q0, avkern, cov, cov_a, cov_e, kern, q21,
     1  qa21, q021, err_q20, qtot21, err_qtot, err_q, kern20, avkern20,
     2  80, 20, nc, .false.)

      call mixratio(p_lvl, rmix, err_q, rmix_pres, err_pres, 81, 15)

      call lodprf(qa21, q021, q21, err_q20, qtot21, err_qtot, rmix_pres,
     1  err_pres, res0, res, kern20, nval_s, tu_a, itrc, dr)
                        !qtot21 should be same as qtot
c
      return
      end