CRTM_Utility.f90

Go to the documentation of this file.

!
! CRTM_Utility 
!
! Module containing CRTM numerical utility routines.
!
! NOTE: This utility is specifically for use with RTSolution codes.
!       Adapted from package of MOM 1991, ASYMTX adapted from DISORT.
!
!
! CREATION HISTORY:
!       Written by:     Quanhua Liu,    Quanhua.Liu@noaa.gov
!                       Yong Han,       Yong.Han@noaa.gov
!                       Paul van Delst, paul.vandelst@noaa.g
!                       08-Jun-2004
!

MODULE crtm_utility 


  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE type_kinds     , ONLY: fp
  USE message_handler, ONLY: display_message, success, failure, warning
  USE crtm_parameters, ONLY: zero, one, two, &
                             point_25, point_5, &
                             pi
  ! Disable implicit typing
  IMPLICIT NONE
  
  
  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE
  PUBLIC :: matinv
  PUBLIC :: double_gauss_quadrature
  PUBLIC :: legendre
  PUBLIC :: legendre_m
  PUBLIC :: asymtx,asymtx_tl,asymtx_ad
  PUBLIC :: asymtx2_ad


  ! -------------------
  ! Procedure overloads
  ! -------------------
  INTERFACE legendre 
    MODULE PROCEDURE legendre_scalar
    MODULE PROCEDURE legendre_rank1
  END INTERFACE 


  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module 
  CHARACTER(*), PARAMETER :: module_version_id = & 
  '$Id: CRTM_Utility.f90 60152 2015-08-13 19:19:13Z paul.vandelst@noaa.gov $'
  ! Numerical small threshold value in Eigensystem
  REAL(fp), PARAMETER :: eigen_threshold = 1.0e-20_fp


CONTAINS


  SUBROUTINE double_gauss_quadrature(NUM, ABSCISSAS, WEIGHTS)
    ! Generates the abscissas and weights for an even 2*NUM point
    ! Gauss-Legendre quadrature.  Only the NUM positive points are returned.
    IMPLICIT NONE
    INTEGER, INTENT(IN) ::  num
    REAL( fp), DIMENSION(:) ::   abscissas, weights
    INTEGER  n, k, i, j, l
    REAL(fp) ::   x, xp, pl, pl1, pl2, dpl
    ! REAL(fp), PARAMETER :: TINY1=3.0E-14_fp
    REAL(fp) :: tiny1
    parameter(tiny1=3.0d-14)
    n = num
    k = (n+1)/2
    DO j = 1, k
      x = cos(pi*(j-point_25)/(n+point_5))
      i = 0
    100  CONTINUE
      pl1 = 1
      pl = x
      DO l = 2, n
        pl2 = pl1
        pl1 = pl
        pl = ( (2*l-1)*x*pl1 - (l-1)*pl2 )/l
      END DO
      dpl = n*(x*pl-pl1)/(x*x-1)
      xp = x
      x = xp - pl/dpl
      i = i+1
      IF (abs(x-xp).GT.tiny1 .AND. i.LT.10) GO TO 100
        abscissas(j) = (one-x)/two
        abscissas(num+1-j) = (one+x)/two
        weights(num+1-j) = one/((one-x*x)*dpl*dpl)
        weights(j)       = one/((one-x*x)*dpl*dpl)
    END DO
    RETURN
  END SUBROUTINE double_gauss_quadrature
!
!
  FUNCTION matinv(A, Error_Status)
    ! --------------------------------------------------------------------
    ! Compute the inversion of the matrix A
    ! Invert matrix by Gauss method
    ! --------------------------------------------------------------------
    IMPLICIT NONE
    REAL(fp), intent(in),dimension(:,:) :: a
    INTEGER, INTENT( OUT ) :: error_status

    INTEGER:: n
    REAL(fp), dimension(size(a,1),size(a,2)) :: b
    REAL(fp ), dimension(size(a,1),size(a,2)) :: matinv 
    REAL(fp), dimension(size(a,1)) :: temp 
    ! Local parameters
    CHARACTER(*), PARAMETER :: routine_name = 'matinv'
    ! - - - Local Variables - - -
    REAL(fp) :: c, d
    INTEGER :: i, j, k, m, imax(1), ipvt(size(a,1))
    ! - - - - - - - - - - - - - -
    error_status = success     
    b = a
    n=size(a,1)
    matinv=a
    ipvt = (/ (i, i = 1, n) /)
    ! Go into loop- b, the inverse, is initialized equal to a above
    DO k = 1,n
    ! Find the largest value and position of that value
      imax = maxloc(abs(b(k:n,k)))
      m = k-1+imax(1)
    ! sigular matrix check
      IF ( abs(b(m,k)).LE.(1.e-40_fp) ) THEN
        error_status = failure
        CALL display_message( routine_name, 'Singular matrix', error_status )                                          
        RETURN
      END IF    
    ! get the row beneath the current row if the current row will
    ! not compute
      IF ( m .ne. k ) THEN
        ipvt( (/m,k/) ) = ipvt( (/k,m/) )
        b((/m,k/),:) = b((/k,m/),:)
      END IF
    ! d is a coefficient - brings the pivot value to one and then is applied
    ! to the rest of the row
      d = 1/b(k,k)
      temp = b(:,k)
      DO j = 1, n
         c = b(k,j)*d
         b(:,j) = b(:,j)-temp*c
         b(k,j) = c
      END DO 
      b(:,k) = temp*(-d)
      b(k,k) = d
    END DO
    matinv(:,ipvt) = b
  END FUNCTION matinv
