c$DEBUG C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C SYMF - SYMMETRIC TOP FITTING PROGRAM with fully diagonalising treatment C of symmetric top splitting contributions: d2,h2,l2.. or h3,l3.. C depending on top symmetry. C C This program is adapted from 'S'-reduction part of ASFITL and fits C prolate type constants. C c Symmetric top frequencies are normally treatable with closed c formulae except for K-doubling near the spherical limit where c full-diagonalisation may be necessary. This program was created to c deal with both the standard and the near-spherical case, as described c in connection with h3 splitting in the near-spherical molecule c trichloroacetonitrile: c c G.Cazzoli and Z.Kisiel, J.Mol.Spectrosc. 159, 96-102 (1993) C c (note that the near spherical oblate top has not been tested) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C - J up to JMAX=250 C - up to NLINES=600 TRANSITIONS C - WEIGHTED FIT C - OPTIONAL CORRELATION COEFFICIENTS AND CONTRIBUTIONS OF C INDIVIDUAL CONSTANTS TO TRANSITION FREQUENCIES (EVALUATED C FROM DERIVATIVES) C - FIT WITH ZERO DEGREES OF FREEDOM (ie. N constants to N C frequencies) IS ALLOWED, BUT THE USER IS WARNED BY RINGING C THE BELL C C C Ver 6.IV.2001 ------- 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 Modification history: C c early 90's: creation for dealing with h3 splitting in C3v tops C 27.02.01: various enhancements to bring up to standard in ASFIT c 6.04.01: interspersed commenting from QSTARK C C_____________________________________________________________________________ C C This program is standard FORTRAN 77 and compiles with standard C 32-bit Fortran compilers. If a compiler has static allocation of C variables as a compilation option then this needs to be used. C Eg. for f77 compilers use the '-static' compilation option. C_____________________________________________________________________________ C C DATA INPUT: C C The data file is best prepared by reediting an old data file. It should C contain: C C Comment C --------------- C number of transitions: NLINES C number of constants: NCON C --------------- C _____ C | C | DEFINITIONS of constants from B.z up to the NCON'th constant C |_____ C C --------------- C Column headers C --------------- C _____ C | C | DEFINITIONS of measured transitions C |_____ C C C In the above change only NLINES, NCON, and the DEFINITIONS. The various c dividing and comment lines should not be altered. C C It is possible to specify an excessively large value of NLINES and let C the program count the actual number of lines on input. C C To NCON variable defines the number of leading constants in the C Hamiltonian that will be specified (it includes the empty line after B). C C To set further constants than have been written to the data set under c current value of NCON it is best to use the CHANGE CONSTANTS option c of the program. C C Measured transitions are specified by quantum numbers J,K, the measured C frequency, measurement error and a flag on whether or not the line is to C be included in the fit. Several variants are possible as shown in a C block of lines from data set for Cl3CCN below: C C 112 3 375918.732000 0.050000 1 (2I7,F20.6,F15.6,I5) C 112 0- 2 375916.973000 0.050000 1 C 112 4 375915.670000 0.050000 1 C 112 5 375914.220000 0.050000 1 C 112 -3 0.000000 0.050000 0 C 112 D 6 375912.777000 0.050000 1 C C The lines are for J=112, the second column specifies K. This should C be positive. If some form of K-splitting is possible this will C correspond to the higher frequency component and negative value should C be specified for the lower frequency component. An entry such as 'D 6' C specifies that a mean of the two splitting components should be C calculated. The entry '0- 2' specifies that mean frequency of components C K=0,1,2 should be calculated - this can be used for fitting lines C consisting of several unresolved components (these are unit C weighted) C c * Empty lines can be placed between transition lines, and will c be translated as spaces in the same places in the printout c * Lines between transition lines beginning with c the '!' or '%' character are treated as comment lines - but only c those beginning with '!' will be echoed in the output C C The program is not completely converted for symmetric top operation - C about half of the modification options are not available and are C marked N/A. C_____________________________________________________________________________ C C The available constants are as follows: C C B.z, B, dummy 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 and it is up to the user not to set any of the inapplicable constants C C I - number of transitions C N - number of energy levels to be used in fit (if errors are present C in line declarations this may be less than 2*I, or if mean C frequencies of line overlaps are evaluated then this is greater C than 2*I) C K1 - number of fitted constants C FR - frequencies of transitions where C (i,1) observed frequencies C (i,2) calculated frequencies C (i,3) obs-calc C ERROR - estimated experimental error in frequency C WEIGHT - weight equal to 1/error**2 C NJA - information on energy levels where C (i,1) J quantum number C (i,2) number of Wang submatrix to which the level belongs =1,2,3,4 for C E+, E-, O+, O- resp. C (i,3) position of the energy level among eigenvalues of the relevant C Wang subm. in ascending order of magnitude C (i,4) number of transition associated with the level, this is negative C for the lower level C IA - pointers to which constants are to be fitted (=1, otherwise 0) C LINFIT - pointers to which lines are to be fitted (+1, otherwise 0) C JK - quantum numbers for transitions where C (i,1), (i,2), (i,3) are J, K-1, K+1 for upper level resp. C (i,4), (i,5), (i,6) are J, K-1, K+1 for lower level resp. C IPR - default number of iterations c c COMBLK - string variable holding all comments concerning transitions c ICOMML(i) - number of transition to be preceded by comment i c ICOMM(i) - position of last character in COMBLK for comment i c ISPACE(i) - number of transition to be preceded by an empty line C_____________________________________________________________________________ C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=600, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1) PARAMETER (maxcbl=2000,maxcom=100,maxspa=100) c COMMON /BLKH/H/BLKR/R COMMON /BLK1/JK,FR,ERROR,WEIGHT,NDEG,LINFIT,A,IA,I, 1 IOUT/BLK2/NJA,B/BLK3/C,D,DD/TBLOCK/COM,T DIMENSION A(NCONST),B(NCONST),C(NCONST),Z(NCONST),IA(NCONST), 1 H(MAXDIM,MAXDIM),R(MAXDIM,MAXDIM),DD(NCONST,NCONST) INTEGER*2 NJA(NNJA,4) DIMENSION FR(NLINES,3),JK(NLINES,6),LINFIT(NLINES),NDEG(NLINES,3) REAL*4 D(NLINES,NCONST),EE,ERROR(NLINES),WEIGHT(NLINES) CHARACTER*6 T(NCONST),TT(NCONST) CHARACTER COM(72),KQUANT*7 CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk c DATA TT/'B.z =','B =',' ',' DJ =',' DJK =',' DK =', 1 ' d1 =',' d2 =','HJ =','HJK =','HKJ =','HK =', 1 'h1 =','h2 =','h3 =', 1 ' LJ =',' LJJK=',' LJK =',' LKKJ=',' LK =',' l1 =', 1 ' l2 =',' l3 =',' l4 =', 1 'PJ =','PJJK =','PJK =','PKJ =','PKKJ =','PK =', 1 'p1 =','p2 =','p3 =','p4 =','p5 ='/ DO 1234 L=1,NCONST T(L)=TT(L) 1234 CONTINUE IOUT=0 C IFLAG=0 CALL DATAIN 400 CALL MODIF(N,IFLAG) IF(IOUT.EQ.1)GOTO 41 C I1=0 S=1.D50 IPR=3 C C C START ITERATION C C 14 WRITE(*,6000) 6000 FORMAT(/' **** CALCULATION IN PROGRESS, CURRENTLY WORKING ON J', 1 ' ='/) I1=I1+1 S1=S DO 15 K=1,I FR(K,2)=0.D0 FR(K,3)=0.D0 DO 15 L=1,NCONST D(K,L)=0.0 15 CONTINUE L1=1 M4=0 WRITE(*,801)NJA(1,1) C C Go through all sorted energy levels C DO 25 K=2,N KK=K-1 IF(NJA(K,1).EQ.NJA(KK,1))GOTO 250 WRITE(*,801)NJA(KK,1) 801 FORMAT(I4,$) 250 IF(NJA(K,1).EQ.NJA(KK,1).AND.NJA(K,2).EQ.NJA(KK,2))GOTO 25 L2=K-1 N1=NJA(KK,2) J1=NJA(KK,1) C C Set up and diagonalize the appropriate Wang matrix for each group of C common levels C CALL ASREDS(M,J1,N1,A) IEGEN=0 CALL HDIAG(M,IEGEN,K1) DO 17 L=L1,L2 M=NJA(L,3) M1=NJA(L,4) M2=ABS(M1) FR(M2,2)=FR(M2,2)+SIGN(1,M1)*H(M,M)/NDEG(M2,1) 17 CONTINUE C C FIND THE DERIVATIVES - since transformation to the space where the Wang C submatrix is diagonal has already been found (eigenvectors), the deri- C vatives for energy levels in this submatrix with respect to a given C constant are obtained by the same similarity transformation of matrix C set up with all other constants set to zero and the required constant C set to 1 (array B) C K1=0 DO 24 J=1,NCONST IF(IA(J).EQ.1)GOTO 18 GOTO 24 18 K1=K1+1 DO 19 L=1,NCONST B(L)=0.D0 19 CONTINUE B(J)=1.D0 CALL ASREDS(M3,J1,N1,B) DO 124 L=L1,L2 M=NJA(L,3) IF(L.NE.L1.AND.M.EQ.M4)GOTO 23 EE=0.0 DO 22 M1=1,M3 DO 122 M2=M1,M3 N2=2 IF(M1.EQ.M2)N2=1 EE=EE+N2*H(M1,M2)*R(M1,M)*R(M2,M) 122 CONTINUE 22 CONTINUE 23 M1=NJA(L,4) M2=ABS(M1) D(M2,K1)=D(M2,K1)+SIGN(1,M1)*EE/NDEG(M2,1) M4=M 124 CONTINUE 24 CONTINUE L1=K 25 CONTINUE C C FIND NEW CONSTANTS - inverse of D'D is evaluated by diagonalization using C inv A = U * inv L * U' where U and L contain the eigenvectors and the C eigenvalues of A. C C To improve conditioning D'D is rescaled to D'D.ij/(D'D.ii D'Djj)**1/2 C (scaling factors in Z) C To save space D'D is placed and diagonalized in the top lh. corner C of H, D'Y is placed in (K1+1)'st column of H and inv D'D is placed C in array DD equivalenced to the rh. side of H C S=0.D0 SW=0.D0 L=0 DO 26 J=1,I FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) IF(LINFIT(J).NE.1)GOTO 26 S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERROR(J))**2 L=L+1 26 CONTINUE N139=L-K1 IF(N139.LT.0)GOTO 48 IF(N139.GT.0)GOTO 810 WRITE(*,89)L,K1,CHAR(7) SDEV=DSQRT(S) SWDEV=DSQRT(SW) GOTO 811 810 SDEV=DSQRT(S/N139) SWDEV=DSQRT(SW/N139) 811 WRITE(*,78)I1,SDEV,SWDEV 78 FORMAT(//' ITERATION STEP NO.',I2,11X,'Deviation ', 1 'of fit =',F14.6/24X,'Weighted deviation of fit =',F14.6) DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(LINFIT(L).NE.1)GOTO 27 E=E+D(L,J)*D(L,K)*WEIGHT(L) 27 CONTINUE H(J,K)=E 28 CONTINUE Z(J)=SQRT(1.D0/H(J,J)) E=0.D0 DO 29 L=1,I IF(LINFIT(L).NE.1)GOTO 29 E=E+D(L,J)*FR(L,3)*WEIGHT(L) 29 CONTINUE H(J,K1+1)=E 30 CONTINUE DO 301 J=1,K1 DO 301 K=J,K1 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 301 CONTINUE IEGEN=0 CALL HDIAG(K1,IEGEN,L3) DO 35 J=1,K1 B(J)=0.D0 DO 34 K=1,K1 E=0.D0 DO 33 L=1,K1 E=E+R(J,L)*R(K,L)/H(L,L) 33 CONTINUE E=E*Z(J)*Z(K) DD(J,K)=E B(J)=B(J)+E*H(K,K1+1) 34 CONTINUE C(J)=SQRT(DD(J,J))*SWDEV 35 CONTINUE C WRITE(*,81) 81 FORMAT(/' RESULTS - constant, change, standard deviation :'/) K=0 DO 40 J=1,NCONST IF(IA(J).EQ.1)GOTO 39 GOTO 40 39 K=K+1 A(J)=A(J)+B(K) IF(J.LE.8)THEN WRITE(*,82)T(J),A(J),B(K),C(K) ELSE WRITE(*,1082)T(J),A(J),B(K),C(K) ENDIF 82 FORMAT(1X,A6,3F23.12) 1082 FORMAT(1X,A6,F27.16,2F23.16) 40 CONTINUE C C Convergence tests and output C 403 WRITE(*,'(/'' Do you want to interrupt the fit (1/0) ? '',$)') READ(*,5,ERR=403)IIFLAG IF(IIFLAG.NE.1)GOTO 6 C 404 WRITE(*,401) 401 FORMAT(1X/ * 9x,'= 0 OBS-CALC DIFFERENCES = 2 DISK OUTPUT'/ * 9x,'=-1 as above, screen at a time = 3 NEXT ITERATION'/ * 9x,'= 1 MODIFICATIONS = 4 EXIT SYMF', * 9X,'... ',$) READ(*,5,ERR=404)IFLAG 5 FORMAT(I5) IF(IFLAG.LT.-1.OR.IFLAG.GT.4)GOTO 404 IF(IFLAG.EQ.1)GOTO 400 IF(IFLAG.EQ.3)GOTO 6 IF(IFLAG.EQ.4)GOTO 99 C c...obs-calc differences c IF(IFLAG.LE.0)THEN WRITE(*,'(1X)') do 108 jj=1,i,24 istart=jj ifin=jj+23 if(ifin.gt.i)ifin=i DO 7 NN=istart,ifin J=JK(NN,4) call label(NN,kquant) IF(LINFIT(NN).NE.1)GOTO 70 WRITE(*,8)NN,J,KQUANT,FR(NN,1),FR(NN,2),FR(NN,3),ERROR(NN) GOTO 7 70 WRITE(*,71)NN,J,KQUANT,FR(NN,1),FR(NN,2),FR(NN,3) 7 CONTINUE if(iflag.eq.-1)then write(*,'('' Press ENTER to continue '',$)') read(*,'(i1)',err=108)j endif 108 continue WRITE(*,78)I1,SDEV,SWDEV GOTO 404 ENDIF 8 FORMAT(1X,I3,'. ',I7,A7,F18.4,2F14.4,F11.4) 71 FORMAT(1X,I3,'. ',I7,A7,F18.4,2F14.4,'--excluded') c c...disk output c IF(IFLAG.EQ.2)THEN IOUT=2 CALL DATOUT(iiflag,SDEV,SWDEV) GOTO 404 ENDIF C c...exit to next iteration c 6 IFLAG=1 E=(S1-S)/S IF(DABS(E).LT.0.001D0)GOTO 41 IF(E.LT.-0.1D0)GOTO 41 IF(I1.LT.IPR)GOTO 14 1900 WRITE(*,450) 450 FORMAT(/' MORE CYCLES ? ',$) READ(*,'(I5)',ERR=1900)L IF(L.NE.1)GOTO 41 I1=0 GOTO 14 C 41 CALL DATOUT(IFLAG,SDEV,SWDEV) C IF(IFLAG.EQ.1)GOTO 400 GOTO 99 48 WRITE(*,89)L,K1,' - ITERATION STOPPED' 89 FORMAT(//' ONLY',I3,' MEASURED TRANSITIONS FOR',I3, 1 ' CONSTANTS',A) GOTO 41 C 99 STOP END C C_____________________________________________________________________________ C SUBROUTINE DATAIN C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=600, NCONST=35, JMAX=250) PARAMETER (maxcbl=2000,maxcom=100,maxspa=100) c COMMON /BLK1/JK,FR,ERROR,WEIGHT,NDEG,LINFIT,A,IA,I,IOUT COMMON /TBLOCK/COM,T DIMENSION JK(NLINES,6),FR(NLINES,3),A(NCONST),IA(NCONST), * LINFIT(NLINES),NDEG(NLINES,3) CHARACTER*6 T(NCONST) CHARACTER*1 COM(72) CHARACTER*25 FILNAM,cdummy*22 REAL*4 ERROR(NLINES),WEIGHT(NLINES) CHARACTER comblk*(maxcbl),line*80 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk nspace=0 ncomm=1 icomm(ncomm)=0 C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | SYMF - FITTING OF SYMMETRIC TOP HAMILTONIAN WITH SPLIT', * 'TING TERMS',T79,'|'/ * ' | (based on prolate type reduction S)', * T79,'|'/ * ' |',76(1H_),'|'/' version 6.IV.2001',T64,'Zbigniew KISIEL'/) write(*,1910)nlines,jmax 1910 format(' Current limits:',i7,' lines and J up to',i4/) 14 FORMAT(I5) 141 format(a22,i5) DO 200 L=1,NCONST A(L)=0.0D0 IA(L)=0.0D0 200 CONTINUE C C...input (only from a data file on disk) C 15 WRITE(*,135) 135 FORMAT(/' NAME OF DATA FILE : ',$) READ(*,'(A)',ERR=15)FILNAM OPEN(2,FILE=FILNAM,ERR=15,status='OLD') READ(2,3500,ERR=15)COM read(2,3500,err=15)cdummy READ(2,141,ERR=15)cdummy,I read(2,141,err=15)cdummy,NCON read(2,3500,err=15)cdummy READ(2,150,ERR=15)(IA(K),A(K),K=1,NCON) read(2,3500,err=15)cdummy read(2,3500,err=15)cdummy read(2,3500,err=15)cdummy c c...transitions c nlin=0 nerr=0 DO 3400 KK=1,I KKK=nlin+1 c 23 read(2,'(a)',err=351,end=51)line if(line(1:1).eq.'!'.or.line(1:1).eq.'%')then if(line(1:1).eq.'%')goto 23 if(ncomm.ge.maxcom)then write(*,35)'Too many comment lines',char(7) goto 23 endif if(icomm(ncomm)+lentrm(line).gt.maxcbl)then write(*,35)'Capacity of the comments buffer exceeded', * char(7) goto 23 endif 35 format(' **** ERROR: ',2a) ncomm=ncomm+1 icomml(ncomm)=kk icomm(ncomm)=icomm(ncomm-1)+lentrm(line) comblk(icomm(ncomm-1)+1:icomm(ncomm))=line(1:lentrm(line)) goto 23 endif if(line(1:3).eq.'end')goto 51 if(lentrm(line).lt.20)then nspace=nspace+1 if(nspace.lt.maxspa)ispace(nspace)=kk if(nspace.gt.1.and.i.eq.ispace(nspace-1))nspace=nspace-1 goto 23 endif backspace(2) c READ(2,'(A)',err=3401,end=51)LINE IF(LINE(12:12).EQ.'-')THEN READ(LINE,3406,err=3401,end=51)J,NDEG(KKK,2),NDEG(KKK,3), * FR(KKK,1),ERROR(KKK),LINFIT(KKK) NDEG(KKK,1)=NDEG(KKK,3)-NDEG(KKK,2)+1 K=NDEG(KKK,2) 3406 FORMAT(I7,I4,1X,I2,F20.6,F15.6,I5) ELSE IF(LINE(11:11).EQ.'D')THEN READ(LINE,3408,err=3401,end=51)J,K,FR(KKK,1),ERROR(KKK), * LINFIT(KKK) 3408 FORMAT(I7,5X,I2,F20.6,F15.6,I5) NDEG(KKK,1)=2 NDEG(KKK,2)=K NDEG(KKK,3)=K ELSE NDEG(KKK,1)=1 READ(LINE,3407,err=3401,end=51)J,K,FR(KKK,1),ERROR(KKK), * LINFIT(KKK) 3407 FORMAT(2I7,F20.6,F15.6,I5) ENDIF JK(KKK,1)=J+1 JK(KKK,4)=J IF(K.GE.0)THEN JK(KKK,3)=JK(KKK,1)-K JK(KKK,6)=JK(KKK,4)-K ELSE K=IABS(K) JK(KKK,3)=JK(KKK,1)-K+1 JK(KKK,6)=JK(KKK,4)-K+1 ENDIF JK(KKK,2)=K JK(KKK,5)=K nlin=nlin+1 goto 3400 3401 IF(LINE(1:10).EQ.'----------')GOTO 3400 351 nerr=nerr+1 if(nerr.eq.1)write(*,'(1x)') write(*,'('' ----> transition error:''/i5,'': '',a)') * kk,LINE(1:72) goto 23 3400 CONTINUE c 51 I=NLIN if(nerr.gt.0)then write(*,'(1x/i5, * '' errors encountered, press ENTER '',$)')nerr read(*,'(i1)',err=3402)nerr endif c 3402 CLOSE(2) 3500 FORMAT(1X,72A1) 150 FORMAT(I5,F30.19) 140 FORMAT(6I5,F20.6,F15.6,I5) DO 1020 K=1,I IF(ERROR(K).EQ.0.)ERROR(K)=1.0 WEIGHT(K)=1./ERROR(K)**2 1020 CONTINUE C 41 RETURN END C C_____________________________________________________________________________ C subroutine label(NN,KQUANT) c IMPLICIT INTEGER*2 (I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NLINES=600, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1) c character KQUANT*7 COMMON /BLK1/JK,FR,ERROR,WEIGHT,NDEG,LINFIT,A,IA,I,IOUT DIMENSION JK(NLINES,6),FR(NLINES,3),A(NCONST),IA(NCONST), 1 LINFIT(NLINES),NDEG(NLINES,3) REAL*4 ERROR(NLINES),WEIGHT(NLINES) c c Place in string KQUANT the appropriate designator for type of K value c associated with transition NN c a/ simple c b/ mean of two near degenerate values c c/ mean of values in a specified range c IF(NDEG(NN,1).EQ.1)THEN KKK=JK(NN,4)-(JK(NN,5)+JK(NN,6)) K=ISIGN(1,KKK)*JK(NN,5) WRITE(KQUANT,3420)K 3420 FORMAT(I7) ELSE IF(NDEG(NN,2).EQ.NDEG(NN,3))THEN WRITE(KQUANT,3461)NDEG(NN,2) 3461 FORMAT(' D ',I2) ELSE WRITE(KQUANT,3421)NDEG(NN,2),NDEG(NN,3) 3421 FORMAT(I4,'-',I2) ENDIF c return end C_____________________________________________________________________________ C SUBROUTINE MODIF(N,IFLAG) C C Modification of data for the fit and identification and sorting of C energy levels: C C N - on output contains number of energy levels identified for fitting C IFLAG - if set to 1 then check of data, listing and reevaluation of NJA C are skipped C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=600, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1) c COMMON /BLK1/JK,FR,ERROR,WEIGHT,NDEG,LINFIT,A,IA,I,IOUT COMMON /BLK2/NJA,B/TBLOCK/COM,T DIMENSION JK(NLINES,6),FR(NLINES,3),A(NCONST),IA(NCONST), 1 B(NCONST),LINFIT(NLINES),NDEG(NLINES,3) CHARACTER*6 T(NCONST) INTEGER*2 NJA(NNJA,4),IINT(6),JJK(6) CHARACTER COM(72),kquant*7 REAL*4 ERROR(NLINES),WEIGHT(NLINES) NQNS=6 C 30 FORMAT(/1X,72A1) 31 FORMAT(/' CURRENT CONSTANTS :'/) 33 FORMAT(1X,A6,F22.12,10X,I1) 34 FORMAT(/' MEASURED TRANSITIONS : ',I3/) 96 FORMAT(1X,I3,'. ',2(I5,2I3),F19.4,F14.4) 97 FORMAT(1X,I3,'. ',2(I5,2I3),F19.4,F14.4,'***INPUT ERROR') 98 FORMAT(1X,I3,'. ',2(I5,2I3),F19.4,F14.4,' EXCL.') C IF(IFLAG.EQ.1)GOTO 11 C C Check of data, evaluation of NJA (identification of energy levels) and C listing C 29 WRITE(*,30)COM WRITE(*,31) DO 32 K=1,NCONST IF(IA(K).EQ.0.AND.A(K).EQ.0.D0)GOTO 32 WRITE(*,33)T(K),A(K),IA(K) 32 CONTINUE WRITE(*,34)I N=0 DO 10 K=1,I FR(K,1)=ABS(FR(K,1)) B(1)=JK(K,1)-JK(K,4) C C Test whether J''-J' = 0, -1 or +1 and J less than allowed maximum C IF(ABS(B(1)).GT.1.D0.OR.JK(K,1).GT.JMAX)GOTO 9 DO 6 L=2,5,3 J=JK(K,L)+JK(K,L+1)-JK(K,L-1) C C Test whether (K-1 + K+1) - J = 0 or 1 C IF(J.LT.0.OR.J.GT.1)GOTO 7 N=N+1 J=J+(JK(K,L)-JK(K,L)/2*2)*2+1 M=J/2*2 NJA(N,1)=JK(K,L-1) NJA(N,2)=J NJA(N,3)=(JK(K,L)-(J-1)/2-(4-J)/2*M)/2+1 NJA(N,4)=K IF(L.EQ.5)NJA(N,4)=-K 6 CONTINUE c c...list the line c jjJ=JK(k,4) call label(k,kquant) IF(LINFIT(k).NE.1)GOTO 1270 WRITE(*,1208)k,jjJ,KQUANT,FR(k,1),ERROR(k) GOTO 10 1270 WRITE(*,1271)k,jjJ,KQUANT,FR(k,1),ERROR(k) GOTO 10 1208 FORMAT(1X,I3,'. ',I7,A7,F18.4,F14.4) 1271 FORMAT(1X,I3,'. ',I7,A7,F18.4,F14.4'--excluded') c c IF(LINFIT(K).NE.1)GOTO 99 c WRITE(*,96)K,(JK(K,M),M=1,NQNS),FR(K,1),ERROR(K) c GOTO 10 c99 WRITE(*,98)K,(JK(K,M),M=1,NQNS),FR(K,1),ERROR(K) c GOTO 10 C C If error in lower level, eliminate also the upper level from NJA C 7 IF(L.EQ.5)GOTO 8 GOTO 9 8 N=N-1 9 WRITE(*,97)K,(JK(K,M),M=1,NQNS),FR(K,1),ERROR(K) FR(K,1)=0.D0 10 CONTINUE C C...set NJA also for other transitions required for frequencies of overlaps C DO 710 K=1,I IF(NDEG(K,1).EQ.1)GOTO 710 JJK(1)=JK(K,1) JJK(4)=JK(K,4) IF(NDEG(K,2).EQ.NDEG(K,3))THEN KKSTRT=NDEG(K,2) ELSE KKSTRT=NDEG(K,2)+1 ENDIF DO 711 KK=KKSTRT,NDEG(K,3) JJK(3)=JJK(1)-KK JJK(6)=JJK(4)-KK JJK(2)=KK JJK(5)=KK IF(NDEG(K,2).EQ.NDEG(K,3))THEN JJK(3)=JJK(3)+1 JJK(6)=JJK(6)+1 ENDIF DO 706 L=2,5,3 J=JJK(L)+JJK(L+1)-JJK(L-1) N=N+1 J=J+(JJK(L)-JJK(L)/2*2)*2+1 M=J/2*2 NJA(N,1)=JJK(L-1) NJA(N,2)=J NJA(N,3)=(JJK(L)-(J-1)/2-(4-J)/2*M)/2+1 NJA(N,4)=K IF(L.EQ.5)NJA(N,4)=-K 706 CONTINUE 711 CONTINUE 710 CONTINUE C C Selection and execution of required modification C 1900 WRITE(*,40) 40 FORMAT(/' MODIFICATIONS ? ',$) READ(*,14,ERR=1900)IFLAG 14 FORMAT(I5) IF(IFLAG.NE.1)GOTO 41 11 WRITE(*,70) 70 FORMAT(/' MODIFICATIONS ARE SELECTED AS FOLLOWS :' 1 //' = 0 LISTING AND EXIT TO CALCULATION',6X, 1 'N/A INSERT LINES'/' = 1 CHANGE CONSTANTS',21X, 1 '= 6 CHANGE THE COMMENT'/' N/A DELETE LINES',25X, 1 '= 7 STORE DATA'/' N/A CHANGE LINES',25X, 1 '= 8 EXCLUDE LINES FROM FIT'/' N/A ADD LINES',28X, 1 '= 9 INCLUDE LINES IN FIT'/43X,'= 10 CHANGE ERROR'/) 72 WRITE(*,71) 71 FORMAT(/' CONTROL PARAMETER = ',$) READ(*,14,ERR=72)ICONTR ICONTR=ICONTR+1 GOTO(29,74, 11, 11, 11, 11,178,179,210,220,300),ICONTR c GOTO(29,74, 75, 76, 77, 500,178,179,210,220,300),ICONTR C C Change of constant C 74 WRITE(*,31) LL=NCONST/2 DO 37 L=1,LL WRITE(*,35)(K,T(K),A(K),IA(K),K=L,L+LL,LL) 37 CONTINUE 35 FORMAT(1X,I2,': ',A6,F26.16,I2,2X,I2,': ',A6,F26.16,I2) IF(2*LL.LT.NCONST)WRITE(*,36)NCONST,T(NCONST),A(NCONST), 1 IA(NCONST) 36 FORMAT(39X,2X,I2,': ',A6,F26.16,I2) 1910 WRITE(*,78) 78 FORMAT(/' NUMBER OF CONSTANT TO BE CHANGED (ENTER=exit) ? ',$) READ(*,14,ERR=1910)NN IF(NN.LT.1.OR.NN.GT.NCONST)GOTO 72 1911 WRITE(*,79)T(NN) 79 FORMAT(/' FITTING PARAMETER, NEW VALUE OF ',A6,' ',$) READ(*,*,ERR=1911)IA(NN),A(NN) GOTO 74 C C Line deletion C 75 WRITE(*,81) 81 FORMAT(/' NUMBER OF LINE TO BE DELETED ? ',$) READ(*,14,ERR=75)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 I=I-1 DO 1 L=NN,I FR(L,1)=FR(L+1,1) LINFIT(L)=LINFIT(L+1) ERROR(L)=ERROR(L+1) WEIGHT(L)=WEIGHT(L+1) DO 1 LL=1,NQNS JK(L,LL)=JK(L+1,LL) 1 CONTINUE GOTO 75 C C Change of lines C 76 WRITE(*,82) 82 FORMAT(/' NUMBER OF LINE TO BE CHANGED ? ',$) READ(*,14,ERR=76)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 1912 WRITE(*,90)NN,(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) 90 FORMAT(/' LINE NO.',I3,' :',2(I5,2I3),F18.4,F14.4/ 1 ' NEW LINE : ',$) READ(*,*,ERR=1912)(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) IF(ERROR(NN).EQ.0.)ERROR(NN)=1. WEIGHT(NN)=1./ERROR(NN)**2 GOTO 76 C C Line addition C 77 I=I+1 WRITE(*,190) 190 FORMAT(/' LINE (LINES) TO BE ADDED :'/) 191 READ(*,*,ERR=1915)(JK(I,L),L=1,NQNS),FR(I,1),ERROR(I) GOTO 1916 1915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 191 1916 IF(JK(I,1).EQ.0.AND.JK(I,4).EQ.0)GOTO 192 LINFIT(I)=1 IF(ERROR(I).EQ.0.)ERROR(I)=1. WEIGHT(I)=1./ERROR(I)**2 I=I+1 GOTO 191 192 I=I-1 GOTO 72 C C Line insertion C 500 WRITE(*,501) 501 FORMAT(/' INSERTIONS ARE TO PRECEDE LINE NO. : ',$) READ(*,14,ERR=500)NINS WRITE(*,502) 502 FORMAT(/' LINE (LINES) TO BE INSERTED :'/) 503 READ(*,*,ERR=1918)(IINT(L),L=1,NQNS),FREQ,TEMP GOTO 1919 1918 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 503 1919 IF(IINT(1).EQ.0.AND.IINT(4).EQ.0)GOTO 72 I=I+1 DO 504 L=I,NINS,-1 FR(L,1)=FR(L-1,1) LINFIT(L)=LINFIT(L-1) ERROR(L)=ERROR(L-1) WEIGHT(L)=WEIGHT(L-1) DO 504 LL=1,NQNS JK(L,LL)=JK(L-1,LL) 504 CONTINUE FR(NINS,1)=FREQ LINFIT(NINS)=1 IF(TEMP.EQ.0.)TEMP=1. ERROR(NINS)=TEMP WEIGHT(NINS)=1./TEMP**2 DO 505 LL=1,NQNS JK(NINS,LL)=IINT(LL) 505 CONTINUE NINS=NINS+1 GOTO 503 C C Change of comment C 178 WRITE(*,'(/'' Old comment:''/1x,72a1)')COM WRITE(*,180) 180 FORMAT(/' New comment:') READ(*,181,ERR=178)COM 181 FORMAT(72A1) GOTO 72 C C Saving of data C 179 IOUT=1 GOTO 200 C C Line exclusion C 210 WRITE(*,211) 211 FORMAT(/' NUMBER OF LINE TO BE EXCLUDED FROM FIT ? ',$) READ(*,14,ERR=210)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=0 GOTO 210 C C Line inclusion C 220 WRITE(*,221) 221 FORMAT(/' NUMBER OF LINE TO BE INCLUDED IN FIT ? ',$) READ(*,14,ERR=220)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=1 GOTO 220 C C Change of error C 300 WRITE(*,301) 301 FORMAT(/' ERROR TO BE CHANGED ON LINE NUMBER ? ',$) READ(*,'(I5)',ERR=300)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 call label(nn,kquant) WRITE(*,302)NN,jk(NN,4),kquant,FR(NN,1),ERROR(NN) 302 FORMAT(/' LINE NO.',I4,' :',I6,a,F18.4,' +-',F14.4/ 1 ' NEW ERROR : ',$) READ(*,*,ERR=300)TEMP IF(TEMP.EQ.0.)TEMP=1. ERROR(NN)=TEMP WEIGHT(NN)=1./TEMP**2 GOTO 300 C C Sort of energy levels in NJA according to J, Wang submatrix and order C of eigenvalue in the submatrix C 41 DO 13 J=1,N-1 K1=J+1 L1=J DO 12 K=K1,N IF(NJA(K,1)-NJA(L1,1))411,401,12 401 IF(NJA(K,2)-NJA(L1,2))411,402,12 402 IF(NJA(K,3)-NJA(L1,3))411,12,12 411 L1=K 12 CONTINUE IF(L1.EQ.J)GOTO 13 DO 400 L=1,4 M=NJA(J,L) NJA(J,L)=NJA(L1,L) NJA(L1,L)=M 400 CONTINUE 13 CONTINUE N=N+1 NJA(N,1)=0 NJA(N,2)=0 NJA(N,3)=0 C 200 RETURN END C C_____________________________________________________________________________ C SUBROUTINE DATOUT(IFLAG,SDEV,SWDEV) C C OUTPUT OF DATA AND RESULTS FROM SYMF - BOTH TO SCREEN AND TO DISK C C IFLAG - on output signals whether further modifications are required C SDEV - standard deviation of the fit C NCONST - total no of fittable constants C also: C IOUT - if set to 1 on input, disk output of data is carried out only C if set to 2, disk output of both data and intermediate results C is enabled C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=600, NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) PARAMETER (maxcbl=2000,maxcom=100,maxspa=100) c COMMON /BLKH/H/BLKR/R COMMON /BLK1/JK,FR,ERROR,WEIGHT,NDEG,LINFIT,A,IA,I,IOUT COMMON /BLK3/C,D,DD/TBLOCK/COM,T INTEGER*2 JK(NLINES,6),IA(NCONST),TPT(NCONST),LINFIT(NLINES), 1 NDEG(NLINES,3),UNITS(NCONST) REAL*8 FR(NLINES,3),A(NCONST),C(NCONST),H(MAXDIM,MAXDIM), 1 R(MAXDIM,MAXDIM),DD(NCONST,NCONST) REAL*4 D(NLINES,NCONST),DER,DIV,ERROR(NLINES),WEIGHT(NLINES) CHARACTER*6 T(NCONST) CHARACTER*7 UNIT,NAMUN(6) CHARACTER*25 FILNAM,conval,erval CHARACTER COM(72),FORMAR(54),KQUANT*7 CHARACTER*38 OUTCON(NCONST) CHARACTER*54 FORMAL,F1,F2 EQUIVALENCE (FORMAL,FORMAR(1)) CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk c NOUT=0 C C...UNITS defines the units in which various constants are to be printed C out, NAMUN contains the names of these units DATA NAMUN/'MHz ','kHz ',' Hz ','mHz ','microHz', * 'nHz '/ DATA UNITS/1,1,1, 2,2,2,2,2, 4,4,4,4,4,4,4, 5,5,5,5,5,5,5,5,5, * 6,6,6,6,6,6,6,6,6,6,6/ C C...F1 - format for first line of printout of contributions, F2 - format C for all following lines for the same transition F1='(1X,I3,3H. ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' F2='(1X,6H ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' IF(IOUT.EQ.1)GOTO 100 C K1=0 DO 9 N=1,NCONST IF(IA(N).NE.1)GOTO 9 K1=K1+1 TPT(K1)=N 9 CONTINUE C IF(IOUT.EQ.2)GOTO 300 C C...Output of results (standard listing) which can be routed to screen C (NOUT=0) or to disk (NOUT=3) C 25 WRITE(NOUT,1)COM 1 FORMAT( T54,'SYMF version 6.IV.2001'/ 1 1X,79(1H-)/1X,72A1/1X,79(1H-)//' Fit of symmetric', 1 ' top Hamiltonian with splitting terms:'//) L=0 DO 2 N=1,NCONST IF(A(N).EQ.0.D0)GOTO 2 IF(IA(N).EQ.0)GOTO 3 L=L+1 WRITE(NOUT,4)T(N),A(N),C(L) GOTO 2 3 WRITE(NOUT,5)T(N),A(N) 2 CONTINUE 4 FORMAT(1X,A6,F30.19,' +-',F25.19) 5 FORMAT(1X,A6,F30.19) WRITE(NOUT,405)SDEV,SWDEV 405 FORMAT(/' Deviation of fit = ',F14.6/ 1 ' Weighted deviation of fit = ',F14.6) c c...TRANSITIONS c WRITE(NOUT,'(1x)') WRITE(NOUT,53) WRITE(NOUT,6) 6 FORMAT(10x,' J'''' K',9x,'Observed',7X,'Calc',10x,'Obs-Calc', 1 6x,'Error') WRITE(NOUT,53) NLF=0 nspace=1 ncomm=2 DO 7 N=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 33 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 33 endif c J=JK(N,4) call label(n,kquant) IF(LINFIT(N).NE.1)GOTO 70 NLF=NLF+1 WRITE(NOUT,8)N,J,KQUANT,FR(N,1),FR(N,2),FR(N,3),ERROR(N) GOTO 7 70 WRITE(NOUT,71)N,J,KQUANT,FR(N,1),FR(N,2),FR(N,3) 7 CONTINUE WRITE(NOUT,53) c 8 FORMAT(1X,I3,'. ',I7,A7,F18.4,2F14.4,F11.4) 71 FORMAT(1X,I3,'. ',I7,A7,F18.4,2F14.4,'--excluded') IF(NOUT.NE.0)WRITE(NOUT,410)SDEV,NLF,SWDEV,NLF-K1 410 FORMAT(/' Sigma = ',F14.6,28X,I4,' Transitions in fit'/ 1 ' Sigma.w = ',F14.6,28X,I4,' Degrees of freedom') C C...CONSTANTS c C...On output to disk additional conversion into customary units C and writing in the convention of error given in units of last quoted C digit in constant value (three digit error) C IF(NOUT.EQ.3)THEN WRITE(NOUT,'(1X)') L=0 DO 1002 N=1,NCONST if(N.eq.3)then outcon(n)=' ' goto 1002 endif UNIT=NAMUN(UNITS(N)) ASCLD=A(N)*10.D0**(3*(UNITS(N)-1)) IF(IA(N).EQ.0)GOTO 1003 C C...value with error L=L+1 ESCLD=C(L)*10.D0**(3*(UNITS(N)-1)) WRITE(CONVAL,'(F25.16)')ASCLD WRITE(ERVAL,'(F25.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR) if(ndcon+nderor.gt.22)ndcon=22-nderor WRITE(OUTCON(N),411)T(N)(1:5),UNIT,CONVAL(1:NDCON), * ERVAL(1:NDEROR) GOTO 1002 C C...value without error, formatting is carried out to insert zero before C the decimal point and to discard digits from rounding errors 1003 IF(ASCLD.NE.0.D0)ASCLD=ASCLD+MOD(ASCLD,1.D-15*ASCLD) WRITE(CONVAL,'(F25.16)')ASCLD IF(ASCLD.EQ.0.0D0)CONVAL=' 0.0 ' IF(CONVAL(8:8).EQ.' ')CONVAL(8:8)='0' IF(CONVAL(8:8).EQ.'-')THEN CONVAL(8:8)='0' CONVAL(7:7)='-' ENDIF c iend=25 do 1104 m=1,25 if(conval(m:m).eq.'0'.and.iend.eq.25)then iend=m-1 goto 1104 endif if(conval(m:m).eq.'.')idot=m if((conval(m:m).ne.'.'.and.conval(m:m).ne.'0') * .and.iend.ne.25)then if(m-2-iend.ge.5)goto 1105 iend=25 goto 1104 endif 1104 continue 1105 if(idot.gt.iend)iend=idot DO 1114 M=25,1,-1 IF(CONVAL(M:M).NE.'0'.AND.CONVAL(M:M).NE.' ')GOTO 1115 IF(CONVAL(M:M).EQ.'0')CONVAL(M:M)=' ' 1114 CONTINUE 1115 WRITE(OUTCON(N),1005)T(N)(1:5),UNIT,CONVAL(1:iend) 1002 CONTINUE c LLLL=NCONST/2 DO 1107 N=1,LLLL WRITE(NOUT,1108)OUTCON(N),OUTCON(N+LLLL) 1107 CONTINUE IF(2*LLLL.NE.NCONST)WRITE(NOUT,1108)' ',OUTCON(NCONST) ENDIF 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C IF(NOUT.EQ.0)GOTO 24 C C...VARIANCE-COVARIANCE matrix (not printed) and correlation coefficients C (lower triangle only) C IF(KCOPY-1.LT.1)GOTO 31 SDEVSQ=SWDEV*SWDEV DO 200 N=1,K1 DO 200 L=N,K1 H(N,L)=DD(N,L)*SDEVSQ 200 CONTINUE DO 201 N=1,K1 DO 201 L=1,N H(N,L)=H(L,N)/(C(L)*C(N)) 201 CONTINUE C WRITE(NOUT,202) 202 FORMAT(/' CORRELATION COEFFICIENTS:') M=1 305 MM=M+7 IF(MM.GT.K1)MM=K1 WRITE(NOUT,301) WRITE(NOUT,230)(T(TPT(L)),L=M,MM) 230 FORMAT(7X,8(3X,A5,1X)) WRITE(NOUT,301) DO 302 N=M,K1 LLLL=MM IF(LLLL.GT.N)LLLL=N WRITE(NOUT,303)T(TPT(N)),(H(N,L),L=M,LLLL) 302 CONTINUE 303 FORMAT(1X,A5,1X,8(F9.4)) IF(MM.GE.K1)THEN GOTO 304 ELSE M=M+8 GOTO 305 ENDIF C C Contributions to frequencies (disk only), F descriptor in output format C is chosen individually for each derivative according to its magnitude C C...values of contributions 304 IF(KCOPY-2.LT.1)GOTO 31 DO 203 N=1,I DO 203 L=1,K1 D(N,L)=D(N,L)*A(TPT(L)) 203 CONTINUE C C...heading WRITE(NOUT,204) 204 FORMAT(/' CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO FREQUENCY:'/) IST=1 ISP=8 390 IF(ISP.GT.K1)ISP=K1 WRITE(NOUT,230)(T(TPT(L)),L=IST,ISP) IF(ISP.LT.K1)THEN IST=IST+8 ISP=ISP+8 GOTO 390 ENDIF WRITE(NOUT,301) 301 FORMAT(1X) C DO 320 N=1,I C FORMAL=F1 IST=1 315 ISP=IST+7 IF(ISP.GT.K1)ISP=K1 314 DO 310 L=IST,ISP C C...determine LLL in format descriptor F9.LLL such that largest number C of significant digits are displayed for the current contribution LLL=5 DER=ABS(D(N,L)) DIV=10. 312 IF(DER/DIV.LT.1.)GOTO 311 DIV=DIV*10. LLL=LLL-1 GOTO 312 311 IF(LLL.LT.0)LLL=0 LLLL=L IF(LLLL.GT.8)LLLL=L-IST+1 FORMAR(13+LLLL*5)=CHAR(LLL+ICHAR('0')) 310 CONTINUE IF(ISP.LE.8)THEN WRITE(NOUT,FORMAL)N,(D(N,L),L=IST,ISP) ELSE WRITE(NOUT,FORMAL)(D(N,L),L=IST,ISP) ENDIF IF(ISP.EQ.K1)GOTO 320 FORMAL=F2 IST=IST+8 GOTO 315 C 320 CONTINUE C C 31 NOUT=0 WRITE(3,'(1X)') write(3,53) CLOSE(3) GOTO 23 c 30 WRITE(NOUT,405)SDEV,SWDEV IFLAG=1 24 IF(IOUT.NE.1.AND.IOUT.NE.2)GOTO 300 500 IOUT=0 GOTO 16 C 300 WRITE(*,22) 22 FORMAT(/' HARD COPY OUTPUT - CHOOSE ONE OF:'// 1 20X,'0 = no hardcopy'/ 1 20X,'1 = standard output'/ 1 20X,'2 = so. + correlation coeff.'/ 1 20X,'3 = so. + cc. + contributions ... ',$) READ(*,15,ERR=300)KCOPY IF(KCOPY.LT.1.OR.KCOPY.GT.3)GOTO 23 NOUT=3 1901 WRITE(*,17) READ(*,'(A)',ERR=1901)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) GOTO 25 C C...output of the data file C 23 WRITE(*,14) 14 FORMAT(/' DATA TO BE WRITTEN ON DISK (1/0) ? ',$) READ(*,15,ERR=23)L 15 FORMAT(I5) IOUT=1 IF(L.NE.1)GOTO 24 100 WRITE(*,17) 17 FORMAT(/' NAME TO BE USED FOR THE OUTPUT FILE : ',$) READ(*,'(A)',ERR=100)FILNAM OPEN(2,FILE=FILNAM,STATUS='UNKNOWN',ERR=100) WRITE(2,51)COM write(2,53) do 101 ncon=nconst,1,-1 if(a(ncon).ne.0.d0)goto 102 101 continue 102 WRITE(2,54)'number of transitions:',I write(2,54)'number of constants :',ncon write(2,53) do 103 k=1,ncon if(k.eq.1)then WRITE(2,150)IA(K),A(K),'B.z = A or C' goto 103 endif if(k.eq.3)then WRITE(2,'(58(1H ))') goto 103 endif WRITE(2,150)IA(K),A(K),T(K)(1:5) 103 continue write(2,55) write(2,56) write(2,55) nspace=1 ncomm=2 DO 3400 KK=1,I if(ispace(nspace).eq.kk)then nspace=nspace+1 write(2,'(1x)') endif c 133 if(icomml(ncomm).eq.kk)then write(2,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 133 endif c J=JK(KK,4) IF(NDEG(KK,1).EQ.1)THEN KKK=JK(KK,4)-(JK(KK,5)+JK(KK,6)) K=ISIGN(1,KKK)*JK(KK,5) WRITE(2,140)J,K,FR(KK,1),ERROR(KK),LINFIT(KK) ELSE IF(NDEG(KK,2).EQ.NDEG(KK,3))THEN WRITE(2,2141)J,NDEG(KK,2),FR(KK,1),ERROR(KK), * LINFIT(KK) ELSE WRITE(2,2140)J,NDEG(KK,2),NDEG(KK,3),FR(KK,1),ERROR(KK), * LINFIT(KK) ENDIF 3400 CONTINUE write(2,53) CLOSE(2) 51 FORMAT(1X,72A1) 52 FORMAT(I5) 53 format(79(1H-)) 54 format(a22,i5) 55 format(58(1H-)) 56 format(' J'''' K frequency error', * ' fit ?') 140 FORMAT(2I7,F20.6,F15.6,I5) 2140 FORMAT(I7,I4,'-',I2,F20.6,F15.6,I5) 2141 FORMAT(I7,' D ',I2,F20.6,F15.6,I5) 150 FORMAT(I5,F30.19,3X,5A) C IOUT=1 GOTO 24 C 16 RETURN END C_____________________________________________________________________________ C C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR) C C Constant and error formating for output C IMPLICIT INTEGER*2 (I-N) CHARACTER*25 CONVAL,ERVAL,CONOUT,EROUT C NERD=2 NDCON=0 NDEROR=0 NDNOTZ=0 ICZERO=0 C C...Go through digits of constant value adding those to output buffer while C checking at the same time digits of the error value, reacting as C necessary. Terminate when either three digits in the error value are C transferred or, if error has more than three digits before the decimal C point, the decimal point is reached. DO 1 N=1,25 NDCON=NDCON+1 CONOUT(NDCON:NDCON)=CONVAL(N:N) C C...ensure that zero precedes the decimal point and use ICZERO to monitor C whether decimal point has been reached IF(CONVAL(N:N).EQ.'.')ICZERO=1 if(ndcon.gt.1.and.n.gt.1)then IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.' ') * CONOUT(NDCON-1:NDCON-1)='0' endif if(ndcon.gt.2.and.n.gt.1)then IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.'-')THEN CONOUT(NDCON-2:NDCON-2)='-' CONOUT(NDCON-1:NDCON-1)='0' ENDIF endif C C...use NDNOTZ (number of digits not zero) to monitor whether significant C digits in constant value have been reached IF(CONVAL(N:N).GE.'1'.AND.CONVAL(N:N).LE.'9') * NDNOTZ=NDNOTZ+1 C C...do not transfer error digit if it is a leading zero, dot or space IF(NDEROR.EQ.0 .AND. (ERVAL(N:N).EQ.' '.OR. * ERVAL(N:N).EQ.'0'.OR.ERVAL(N:N).EQ.'.') )GOTO 1 C C...exit if error larger than value and decimal point reached IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 C C...do not transfer the dot in error value IF(ERVAL(N:N).EQ.'.')GOTO 1 C C...transfer error digits until three or, if more, the first significant C digit in value is reached NDEROR=NDEROR+1 EROUT(NDEROR:NDEROR)=ERVAL(N:N) IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 1 CONTINUE 2 CONVAL(1:NDCON)=CONOUT(1:NDCON) if(nderor.lt.6.and.nderor.ge.1)then ERVAL(1:NDEROR)=EROUT(1:NDEROR) else if(nderor.eq.0)then nderor=1 erval(1:1)='*' endif if(nderor.ge.6)then nderor=5 ERVAL(1:5)='*****' endif endif C RETURN END C_____________________________________________________________________________ C C SUBROUTINE ASREDS(I,J,N,A) C C Evaluation of a Wang submatrix for reduction S, where: C C I - on output contains the dimensionality of te calculated subm. C J - J quantum number for the subm. C N - number of the required subm. = 1,2,3,4 for E+,E-,O+,O- resp. C A - array containing rotational and cd. constants in usual order C C Evaluated submatrix is placed in H C IMPLICIT REAL*8(A-H,O-Z) INTEGER*2 I,J,N PARAMETER (NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) c COMMON /BLKH/H/BLKR/U DIMENSION A(NCONST),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C RJ=J Q=(RJ+1.d0)*RJ QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C C=0.5D0*(A(2)+A(2)) 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(2))+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 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 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 M=K*K FK2=M 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) 6 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 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 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_____________________________________________________________________________ C SUBROUTINE HDIAG(N,IEGEN,NR) C C EIGENVALUES AND EIGENVECTORS OF A SQUARE, SYMMETRIC ARRAY C C N - dimensionality of array to be diagonalized C IEGEN - if =0 then calculate eigenvectors, if not then do not C NR - on output contains the number of iterations that were carried out C C Matrix to be diagonalized is placed in H C Eigenvalues are output on the diagonal of H, in ascending order C Eigenvectors are output in U in order of eigenvalues C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (JMAX=250, MAXDIM=JMAX/2+1) COMMON /BLKH/H/BLKR/U DIMENSION H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),X(MAXDIM),IQ(MAXDIM) C IF(IEGEN.NE.0)GOTO 15 DO 14 I=1,N DO 14 J=1,N IF(I-J.NE.0)GOTO 12 U(I,J)=1.D0 GOTO 14 12 U(I,J)=0.D0 14 CONTINUE 15 NR=0 IF(N-1.LE.0)GOTO 1000 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)).GT.0.D0)GOTO 30 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE 40 DO 70 I=1,NMI1 IF(I-1.LE.0)GOTO 60 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,148 148 CONTINUE NR=NR+1 130 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 HTEMP=SIGN(1.D0,-SINE)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 CONTINUE DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV.EQ.0)GOTO 240 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) 250 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)).GT.0.D0)GOTO 320 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)).GE.0.D0)GOTO 390 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)).GE.0.D0)GOTO 530 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)).GE.0.D0)GOTO 450 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)).GE.0.D0)GOTO 530 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)).GE.0.D0)GOTO 500 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)).GE.0.D0)GOTO 530 510 X(JPIV)=ABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF(IEGEN.NE.0)GOTO 40 DO 550 I=1,N HTEMP=U(I,IPIV) U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV) U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV) 550 CONTINUE GOTO 40 1000 CONTINUE C RETURN END c_____________________________________________________________________________ C C Actual length of a string (equivalent to LEN_TRIM) C integer*2 function lentrm(carg) character carg*(*) integer*2 nn,n c nn=len(carg) do 1 n=nn,1,-1 if(carg(n:n).gt.' ')goto 2 1 continue 2 lentrm=n c return end C_____________________________________________________________________________ C_____________________________________________________________________________