c2--567--1---------2---------3---------4---------5---------6---------712
      Subroutine update(x, y, kern, y_obs, x_a, s_a, s_e, x_u,
     1   kern_aux, s_aux, n, m, pindex, pivot)
c
c     -- Update/Update ozone profile using C. Rodgers' method/N. Nath
c        /May '01
c
c     -- Function:  This routine implements the Rodgers solution for the
c        updated ozone profile, symbolized by:
c
c              x(i+1) = x_a + S_a*K_t(i) * (K(i)*S_a*K_t(i) + S_e)**-1 *
c                         ((y - F(x(i)) + K(i)*(x(i) - x_a))
c
c        where i+1 and i refer to the current and previous iterations, 
c        x is the profile vector (of dimension n), y is the measurement 
c        vector (of dimension m < n), x_a is the a-priori profile 
c        vector, F(x) is the model vector simulating y, K is the kernel
c        or the weighting function matrix (del(F)/del(x), of dimension
c        m by n), 't' stands for transpose, and S_e and S_a are the 
c        covariance matrices connected with the measurement and the 
c        a-priori profile (of dimensions m by m and n by n),
c        respectively.  Both S_a and S_e are symmetric.
c
c        The LU factorization of the matrix, K(i)*S_a*K_t(i) + S_e, is
c        used in the evaluation of the above expression.
c
c        On convergence, the covariance of the solution can be readily
c        computed from:
c
c              S = S_a - S_a*K_t(i)*(K(i)*S_a*K_t(i) + S_e)**-1*K(i)*S_a
c
c        by using the matrix K(i)*S_a and the LU-decomposed version of 
c        K(i)*S_a*K_t(i) + S_e, which are provided as output.  Similar
c        comment applies for computation of the averaging kernel, 
c        i.e., S_a*K_t(i)*(K(i)*S_a*K_t(i) + S_e)**-1*K(i).  (In either
c        case, remember that the operand on which the inverse works is
c        overwritten in the LU approach.)
c
c     -- Notes:  Cholesky decomposition, QR decomposition, and singular
c        value decomposition are possible alternatives to LU decom-
c        position; the first two for efficiency, and SVD to handle 
c        extreme ill-conditioning.  With inclusion of a-priori cons-
c        traint, the inverse problem is already regularized; thus, SVD 
c        should not be needed.
c
            implicit none
c
c     -- Input & Output Variables:
c
            real        x(n), y(m), kern(m,n),
                              !profile to be updated, and the corres-                              
                              !  ponding values of the vector in the 
                              !  measurement space and the weighting 
                              !  function matrix (input)
     1                  y_obs(m),
                              !observed vector in the measurement space
                              !  (input)
     2                  x_a(n), s_a(n,n), s_e(m,m)
                              !a-priori profile, a-priori covariance 
                              !  matrix, & measurement-error covariance
                              !  matrix (input)

            real        x_u(n), kern_aux(m,n), s_aux(m,m)
                              !updated profile, kern*s_a & LU-factored
                              !  form of kern*s_a*kern_t + s_e (output)

            integer     n, m  !dimensions of x and y (input)

            integer     pindex(m)
                              !row permutation index of kern*s_a*kern_t
                              !  (output)

            logical     pivot !true <-> partial pivoting, false <-> no
                              !  pivoting (input)
c
c     -- Local Variables:
c
            integer     ml
            parameter   (ml = 10)

            real        diff(ml)
c >>>
c            real        s_auxc(ml,ml), s_auxcinv(ml,ml), prod(ml,ml)
c >>>
            integer     i, j, k, l
c
c     -- Ref:  C. Rodgers, Inverse Methods for Atmospheric Sounding, 
c        World Scientific, 2000.
c
c           Compute K_aux = K*S_a
c
      do i = 1,n              !organization of indices:  use i & j in 
                              !  the space of dimension n; k & l in 
                              !  the space of dimension m; and remember
                              !  [kern] = (m,n), [s_a] = (n,n), [s_e] =
                              !  (m,m), [y] = m, [x] = n 
         do k = 1,m
            kern_aux(k,i) = 0.0
            do j = 1,n
               kern_aux(k,i) = kern_aux(k,i) + kern(k,j)*s_a(j,i)
            end do
         end do
      end do
c
c           Compute S_aux = K*S_a*K_t + S_e
c                 !remember symmetry of S_aux
c
      do l = 1,m
         do k = 1,l
            s_aux(k,l) = s_e(k,l)
            do j = 1,n
               s_aux(k,l) = s_aux(k,l) + kern_aux(k,j)*kern(l,j)
            end do
            if (k .ne. l) s_aux(l,k) = s_aux(k,l)
         end do
      end do
c
c >>>       Copy s_aux to s_auxc to check correctness of inversion
c
c      do l = 1,m
c         do k = 1,m
c            s_auxc(k,l) = s_aux(k,l)
c         end do
c      end do
c  >>>
c
c           Factor S_aux into lower (L) and upper (U) triangular 
c           matrices, & overwrite S_aux with its decomposed form
c                 !decomposition: pindex*S_aux = L*U, where pindex is
c                 !  a permutation symbol connected with partial 
c                 !  pivoting, & diagonal elements of L are all 1,
c                 !  which are not written
c
      call factor(s_aux, pindex, m, pivot)
c
c >>>       
c
c      do l = 1,m              !define unit matrix
c         do k = 1,m
c            s_auxcinv(k,l) = 0.0
c         end do
c         s_auxcinv(l,l) = 1.0
c      end do
c      call inverse(s_aux, s_auxcinv, pindex, m, m, pivot)
c                              !find inverse of s_auxc; first parameter
c                              !  is s_aux
c      do l = 1,m              !check if s_auxc*s_auxcinv = 1
c         do k = 1,m
c            prod(k,l) = 0.0
c            do j = 1,m
c               prod(k,l) = prod(k,l) + s_auxc(k,j)*s_auxcinv(j,l)
c            end do
c         end do
c      end do
c      write (*, *) "prod: ", prod
c >>>
c
c           Compute S_aux**-1*((y - F(x) + K*(x - x_a))
c                 !the vector ((y - F(x) + K*(x - x_a)) is overwritten
c                 !  in the process
c
      do l = 1,m
         diff(l) = y_obs(l) - y(l)
         do j = 1,n
            diff(l) = diff(l) + kern(l,j)*(x(j) - x_a(j))
         end do
      end do
      call inverse(s_aux, diff, pindex, m, 1, pivot)
c
c           Compute updated x
c
      do i = 1,n
         x_u(i) = x_a(i)
         do l = 1,m
            x_u(i) = x_u(i) + kern_aux(l,i)*diff(l)
         end do
      end do
c
      return
      end