!
!
  SUBROUTINE legendre_scalar(MOA,ANG,PL)
    !  calculating Legendre polynomial using recurrence relationship
    IMPLICIT NONE
    INTEGER, INTENT(in) :: MOA
    REAL(fp), intent(in) :: ANG
    REAL(fp), intent(out), dimension(0:MOA):: PL
    INTEGER :: j
    pl(0)=one
    pl(1)=ang
    IF ( moa.GE.2 ) THEN
      DO j = 1, moa -1
        pl(j+1)=REAL(2*J+1,fp)/REAL(J+1,fp)*ANG*PL(J) &
          -REAL(J,fp)/REAL(J+1,fp)*PL(J-1)
      END DO
    END IF
    RETURN
  END SUBROUTINE legendre_scalar
!
!
   SUBROUTINE legendre_rank1(MOA,NU,ANG,PL)
     !  calculating Legendre polynomial using recurrence relationship
     IMPLICIT NONE
     INTEGER, intent(in):: MOA,NU
     REAL(fp), INTENT(in), dimension(NU):: ANG
     REAL(fp), INTENT(out),dimension(0:MOA,NU):: PL
     REAL(fp) :: TE
     INTEGER :: i,j
     DO i = 1,nu
       te=ang(i)
       pl(0,i)=one
       pl(1,i)=te
       IF(moa.GE.2) THEN
         DO j = 1, moa-1
           pl(j+1,i)=REAL(2*J+1, fp)/REAL(J+1, fp)*TE*PL(J,i) &
           -REAL(J, fp)/REAL(J+1, fp)*PL(J-1,i)
         END DO
       END IF
     END DO
     RETURN
   END SUBROUTINE legendre_rank1
!       
!
   SUBROUTINE legendre_m(MF,N,U,PL)
     IMPLICIT NONE
     REAL(fp), INTENT(IN) :: u
     INTEGER, INTENT(IN) :: n,mf
     REAL(fp), DIMENSION(0:N) :: pl 
     REAL(fp) :: f
     INTEGER:: j
     IF ( mf.GT.n ) THEN
       pl = 0.0_fp
     ELSE
       IF ( mf .EQ. 0 ) THEN
         f = 1.0_fp
         pl(0) = 1.0_fp
       ELSE
         f=sqrt(gamma2(2*mf))/gamma2(mf)/(2**mf)
         pl(mf)=f*sqrt((1.0_fp-u*u)**mf)
       END IF
       IF ( n.GT.mf )  pl(mf+1)=u*sqrt(2*mf+1.0_fp)*pl(mf)
       IF( n.gt.(mf+1) ) THEN
         DO j=mf+1,n-1
           pl(j+1)=(float(2*j+1)*u*pl(j) &
           -sqrt(float(j*j-mf*mf))*pl(j-1))/sqrt(float((j+1)**2-mf*mf))
         END DO
       END IF
     END IF     
   END SUBROUTINE legendre_m
!
!
   FUNCTION gamma2(N)
     ! Compute gamma function value
     IMPLICIT NONE
     INTEGER, INTENT(IN) :: n
     INTEGER :: i
     REAL(fp) :: gamma2
       
     IF ( n.lt.0 ) THEN
       CALL display_message('gamma2','Invalid input', warning)
       gamma2 = 1.0_fp
     ELSE
       gamma2=1.0_fp
       IF( n.gt.0 ) THEN
         DO i = 1 , n
         gamma2=gamma2*float(i)
         END DO
        END IF
     END IF
   END FUNCTION gamma2
!
!
   SUBROUTINE  asymtx( AAD, M, IA, IEVEC, &
                            EVECD, EVALD, IER)

