C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C ASROT - ASYMMETRIC ROTOR ENERGY LEVELS, TRANSITIONS, LINE STRENGTHS AND C INTENSITIES C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C - J UP TO 300 AND UP TO 92104 LINES IN PREDICTION C - WATSON'S REDUCED HAMILTONIAN (repr Ir, reductions 'A' and 'S') C UP TO AND INCLUDING DECADIC TERMS C - ALL TYPES OF TRANSITIONS AND BROAD RANGE OF SELECTION RULES ARE C SUPPORTED, CLEAR INDICATION OF TRANSITION TYPE IS GIVEN IN THE C PRINTOUT C - LINE INTENSITIES (MAXIMUM ABSORPTION COEFFICIENTS) ARE CALCULATED C FOR A SPECIFIED TEMPERATURE, VALUES ASSUMED FOR OTHER KEY C PARAMETERS ARE PRINTED AT THE END OF THE FREQUENCY LISTING C C C This program, though much developed, is still recognisably descended C from ASROT of P.J.Mjoberg et al (ca 1972) - handling of Wang C factorisation is common with that used in ASFIT. C C C ver. 2.06.2015 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c 14.05.95: Warning on reaching the maximum number of lines c 13.04.97: Overhaul of line packing/unpacking c 14.04.97: New line sorting code c 15.04.97: New exit for too many lines c 29.05.97: More screen information+some mods. c 4.09.03: Increase to MAXDIM=91204 C 21.03.14: Code tweaks for compilation with contemporary FORTRAN compilers C 2.06.15: Small mods for Linux compilation C C This program is standard F77 and this version requires about 5 Mb of C memory. C C_____________________________________________________________________________ C C ASROT can be used completely interactively, questions being asked by the C program being (hopefully) self explanatory. C C In cases where many constants are to be put in, online typing-in C may not be preferred and it is possible to read these from an ASFIT C produced or manually edited data file of constants with the following C format: C C Line.1: NCON = 8 C Line.2: A = 48553.1122 C . B = 5010.04024 C . C = 4529.89707 C DJ = .00230008 C DJK = -.0691118 C DK = 1.74704 C dJ = .000477929 C dK = .029516 C C This file is read in under a fixed (A7,F30) FORMAT, the alphanumeric C field being a dummy field for constant identification. NCON can be set C as required, but no constants can be skipped i.e. they have to adhere to C the following sequence: C C 'A' reduction: C C A, B, C C DJ, DJK, DK, dJ, dK C HJ, HJK, HKJ, HK, hJ, hJK, hK C LJ, LJJK, LJK, LKKJ, LK, lJ, lJK, lKJ, lK C PJ, PJJK, PJK, PKJ, PKKJ, PK, pJ, pJJK, pJK, pKKJ, pK C C C 'S' reduction: C C A, B, C C DJ, DJK, DK, d1, d2 C HJ, HJK, HKJ, HK, h1, h2, h3 C LJ, LJJK, LJK, LKKJ, LK, l1, l2, l3, l4 C PJ, PJJK, PJK, PKJ, PKKJ, PK, p1, p2, p3, p4, p5 C C C - On interactive input the cutoff in constants is specified by means C of the order of the Hamiltonian eg. 4 for up to and including quartic, C 6 - sextic etc. Any other values are truncated to their nearest lower C multiple of 2. C C - Various printout combinations are possible (see appropriate question C in the interactive mode), options 1 (frequencies only)+sorting+discard, C and option 3 (energy levels only) being recommended for routine use C C - ASROT has built in detection of close K-1, or K+1 doublets on sorting; C when these are closer than 0.001 in frequency (FRESOL in SORT0) an C entry for a doublet is printed - this has a D after the frequency and C only the common K quantum number is printed. C C - ASROT writes identified transitions to a temporary file, each line C taking 24 bytes. For J=300, i.e. 22801 lines this implies C a maximum of ca 550Kbytes for this file. The ASCII output C file may be considerably larger, up to over 1.8 Mb C C - if ASROT is run under options: 'frequencies only, lines to be sorted, C unsorted lines to be discarded' then the results can be viewed C graphically with ASCP_L (or ASCP). C C - On repetitive use of ASROT for similar predictions interactive use C becomes tedious and many systems allow a pipeline mechanism such as C ASROT reduction 'A', IRED=1 -> reduction 'S' C IRED=0 IF(IP.LT.0)THEN IRED=1 IP=IABS(int(IP)) DO 3400 I=1,NCONST ISTP(I)=ISTPS(I) 3400 CONTINUE ENDIF C 2901 WRITE(*,4301) 4301 FORMAT(/' JMIN = ',$) READ(*,*,ERR=2901)JMIN 2903 WRITE(*,4302) 4302 FORMAT(' JRMAX = ',$) READ(*,*,ERR=2903)JRMAX 2904 WRITE(*,4303) 4303 FORMAT(' JQMAX = ',$) READ(*,*,ERR=2904)JQMAX 200 WRITE(*,201) 201 FORMAT(' K-1.MAX = ',$) READ(*,*,ERR=200)KMMAX C 2905 WRITE(*,4304) 4304 FORMAT(/' LINES TO BE SORTED (1=YES, 0=NO) ? ',$) READ(*,40,ERR=2905)IS IF(IS.LT.0.OR.IS.GT.1)GOTO 2905 ISON=0 IF(IP.NE.1.OR.IS.NE.1)GOTO 4390 2906 WRITE(*,4391) 4391 FORMAT(/' ARE UNSORTED FREQUENCIES TO BE DISCARDED ? ',$) READ(*,40,ERR=2906)ISON IF(ISON.LT.0.OR.ISON.GT.1)GOTO 2906 4390 WRITE(*,4305) 4305 FORMAT(/' COMMENT :') READ(*,4000,ERR=4390)COMENT 4000 FORMAT(72A1) WRITE(4,42)COMENT WRITE(4,43)namred(ired+1) DO 5848 I=1,NCONST A(I)=0.D0 5848 CONTINUE C 2907 WRITE(*,6853) 6853 FORMAT(/' ARE ROTATIONAL CONSTANTS AVAILABLE ON DISK ? ',$) READ(*,40,ERR=2907)I IF(I.NE.1)GOTO 6854 2908 WRITE(*,6855) 6855 FORMAT(/' NAME OF FILE CONTAINING ROTATIONAL CONSTANTS : ',$) 2909 READ(*,'(A)',ERR=2908)FILROT OPEN(3,FILE=FILROT,STATUS='OLD',ERR=2908) READ(3,'(A)')LINE READ(LINE,6866,ERR=2950)BLUFF,NCON 6866 FORMAT(A7,I5) GOTO 2951 C C...Constants from a .ASL file C 2950 READ(3,'(A)')LINE NCON=0 DO 2953 I=1,NCONST READ(3,'(I5,F30.19)',ERR=6851)K,A(I) NCON=NCON+1 2953 CONTINUE GOTO 6851 C C...Constants from a .CON file C 2951 DO 6856 I=1,NCON 2911 READ(3,6867,ERR=2910)BLUFF,A(I) 6867 FORMAT(A7,F30.19) GOTO 6856 2910 WRITE(*,'(1H+,A)')CHAR(7) GOTO 2911 6856 CONTINUE GOTO 6851 C 6854 WRITE(*,6849) 6849 FORMAT(/' Order of Hamiltonian (eg. 2, 4, 6 etc.) = ',$) READ(*,40,ERR=6854)IORDER IORDER=(IORDER/2)*2 NCON=0 DO 6850 I=IORDER,2,-2 NCON=NCON+I+1 6850 CONTINUE WRITE(*,4306) 4306 FORMAT(/' ROTATIONAL AND CENTRIFUGAL DISTORTION CONSTANTS /MHz', 1 ' :'/) DO 5000 I=1,NCONST WRITE(*,5001)ISTP(I) 2916 READ(*,41,ERR=2915)A(I) IF(I.GE.NCON)GOTO 6851 GOTO 5000 2915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 2916 5000 CONTINUE 5001 FORMAT(20X,A6,' ',$) C 6851 WRITE(*,7005) 7005 FORMAT(//' CALCULATION WILL BE CARRIED OUT WITH FOLLOWING ', 1 'CONSTANTS:'/) DO 6870 I=1,NCON IF(A(I).EQ.0.D0)GOTO 6870 WRITE(*,44)ISTP(I),A(I) WRITE(4,44)ISTP(I),A(I) 6870 CONTINUE C IF(IP.EQ.3)GOTO 5 WRITE(*,4307) 4307 FORMAT(/' DIPOLE MOMENT COMPONENTS /D :'/) DO 5002 I=1,3 2918 WRITE(*,5001)ISTP(I+NCONST) READ(*,41,ERR=2917)D(I) GOTO 5002 2917 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 2918 5002 CONTINUE 2919 WRITE(*,5003) 5003 FORMAT(/' CUTOFF LINE STRENGTH = ',$) READ(*,41,ERR=2919)EPS 3900 WRITE(*,3901) 3901 FORMAT(' ROTATIONAL TEMPERATURE /K = ',$) READ(*,41,ERR=3900)TEMP 2920 WRITE(*,4308) 4308 FORMAT(' LOW FREQUENCY LIMIT /MHz = ',$) READ(*,41,ERR=2920)FMIN 2921 WRITE(*,4309) 4309 FORMAT(' HIGH FREQUENCY LIMIT /MHz = ',$) READ(*,41,ERR=2921)FMAX WRITE(4,50) WRITE(4,44)(ISTP(I+NCONST),D(I),I=1,3) C 5 JMIN=ABS(JMIN) JRMAX=MIN(ABS(JRMAX),maxj) JQMAX=MIN(ABS(JQMAX),maxj) JMAX=MAX(JRMAX,JQMAX) IF(KMMAX.EQ.0)KMMAX=JMAX K=1 IF(JMIN.GT.JMAX)GOTO 28 K=2 IF(A(1).LT.A(2).OR.A(2).LT.A(3))GOTO 28 IF(IP.NE.3)GOTO 9 WRITE(4,43)namred(ired+1) GOTO 14 9 EPS=DABS(EPS) FMIN=DABS(FMIN) FMAX=DABS(FMAX) K=4 IF(FMAX.LT.FMIN)GOTO 29 RJR=0.D0 RJQ=0.D0 DO 13 I=1,3 RJR=RJR*16.D0 RJQ=RJQ*4.D0 IF(D(4-I).EQ.0.D0)GOTO 13 RJR=RJR+15.D0 RJQ=RJQ+12288.D0 13 CONTINUE 14 M=JMAX/2+1 DO 15 I=1,M DO 15 K=1,M H(I,K)=0.D0 15 CONTINUE WRITE(*,2300) 2300 FORMAT(/' ***** CALCULATION IN PROGRESS, NOW WORKING ON J = '//) JQMIN=JMIN IF(JQMAX.GT.JMIN)GOTO 7006 JQMAX=0 JQMIN=0 7006 WRITE(4,46)JMIN,JRMAX,FMIN,FMAX,JQMIN,JQMAX,EPS,KMMAX IF(ISON.EQ.0)WRITE(4,4392) JJMIN=JMIN+1 JJMAX=JMAX+1 CLOSE(4) OPEN(4,FILE=FILNAM,ACCESS='APPEND') C C DO 24 JJ=JJMIN,JJMAX J=JJ-1 WRITE(*,2400)J 2400 FORMAT(1H+,I4,$) Intel/Windows c2400 FORMAT( I4,$) gfortran/Linux L=15 IF(IRED.NE.1)THEN CALL ASROHA(NCON) ELSE CALL ASROHS ENDIF IF(IP.LE.1)GOTO 19 IF(ISON.EQ.0)WRITE(4,50) DO 118 N=1,4 L=N-1 M=N/2*2 K=L/2+(4-N)/2*M L=J-K+M-L M=(J-K+2)/2 NN=(J-J/2*2)*4+N IF(M.LT.1)GOTO 118 DO 18 I=1,M WRITE(4,47)J,K,L,E(NN,I) K=K+2 L=L-2 18 CONTINUE 118 CONTINUE IF(IP.EQ.3)GOTO 24 19 RJS=1.D-6 IF(J.EQ.JMIN)GOTO 21 IF(J.GT.JRMAX)GOTO 21 RJS=RJS+RJR 21 IF(J.GT.JQMAX)GOTO 23 RJS=RJS+RJQ 23 CALL DIPTRA(DIDUTY,RJS,ISON,KMMAX,iexit) if(iexit.eq.1)goto 244 24 CONTINUE C C C...IP = 1/2/3 calculate freq/freq+el/el 244 IF(IP.EQ.3)GOTO 2 C C...IS = 0/1 do not sort/sort frequencies 25 IF(IS.NE.0)CALL SORT0(FV,SIGMA,GI,HALFW,ABUND,TEMP) GOTO 31 C C...Various error reports C 28 WRITE(*,48)K,JMIN,JMAX IF(IP.EQ.3)GOTO 2 GOTO 30 29 WRITE(*,49)K,EPS,FMIN,FMAX 30 IF(IS.NE.1.OR.NT.GT.1)GOTO 25 C 31 WRITE(*,50) C C 2 STOP END C C_____________________________________________________________________________ C SUBROUTINE ASROHA(NCON) C C ASYMMETRIC ROTOR ENERGY LEVELS FOR A GIVEN J - Reduction 'A' C C The coding is set up to confer a slight speed advantage when constants C do not go beyond quartic. C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=300, MAXDIM=MAXJ/2+1) C COMMON /ASRO/FMIN,FMAX,A(NCONST),EA,EB,EE,EPS, 1 S(3),E(8,MAXDIM),H(MAXDIM,MAXDIM),V(MAXDIM),R, 1 ID(MAXDIM),I,J,L,NT,IP COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),R(MAXDIM,MAXDIM) INTEGER*4 MM,NT C C...various ancillary coefficients Q=J Q=Q*(J+1) C=(A(2)+A(3))*0.5D0 B=(-A(4)*Q+C)*Q C=-A(5)*Q+A(1)-C D=(A(2)-A(3))*0.25D0-A(7)*Q F=-A(6) G=-A(8) IF(NCON.LE.8)GOTO 50 QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q B=B+QQQ*A(9)+QQQQ*A(16)+QQQQQ*A(25) C=C+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ F=F+A(11)*Q+A(18)*QQ+A(27)*QQQ FK6=A(12)+A(19)*Q+A(28)*QQ FK8=A(20)+A(29)*Q FK10=A(30) D=D+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ G=G+A(14)*Q+A(22)*QQ+A(32)*QQQ FKK4=0.5D0*(A(15)+A(23)*Q+A(33)*QQ) FKK6=0.5D0*(A(24)+A(34)*Q) FKK8=0.5D0*A(35) C C...Process the four Wang matrices for a given J C 50 DO 12 N=1,4 LN=L-2*(L/2) L=L/2 IF(LN.EQ.0)GOTO 12 K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K IF(KE.LT.KA)GOTO 12 I=0 KA=KA+1 KE=KE+1 C C...go over K (the loop variable KK is K+1 to avoid problems for K=0) DO 7 KK=KA,KE,2 K=KK-1 I=I+1 MM=K MM=MM*K FLTA=MM IF(NCON.GT.8)GOTO 51 H(I,I)=(FLTA*F+C)*FLTA+B GOTO 52 51 H(I,I)=((((FK10*FLTA+FK8)*FLTA+FK6)*FLTA+F)*FLTA+C)*FLTA+B RMP2=(K+2)**2 SRMP2=RMP2*RMP2 SKSQ=FLTA*FLTA FLT4=SKSQ+SRMP2 FLT6=SKSQ*FLTA+SRMP2*RMP2 FLT8=SKSQ*SKSQ+SRMP2*SRMP2 52 IF(KK.GE.KE)GOTO 7 MM=MM+K+K+1 FLTA=MM+1 FLTB=-2*MM FLTC=MM FLTC=(FLTC-1.D0)*FLTC IF(NCON.GT.8)GOTO 60 H(I,I+1)=(D+FLTA*G)*SQRT((FLTB+Q)*Q+FLTC) GOTO 61 60 H(I,I+1)=(D+FLTA*G+FLT4*FKK4+FLT6*FKK6+FLT8*FKK8)* 1 SQRT((FLTB+Q)*Q+FLTC) 61 H(I+1,I)=H(I,I+1) 7 CONTINUE C C...touch up the top lh. corner elements IF(2-N)8,10,9 8 FLTA=-2*N+7 IF(NCON.LE.8)GOTO 55 H(1,1)=FLTA*(D+G+2.D0*(FKK4+FKK6+FKK8))*Q+H(1,1) GOTO 56 55 H(1,1)=FLTA*(D+G)*Q+H(1,1) 56 GOTO 10 9 H(1,2)=H(1,2)*1.414213562373095D0 H(2,1)=H(1,2) C C...diagonalise the Wang matrix, distribute eigenvalues and eigenvectors 10 CALL HDIAG0 LN=(J-J/2*2)*4+N DO 111 K=1,I E(LN,K)=H(K,K) DO 11 M=1,I T(LN,K,M)=R(K,M) H(K,M)=0.0D0 11 CONTINUE 111 CONTINUE 12 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ASROHS C C ASYMMETRIC ROTOR ENERGY LEVELS FOR A GIVEN J - Reduction 'S' C C I - on output contains the dimensionality of the calculated submatrix C J - J quantum number for the submatrix C N - number of the required submatrix = 1,2,3,4 for E+,E-,O+,O- resp. C A - array containing rotational and c.d. constants in usual order C C Evaluated submatrix is placed in H, then diagonalised and eigenvalues C transferred to E, eigenvectors to T C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=300, MAXDIM=MAXJ/2+1) C COMMON /ASRO/FMIN,FMAX,A(NCONST),EA,EB,EE,EPS, 1 S(3),E(8,MAXDIM),H(MAXDIM,MAXDIM),V(MAXDIM),R, 1 ID(MAXDIM),I,J,L,NT,IP COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),R(MAXDIM,MAXDIM) INTEGER*4 MM,NT C C...ancillary coefficients Q=J Q=Q*(J+1) QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C C=0.5D0*(A(2)+A(3)) DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ DK4=-A(6)+A(11)*Q+A(18)*QQ+A(27)*QQQ DK6=A(12)+A(19)*Q+A(28)*QQ DK8=A(20)+A(29)*Q DK10=A(30) C OKK2=0.25D0*(A(2)-A(3))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ OKK4=A(8)+A(14)*Q+A(22)*QQ+A(32)*QQQ OKK6=A(15)+A(23)*Q+A(33)*QQ OKK8=A(24)+A(34)*Q OKK10=A(35) C C...Process the four Wang matrices for the specified J C C (KA and KE are the lowest and highest K in the submatrix: C KA = 0, 2, 1, 1 for E+, E-, O+, O- resp.) C DO 120 N=1,4 K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K M=(KE-KA)/2+1 C DO 1 I=1,M DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE H(1,1)=0.D0 H(1,2)=0.D0 C C For summing over K loop variable KK equal to K+1 is used to avoid C loop count from zero for E+ C Variable I is a pointer to elements in the Wang submatrix C I=0 KA1=KA+1 KE1=KE+1 DO 3 KK=KA1,KE1,2 K=KK-1 I=I+1 MM=K MM=MM*K FK2=MM c c... H(I,I)=((((DK10*FK2+DK8)*FK2+DK6)*FK2+DK4)*FK2+DK2)*FK2 * +DK0 c c... IF(K.GE.KE)GOTO 3 H(I,I+1)=OKK2*DSQRT(FPLUS(J,K,2)) H(I+1,I)=H(I,I+1) c c... IF(K.GE.KE-2)GOTO 3 H(I,I+2)=OKK4*DSQRT(FPLUS(J,K,4)) H(I+2,I)=H(I,I+2) c c... IF(K.GE.KE-4)GOTO 3 H(I,I+3)=OKK6*DSQRT(FPLUS(J,K,6)) H(I+3,I)=H(I,I+3) c c... IF(K.GE.KE-6)GOTO 3 H(I,I+4)=OKK8*DSQRT(FPLUS(J,K,8)) H(I+4,I)=H(I,I+4) c c... IF(K.GE.KE-8)GOTO 3 H(I,I+5)=OKK10*DSQRT(FPLUS(J,K,10)) H(I+5,I)=H(I,I+5) c 3 CONTINUE C C Final touch up: 4 -> O+ and O-, 12 -> E-, 5 -> E+ C IF(2-N)4,12,5 c c...O+ and O-, FLTA = +1 for O+, -1 for O- 4 FLTA=-2*N+7 H(1,1)=FLTA*OKK2*DSQRT(FMINUS(J,1,2))+H(1,1) IF(KE.LT.3)GOTO 6 H(1,2)=FLTA*OKK4*DSQRT(FMINUS(J,3,4))+H(1,2) H(2,1)=H(1,2) IF(KE.LT.5)GOTO 6 H(2,2)=FLTA*OKK6*DSQRT(FMINUS(J,3,6))+H(2,2) H(1,3)=FLTA*OKK6*DSQRT(FMINUS(J,5,6))+H(1,3) H(3,1)=H(1,3) IF(KE.LT.7)GOTO 6 H(1,4)=FLTA*OKK8*DSQRT(FMINUS(J,7,8))+H(1,4) H(4,1)=H(1,4) H(2,3)=FLTA*OKK8*DSQRT(FMINUS(J,5,8))+H(2,3) H(3,2)=H(2,3) IF(KE.LT.9)GOTO 6 H(1,5)=FLTA*OKK10*DSQRT(FMINUS(J,9,10))+H(1,5) H(5,1)=H(1,5) H(2,4)=FLTA*OKK10*DSQRT(FMINUS(J,7,10))+H(2,4) H(4,2)=H(2,4) H(3,3)=FLTA*OKK10*DSQRT(FMINUS(J,5,10))+H(3,3) GOTO 6 c c...root(2) elements of E+ 5 IF(KE.LT.2)GOTO 6 K=4 IF(KE.LT.6)K=1+KE/2 DO 10 KK=2,K H(1,KK)=H(1,KK)*1.414213562373095D0 H(KK,1)=H(1,KK) 10 CONTINUE c c...E- and remaining elements of E+, FLTA = +1 for E+ and -1 for E- 12 IF(KE.LT.2)GOTO 6 FLTA=-2*N+3 K=2 IF(N.EQ.2)K=1 H(K,K)= FLTA*OKK4*DSQRT(FMINUS(J,2,4))+H(K,K) IF(KE.LT.4)GOTO 6 H(K,K+1)=FLTA*OKK6*DSQRT(FMINUS(J,4,6))+H(K,K+1) H(K+1,K)=H(K,K+1) IF(KE.LT.6)GOTO 6 H(K,K+2)=FLTA*OKK8*DSQRT(FMINUS(J,6,8))+H(K,K+2) H(K+2,K)=H(K,K+2) H(K+1,K+1)=FLTA*OKK8*DSQRT(FMINUS(J,4,8))+H(K+1,K+1) IF(KE.LT.8)GOTO 6 H(K,K+3)=FLTA*OKK10*DSQRT(FMINUS(J,8,10))+H(K,K+3) H(K+3,K)=H(K,K+3) H(K+1,K+2)=FLTA*OKK10*DSQRT(FMINUS(J,6,10))+H(K+1,K+2) H(K+2,K+1)=H(K+1,K+2) C C...diagonalise the Wang matrix, distribute eigenvalues and eigenvectors 6 CALL HDIAG0 LN=(J-J/2*2)*4+N DO 111 K=1,I E(LN,K)=H(K,K) DO 112 M=1,I T(LN,K,M)=R(K,M) H(K,M)=0.0D0 112 CONTINUE 111 CONTINUE 120 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FPLUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) integer*4 kdelta C C Values of f for matrix elements : C C dK-1 dK dK C + ___ ___ ___ C f (J,K) = | | (J-K-l) | | (J+K+l) = | | [J(J+1)-(K+l)(K+l-1)] C dK l=0 l=1 l=1 C FPLUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FPLUS=FPLUS*(RJ-RK-L) 1 CONTINUE DO 2 L=1,KDELTA FPLUS=FPLUS*(RJ+RK+L) 2 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FMINUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) integer*4 k,kdelta C C Values of f for matrix elements : C C dK-1 dK dK C - ___ ___ ___ C f (J,K) = | | (J+K-l) | | (J-K+l) = | | [J(J+1)-(K-l)(K-l+1)] C dK l=0 l=1 l=1 C FMINUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FMINUS=FMINUS*(RJ+RK-L) 1 CONTINUE DO 2 L=1,KDELTA FMINUS=FMINUS*(RJ-RK+L) 2 CONTINUE C RETURN END C_____________________________________________________________________________ C SUBROUTINE HDIAG0 C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, MAXJ=300, MAXDIM=MAXJ/2+1) C COMMON /ASRO/FMIN,FMAX,A(NCONST),EA,EB,EE,SPE, 1 D(3),E(8,MAXDIM),H(MAXDIM,MAXDIM),X(MAXDIM),U, 1 IQ(MAXDIM),N,M,L,NT,IP COMMON /BIG/T C REAL*4 T(8,MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),RSIN,RCOS,RHTEMP INTEGER*4 NT C IF(IP.EQ.3)GOTO 15 DO 14 I=1,N DO 14 J=1,N IF(I-J)12,11,12 11 U(I,J)=1. GOTO 14 12 U(I,J)=0. 14 CONTINUE 15 NR=0 IF(N-1)1000,1000,17 17 NMI1=N-1 DO 30 I=1,NMI1 X(I)=0.D0 IPL1=I+1 DO 30 J=IPL1,N IF(X(I)-ABS(H(I,J)))20,20,30 20 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE 40 DO 70 I=1,NMI1 IF(I-1)60,60,45 45 IF(XMAX-X(I))60,70,70 60 XMAX=X(I) IPIV=I JPIV=IQ(I) 70 CONTINUE EPS=1.0D-12 IF(XMAX-EPS)1000,1000,120 120 NR=NR+1 IF(H(IPIV,IPIV)-H(JPIV,JPIV))140,141,142 140 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 -SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) GOTO 150 141 TANG=1.D0 GOTO 150 142 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 +SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) 150 COSINE=1.D0/SQRT(1.D0+TANG**2) SINE=TANG*COSINE HII=H(IPIV,IPIV) H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.D0*H(IPIV,JPIV)+TANG* 1 H(JPIV,JPIV))) H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.D0*H(IPIV,JPIV)- 1 TANG*HII)) H(IPIV,JPIV)=0.D0 IF(H(IPIV,IPIV)-H(JPIV,JPIV))153,153,152 152 HTEMP=H(IPIV,IPIV) H(IPIV,IPIV)=H(JPIV,JPIV) H(JPIV,JPIV)=HTEMP SUNO=-SINE HTEMP=SIGN(1.D0,SUNO)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV)230,240,230 230 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) HTEMP=H(I,K) H(I,K)=0.D0 IPL1=I+1 X(I)=0.D0 DO 320 J=IPL1,N IF(X(I)-ABS(H(I,J)))300,300,320 300 X(I)=ABS(H(I,J)) IQ(I)=J 320 CONTINUE H(I,K)=HTEMP 350 CONTINUE X(IPIV)=0.D0 X(JPIV)=0.D0 DO 530 I=1,N IF(I-IPIV)370,530,420 370 HTEMP=H(I,IPIV) H(I,IPIV)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(I)-ABS(H(I,IPIV)))380,390,390 380 X(I)=ABS(H(I,IPIV)) IQ(I)=IPIV 390 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(I)-ABS(H(I,JPIV)))400,530,530 400 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 420 IF(I-JPIV)430,530,480 430 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(IPIV)-ABS(H(IPIV,I)))440,450,450 440 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 450 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(IPIV)-ABS(H(I,JPIV)))460,530,530 460 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 480 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(JPIV,I) IF(X(IPIV)-ABS(H(IPIV,I)))490,500,500 490 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 500 H(JPIV,I)=-SINE*HTEMP+COSINE*H(JPIV,I) IF(X(JPIV)-ABS(H(JPIV,I)))510,530,530 510 X(JPIV)=ABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF(IP.EQ.3)GOTO 40 RSIN=SINE RCOS=COSINE DO 550 I=1,N RHTEMP=U(I,IPIV) U(I,IPIV)=RCOS*RHTEMP+RSIN*U(I,JPIV) U(I,JPIV)=-RSIN*RHTEMP+RCOS*U(I,JPIV) 550 CONTINUE GOTO 40 1000 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE DIPTRA(DIDUTY,RJS,ISON,KMMAX,iexit) C C Calculate and determine acceptable transitions between levels currently C being worked on. Write these out to an intermediate file for later C sorting. Space saving coding of quantum numbers from the original ASROT C is still retained (though this, as well as writing to an intermediate c file could be dispensed with as the days of memory limitations are gone) C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) INTEGER*4 MAXLIN,NT,NREC PARAMETER (NCONST=35, MAXJ=300, MAXDIM=MAXJ/2+1, * MAXLIN=MAXDIM*MAXDIM*4) C COMMON /ASRO/FMIN,FMAX,A(NCONST),EA,EB,EE,EPS, 1 D(3),E(8,MAXDIM),H(MAXDIM,MAXDIM),V(MAXDIM),U, 1 ID(MAXDIM),I,J,JS,NT,IP COMMON /DEF/NREC COMMON /BIG/T C REAL*8 DIDUTY(18) REAL*4 T(8,MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),EASP,ELOW C 90 FORMAT(/1X) 91 FORMAT(1X,3I4,4X,3I4,I9,F21.6,F20.6) C C GF(K)=SIGN(1,L-N3+1)*SQRT(DBLE(REAL(((K+1)*N3+J)*(K*N3+J) 1 *(1+1/(K+1))))) NF(K)=2*K-(K+1)/2*4+1 JP=0 IRQ=1 ID(1)=J I=(J/2*2-J)*2+1 ID(12)=(I+1)*2 ID(13)=(1-I)*2 DO 80 I=1,18 IF(I.NE.13)GOTO 3 JP=0 IRQ=0 3 K=MOD(RJS,2.D0) RJS=(RJS-K)*0.5D0 IF(K.EQ.0)GOTO 80 XXX=DIDUTY(I) DO 5 K=3,11 XXXX=MOD(XXX,8.D0) ID(K)=XXXX+1.D-6 XXX=(XXX-XXXX)*.125D0 5 CONTINUE IS=XXX+1.D-6 M1=(J+2-IRQ-ID(8))/2 IF(M1.EQ.0)GOTO 80 M2=(J+2-ID(9))/2 IF(M2.EQ.0)GOTO 80 I1=ID(5)+ID(13-IRQ) I2=ID(6)+ID(13) L=ID(4)-1 N=ID(3) ID(4)=J-IRQ GOTO(10,20,30,40,50),IS 10 M=(J-1)/2*2 NN=N+1 MM=M+1 DO 11 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=NF(K) H(N1,N2)=GF(K) 11 CONTINUE GOTO 60 20 M=J/2*2-1 NN=N+1 MM=M+1 DO 21 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=-NF(K) H(N2,N1)=GF(K) 21 CONTINUE GOTO 60 30 DO 31 K=1,M1 G=J G=G*J-((K-1)*2+ID(8))**2 H(K,K)=2.D0*SQRT(G) 31 CONTINUE GOTO 60 40 M=J-1 NN=N+1 MM=M+1 DO 41 KKK=NN,MM K=KKK-1 N1=(K-2*N+3)/2 N2=K/2+1 N3=NF(K)+L+1 DBLA=REAL(J+1)*J-REAL(K+1)*K DBLA=DBLA*(1+1/(K+1)) H(N1,N2)=SIGN(1,int(N3))*SQRT(DBLA) 41 CONTINUE GOTO 60 50 DO 51 K=1,M2 N1=K-N+1 H(N1,K)=((K-1)*2+ID(9))*2 51 CONTINUE 60 G=((1-IRQ)*(2*J+1)/(J+1.D0)+IRQ)*0.25D0/J DO 169 K=1,M2 N1=0 DO 69 L=1,M1 EE=E(I2,K)-E(I1,L) EB=ABS(EE) IF(EB.LT.FMIN.OR.EB.GT.FMAX)GOTO 69 ELOW=E(I2,K) IF(E(I1,L).LT.E(I2,K))ELOW=E(I1,L) ELOW=ELOW*3.3356409E-5 IF(N1.NE.0)GOTO 64 N1=1 DO 63 M=1,M1 V(M)=0.D0 DO 63 N=1,M2 V(M)=T(I2,N,K)*H(M,N)+V(M) 63 CONTINUE 64 EA=0.D0 DO 65 M=1,M1 EA=T(I1,M,L)*V(M)+EA 65 CONTINUE N2=ID(7) EA=EA**2*G IF(EA.LT.EPS)GOTO 69 N2=(L-1)*2 N3=N2+ID(8) ID(5)=N3 IF(ID(5).GT.KMMAX)GOTO 69 ID(6)=ID(4)-N3+ID(10) N2=(K-1)*2 N3=N2+ID(9) ID(2)=N3 IF(ID(2).GT.KMMAX)GOTO 69 ID(3)=J-N3+ID(11) IF(JP.NE.0)GOTO 68 IF(ISON.EQ.0)WRITE(4,90) 68 JP=1 IF(ISON.EQ.0)WRITE(4,91)(ID(KSQ),KSQ=1,7),EE,EA EASP=EA IF(NT.GE.MAXLIN)THEN if(iexit.eq.0)write(*,'(1x//)') WRITE(*,'('' J='',I3,'' <- maximum of'',i6, * '' lines has been reached'',a)')j,maxlin,char(7) iexit=1 GOTO 269 ENDIF C C...Compression of quantum numbers is simplified from previous version to c avoid problems on possible increase in J. The same amount of space c (8 bytes) is now used to store quantum numbers of the lower state and c packed changes in quantum numbers c C (DeltaK-1+10)*1000 + (DeltaK+1+10) + DeltaJ + 1 c c so that a.R0,1 is 10112, bQ1,-1 is 11091. There might be problems if c DeltaK exceeds +-9 c C ID(1),ID(2),ID(3) - J,K-1,K+1 for upper state C ID(4),ID(5),ID(6) - J,K-1,K+1 for lower state C NT=NT+1 IF(EE.GT.0.0d0)THEN ID(21)=ID(4) ID(22)=ID(5) ID(23)=ID(6) ID(24)=(ID(2)-ID(5)+10)*1000+(ID(3)-ID(6)+10)*10+ * (ID(1)-ID(4)+1) ELSE ID(21)=ID(1) ID(22)=ID(2) ID(23)=ID(3) ID(24)=(ID(5)-ID(2)+10)*1000+(ID(6)-ID(3)+10)*10+ * (ID(4)-ID(1)+1) ENDIF NREC=NT+1 WRITE(2,REC=NREC)ID(21),ID(22),ID(23),ID(24), * EB,EASP,ELOW cdump WRITE(7,'(3I3,2X,3I3,2X,3I3,I12,F20.2)') cdump * (ID(KKK),KKK=1,6),(ID(KKK),KKK=21,24),EE 69 CONTINUE WRITE(2,REC=1)NT 169 CONTINUE c 269 WRITE(2,REC=1)NT H(1,1)=0.D0 N3=J/2+1 DO 70 K=2,N3 H(K,K)=0.D0 H(K,K-1)=0.D0 H(K-1,K)=0.D0 70 CONTINUE 80 CONTINUE C C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SORT0(FV,SIGMA,GI,HALFW,ABUND,TEMP) C C Read and decode transitions stored in an intermediate SCRATCH file and C sort these in ascending order of frequency. C Arrays unpacked to are F,EL,S,NQ,ITRANS,M C C NQ(I,1),NQ(I,2),NQ(I,3) - J,K-1,K+1 of the upper state C NQ(I,4),NQ(I,5),NQ(I,6) - J,K-1,K+1 of the lower state C NQ(I,7) - type of transition flag = 1,2,3 for a,b,c C NQ(I,8) - Delta K-1 C NQ(I,9) - Delta K+1 C NQ(I,10) - (Delta J) + 2 C C To save a bit of space (for maxdim=151, maxlin=91204): c c - the eigenvector storage array R*4 T(8,maxdim,maxdim) in /BIG/ C is here reused as R*8 F(maxlin) c c - two other big arrays R*8 H(MAXDIM,MAXDIM)=R*4(maxlin/2) in /ASRO/ c and the eigenvector array R*4 R(MAXDIM,MAXDIM)=I*2(maxlin/2) c in /ASRO/ are no longer reused C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 MAXDIM,maxlin c c MAXJ=300: MAXDIM=151, MAXLIN=91204 c PARAMETER (NCONST=35, MAXJ=300, MAXDIM=MAXJ/2+1, * MAXLIN=MAXDIM*MAXDIM*4) COMMON /ASRO/FMIN,FMAX,A(NCONST),EA,EB,EE,EPS, 1 D(3),E(8,MAXDIM),H(MAXDIM,MAXDIM),V(MAXDIM),R, 1 ID,B,C,LL,N,IP COMMON /DEF/NREC COMMON /BIG/F C REAL*4 R(MAXDIM,MAXDIM) REAL*8 F(maxlin),WWK INTEGER*2 B,C,NQ(maxlin,10),ID(MAXDIM),IP,LL INTEGER*4 ITRANS(maxlin),M(MAXLIN),N,NREC,I, * LSORT,IRSORT,ISORT,JSORT,IIPT,NOUTP CHARACTER DIPOLE(3),TYPETR(3) REAL*4 EL(maxlin),ALFA,FV,SIGMA,GI,HALFW,ABUND,TEMP, 1 CONST,CONST1,CONST3,FREQ,S(maxlin) c DATA DIPOLE/'a','b','c'/TYPETR/'P','Q','R'/ FRESOL=0.001D0 C 20 FORMAT(/1x,78(1H-)/' LINES IN ORDER OF ASCENDING FREQUENCY:', 1 10X,I5,' lines in temp. file' 1 //' FREQUENCY', 1 4X,'Alpha max',4X,'J K K J K K',4X,'Branch',4X, 1 'Line',8X,'Lower'/34X,'-1 +1 -1 +1',12X,'Strength', 1 4X,'Level'/4X,'/MHz',9X,'/cm-1',52X,'cm-1'/) 2500 FORMAT(//' *****',i6,' LINES ARE BEING SORTED') 2501 FORMAT(1X) C C C...Read unsorted frequencies: NQ(i,4 to 6) - quantum numbers of lower C state, NQ(i,7) - packed changes in quantum numbers for the transition, C F - frequencies, S - line strengths, EL - energy of the lower level C READ(2,REC=1)N WRITE(*,2500)N DO 200 I=1,N NREC=I+1 READ(2,REC=NREC)NQ(i,4),NQ(i,5),NQ(i,6),NQ(i,7), * F(I),S(I),EL(I) cdump WRITE(7,'(3I3,I7,F20.2)')NQ(i,4),NQ(i,5),NQ(i,6),NQ(i,7), cdump * F(I) 200 CONTINUE c C Unpack quantum number changes, calculate quantum numbers of upper state, c dipole type and generate transition type index ITRANS C - this is used for identification of unresolved K+1 or K-1 doublets C DO 10 I=1,N C NQ(I,1)=NQ(I,4)+MOD(NQ(I,7),10)-1 NQ(I,2)=NQ(I,5)+NQ(I,7)/1000-10 NQ(I,3)=NQ(I,6)+MOD(NQ(I,7),1000)/10-10 C NQ(I,8)=NQ(I,2)-NQ(I,5) NQ(I,9)=NQ(I,3)-NQ(I,6) NQ(I,10)=NQ(I,1)-NQ(I,4)+2 C c...(DeltaK-1,DeltaK+1) is (e,o) for a-type c (o,o) for b-type c (o,e) for c-type NQ(I,7)=MOD(IABS(int(NQ(I,8))),2)-MOD(IABS(int(NQ(I,9))),2)+2 C ITRANS(I)=(NQ(I,7)*10+NQ(I,10))*10 ITRANS(I)=( (ITRANS(I)+(NQ(I,8)+50))*100+NQ(I,9) )+50 cdump WRITE(7,'(3I3,2X,3I3,I5,I3,2I3,F20.2)') cdump * (nq(I,KKK),KKK=1,10),F(I) 10 CONTINUE c C C...Sorting loop: c c The code is adapted from rhe SORT2 'heapsort' routine from Numerical c Recipes: quantities in vector F are sorted in ascending order of magnitude c and the vector M of pointers to original positions of sorted quantities c is rearranged accordingly. NQ, S and EL remain unsorted C WRITE(4,20)N DO 1 I=1,N M(I)=I 1 CONTINUE c LSORT=N/2+1 IRSORT=N 50 CONTINUE IF(LSORT.GT.1)THEN LSORT=LSORT-1 WWK=f(LSORT) IIPT=m(LSORT) ELSE WWK=f(IRSORT) IIPT=m(IRSORT) f(IRSORT)=f(1) m(IRSORT)=m(1) IRSORT=IRSORT-1 IF(IRSORT.EQ.1)THEN f(1)=WWK m(1)=IIPT GOTO 52 ENDIF ENDIF ISORT=LSORT JSORT=LSORT+LSORT 51 IF(JSORT.LE.IRSORT)THEN IF(JSORT.LT.IRSORT)THEN IF(f(JSORT).LT.f(JSORT+1))JSORT=JSORT+1 ENDIF IF(WWK.LT.f(JSORT))THEN f(ISORT)=f(JSORT) m(ISORT)=m(JSORT) ISORT=JSORT JSORT=JSORT+JSORT ELSE JSORT=IRSORT+1 ENDIF GO TO 51 ENDIF f(ISORT)=WWK m(ISORT)=IIPT GO TO 50 52 CONTINUE c C C...Calculate intensities, establish types of transitions and C write out sorted transitions C CONST=3.85E-14*DSQRT(A(1)*A(2)*A(3))*FV*GI*SIGMA/(HALFW 1 *TEMP**2.5)*ABUND CONST1=0.5/(29979*0.69497) IOMIT=0 NOUTP=1 C C DO 6 I=1,N IF(IOMIT.EQ.1)GOTO 11 J=M(I) JP1=M(I+1) C C...absorption coefficient CONST3=-EL(J)/(0.69497*TEMP) IF(CONST3.LT.-75.)GOTO 6 FREQ=F(I) ALFA=CONST*SNGL(D(NQ(J,7)))**2*S(J)*FREQ**2*EXP(CONST3) 1 *(1.-CONST1*FREQ/TEMP) IF(ALFA.LT.1.E-12)GOTO 6 C C...check for unresolved doublets IF(I.EQ.N)GOTO 12 IF(DABS(F(I)-F(I+1)).GT.FRESOL)GOTO 12 IF(ITRANS(J).NE.ITRANS(JP1))GOTO 12 IF((NQ(J,1).NE.NQ(JP1,1)).OR.(NQ(J,4).NE.NQ(JP1,4)))GOTO 12 IF((NQ(J,2).EQ.NQ(JP1,2)).AND.(NQ(J,5).EQ.NQ(JP1,5)))GOTO 14 IF((NQ(J,3).EQ.NQ(JP1,3)).AND.(NQ(J,6).EQ.NQ(JP1,6)))GOTO 15 GOTO 12 C C...doublet with K-1's the same 14 ALFA=ALFA*2. F(I)=0.5D0*(F(I)+F(I+1)) IOMIT=1 WRITE(4,17)F(I),ALFA,NQ(J,1),NQ(J,2),NQ(J,4),NQ(J,5), 1 DIPOLE(NQ(J,7)), 1 TYPETR(NQ(J,10)),NQ(J,8),NQ(J,9),S(J),EL(J) 17 FORMAT(F14.6,1PE10.2,' D',I3,',',I3,',',I8,',',I3,',', 1 4X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C C...doublet with K+1's the same 15 ALFA=ALFA*2. F(I)=0.5D0*(F(I)+F(I+1)) IOMIT=1 WRITE(4,16)F(I),ALFA,NQ(J,1),NQ(J,3),NQ(J,4),NQ(J,6), 1 DIPOLE(NQ(J,7)), 1 TYPETR(NQ(J,10)),NQ(J,8),NQ(J,9),S(J),EL(J) 16 FORMAT(F14.6,1PE10.2,' D',I3,', ,',I3,I5,', ,',I3, 1 1X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C C...single transition 12 WRITE(4,21)F(I),ALFA,(NQ(J,K),K=1,6),DIPOLE(NQ(J,7)), 1 TYPETR(NQ(J,10)),NQ(J,8),NQ(J,9),S(J),EL(J) 21 FORMAT(F14.6,1PE10.2,2X,2(I3,','),I3,I5,2(',',I3), 1 1X,A1,'.',A1,I2,',',I2,0PF10.6,F10.3) IF((NOUTP/10)*10.EQ.NOUTP)WRITE(4,2501) NOUTP=NOUTP+1 GOTO 6 C 11 IOMIT=0 C 6 CONTINUE C C C...Concluding information C WRITE(*,2502)NOUTP-1 2502 FORMAT(i12,' LINES WRITTEN TO FILE') WRITE(4,666)NOUTP-1,FV,SIGMA,GI,HALFW,ABUND,TEMP 666 FORMAT(1X/49X,I5,' lines printed' 1 //' ABSORPTION COEFFICIENTS WERE CALCULATED WITH:'// 1 7X,' Vibrational population factor, F.v = ',F9.4/ 1 19X,' Symmetry number, sigma = ',F9.4/ 1 2X,' Reduced nuclear statistical weight, g.I = ',F9.4/ 1 ' Pressure broadening parameter /MHz Torr-1 = ',F9.4/ 1 11X,' Isotopic abundance factor, i.c = ',F9.4/ 1 25X,' Temperature, T/K = ',F9.4/1x,78(1H-)) C RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________ C C Routines in the order in this listing: C C MP C ASROHA - A-reduction energy levels C ASROHS - S-reduction energy levels C FPLUS - evaluate the Fplus function for the Hamiltonian C FMINUS - evaluate the Fminus function for the Hamiltonian C HDIAG0 - diagonalise the Hamiltonian matrix C DIPTRA - determine and calculate transitions between energy levels C SORT0 - read from scratch file, decode and sort transitions C C_____________________________________________________________________________ C_____________________________________________________________________________