<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
!
! 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
!
<A NAME='CRTM_UTILITY'><A href='../../html_code/crtm/CRTM_Utility.f90.html#CRTM_UTILITY' TARGET='top_target'><IMG SRC="../../gif/bar_purple.gif" border=0></A>
MODULE CRTM_UTILITY 4,7
! -----------------
! 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
! -------------------
<A NAME='LEGENDRE'><A href='../../html_code/crtm/CRTM_Utility.f90.html#LEGENDRE' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
INTERFACE Legendre
MODULE PROCEDURE
MODULE PROCEDURE
END INTERFACE
! -----------------
! Module parameters
! -----------------
! Version Id for the module
CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
'$Id: CRTM_Utility.f90 29405 2013-06-20 20:19:52Z paul.vandelst@noaa.gov $'
! Numerical small threshold value in Eigensystem
REAL(fp), PARAMETER :: EIGEN_THRESHOLD = 1.0e-20_fp
CONTAINS
<A NAME='DOUBLE_GAUSS_QUADRATURE'><A href='../../html_code/crtm/CRTM_Utility.f90.html#DOUBLE_GAUSS_QUADRATURE' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE DOUBLE_GAUSS_QUADRATURE(NUM, ABSCISSAS, WEIGHTS) 1
! 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
!
!
<A NAME='MATINV'><A href='../../html_code/crtm/CRTM_Utility.f90.html#MATINV' TARGET='top_target'><IMG SRC="../../gif/bar_green.gif" border=0></A>
FUNCTION matinv(A, Error_Status) 11,1
! --------------------------------------------------------------------
! 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
!
!
<A NAME='LEGENDRE_SCALAR'><A href='../../html_code/crtm/CRTM_Utility.f90.html#LEGENDRE_SCALAR' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE Legendre_Scalar(MOA,ANG,PL) 1
! 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
!
!
<A NAME='LEGENDRE_RANK1'><A href='../../html_code/crtm/CRTM_Utility.f90.html#LEGENDRE_RANK1' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE Legendre_Rank1(MOA,NU,ANG,PL) 1
! 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
!
!
<A NAME='LEGENDRE_M'><A href='../../html_code/crtm/CRTM_Utility.f90.html#LEGENDRE_M' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE Legendre_M(MF,N,U,PL) 2
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
!
!
<A NAME='GAMMA2'><A href='../../html_code/crtm/CRTM_Utility.f90.html#GAMMA2' TARGET='top_target'><IMG SRC="../../gif/bar_green.gif" border=0></A>
FUNCTION gamma2(N),1
! 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
!
!
<A NAME='ASYMTX'><A href='../../html_code/crtm/CRTM_Utility.f90.html#ASYMTX' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE ASYMTX( AAD, M, IA, IEVEC, & 1
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
!
!
<A NAME='ASYMTX_TL'><A href='../../html_code/crtm/CRTM_Utility.f90.html#ASYMTX_TL' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE ASYMTX_TL( nZ, V, VAL, A_TL, V_TL, VAL_TL, Error_Status) 1,2
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
!
!
<A NAME='ASYMTX_AD'><A href='../../html_code/crtm/CRTM_Utility.f90.html#ASYMTX_AD' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE ASYMTX_AD( nZ, V, VAL, V_AD, VAL_AD, A_AD, Error_Status) 1,2
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
!
!
<A NAME='ASYMTX2_AD'><A href='../../html_code/crtm/CRTM_Utility.f90.html#ASYMTX2_AD' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
SUBROUTINE ASYMTX2_AD( COS_Angle,n_streams, nZ, V, VAL, V_AD, VAL_AD, A_AD, Error_Status),1
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
!