!    =======  D O U B L E    P R E C I S I O N    V E R S I O N  ======

!       Solves eigenfunction problem for real asymmetric matrix
!       for which it is known a priori that the eigenvalues are real.

!       This is an adaptation of a subroutine EIGRF in the IMSL
!       library to use real instead of complex arithmetic, accounting
!       for the known fact that the eigenvalues and eigenvectors in
!       the discrete ordinate solution are real.  Other changes include
!       putting all the called subroutines in-line, deleting the
!       performance index calculation, updating many DO-loops
!       to Fortran77, and in calculating the machine precision
!       TOL instead of specifying it in a data statement.

!       EIGRF is based primarily on EISPACK routines.  The matrix is
!       first balanced using the parlett-reinsch algorithm.  Then
!       the Martin-Wilkinson algorithm is applied.

!       References:
!          Dongarra, J. and C. Moler, EISPACK -- A Package for Solving
!             Matrix Eigenvalue Problems, in Cowell, ed., 1984:
!             Sources and Development of Mathematical Software,
!             Prentice-Hall, Englewood Cliffs, NJ
!         Parlett and Reinsch, 1969: Balancing a Matrix for Calculation
!             of Eigenvalues and Eigenvectors, Num. Math. 13, 293-304
!         Wilkinson, J., 1965: The Algebraic Eigenvalue Problem,
!             Clarendon Press, Oxford

!   I N P U T    V A R I A B L E S:

!        AAD  :  input asymmetric matrix, destroyed after solved
!        M    :  order of  A
!       IA    :  first dimension of  A
!    IEVEC    :  first dimension of  EVECD

!   O U T P U T    V A R I A B L E S:

!       EVECD :  (unnormalized) eigenvectors of  A
!                   ( column J corresponds to EVALD(J) )

!       EVALD :  (unordered) eigenvalues of  A ( dimension at least M )

!       IER   :  if .NE. 0, signals that EVALD(IER) failed to converge;
!                   in that case eigenvalues IER+1,IER+2,...,M  are
!                   correct but eigenvalues 1,...,IER are set to zero.
      IMPLICIT NONE
      LOGICAL            bad_status
      CHARACTER(len=200) :: message

!   S C R A T C H   V A R I A B L E S:

!       WKD    :  WORK AREA ( DIMENSION AT LEAST 2*M )
!+---------------------------------------------------------------------+

!  include file for dimensioning


!!      INCLUDE '../includes/VLIDORT.PARS'

!  input/output arguments

      INTEGER :: m, ia, ievec, ier
      INTEGER, PARAMETER :: maxstrmstks = 30
      REAL( fp), PARAMETER :: c1=0.4375_fp,c2=0.5_fp,c3=0.75_fp,c4=0.95_fp,c5=16.0_fp,c6=256.0_fp
      REAL( fp), DIMENSION(:,:) :: aad, evecd
      REAL( fp), DIMENSION(:) :: evald
      REAL( fp) :: wkd(4*maxstrmstks), nor_factor
      

