c2--567--1---------2---------3---------4---------5---------6---------712
c     Program apriori
c
c     -- Function:  Given latitude and day of the year, this program
c        determines and outputs the 21-layer apriori ozone profile used
c        in SBUV Version 8 ozone retrieval package
c     PROGRAM CURRENT AS OF MAY 2004.
c
c     -- Input & Output Variables
c
            real     lat   !latitude (degree) (input)
            integer  day   !day of the year (input)

            real     qu(13)
                     !ozone amounts in Umkehr layers (12:0) (DU)
                     !  (output)

            real     q(81), q21(21), qtot
                     !ozone amounts in fine layers (1:81), 21-layer
                     !  scheme (1:21), and in total column, with
                     !  layer 1 at the bottom of the atmosphere
                     !  (DU) (output) 
            real     presly(22)
                     ! Pressure at 21 layer boundaries (atmospheres)

           data presly / 1.,  0.631,   0.398,   0.251,   0.158,
     1                  .1,  0.0631,  0.0398,  0.0251,  0.0158,
     2                 .01, 0.00631, 0.00398, 0.00251, 0.00158,
     3               0.001,0.000631,0.000398,0.000251,0.000158,
     4              0.0001,0./

c23456789112345678921234567893123456789412345678951234567896123456789712
c
c
c     -- Local Variables
c
            real*4   qu_ref(13,18,12)
                     !climatological ozone amounts in 13 Umkehr layers
                     !  arranged from layer 12 to layer 0, in 18 10-degree
                     !  latitude bins from south to north, and
                     !  12 months, January to December (DU) (input)
            real     cntr(18)
            data     cntr  /-85,-75,-65,-55,-45,-35,-25,-15,-5,
     1                        5,15,25,35,45,55,65,75,85/
            real     mon, fm, fl
            integer  m1, m2, l1, l2, iu, i, i1

            real     ptop, slope, lam1, lam2
            integer  nu, n
            logical  normal
            parameter   (ptop = 1e-4, lam1 = 0.0, lam2 = 0.0,
     1                      nu = 13, n = 81, normal = .false.)
c
c           Read Umkehr-layer ozone climatology data
c
      open (unit=1, file='climat.txt', status='old', form='formatted')
      read (1, '(13f7.3)') qu_ref
      close (1)
c
c           Input latitude and day
c
      write (*, *)
      write (*, *) 'SBUV V8 Apriori'
      write (*, *)
      write (*, *)
     1 'When prompted, input latitude (degree) and day of the year ',
     2 '(in integer) and press enter:  The output is latitude- and ',
     3 'season-dependent 21-layer apriori ozone profile and total ',
     4 'column ozone, in DU, with layer bottoms at 1.0, 0.631, 0.398, ',
     5 '0.251, 0.158, 0.1, 0.0631, 0.0398, 0.0251, 0.0158, 0.01, ',
     6 '0.00631, 0.00398, 0.00251, 0.00158, 0.001, 0.000631, ',
     7 '0.000398, 0.000251, 0.000158, 0.0001 atmos, and top layer ',
     8 'extending to infinity.                       '
      write (*, *)
      write (*, *)
     1 'For another latitude and day, repeat the process; to end ',
     2 'the program, press cntl-c.                   '
      write (*, *)

  100 write (*,*) 'lat, day ?'
      read (*,*) lat, day
      if ((lat .lt. -90.0) .or. (lat .gt. 90.0)) then
         write (*, *) 'lat must lie between -90.0 and 90.0'
         go to 100
      end if
      if ((day .lt. 1) .or. (day .gt. 366)) then
         write (*, *) 'day should be an integer between 1 and 366'
         go to 100
      end if
c
c           Perform bi-linear interpolation in latitude and month
c
      l1 = int((lat + 85.0)/10.0) + 1
      if (l1 .le. 0) l1 = 1
      if (l1 .ge. 18) l1 = 17
      l2 = l1 + 1
      fl = (lat - cntr(l1))/(cntr(l2) - cntr(l1))

      mon = (day + 15.25)/30.5
      m1 = int(mon)
      m2 = m1 + 1
      fm = mon - m1
      if (m1 .eq. 0) m1 = 12
      if (m2 .eq. 13) m2 = 1

      do iu = 1,13
         qu(iu)=(1.0 - fl)*(1.0 - fm)*qu_ref(iu,l1,m1)
     1      + fl*(1.0 - fm)*qu_ref(iu,l2,m1)
     2      + (1.0 - fl)*fm*qu_ref(iu,l1,m2)
     3      + fl*fm*qu_ref(iu,l2,m2)
      end do
c
c           Perform interpolation for 81-layer profile
c
      call interpol_qo3(qu, ptop, q, slope, nu, n, lam1, lam2, normal)
c
c           Convert 81-layer profile to 21-layer profile
c
      qtot = 0.0
      do i = 1,20
         i1 = (i-1)*4 + 1
         q21(i) = q(i1) + q(i1+1) + q(i1+2) + q(i1+3)
         qtot = qtot + q21(i)
      end do
      q21(21) = q(81)
      qtot = qtot + q21(21)
c
c           Output 21-layer profile
c
      write (*, *) 'profile (bottom-to-top):'
      write (*, '(2x, 10f7.3)') q21
      write (*, *) 'column:'
      write (*, '(2x, f7.3)') qtot
      write (*, *)
      go to 100
c
      end

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 are
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     -- 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)

            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
         do i = 1,n-1
            devn = abs(1.0 - qsum(i)/q(i))
            devn2 = devn2 + devn**2
         end do
         devn = sqrt(devn2/(n-1))
         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
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

