c$debug C----------------------------------------------------------------------------- c c P I S O R T - PIckett Sort (Sorting of the .OUT file from H.M.Pickett's C program SPCAT) C c The purpose of this program is to produce a data file c which can be displayed in stick diagram form by c programs ASC and ASCP (the output is compatible with c that from program ASROT) c C ver.16b.IV.2003 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c c Sorting of output consisting of six quantum numbers per energy level. c c - the 12 quantum numbers of .OUT file can optionally be redirected into c any order in the .ASR file, the redirection information is read from c file PISORT.RED c PISORT always looks for this file (in the default directory), if it c is not found then default ordering is assumed, ie. same as in .OUT c file c c - If the first three quantum numbers for each state are the c asymmetric top quantum numbers then PISORT will correctly assign c the dipole transition flag: a,b, or c c c - for sorting of transitions into vibrational states it is assumed that c the fourth quantum number for each state is v and can take c one of three values c c - intensity calculated by SPCAT (CALCAT) for 300 K can be recalculated c for a different temperature and into alpha.max at the same conditions c of calculations as empoyed by ASROT c c - the upper limit for setting MAXLIN depends on the particular FORTRAN c being used, specifically the limits imposed on size of character c arrays c C----------------------------------------------------------------------------- c C Modification history: C c 28.09.96: recalculation of intensity for temperature different from 300K c 2.01.97: redirection of quantum numbers c 24.02.97: status='OLD' on input and steps against overoptimisation of nwrit c 14.03.97: decoding of notation for qn's smaller than -9 c 23.11.97: scrapping of all INTEGER*2 for higher MAXLIN c 27.11.98: decoding of notation for qn's>=100 c 24.10.99: check against 1999 C-version of SPCAT + small mods c 10.02.01: modification to work with the 12a2 compressed quantum number c field beginning in the 76th column of .OUT files c 23.10.01: modifications to work with millions of lines c 24.10.02: filtering out of weakest lines + linestrength output for .ASR c 16.04.03: addition of astrophysical output format c C----------------------------------------------------------------------------- c parameter (maxlin=1000000) c character*108 L(maxlin),line,linem,linea real*8 f(maxlin),scut,slin,strlin, * slimit,elimit,ferror,elevel,dipb,freql integer i(maxlin),nstart,nlines integer iredir(12) character*30 filin,filout character cflag*2,lstr*13,level*10 common /lbuf/L common /freq/f common /point/i common /redir/iredir c scut=1.d-15 c WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | PISORT - Sorting of .OUT file from Pickett''s ', * 'program SPCAT',T79,'|'/ * ' |',76(1H_),'|'/' version 16b.IV.2003',T64,'Zbigniew KISIEL'/) write(*,'(1x/'' - compiled for maximum of'',i9,'' lines'')') * maxlin c C...open redirection file PISORT.RED, if found c open(4,file='PISORT.RED',status='OLD',err=50) 52 read(4,'(a)',err=50,end=50)line if(line(1:1).eq.';')goto 52 backspace(4) read(4,*,err=50)iredir goto 2 c 50 write(*,'(1x/'' - PISORT.RED not found, quantum numbers will'', * '' be in the default order''/ * '' ie. 1,2,3,7,8,9,4,5,6,10,11,12'')') do 51 n=1,12 iredir(n)=n 51 continue do 551 n=4,9 if(n.lt.7)iredir(n)=n+3 if(n.ge.7)iredir(n)=n-3 551 continue c c...open .OUT file c 2 write(*,1)' Name of input .OUT file (without extension):' 1 format(1x//1x,a,' ',$) read(*,'(a)',err=2)filin do 100 n=30,1,-1 if(filin(n:n).ne.' '.and.ichar(filin(n:n)).ne.0)then filin=filin(1:n)//'.out' goto 101 endif 100 continue 101 write(*,'(1x/'' R E A D I N G: '',a)')filin open(3,file=filin,status='OLD',err=2) c c...Loop for fishing out spectral lines from the .OUT file c nignor=0 53 nlines=0 5 read(3,'(a)',end=6,err=5)line if(line(9:9).ne.'.'.or.line(18:18).ne.'.')goto 5 nlines=nlines+1 L(nlines)=line read(line(1:13),'(f13.4)',err=8)f(nlines) i(nlines)=nlines if(nlines.eq.maxlin)goto 6 goto 5 8 nlines=nlines-1 nignor=nignor+1 goto 5 c 6 write(*,7)nlines,nignor 7 format(1x/i9,' lines read in'/ * i9,' lines appearing to begin with frequency ignored'/) close(3) c c...frequency sort c write(*,'(1x/'' S O R T I N G'')') nstart=1 call sorth(nstart,nlines) c c...output options c 14 write(*,15) 15 format(1x/' OUTPUT OPTIONS -- - exit'/ * ' aa - all'/ * ' 00 - state 0'/ * ' 11 - state 1'/ * ' 22 - state 2'/ * ' 10 - inter state 1<-0'/ * ' 01 - inter state 0<-1'// * ' . . . . . ',$) read(*,'(a2)')cflag if(cflag.ne.'00'.and.cflag.ne.'11'.and.cflag.ne.'22'.and. * cflag.ne.'--'.and.cflag.ne.'10'.and.cflag.ne.'01'.and. * cflag.ne.'aa'.and.cflag.ne.'AA')goto 14 if(cflag.eq.'--')goto 9 c 22 write(*,21) 21 format(1x/21x,'-1 = astrophysical format'/ * 21x,' 0 = unaltered .OUT format'/ * 21x,' 1 = ASROT format'// * 21x,' . . . . . ',$) read(*,'(i5)',err=22)iasrot c if(iasrot.lt.-1.or.iasrot.gt.1)goto 22 c if(iasrot.eq.1)then 26 write(*,25) 25 format(1x// * ' T E M P E R A T U R E /K for the output spectrum:' * //10x,'- it is assumed that the calculation was for 300K' * /10x,'- Temp <= 0 specifies no temperature conversion'// * 10x,'. . . . . ',$) read(*,'(F15.5)',err=26)temp endif if(temp.eq.0.d0)temp=300.d0 c c...output headers c 3 write(*,1)' Name of output file:' read(*,'(a)',err=3)filout open(4,file=filout,err=3,status='unknown') c if(iasrot.eq.-1)then write(4,210)filin endif 210 format(79(1h-)/' Converted from the SPCAT output file: ',a/79(1h-) * //3x,'Upper Lower Frequency Error', * ' Linestr. G Upper' * / 3x,'quanta quanta ', * ' * muG**2 Level' * /29x, '(MHz) (MHz)', * ' (D**2) (K) '/) C if(iasrot.eq.0)then write(4,200)filin write(4,203) endif 200 format(108(1h-)/' Unchanged lines echoed from the SPCAT ', * 'output file: ',a// * ' sorted according to increasing frequency'/108(1h-)/) 203 format(1x/' FREQUENCY LStr.DIP**2 DR ', * ' E.low GUP QNFORM Upper Lower'/ * ' EST.ERROR. LGSTR. | /cm-1 ', * ' | TAG | quanta quanta'/ * T46,'| | | |'/ * T46,'| | | | 1 2 3 4 5 6', * ' 1 2 3 4 5 6'/) c if(iasrot.eq.1)then write(4,201)filin,temp write(4,202) endif 201 format(93(1h-)/' Converted from the SPCAT output file: ',a/93(1h-) * //' Temperature converted from assumed 300K to',f7.1,'K'// * ' Peak absorption coefficient calculated assuming: ', * ' F.v = 0.95'/T43,'Delta.nu at 1 Torr = 20 MHz' * /93(1h-)/) 202 format(1x/' FREQUENCY abs.coeff. upper <- lower', * ' G upper <- lower E.low Linestr.'/ * ' /MHz /cm-1 quanta quanta', * ' quanta quanta /cm-1 * muG**2'/) c nwrit=0 nwrit1=0 nwrit2=0 nwrit3=0 nwrit4=0 nwrit5=0 C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...deal with astrophysical output by using ASROT conversion to decode quantum C numbers. C C ASROT columns: C QNS 27-50 and 30,34,43,47 are comas C frequency 1-12 C dipole flag 52 C lower level 72-79 C C OUT columns: C error 14-22 C if(iasrot.eq.-1)then slimit=1.d0 elimit=1000.d0 write(*,45)slimit,elimit 45 format(1x//' Limits on output:'// * f15.2,' on Linestrength*dip**2'/ * f15.2,'K on Upper level energy'/) do 41 n=1,nlines ii=i(n) line=L(ii) linem=L(ii) call casrot(line,temp) line(30:30)=' ' line(34:34)=' ' line(43:43)=' ' line(47:47)=' ' read(L(ii)(23:35),'(e13.6)')strlin read(Line(72:79),'(f8.3)')elevel read(Line( 1:12),'(f12.4)')freql elevel=elevel+freql/29979.2458d0 elevel=elevel*1.43877d0 if(strlin.lt.slimit)goto 41 write(lstr,'(F10.3)')strlin if(elevel.gt.elimit)goto 41 write(level,'(F10.3)')elevel c write(4,'(a)')line(27:50)//' '//line(1:12)//L(ii)(14:22)// * lstr(1:10)//' '//line(51:52)//level nwrit=nwrit+1 41 continue write(4,45)slimit,elimit goto 40 endif C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...deal with .OUT and ASROT output C c c...first check whether quantum numbers are in the old (three character fields) C or new (two character fields) format c nva=84 nvb=102 if(l(1)(84:84).eq.' ')then nva=83 nvb=95 endif c do 12 n=1,nlines ii=i(n) line=L(ii) if(iasrot.ne.1)then read(line(36:44),'(f9.4)')slin slin=10.d0**slin endif c if(iasrot.eq.1)then call casrot(line,temp) read(line(17:24),'(e8.2)')slin read(L(ii)(23:35),'(e13.6)')strlin write(lstr,'(F13.7)')strlin endif c if(cflag.eq.'aa'.or.cflag.eq.'AA')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit=nwrit+1 goto 12 endif c if(cflag.eq.'22'.and.L(ii)(nva:nva).eq.'2'.and. * L(ii)(nvb:nvb).eq.'2')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit1=nwrit1+1 goto 12 endif c if(cflag.eq.'11'.and.L(ii)(nva:nva).eq.'1'.and. * L(ii)(nvb:nvb).eq.'1')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit2=nwrit2+1 goto 12 endif c if(cflag.eq.'00'.and.L(ii)(nva:nva).eq.'0'.and. * L(ii)(nvb:nvb).eq.'0')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit3=nwrit3+1 goto 12 endif c if(cflag.eq.'01'.and.L(ii)(nva:nva).eq.'0'.and. * L(ii)(nvb:nvb).eq.'1')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit4=nwrit4+1 goto 12 endif c if(cflag.eq.'10'.and.L(ii)(nva:nva).eq.'1'.and. * L(ii)(nvb:nvb).eq.'0')then if(iasrot.eq.1)then if(slin.gt.scut)write(4,'(a)')line(1:79)//lstr else if(slin.gt.scut)write(4,'(a)')line endif if(slin.gt.scut)nwrit5=nwrit5+1 goto 12 endif c 12 continue C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c C...concluding output lines c 40 nwrit=nwrit+nwrit1+nwrit2+nwrit3+nwrit4+nwrit5 C write(*,20)nwrit if(iasrot.eq.-1)write(4,35)nwrit if(iasrot.eq. 0)write(4,31)nwrit if(iasrot.eq. 1)write(4,30)nwrit 20 format(1x/1x,i8,' lines written to the output file'/) 30 format(1x/1x,i8,' lines written to this file'//93(1H-)/93(1h-)) 31 format(1x/1x,i8,' lines written to this file'//108(1H-)/108(1h-)) 35 format(1x/1x,i8,' lines written to this file'//79(1H-)/79(1h-)) C close(4) c goto 14 c 9 stop end C_____________________________________________________________________________ C subroutine casrot(line,temp) c c Routine to convert a line of .OUT output into a line equivalent c to .ASR output c c typical ASROT line is: c c 25277.238471 0.99E-07 116, 7,110 115, 6,109 b.R 1, 1 99.048770 543.614 c c but the information past the dipole letter is not read by ASC. c PISORT fills out the lower level energy but uses the rest of the c field past the dipole flag for encoding the remaining six PICKETT quantum c numbers eg. for interstate v2=1<-v3=1 transition with both Br and N coupling c in the F1,F scheme: c C117444.1572 2.86E-04 16, 1, 15 15, 2, 14 b 1,16,15-- 0,15,14 40.144 c c TEMP = temperature (K) c character*108 line,linasr character confld*36 integer iredir(12),nq(12) character*2 cdip real*8 s common /redir/iredir c linasr(1:12)=line(2:13) linasr(13:14)=' ' c c...Convert the SPCAT intensity, which is in units of (nm**2 MHz): c c 1/ to alpha.max (cm-1) calculated at conditions used in the ASROT c calculation c 2/ to a different temperature (it is assumed that in the .int file the c value 5.34E+6 (300**3/(A B C) )**(1/2) has been entered for Q.r c (A,B,C in MHz), and that the SPCAT alculation was for the default c temperature of 300K c c The connection between the SPCAT intensity and alpha.max is given in the c documentation of the jpl catalog and is: c alpha.max /cm-1 = 102.458 I(SPCAT)/Delta.nu c where Delta.nu is the HWHH in MHz at 1 Torr. Since ASROT uses Delta.nu=20, c and F.v=0.95 then c alpha.max (ASROT) = 4.87 I(SPCAT) c c The temperature conversion is according to eq.7.73 3rd.ed G&C, but with c accounting for the (1-hv/2kT) term dropped in G&C from Eq.3.27 c read(line(36:44),'(f9.4)',err=110)s s=10.d0**s c if(nint(temp).ne.300.and.temp.gt.0.)then read(line(48:57),'(F8.3)',err=112)EROT read(line(1:13),'(F13.4)',err=112)F fnk=F/(29979.2458d0*0.6950356D0) s=S*(1558845.7D0/TEMP**2.5)* * dexp(-erot/(0.6950356D0*TEMP))/dexp(-erot/208.51068D0)* * (1.d0-0.5d0*fnk/temp)/(1.d0-0.5d0*fnk/300.d0) s=s*4.87d0 endif c write(linasr(15:24),'(1PE10.2)')s c c...The twelve quantum numbers are in 12a3 format in columns: c c 73:75 76:78 79:81 82:84 85:87 88:90 c 91:93 94:96 97:99 100:102 103:105 106:108 c c c...In later versions the twelve quantum numbers are in 12a2 format in c columns: c c 76:77 78:79 80:81 82:83 84:85 86:87 c 88:89 90:91 92:93 94:95 96:97 98:99 c c These versions can be identified by testing the character in column 75 - c it should be a space, whereas in previous versions it is the least c significant character of the first quantum number. If this case is c encountered then pad out quantum numbers to the earlier 12a3 format c beginning in column 73: c if(line(75:75).eq.' ')then confld=line(75:77)//' '//line(78:79)//' '//line(80:81) confld=confld(1:9)//' '//line(82:83)//' '//line(84:85)// * ' '//line(86:87) confld=confld(1:18)//' '//line(88:89)//' '//line(90:91)// * ' '//line(92:93) confld=confld(1:27)//' '//line(94:95)//' '//line(96:97)// * ' '//line(98:99) line=line(1:72)//confld endif c c...decode two digit compression of negative 2-digit qn's c decode two digit compression of positive three digit qn's c call check(line,73,75) call check(line,76,78) call check(line,79,81) call check(line,82,84) call check(line,85,87) call check(line,88,90) call check(line,91,93) call check(line,94,96) call check(line,97,99) call check(line,100,102) call check(line,103,105) call check(line,106,108) c read(line(73:108),'(12i3)',err=100)(nq(i),i=1,12) goto 101 c 100 write(*,111)'Cannot read quantum numbers from:', * 73,108,line(73:108),line(1:60),line(61:108) 111 format(1x//' **** ERROR: ',a// * ' Columns ',I2,':',I3,'-->',a//' in line:'// * ' Columns 1: 60-->',a/ * ' Columns 61:108-->',a// * 18x,6(10H....,....I)//) stop c 110 write(*,111)'Cannot read intensity from:', * 36,44,line(36:44),line(1:60),line(61:108) stop c 112 write(*,111)'Cannot read lower state energy from line:', * 49,56,line(49:56),line(1:60),line(61:108) stop c 101 write(linasr(25:71),2)(nq(iredir(i)),i=1,12) c * nq1l,nq2l,nq3l 2 format(2(i5,',',i3,',',i3),' ',i2,',',i2,',',i2, * '--',i2,',',i2,',',i2) c kdelm=nq(2)-nq(8) c kmu-kml kdelp=nq(3)-nq(9) c kpu-kpl im=iabs(mod(kdelm,2)) ip=iabs(mod(kdelp,2)) if(im.eq.0.and.ip.eq.1)cdip=' a' if(im.eq.1.and.ip.eq.1)cdip=' b' if(im.eq.1.and.ip.eq.0)cdip=' c' linasr(51:52)=cdip c read(line(49:56),'(f8.3)',err=112)elower linasr(72:79)=line(49:56) c line=linasr c return end c C_____________________________________________________________________________ C subroutine check(line,i,j) c c...routine to decode quantum numbers: c -- two digit negative coded a1 for -11, b1 for -21 etc., up to i1 c for -91 c -- three digit positive coded A0 for 100, B0 for 110 etc., up to N0 c for 240 c into standard I3 notation c character line*108,cdig c read(line(i:j),'(i3)',err=1)k return 1 read(line(i+1:i+1),'(a)')cdig c if( .not.(cdig.ge.'a'.and.cdig.le.'i') .and. * .not.(cdig.ge.'A'.and.cdig.le.'N') )then write(*,111)'Quantum number checks failed on:', * i,j,line(i:j),line(1:60),line(61:108) 111 format(1x//' **** ERROR: ',a// * ' Columns ',I2,':',I3,'-->',a//' in line:'// * ' Columns 1: 60-->',a/ * ' Columns 61:108-->',a// * 18x,6(10H....,....I)//) stop endif c read(line(j:j),'(i1)')k if(cdig.ge.'a'.and.cdig.le.'i')then ntens=ichar(cdig)-ichar('a')+1 k=-(k+10*ntens) endif if(cdig.ge.'A'.and.cdig.le.'N')then ntens=ichar(cdig)-ichar('A')+1 k=(k+(9+ntens)*10) endif write(line(i:j),'(i3)')k c return end C_____________________________________________________________________________ c SUBROUTINE SORTH(NSTART,N) c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c parameter (maxlin=1000000) COMMON /FReq/WK COMMON /point/IPT INTEGER IPT(maxlin) REAL*8 WK(maxlin),WWK C L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________