!  local variables (explicit declaration

      LOGICAL           noconv, notlas
      INTEGER :: i, j, l, k, kkk, lll,n, n1, n2, in, lb, ka, ii
!      DOUBLE PRECISION  TOL, DISCRI, SGN, RNORM, W, F, G, H, P, Q, R
      REAL(fp) :: tol, discri, sgn, rnorm, w, f, g, h, p, q, r
      REAL(fp) :: repl, col, row, scale, t, x, z, s, y, uu, vv
!
      ier = 0
      bad_status = .false.
      message = ' '

!       Here change to bypass D1MACH:
      tol = 1.0d-12
!        TOL = D1MACH(4)
      IF ( m.LT.1 .OR. ia.LT.m .OR. ievec.LT.m ) THEN
        message = 'ASYMTX--bad input variable(s)'
        bad_status = .true.
        RETURN
      ENDIF
!                           ** HANDLE 1X1 AND 2X2 SPECIAL CASES
      IF ( m.EQ.1 )  THEN
         evald(1) = aad(1,1)
         evecd(1,1) = 1.0d0
         RETURN
      ELSE IF ( m.EQ.2 )  THEN
         discri = ( aad(1,1) - aad(2,2) )**2 + 4.0d0*aad(1,2)*aad(2,1)
         IF ( discri.LT.zero ) THEN
           message = 'ASYMTX--COMPLEX EVALS IN 2X2 CASE'
           bad_status = .true.
           RETURN
         ENDIF
         sgn = one
         IF ( aad(1,1).LT.aad(2,2) )  sgn = - one
         evald(1) = 0.5d0*( aad(1,1) + aad(2,2) + sgn*sqrt(discri) )
         evald(2) = 0.5d0*( aad(1,1) + aad(2,2) - sgn*sqrt(discri) )
         evecd(1,1) = one
         evecd(2,2) = one
         IF ( aad(1,1).EQ.aad(2,2) .AND. &
               (aad(2,1).EQ.zero.OR.aad(1,2).EQ.zero) ) THEN
            rnorm = abs(aad(1,1))+abs(aad(1,2))+ &
                      abs(aad(2,1))+abs(aad(2,2))
            w = tol * rnorm
            evecd(2,1) = aad(2,1) / w
            evecd(1,2) = - aad(1,2) / w
         ELSE
            evecd(2,1) = aad(2,1) / ( evald(1) - aad(2,2) )
            evecd(1,2) = aad(1,2) / ( evald(2) - aad(1,1) )
         ENDIF
         RETURN
      END IF
!                                        ** INITIALIZE OUTPUT VARIABLES
      DO 20 i = 1, m
         evald(i) = zero
         DO 10 j = 1, m
            evecd(i,j) = zero
10       CONTINUE
         evecd(i,i) = one
20    CONTINUE
!                  ** BALANCE THE INPUT MATRIX AND REDUCE ITS NORM BY
!                  ** DIAGONAL SIMILARITY TRANSFORMATION STORED IN WK;
!                  ** THEN SEARCH FOR ROWS ISOLATING AN EIGENVALUE
!                  ** AND PUSH THEM DOWN
      rnorm = zero
      l  = 1
      k  = m

30    kkk = k
         DO 70  j = kkk, 1, -1
            row = zero
            DO 40 i = 1, k
               IF ( i.NE.j ) row = row + abs( aad(j,i) )
40          CONTINUE
            IF ( row.EQ.zero ) THEN
               wkd(k) = j
               IF ( j.NE.k ) THEN
                  DO 50 i = 1, k
                     repl   = aad(i,j)
                     aad(i,j) = aad(i,k)
                     aad(i,k) = repl
50                CONTINUE
                  DO 60 i = l, m
                     repl   = aad(j,i)
                     aad(j,i) = aad(k,i)
                     aad(k,i) = repl
60                CONTINUE
               END IF
               k = k - 1
               GO TO 30
            END IF
70       CONTINUE
!                                     ** SEARCH FOR COLUMNS ISOLATING AN
!                                       ** EIGENVALUE AND PUSH THEM LEFT
80    lll = l
         DO 120 j = lll, k
            col = zero
            DO 90 i = l, k
               IF ( i.NE.j ) col = col + abs( aad(i,j) )
90          CONTINUE
            IF ( col.EQ.zero ) THEN
               wkd(l) = j
               IF ( j.NE.l ) THEN
                  DO 100 i = 1, k
                     repl   = aad(i,j)
                     aad(i,j) = aad(i,l)
                     aad(i,l) = repl
100               CONTINUE
                  DO 110 i = l, m
                     repl   = aad(j,i)
                     aad(j,i) = aad(l,i)
                     aad(l,i) = repl
110               CONTINUE
               END IF
               l = l + 1
               GO TO 80
            END IF
120      CONTINUE
!                           ** BALANCE THE SUBMATRIX IN ROWS L THROUGH K
      DO 130 i = l, k
         wkd(i) = one
130   CONTINUE

140   noconv = .false.
         DO 200 i = l, k
            col = zero
            row = zero
            DO 150 j = l, k
               IF ( j.NE.i ) THEN
                  col = col + abs( aad(j,i) )
                  row = row + abs( aad(i,j) )
               END IF
150         CONTINUE
            f = one
            g = row / c5
            h = col + row
160         IF ( col.LT.g ) THEN
               f   = f * c5
               col = col * c6
               GO TO 160
            END IF
            g = row * c5
170         IF ( col.GE.g ) THEN
               f   = f / c5
               col = col / c6
               GO TO 170
            END IF
!                                                         ** NOW BALANCE
            IF ( (col+row)/f .LT. c4*h ) THEN
               wkd(i)  = wkd(i) * f
               noconv = .true.
               DO 180 j = l, m
                  aad(i,j) = aad(i,j) / f
180            CONTINUE
               DO 190 j = 1, k
                  aad(j,i) = aad(j,i) * f
190            CONTINUE
            END IF
200      CONTINUE

      IF ( noconv ) GO TO 140
!                                  ** IS -A- ALREADY IN HESSENBERG FORM?
      IF ( k-1 .LT. l+1 ) GO TO 350
!                                   ** TRANSFER -A- TO A HESSENBERG FORM
      DO 290 n = l+1, k-1
         h        = zero
         wkd(n+m) = zero
         scale    = zero
!                                                        ** SCALE COLUMN
         DO 210 i = n, k
            scale = scale + abs(aad(i,n-1))
210      CONTINUE
         IF ( scale.NE.zero ) THEN
            DO 220 i = k, n, -1
               wkd(i+m) = aad(i,n-1) / scale
               h = h + wkd(i+m)**2
220         CONTINUE
            g = - sign( sqrt(h), wkd(n+m) )
            h = h - wkd(n+m) * g
            wkd(n+m) = wkd(n+m) - g
!                                                 ** FORM (I-(U*UT)/H)*A
            DO 250 j = n, m
               f = zero
               DO 230  i = k, n, -1
                  f = f + wkd(i+m) * aad(i,j)
230            CONTINUE
               DO 240 i = n, k
                  aad(i,j) = aad(i,j) - wkd(i+m) * f / h
240            CONTINUE
250         CONTINUE
!                                    ** FORM (I-(U*UT)/H)*A*(I-(U*UT)/H)
            DO 280 i = 1, k
               f = zero
               DO 260  j = k, n, -1
                  f = f + wkd(j+m) * aad(i,j)
260            CONTINUE
               DO 270 j = n, k
                  aad(i,j) = aad(i,j) - wkd(j+m) * f / h
270            CONTINUE
280         CONTINUE
            wkd(n+m)  = scale * wkd(n+m)
            aad(n,n-1) = scale * g
         END IF
290   CONTINUE

      DO 340  n = k-2, l, -1
         n1 = n + 1
         n2 = n + 2
         f  = aad(n1,n)
         IF ( f.NE.zero ) THEN
            f  = f * wkd(n1+m)
            DO 300 i = n2, k
               wkd(i+m) = aad(i,n)
300         CONTINUE
            IF ( n1.LE.k ) THEN
               DO 330 j = 1, m
                  g = zero
                  DO 310 i = n1, k
                     g = g + wkd(i+m) * evecd(i,j)
310               CONTINUE
                  g = g / f
                  DO 320 i = n1, k
                     evecd(i,j) = evecd(i,j) + g * wkd(i+m)
320               CONTINUE
330            CONTINUE
            END IF
         END IF
340   CONTINUE

350   CONTINUE
      n = 1
      DO 370 i = 1, m
         DO 360 j = n, m
            rnorm = rnorm + abs(aad(i,j))
360      CONTINUE
         n = i
         IF ( i.LT.l .OR. i.GT.k ) evald(i) = aad(i,i)
370   CONTINUE
      n = k
      t = zero
!                                         ** SEARCH FOR NEXT EIGENVALUES
380   IF ( n.LT.l ) GO TO 530
      in = 0
      n1 = n - 1
      n2 = n - 2
!                          ** LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
390   CONTINUE
      DO 400 i = l, n
         lb = n+l - i
         IF ( lb.EQ.l ) GO TO 410
         s = abs( aad(lb-1,lb-1) ) + abs( aad(lb,lb) )
         IF ( s.EQ.zero ) s = rnorm
         IF ( abs(aad(lb,lb-1)) .LE. tol*s ) GO TO 410
400   CONTINUE

410   x = aad(n,n)
      IF ( lb.EQ.n ) THEN
!                                        ** ONE EIGENVALUE FOUND
         aad(n,n)  = x + t
         evald(n) = aad(n,n)
         n = n1
         GO TO 380
      END IF

      y = aad(n1,n1)
      w = aad(n,n1) * aad(n1,n)
      IF ( lb.EQ.n1 ) THEN
!                                        ** TWO EIGENVALUES FOUND
         p = (y-x) * c2
         q = p**2 + w
         z = sqrt( abs(q) )
         aad(n,n) = x + t
         x = aad(n,n)
         aad(n1,n1) = y + t
!                                        ** REAL PAIR
         z = p + sign(z,p)
         evald(n1) = x + z
         evald(n)  = evald(n1)
         IF ( z.NE.zero ) evald(n) = x - w / z
         x = aad(n,n1)
!                                  ** EMPLOY SCALE FACTOR IN CASE
!                                  ** X AND Z ARE VERY SMALL
         r = sqrt( x*x + z*z )
         p = x / r
         q = z / r
!                                             ** ROW MODIFICATION
         DO 420 j = n1, m
            z = aad(n1,j)
            aad(n1,j) = q * z + p * aad(n,j)
            aad(n,j)  = q * aad(n,j) - p * z
420      CONTINUE
!                                             ** COLUMN MODIFICATION
         DO 430 i = 1, n
            z = aad(i,n1)
            aad(i,n1) = q * z + p * aad(i,n)
            aad(i,n)  = q * aad(i,n) - p * z
430      CONTINUE
!                                          ** ACCUMULATE TRANSFORMATIONS
         DO 440 i = l, k
            z = evecd(i,n1)
            evecd(i,n1) = q * z + p * evecd(i,n)
            evecd(i,n)  = q * evecd(i,n) - p * z
440      CONTINUE

         n = n2
         GO TO 380
      END IF

      IF ( in.EQ.30 ) THEN
!                    ** NO CONVERGENCE AFTER 30 ITERATIONS; SET ERROR
!                    ** INDICATOR TO THE INDEX OF THE CURRENT EIGENVALUE
         ier = n
         GO TO 670
      END IF
!                                                          ** FORM SHIFT
      IF ( in.EQ.10 .OR. in.EQ.20 ) THEN
         t = t + x
         DO 450 i = l, n
            aad(i,i) = aad(i,i) - x
450      CONTINUE
         s = abs(aad(n,n1)) + abs(aad(n1,n2))
         x = c3 * s
         y = x
         w = - c1 * s**2
      END IF

      in = in + 1
!                ** LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL ELEMENTS

      DO 460 j = lb, n2
         i = n2+lb - j
         z = aad(i,i)
         r = x - z
         s = y - z
         p = ( r * s - w ) / aad(i+1,i) + aad(i,i+1)
         q = aad(i+1,i+1) - z - r - s
         r = aad(i+2,i+1)
         s = abs(p) + abs(q) + abs(r)
         p = p / s
         q = q / s
         r = r / s
         IF ( i.EQ.lb ) GO TO 470
         uu = abs( aad(i,i-1) ) * ( abs(q) + abs(r) )
         vv = abs(p)*(abs(aad(i-1,i-1))+abs(z)+abs(aad(i+1,i+1)))
         IF ( uu .LE. tol*vv ) GO TO 470
460   CONTINUE

470   CONTINUE
      aad(i+2,i) = zero
      DO 480 j = i+3, n
         aad(j,j-2) = zero
         aad(j,j-3) = zero
480   CONTINUE

!             ** DOUBLE QR STEP INVOLVING ROWS K TO N AND COLUMNS M TO N

      DO 520 ka = i, n1
         notlas = ka.NE.n1
         IF ( ka.EQ.i ) THEN
            s = sign( sqrt( p*p + q*q + r*r ), p )
            IF ( lb.NE.i ) aad(ka,ka-1) = - aad(ka,ka-1)
         ELSE
            p = aad(ka,ka-1)
            q = aad(ka+1,ka-1)
            r = zero
            IF ( notlas ) r = aad(ka+2,ka-1)
            x = abs(p) + abs(q) + abs(r)
            IF ( x.EQ.zero ) GO TO 520
            p = p / x
            q = q / x
            r = r / x
            s = sign( sqrt( p*p + q*q + r*r ), p )
            aad(ka,ka-1) = - s * x
         END IF
         p = p + s
         x = p / s
         y = q / s
         z = r / s
         q = q / p
         r = r / p
!                                                    ** ROW MODIFICATION
         DO 490 j = ka, m
            p = aad(ka,j) + q * aad(ka+1,j)
            IF ( notlas ) THEN
               p = p + r * aad(ka+2,j)
               aad(ka+2,j) = aad(ka+2,j) - p * z
            END IF
            aad(ka+1,j) = aad(ka+1,j) - p * y
            aad(ka,j)   = aad(ka,j)   - p * x
490      CONTINUE
!                                                 ** COLUMN MODIFICATION
         DO 500 ii = 1, min0(n,ka+3)
            p = x * aad(ii,ka) + y * aad(ii,ka+1)
            IF ( notlas ) THEN
               p = p + z * aad(ii,ka+2)
               aad(ii,ka+2) = aad(ii,ka+2) - p * r
            END IF
            aad(ii,ka+1) = aad(ii,ka+1) - p * q
            aad(ii,ka)   = aad(ii,ka) - p
500      CONTINUE
!                                          ** ACCUMULATE TRANSFORMATIONS
         DO 510 ii = l, k
            p = x * evecd(ii,ka) + y * evecd(ii,ka+1)
            IF ( notlas ) THEN
               p = p + z * evecd(ii,ka+2)
               evecd(ii,ka+2) = evecd(ii,ka+2) - p * r
            END IF
            evecd(ii,ka+1) = evecd(ii,ka+1) - p * q
            evecd(ii,ka)   = evecd(ii,ka) - p
510      CONTINUE

520   CONTINUE
      GO TO 390
!                     ** ALL EVALS FOUND, NOW BACKSUBSTITUTE REAL VECTOR
530   CONTINUE
      IF ( rnorm.NE.zero ) THEN
         DO 560  n = m, 1, -1
            n2 = n
            aad(n,n) = one
            DO 550  i = n-1, 1, -1
               w = aad(i,i) - evald(n)
               IF ( w.EQ.zero ) w = tol * rnorm
               r = aad(i,n)
               DO 540 j = n2, n-1
                  r = r + aad(i,j) * aad(j,n)
540            CONTINUE
               aad(i,n) = - r / w
               n2 = i
550         CONTINUE
560      CONTINUE
!                      ** END BACKSUBSTITUTION VECTORS OF ISOLATED EVALS

         DO 580 i = 1, m
            IF ( i.LT.l .OR. i.GT.k ) THEN
               DO 570 j = i, m
                  evecd(i,j) = aad(i,j)
570            CONTINUE
            END IF
580      CONTINUE
!                                   ** MULTIPLY BY TRANSFORMATION MATRIX
         IF ( k.NE.0 ) THEN
            DO j = m, l, -1
               DO i = l, k
                  z = zero
                  DO 590 n = l, min0(j,k)
                     z = z + evecd(i,n) * aad(n,j)
590               CONTINUE
                  evecd(i,j) = z
               END DO
            END DO
         END IF

      END IF

      DO i = l, k
         DO j = 1, m
            evecd(i,j) = evecd(i,j) * wkd(i)
         END DO
      END DO
!                           ** INTERCHANGE ROWS IF PERMUTATIONS OCCURRED
      DO 640  i = l-1, 1, -1
         j = int(wkd(i))
         IF ( i.NE.j ) THEN
            DO 630 n = 1, m
               repl       = evecd(i,n)
               evecd(i,n) = evecd(j,n)
               evecd(j,n) = repl
630         CONTINUE
         END IF
640   CONTINUE

      DO 660 i = k+1, m
         j = int(wkd(i))
         IF ( i.NE.j ) THEN
            DO 650 n = 1, m
               repl       = evecd(i,n)
               evecd(i,n) = evecd(j,n)
               evecd(j,n) = repl
650         CONTINUE
         END IF
660   CONTINUE
!
  670 CONTINUE
!
!  normalizing eigenvector
      DO j = 1, m
        nor_factor = zero
        DO i = 1, m
        nor_factor = nor_factor + evecd(i,j)**2
        END DO
        nor_factor = sqrt(nor_factor)
        evecd(:,j) = evecd(:,j)/nor_factor
      END DO
      RETURN
   END SUBROUTINE  asymtx
!
!
     SUBROUTINE asymtx_tl( nZ, V, VAL, A_TL, V_TL, VAL_TL, Error_Status)
       IMPLICIT NONE
       INTEGER :: nz,error_status
       REAL(fp), DIMENSION(:,:) :: v, a_tl, v_tl
       REAL(fp), DIMENSION(:) :: val, val_tl
       REAL(fp), DIMENSION(nZ,nZ) :: v_int, a1_tl
       REAL(fp) :: b_tl
       INTEGER :: i,j
       CHARACTER(*), PARAMETER :: routine_name = 'ASYMTX_TL'
       CHARACTER(256) :: message
     !    
       v_int = matinv( v(1:nz,1:nz), error_status )
       IF( error_status /= success  ) THEN
         WRITE( message,'("Error in matrix inversion matinv( V(1:nZ,1:nZ), Error_Status ) ")' ) 
         CALL display_message( routine_name, &                                                    
                               trim(message), &                                                   
                               error_status )                                          
         RETURN                                                                                    
       END IF              
       ! part 1
       !!  A1_TL = ZERO
       a1_tl = matmul( v_int, matmul(a_tl, v) )
       
       ! part 2
       !!  A1_TL = ZERO
       !!  V_TL = ZERO
       DO i = 1, nz
         val_tl(i) = a1_tl(i,i)
         DO j = 1, nz
          IF( i /= j ) THEN
            if( abs(val(j)-val(i)) > eigen_threshold ) then
              v_tl(i,j) = a1_tl(i,j)/(val(j)-val(i))
            else
              v_tl(i,j) = zero
            end if
          ELSE
            v_tl(i,j) = zero
          END IF
         END DO
       END DO
       
       ! part 3
       !!  A1_TL = ZERO
       a1_tl = matmul( v, v_tl) 
       !!  V_TL = ZERO   
       ! part 4
       DO i = 1, nz  
         b_tl = zero
         DO j = 1, nz
           b_tl = b_tl - v(j,i)*a1_tl(j,i)
         END DO
!
         DO j = 1, nz
           v_tl(j,i) = a1_tl(j,i) + v(j,i)*b_tl
         END DO
       END DO

       RETURN
       END SUBROUTINE asymtx_tl
!
!
     SUBROUTINE asymtx_ad( nZ, V, VAL, V_AD, VAL_AD, A_AD, Error_Status)
       IMPLICIT NONE
       INTEGER :: nz,error_status
       REAL(fp), DIMENSION(:,:) :: v, v_ad, a_ad
       REAL(fp), DIMENSION(:) :: val, val_ad
       INTEGER :: i,j,k
       REAL(fp), DIMENSION(nZ,nZ) :: v_int, a1_ad
       CHARACTER(*), PARAMETER :: routine_name = 'ASYMTX_AD'
       CHARACTER(256) :: message
     !
       a1_ad = zero
       DO i = 1, nz
         DO j = 1, nz
         DO k = 1, nz
          IF( k /= j ) THEN
            if( abs(val(j)-val(k)) > eigen_threshold ) &
             a1_ad(k,j) = a1_ad(k,j) + v(i,k)*v_ad(i,j)/(val(j)-val(k))
          END IF
         END DO
         END DO
         a1_ad(i,i) = a1_ad(i,i) + val_ad(i)
       END DO       
       v_int = matinv( v(1:nz,1:nz), error_status )
       IF( error_status /= success  ) THEN
         WRITE( message,'("Error in matrix inversion matinv( V(1:nZ,1:nZ), Error_Status ) ")' ) 
         CALL display_message( routine_name, &                                                    
                               trim(message), &                                                   
                               error_status )                                          
         RETURN                                                                                    
       END IF           
       
       a_ad = matmul( transpose(v_int), matmul(a1_ad, transpose(v) ) )     
       RETURN
     END SUBROUTINE asymtx_ad
!
!
     SUBROUTINE asymtx2_ad( COS_Angle,n_streams, nZ, V, VAL, V_AD, VAL_AD, A_AD, Error_Status)
       IMPLICIT NONE
       INTEGER :: nz,error_status
       REAL(fp), DIMENSION(:,:) :: v, v_ad, a_ad
       REAL(fp), DIMENSION(:) :: val, val_ad
       INTEGER :: i,j,n_streams,n
       REAL(fp) :: cos_angle,d0,d1
       REAL(fp), DIMENSION(nZ,nZ) :: v_int, a1_ad
       REAL(fp) :: b_ad
     !
       n = 0
       IF( nz /= n_streams ) THEN
       ! additional stream for satellite viewing angle is used.
          d0 = one/cos_angle
          d0 = d0 * d0
          d1 = 100000.0_fp
          DO i = 1, nz
            if( abs( d0 - val(i) ) < d1 .and. abs(v(2,i)) < eigen_threshold ) then
              n = i
              d1 = abs( d0 - val(i) )
            end if
          END DO
       END IF 
       IF( n /= 0 ) val_ad(n) = zero 

     ! using normoalization condition
       b_ad = zero
       !! A1_AD = ZERO
       ! part 4
       DO i = nz, 1, -1 
         DO j = nz, 1, -1
           b_ad = b_ad + v(j,i)*v_ad(j,i)
           a1_ad(j,i) = v_ad(j,i)
           !! A1_AD(j,i) = A1_AD(j,i) + V_AD(j,i)
         END DO

         DO j = nz, 1, -1
           a1_ad(j,i) = a1_ad(j,i) - v(j,i)*b_ad
         END DO
         b_ad = zero
       END DO

       ! part 3
       !! V_AD = ZERO
       !! V_AD = V_AD + matmul( transpose(V), A1_AD)        
       v_ad = matmul( transpose(v), a1_ad) 

       ! part 2
       a1_ad = zero
       DO i = nz, 1, -1
         DO j = nz, 1, -1
          IF( i /= j .and. j /= n ) THEN
            if( abs(val(j)-val(i)) > eigen_threshold ) then
              a1_ad(i,j) = v_ad(i,j)/(val(j)-val(i))
            else
              v_ad(i,j) = zero
            end if
          ELSE
            v_ad(i,j) = zero
          END IF
         END DO
         a1_ad(i,i) = a1_ad(i,i) + val_ad(i)
       END DO

       v_int = matinv( v(1:nz,1:nz), error_status )
       ! part 1
       !! A_AD = ZERO
       a_ad = matmul( transpose(v_int), matmul(a1_ad, transpose(v) ) )

       RETURN
     END SUBROUTINE asymtx2_ad
!
!
 END MODULE crtm_utility 
!