c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine spline(x0, y0, n0, x, y, n, lam1, lamn0)
c
c     -- Spline/Spline interpolation/N. Nath/April '01.  (Adapted from
c        a previous version, this routine has the following features:
c        (1) it incorporates different types of end conditions; (2) the
c        argument array {x0} may be an increasing or a decreasing
c        sequence; (3) {x} need not be ordered; (4) for x lying outside
c        {x0}, the cubic in the nearest interval is used; this feature
c        should be used with extreme caution.)
c
c     -- Function:  Performs a cubic spline interpolation for function
c        values (y) at an array of arguments (x) of dimension n, with
c        different end parameters lam1 & lamn0.  x0 & y0 (of dimension
c        n0) are the input argument & associated function array.  lam1
c        [alternately lamn0] may take the following values:  0.0, -1.0,
c        -2.0, 1.0, which correspond respectively to y20(1) = (0.0, 0.5,
c        1.0)*y20(2), & y10(1) = 0.0 [y20(n0) = (0.0, 0.5, 1.0)*
c        y20(n0-1), & y10(n0) = 0.0]; y20 and y10 are respectively the
c        second and first derivative at x0.
c
c     -- Usage Notes:  (1) Choice of lam1 and lamn0 depends on the
c        knowledge of the function.  If nothing is known, use one of the
c        first two values.  Note also that the second & third values
c        correspond respectively to a single cubic in the domain:
c        (2.0*x0(1) - x0(2), x0(2)), with y20 at the left boundary =
c        0.0, & a quadratic in the domain (x0(1), x0(2)), and similarly
c        for the other end.  (2) x0 must be monotonically increasing or
c        decreasing.  (3) n0 must be > 2; if it exceeds 256, change the
c        value of n0max.
c
c     -- Input & Output Variables:
c
            real        x0(n0), y0(n0),
                              !base argument & function array (input),
     1                  x(n), y(n),
                              !argument array at which function values
                              !  are desired (input), & interpolated
                              !  function array (output)
     2                  lam1, lamn0
                              !end parameters at x0(1) & x0(n0) (input)

            integer     n0, n !dimensions of (x0, y0) & (x, y) (input)
c
c     -- Local Variables:
c
            integer     n0max
            parameter   (n0max = 256)
c
            real        y20(n0max) !second derivative array

            real        h(n0max), a(n0max), b(n0max), d(n0max),
     1                  phi, lam, dy0, dy00, c, c1, f

            integer     n01, l, u, m, i, j, k, i1

            logical     reverse
c
c     -- Subprogram Called:  order
c
c     -- Ref.:  Ahlberg et al, Theory Of Splines & Their Applications,
c        Academic Press, 1967; Hamming, Numerical Methods, Mcgraw-Hill,
c        1973.
c
      phi(lam) = lam*(lam + 1.0)*(lam + 2.0)
c
c           Order x0 in an increasing sequence; change x accordingly
c                 !note:  this ordering, viz., x0 -> -x0 while y0 -> y0,
c                 !  only changes the signs of the odd derivatives; thus
c                 !  y20 remains unchanged
c
      reverse =  (x0(1) .gt. x0(n0))
      if (reverse) call order(x0, n0, x, n)
c
c           Determine y20 by solving tri-diagonal linear equations
c
c           o  Define h, a, & d (quantities independent of y0)
c
      n01 = n0 - 1
      do i = 1,n01
         h(i) = x0(i+1) - x0(i)
      end do
      a(1) = -lam1/2.0
      do i = 2,n01
         i1 = i - 1
         d(i) = h(i1)*a(i1) + 2.0*(h(i1) + h(i))
         a(i) = -h(i)/d(i)
      end do
      d(n0) = h(n01)*(lamn0*a(n01) + 2.0)
c
c           o  Define b, & complete computations
c
      dy0 = (y0(2) - y0(1))/h(1)
      b(1) = (phi(lam1)*dy0)/h(1)/2.0
      do i = 2,n01
         dy00 = dy0
         dy0 = (y0(i+1) - y0(i))/h(i)
         b(i) = (6.0*(dy0 - dy00) - h(i-1)*b(i-1))/d(i)
      end do
      b(n0) = (-phi(lamn0)*dy0 - lamn0*h(n01)*b(n01))/d(n0)
      y20(n0) = b(n0)
      k = n0
      do i = 2,n0
         k = k - 1
         y20(k) = a(k)*y20(k+1) + b(k)
      end do
c
c           Interpolate using binary domain search
c
      do j = 1,n
         l = 1
         u = n0
         m = (l + u)/2
         do while (m .gt. l)
            if (x(j) .lt. x0(m)) then
                  u = m
               else
                  l = m
            end if
            m = (l + u)/2
         end do
         c = (x0(l+1) - x(j))/h(l)
         c1 = 1.0 - c
         f = (h(l)**2)/6.0
         y(j) = c*(y0(l) - f*y20(l)*(1.0 - c**2))
     1            + c1*(y0(l+1) - f*y20(l+1)*(1.0 - c1**2))
      end do
c
c           Reorder x0 and x
c
      if (reverse) call order(x0, n0, x, n)
c
      return
      end

c     ---
      Subroutine order(x0, n0, x, n)
c
c     -- Order/Adjunct to Spline/N. Nath/April '01
c
c     -- Input & Output Variables:
c
            real        x0(n0), x(n)

            integer     n0, n
c
c     -- Local Variables:
c
            integer     i, j
c
c           Change signs of x0 and x
c
      do i = 1,n0
         x0(i) = -x0(i)
      end do
      do j = 1,n
         x(j) = -x(j)
      end do
c
      return
      end


