C----------------------------------------------------------------------------- c c P M I X - Extraction of data from an EGY output file for a plot of c mixing coefficients C C----------------------------------------------------------------------------- C C C Version 30a.X.2007 ----- Zbigniew KISIEL ----- C 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 C Modification history: C C 3.06.06: creation C 18.04.07: mods C 30.10.07: adaptation for change in EGY format C C=========================================================================== c c c CM(i,j,k) - matrix storing mixing coefficients where c i = j+1 c j = tau+j+1 c k = v+1, v being the SPFIT vibrational identifier c c parameter (maxj=300, maxtau=2*maxj+1,maxv=4,maxcut=7) c real*4 cm(maxtau,maxtau,maxv),cutlev(maxcut) integer*4 jmax(maxv),nlev(maxv) character line*100,filin*50,ivc c c data cutlev/0.01, 0.02, 0.05, 0.1, 0.2, 0.4, 1.0/ c WRITE(*,3344) 3344 FORMAT(1X//' ',76('_')/' |',T79,'|'/ * ' | PMIX - Extraction mixing coefficient data from EGY ', * 'output file',T79,'|'/ * ' | for use in a GLE diagnostic plot',T79,'|'/ * ' |',76('_'),'|'/' version 30a.X.2007',T64,'Zbigniew KISIEL'/) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c Read the cutoff levels for 1-Pmix from file PMIX.INP c OPEN (1,file='pmix.inp',status='old') 50 read(1,'(a)')line if(line(1:1).eq.'!')then goto 50 else backspace(1) endif read(1,*)(cutlev(i),i=1,maxcut) close(1) c write(*,51)(cutlev(i),i=1,maxcut) 51 format(1x//' Cutoff levels:'//5(10F7.4)) c c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...Read the EGY file, assuming that leading quantum numbers are three c asymmetric rotor qns followed by the vibrational identifier c 2 write(*,1)' Name of the EGY file:' 1 format(1x//1x,a,' ',$) read(*,'(a)',err=2)filin open(3,file=filin,err=2,status='old') c c older: c c364 122 1682.981020 0.276285 0.644081 362: 90 62 28 2 c c 2007: c c 1 2 21.008525 0.006859 0.863177 4: 0 0 0 1 c 5 read(3,'(a)',end=6,err=5)line if(line(62:62).eq.':')read(line(63:74),'(4i3)',err=5)j,ka,kc,iv if(line(64:64).eq.':')read(line(65:76),'(4i3)',err=5)j,ka,kc,iv c if(j.gt.maxj)then write(*,3)j,maxj 3 format(1x//'**** ERROR: data contains j=',i4/ * ' current limit is j=',i4//) stop endif c if(iv+1.gt.maxv)then write(*,4)iv,maxv-1 4 format(1x//'**** ERROR: data contains v=',i4/ * ' current limit is v=',i4//) stop endif c if(line(62:62).eq.':')read(line(46:56),'(f11.6)',err=5)c if(line(64:64).eq.':')read(line(48:58),'(f11.6)',err=5)c itau=j+1+ka-kc cm(j+1,itau,iv+1)=c nlev(iv+1)=nlev(iv+1)+1 if(j.gt.jmax(iv+1))jmax(iv+1)=j if(iv.gt.ivmax)ivmax=iv c goto 5 c 6 write(*,'(1x)') write(*,7)((i-1,nlev(i),jmax(i)),i=1,ivmax+1) 7 format(' state v=',i2,': no of levels=',i6,' maxJ=',i3) close(3) c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...Generate data files for the plot c do 11 n=1,ivmax+1 c write(ivc,'(i1)')n-1 open(11,file='pm_1_v'//ivc//'.out',status='unknown') open(12,file='pm_2_v'//ivc//'.out',status='unknown') open(13,file='pm_3_v'//ivc//'.out',status='unknown') open(14,file='pm_4_v'//ivc//'.out',status='unknown') open(15,file='pm_5_v'//ivc//'.out',status='unknown') open(16,file='pm_6_v'//ivc//'.out',status='unknown') open(17,file='pm_7_v'//ivc//'.out',status='unknown') write(11,15) write(12,15) write(13,15) write(14,15) write(15,15) write(16,15) write(17,15) 15 format('!'/'! J Ka Kc 1-Pmix Pmix'/'!'/ * ' * * * * *'/'!') do 10 nn=1,jmax(n)+1 j=nn-1 do 10 nnn=1,2*j+1 dcm=1.0-cm(nn,nnn,n) itau=nnn-j-1 kc=(j-itau+1)/2 ka=itau+kc c if(dcm.le.cutlev(1))then write(11,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(2))then write(12,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(3))then write(13,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(4))then write(14,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(5))then write(15,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(6))then write(16,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c if(dcm.le.cutlev(7))then write(17,30)j,ka,kc,dcm,cm(nn,nnn,n) goto 10 endif c 10 continue 30 format(3i5,2f15.6) c close(11) close(12) close(13) close(14) close(15) close(16) close(17) c 11 continue c 9 stop end C_____________________________________________________________________________ C_____________________________________________________________________________