c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine interpol_qo3(q, p_top, q_fine, slope, n, n_fineplus,
     1   lam1, lamn, normal)
c
c     -- Interpol_qo3/Interpolation of layer ozone from Umkehr to fine
c        layers/N. Nath/June '01
c
c     -- Function:  This routine interpolates layer ozone amount (q)
c        data from Umkehr layers to fine layers, as a preparation for 
c        computing back-scattered solar radiation from the SBUV 
c        experiment.  The fine layers, like Umkehr layers, are of equal 
c        delta[log(p)], where p is the pressure, and the output data is
c        in the same unit as Umkehr data (DU).  The Umkehr data may be
c        given from bottom to top of the atmosphere (i.e., starting with
c        layer 0), or in the reverse order.  The fine layer data are 
c        from bottom to top.  Spline interpolation of log(q) is used for
c        the lower layers, and linear interpolation is used for the top
c        layers starting at the middle of the next-to-top Umkehr layer.
c        Following the interpolation, the amount in the topmost layer,
c        which extends to infinity, is calculated based on the assumed
c        exponential drop of q (same as linear drop of log(q)) in top 
c        fine layers.  (Check if q's are positive before entering this 
c        routine.)
c
c     -- Note:  Ozone partial pressure (atmos) is related to the output
c        q (DU) and dz (cm) by: q/1000 = p_o3*(1 + z_phys/r)**2*dz, 
c        where z_phys is the physical height corresponding to z and r is
c        the radius of the Earth.
c
            implicit none
c
c     -- Input & Output Variables:
c
            real        q(n), p_top
                              !layer ozone amount in Umkehr layers (DU),
                              !  atmospheric pressure at top of the 
                              !  uppermost fine layer (atmos) (input)

            real        q_fine(n_fineplus), slope
                              !interpolated ozone amount in fine layers
                              !  (DU) (includes amount in topmost layer,
                              !  an extension of fine layers extending
                              !  to infinity--computed at the end after
                              !  interpolation), slope of log[q_fine]
                              !  vs. z at the top layers (output)

            integer     n, n_fineplus
                              !number of Umkehr layers & fine layers
                              !  (including the layer extending to 
                              !  infinity) (input)

            real        lam1, lamn
                              !end parameters for spline interpolation
                              !  (input)

            logical     normal!if true, Umkehr layer data are ordered
                              !  from 0 upward; if false, the opposite
                              !  (input)
c
c     -- Local Variables:
c
            integer     nl, nl_fine, ns
                              !maximum values of n and n_fine, & number
                              !  of auxiliary q's 
            parameter   (nl = 13, nl_fine = 81, ns = 8)
c
            real        h     !pressure scale height at z = 0 & standard 
                              !  temperature (cm) (i.e., p proportional
                              !  to exp(-z/h))
            parameter   (h = 7.996e5)

            real        q0(nl), z(nl), dz, z_fine(nl_fine), dz_fine
                              !q save array, z array & dz for Umkehr 
                              !  layers, z array & dz for fine layers

            real        logqs(nl), logq_fine(nl_fine)
                              !log(scaled-q), log(q_fine)

            real        r     !ratio of successive q_fine's at the top

            real        z_aux(nl*ns), logq_aux(nl*ns), q_aux(nl*ns),
     1                  dz_aux, qsum(nl), sum
                              !quantities needed for auxiliary computns

            real        devn2, devn, tol
            parameter   (tol = 1e-4)

c            real        maxdevn
c            integer     imax

            integer     n_fine, n_spline, i, is, n_aux, irpt
c
c           Order data from bottom to top
c
      if (.not. normal) then
         do i = 1,n
            q0(i) = q(i)
         end do
         do i = 1,n
            q(i) = q0(n+1-i)
         end do
      end if
c
c           Prepare for interpolation
c
      dz = h*log(2.0)         !define dz & z array
      do i = 1,n
         z(i) = dz/2.0 + (i - 1)*dz
      end do
c
      n_fine = n_fineplus - 1 !number of actual fine layers

      dz_fine = - h*log(p_top)/n_fine
                              !define dz_fine & z_fine array
      do i = 1,n_fine
         z_fine(i) = dz_fine/2.0 + (i - 1)*dz_fine
      end do

      n_spline = 0            !count how many z_fine's are below z(n-1),
                              !  the extent of spline interpolation 
      do i = 1,n_fine
         if (z_fine(i) .gt. z(n-1)) go to 100
         n_spline = n_spline + 1
      end do
  100 continue
c
c           Scale q to fine layers & convert to log
c
      do i = 1,n-1            !note: layer n being excluded (see later)
         logqs(i) = log(q(i)*dz_fine/dz)
      end do
c
c           o  Iterate on the scaling so the layer ozone amount is
c              unchanged
c
      dz_aux = dz/ns
      n_aux = (n - 1)*ns      !note: layer n being excluded
      do i = 1,n_aux
         z_aux(i) = dz_aux/2.0 + (i - 1)*dz_aux
      end do
      do irpt = 1,8           !four iterations should suffice
         call spline(z, logqs, n-1, z_aux, logq_aux, n_aux, lam1, lamn)
                              !first step in finding how the choice of 
                              !  qs affects q in Umkehr layers, as a 
                              !  result of spline interpolation
         do i = 1,n_aux
            q_aux(i) = exp(logq_aux(i))
         end do
         do i = 1,n-1
            sum = 0.0
            do is = 1,ns
               sum = sum + q_aux((i-1)*ns + is)
            end do
            qsum(i) = sum*dz_aux/dz_fine
                              !qsum(i) should be close to q(i)
            logqs(i) = logqs(i) + log(q(i)/qsum(i))
                              !correction so in the next step qsum(i)
                              !  will be closer to q(i)
         end do
         devn2 = 0.0
c         imax = 0
c         maxdevn = 0.0
         do i = 1,n-1
            devn = abs(1.0 - qsum(i)/q(i))
            devn2 = devn2 + devn**2
c            if (devn .gt. maxdevn) then
c               imax = i
c               maxdevn = devn
c            end if
         end do
         devn = sqrt(devn2/(n-1))
c         write (*, *) qsum
c         write (*, *) "devn, imax, maxdevn: ", devn, imax, maxdevn
         if (devn .lt. tol) go to 200
      end do
  200 continue
c
c           Perform spline interpolation for lower layers
c
      call spline(z, logqs, n-1, z_fine, logq_fine, n_spline, lam1, 
     1  lamn)
                              !note: q(n) is excluded, as that includes
                              !  ozone up to pressure zero and is not 
                              !  representative of the layer itself
                              !extrapolation allowed at both ends with
                              !  end parameter -1, though on the upper
                              !  end actual extrapolation is not done
c
c           Perform linear interpolation for top layers
c
      slope = (logqs(n-1) - logqs(n-2))/(z(n-1) - z(n-2))
                              !note: q(n) is excluded
      do i = n_spline+1,n_fine
         logq_fine(i) = logqs(n-1) + slope*(z_fine(i) - z(n-1))
      end do
c
c           Convert log(q_fine) to q_fine
c
      do i = 1,n_fine
         q_fine(i) = exp(logq_fine(i))
      end do
c
c           Compute q at the topmost layer
c
      r = q_fine(n_fine)/q_fine(n_fine-1)
      q_fine(n_fineplus) = r*q_fine(n_fine)/(1.0 - r)
                              !constant ratio r in successive layers,
                              !  with (1 - r)**-1 coming from the 
                              !  series 1 + r + r**2 + ...
c
c           Return original data
c
      if (.not. normal) then
         do i = 1,n
            q(i) = q0(i)
         end do
      end if
c
      return
      end