c subroutine ohrad(kmbot,kmtop,mvl,nkm,TA,SO2,SO,SN2,SH,SO3,SHO2, + SOH0,SM,ideltav,dohv,bohv,iwatts,spval,iprint) c c B. Foster 7/19/95: c c Input: c kmbot,kmtop = bottom, top altitudes (km) (30->120km) c nkm = number of altitudes (120-30+1km = 91) c mvl = number of vibrational levels including ground (10) c TA = neutral temp (k) (time-gcm) c SO2,SO,SN2,SH,SO3,SHO2,SOHO,SM = number densities from time-gcm c o2, o, n2, h, o3, ho2, oh ,rho c ideltav = iteration flag (max 6) c iprint = 1 -> write output to stdout c 9/11/96: iwatts=1 -> bohv units in watts/cm3, else photons/cm3-sec c Output: c dohv(mvl,nkm) = 10 vib levels for each altitude (cm-3) c bohv(60,nkm) = brightness (kilo Raleighs) at each band (see aoh) c C ... CALCULATES C 1) THE NUMBER DENSITY FOR THE FIRST 9 VIBRATIONAL LEVELS C ... OF THE OH. c ..................................................................... c c /~~~~~\ c { @ @ } c +=======================oooo=====o=====oooo=========================+ c * * c * author: * c * U. B. Makhlouf * c * SDL/Stewart Radiance Laboratory * c * 139 the Great Rd. * c * Bedford, MA 01730 * c * * c * e-mail: makhlouf@pldac.plh.af.mil * c * tel: (617) 377-4203 * c * * c ********************************************************************* DOUBLE PRECISION DET INTEGER DELTAV real dohv(mvl,nkm),bohv(60,nkm) C C ... NVL: IS # OF VIBRATIONAL LEVELS CALCULATED PLUS GROUND parameter(nvl=10) DIMENSION SH(nkm),SO3(nkm),SO(nkm),SHO2(nkm),SN2(nkm),SM(nkm) DIMENSION AOH(nvl,6),ENERGY(nvl),ENERGYD(nvl),ALT(nkm) DIMENSION SOH0(nkm) DIMENSION SO2(nkm),TA(nkm),FVL(nvl),PROHVS(nkm) DIMENSION PROHV(nkm),FVLS(14),FAC1(nvl),FAC0(nvl) DIMENSION Q1O2(nvl),Q1N2(nvl),CLO(nvl),RQ1O2(nvl),RQ1N2(nvl) DIMENSION AA(nvl,nvl),ALU(nvl,nvl),B(nvl),IPS(nvl),X(nvl) c DIMENSION A2(nvl,nvl),X2(nvl),FAC9(nvl) DIMENSION A2(nvl,nvl),FAC9(nvl) DIMENSION AAO2(nvl),EEO2(nvl),AAN2(nvl),EEN2(nvl) DIMENSION AAO(nvl),EEO(nvl),REE(nvl),FRAQ(nvl,nkm) c DIMENSION RAAO2(nvl),RAAN2(nvl),DOHV2(nvl,nkm) DIMENSION RAAO2(nvl),RAAN2(nvl) DIMENSION RK1(nkm),RK2(nkm),RK3(nkm),RK4(nkm),RK5(nkm) DIMENSION RK6(nkm) DIMENSION FQV(nvl,6),XMQO2(nvl,nvl),RMQO2(nvl,nvl),SQV(nvl), + KIM(nvl) C c COMMON/A1/ AOH,ENERGYD,ILA,IHA,NTG c COMMON/B1/ PI,BK,C,HP,DELTAV C C ... FRACTIONAL VIBRATIONAL LEVEL TAKEN FROM 'KLENERMAN & SMITH' DATA FVL/6*0.0,0.08,0.17,0.27,0.48/ C C ... FVLS IS FVL FOR THE SECONDARY REACTION C......... KAYE ............................ DATA FVLS/0.52,0.34,0.13,.01,10*0.0/ C DATA FAC1/0.0,9*1.0/ DATA FAC0/1.0,9*0.0/ DATA FAC9/9*1.0,0.0/ c c B. Foster: use data statements instead of read(5) (see end of this file) c data energy + /00000.00, 03568.50, 06971.37, 10210.57, 13287.19, 16201.33, + 18951.90, 21536.25, 23949.83, 26185.79/ c c aoh are Einstein coeffs as follows: c c 0 1-0 2-1 3-2 4-3 5-4 6-5 7-6 8-7 9-8 c 0 0 2-0 3-1 4-2 5-3 6-4 7-5 8-6 9-7 c 0 0 0 3-0 4-1 5-2 6-3 7-4 8-5 9-6 c 0 0 0 0 4-0 5-1 6-2 7-3 8-4 9-5 c 0 0 0 0 0 5-0 6-1 7-2 8-3 9-4 c 0 0 0 0 0 0 6-0 7-1 8-2 9-3 c data aoh + /0.000,22.74,30.43,28.12,20.30,11.05,4.000,2.340,8.600,23.72, + 0.000,0.000,15.42,40.33,69.77,99.42,125.6,145.1,154.3,148.9, + 0.000,0.000,0.000,2.032,7.191,15.88,27.94,42.91,59.98,78.64, + 0.000,0.000,0.000,0.000,0.299,1.315,3.479,7.165,12.68,19.94, + 0.000,0.000,0.000,0.000,0.000,0.051,0.274,0.847,2.007,4.053, + 0.000,0.000,0.000,0.000,0.000,0.000,0.010,0.063,0.230,0.620/ c data aao2 /0.0000000, 1.300E-13, 2.700E-13, 5.200E-13, 8.800E-13, + 1.700E-12, 3.000E-12, 5.420E-12, 9.810E-12, 1.700E-11/ data aan2 /0.0000000, 5.757E-15, 1.000E-14, 1.737E-14, 3.017E-14, + 5.241E-14, 9.103E-14, 1.581E-13, 2.746E-13, 4.770E-13/ data aao /3.900E-11, 10.50E-11, 2.500E-10, 2.500E-10, 2.500E-10, + 2.500E-10, 2.500E-10, 2.500E-10, 2.500E-10, 2.500E-10/ data raao2 /1.300E-13, 2.700E-13, 5.200E-13, 8.800E-13, 1.700E-12, + 3.000E-12, 5.420E-12, 9.810E-12, 1.700E-11, 0.0000000/ data raan2 /5.757E-15, 1.000E-14, 1.737E-14, 3.017E-14, 5.241E-14, + 9.103E-14, 1.581E-13, 2.746E-13, 4.770E-13, 0.0000000/ c ree(:) = 0. eeo2(:) = 0. een2(:) = 0. eeo(:) = 0. aa(:,:) = 0. alu(:,:) = 0. b(:) = 0. xmqo2(:,:) = 0. rmqo2(:,:) = 0. bohv(:,:) = 0. c c write(6,"('ohrad: ta=',/(6e12.4))") ta c write(6,"('ohrad: so2=',/(6e12.4))") so2 c write(6,"('ohrad: so=',/(6e12.4))") so c write(6,"('ohrad: sn2=',/(6e12.4))") sn2 c write(6,"('ohrad: sh=',/(6e12.4))") sh c write(6,"('ohrad: so3=',/(6e12.4))") so3 c write(6,"('ohrad: sho2=',/(6e12.4))") sho2 c write(6,"('ohrad: soh0=',/(6e12.4))") soh0 c write(6,"('ohrad: sm=',/(6e12.4))") sm c C ... SPECIES # AS ARRANGED IN THE MATRICES C ... 1-10 OH(V=0-9) C ... 11 O3 C ... 12 O C ... 13 H C ... 14 HO2 C C ... ILA: IS THE LOWER ALTITUDE, IHA: IS THE HIGHER ALTITUDE ILA=30 IHA=120 NPTS=(IHA-ILA)+1 c c Check input arg list: c if (nvl.ne.mvl.or.ila.ne.kmbot.or.iha.ne.kmtop.or.npts.ne.nkm) + then write(6,"(/72('-')/,'>>> ohrad: kmbot, kmtop, nkm, mvl = ',4i4, + /' but should be = ',4i4,/72('-')/)") + kmbot,kmtop,nkm,mvl, ila,iha,npts,nvl stop 'ohrad' endif if (ideltav.lt.1.or.ideltav.gt.6) then write(6,"(/72('-'),'>>> ohrad: bad ideltav = ',i4, + ' must be 1 -> 6',/72('-')/)") ideltav stop 'ohrad' endif deltav = ideltav c C ... NTG # OF TANGENT HEIGHTS NTG=41 INDEX=0 C C ... If XMMQ=0 or 1 a '1' invokes a multi quanta quenching by O2 XMMQ=1.0 C RORIG=1.-XMMQ CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PI=3.1415927 BK=1.3805E-16 C=2.997925E+10 HP=6.6256E-34 C2K=1.438818 C......C2K UNITS IS (DEGREES KELVIN.CM) C......IO...IO...IO...IO...IO...IO...IO...IO...IO...IO...IO... C c OPEN(3,FILE=' INSERT FILE1', c 1 status='readonly') c REWIND(3) c READ(3,11) (SO(I),SH(I),SO3(I),SHO2(I),SOH0(I),I=1,NPTS) c11 FORMAT(5E10.3) C c OPEN(4,FILE='INSERT FILE2',status='readonly') c REWIND(4) c DO 10 I=1,NPTS c J=NPTS+1-I c10 READ(4,12) TA(J),SN2(J),SO2(J),SM(J) c12 FORMAT(4E12.4) c210 CONTINUE C c OPEN(5,FILE='oh_qrates.lowe',status='readonly') c REWIND(5) c DO 18 KV=NVL,1,-1 c READ(5,*) (AOH(KV,J),J=1,6),ENERGY(KV) c18 CONTINUE c READ(5,*) (AAO2(J),AAN2(J),AAO(J),RAAO2(J),RAAN2(J), J=NVL,1,-1) C c OPEN(6,FILE=INSERT FILE3) C c PRINT 7 c7 FORMAT('ENTER DELTA v (MAX 6): ',$) c READ(*,*) DELTAV C C......IO...IO...IO...IO...IO...IO...IO...IO...IO...IO...IO... C C**************************************************************** C ... AOH'S ARE THE EINSTEIN COEFFICIENTS, 1 FOR SINGLE QUANTA, 2 FOR C ... DOUBLE QUANTA,... ETC. C ... THE AOH ARE READ FROM THE OH_QRATES FILE IN THE FOLLOWING ORDER C ... FROM THE FIRST 9x6 MATRIX: C ... 9-8 9-7 9-6 9-5 9-4 9-3 C ... 8-7 8-6 8-5 8-4 8-3 8-2 C ... 7-6 7-5 7-4 7-3 7-2 7-1 C ... 6-5 6-4 6-3 6-2 6-1 6-0 C ... 5-4 5-3 5-2 5-1 5-0 C ... 4-3 4-2 4-1 4-0 C ... 3-2 3-1 3-0 C ... 2-1 2-0 C ... 1-0 C ... FOR THE QUENCHING RATES THE SHARC FORMAT HAS BEEN ADOPTED, WHERE C ... THE QUENCHING RATES ARE WRITTEN IN THE FORM: C ... A*T**B*EXP(-C/T**(1/3))*EXP(-E/T) C ... FOR OUR CASE THE B's AND THE C's ARE ZERO's, SO WE ONLY READ THE A's C ... AND THE E's C ... AAO2 IS THE A FOR THE O2 QUENCHING C ... REE IS THE E FOR THE REVERSE QUENCHING REACTION C ... Q1O2 ARE THE QUENCHING RATE DUE TO O2, AND C ... RQ1O2 ARE THE REVERSE RATE OF THE REACTION DUE TO DETAILED BALANCE C *********************************************************************** REE(9)=0.0 EEO2(9)=0.0 EEN2(9)=0.0 EEO(9)=0.0 DO 23 KV=NVL-1,1,-1 EEO2(KV)=0.0 EEN2(KV)=0.0 EEO(KV)=0.0 REE(KV)=(ENERGY(KV+1)-ENERGY(KV))*C2K 23 CONTINUE DO 14 K=NVL,DELTAV+1,-1 ENERGYD(K)=ENERGY(K)-ENERGY(K-DELTAV) ENERGYD(K)=ENERGYD(K)*HP*C 14 CONTINUE if (iprint.gt.0) WRITE(6,44) 44 FORMAT('C',2X,'ALT.',3X,'OH(0)',5X,'OH(1)',5X,'OH(2)',5X,'OH(3)', 1 5X,'OH(4)',5X,'OH(5)',5X,'OH(6)',5X,'OH(7)',5X,'OH(8)',5X, 2 'OH(9)') aa(:,:) = 0. alu(:,:) = 0. b(:) = 0. 46 DO 888 L=1,NPTS c c Return spval if any inputs are spval (not available from time-gcm c usually because near bottom boundary at 30-35km): c if (ta(l).eq.spval.or.so2(l).eq.spval.or.so(l).eq.spval.or. + sn2(l).eq.spval.or.sh(l).eq.spval.or.so3(l).eq.spval.or. + sho2(l).eq.spval.or.soh0(l).eq.spval.or.sm(l).eq.spval) then dohv(:,l) = spval bohv(:,l) = spval if (iprint.gt.0) + write(6,"('ohrad skipping L=',i2,' because of spval')") l goto 151 endif INDEX=INDEX+1 ALT(L)=FLOAT(ILA+L-1)*1.E5 RK1(L)=1.4E-10*EXP(-470./TA(L)) RK2(L)=3.0E-11*EXP(200./TA(L)) RK3(L)=6.0E-34*(300./TA(L))**2.3 RK4(L)=5.7E-32*(300./TA(L))**1.6 RK5(L)=1.1E-14 RK6(L)=1.6E-12*EXP(-1140./TA(L)) RK7=4.8E-11 RK8=4.2E-12 C ... PROHVS: SECONDARY REACTION FOR PRODUCTION OF OH* <= 6 PROHVS(L)=RK2(L)*SHO2(L)*SO(L) PROHV(L)=RK1(L)*SH(L)*SO3(L) C..................................................................... C ... LOSS TERMS DO 50 KV=NVL,1,-1 Q1O2(KV)=AAO2(KV)*EXP(-EEO2(KV)/TA(L)) Q1N2(KV)=AAN2(KV)*EXP(-EEN2(KV)/TA(L)) CLO(KV)=AAO(KV)*EXP(-EEO(KV)/TA(L)) RQ1O2(KV)=RAAO2(KV)*EXP(-REE(KV)/TA(L)) RQ1N2(KV)=RAAN2(KV)*EXP(-REE(KV)/TA(L)) AA(KV,KV)=-(AOH(KV,1)+AOH(KV,2)+AOH(KV,3)+AOH(KV,4)+AOH(KV,5) 1 +AOH(KV,6)+Q1O2(KV)*SO2(L)+Q1N2(KV)*SN2(L)+CLO(KV)*SO(L) 2 +RQ1O2(KV)*SO2(L)+RQ1N2(KV)*SN2(L)+FAC0(KV)*(RK6(L)*SO3(L) 3 +RK7*SHO2(L)+2.*RK8*SOH0(L))) C ... A2: IS A WITHOUT QUENCHING A2(KV,KV)=-(AOH(KV,1)+AOH(KV,2)+AOH(KV,3)+AOH(KV,4)+AOH(KV,5) 1 +AOH(KV,6)+FAC0(KV)*(RK6(L)*SO3(L)+RK7*SHO2(L))) C B(KV)=-PROHV(L)*FVL(KV)-PROHVS(L)*FVLS(KV)-RK5(L)*SO3(L)*SHO2(L) 1 *FAC0(KV) 50 CONTINUE C DO 53 KV=NVL,2,-1 IF(INDEX.GE.2) GOTO 66 SQV(KV)=0.0 IF(KV-7) 57,59,59 57 KIM(KV)=KV-1 GOTO 63 59 KIM(KV)=6 GOTO 63 63 DO 65 JV=1,KIM(KV) FQV(KV,JV)=EXP(-FLOAT(JV-1)) SQV(KV)=FQV(KV,JV)+SQV(KV) 65 CONTINUE 66 DO 67 JV=1,KIM(KV) XMQO2(KV-JV,KV)=Q1O2(KV)*FQV(KV,JV)/SQV(KV) RMQO2(KV-JV,KV)=RQ1O2(KV)*FQV(KV,JV)/SQV(KV) 67 CONTINUE 53 CONTINUE C ... PRODUCTION TERMS DO 60 KV=1,NVL-1 AA(KV,KV+1)=AOH(KV+1,1)+XMQO2(KV,KV+1)*SO2(L)*XMMQ 1 +Q1O2(KV+1)*SO2(L)*RORIG+Q1N2(KV+1)*SN2(L) A2(KV,KV+1)=AOH(KV+1,1) 60 CONTINUE DO 70 KV=1,NVL-2 AA(KV,KV+2)=AOH(KV+2,2)+XMQO2(KV,KV+2)*SO2(L)*XMMQ A2(KV,KV+2)=AOH(KV+2,2) 70 CONTINUE DO 80 KV=1,NVL-3 AA(KV,KV+3)=AOH(KV+3,3)+XMQO2(KV,KV+3)*SO2(L)*XMMQ A2(KV,KV+3)=AOH(KV+3,3) 80 CONTINUE DO 90 KV=1,NVL-4 AA(KV,KV+4)=AOH(KV+4,4)+XMQO2(KV,KV+4)*SO2(L)*XMMQ A2(KV,KV+4)=AOH(KV+4,4) 90 CONTINUE DO 100 KV=1,NVL-5 AA(KV,KV+5)=AOH(KV+5,5)+XMQO2(KV,KV+5)*SO2(L)*XMMQ A2(KV,KV+5)=AOH(KV+5,5) 100 CONTINUE DO 110 KV=1,NVL-6 AA(KV,KV+6)=AOH(KV+6,6)+XMQO2(KV,KV+6)*SO2(L)*XMMQ A2(KV,KV+6)=AOH(KV+6,6) 110 CONTINUE DO 120 KV=1,NVL-7 AA(KV,KV+7)=0.0 A2(KV,KV+7)=0.0 120 CONTINUE DO 125 KV=1,NVL-8 AA(KV,KV+8)=0.0 A2(KV,KV+8)=0.0 125 CONTINUE AA(1,10)=0.0 A2(1,10)=0.0 DO 130 KV=2,NVL AA(KV,KV-1)=RMQO2(KV-1,KV)*SO2(L)+RQ1N2(KV-1)*SN2(L) 130 CONTINUE DO 140 KV=3,NVL DO 141 I=1,KV-2 AA(KV,I)=RMQO2(I,KV)*SO2(L)*XMMQ 141 CONTINUE 140 CONTINUE DO 155 KV=2,NVL DO 156 I=1,KV-1 156 A2(KV,I)=0.0 155 CONTINUE DO 145 KV=1,NVL FRAQ(KV,L)=(Q1O2(KV)*SO2(L)+Q1N2(KV)*SN2(L)+CLO(KV)*SO(L)+ 1 RQ1O2(KV)*SO2(L)+RQ1N2(KV)*SN2(L)+FAC0(KV)*(RK6(L)* 2 SO3(L)+RK7*SHO2(L)))/(-AA(KV,KV)) 145 CONTINUE C ... CALLING LINEAR EQUATIONS SOLVING SUBROUTINES CALL DECOMP(NVL,NVL,AA,ALU,DET,IPS,IER) CALL SOLVE(NVL,NVL,ALU,B,IPS,X) DO 150 KV=1,NVL DOHV(KV,L)=X(KV) C ... DOHV2 : IS THE # DENSITY WITHOUT QUENCHING 150 continue c150 DOHV2(KV,L)=X2(KV) c c btf 9/11/96: fix units to be either photons/cm3-sec or watts/cm3 c iwatts=0: mol/cm3 * ph/mol-s = ph/cm3-s c iwatts=1: mol/cm3 * ph/mol-s * J/ph = J/cm3-s = watts/cm3 c do kv=1,nvl do iv=1,6 if (kv-1-iv.ge.0) then ! valid band range bohv((iv-1)*nvl+kv,l) = dohv(kv,l)*aoh(kv,iv) ! photons/cm3-sec if (iwatts.gt.0) ! watts/cm3 + bohv((iv-1)*nvl+kv,l) = bohv((iv-1)*nvl+kv,l)*energy(kv) endif enddo enddo 151 continue C WRITE(6,55) ALT(L)/1.E+5,(DOHV(KV,L)/DOHV(4,L),KV=1,NVL) C WRITE(7,55) (FRAQ(KV,L),KV=1,NVL) C WRITE(8,55) (DOHV2(KV,L),KV=1,NVL) if (iprint.gt.0) WRITE(6,55) ALT(L)/1.E+5,(DOHV(KV,L),KV=1,NVL) 55 FORMAT(F6.1,10E11.3) C...................................................................... C DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD 888 CONTINUE c do i=1,nvl c write(6,"('ohrad: vib level i=',i2,' dohv(i,:)=',/(6e12.4))") c + i,(dohv(i,k),k=1,npts) c enddo C return END C LINEQSV FROM NSSL 07/29/80 SUBROUTINE DECOMP (N,NDIM,A,LU,D,IPS,IER) C C PACKAGE LINEQSV C C C LATEST REVISION JANUARY 1979 C C**** ***** Sept. 1, 1983 revised setting D as double precision number c**** it is never really used ,but otherwise it can overflow c** C PURPOSE LINEQSV IS A PACKAGE OF THREE C SUBROUTINES--DECOMP, SOLVE AND IMPRUV--WHICH C MAY BE USED TO SOLVE A SYSTEM OF LINEAR C EQUATIONS BY GAUSSIAN ELIMINATION WITH PARTIAL C PIVOTING AND (OPTIONALLY) WITH ITERATIVE C IMPROVEMENT. C C ACCESS CARDS *FORTRAN,S=ULIB,N=LINEQSV C C C USAGE TO SOLVE A SYSTEM OF EQUATIONS, AX = B, THE C USER FIRST CALLS DECOMP BY C CALL DECOMP (N,NDIM,A,LU,D,IPS,IER) C WHICH DECOMPOSES A INTO THE PRODUCT OF A UNIT C LOWER TRIANGULAR MATRIX AND AN UPPER TRIANGULAR C MATRIX. AFTER CHECKING D TO INSURE THAT THE C DETERMINANT OF A IS NOT ZERO, THE USER THEN C FOLLOWS THIS CALL WITH A CALL TO SOLVE BY C CALL SOLVE (N,NDIM,LU,B,IPS,X) C WHICH USES THE DECOMPOSITION TO FIND THE C SOLUTION X. IF DESIRED, THE USER MAY THEN TRY C TO IMPROVE THE ACCURACY OF THE SOLUTION BY C CALLING IMPRUV BY C CALL IMPRUV (N,NDIM,A,LU,B,IPS,X,R,DIGITS, C IER) C IF THE USER WANTS TO SOLVE SEVERAL LINEAR C SYSTEMS OF EQUATIONS WITH THE SAME COEFFICIENT C MATRIX A, ONLY ONE CALL TO DECOMP IS NECESSARY. C THIS IS THEN FOLLOWED BY A CALL TO SOLVE FOR C EACH DIFFERENT B VECTOR. C . EACH OF THE ABOVE ROUTINES IS WRITTEN UP IN C DETAIL BELOW. C C ENTRY POINTS DECOMP, SOLVE, IMPRUV C C SPACE REQUIRED 723 (OCTAL) = 467 (DECIMAL) LOCATIONS C C SUBROUTINE DECOMP (N,NDIM,A,LU,D,IPS,IER) C C C DIMENSION OF A(NDIM,N),LU(NDIM,N),IPS(N) C ARGUMENTS C C PURPOSE DECOMPOSE MATRIX INTO THE PRODUCT OF UPPER AND C UNIT LOWER TRIANGULAR MATRICES USING GAUSSIAN C ELIMINATION WITH PARTIAL PIVOTING. C C USAGE CALL DECOMP (N,NDIM,A,LU,D,IPS,IER) C C ARGUMENTS C C ON INPUT N C ORDER OF MATRIX A. C C NDIM C DECLARED FIRST (ROW) DIMENSIONS OF A AND LU C IN THE CALLING PROGRAM. C C A C MATRIX TO BE DECOMPOSED. C C ON OUTPUT LU C LOWER AND UPPER TRIANGULAR MATRICES OF A, C WHERE C LU(I,J), FOR I .LE. J, C IS THE UPPER TRIANGULAR MATRIX WITH ROW C INTERCHANGES COMPLETED, AND C LU(I,J), FOR I .GT. J, C IS THE NEGATIVE OF LOWER TRIANGULAR MATRIX C EXCLUDING DIAGONAL ELEMENTS, BUT WITH ROW C INTERCHANGES PARTIALLY COMPLETED. C C IF MATRIX A NEED NOT BE SAVED, THEN THE C MATRIX A MAY BE USED AS THE DECOMPOSITION C MATRIX LU IN THE CALL, I.E., C CALL DECOMP (N,NDIM,A,A,D,IPS,IER) C NOTE: IF SUBROUTINE IMPRUV WILL BE USED, A C MUST BE RETAINED. C C D C DETERMINANT OF A. C C IPS C VECTOR (OF DIMENSION AT LEAST N) STORING ROW C PERMUTATIONS RESULTING FROM PIVOTING, WHERE C IPS(K) = KTH PIVOT ROW AFTER ROW INTERCHANGES C FROM PREVIOUS PIVOTS. C C IER C ERROR FLAG USED BY STANDARD ERROR MESSAGE C ROUTINE, ULIBER, AND RETURNED TO USER. C = 0, NO ERROR. C = 32, DETERMINANT OF A EQUALS ZERO, HENCE C DECOMPOSITION IS IMPOSSIBLE. A C MESSAGE IS PRINTED. C C C NOTE ULIBER IS CALLED TO OUTPUT ERROR MESSAGES. C C C C C COMMON BLOCKS NONE C C I/O AN NCAR RESIDENT ROUTINE, ULIBER, IS USED TO C OUTPUT ERROR MESSAGES. C C PRECISION SINGLE C C REQUIRED ULIB NONE C ROUTINES C C SPECIALIST R. K. SATO, NCAR, BOULDER, COLORADO 80303 C C LANGUAGE FORTRAN C C HISTORY TRANSCRIBED FROM A PAPER BY C. B. MOLER, C LINEAR EQUATION SOLVER, COMM. ACM 15 C (APRIL, 1972), 74. C C double precision D REAL LU DIMENSION A(NDIM,NDIM) ,LU(NDIM,NDIM) , 1 IPS(NDIM) C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR cc** CALL Q8QST4( 4HNSSL , 7HLINEQSV ,6HDECOMP ,10HVERSION 10) IER = 0 C C INITIALIZE ARRAYS C C DO 101 J=1,N IPS(J) = J 101 CONTINUE C DO 103 J=1,N DO 102 I=1,N LU(I,J) = A(I,J) 102 CONTINUE 103 CONTINUE 104 D = 1.0D0 C C LU DECOMPOSITION USING GAUSSIAN ELIMINATION WITH PARTIAL PIVIOTING. C DO 110 K=1,N IF (K .EQ. N) GO TO 109 KP1 = K+1 M = K C C FIND PIVOT ELEMENT. C DO 105 I=KP1,N IF (ABS(LU(I,K)) .GT. ABS(LU(M,K))) M = I 105 CONTINUE C C SAVE ROW PERMUTATIONS. C IPS(K) = M IF (M .NE. K) D = -D T = LU(M,K) LU(M,K) = LU(K,K) LU(K,K) = T IF (T .EQ. 0.0) GO TO 112 C C CALCULATE AND STORE ELEMENTS OF FACTORS, L AND U. C DO 106 I=KP1,N LU(I,K) = -LU(I,K)/T 106 CONTINUE DO 108 J=KP1,N T = LU(M,J) LU(M,J) = LU(K,J) LU(K,J) = T IF (T .EQ. 0.0) GO TO 108 DO 107 I=KP1,N LU(I,J) = LU(I,J)+LU(I,K)*T 107 CONTINUE 108 CONTINUE 109 IF (LU(K,K) .EQ. 0.0) GO TO 112 110 CONTINUE C C COMPUTE DETERMINANT OF MATRIX A. C DO 111 K=1,N D = D*LU(K,K) 111 CONTINUE RETURN 112 IER = 32 D = 0.0 CALL ULIBER_old (IER,26H SINGULAR MATRIX IN DECOMP,26) RETURN END SUBROUTINE SOLVE (N,NDIM,LU,B,IPS,X) C C C DIMENSION OF LU(NDIM,N),B(N),IPS(N),X(N) C ARGUMENTS C C PURPOSE SOLVES FOR VECTOR X OF UNKNOWNS OF THE LINEAR C SYSTEM AX = B USING THE LU DECOMPOSITION OF A C FROM SUBROUTINE DECOMP. C C USAGE CALL SOLVE (N,NDIM,LU,B,IPS,X) C C ARGUMENTS C C ON INPUT N C ORDER OF MATRIX A (NUMBER OF EQUATIONS). C C NDIM C DECLARED FIRST (ROW) DIMENSION OF LU IN THE C CALLING PROGRAM. C C LU C LU DECOMPOSITION MATRIX FROM SUBROUTINE C DECOMP. C C B C RIGHT HAND SIDE VECTOR (OF DIMENSION AT C LEAST N) OF CONSTANTS. C C IPS C VECTOR (OF DIMENSION AT LEAST N) FROM C SUBROUTINE DECOMP STORING ROW PERMUTATIONS C RESULTING FROM PIVOTING, WHERE IPS(K) = KTH C PIVOT ROW AFTER ROW INTERCHANGES FROM C PREVIOUS PIVOTS. C C ON OUTPUT X C SOLUTION VECTOR (OF DIMENSION AT LEAST N). C IF VECTOR B NEED NOT BE SAVED, THEN THE C VECTOR B MAY BE USED AS THE SOLUTION VECTOR C IN THE CALL, I.E., C CALL SOLVE (N,NDIM,LU,B,IPS,B) C NOTE: IF SUBROUTINE IMPRUV WILL BE USED, B C MUST BE RETAINED. C C NOTE AN NCAR SYSTEM ROUTINE, LOC, IS USED TO CHECK C EQUIVALENCE OF ARGUMENTS B AND X BY CHECKING C THEIR ADDRESSES TO ELIMINATE UNNECESSARY C TRANSFER OF B TO X WHEN POSSIBLE. C C COMMON BLOCKS NONE C C I/O NONE C C PRECISION SINGLE C C REQUIRED ULIB NONE C ROUTINES C C SPECIALIST R. K. SATO, NCAR, BOULDER, COLORADO 80303 C C LANGUAGE FORTRAN C C HISTORY TRANSCRIBED FROM A PAPER BY C. B. MOLER, C LINEAR EQUATION SOLVER, COMM. ACM 15 C (APRIL, 1972), 74. C C REAL LU DIMENSION LU(NDIM,NDIM) ,B(NDIM) ,IPS(NDIM) , 1 X(NDIM) C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR cc*** CALL Q8QST4( 4HNSSL , 7HLINEQSV ,5HSOLVE ,10HVERSION 10) C C C C C DO 101 K=1,N X(K) = B(K) 101 CONTINUE 102 IF (N .EQ. 1) GO TO 107 NM1 = N-1 DO 104 K=1,NM1 KP1 = K+1 M = IPS(K) T = X(M) X(M) = X(K) X(K) = T DO 103 I=KP1,N X(I) = X(I)+LU(I,K)*T 103 CONTINUE 104 CONTINUE DO 106 KB=1,NM1 KM1 = N-KB K = KM1+1 X(K) = X(K)/LU(K,K) T = -X(K) DO 105 I=1,KM1 X(I) = X(I)+LU(I,K)*T 105 CONTINUE 106 CONTINUE 107 X(1) = X(1)/LU(1,1) RETURN END SUBROUTINE IMPRUV (N,NDIM,A,LU,B,IPS,X,R,DIGITS,IER) C C C DIMENSION OF A(NDIM,N),LU(NDIM,N),B(N),IPS(N),X(N),R(N) C ARGUMENTS C C PURPOSE ITERATIVE IMPROVEMENT OF THE SOLUTION OBTAINED C FROM SUBROUTINE SOLVE. C C USAGE CALL IMPRUV (N,NDIM,A,LU,B,IPS,X,R,DIGITS,IER) C C ARGUMENTS C C ON INPUT N C ORDER OF MATRIX A (NUMBER OF EQUATIONS). C C NDIM C DECLARED FIRST (ROW) DIMENSIONS OF A AND LU C IN THE CALLING PROGRAM. C C A C ORIGINAL COEFFICIENT MATRIX. C C LU C LU DECOMPOSITION MATRIX OF A FROM SUBROUTINE C DECOMP. C C B C ORIGINAL RIGHT HAND SIDE VECTOR (OF DIMENSION C AT LEAST N) OF CONSTANTS. C C IPS C VECTOR (OF DIMENSION AT LEAST N) FROM C SUBROUTINE DECOMP STORING ROW PERMUTATIONS C RESULTING FROM PIVOTING, WHERE IPS(K) = KTH C PIVOT ROW AFTER ROW INTERCHANGES FROM C PREVIOUS PIVOTS. C C X C SOLUTION VECTOR (OF DIMENSION AT LEAST N) C FROM SUBROUTINE SOLVE. C C R C WORK ARRAY (OF DIMENSION AT LEAST N) USED IN C THE COMPUTATION OF RESIDUALS. C C ON OUTPUT X C SOLUTION VECTOR AFTER ITERATIVE IMPROVEMENT. C C DIGITS C APPROXIMATE NUMBER OF CORRECT DIGITS IN X. C C IER C ERROR FLAG USED BY STANDARD ERROR MESSAGE C ROUTINE ULIBER AND RETURNED TO USER. C = 0, NO ERROR. C = 1, NO CONVERGENCE IN THE ITERATION. C C COMMON BLOCKS NONE C C I/O AN NCAR RESIDENT ROUTINE, ULIBER, IS CALLED TO C OUTPUT THE ERROR MESSAGE. C C PRECISION DOUBLE PRECISION IS USED IN THE CALCULATION OF C THE RESIDUALS. C C REQUIRED ULIB NONE C ROUTINES C C SPECIALIST R. K. SATO, NCAR, BOULDER, COLORADO 80303 C C LANGUAGE FORTRAN C C HISTORY TRANSCRIBED FROM FORSYTHE AND MOLER, COMPUTER C SOLUTIONS OF LINEAR ALGEBRAIC SYSTEMS, 1967. C C C C DIMENSION A(NDIM,NDIM) ,LU(NDIM,NDIM) , 1 B(NDIM) ,X(NDIM) ,IPS(NDIM) ,R(NDIM) REAL LU DOUBLE PRECISION SUM C C EPS AND ITMAX ARE MACHINE DEPENDENT CONSTANTS. C EPS = LARGEST FLOATING POINT NUMBER SUCH THAT 1.0 + EPS = 1.0 C ITMAX = NUMBER OF BINARY BITS IN THE MANTISSA OF A FLOATING POINT C NUMBER AND IS THE MAXIMUM NUMBER OF ITERATIONS. DATA EPS,ITMAX/1.E-15,48/ C THE FOLLOWING CALL IS FOR GATHERING STATISTICS ON LIBRARY USE AT NCAR c*** CALL Q(QSTR( 4HNSSL , 7HLINEQSV ,6HIMPRUV ,10HVERSION 10) C IER = 0 XNORM = 0.0 DO 101 I=1,N XNORM = AMAX1(XNORM,ABS(X(I))) 101 CONTINUE IF (XNORM .GT. 0.0) GO TO 102 DIGITS = -ALOG10(EPS) GO TO 108 102 DO 107 ITER=1,ITMAX C C CALCULATION OF DOUBLE PRECISION RESIDUALS. C DO 104 I=1,N SUM = 0.0 DO 103 J=1,N SUM = SUM+A(I,J)*X(J) 103 CONTINUE SUM = B(I)-SUM R(I) = SUM 104 CONTINUE CALL SOLVE (N,NDIM,LU,R,IPS,R) DXNORM = 0.0 DO 105 I=1,N T = X(I) X(I) = X(I)+R(I) DXNORM = AMAX1(DXNORM,ABS(X(I)-T)) 105 CONTINUE IF (ITER .GT. 1) GO TO 106 DIGITS = -ALOG10(AMAX1(DXNORM/XNORM,EPS)) 106 IF (DXNORM .LE. EPS*XNORM) GO TO 108 107 CONTINUE C C ITERATION DID NOT CONVERGE C IER = 1 CALL ULIBER_old (IER, + 54H NO CONVERGENCE IN IMPRUV. MATRIX IS NEARLY SINGULAR. , + 54) 108 CONTINUE RETURN C C REVISION HISTORY--- C C JUNE 1977 REMOVED REFERENCE TO THE LOC FUNCTION C TO ENHANCE PORTABILITY. C C JANUARY 1978 DELETED REFERENCES TO THE *COSY CARDS, MOVED C THE REVISION HISTORIES TO APPEAR BEFORE THE C FINAL END CARD, AND MOVED THE INITIAL COMMENT C CARDS TO APPEAR AFTER THE FIRST SUBROUTINE CARD C JANUARY 1979 REMOVED COMMENTS ABOUT LOC FUNCTION. C----------------------------------------------------------------------- END C C C C C C C C TRUNC FROM PORTLIB 07/29/80 FUNCTION TRUNC(X) C C TRUNC IS A PORTABLE FORTRAN FUNCTION WHICH TRUNCATES A VALUE TO THE C MACHINE SINGLE PRECISION WORD SIZE, REGARDLESS OF WHETHER LONGER C PRECISION INTERNAL REGISTERS ARE USED FOR FLOATING POINT ARITHMETIC IN C COMPUTING THE VALUE INPUT TO TRUNC. THE METHOD USED IS TO FORCE A C STORE INTO MEMORY BY USING A COMMON BLOCK IN ANOTHER SUBROUTINE. C COMMON /VALUE/ V CALL STORES(X) TRUNC=V RETURN END SUBROUTINE STORES(X) C C FORCES THE ARGUMENT VALUE X TO BE STORED IN MEMORY LOCATION V. C COMMON /VALUE/ V V=X RETURN END C C C ULIBER FROM PORTLIB 07/29/80 SUBROUTINE ULIBER_old (IERR,MESSG,NMESSG) C C DIMENSION OF MESSG(1+(NMESSG-1)/NCPWD), WHERE NCPWD IS THE C ARGUMENTS NUMBER OF CHARACTERS THAT CAN BE STORED IN ONE C INTEGER WORD. C C LATEST REVISION APRIL 1977. C C PURPOSE PRINTS AN ERROR MESSAGE. C C USAGE CALL ULIBER (IERR,MESSG,NMESSG) C C ARGUMENTS C C ON INPUT IERR C ERROR NUMBER. IF IERR .LT. 32, THE ERROR IS C CONSIDERED TO BE NON-FATAL AND ULIBER RETURNS C AFTER PRINTING THE ERROR MESSAGE. IF IERR C .GE. 32 THE ERROR IS CONSIDERED FATAL, AND C ULIBER STOPS AFTER PRINTING THE MESSAGE. C MESSG C THE ERROR MESSAGE TO BE PRINTED. C NMESSG C THE LENGTH OF THE MESSAGE, IN CHARACTERS. C DIMENSION MESSG(1) INTEGER ERUNIT,PRUNIT C C ERUNIT SHOULD BE SET TO THE LOCAL UNIT NUMBER FOR ERROR MESSAGES. C DATA ERUNIT/0/,PRUNIT/6/ C IF (ERUNIT .EQ. 0) WRITE(PRUNIT,1000) C 1000 FORMAT(50H1ULIBER, A SUBROUTINE TO PRINT ERROR MESSAGES, HAS/ 1 50H BEEN CALLED, BUT NO LOCAL IMPLEMENTATION OF / 2 50H ULIBER HAS BEEN PROVIDED. TO PROPERLY IMPLEMENT / 3 50H THIS SUBROUTINE FOR THE LOCAL MACHINE AND / 4 50H ENVIRONMENT, PLEASE SEE THE COMMENTS IN ULIBER. /) C C REPLACE THE IF STATEMENT AND FORMAT STATEMENT ABOVE WITH THE FOLLOWING C CODE WHERE $NCPWD SHOULD BE REPLACED WITH THE NUMBER OF CHARACTERS C THAT CAN BE STORED IN ONE INTEGER WORD, AND $N SHOULD BE REPLACED BY C THE VALUE OF 1+(79/$NCPWD). C C MM=1+(NMESSG-1)/$NCPWD C WRITE(ERUNIT,1000) IERR,(MESSG(I),I=1,MM) C1000 FORMAT(18H ****ERROR NUMBER ,I5,25H,ERROR MESSAGE FOLLOWS.../ C 1 ($N A $NCPWD )) C 10 IF (IERR .GE. 32) STOP RETURN END C c oh_qrates.lowe c 23.72 148.9 78.64 19.94 4.053 0.620 26185.79 c 8.600 154.3 59.98 12.68 2.007 0.230 23949.83 c 2.340 145.1 42.91 7.165 0.847 0.063 21536.25 c 4.000 125.6 27.94 3.479 0.274 0.010 18951.90 c 11.05 99.42 15.88 1.315 0.051 0.000 16201.33 c 20.30 69.77 7.191 0.299 0.000 0.000 13287.19 c 28.12 40.33 2.032 0.000 0.000 0.000 10210.57 c 30.43 15.42 0.000 0.000 0.000 0.000 06971.37 c 22.74 0.000 0.000 0.000 0.000 0.000 03568.50 c 0.000 0.000 0.000 0.000 0.000 0.000 00000.00 c 1.700E-11 4.770E-13 2.500E-10 0.0000000 0.0000000 c 9.810E-12 2.746E-13 2.500E-10 1.700E-11 4.770E-13 c 5.420E-12 1.581E-13 2.500E-10 9.810E-12 2.746E-13 c 3.000E-12 9.103E-14 2.500E-10 5.420E-12 1.581E-13 c 1.700E-12 5.241E-14 2.500E-10 3.000E-12 9.103E-14 c 8.800E-13 3.017E-14 2.500E-10 1.700E-12 5.241E-14 c 5.200E-13 1.737E-14 2.500E-10 8.800E-13 3.017E-14 c 2.700E-13 1.000E-14 2.500E-10 5.200E-13 1.737E-14 c 1.300E-13 5.757E-15 10.50E-11 2.700E-13 1.000E-14 c 0.0000000 0.0000000 3.900E-11 1.300E-13 5.757E-15 c