SUBROUTINE CROUT(A,N,INDEX,B) C THIS VERSION OF MATRIX DOESN'T USE THE LAPACK DGET ROUTINES C IT CAME FROM KARL HOPPEL C INCLUDE 'PARAMETERS.INC' DOUBLE PRECISION A(MB,MB) , B(MB) DIMENSION INDEX(MB) C THIS GENERIC SUBROUTINE IS USED IN CONJUCTION WITH SUBROUTINE LU TO SOLVE A SYSTEM C OF EQUATIONS. THE SOLUTION VECTOR OVERWRITES THE RIGHT HAND SIDE VECTOR B. C THIS ROUTINE ACCOUNTS FOR THE POSSIBILITY THAT B WILL BEGIN WITH MANY ZERO C ELEMENTS THUS REDUCING THE NUMBER OF CALCULATIONS WHICH NEED TO BE DONE. II = 0 DO I = 1 , N L = INDEX(I) SUM = B(L) B(L) = B(I) IF (II .NE. 0)THEN DO J = II , I-1 SUM = SUM - A(I,J)*B(J) ENDDO ELSE IF (SUM .NE. 0.0) II = I ENDIF B(I) = SUM ENDDO C NOW GET THE SOLUTION VECTOR DO I = N , 1 , -1 SUM = B(I) DO J = I+1 , N SUM = SUM - A(I,J)*B(J) ENDDO B(I) = SUM/A(I,I) ENDDO RETURN END C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- SUBROUTINE LU(A,N,INDEX) INCLUDE 'PARAMETERS.INC' PARAMETER (EPSILON = 1.0e-20,NMAX = 300) DOUBLE PRECISION SCALE(NMAX) , A(MB,MB) DIMENSION INDEX(MB) C THIS GENERIC ROUTINE COMPUTES THE L AND U TRIANGULAR MATRICES EQUIVALENT TO THE C A MATRIX, SUCH THAT LU = A. THESE MATRICES ARE RETURNED IN THE SPACE OF A. THE C U MATRIX HAS ONES ON ITS DIAGONAL. THAT IS, THIS SUBROUTINE RETURNS C L(1,1) = A(1,1) , U(1,1) = 1 , L(2,1) = A(2,1) , U(1,2) = A(1,2) , ETC. C PARTIAL PIVOTING IS USED TO GIVE MAXIMUM VALUED ELEMENTS ON THE DIAGONAL OF L. C THE ORDER OF THE ROWS AFTER PIVOTING IS RETURNED IN THE VECTOR INDEX. C THIS SHOULD BE USED TO REORDER THE RIGHT HAND SIDE VECTOR BEFORE SOLVING Ax = b. C FIND LARGEST ELEMENT IN EACH ROW DO I = 1 , N AAMAX = 0.0D0 DO J = 1 , N c AAMAX = AMAX1(AAMAX,ABS(A(I,J))) AAMAX = dMAX1(AAMAX,dABS(A(I,J))) ENDDO IF (AAMAX .EQ. 0.0) THEN PAUSE '**SINGULAR MATRIX** !!' ENDIF SCALE(I) = 1.0D0/AAMAX ! SCALING FOR ROW I ENDDO DO J = 1 , N DO I = 1 , J-1 SUM = A(I,J) DO K = 1 , I-1 SUM = SUM - A(I,K)*A(K,J) ENDDO A(I,J) = SUM ENDDO AAMAX = 0.0 C COMPUTE A COLUMN OF L DO I = J , N SUM = A(I,J) DO K = 1 , J-1 SUM = SUM - A(I,K)*A(K,J) ENDDO A(I,J) = SUM TEMP = SCALE(I)*ABS(SUM) IF (TEMP .GE. AAMAX) THEN IMAX = I AAMAX = TEMP ENDIF ENDDO C INTERCHANGE ROWS IF NEED BE IF (J .NE. IMAX)THEN DO K = 1 , N TEMP = A(IMAX,K) A(IMAX,K) = A(J,K) A(J,K) = TEMP ENDDO SCALE(IMAX) = SCALE(J) ENDIF INDEX(J) = IMAX IF(A(J,J) .EQ. 0.0) A(J,J) = EPSILON C COMPUTE A ROW OF U IF(J .NE. N)THEN TEMP = 1.0D0/A(J,J) DO I = J+1 , N A(I,J) = A(I,J)*TEMP ENDDO ENDIF ENDDO RETURN END C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- SUBROUTINE THOMAS INCLUDE 'PARAMETERS.INC' COMMON & /ONED7/ A(MB,NALT) , C(MB,NALT) & /ONED9/ DCONC(MB,NALT) & /ONED13/ ADELTA(MB,MB,NALT) , CDELTA(MB,MB,NALT) & /ONED17/ NTRANS & /ONED30/ B(MB,NALT) & /ONED31/ BDELTA(MB,MB,NALT) & /ONED39/ RESID(MB,NALT) & /ONED40/ MTRANS(MB) DOUBLE PRECISION GAMMA(MB,MB,NALT) , BTEMP(MB,MB) , BORIG(MB,MB) , & R(MB) , RTEMP(MB) , TEMP(MB,NALT) , SUMR DIMENSION INDEX(MB) C THIS MODULE USES THE THOMAS ALGORITHM TO SOLVE A SET OF EQUATIONS WHICH ARE C TRIDIAGONAL IN THE ALTITUDE VARIABLE. THE ELEMENTS OF THE TRIDIAGONAL SYSTEM C ARE THEMSELVES MATRICES. THUS, THE MULTIPLICATIONS AND DIVISIONS WHICH ARE FOUND C IN THE THOMAS ALGORITHM MUST BE INTERPRETED AS MULTRIX MULTIPLICATION AND DIVISION C (INVERSION). A RESIDUAL CORRECTION PROCEEDURE IS APPLIED TO THE COMPUTED C SOLUTION IN ORDER TO INCREASE THE ACCURACY. C THE LU DECOMPOSITION IS USED TO PERFORM THE MATRIX MULTIPLICATIONS AND C MATRIX SOLVES WHICH APPEAR IN THE THOMAS ALGORITHM. C--------------------- LIST OF RELEVANT VARIABLES ------------------------------- C BTEMP = BDELTA AT A SPECIFIC ALTITUDE LEVEL. IS PASSED TO SUBROUTINE LU AND C THE LU DECOMPOSITION IS PERFORMED ON IT. C BORIG = BTEMP BEFORE LU DECOMPOSITION IS PERFORMED. IS USED IN THE RESIDUAL C CORRECTION PROCEDURE. IT IS NEEDED BECAUSE SUBROUTINE LU CHANGES C THE FORM OF BTEMP. C RTEMP = RESID AT A SPECIFIC ALTITUDE LEVEL. IS THE RIGHT HAND SIDE VECTOR C IN THE CALL TO CROUT. THE SOLUTION OF THE MATRIX EQUATION IS C RETURNED IN RTEMP. C R = RESID - BDELTA*DCONC (I.E. RESIDUAL) THIS IS NOT ZERO DUE TO C THE INACCURACIES WHICH CREEP IN WHEN SOLVING BDELTA*DCONC = RESID C R IS USED AS A RESIDUAL CORRECTION IN ORDER TO INCREASE THE C ACCURACY OF THE COMPUTED SOLUTION. C-------------------------------------------------------------------------------- C SET UP MATRIX ELEMENTS C BDELTA(1) AND BDELTA(NALT) ARE COMPUTED IN SUBROUTINE BNDRY C BDELTA(K) , K = 2 , NALT-1 IS COMPUTED IN SUBROUTINE BMATRIX DO K = 2 , NALT-1 DO M2 = 1 , NTRANS DO M1 = 1 , NTRANS ADELTA(M1,M2,K) = 0.0D0 CDELTA(M1,M2,K) = 0.0D0 ! CDELTA(M1,M2,1) IS SET IN SUBROUTINE BNDRY ENDDO ADELTA(M2,M2,K) = A(M2,K) CDELTA(M2,M2,K) = C(M2,K) ENDDO ENDDO DO M2 = 1 , NTRANS CDELTA(M2,M2,NALT) = C(M2,NALT) ! ADELTA(NALT) IS SET IN BNDRY ENDDO C NOW PERFORM THOMAS ALGORITHM C LEVEL 1 ELEMENTS-----: BTEMP = BDELTA(1) ; DCONC(1) = (BTEMP^-1)*DCONC(1) DO M2 = 1 , NTRANS DO M1 = 1 , NTRANS BTEMP(M1,M2) = BDELTA(M1,M2,1) ! TRANSFER BDELTA TO AN NTRANSxNTRANS ! ARRAY BTEMP BORIG(M1,M2) = BTEMP(M1,M2) ENDDO RTEMP(M2) = RESID(M2,1) ENDDO CALL LU(BTEMP,NTRANS,INDEX) CALL CROUT(BTEMP,NTRANS,INDEX,RTEMP) ! SOLVE FOR DCONC DO M = 1 , NTRANS DCONC(MTRANS(M),1) = RTEMP(M) ENDDO C---------------------- APPLY RESIDUAL CORRECTION --------------------------------- DO M = 1 , NTRANS SUMR = -RESID(M,1) DO M1 = 1 , NTRANS SUMR = SUMR + DBLE(BDELTA(M,M1,1))*DBLE(DCONC(MTRANS(M1),1)) ENDDO RTEMP(M) = SUMR ENDDO CALL CROUT(BTEMP,NTRANS,INDEX,RTEMP) ! SOLVE FOR RESIDUAL DO M = 1 , NTRANS MVARY = MTRANS(M) DCONC(MVARY,1) = DCONC(MVARY,1) - RTEMP(M) ! MORE ACCURATE VALUE OF DCONC(M,1) ENDDO C---------------------------------------------------------------------------------- C REMAINING LEVELS. GENERAL ALGORITHM IS: C GAMMA(K) = (BTEMP^-1)*CDELTA(K-1) C BETA = BDELTA(K) - ADELTA(K)*GAMMA(K) C DCONC(K) = (BTEMP^-1)*(RESID(K) - ADELTA(K)*DCONC(K-1)) DO K = 2 , NALT DO M1 = 1 , NTRANS DO M2 = 1 , NTRANS RTEMP(M2) = CDELTA(M2,M1,K-1) ! WE WANT TO SOLVE FOR GAMMA (A MATRIX - ! NOT A VECTOR) SO THE RTEMP IS A COLUMN ! OF CDELTA ENDDO CALL CROUT(BTEMP,NTRANS,INDEX,RTEMP) ! SOLVE FOR RTEMP (LU DECOMPOSITION ! OF BTEMP HAS ALREADY BEEN DONE) C--------------------------- RESIDUAL CORRECTION ------------------------------------ DO M3 = 1 , NTRANS SUMR = -CDELTA(M3,M1,K-1) DO M4 = 1 , NTRANS SUMR = SUMR + DBLE(BORIG(M3,M4))*DBLE(RTEMP(M4)) ENDDO R(M3) = SUMR ENDDO CALL CROUT(BTEMP,NTRANS,INDEX,R) ! SOLVE FOR RESIDUAL DO M3 = 1 , NTRANS RTEMP(M3) = RTEMP(M3) - R(M3) ENDDO DO M2 = 1 , NTRANS GAMMA(M2,M1,K) = RTEMP(M2) ! IMPROVED VALUE OF GAMMA ENDDO C-------------------------------------------------------------------------------------- ENDDO DO M1 = 1 , NTRANS SUM = 0.0D0 DO M2 = 1 , NTRANS MVARY2 = MTRANS(M2) SUMG = 0.0D0 DO M3 = 1 , NTRANS SUMG = SUMG + ADELTA(M1,M3,K)*GAMMA(M3,M2,K) ! ADELTA*GAMMA ENDDO BTEMP(M1,M2) = BDELTA(M1,M2,K) - SUMG BORIG(M1,M2) = BTEMP(M1,M2) SUM = SUM + ADELTA(M1,M2,K)*DCONC(MVARY2,K-1) ! ADELTA*DCONC ENDDO TEMP(M1,K) = RESID(M1,K) - SUM ! DCONC(K) - ADELTA(K)*DCONC(K-1) ENDDO C (BTEMP^-1)*(RESID(K)-ADELTA(K)*DCONC(K-1)) DO M = 1 , NTRANS RTEMP(M) = TEMP(M,K) ENDDO CALL LU(BTEMP,NTRANS,INDEX) CALL CROUT(BTEMP,NTRANS,INDEX,RTEMP) ! REURNS DCONC IN PLACE OF RTEMP C--------------------------------- RESIDUAL CORRECTION ----------------------------------- DO M = 1 , NTRANS SUMR = -TEMP(M,K) DO M1 = 1 , NTRANS SUMR = SUMR + DBLE(BORIG(M,M1))*DBLE(RTEMP(M1)) ENDDO R(M) = SUMR ENDDO CALL CROUT(BTEMP,NTRANS,INDEX,R) ! SOLVE FOR RESIDUAL DO M = 1 , NTRANS RTEMP(M) = RTEMP(M) - R(M) ENDDO DO M = 1 , NTRANS DCONC(MTRANS(M),K) = RTEMP(M) ENDDO C------------------------------------------------------------------------------------------ ENDDO C NOW BACKSOLVE. ALGORITHM IS: DCONC(K) = DCONC(K) - GAMMA(K+1)*DCONC(K+1) DO K = NALT-1 , 1 , -1 DO M1 = 1 , NTRANS SUM = 0.0D0 DO M2 = 1 , NTRANS SUM = SUM + GAMMA(M1,M2,K+1)*DCONC(MTRANS(M2),K+1) ENDDO MVARY1 = MTRANS(M1) DCONC(MVARY1,K) = DCONC(MVARY1,K) - SUM ENDDO ENDDO RETURN END C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- C--------------------------------------------------------------------------------- SUBROUTINE JACOB INCLUDE 'PARAMETERS.INC' COMMON & /ONED16/ DNET(MB,MB,NALT) & /ONED17/ NTRANS & /ONED18/ NFORM(MB) , NDEST(MB) , & IFORM(IB,MB) , IDEST(IB,MB) & /ONED19/ KFORM(IB,MB) , KDEST(IB,MB) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /ONED40/ MTRANS(MB) & /ONED42/ NMOL & /BLOCK3/ NOREAC(IB) ! FROM CHEM DOUBLE PRECISION DF(MB,MB,NALT) , DD(MB,MB,NALT) C CALCULATES JACOBIAN OF PRODUCTION - LOSS VECTOR DO K = 1 , NALT DO M = 1 , NMOL DO M1 = 1 , NTRANS DF(M1,M,K) = 0.0D0 DD(M1,M,K) = 0.0D0 ENDDO ENDDO ENDDO DO K = 1 , NALT DO M = 1 , NTRANS MV = MTRANS(M) C CALCULATE PARTIAL DERIVATIVE CONTRIBUTIONS FROM FORMATION REACTIONS. NF = NFORM(MV) ! NUMBER OF REACTIONS IN WHICH MOLECULE MV IS FORMED DO I = 1 , NF IFRM = IFORM(I,MV) ! REACTION # CORRESPONDING TO THE ITH ! REACTION IN WHICH MV IS FORMED NRCT = NOREAC(IFRM) DO J = 1 , NRCT M2 = MOL(J,IFRM) DF(M,M2,K) = DF(M,M2,K) + KOEFF(J,IFRM)*KFORM(IFRM,MV)* & DERIV(IFRM,M2,K) ENDDO ENDDO C CALCULATE PARTIAL DERIVATIVE CONTRIBUTIONS FROM DESTRUCTION REACTIONS. ND = NDEST(MV) ! NUMBER OF REACTIONS IN WHICH MOLECULE MV IS DESTROYED DO I = 1 , ND IDST = IDEST(I,MV) ! REACTION # CORRESPONDING TO THE ITH ! REACTION IN WHICH MV IS DESTROYED NRCT = NOREAC(IDST) DO J = 1 , NRCT M2 = MOL(J,IDST) DD(M,M2,K) = DD(M,M2,K) + KOEFF(J,IDST)*KDEST(IDST,MV)* & DERIV(IDST,M2,K) ENDDO ENDDO C CALCULATE NET PARTIAL DERIVATIVES. DO M2 = 1 , NMOL DNET(M,M2,K) = DF(M,M2,K) - DD(M,M2,K) ENDDO ENDDO ENDDO RETURN END C-------------------------------------------------------------------------------- C-------------------------------------------------------------------------------- C-------------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DERIV(IREACT,M,K) INCLUDE 'PARAMETERS.INC' COMMON & /BLOCK7/ RATE(IB,NALT) & /ONED20/ RCONST(IB,NALT) & /BLOCK18/ CONC(MB,NALT) & /ONED27/ MOL(6,IB) , KOEFF(6,IB) & /BLOCK3/ NOREAC(IB) ! FROM CHEM C CALCULATES THE PARTIAL DERIVATIVE OF THE PRODUCTION/LOSS RATE. SINCE THE C GENERAL FORM FOR THE PROD/LOSS RATE IS GIVEN BY C K*([C1]**K1)*([C2]**K2)*([C3]**K3) THE DERIVATIVE CAN BE EASILY COMPUTED IN C TERMS OF THE PRODUCTION/LOSS RATE AS LONG AS [C] IS NOT ZERO. IF [C] IS C ZERO, THEN THE RESULTING EXPRESSION WOULD BE OF THE FORM ZERO/ZERO. IN THAT CASE C THE DERIVATIVE MUST BE CALCULATED EXPLICITLY USING THE [C]'s AND K's. IF (CONC(M,K) .GT. TINY) THEN DERIV = RATE(IREACT,K)/CONC(M,K) ELSE DERIV = RCONST(IREACT,K) NRCT = NOREAC(IREACT) DO J = 1 , NRCT IF (M .NE. MOL(J,IREACT)) THEN DERIV = DERIV*(CONC(MOL(J,IREACT),K)**KOEFF(J,IREACT)) ELSE IF( KOEFF(J,IREACT) .GT. 1 ) THEN DERIV = DERIV* & (CONC(MOL(J,IREACT),K)**(KOEFF(J,IREACT)-1)) ENDIF ENDIF ENDDO ENDIF RETURN END