C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C ASFIT - ASYMMETRIC ROTOR FITTING PROGRAM, using Watson's reduced C Hamiltonian with terms up to decadic. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Principal features: C C - J UP TO JMAX=250 C - ALL CONSTANTS IN WATSON'S REDUCED HAMILTONIAN C UP TO AND INCLUDING DECADIC CAN BE FITTED FOR REDUCTIONS C 'A' AND 'S' AND REPRESENTATIONS 'I.r' AND 'III.r'. C - UP TO NLINES=15000 TRANSITIONS C - WEIGHTED FIT C - OPTIONAL CORRELATION COEFFICIENTS, MOST INFLUENTIAL TRANSITIONS C AND CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO TRANSITION C FREQUENCIES (EVALUATED FROM DERIVATIVES) C - FIT WITH ZERO DEGREES OF FREEDOM (ie. N constants to N C frequencies) IS ALLOWED, BUT THE USER IS WARNED BY THE BELL C - OPTION TO FIT PEAK FREQUENCIES OF BLENDED LINES C - GENERATION OF QUANTUM NUMBERS OF COMPONENT TRANSITIONS FOR C MANY BAND TYPES C - ON-LINE DATA MODIFICATION (constants, transitions) DURING C THE FIT C - POSSIBILITY OF ANNOTATING EACH LINE WITH AN 8-BYTE COMMENT as C well as allowance for spacing/comment lines between C declared transitions C C This program, though much developed, owes much to the excellent ASFIT C of P.J.Mjoberg et al (ca 1972), which also seems to have carried C code written in E.B.Wilson's group. C C ver. 22.I.2024 ----- 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 1979: port from IBM-360 to PDP-11 C 1980's: MERA-60 -> HP-150 -> PC's -> PRIME, SUN, SG C 1.06.95: allowance for annotations at ends of lines C 11.11.95: modification of data file creation from keyboard; C generation of oblate-type bands according to n=l/m generalization C 15.11.95: addition of deposition format for frequency table C 11.03.96: nomenclature I(+)(-),II(+)(-) and n=l/m scheme for II(-) C 12.03.96: modified line change option C 28.03.96: further flexibility in printout of derivatives C 10.04.96: explicit declaration of the handedness of representations used C 15.04.96: flagging of non-excluded lines with o-c > 3*exptal error C 11.04.97: input file='old', jmax=250, elimination of excl/annotations clash C 16.04.97: change in .ASF format, in output + debugging C 30.03.00: debugging of output C 7.04.01: overhaul of input/output C 10.01.02: addition of hjj/DFBETAS tests + major overhaul/debugging C 4.04.02: fitting of blended transitions C 4.01.03: various modifications to input and output C 7.04.03: .LIN output for SPFIT (preserving some comments for use by PIFORM) C 3.09.03: elimination of bug in fitting blended lines, line generation by C cloning + many incremental mods, including systematisation of C end-of-line comments and switch to 8-byte length C 12.08.04: modified line counting in data file, .PAR output + other mods C 3.12.04: modified output formats and data checks C 19.01.06: detection of potential split blends in dataset C 14.04.08: mods to SPFIT output C 6.06.08: option of symmetric top quantisation to SPFIT output C 1.12.10: correction of DFBETAS bug for data with blended transitions C 29.05.15: debugging and compatibility with gfortran C 3.09.18: NLINES up to 15000 and associated formatting changes C 27.12.18: rounding of formatted constants and errors using modified C CONFOR and unmodified CONFIX from PIFORM C 21.01.19: allowance for 50 character end of line annotations C 19.05.20: Cleaning up and debugging code closer to Fortran2018 C but COMMON and EQUIVALENCE still remaining C 22.01.24: Debugging and clarifying SPFIT output C C_____________________________________________________________________________ C C This program is standard F77 and compiles accordingly (tested with C quite a few FORTRANs, although not necessarily this particular version). C Any metacommands at the top of the listing beginning with a $ are for C MSFortran only and should be discarded on transfer to a different Fortran. C C Intel Visula Fortran compilation: C C ifort -nopdbfile -nodebug -traceback -arch:IA32 -O3 -Qsave C -fpscomp:filesfromcmd asfit.for C gfortran compilation: C C gfortran -fno-automatic -fbacktrace -O2 asfit.for -o asfit C C_____________________________________________________________________________ C C DATA INPUT: C C A new data file is best generated by editing a previous one, which C should be self explanatory. The separating lines written by ASFIT into C the data set define the columns where the various numbers are to C be written (in right-justified form). C C The program provides its own limited editing facilities for C generating or modifying a data set. Generation from scratch is now rarely C used but the editing facility is still invaluable for making C modifications to the data as a result of a fit, and the program C allows full cyclic operation of C fitting->inspection of results->modification->fitting C C The following input conventions apply: C C - 1/0 convention for YES/NO (usually the answer is read in in fixed C format so an suffices for NO) C - constants are to be supplied in the form: 1,constant value C or: 0,constant value C depending on whether the constant is to be fitted or not C C - transitions are to be specified with the three rotational quantum C numbers J,K-1,K+1 of the upper state first, followed by the q.n's C of the lower state, the frequency, its error and fitting parameter: C 3,1,3,2,1,2,8669.6161,0.001,1 C the fitting parameter can be one of: C 0=do not fit C 1=fit C 2=fit blends, in which case intensity weighted mean frequency C of several calculated transitions is fitted to a single C measured frequency of the blended line. In this case C relative intensities of blended components follow C the fitting parameter, for example for unresolved C K2 doublet suspected of splitting greater than the frequency C error: C 3,2,2,2,2,1,8669.6161,0.05,2,1.0 C 3,2,1,2,2,0,8669.6161,0.05,2,1.0 C In the file output each component will have the same o-c C frequency, but also its individual relative intensity and C calculated frequency C C - termination of frequency input loop (including insert and C add modes) is by inputting 0's for every parameter eg. C 0,0,0,0,0,0,0,0 or ,,,,,,,, C C - exit from most modification options is by just pressing C the key C C If ASFIT is used to generate a data set from scratch (which may C sometimes be faster) use option 7 to save intermediate stages as well C as just before the first fitting cycle. C C * The number of transitions specified in the data file only serves as C the count of the number of lines at the moment of writing this file C by ASFIT. C On input the program reads up to the end of the data file, C IRRESPECTIVE of the number of transitions in the third line of file. 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_____________________________________________________________________________ C C The available 'A' reduction constants are as follows: 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 The available 'S' reduction constants are as follows: 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 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) C IRED - if set to +-1 then reduction S, otherwise (+-2) reduction A C (positive values representation I.r, negative representation III.r) C K1 - number of fitted constants C FR - information on transitions where C (i,1) observed frequencies C (i,2) calculated frequencies C (i,3) obs-calc C (i,4) relative intensity for blend calculation 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 as to what to do with specified transitions C =0 do not fit, predict only C =1 fit C =2 fit intensity weighted mean of several transitions to a single C experimental frequency 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 DESCR - string array holding up to 50-byte end-of-line annotations C COMBLK - string variable holding all between-the-line comments C 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, the total C number of spacing lines is NSPCES C_____________________________________________________________________________ C C IMPLICIT REAL(8) (A-H,O-Z) c PARAMETER (NLINES=15000, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1, lpersc=25) PARAMETER (maxcbl=30000,maxcom=1000,maxspa=1000) C COMMON /BLKH/H * /BLKR/R * /SORTCC/WORK,IPT COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR * /BLK2/NJA,B * /BLK3/C,D,DD,HJJ * /TBLOCK/COM,T,pinums DIMENSION A(NCONST),B(NCONST),C(NCONST),Z(NCONST), * H(MAXDIM,MAXDIM),R(MAXDIM,MAXDIM),DD(NCONST,NCONST), * IA(NCONST),LINFIT(NLINES),JK(NLINES,6),IPT(NLINES), * NJA(NNJA,4),HJJ(NLINES) DIMENSION FR(NLINES,4),WORK(NLINES) REAL(8) ERRORL(NLINES),WEIGHT(NLINES),D(NLINES,NCONST) character(50) descr(nlines) CHARACTER(6) T(NCONST),TT(NCONST) CHARACTER COM(72) CHARACTER comblk*(maxcbl),pinums(nconst)*3 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c DATA TT/'A =','B =','C =',' DJ =',' DJK =',' DK =', 1 ' dJ =',' dK =','HJ =','HJK =','HKJ =','HK =', 1 'hJ =','hJK =','hK =', 1 ' LJ =',' LJJK=',' LJK =',' LKKJ=',' LK =',' lJ =', 1 ' lJK =',' lKJ =',' lK =', 1 'PJ =','PJJK =','PJK =','PKJ =','PKKJ =','PK =', 1 'pJ =','pJJK =','pJK =','pKKJ =','pK ='/ DATA PINUMS/'100','200','300',' 2',' 11',' 20', 1 '401','410',' 3',' 12',' 21',' 30', 1 '402','411','420', 1 ' 4',' 13',' 22',' 31',' 40','403', 1 '412','421','430', 1 ' 5',' 14',' 23',' 32',' 41',' 50', 1 '404','413','422','431','440'/ DO 1234 L=1,NCONST T(L)=TT(L) 1234 CONTINUE IOUT=0 C IFLAG=0 CALL DATAIN ! <--- DATAIN 400 CALL MODIF(N,IFLAG) ! <--- MODIF IF(IOUT.EQ.1)GOTO 41 C I1=0 S=1.D50 IPR=3 C C C START ITERATION C C 14 WRITE(*,6000)SDEV,SWDEV 6000 FORMAT(/' CALCULATION CURRENTLY ON J ---',14X, * 'Last deviations:',F10.4,' (',F6.3,')'/31X,'|'/) I1=I1+1 S1=S DO 1015 K=1,I ! Fortran2018 FR(K,2)=0.D0 FR(K,3)=0.D0 DO 15 L=1,NCONST D(K,L)=0.0 15 CONTINUE 1015 CONTINUE ! Fortran2018 L1=1 M4=0 WRITE(*,801,advance='NO')NJA(1,1) C C Go through all sorted energy levels C njs=0 DO 25 K=2,N KK=K-1 IF(NJA(K,1).EQ.NJA(KK,1))GOTO 250 c njs=njs+1 if((njs/15)*15.eq.njs)write(*,'(1x)') WRITE(*,801,advance='NO')NJA(KK,1) ! counter 801 FORMAT(1x,I4) c 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 diagonalise the appropriate Wang matrix for each group of C common levels C IF(IABS(int(IRED)).EQ.1)THEN CALL ASREDS(M,J1,N1,A,IRED) ! <--- ASREDS ELSE CALL ASREDA(M,J1,N1,A,IRED) ! <--- ASREDA ENDIF IEGEN=0 CALL HDIAG(M,IEGEN,K1) ! <--- HDIAG CALL HSORT(M) ! <--- HSORT DO 17 L=L1,L2 M=NJA(L,3) M1=NJA(L,4) M2=ABS(M1) FR(M2,2)=FR(M2,2)+SIGN(1,INT(M1))*H(M,M) 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 IF(IABS(INT(IRED)).EQ.1)THEN CALL ASREDS(M3,J1,N1,B,IRED) ! <--- ASREDS ELSE CALL ASREDA(M3,J1,N1,B,IRED) ! <--- ASREDA ENDIF DO 124 L=L1,L2 M=NJA(L,3) IF(L.NE.L1.AND.M.EQ.M4)GOTO 23 EE=0.0D0 DO 22 M1=1,M3 DO 122 M2=M1,M3 N2=2 IF(M1.EQ.M2)N2=1 EE=EE+dble(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,INT(M1))*EE M4=M 124 CONTINUE 24 CONTINUE L1=K 25 CONTINUE C C FIND NEW CONSTANTS - inverse of D'D is evaluated by diagonalisation 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 C D'D is placed and diagonalised in the top lh. corner of H, C D'Y is placed in (K1+1)'st column of H C inv(D'D) is placed in array DD C S=0.D0 SW=0.D0 L=0 DO 26 J=1,I C C...unfitted line IF(LINFIT(J).eq.0)then FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) endif C C...single line if(linfit(j).eq.1)then FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERRORL(J))**2 L=L+1 endif C C...blend: deal with the whole blend at the first component if(linfit(j).eq.2)then if(j.eq.1.or.(j.gt.1.and. * ( (fr(J-1,1).ne.fr(j,1)).or.linfit(j-1).ne.2 ) ) )then jfirst=j jlast=j over=fr(j,4)*fr(j,2) below=fr(j,4) do 526 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 527 jlast=jjj over=over+fr(jjj,4)*fr(jjj,2) below=below+fr(jjj,4) 526 continue 527 fr(j,3)=fr(j,1)-over/below do 528 jjj=jfirst+1,i fr(jjj,3)=fr(j,3) if(jjj.eq.jlast)goto 529 528 continue 529 S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERRORL(J))**2 L=L+1 endif endif 26 CONTINUE C NINFIT=L 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(*,'(1x)') 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) C C...in the multiplications below also deal with effect of the whole blend C on encountering the first component (and successive components the same C blend are then ignored) C DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(LINFIT(L).eq.0)GOTO 27 IF(LINFIT(L).eq.1)E=E+D(L,J)*D(L,K)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*D(L,K)*WEIGHT(L) below=FR(L,4) do 530 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 531 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*D(jjj,K)*WEIGHT(jjj) below=below+fr(jjj,4) 530 continue 531 E=E+over/below endif ENDIF 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).eq.0)GOTO 29 IF(LINFIT(L).eq.1)E=E+D(L,J)*FR(L,3)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*FR(L,3)*WEIGHT(L) below=FR(L,4) do 532 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 533 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*FR(jjj,3)*WEIGHT(jjj) below=below+fr(jjj,4) 532 continue 533 E=E+over/below endif ENDIF 29 CONTINUE H(J,K1+1)=E 30 CONTINUE C DO 1301 J=1,K1 ! Fortran2018 DO 301 K=J,K1 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 301 CONTINUE 1301 CONTINUE ! Fortran2018 C IEGEN=0 CALL HDIAG(K1,IEGEN,L3) ! <--- HDIAG C 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) ? '')', * advance='NO') READ(*,5,ERR=403)IFLAG IF(IFLAG.NE.1)GOTO 6 C 404 WRITE(*,401,advance='NO') 401 FORMAT(1X/ * 9x,'= 0 OBS-CALC + diagnostics = 2 DISK OUTPUT'/ * 9x,'=-1 as above, screen at a time = 3 NEXT ITERATION'/ * 9x,'= 1 MODIFICATIONS = 4 EXIT ASFIT', * 8X,'.... ') 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 lmax=nlines IF(IFLAG.LT.0)lmax=lpersc-1 WRITE(*,'(1X/1x,78(1H-)/1x,t17,''TRANSITION'', * T37,''Frequency'',T52,''Obs-Calc'',T63,''Error'', * T72,''Notes''/1x,78(1H-))') nspace=1 ncomm=2 lcount=3 do 108 j=1,i c if(ispace(nspace).eq.j)then nspace=nspace+1 write(*,'(1x)') lcount=lcount+1 if(lcount.eq.lmax)then write(*,516,advance='NO') read(*,'(i1)',err=515)jj if(jj.eq.1)goto 507 515 lcount=0 endif endif 516 format( * ' Press to continue, 1 to stop .... ') c 501 if(icomml(ncomm).eq.j)then write(*,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 lcount=lcount+1 if(lcount.eq.lmax)then write(*,516,advance='NO') read(*,'(i1)',err=506)jj if(jj.eq.1)goto 507 506 lcount=0 goto 501 endif goto 501 endif c IF(LINFIT(J).gt.0)then if(abs(fr(j,3)).le.3.d0*errorl(j))then if(linfit(j).eq.2) * WRITE(*,93)J,' ',(JK(J,L),L=1,6),' ',FR(J,1),' B', * FR(J,3),ERRORL(J),DESCR(J)(1:10) if(linfit(j).eq.1) * WRITE(*,9)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3),ERRORL(J), * DESCR(J)(1:10) else if(linfit(j).eq.2) * WRITE(*,93)J,'?',(JK(J,L),L=1,6),' ?',FR(J,1),' B', * FR(J,3),ERRORL(J),'>3*Error' if(linfit(j).eq.1) * WRITE(*,92)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3),ERRORL(J), * '>3*Error' endif else WRITE(*,91)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3), * '--excl' endif lcount=lcount+1 if(lcount.eq.lmax)goto 500 goto 108 c 500 write(*,516,advance='NO') read(*,'(i1)',err=505)jj if(jj.eq.1)goto 507 505 lcount=0 c 108 continue if(iflag.lt.0.and.lcount.gt.lmax-12)then write(*,516,advance='NO') read(*,'(i1)',err=507)jj endif c 507 WRITE(*,78)I1,SDEV,SWDEV c c...most influential lines CALL hjjset(K1,hjjav) ! <--- hjjset WRITE(*,101) nhj=i-5 if(nhj.lt.1)nhj=1 write(*,102)(ipt(L),hjj(ipt(L)),L=I,nhj,-1) c c...worst fitting lines write(*,105) do 106 L=1,I if(LINFIT(L).gt.0)then work(L)=abs(fr(L,3)/errorl(L)) if(L.ne.1)then if( (linfit(L).eq.2.and.linfit(L-1).eq.2) .and. * (fr(L,1).eq.fr(L-1,1)) )work(L)=0.d0 endif else WORK(L)=0.d0 endif ipt(L)=L 106 continue int21=1 call sortc(int21,I) ! <--- sortc nhj=I-5 if(nhj.lt.1)nhj=1 write(*,110)(ipt(L),FR(ipt(L),3)/errorl(ipt(L)),L=I,nhj,-1) c GOTO 404 ENDIF 9 FORMAT(1X,I5,'. ', 2(I5,2I3), F16.4, F14.4,F8.4,1x,a) 91 FORMAT(1X,I5,'.--',I4,2I3,I5,2I3,'--',F14.4, F14.4,a) 92 FORMAT(1X,I5,'.?', 2(I5,2I3),' ?',F14.4, F14.4,F8.4,1x,a) 93 FORMAT(1X,I5,'.',A, 2(I5,2I3),A, F14.4,a,F12.4,F8.4,1x,a) 101 FORMAT(1x/' MOST INFLUENTIAL TRANSITIONS (line number=hjj):') 102 FORMAT(1x,6(I7,'=',F5.3)) 105 FORMAT(' WORST FITTING LINES (line number=(o-c)/error):') 110 FORMAT(1x,6(I7,'=',F5.1)) c c...disk output c IF(IFLAG.EQ.2)THEN IOUT=2 CALL DATOUT(SDEV,SWDEV,NINFIT) ! <--- DATOUT 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,advance='NO') 450 FORMAT(/' MORE CYCLES ? ') READ(*,'(I5)',ERR=1900)L IF(L.NE.1)GOTO 41 I1=0 GOTO 14 C 41 CALL DATOUT(SDEV,SWDEV,NINFIT) ! <--- DATOUT C IF(IFLAG.EQ.1)GOTO 400 GOTO 99 48 WRITE(*,89)L,K1,' - ITERATION STOPPPED' 89 FORMAT(//' ONLY',I3,' MEASURED TRANSITIONS FOR',I3, 1 ' CONSTANTS',A) GOTO 41 C 99 STOP END C C_____________________________________________________________________________ C SUBROUTINE DATAIN IMPLICIT REAL(8) (A-H,O-Z) PARAMETER (NLINES=15000, NCONST=35, JMAX=250) PARAMETER (maxcbl=30000,maxcom=1000,maxspa=1000) C COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR * /TBLOCK/COM,T,PINUMS DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST), * ERRORL(NLINES),WEIGHT(NLINES), * IA(NCONST),LINFIT(NLINES) character(50) descr(nlines) CHARACTER(6) T(NCONST) CHARACTER COM(72) CHARACTER FILNAM*25,cdummy*22,line*140,errmes*50 CHARACTER comblk*(maxcbl),pinums(nconst)*3 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c IRED=2 C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | ASFIT - FITTING OF ASYMMETRIC TOP ROTATIONAL SPECTRUM ', * T79,'|'/ * ' | (reduction A or S, all terms up to decadic)', * T79,'|'/ * ' |',76(1H_),'|'/' version 22.I.2024',T64,'Zbigniew KISIEL'/) write(*,1910)nlines,jmax 1910 format(' Current limits:',i7,' lines and J up to',i4) c 1900 WRITE(*,12,advance='NO') 12 FORMAT(//' DATA TO BE CREATED (+-1), OR READ FROM DISK (+-2,+-3):' 1 //10x,'positive values - reduction A'/ 1 10x,'negative values - reduction S'/ 1 10x,' +-2 - prolate-type representation I.r'/ 1 10x,' +-3 - oblate-type representation III.r'// 1 10x,' ..... ') READ(*,14,ERR=1900)IFLAG 14 FORMAT(2I5) C C...IRED=2 -> reduction 'A', IRED=1 -> reduction 'S' C IF(IFLAG.GT.0)GOTO 31 IRED=1 if(iflag.eq.-3)ired=-1 T(7)= ' d1 =' T(8)= ' d2 =' T(13)='h1 =' T(14)='h2 =' T(15)='h3 =' T(21)='l1 =' T(22)='l2 =' T(23)='l3 =' T(24)='l4 =' T(31)='p1 =' T(32)='p2 =' T(33)='p3 =' T(34)='p4 =' T(35)='p5 =' DO 200 L=21,NCONST A(L)=0.0D0 IA(L)=0 200 CONTINUE c 31 if(iflag.eq.3)ired=-2 IF(ABS(IFLAG).gt.1)GOTO 15 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C c c...Input from keyboard (IFLAG=+-1) c 1901 WRITE(*,16) 16 FORMAT(/' COMMENT :') READ(*,17,ERR=1901)COM 17 FORMAT(72A1) 142 write(*,40,advance='NO') 40 format(/' NUMBER OF CONSTANTS TO BE SPECIFIED : ') read(*,'(i2)',err=41)nfin if(nfin.lt.1.or.nfin.gt.nconst)goto 142 WRITE(*,18) 18 FORMAT(/' FITTING PARAMETER, INITIAL VALUE OF CONSTANT : ' 1 /) DO 19 L=1,nfin 1903 WRITE(*,20,advance='NO')T(L) 20 FORMAT(1X,A6,2X) READ(*,*,ERR=1902)IA(L),A(L) GOTO 19 1902 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1903 19 CONTINUE I=0 WRITE(*,22) 22 FORMAT(/' LINE ASSIGNMENTS, i.e. six comma separated qn''s, ', 1 'frequency, error: '/ 2 ' (terminate with -1,,,,,,,,)'/) 24 I=I+1 1906 WRITE(*,'('' Line'',i5,'': '')',advance='NO')I READ(*,*,ERR=1904)(JK(I,K),K=1,6),FR(I,1),ERRORL(I) descr(i)=' ' IF(ERRORL(I).EQ.0.0d0)ERRORL(I)=1.0d0 WEIGHT(I)=1./ERRORL(I)**2 GOTO 1905 1904 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1906 1905 IF(JK(I,1).GE.0)GOTO 24 I=I-1 DO 70 L=1,I LINFIT(L)=1 70 CONTINUE GOTO 41 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C C...input from disk (some features are added so that in DOS compilation there C are no problems with reading UNIX files and files from botched UNIX/DOS C conversion ) C 15 WRITE(*,135,advance='NO') 135 FORMAT(/' NAME OF DATA FILE : ') READ(*,'(A)',ERR=15)FILNAM OPEN(2,FILE=FILNAM,ERR=153,STATUS='OLD') c errmes='Bad comment line' READ(2,'(a)',err=152,end=152)line READ(line,3500,ERR=151)COM c errmes='Cannot skip to the number of lines' READ(2,'(a)',err=152,end=152)line 3500 FORMAT(1X,72A1) c errmes='Bad number of transitions' READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 if(len_trim(line).le.27)line=line(1:len_trim(line))//' ' READ(line(1:27),141,ERR=151)cdummy,I if(I.le.0)then errmes='Zero or negative number of transitions' goto 151 endif c errmes='Bad number of constants' READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 if(len_trim(line).le.27)line=line(1:len_trim(line))//' ' read(line,141,err=151)cdummy,NCON if(ncon.le.0.or.ncon.gt.nconst)then errmes='Zero, negative or excessive number of constants' goto 151 endif 141 format(a22,i5) errmes='Cannot skip to the list of constants' READ(2,'(a)',err=152,end=152)line c if(ncon.eq.0)ncon=NCONST do 159 k=1,ncon write(errmes,'(a,i3)')'Bad fitting flag or constant number',k READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 READ(line,150,ERR=151)IA(K),A(K) if(ia(k).lt.0.or.ia(k).gt.1)ia(k)=0 159 continue 150 FORMAT(I5,F30.19) c c...compatibility with files without transitions header c read(2,'(a)',err=151)cdummy if(cdummy(2:11).eq.'----------')then read(2,3500,err=151)cdummy read(2,3500,err=151)cdummy else backspace(2) endif c c...input transitions c nlin=0 nspace=1 ncomm=1 icomm(ncomm)=0 do 3400 L=1,NLINES kkk=nlin+1 c 23 read(2,'(a)',err=241,end=51)line if(icheck(line).ne.0)goto 157 if(len_trim(line).le.70) * line=line(1:len_trim(line))//' ' c c...deal with a between-the-line comment 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)+len_trim(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)=L icomm(ncomm)=icomm(ncomm-1)+len_trim(line) comblk(icomm(ncomm-1)+1:icomm(ncomm))=line(1:len_trim(line)) goto 23 endif if(line(1:3).eq.'end')goto 51 if(line(6:10).eq.'-----')goto 23 c c...deal with a separating line if(len_trim(line).lt.20)then nspace=nspace+1 if(nspace.lt.maxspa)ispace(nspace)=L if(nspace.gt.1.and.i.eq.ispace(nspace-1))nspace=nspace-1 goto 23 endif c c...read actual information READ(LINE,140,ERR=241)(JK(kkk,M),M=1,6),FR(kkk,1),ERRORL(kkk), * LINFIT(kkk) c read(line(71:80),'(f10.5)',err=3405)fr(kkk,4) if(linfit(kkk).ne.2)fr(kkk,4)=0.d0 goto 3402 3405 fr(kkk,4)=0.d0 c 3402 if(len_trim(line).gt.80)then read(line(81:),'(a)')descr(kkk) else descr(kkk)='' endif c if(linfit(kkk).lt.0.or.linfit(kkk).gt.2)linfit(kkk)=0 140 FORMAT(6I5,F20.6,F15.6,I5) if(linfit(kkk).eq.2.and.fr(kkk,4).eq.0.d0)fr(kkk,4)=1.d0 nlin=nlin+1 3400 CONTINUE 51 I=NLIN if(nlin.eq.nlines)then write(*,3420,advance='NO')nlines 3420 FORMAT(1x/' ***** WARNING: The current maximum of',i7, * ' on the number of lines'/ * ' has been reached.'/ * ' Some lines may not have been read ', * 'and may be lost if'/' the data file is ', * 'overwritten during the fit.'// * ' Program can continue but you are advised to rectify this'// * /16x,'Press ENTER to continue ') read(*,'(i5)',err=1057)k 1057 continue endif nspces=nspace c CLOSE(2) c c...Check whether blends are specified correctly c nberr=0 nblend=0 bfreq=0.d0 DO 1050 k=1,I c if(linfit(k).ne.2)then if(nblend.eq.0)goto 1050 if(nblend.eq.1)then nberr=nberr+1 if(nberr.eq.1)write(*,'(2x)') write(*,1051)k-1,fr(k-1,1) 1051 format(' **** WARNING: line',i6,F16.5,' is only a single', * ' blend component') nblend=0 bfreq=0.d0 goto 1050 endif if(nblend.ge.2)then nblend=0 bfreq=0.d0 goto 1050 endif endif c if(linfit(k).eq.2.and.nblend.eq.0)then nblend=1 bfreq=fr(k,1) goto 1050 endif c if(linfit(k).eq.2.and.nblend.ne.0)then if(fr(k,1).eq.bfreq)then nblend=nblend+1 goto 1050 else if(nblend.eq.1)then write(*,1051)k-1,fr(k-1,1) nberr=nberr+1 nblend=1 bfreq=fr(k,1) goto 1050 else nblend=1 bfreq=fr(k,1) goto 1050 endif endif endif c 1050 continue if(nberr.gt.0)then write(*,1052,advance='NO') 1052 format(1x/' ----> ADVICE: Program will continue but you mig', * 'ht want to check the dataset'/15x,'Press ENTER to continue ') read(*,'(i5)',err=1053)k 1053 continue endif c c...check for duplicated lines c ndupl=0 DO 1060 kk=1,i-1 if(linfit(kk).eq.0)goto 1060 DO 1062 k=kk+1,i if(linfit(k).eq.0)goto 1062 do 1061 kkk=1,6 if(jk(kk,kkk).ne.jk(k,kkk))goto 1062 1061 continue ndupl=ndupl+1 if(ndupl.eq.1)write(*,'(2x)') write(*,'('' **** WARNING: duplicated lines ''/ * 9x,i5,'':'',3i4,2x,3i4,f16.4,f10.4/ * 9x,i5,'':'',3i4,2x,3i4,f16.4,f10.4)') * k,(jk(k,kkk),kkk=1,6),fr(k,1),errorl(k), * kk,(jk(kk,kkk),kkk=1,6),fr(kk,1),errorl(kk) 1062 CONTINUE 1060 CONTINUE if(ndupl.gt.0)then write(*,1052) read(*,'(i5)',err=1063)k 1063 continue endif c c...check for split blends c nsplit=0 DO 1055 k=1,I-1 if(fr(k,1).eq.0.d0)goto 1055 if(linfit(k).eq.0)goto 1055 berr=errorl(k) C DO 1054 kk=k+1,I if(fr(kk,1).eq.0.d0)goto 1054 if(linfit(kk).eq.0)goto 1054 berra=dsqrt(errorl(kk)**2+berr**2) c if(abs(fr(k,1)-fr(kk,1)).le.berra)then C IF(fr(k,1).eq.fr(kk,1).and.kk-k.eq.1)goto 1054 idiff=0 do 1056 kkk=k+1,kk-1 if(abs(fr(kkk,1)-fr(kk,1)).gt.berra)idiff=idiff+1 1056 continue if(idiff.eq.0)goto 1054 c nsplit=nsplit+1 if(nsplit.eq.1)write(*,'(2x)') write(*,1058) k, (jk(k,j), j=1,6),fr(k,1),errorl(k), * kk,(jk(kk,j),j=1,6),fr(kk,1),errorl(kk) endif c 1054 CONTINUE 1055 CONTINUE 1058 format(' **** WARNING - possible split blend: '/ * 9x,i5,':',3i4,2x,3i4,f16.4,f10.4/ * 9x,i5,':',3i4,2x,3i4,f16.4,f10.4) c if(nsplit.gt.0)then write(*,1052) read(*,'(i5)',err=1059)k 1059 continue endif c c...weights c DO 1020 K=1,I IF(ERRORL(K).EQ.0.0d0)ERRORL(K)=1.0d0 WEIGHT(K)=1./ERRORL(K)**2 1020 CONTINUE GOTO 41 C 153 write(*,'(1x/ '' **** ERROR: nonexistent file, try again'')') goto 15 c 152 write(*,'(1x//'' **** ERROR: '',a// * '' **** END OF FILE encountered on line:''/1x,a//)') * errmes,line(1:len_trim(line)) stop c 151 write(*,'(1x//'' ***** ERROR: '',a// * '' reading line ---''/17x,''|''/1x,a//)') * errmes,line(1:len_trim(line)) if(len_trim(line).eq.0)write(*,'('' (This is an empty line)''//)') stop c 157 write(*,'(1x//'' ***** ERROR: found an embedded non-FORTRAN '', * ''character with ASCII code ='',i4/14x, * ''- please edit this character out of the data set''// * '' reading line ---''/17x,''|''/1x,a//)') * icheck(line),line(1:len_trim(line)) stop c 241 WRITE(*,'(1X/'' ***** ERROR: bad transition number'',I5// * '' reading line ---''/17x,''|''/1x,a//)') * nlin+1,line(1:len_trim(line)) if(len_trim(line).eq.0)write(*,'('' (This is an empty line)''//)') STOP c WRITE(*,'(1X/ * '' ***** ERROR: bad intensity of blend component: '',A// * '' reading line ---''/17x,''|''/1x,a)') * descr(kkk),line STOP C 41 RETURN END C 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) c PARAMETER (NLINES=15000, JMAX=250, NCONST=35, NCONSTH=17, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1, lpersc=25) PARAMETER (maxcbl=30000,maxcom=1000,maxspa=1000) C COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR * /BLK2/NJA,B * /TBLOCK/COM,T,PINUMS COMMON /BLKH/H * /BLKR/R DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST),IINT(6), * B(NCONST),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),ERRORL(NLINES),WEIGHT(NLINES), * IA(NCONST),LINFIT(NLINES),NJA(NNJA,4),IWORK(MAXDIM) character(50) descr(nlines) CHARACTER(6) T(NCONST) CHARACTER COM(72) character wang(4)*2,WANGL*8,pinums(nconst)*3 EQUIVALENCE (WANGL,WANG(1)) CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c NQNS=6 write(wangl,'(a)')'E+E-O+O-' 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 NBADL=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 C C Test whether (K-1 + K+1) - J = 0 or 1 C J=JK(K,L)+JK(K,L+1)-JK(K,L-1) IF(J.LT.0.OR.J.GT.1)GOTO 7 N=N+1 NJA(N,1)=JK(K,L-1) NJA(N,4)=K IF(L.EQ.5)NJA(N,4)=-K C C...Wang assignment for representation I.r, old compact coding IF(IRED.GT.0)THEN J=J+(JK(K,L)-JK(K,L)/2*2)*2+1 M=J/2*2 NJA(N,2)=J NJA(N,3)=(JK(K,L)-(J-1)/2-(4-J)/2*M)/2+1 C C...Wang assignment for representation III.r ELSE J=JK(K,L)-JK(K,L)/2*2 IF(JK(K,L-1)/2*2.NE.JK(K,L-1))J=1-J J=J+(JK(K,L+1)-JK(K,L+1)/2*2)*2+1 NJA(N,2)=J IF(JK(K,L-1)-JK(K,L-1)/2*2.EQ.1)THEN IPLUS=J+1 IF(J.GE.3)IPLUS=4-J ELSE IPLUS=J-1 IF(J.GE.3)IPLUS=6-J ENDIF NJA(N,3)=1+(JK(K,L)-JK(K,L+1)+JK(K,L-1)-IPLUS)/4 ENDIF 6 CONTINUE 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 NBADL=NBADL+1 if(NBADL.eq.1)write(*,797) 797 FORMAT(1x/' The following lines have errors in quantum ', * 'numbers and have been'/ * ' automatically excluded from fit:'/) WRITE(*,97)K,(JK(K,M),M=1,NQNS),FR(K,1),ERRORL(K) LINFIT(K)=0 10 CONTINUE c 30 FORMAT(/1X,72A1) 31 FORMAT(/' CURRENT CONSTANTS :'/) 33 FORMAT(1X,A6,F28.18,I5) 34 FORMAT(/' TRANSITIONS IN DATA SET:',I6) 97 FORMAT(1X,'Line',I5,': ',2(I5,2I4),F15.4,F10.4, * ' <---INPUT ERROR') c if(nbadl.gt.0)then write(*,1052,advance='NO') 1052 format(1x/' **** WARNING: Program will continue but it is adv', * 'ised to correct the dataset'/15x,'Press ENTER to continue ') read(*,'(i5)',err=1053)k 1053 continue endif C C Listing Selection and execution of required modification C 1900 WRITE(*,40,advance='NO') 40 FORMAT(/' SELECT AN OPTION: = -1 list lines page at a time' * /24x, ' = 0 proceed to fit'/ * 24x, ' = 1 list lines'/ * 24x, ' = 2 modifications .... ') READ(*,14,ERR=1900)IFLAG 14 FORMAT(I5) IF(IFLAG.EQ.0)GOTO 41 c IF(iabs(int(IFLAG)).eq.1)then lmax=nlines+maxcom+maxspa IF(IFLAG.LT.0)lmax=lpersc-1 WRITE(*,'(1X/1x,78(1H-)/1x,t19,''TRANSITION'', * T39,''Frequency'',T53,''Error'', * T65,''Notes''/1x,78(1H-))') nspace=1 ncomm=2 lcount=3 do 708 j=1,i c if(ispace(nspace).eq.j)then nspace=nspace+1 write(*,'(1x)') lcount=lcount+1 if(lcount.eq.lmax)then write(*,716,advance='NO') read(*,'(i1)',err=711)jj if(jj.eq.1)goto 1900 711 lcount=0 endif endif 716 format( * ' Press to continue, 1 to stop .... ') c 701 if(icomml(ncomm).eq.j)then write(*,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 lcount=lcount+1 if(lcount.eq.lmax)then write(*,716,advance='NO') read(*,'(i1)',err=706)jj 706 lcount=0 goto 701 endif goto 701 endif c IF(linfit(J).eq.2) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERRORL(J), * '<==Blended line' IF(LINFIT(J).eq.1) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERRORL(J), * DESCR(J)(1:10) IF(linfit(J).eq.0) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERRORL(J), * '<--Excluded line' 790 FORMAT(1X,I5,': ',2(I5,2I4),F15.4,F10.4,2x,A) lcount=lcount+1 if(lcount.eq.lmax)goto 700 goto 708 c 700 write(*,716,advance='NO') read(*,'(i1)',err=705)jj if(jj.eq.1)goto 1900 705 lcount=0 c 708 continue if(iflag.lt.0.and.lcount.gt.lmax-12)then write(*,716,advance='NO') read(*,'(i1)',err=707)jj endif 707 goto 1900 endif c if(Iflag.ne.2)GOTO 1900 C C...M O D I F I C A T I O N M E N U C 11 WRITE(*,70) 70 FORMAT(/' MODIFICATIONS ARE SELECTED AS FOLLOWS :' 1 //' = 0 LISTING AND EXIT TO CALCULATION',T42, 1 '= 6 CHANGE THE COMMENT'/ 1 ' = 1 CHANGE CONSTANTS',T42,'= 7 STORE DATA'/ 1 ' = 2 DELETE LINES',T42,'= 8 EXCLUDE LINES FROM FIT'/ 1 ' = 3 CHANGE LINES',T42,'= 9 INCLUDE LINES IN FIT'/ 1 ' = 4 ADD LINES',T42,'= 10 CHANGE ERROR'/ 1 ' = 5 INSERT LINES',T42,'= 11 GENERATE LINES'/ 1 1X,T42,'= 12 DUMP EIGENVALUES') 72 WRITE(*,71,advance='NO') 71 FORMAT(/' CONTROL PARAMETER (555=EXIT): ') READ(*,14,ERR=72)ICONTR if(icontr.eq.555)stop IF(ICONTR.LT.0.OR.ICONTR.GT.12)GOTO 11 C if(icontr.eq. 0)goto 29 ! exit if(icontr.eq. 1)goto 74 ! change constant if(icontr.eq. 2)goto 75 ! delete line if(icontr.eq. 3)goto 76 ! change line if(icontr.eq. 4)goto 77 ! add line if(icontr.eq. 5)goto 500 ! insert line if(icontr.eq. 6)goto 178 ! change comment if(icontr.eq. 7)goto 179 ! store data if(icontr.eq. 8)goto 210 ! exclude line if(icontr.eq. 9)goto 220 ! include line if(icontr.eq.10)goto 300 ! change error if(icontr.eq.11)goto 650 ! generate if(icontr.eq.12)goto 630 ! dump c GOTO 72 C C...old: 0 1 2 3 4 5 6 7 8 9 10 11 12 c GOTO(29,74,75,76,77,500,178,179,210,220,300,650,630),ICONTR+1 C C Change of constant (with 1) C 74 WRITE(*,31) LL=NCONSTH 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,advance='NO') 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,advance='NO')T(NN) 79 FORMAT(/' FITTING PARAMETER, NEW VALUE OF ',A6,' ') READ(*,*,ERR=1911)IA(NN),A(NN) GOTO 74 C C Line deletion (with 2) C 75 WRITE(*,81,advance='NO') 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) FR(L,4)=FR(L+1,4) LINFIT(L)=LINFIT(L+1) ERRORL(L)=ERRORL(L+1) WEIGHT(L)=WEIGHT(L+1) DESCR(L)=DESCR(L+1) DO 111 LL=1,NQNS JK(L,LL)=JK(L+1,LL) 111 CONTINUE 1 CONTINUE c do 150 L=1,maxspa if(ispace(L).ge.NN.and.ispace(L).gt.1)ispace(L)=ispace(L)-1 150 continue do 151 L=1,maxcom if(icomml(L).ge.NN.and.icomml(L).gt.1)icomml(L)=icomml(L)-1 151 continue c GOTO 75 C C Change of lines (with 3) C 76 WRITE(*,82,advance='NO') 82 FORMAT(/ * ' NUMBER OF LINE TO BE CHANGED (only frequency if -ve)? ') READ(*,14,ERR=76)NN IF(iabs(int(NN)).LT.1.OR.iabs(int(NN)).GT.I)GOTO 72 if(nn.lt.0)then nn=iabs(int(nn)) isele=1 else isele=0 endif 1912 WRITE(*,90)NN,(JK(NN,L),L=1,NQNS),FR(NN,1),ERRORL(NN) 90 FORMAT(/ ' LINE NO.',I3,' :',2(I5,2I3),F18.4,F14.4) if(isele.eq.0)then write(*,1913,advance='NO') 1913 format(' NEW LINE : ') READ(*,*,ERR=1912)(JK(NN,L),L=1,NQNS),FR(NN,1),ERRORL(NN) IF(ERRORL(NN).EQ.0.0d0)ERRORL(NN)=1.0d0 WEIGHT(NN)=1./ERRORL(NN)**2 else write(*,1914,advance='NO') 1914 format(23x,' NEW FREQUENCY : ') READ(*,*,ERR=1912)FR(NN,1) endif GOTO 76 C C Line addition (with 4) C 77 I=I+1 WRITE(*,190) 190 FORMAT(/' LINE (LINES) TO BE ADDED - qn''s, freq, error:'/ * ' (terminate with -1,,,,,,,,)'/) 191 WRITE(*,'('' Line'',i5,'': '')',advance='NO')I READ(*,*,ERR=1915)(JK(I,L),L=1,NQNS),FR(I,1),ERRORL(I) GOTO 1916 1915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 191 1916 IF(JK(I,1).LT.0)GOTO 192 LINFIT(I)=1 FR(I,4)=0.d0 IF(ERRORL(I).EQ.0.)ERRORL(I)=1. WEIGHT(I)=1./ERRORL(I)**2 DESCR(I)=' ' I=I+1 GOTO 191 192 I=I-1 GOTO 72 C C Line insertion (with 5) C 500 WRITE(*,501,advance='NO') 501 FORMAT(/' INSERTIONS ARE TO PRECEDE LINE NO. : ') READ(*,14,ERR=500)NINS WRITE(*,502) 502 FORMAT(/' LINE (LINES) TO BE INSERTED - qn''s, freq, error:'/ * ' (terminate with -1,,,,,,,,)'/) 503 WRITE(*,'('' Line'',i5,'': '')',advance='NO')NINS READ(*,*,ERR=1918)(IINT(L),L=1,NQNS),FREQ,TEMP GOTO 1919 1918 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 503 1919 IF(IINT(1).lt.0)GOTO 72 I=I+1 DO 1504 L=I,NINS,-1 ! Fortran2018 FR(L,1)=FR(L-1,1) FR(L,4)=FR(L-1,4) LINFIT(L)=LINFIT(L-1) ERRORL(L)=ERRORL(L-1) WEIGHT(L)=WEIGHT(L-1) DESCR(L)=DESCR(L-1) DO 504 LL=1,NQNS JK(L,LL)=JK(L-1,LL) 504 CONTINUE 1504 CONTINUE ! Fortran2018 FR(NINS,1)=FREQ FR(NINS,4)=0.d0 DESCR(NINS)=" " LINFIT(NINS)=1 IF(TEMP.EQ.0.0d0)TEMP=1.0d0 ERRORL(NINS)=TEMP WEIGHT(NINS)=1./TEMP**2 DO 505 LL=1,NQNS JK(NINS,LL)=IINT(LL) 505 CONTINUE NINS=NINS+1 c do 152 L=1,maxspa if(ispace(L).ge.Nins-1)ispace(L)=ispace(L)+1 152 continue do 153 L=1,maxcom if(icomml(L).ge.Nins-1)icomml(L)=icomml(L)+1 153 continue c GOTO 503 C C Change of comment (with 6) C 178 WRITE(*,180)COM 180 FORMAT(/' CURRENT COMMENT :'/1X,A//' NEW COMMENT :') READ(*,181,ERR=178)COM 181 FORMAT(72A1) GOTO 72 C C Saving of data (with 7) C 179 IOUT=1 GOTO 200 C C Line exclusion (with 8) C 210 WRITE(*,211,advance='NO') 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 (with 9) C 220 WRITE(*,221,advance='NO') 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 (with 10) C 300 WRITE(*,301,advance='NO') 301 FORMAT(/' ERROR TO BE CHANGED ON LINE NUMBER ? ') READ(*,'(i5)',ERR=300)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 WRITE(*,302,advance='NO')NN,(JK(NN,K),K=1,NQNS), * FR(NN,1),ERRORL(NN) 302 FORMAT(/' LINE NO.',I3,' :',2(I5,2I3),F18.4,' +-',F14.4/ 1 ' NEW ERROR : ') READ(*,'(f20.10)',ERR=300)TEMP IF(TEMP.EQ.0.0d0)TEMP=1.0d0 ERRORL(NN)=TEMP WEIGHT(NN)=1./TEMP**2 GOTO 300 C C Generation of lines of a given type for prediction or as templates C for input (with 11) C 650 write(*,651,advance='NO') 651 format(/' SELECT 0 = generation according to band rules ', * '(and exit)'/ * ' 1 = generation by cloning a given line'// * ' ..... ') READ(*,14,ERR=650)ICLONE if(iclone.eq.0)goto 600 if(iclone.ne.1)goto 650 c c...clones of a supplied line c iold=i CALL GENER(ICLONE,IBAND,IBOTH) ! <--- gener goto 652 c C...bands c 600 WRITE(*,601,advance='NO') 601 FORMAT(/' SELECT TYPE OF BAND (-ve parameter specifies only ', * 'one of two like components):'// * T10,'0 = EXIT this option'/ * T10,'1 = type- I(-) a-R.0,1 (various K-1 for given J)'/ * T9,' 2 = type- I(+) c-R.1,0 (various K+1 for given J)'/ * T10,'3 = type-II(+) a-R.0,1 \ J decreases down ', * 'the band,',/ * T10,'4 = type-II(+) b-R.1,1 + b.R.-1,1 / (K+1 large and', * ' degenerate)'/ * T10,'5 = type-II(-) b-R.1,1 + b.R.-1,1 \ J increases down ', * 'the band,'/ * T10,'6 = type-II(-) c-R.1,0 / (K-1 large and', * ' degenerate)'/ * T10,'7 = a-Q.1,0'/ * T10,'8 = b-Q.1,-1 + b-Q.-1,1'/ * T10,'9 = c-Q.1,0 + c-Q.1,-2'/ * T9,'10 = c-R.-1,2 '/ * T40,'.... ') IBOTH=1 READ(*,14,ERR=600)IBAND IF(IBAND.EQ.0)GOTO 72 IF(IBAND.LT.0)THEN IBOTH=0 IBAND=IABS(INT(IBAND)) ENDIF IF(IBAND.LT.1.OR.IBAND.GT.11)GOTO 600 IOLD=I C CALL GENER(ICLONE,IBAND,IBOTH) ! <--- gener C 652 WRITE(*,606)I-IOLD,I 606 FORMAT(1X/I5,' lines added, total number of lines is now:',I5) GOTO 650 C C...Dump of energy levels for a specified J+1<-J transition (with 12) C C NOTE: K-1,K+1 parity: C C J.even J.odd J.even J.odd C E+ I: ee eo III: ee oe C E- eo ee oe ee C O+ oo oe eo oo C O- oe oo oo eo C C K-1 = ( Tau + J+1 ) /2 C K+1 = ( J+1 - Tau ) /2 C C (from Table 7.5 of Gordy and Cook, p.243 third ed. these parities are C for I.l and III.l and in that table entries for I.r and III.r differ C by reversal of O+ and O- lines. --To be resolved!) C 630 NEWDMP=1 636 WRITE(*,631,advance='NO') 631 FORMAT(1X/' J for J+1<-J transition: ') READ(*,14,ERR=636)J IF(J.EQ.0.or.J.eq.555.or.j-1.gt.jmax)GOTO 72 OPEN(1,FILE='DUMP.ASF',POSITION='APPEND') IF(NEWDMP.EQ.1)THEN WRITE(1,30)COM WRITE(1,31) LL=NCONSTH DO 635 L=1,LL WRITE(1,35)(K,T(K),A(K),IA(K),K=L,L+LL,LL) 635 CONTINUE IF(2*LL.LT.NCONST)WRITE(1,36)NCONST,T(NCONST),A(NCONST), * IA(NCONST) ENDIF C IEGEN=1 write(*,634)j 634 format(' ---- Now working on J=',i4) do 640 nwang=1,4 if(IABS(INT(ired)).ne.1)then call asreda(l,j,nwang,a,IRED) ! <--- asreda else call asreds(l,j,nwang,a,IRED) ! <--- asreds endif CALL HDIAG(L,IEGEN,NIT) ! <--- hdiag CALL HSORT(L) ! <--- hsort WRITE(1,632)J,wang(nwang) L=0 IF(IRED.GT.0)THEN IPLUS=5-NWANG IF(NWANG.EQ.1)IPLUS=0 ELSE IF(J-J/2*2.EQ.1)THEN IPLUS=NWANG+1 IF(NWANG.GE.3)IPLUS=4-NWANG ELSE IPLUS=NWANG-1 IF(NWANG.GE.3)IPLUS=6-NWANG ENDIF ENDIF DO 642 LL=-J+IPLUS,J,4 L=L+1 IWORK(L)=LL 642 CONTINUE WRITE(1,633)(J,(J+1+IWORK(LL))/2,(J+1-IWORK(LL))/2, * H(LL,LL),LL=1,L) 640 continue 632 FORMAT(1X/' J=',I3,', Wang matrix ',A,':') 633 FORMAT(3(I6,',',I2,',',I2,F14.3)) C J=J+1 write(*,634)j do 641 nwang=1,4 if(IABS(INT(ired)).ne.1)then call asreda(l,j,nwang,a,IRED) ! <--- asreda else call asreds(l,j,nwang,a,IRED) ! <--- asreds endif CALL HDIAG(L,IEGEN,NIT) ! <--- hdiag CALL HSORT(L) ! <--- hsort WRITE(1,632)J,wang(nwang) L=0 IF(IRED.GT.0)THEN IPLUS=5-NWANG IF(NWANG.EQ.1)IPLUS=0 ELSE IF(J-J/2*2.EQ.1)THEN IPLUS=NWANG+1 IF(NWANG.GE.3)IPLUS=4-NWANG ELSE IPLUS=NWANG-1 IF(NWANG.GE.3)IPLUS=6-NWANG ENDIF ENDIF DO 643 LL=-J+IPLUS,J,4 L=L+1 IWORK(L)=LL 643 CONTINUE WRITE(1,633)(J,(J+1+IWORK(LL))/2,(J+1-IWORK(LL))/2, * H(LL,LL),LL=1,L) 641 continue CLOSE(1) NEWDMP=0 GOTO 636 C C------------------------------------- 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 c IF(NJA(K,1)-NJA(L1,1))411,401,12 IF(NJA(K,1)-NJA(L1,1).ne.0)then ! Fortran2018 IF(NJA(K,1)-NJA(L1,1).lt.0)goto 411 ! Fortran2018 goto 12 ! Fortran2018 ENDIF ! Fortran2018 c401 IF(NJA(K,2)-NJA(L1,2))411,402,12 IF(NJA(K,2)-NJA(L1,2).ne.0)then ! Fortran2018 old label 401 IF(NJA(K,2)-NJA(L1,2).lt.0)goto 411 ! Fortran2018 goto 12 ! Fortran2018 ENDIF ! Fortran2018 c402 IF(NJA(K,3)-NJA(L1,3))411,12,12 IF(NJA(K,3)-NJA(L1,3).ge.0)goto 12 ! Fortran2018 old label 402 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 GENER(ICLONE,IBAND,IBOTH) C C GENERATION OF TRANSITIONS IN BANDS WHERE C C ICLONE =1 lines are to be cloned from supplied quantum numbers, in which C case IBAND and IBOTH do not apply C =0 lines are to be generated from band rules C IBAND = band type as below or in header in MODIF c IBOTH=1 both transitions which may become degenerate for each line to be C calculated C IBOTH=0 only one of the like transitions to be calculated C C Generated transitions will be assigned the error carried by the last C transition in the data set C C IMPLICIT REAL(8) (A-H,O-Z) c PARAMETER (NLINES=15000, NCONST=35, JMAX=250) PARAMETER (maxcbl=30000,maxcom=1000,maxspa=1000) C COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR COMMON /annot/ispace,icomm,icomml,comblk,nspces character(50) descr(nlines) DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST), * ERRORL(NLINES),WEIGHT(NLINES), * IA(NCONST),LINFIT(NLINES) CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) C C C___LINE CLONING C iold=i if(iclone.eq.0)goto 700 702 write(*,701,advance='NO') 701 format(1x/' Type in six quantum numbers for the clone template'/ * ' (comma separated)'// * ' ..... ') read(*,*,err=702)ju,kau,kcu,jl,kal,kcl if(iabs(int(ju-jl)).gt.1)goto 702 n=iabs(int(ju-(kau+kcu))) if(n.lt.0.or.n.gt.1)goto 702 n=iabs(int(jl-(kal+kcl))) if(n.lt.0.or.n.gt.1)goto 702 kaddu=kau+kcu-ju kaddl=kal+kcl-jl c 704 write(*,703,advance='NO') 703 format(1x/' Select 1 = fix J'/ * ' 2 = fix K-1'/ * ' 3 = fix K+1'// * ' (-ve value selects both possible', * ' components)'// * ' ..... ') read(*,*,err=704)nfix iboth=0 if(nfix.lt.0)then nfix=-nfix iboth=1 endif if(nfix.lt.1.or.nfix.gt.3)goto 704 c 706 write(*,705,advance='NO') 705 format(1x/' Shift in quantum number relative to this transition'/ * ' (comma separated)'// * ' ..... ') read(*,*,err=706)minqn,maxqn if(minqn.gt.maxqn)goto 706 c if(iboth.eq.1)then if(nfix.eq.1)then 716 write(*,715,advance='NO') 715 format(1x/' Select -1 = prolate (K-1) doubling'/ * ' 1 = oblate (K+1) doubling'// * ' ..... ') read(*,*,err=716)kdubl if(iabs(int(kdubl)).ne.1)goto 716 else if(nfix.eq.2)kdubl=-1 if(nfix.eq.3)kdubl=1 endif endif c C...Generation loop C do 707 k=minqn,maxqn c c...fixed J c if(nfix.eq.1)then jjl=jl kkl=kal+k if(kkl.gt.jjl)goto 707 kkkl=jjl-kkl+kaddl if(kkkl.gt.jjl)goto 707 jju=ju kku=kau+k if(kku.gt.jju)goto 707 kkku=jju-kku+kaddu if(kkku.gt.jju)goto 707 endif c c...fixed K-1 c if(nfix.eq.2)then jju=ju+k if(jju.lt.0)goto 707 jjl=jl+k if(jjl.lt.0)goto 707 kku=kau kkl=kal kkku=jju-kku+kaddu kkkl=jjl-kkl+kaddl endif c c...fixed K+1 c if(nfix.eq.3)then jju=ju+k if(jju.lt.0)goto 707 jjl=jl+k if(jjl.lt.0)goto 707 kkku=kcu kkkl=kcl kku=jju-kkku+kaddu kkl=jjl-kkkl+kaddl endif c if(kkl .lt.0.or.kku .lt.0)goto 707 if(kkkl.lt.0.or.kkku.lt.0)goto 707 i=i+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(I-1) WEIGHT(I)=1.D0/ERRORL(I)**2 DESCR(I)=' ' LINFIT(I)=0 c JK(I,1)=jju JK(I,4)=jjl JK(I,2)=kku JK(I,5)=kkl JK(I,3)=kkku JK(I,6)=kkkl c if(iboth.eq.1)then if(kdubl.eq.-1.and.( kkl.eq.0.or. kku.eq.0))goto 707 if(kdubl.eq. 1.and.(kkkl.eq.0.or.kkku.eq.0))goto 707 c if(kdubl.eq.-1)then if(kku+kkku.eq.jju)then kkku=kkku+1 else kkku=kkku-1 endif if(kkl+kkkl.eq.jjl)then kkkl=kkkl+1 else kkkl=kkkl-1 endif else if(kku+kkku.eq.jju)then kku=kku+1 else kku=kku-1 endif if(kkl+kkkl.eq.jjl)then kkl=kkl+1 else kkl=kkl-1 endif endif c i=i+1 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(I-1) WEIGHT(I)=1.D0/ERRORL(I)**2 DESCR(I)=' ' LINFIT(I)=0 JK(I,1)=jju JK(I,4)=jjl JK(I,2)=kku JK(I,5)=kkl JK(I,3)=kkku JK(I,6)=kkkl endif c 707 continue return C C C___GENERATION ACCORDING TO BAND RULES C C...a-R.0,1, type-I(-): C 700 IERROR=I IF(IBAND.EQ.1)THEN 603 WRITE(*,610,advance='NO')'Type-I(-) a-R.0,1: J' 610 FORMAT(1X/1X,8(1H-),T10,A,' = ') READ(*,14,ERR=603)JBAND 14 FORMAT(I5) IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 603 604 WRITE(*,602,advance='NO')' K-1 range (lower,upper)' 602 FORMAT(1X/T10,A,' = ') READ(*,*,ERR=604)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 604 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 604 IF(KMHIGH.GT.JBAND)GOTO 604 DO 1605 K=KMLOW,KMHIGH ! Fortran2018 DO 605 KK=0,1 IF((K.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 605 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=JBAND+1 JK(I,4)=JBAND JK(I,2)=K JK(I,5)=K JK(I,3)=JK(I,1)-JK(I,2)+KK JK(I,6)=JK(I,4)-JK(I,5)+KK FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 DESCR(I)=' ' LINFIT(I)=0 605 CONTINUE 1605 CONTINUE ! Fortran2018 ENDIF C C...c-R.1,0, type-I(+): C IF(IBAND.EQ.2)THEN 653 WRITE(*,610)'Type-I(+) c-R.1,0: J range (lower,upper)' READ(*,*,ERR=653)JBLOW,JBHIGH IF(JBHIGH.LT.JBLOW)GOTO 653 IF(JBLOW.LT.0)GOTO 653 IF(JBHIGH.GT.JMAX)GOTO 653 654 WRITE(*,602)' K+1 range (lower,upper)' READ(*,*,ERR=654)KPLOW,KPHIGH IF(KPHIGH.LT.KPLOW)GOTO 654 IF(KPLOW.LT.0.OR.KPLOW.GT.JBHIGH)GOTO 654 IF(KPHIGH.GT.JBHIGH)GOTO 654 DO 1655 JBAND=JBLOW,JBHIGH ! Fortran2018 DO 655 K=KPLOW,KPHIGH IF(K.GT.JBAND)GOTO 655 DO 656 KK=0,1 IF((K.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 656 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=JBAND+1 JK(I,4)=JBAND JK(I,3)=K JK(I,6)=K JK(I,2)=JK(I,1)-JK(I,3)+KK JK(I,5)=JK(I,4)-JK(I,6)+KK FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' 656 CONTINUE 655 CONTINUE 1655 CONTINUE ! Fortran2018 ENDIF C C...a-R.0,1 and b-R, type-II(+): bands are formed when (A+B)=(l/m)(2C), where C l and m are integers such that (l/m)>1 C IF(IBAND.EQ.3.OR.IBAND.EQ.4)THEN 616 WRITE(*,610) 'Type-II(+): l and m defining n=(l/m)' read(*,*,err=616)LBAND,MBAND if(lband.lt.2.or.mband.lt.1)goto 616 if(lband.le.mband)goto 616 IF(IBAND.EQ.3)THEN 611 WRITE(*,610)'Type-II(+) a-R.0,1: J of leading line' READ(*,14,ERR=611)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 611 ELSE 615 WRITE(*,610)'Type-II(+) b-R: J of leading line' READ(*,14,ERR=615)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 615 ENDIF isband=0 if(mband.gt.1)then write(*,'(1x/'' Subbands are possible for following K-1: '', * 15i3)')(kk,kk=0,mband-1) 630 write(*,602)' K-1 for leading line in subband' read(*,14,err=630)isband if(isband.lt.0.or.isband.gt.mband-1)goto 630 endif c 612 WRITE(*,602) ' K-1 range (lower,upper)' READ(*,*,ERR=612)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 612 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 612 IF(KMHIGH.GT.JBAND)GOTO 612 jline=jband-(mband-lband) DO 6614 K=ISBAND,KMHIGH,MBAND jline=jline+(mband-lband) if(k.lt.kmlow)goto 6614 if(jline-k.lt.0)goto 6614 DO 614 KK=0,1 IF(IBOTH.EQ.0.AND.KK.EQ.1)GOTO 6614 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,4)=jline JK(I,1)=JK(I,4)+1 IF(IBAND.EQ.3)THEN ! aR 0,1 JK(I,2)=K+KK JK(I,5)=K+KK JK(I,6)=JLINE-K JK(I,3)=JK(I,6)+1 ELSE ! bR 1,1 JK(I,5)=K+KK ! bR 1,-1 JK(I,2)=K+1-KK JK(I,6)=JLINE-K JK(I,3)=JK(I,6)+1 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' 614 CONTINUE 6614 CONTINUE ENDIF C C...b-R and c.R.1,0, type-II(-): bands are formed when 2A=(l/m)(B+C), where C l and m are integers such that (l/m)>1 C IF(IBAND.EQ.5.OR.IBAND.EQ.6)THEN 516 WRITE(*,610) 'Type-II(-): l and m defining n=(l/m)' read(*,*,err=516)LBAND,MBAND if(lband.lt.2.or.mband.lt.1)goto 516 if(lband.le.mband)goto 516 c IF(IBAND.EQ.5)THEN 511 WRITE(*,610)'Type-II(-) b-R: J of leading line' READ(*,14,ERR=511)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 511 ELSE 515 WRITE(*,610)'Type-II(-) c-R.1,0: J of leading line' READ(*,14,ERR=515)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 515 ENDIF c write(*,'(1x/'' Subbands are possible for following K+1: '', * 15i3)')(kk,kk=0,lband-1) 530 write(*,602)' K-1 for leading line in subband' read(*,14,err=530)isband if(isband.lt.0.or.isband.gt.lband-1)goto 530 c 512 WRITE(*,602) ' K+1 range (lower,upper)' READ(*,*,ERR=512)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 512 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 512 IF(KMHIGH.GT.JBAND)GOTO 512 jline=jband-(lband-mband) DO 5514 K=isband,KMHIGH,lband jline=jline+(lband-mband) if(k.lt.kmlow)goto 5514 if(jline-k.lt.0)goto 5514 DO 514 KK=0,1 IF(IBOTH.EQ.0.AND.KK.EQ.1)GOTO 514 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,4)=jline JK(I,1)=JK(I,4)+1 IF(IBAND.EQ.6)THEN ! cR 1,0 JK(I,6)=K+KK JK(I,3)=K+KK JK(I,5)=JLINE-K JK(I,2)=JK(I,5)+1 ELSE ! bR 1,1 JK(I,6)=K+KK ! bR 1,-1 JK(I,3)=K+1-KK JK(I,5)=JLINE-K JK(I,2)=JK(I,5)+1 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' 514 CONTINUE 5514 CONTINUE ENDIF C C...a-Q.1,0 C IF(IBAND.EQ.7)THEN WRITE(*,'(1X/'' ***** SORRY, NOT YET IMPLEMENTED''/)') ENDIF C C...b-Q.1,-1 + b-Q.-1,1 C IF(IBAND.EQ.8)THEN 607 WRITE(*,610)'b-Q.1,-1 (b-Q.-1,1): K-1' READ(*,14,ERR=607)KBAND IF(KBAND.LT.0.OR.KBAND.GT.JMAX)GOTO 607 608 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=608)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 608 IF(JLOW.LT.KBAND)GOTO 608 IF(JHIGH.GT.JMAX)GOTO 608 DO 1609 J=JLOW,JHIGH ! Fortran2018 DO 609 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 609 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J JK(I,4)=J JK(I,2)=KBAND+1 JK(I,5)=KBAND JK(I,6)=JK(I,4)-JK(I,5)+KK JK(I,3)=JK(I,6)-1 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' 609 CONTINUE 1609 CONTINUE ! Fortran2018 ENDIF C C...c-Q.1,0 + c-Q.1,-2 C IF(IBAND.EQ.9)THEN 670 WRITE(*,610)'c-Q.1,0 (c-Q.1,-2): K-1' READ(*,14,ERR=670)KBAND IF(KBAND.LT.0.OR.KBAND.GT.JMAX)GOTO 670 671 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=671)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 671 IF(JLOW.LT.KBAND.OR.JLOW.EQ.0)GOTO 671 IF(JHIGH.GT.JMAX)GOTO 671 DO 1672 J=JLOW,JHIGH ! Fortran2018 DO 672 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 672 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J JK(I,4)=J JK(I,2)=KBAND+1 JK(I,5)=KBAND IF(KK.EQ.0)THEN JK(I,6)=JK(I,4)-JK(I,5) JK(I,3)=JK(I,6) ELSE JK(I,6)=JK(I,4)-JK(I,5)+1 JK(I,3)=JK(I,6)-2 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' IF(JK(I,2)+JK(I,3).GT.JK(I,1)+1)I=I-1 IF(JK(I,2).LT.0)I=I-1 672 CONTINUE 1672 CONTINUE ! Fortran2018 ENDIF C C...c-R.-1,2 C IF(IBAND.EQ.10)THEN 650 WRITE(*,610)'c-R.-1,2: K-1' READ(*,14,ERR=650)KBAND IF(KBAND.LT.1.OR.KBAND.GT.JMAX-2)GOTO 650 651 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=651)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 651 IF(JLOW.LT.KBAND)GOTO 651 IF(JHIGH.GT.JMAX)GOTO 651 DO 1652 J=JLOW,JHIGH ! Fortran2018 DO 652 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 652 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J+1 JK(I,4)=J JK(I,2)=KBAND+KK-1 JK(I,5)=KBAND+KK JK(I,6)=JK(I,4)-JK(I,5)+KK JK(I,3)=JK(I,6)+2 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERRORL(I)=ERRORL(IERROR) WEIGHT(I)=1.D0/ERRORL(I)**2 LINFIT(I)=0 DESCR(I)=' ' IF(JK(I,2)+JK(I,3).GT.JK(I,1)+1)I=I-1 IF(JK(I,2).LT.0)I=I-1 652 CONTINUE 1652 CONTINUE ! Fortran2018 ENDIF C RETURN END C C_____________________________________________________________________________ C SUBROUTINE DATOUT(SDEV,SWDEV,NINFIT) C C OUTPUT OF DATA AND RESULTS FROM ASFIT - BOTH TO SCREEN AND TO DISK C C SDEV - standard deviation of the fit C SWDEV - deviation per unit weight C NINFIT - number of different frequency transitions in fit C NCONST - total no of fittable constants C also: C IOUT - if set to 1 on input, then only disk output of data is carried out C if set to 2, then output menu is provided for selection of C disk output of data and/or intermediate results C K - output mode which is determined from the menu at the bottom of C the routine C =0 no output C =1 standard output (including c.c., worst and most infl.transitions) C =2 above + c.c. distribution plot C =3 above + contributions or hjj/DFBETAS C =4 constants file for ASROT C =5 frequency listing in deposition format C =6 .LIN file for SPFIT C IMPLICIT REAL(8) (A-H,O-Z) INTEGER(4) NERD,NDCON,NDEROR,iend c PARAMETER (NLINES=15000,JMAX=250, MAXDIM=JMAX/2+1, * lpersc=25, nerd=2, NCONST=35, NCONSTH=17) PARAMETER (maxcbl=30000,maxcom=1000,maxspa=1000) C COMMON /BLKH/H * /BLKR/R * /SORTCC/WORK,IPT COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR * /BLK3/C,D,DD,HJJ * /TBLOCK/COM,T,PINUMS INTEGER(4) IA(NCONST),LINFIT(NLINES),JK(NLINES,6),IPT(NLINES), * TPT(NCONST),UNITS(NCONST),NCC(20), * namcd(10),tptinv(nconst) REAL(8) FR(NLINES,4),A(NCONST),C(NCONST),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),DD(NCONST,NCONST),WORK(NLINES), * DFBETA(NCONST),PICON(NCONST),ERRORL(NLINES), * WEIGHT(NLINES),D(NLINES,NCONST),HJJ(NLINES) c REAL(4) DIV character(50) descr(nlines) character(8) tdescr,ctemp CHARACTER T(NCONST)*6,linout*500,lincom*100,linrel*10,lannot*100 CHARACTER(25) FILNAM CHARACTER(27) conval,erval CHARACTER(40) OUTCON(NCONST) CHARACTER(7) UNIT,NAMUN(6) CHARACTER COM(72),FORMAR(54),NAMRED,namrep*5,cdum1*20,cdum2*20 CHARACTER(54) FORMAL,F1,F2 EQUIVALENCE (FORMAL,FORMAR(1)) CHARACTER comblk*(maxcbl),csymb,pinums(nconst)*3,vv*2 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces 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, 3,3,3,3,3,3,3, 4,4,4,4,4,4,4,4,4, * 5,5,5,5,5,5,5,5,5,5,5/ c iexdep=0 NOUT=0 C C...F1 - format for first line of printout of contributions, F2 - format C for all following lines for the same transition F1='( I4,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)' namrep='I.r ' NAMRED='A' IF(iabs(int(IRED)).EQ.1)NAMRED='S' IF(IRED.LT.0)NAMREP='III.r' IF(IOUT.EQ.1)GOTO 100 C K1=0 DO 9 NN=1,NCONST IF(IA(NN).NE.1)GOTO 9 K1=K1+1 TPT(K1)=NN TPTINV(nn)=k1 9 CONTINUE C C...go to the output menu C IF(IOUT.EQ.2)GOTO 300 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...Output of results (standard listing) which can be routed to screen C (NOUT=0) or to disk (NOUT=3) C 25 if(k.ne.5)then WRITE(NOUT,1)COM,NAMRED,NAMREP else write(nout,61)com,namred,namrep endif 1 FORMAT( T54,'ASFIT version 22.I.2024'/ * 79(1H_)//1x,72A1/79(1H_)//' Fit of Watson''s ', * 'Hamiltonian --> reduction-',A,', representation-',A,':'//) 61 FORMAT(79(1H_)//1X,72A1/79(1H_)//' Obs-calc differences', * ' are from fit of Watson''s reduced Hamiltonian in'/33x, * '---> reduction-',A,', representation-',A,'<---'/) L=0 IF(K.NE.5)THEN 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) WRITE(NOUT,'(1x)') if(nout.eq.0)then 902 write(*,901,advance='NO') 901 format(' F I T C O M P L E T E D: ', 1 '0 = continue'/ 1 29x,' 1 = scroll through fitted lines'/ 1 29x,'-1 = as above, screen at a time .... ') read(*,'(I2)',err=902)nscrol if(iabs(int(nscrol)).ne.1)goto 300 WRITE(NOUT,'(1x)') WRITE(NOUT,53) WRITE(NOUT,6) 6 FORMAT(17X,'TRANSITION',T42,'Observed',T56,'Obs-Calc',T68, 1 'Error') WRITE(NOUT,53) else WRITE(NOUT,553) WRITE(NOUT,666) 666 FORMAT(17X,'TRANSITION',T42,'Observed',T56,'Obs-Calc',T68, 1 'Error Note Blend components') WRITE(NOUT,553) endif ELSE WRITE(NOUT,67) 67 FORMAT(/10X,'TRANSITION',T35,'Observed',T47,'Obs-Calc',T59, 1 'Error Note wt for'/T76,'blend'/) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...TRANSITIONS, listed for all k but: C C if k=5 then deposition type printout C if NOUT=0 and nscrol=-1 then scrolled page at a time to screen c nspace=1 ncomm=2 if(iabs(int(nscrol)).eq.1)then lmax=nlines+maxcom+maxspa IF(nscrol.LT.0)lmax=lpersc-1 lcount=3 endif DO 7 N=1,I lendes=len_trim(descr(n)) if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') lcount=lcount+1 if(nout.eq.0.and.lcount.eq.lmax)then write(*,716,advance='NO') read(*,'(i1)',err=707)jj if(jj.eq.1)goto 709 707 lcount=0 endif endif 716 format( * ' Press to continue, 1 to stop .... ') c 33 if(icomml(ncomm).eq.n)then if(nout.eq.0)then write(NOUT,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) else write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) endif ncomm=ncomm+1 lcount=lcount+1 if(nout.eq.0.and.lcount.eq.lmax)then write(*,716,advance='NO') read(*,'(i1)',err=706)jj if(jj.eq.1)goto 709 706 lcount=0 goto 33 endif goto 33 endif c np=n if(n/10000.ge.1)np=n-(n/10000)*10000 c c...transition in fit c if(len_trim(descr(n)).le.1.and.descr(n)(1:1).eq.' ')then lannot=' ' else lannot=' /'//descr(n) endif c IF(LINFIT(N).gt.0)then c c...non-deposition if(k.ne.5)then if(abs(fr(n,3)).le.3.d0*errorl(n))then if(nout.ne.0)then ! disk output if(linfit(n).eq.1) * WRITE(NOUT,81)NP,' ',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERRORL(N),lannot(1:len_trim(lannot)) if(linfit(n).eq.2) * WRITE(NOUT,181)NP,'B',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERRORL(N),lannot(1:len_trim(lannot)),fr(n,4),FR(n,2) else ! screen output WRITE(NOUT,721)N,' ',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERRORL(N) endif else csymb='?' if(nout.ne.0)then ! disk output if(linfit(n).eq.1) * WRITE(NOUT,81)NP,csymb,(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERRORL(N),' >3*Err ', * lannot(1:len_trim(lannot)) if(linfit(n).eq.2) * WRITE(NOUT,181)NP,csymb,(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERRORL(N),' >3*Err ',fr(n,4),FR(n,2) else ! screen output tdescr='>3*Err' if(abs(fr(n,3)).le.3.d0*errorl(n))tdescr='<=Blnd' WRITE(NOUT,721)N,csymb,(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERRORL(N),tdescr(1:6) endif endif c c...deposition else ndig=10 if(errorl(n).gt.0.01d0)then write(cdum1,'(f10.3)')fr(n,3) write(cdum2,'(f10.3)')errorl(n) else write(cdum1,'(f10.4)')fr(n,3) write(cdum2,'(f10.4)')errorl(n) endif call lzero(cdum1,ndig) ! <--- lzero call lzero(cdum2,ndig) ! <--- lzero if(fr(n,4).eq.0.d0)then if(errorl(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),lannot(1:len_trim(lannot)) else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),lannot(1:len_trim(lannot)) endif else if(errorl(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),lannot(1:len_trim(lannot)),fr(n,4) else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),lannot(1:len_trim(lannot)),fr(n,4) endif endif endif c c...transition excluded from fit c c...non-deposition else if(k.ne.5)then if(nout.ne.0)then WRITE(NOUT,81)NP,'E',(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERRORL(n),'--excl ',lannot(1:len_trim(lannot)) else WRITE(NOUT,721)N,'E',(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERRORL(n),'--excl' endif else c c...deposition (optional) if(iexdep.eq.1)then ndig=10 if(errorl(n).gt.0.01d0)then write(cdum1,'(f10.3)')fr(n,3) write(cdum2,'(f10.3)')errorl(n) else write(cdum1,'(f10.4)')fr(n,3) write(cdum2,'(f10.4)')errorl(n) endif call lzero(cdum1,ndig) ! <--- lzero call lzero(cdum2,ndig) ! <--- lzero if(fr(n,4).eq.0.d0)then if(errorl(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'---excl ' else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'--excl ' endif else if(errorl(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'---excl ',fr(n,4) else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'--excl ',fr(n,4) endif endif endif c endif endif lcount=lcount+1 if(NOUT.eq.0.and.lcount.eq.lmax)goto 700 goto 7 c 700 write(*,716,advance='NO') read(*,'(i1)',err=705)jj if(jj.eq.1)goto 709 705 lcount=0 c 7 CONTINUE if(NOUT.eq.0.and.lcount.gt.lmax-12)then write(*,'('' Press ENTER to continue '')',advance='NO') read(*,'(i1)',err=709)jj endif c 709 if(k.ne.5)then if(nout.eq.0)then WRITE(NOUT,53) else WRITE(NOUT,553) endif endif c 81 FORMAT(I4,'.',A,2(I5,2I4,1x),F15.4,F14.4,F9.4,2a) 181 FORMAT(I4,'.',A,2(I5,2I4,1x),F15.4,' B',F12.4,F9.4,a, * 1PE8.1,0PF14.4) 721 FORMAT(I5,'.',A,2(I5,2I4,1x),F15.4,F14.4,F9.4,2a) 68 FORMAT(1X,2(i3,','),i3,' <-',2(I3,','),i3,F17.3,3a,1PE10.2) 168 FORMAT(1X,2(i3,','),i3,' <-',2(I3,','),i3,F18.4,3a,1PE9.2) IF(NOUT.NE.0.and.k.ne.5)WRITE(NOUT,410) * SDEV,NINFIT,SWDEV,NINFIT-K1 410 FORMAT(/' Sigma = ',F14.6,27X,I5,' Distinct lines in fit'/ 1 ' Sigma.w = ',F14.6,27X,I5,' Degrees of freedom') C C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...CONSTANTS C C...On output to disk additional conversion of constants 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.and.k.ne.5)THEN WRITE(NOUT,'(1X)') L=0 DO 1002 N=1,NCONST 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 (rounding as coded tends in rare cases to go into c an infinite loop, hence the NITER cludge) c niter=0 L=L+1 ESCLD=C(L)*10.D0**(3*(UNITS(N)-1)) 417 WRITE(CONVAL,'(F27.16)')ASCLD WRITE(ERVAL,'(F27.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR,NERD) ! <--- confor c c...rounding of the parameter value, if necessary if(conval(ndcon+1:ndcon+1).eq.'5'.or. * conval(ndcon+1:ndcon+1).eq.'6'.or. * conval(ndcon+1:ndcon+1).eq.'7'.or. * conval(ndcon+1:ndcon+1).eq.'8'.or. * conval(ndcon+1:ndcon+1).eq.'9')then c c write(*,'(1x,2a,i8,f27.16)')conval,conval(1:ndcon), DEBUG c * ndcon,10.d0**(-(ndcon-10)) DEBUG ascld=ascld+ * dsign(1.d0,ascld)*10.d0**(-(ndcon-10)) niter=niter+1 if(niter.le.5)goto 417 endif c write(*,'(1x,a,14x,a,i8,f27.16)')'Final_fitted:', DEBUG c * conval(1:ndcon),ndcon,10.d0**(-(ndcon-10)) DEBUG c c 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 c 1003 call CONFIX(a(n),ascld,conval,iend) ! <----- WRITE(OUTCON(N),1005)T(N)(1:5),UNIT,CONVAL(1:iend) c 1002 CONTINUE c LLLL=NCONSTH 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,A40,1X,A40) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C IF(NOUT.EQ.0)GOTO 24 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C Variance covariance matrix (not printed) and correlation coefficients C (lower triangle only). Also mean values of correlation coefficients. C IF(K.LT.1.or.k.eq.5)GOTO 31 SDEVSQ=SWDEV*SWDEV DO 200 N=1,K1 DO 1200 L=N,K1 ! Fortran2018 H(N,L)=DD(N,L)*SDEVSQ 1200 CONTINUE ! Fortran2018 200 CONTINUE DO 201 N=1,K1 DO 1201 L=1,N ! Fortran2018 DENOM=C(L)*C(N) IF(DENOM.EQ.0.D0)DENOM=1.D-50 H(N,L)=H(L,N)/DENOM 1201 CONTINUE ! Fortran2018 201 CONTINUE C WRITE(NOUT,202) 202 FORMAT(/' CORRELATION COEFFICIENTS, C.ij:') DO 3035 N=1,20 NCC(N)=0 3035 CONTINUE NCORCO=0 CORCAV=0.D0 CCSUM=0.0D0 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) DO 3030 L=M,LLLL IF(L.NE.N)THEN NCORCO=NCORCO+1 CORCAV=CORCAV+DABS(H(N,L)) CCSUM=CCSUM+H(N,L) NNN=int( DABS(H(N,L))/0.05D0+1.D0 ) IF(NNN.EQ.21)NNN=20 NCC(NNN)=NCC(NNN)+1 IF(NCC(NNN).GT.67)NCC(NNN)=67 ENDIF 3030 CONTINUE 302 CONTINUE 303 FORMAT(1X,A5,1X,8(F9.4)) IF(MM.GE.K1)THEN GOTO 3031 ELSE M=M+8 GOTO 305 ENDIF 3031 IF(NCORCO.NE.0)THEN WRITE(NOUT,3032)CORCAV/NCORCO,CCSUM/NCORCO 3032 FORMAT(1X/' Mean value of |C.ij|, i.ne.j =', * F9.4/ * ' Mean value of C.ij, i.ne.j =', * F9.4/) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C...bar graph of distribution of correlation coefficients (disk only) IF(K.eq.2.or.k.eq.3)THEN WRITE(NOUT,'(1X/'' Distribution of values of |C.ij|, i.ne.j:'' * //10X,6(8X,I2)/10X,6(''....,....|''),''....,..'' * )')(N,N=10,60,10) DO 3036 N=1,20 write(ctemp,'(f8.2)')N*0.05 ! f4.2 simulation WRITE(NOUT,3037)NCC(N),ctemp(5:8),('*',LLLL=1,NCC(N)) 3036 CONTINUE 3037 FORMAT(1X,I2,' <',A,1X,67A1) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...MOST INFLUENTIAL TRANSITIONS and WORST FITTING LINES C C Most infl. trans. are defined by magnitude of diagonal elements hjj of C H = J (J'J)**(-1) J' C C here the design (derivative) matrix J is D, and (J'J)**(-1) is DD, C in ZK thesis J is X C IF(NOUT.EQ.3.and.k.ne.5)THEN call hjjset(K1,hjjav) ! <--- hjjset c hjjcut=0.5 write(nout,'(79(1h_))') IF(work(i).lt.hjjcut)then write(nout,609)hjjcut goto 608 else write(nout,610)hjjcut endif 609 format(1x/' There are no influential transitions with hjj >', * f6.3) 610 FORMAT(1x/' MOST INFLUENTIAL TRANSITIONS:'/1x, * T43,'Observed',T56,'Obs-Calc',T67,'hjj >',F6.3/) do 607 L=I,1,-1 LL=ipt(L) if(linfit(LL).eq.0)goto 607 if(work(L).lt.hjjcut)goto 608 write(nout,606)LL,(JK(LL,LLL),LLL=1,6),FR(LL,1),FR(LL,3), * hjj(ll) 607 continue 606 FORMAT(I5,'. ',2(I5,2I4,1x),F15.4,F14.4,F9.4) 608 write(nout,611)hjjav 611 format(1x/1x,T47,'mean of all values =',f7.4) c c...worst fitting lines write(NOUT,105) do 106 L=1,I if(LINFIT(L).gt.0)then work(L)=abs(fr(L,3)/errorl(L)) if(L.ne.1)then if( (linfit(L).eq.2.and.linfit(L-1).eq.2) .and. * (fr(L,1).eq.fr(L-1,1)) )work(L)=0.d0 endif else WORK(L)=0.d0 endif ipt(L)=L 106 continue int21=1 call sortc(int21,I) ! <--- sortc c if(i.le.5)then nhj=1 else nhj=i-5 endif c write(NOUT,110)(ipt(L),FR(ipt(L),3)/errorl(ipt(L)),L=I,nhj,-1) 105 FORMAT(1x/' WORST FITTING LINES (line number=(o-c)/error):') 110 FORMAT(1x,6(I7,'=',F5.1)) c ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C Contributions to frequencies (disk only) or hjj/DFBETAS values for C transitions C IF(K-2.LT.1)GOTO 31 ! old label 304 C C...heading and options C 501 write(*,500,advance='NO') 500 format(1x/ * 8x,'-1 = values of hjj/DFBETAS for transitions'/ * 8x,' 0 = All contributions to transition frequencies'/ * 8x,' 1 = Only contributions from selected constants'// * 8x,'.... ') read(*,'(i5)',err=501)iconse c c...hjj/DFBETAS - for exposition see: C C J.Demaison, J.Cosleou, R.Bocquet, and A.G.Lesarri, C J.Mol.Spectrosc. 167,400-418(1994) c if(iconse.eq.-1)then c c..Transitions and HJJ's c call hjjset(K1,hjjav) ! <--- hjjset write(nout,'(79(1h_))') WRITE(NOUT,350) nspace=1 ncomm=2 write(*,'(1x/7x,'' '')',advance='NO') do 351 n=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 355 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 355 endif C if(linfit(N).eq.0)goto 351 write(NOUT,352)n,(jk(n,l),l=1,6),fr(n,1),fr(n,3), * errorl(n),fr(n,3)/errorl(n),hjj(n) 351 continue 350 FORMAT(//' FITTED TRANSITIONS, SCALED DEVIATIONS, HJJ''S:'// * 1x,t35,'obs',t46,'obs-calc',t57,'error ', * '(o-c)/Err',t77,'hjj'/) 352 FORMAT(i5,'.',2x,3i3,i4,2I3,F14.4,F12.4,f8.4,F9.3,F10.5) c C...HJJ's and DFBETAS C C The mode of DFBETAS calculation used here i.e. repetition of the whole C least squares fit for each line in turn is probably the slowest possible, C but it is a rare operation and is fast enough anyway - the Demaison C paper cites references, which probably contain a much faster route. C nspace=1 ncomm=2 write(nout,370) write(nout,376)(T(TPT(L))(1:5),L=1,K1) do 371 n=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 375 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 375 endif if(n.eq.(n/30)*30)write(nout,376)(T(TPT(L))(1:5),L=1,K1) c if(linfit(N).eq.0)goto 371 if(n.gt.1.and.(fr(n-1,1).eq.fr(n,1)) )goto 371 WRITE(*,802,advance='NO') ! counter 802 FORMAT('.') if((n/50)*50.eq.n)write(*,'(1x)') C CALL DFBET(DFBETA,K1,N) ! <--- dfbet C KFIN=K1 IF(KFIN.GT.11)KFIN=11 write(NOUT,373)n,(jk(n,l),l=1,6),fr(n,1), * hjj(n),(dfbeta(L),L=1,KFIN) 377 IF(KFIN.LT.K1)THEN KST=KFIN+1 KFIN=KST+10 IF(KFIN.GT.K1)KFIN=K1 write(NOUT,374)(dfbeta(L),L=KST,KFIN) GOTO 377 ENDIF 371 continue write(*,'(1x)') goto 31 endif 370 FORMAT(//' HJJ''S AND DFBETAS:') 376 FORMAT(1x/46x,'hjj ',11(1x,a)/52x,11(1x,a)/52x,11(1x,a)/ * 52x,11(1x,a)/) 373 FORMAT(i5,'.',2x,3i3,i4,2I3,F14.4,F10.5,11F6.2) 374 FORMAT(51x,11F6.2) c c...continuation of selection for contributions c if(iconse.lt.0.or.iconse.gt.1)goto 501 if(iconse.eq.1)then 502 write(*,503,advance='NO') 503 format(1x/' Please specify up to 10 comma separated numbers', * ' of constants for contributions'/' (smaller number of constants' * ,' is to be padded out with commas)'/10x,' ') do 515 n=1,10 namcd(n)=0 515 continue read(*,*,err=502)namcd noutc=0 do 504 n=1,10 if(namcd(n).ne.0)then noutc=noutc+1 if(namcd(n).lt.0.or.namcd(n).gt.nconst. * or.ia(namcd(n)).ne.1)namcd(n)=tpt(1) else goto 505 endif 504 continue 505 write(*,506)(t(namcd(n))(1:5),n=1,noutc) 506 format(1x/' Derivatives will be output for following', * ' constants:'/11x,10a) 508 write(*,507,advance='NO') 507 format(30x,' OK (1/0) ? ') read(*,'(i5)',err=508)iok if(iok.ne.1)goto 502 endif C C...CONTRIBUTIONS - the F descriptor in output format C is chosen individually for each derivative according to its magnitude C C...values of contributions DO 203 N=1,I DO 1203 L=1,K1 ! Fortran2018 D(N,L)=D(N,L)*A(TPT(L)) 1203 CONTINUE ! Fortran2018 203 CONTINUE c write(nout,'(79(1h_))') WRITE(NOUT,204) 204 FORMAT(/' CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO FREQUENCY:'/) if(iconse.eq.1)then write(nout,509)(t(namcd(n))(1:5),n=1,noutc) 509 format(1x,t33,'obs',t46,'o-c IFIT',2x,10(3x,a,3x)) write(nout,'(1x)') do 511 n=1,i write(NOUT,510)n,(jk(n,l),l=1,6),fr(n,1),fr(n,3),linfit(n), * (d(n,tptinv(namcd(l))),l=1,noutc) 511 continue goto 514 endif 510 FORMAT(i5,'.',2x,3i3,i4,2I3,2F12.3,i2,' ',10F11.3) c 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) nspace=1 ncomm=2 DO 320 N=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 380 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 380 endif C FORMAL=F1 IST=1 315 ISP=IST+7 IF(ISP.GT.K1)ISP=K1 DO 310 L=IST,ISP ! old label 314 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...restore original derivatives C 514 DO 512 N=1,I DO 1512 L=1,K1 ! Fortran2018 D(N,L)=D(N,L)/A(TPT(L)) 1512 CONTINUE ! Fortran2018 512 CONTINUE c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...end of results output c 31 if(k.eq.5)then write(nout,'(79(1H_)/79(1H_))') else if(k.le.3)write(nout,'(79(1h_))') endif NOUT=0 CLOSE(3) GOTO 23 c WRITE(NOUT,405)SDEV,SWDEV ! old label 30 24 IF(IOUT.NE.1.AND.IOUT.NE.2)GOTO 300 IOUT=0 GOTO 16 C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C O U T P U T M E N U C 300 WRITE(*,22,advance='NO') 22 FORMAT(/' OUTPUT OF RESULTS - CHOOSE ONE OF:'// 1 9X,'0 = no output'/ 1 9X,'1 = standard output'/ 1 9X,'2 = as above + distribution plot of corr. coeff.'/ 1 9X,'3 = as above + freq. contributions or hjj/DFBETAS'/ 1 9X,'4 = constants file for ASROT'/ 1 9X,'5 = frequency listing in deposition format'/5x, 1'+6 -6 = generate .LIN and .PAR files for SPFIT'/ 1 9x,' (+6 for asymmetric,-6 for symmetric quantisation)'//5x, 1' .... ') READ(*,15,ERR=300)KINDX K=iabs(int(KINDX)) IF(K.LT.1.OR.K.GT.6)GOTO 23 NOUT=3 1901 if(K.ne.6)then WRITE(*,17,advance='NO') READ(*,'(A)',ERR=1901)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) endif if(k.eq.5)then 401 write(*,400,advance='NO') 400 FORMAT(/' ARE EXCLUDED LINES ALSO TO BE LISTED (1/0): ') read(*,'(I5)',err=401)iexdep if(iexdep.lt.0.or.iexdep.gt.1)goto 401 endif C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...output of a constants file (with 4) C IF(K.EQ.4)THEN DO 6869 LL=NCONST,1,-1 IF(A(LL).NE.0)GOTO 6870 6869 CONTINUE LL=NCONST 6870 WRITE(3,6866)LL LLL=0 DO 6867 L=1,LL ERCON=0.D0 IF(IA(L).EQ.1)THEN LLL=LLL+1 ERCON=C(LLL) ENDIF WRITE(3,6868)T(L),A(L),ERCON 6867 CONTINUE 6866 FORMAT(' NCON =',I5) 6868 FORMAT(1X,A6,2F30.19) WRITE(3,'('' |''/3x,22(1H-),4A)')' Reduction-',NAMRED, * ', Representation-',NAMREP GOTO 31 ENDIF C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...output of a .LIN file for SPFIT (with 6) asymmetric top quantisation c (with -6) symmetric top quantisation C IF(K.EQ.6)THEN 1401 write(*,1400,advance='NO') 1400 FORMAT(/' Generic name for .LIN and .PAR files: ') READ(*,'(A)',ERR=1401)FILNAM filnam=filnam(1:len_trim(filnam))//'.lin' OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1401) nspace=1 ncomm=2 ivnum=-1 7005 write(*,7004,advance='NO') 7004 format(1x/ ' Add a vibrational quantum number (0/1) ? ') vv='00' read(*,'(i5)',err=7005)kk if(kk.eq.1)then 7007 write(*,7006,advance='NO') 7006 format( ' Value of the vibrational quantum number = ') read(*,'(i5)')m if(m.lt.0.or.m.gt.9)goto 7007 ivnum=m write(vv,'(2i1)')ivnum,ivnum endif c c...Loop over all lines c do 7000 kk=1,I linout='' if(ispace(nspace).eq.kk)then nspace=nspace+1 linout(1:2)=' !' endif c c...Transfer between the line comments: COMBLK -> LINCOM -> LINOUT c 7001 if(icomml(ncomm).eq.kk)then icbeg=icomm(ncomm-1)+1 icend=icomm(ncomm) if(icend-icbeg+1.ge.100)then write(*,7025)'***** WARNING: '// * 'Truncating a long .LIN annotation:', * comblk(icbeg:icend) 7025 format(1x/1x,a/1x,a/1x) icend=icbeg+99 endif write(lincom,'(a)')comblk(icbeg:icend) ncomm=ncomm+1 if(len_trim(linout).gt.0)then linout=linout(1:len_trim(linout))// * lincom(1:len_trim(lincom)) else linout=' '//lincom(1:len_trim(lincom)) endif goto 7001 endif c c...Append end of line annotations from DESCR to LINOUT c if(fr(kk,4).ne.0.d0)then write(linrel,'(1PE10.2)')fr(kk,4) else linrel=' 1.00000' endif lendes=len_trim(linout) if(lendes.eq.1.and.linout(l:1).le.' ')lendes=0 lendsc=len_trim(descr(kk)) if(lendsc.eq.1.and.descr(kk)(l:1).le.' ')lendsc=0 if(lendsc.gt.0)then if(lendes.eq.0)then linout=' #'//descr(kk)(1:lendsc) else linout=linout(1:lendes)//' #'//descr(kk)(1:lendsc) endif lendes=len_trim(linout) endif c errprn=errorl(kk) if(linfit(kk).eq.0)errprn=errorl(kk)+900.d0 c c...Output a .LIN line for standard asymmetric top quantisation c if(kindx.eq.6)then if(lendes.eq.0)then if(ivnum.lt.0)WRITE(3,7002)(JK(kK,M),M=1,6),FR(kK,1), * ERRprn,linrel if(ivnum.ge.0)WRITE(3,7003)(JK(kK,M),M=1,3),ivnum, * (JK(kk,m),m=4,6),ivnum,FR(kK,1),ERRprn, * linrel else if(ivnum.lt.0)WRITE(3,7002)(JK(kK,M),M=1,6),FR(kK,1), * ERRprn,linrel,linout(1:lendes) if(ivnum.ge.0)WRITE(3,7003)(JK(kK,M),M=1,3),ivnum, * (jk(kk,m),m=4,6),ivnum,FR(kK,1),ERRprn, * linrel,linout(1:lendes) endif endif c c...Output a .LIN line for symmetric top quantisation, the K that is used is c decided on the basis of NAMREP c if(kindx.eq.-6)then c if(namrep.eq.'III.r')then ksymu=jk(kk,3) if(((jk(kk,1)/2)*2.ne.jk(kk,1)).and. * (jk(kk,2)+jk(kk,3).eq.jk(kk,1)))ksymu=-jk(kk,3) if(((jk(kk,1)/2)*2.eq.jk(kk,1)).and. * (jk(kk,2)+jk(kk,3).eq.jk(kk,1)+1))ksymu=-jk(kk,3) ksyml=jk(kk,6) if(((jk(kk,4)/2)*2.ne.jk(kk,4)).and. * (jk(kk,5)+jk(kk,6).eq.jk(kk,4)))ksyml=-jk(kk,6) if(((jk(kk,4)/2)*2.eq.jk(kk,4)).and. * (jk(kk,5)+jk(kk,6).eq.jk(kk,4)+1))ksyml=-jk(kk,6) endif c if(namrep.eq.'I.r ')then ksymu=jk(kk,2) if(((jk(kk,1)/2)*2.ne.jk(kk,1)).and. * (jk(kk,2)+jk(kk,3).eq.jk(kk,1)))ksymu=-jk(kk,2) if(((jk(kk,1)/2)*2.eq.jk(kk,1)).and. * (jk(kk,2)+jk(kk,3).eq.jk(kk,1)+1))ksymu=-jk(kk,2) ksyml=jk(kk,5) if(((jk(kk,4)/2)*2.ne.jk(kk,4)).and. * (jk(kk,5)+jk(kk,6).eq.jk(kk,4)))ksyml=-jk(kk,5) if(((jk(kk,4)/2)*2.eq.jk(kk,4)).and. * (jk(kk,5)+jk(kk,6).eq.jk(kk,4)+1))ksyml=-jk(kk,5) endif c c...Preliminaries completed: now the actual .LIN line c if(lendes.eq.0)then if(ivnum.lt.0)WRITE(3,7011)JK(kk,1),ksymu,jk(kk,4),ksyml, * FR(kK,1),ERRprn,linrel if(ivnum.ge.0)WRITE(3,7002)JK(kk,1),ksymu,ivnum, * JK(kk,4),ksyml,ivnum,FR(kK,1),ERRprn, * linrel else if(ivnum.lt.0)WRITE(3,7011)JK(kk,1),ksymu,jk(kk,4),ksyml, * FR(kK,1),ERRprn,linrel,linout(1:lendes) if(ivnum.ge.0)WRITE(3,7002)JK(kk,1),ksymu,ivnum, * JK(kk,4),ksyml,ivnum,FR(kK,1),ERRprn, * linrel,linout(1:lendes) endif c endif c 7000 continue c 7011 format(4i3,' 0 0 0 0 0 0 0 0',f15.5,f12.5,2a) 7002 format(6i3,' 0 0 0 0 0 0',f15.5,f12.5,2a) 7003 format(8i3,' 0 0 0 0',f15.5,f12.5,2a) close(3) c c...output of a .PAR file for SPFIT C k=len_trim(filnam) filnam=filnam(1:k-4)//'.par' ispind=1 if(kindx.eq.-6)ispind=-1 if(ivnum.eq.-1)ivnum=0 OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) c write(3,'(72a1)')com do 1405 ncon=nconst,1,-1 if(a(ncon).ne.0.d0)goto 1406 1405 continue 1406 write(3,1407)ncon,i,5,0, * ' 0.0000E+000 1.0000E+002 1.0000E+000 1.0000000000' 1407 format(i5,i6,2i5,a) if(namred.eq.'A')then ! A-reduction if(namrep.eq.'I.r ')then ! prolate write(3,'(a,2i3,a)') * 'a ',ispind,ivnum+1,' 0, , , , , , , ,,,,,' else ! oblate write(3,'(a,2i3,a)') * 'a ',ispind,-(ivnum+1),' 0, , , , , , , ,,,,,' endif else ! S-reduction if(namrep.eq.'I.r ')then ! prolate write(3,'(a,2i3,a)') * 's ',ispind,ivnum+1,' 0, , , , , , , ,,,,,' else ! oblate write(3,'(a,2i3,a)') * 's ',ispind,-(ivnum+1),' 0, , , , , , , ,,,,,' endif endif c do 1408 kk=1,ncon picon(kk)=a(kk) if(kk.ge.4.and.kk.le.8)picon(kk)=-a(kk) if((kk.eq.7.or.kk.eq.8).and.namred.eq.'S')picon(kk)=a(kk) if(namrep.eq.'III.r')then if(kk.eq. 7)picon(kk)=-picon(kk) if(kk.eq. 8)picon(kk)=-picon(kk) if(kk.eq.13)picon(kk)=-picon(kk) if(kk.eq.14)picon(kk)=-picon(kk) if(kk.eq.15)picon(kk)=-picon(kk) if(kk.eq.21)picon(kk)=-picon(kk) if(kk.eq.22)picon(kk)=-picon(kk) if(kk.eq.23)picon(kk)=-picon(kk) if(kk.eq.24)picon(kk)=-picon(kk) if(kk.eq.31)picon(kk)=-picon(kk) if(kk.eq.32)picon(kk)=-picon(kk) if(kk.eq.33)picon(kk)=-picon(kk) if(kk.eq.34)picon(kk)=-picon(kk) if(kk.eq.35)picon(kk)=-picon(kk) endif if(namred.eq.'S'.and.kk.eq. 8)pinums(kk)='500' if(namred.eq.'S'.and.kk.eq.14)pinums(kk)='501' if(namred.eq.'S'.and.kk.eq.15)pinums(kk)='600' if(namred.eq.'S'.and.kk.eq.22)pinums(kk)='502' if(namred.eq.'S'.and.kk.eq.23)pinums(kk)='601' if(namred.eq.'S'.and.kk.eq.24)pinums(kk)='700' if(namred.eq.'S'.and.kk.eq.32)pinums(kk)='503' if(namred.eq.'S'.and.kk.eq.33)pinums(kk)='602' if(namred.eq.'S'.and.kk.eq.34)pinums(kk)='701' if(namred.eq.'S'.and.kk.eq.35)pinums(kk)='800' 1408 continue c do 1402 kk=1,ncon if(a(kk).eq.0.d0.and.ia(kk).eq.0)goto 1402 itstrt=1 if(t(kk)(1:1).eq.' ')itstrt=2 if(ia(kk).eq.0)then parwt=1.D-35 else parwt=1.D+20 endif if(kk.ge.4.and.kk.le.8)then if((kk.eq.7.or.kk.eq.8).and.namred.eq.'S')then write(3,1403)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) else write(3,1404)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) endif else write(3,1403)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) endif 1402 continue write(3,1410) close(3) NOUT=0 c 1403 format(11x,a5,1PE24.15,E16.8,' /',a) 1404 format(11x,a5,1PE24.15,E16.8,' /-',a) 1410 format('!'/'!',78(1H-)/'!'/ * '! Generated automatically by ASFIT'/'!'/ * '!',78(1H-)/) c write(*,'(1x/'' SUCCESSFUL OUTPUT OF: '',a/18x,''AND: '',a/)') * filnam(1:k-4)//'.lin',filnam(1:len_trim(filnam)) GOTO 23 ENDIF c GOTO 25 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...Output of an .ASF data file C 23 WRITE(*,14,advance='NO') 14 FORMAT(/' DATA FILE TO BE WRITTEN TO DISK (1/0) ? ') READ(*,15,ERR=23)L 15 FORMAT(I5) IOUT=1 IF(L.NE.1)GOTO 24 100 WRITE(*,17,advance='NO') 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,'(22(1H-),5(1H+),61(1H-))') 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,'(5(1H-),30(1H+),53(1H-))') do 103 k=1,ncon if(k.eq.1)then WRITE(2,150)IA(K),A(K),T(K)(1:5),' Reduction-', * namred,', Representation-',namrep else WRITE(2,150)IA(K),A(K),T(K)(1:5) endif 103 continue c write(2,55) write(2,56) write(2,554) 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 lendes=len_trim(descr(kk)) if(lendes.eq.0)then if(fr(kk,4).eq.0.d0)then WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERRORL(kK),LINFIT(kK) else WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERRORL(kK),LINFIT(kK), * FR(kk,4) endif else if(fr(kk,4).eq.0.d0)then WRITE(2,141)(JK(kK,M),M=1,6),FR(kK,1),ERRORL(kK),LINFIT(kK), * DESCR(kK)(1:lendes) else WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERRORL(kK),LINFIT(kK), * fr(kk,4),DESCR(kK)(1:lendes) endif endif 3400 continue write(2,554) write(2,55) CLOSE(2) c 51 FORMAT(1X,72A1) 53 format(79(1H-)) 553 format(102(1H-)) 554 format(3(5(1H+),5(1H-)),20(1H+),15(1H-),5(1H+),10(1H-),8(1H+)) 54 format(a22,i5) 55 format(88(1H-)) 56 format(3x,'J'' K-1'' K+1'' J K-1 K+1 frequency', * ' error fit? blend Note') 140 FORMAT(6I5,F20.6,F15.6,I5,1PE9.2,1x,A) c140 FORMAT(6I5,F20.6,F15.6,I5,F9.4,1x,A) 141 FORMAT(6I5,F20.6,F15.6,I5,10x,A) 150 FORMAT(I5,F30.19,3X,5A) C IOUT=1 GOTO 24 C 16 RETURN END C C_____________________________________________________________________________ C SUBROUTINE DFBET(Z,K1,N) C C DFBETAS for K1 fitted constants associated with the N'th transition C (the K1 calculated values are passed back in array Z) C C Changes in constants relative to last iteration are in B, their errors C are in C, D'Y is in K1+1'st column of H, (D'D) is set up in top-left C corner of H and inv(D'D) placed in DD C IMPLICIT REAL(8) (A-H,O-Z) PARAMETER (NLINES=15000, NCONST=35, JMAX=250, MAXDIM=JMAX/2+1, * NNJA=2*NLINES+1) COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR COMMON /BLKH/H * /BLKR/R * /BLK2/NJA,B * /BLK3/C,D,DD,HJJ INTEGER(4) IA(NCONST),LINFIT(NLINES),JK(NLINES,6),NJA(NNJA,4) REAL(8) FR(NLINES,4),A(NCONST),C(NCONST),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),DD(NCONST,NCONST),Z(NCONST), * B(NCONST),BNEW(NCONST),ERRORL(NLINES),WEIGHT(NLINES), * D(NLINES,NCONST),HJJ(NLINES) character(50) descr(nlines) C if(linfit(N).eq.0)goto 300 C C...Evaluate inv(D'D) from the old derivative matrix in D, but without C line N - inversion is through diagonalisation as in M/P C C array Z is used: C C 1. to hold scaling constants improving conditioning in evaluation of C inv(D'D) C 2. for DFBETAS associated with line N C DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(L.eq.N)GOTO 27 IF(LINFIT(L).eq.0)GOTO 27 IF(LINFIT(L).eq.1)E=E+D(L,J)*D(L,K)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*D(L,K)*WEIGHT(L) below=FR(L,4) do 530 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 531 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*D(jjj,K)*WEIGHT(jjj) below=below+fr(jjj,4) 530 continue 531 E=E+over/below endif ENDIF 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(L.eq.N)GOTO 29 IF(LINFIT(L).eq.0)GOTO 29 IF(LINFIT(L).eq.1)E=E+D(L,J)*FR(L,3)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*FR(L,3)*WEIGHT(L) below=FR(L,4) do 532 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 533 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*FR(jjj,3)*WEIGHT(jjj) below=below+fr(jjj,4) 532 continue 533 E=E+over/below endif ENDIF 29 CONTINUE H(J,K1+1)=E 30 CONTINUE C DO 301 J=1,K1 DO 1301 K=J,K1 ! Fortran2018 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 1301 CONTINUE ! Fortran2018 301 CONTINUE C IEGEN=0 CALL HDIAG(K1,IEGEN,L3) ! <--- hdiag C DO 35 J=1,K1 BNEW(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 BNEW(J)=BNEW(J)+E*H(K,K1+1) 34 CONTINUE 35 CONTINUE C C...DFBETAS (from difference of corrections B and BNEW to constants from C previous least squares iteration) C 300 if(LINFIT(N).eq.0)then do 360 L=1,K1 Z(L)=0.0d0 360 continue return else do 361 L=1,K1 Z(L)=( B(L)-BNEW(L) )/C(L) 361 CONTINUE c endif C RETURN END C_____________________________________________________________________________ C subroutine hjjset(K1,hjjav) C C Calculation of the HJJ coefficients describing the influence of transition C J on the fit: C transitions with hjj > 0.5 should be used with care C transitions with hjj > 0.9 strongly skew the fit and should really only C be tolerated in intermediate stages of constant determination C IMPLICIT REAL(8) (A-H,O-Z) PARAMETER (NLINES=15000, NCONST=35) COMMON /BLK1/JK,FR,ERRORL,WEIGHT,LINFIT,A,IA,I,IOUT,IRED,DESCR * /BLK3/C,D,DD,HJJ * /SORTCC/WORK,IPT INTEGER(4) IA(NCONST),LINFIT(NLINES),JK(NLINES,6),IPT(NLINES) REAL(8) FR(NLINES,4),A(NCONST),C(NCONST),DD(NCONST,NCONST), * WORK(NLINES),ERRORL(NLINES),WEIGHT(NLINES), * D(NLINES,NCONST),HJJ(NLINES) character(50) descr(nlines) C C...set up hjj's = diagonal elements of H = D inv(D'D) D' C do 600 L=1,I if(linfit(L).eq.0)goto 600 do 601 LL=1,K1 work(LL)=0.0 do 602 LLL=1,K1 work(LL)=work(LL)+D(L,LLL)*DD(LLL,LL) 602 continue 601 continue hjj(L)=0.0 do 603 LL=1,K1 hjj(L)=hjj(L)+work(LL)*D(L,LL) 603 continue hjj(L)=hjj(L)*weight(L) 600 continue C C...sort hjj's according to magnitude C hjjcut=0.0 hjjav=0.0 nhjj=0 do 605 L=1,I if(linfit(L).gt.0)then work(L)=hjj(L) hjjav=hjjav+hjj(L) nhjj=nhjj+1 else work(L)=0.0d0 endif ipt(L)=L 605 continue if(nhjj.gt.0)hjjav=hjjav/nhjj int21=1 call sortc(int21,i) ! <--- sortc c return end C C_____________________________________________________________________________ C subroutine lzero(number,ndig) character(20) number c c Addition of a leading zero to a number in a string c do 1 n=1,ndig-2 if(number(n:n+1).eq.' .')number(n:n+1)='0.' if(number(n:n+2).eq.' -.')number(n:n+2)='-0.' 1 continue c return end C_____________________________________________________________________________ C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR,NERD) C C Constant and error formatting for output c c CONVAL - String containing the constant value. This is to be c input in F format and will be replaced on output by the result c string of length extending to the last digit of the error c ERVAL - String containing the error value. This is to be c input in F format and will be replaced on output by the result c string, which does not contain the decimal point and is meant c to be included in brackets c NDCON - The number of digits in the CONVAL string (inclusive of any c leading zeros) c NDEROR - The number of digits in the ERVAL string (just the significant c digits) and is either equal to NERD, or is larger if there c are more significant digits than NERD before the decimal point. c NERD - the number of desired error digits, which is to be set on input. c c Both NDCON and NDEROR are generated on output C CHARACTER(27) CONVAL,ERVAL,CONOUT,EROUT C NDEROR=0 NDNOTZ=0 ICZERO=0 erout='???????????' 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. C Terminate when either NERD digits in the error value are C transferred or, if error has more digits before the decimal C point, the decimal point is reached. C conout(1:2)=conval(1:2) ndcon=2 DO 1 N=3,27 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(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.' ') * CONOUT(NDCON-1:NDCON-1)='0' 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 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 NERD 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 c c...output string with fitted value (rounding is to be carried out externally) c 2 CONVAL(1:NDCON)=CONOUT(1:NDCON) c c...output string with error (rounded if necessary) c if(nderor.gt.0)then c write(*,'(1x,2a,5x,a)') DEBUG c * 'ERROR: ',erval,erout(1:nderor) DEBUG c if(erval(ndcon+1:ndcon+1).eq.'5'.or. * erval(ndcon+1:ndcon+1).eq.'6'.or. * erval(ndcon+1:ndcon+1).eq.'7'.or. * erval(ndcon+1:ndcon+1).eq.'8'.or. * erval(ndcon+1:ndcon+1).eq.'9')then read(erout(1:nderor),*)ieror net=int( dlog10(dble(ieror)) ) ieror=ieror+1 net1=int( dlog10(dble(ieror)) ) if(net1-net.gt.0)nderor=nderor+1 write(erout,'(i12)')ieror erval(1:nderor)=erout(12-nderor+1:12) else ERVAL(1:NDEROR)=EROUT(1:NDEROR) endif endif C RETURN END C C_____________________________________________________________________________ c subroutine confix(const,ascld,conval,iend) c C Fixed constant formatting for output c implicit real(8) (a-h,o-z) character conval*27 C if(abs(const).le.1.d-17)ascld=0.d0 IF(ASCLD.NE.0.D0)ASCLD=ASCLD+MOD(ASCLD,1.D-15*ASCLD) WRITE(CONVAL,'(F27.16)')ASCLD c c...deal with leading zero IF(ASCLD.EQ.0.0D0)CONVAL=' 0.0 ' IF(CONVAL(10:10).EQ.' ')CONVAL(10:10)='0' IF(CONVAL(10:10).EQ.'-')THEN CONVAL(10:10)='0' CONVAL(9:9)='-' ENDIF c c...set leading square bracket do 2115 m=1,27 if(CONVAL(M:M).NE.' ')GOTO 2116 2115 continue 2116 if(m.gt.8)then conval(8:8)='[' else if(m.gt.1.and.m.le.8)then conval(m-1:m-1)='[' endif c c...blank at the end so that number processed is within 15 digit accuracy if(ascld.ne.0.d0)then mmm=int( alog10(real(abs(ascld))) ) if(mmm.ge.-1)mmm=mmm+1 do 2117 m=27,27-mmm,-1 conval(m:m)='0' 2117 continue endif c c...find last non-zero digit DO 1114 M=27,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 C C...terminating square bracket C 1115 iend=m+1 if(iend.gt.27)iend=27 conval(iend:iend)=']' c return end c C_____________________________________________________________________________ C SUBROUTINE ASREDA(I,J,N,A,IRED) C C Evaluation of a Wang submatrix for reduction A, where: C C I - on output contains the dimensionality of the 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) PARAMETER (NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) INTEGER(4) JJ1,K,KK COMMON /BLKH/H * /BLKR/U DIMENSION A(NCONST),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C JJ1=J JJ1=JJ1*(JJ1+1) Q=JJ1 QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C IF(IRED.GT.0)THEN C=0.5D0*(A(2)+A(3)) DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ELSE C=0.5D0*(A(2)+A(1)) DK2=A(3)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ENDIF DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ 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 IF(IRED.GT.0)THEN OK0=0.25D0*(A(2)-A(3))-A(7)*Q+A(13)*QQ+A(21)*QQQ+ * A(31)*QQQQ ELSE OK0=0.25D0*(A(1)-A(2))-A(7)*Q+A(13)*QQ+A(21)*QQQ+ * A(31)*QQQQ ENDIF OK2=0.5D0*(-A(8)+A(14)*Q+A(22)*QQ+A(32)*QQQ) OK4=0.5D0*(A(15)+A(23)*Q+A(33)*QQ) OK6=0.5D0*(A(24)+A(34)*Q) OK8=0.5D0*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 1001 I=1,M ! Fortran2018 DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE 1001 CONTINUE ! Fortran2018 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 FK2=K**2 H(I,I)=((((DK10*FK2+DK8)*FK2+DK6)*FK2+DK4)*FK2+DK2)*FK2 * +DK0 IF(K.LT.KE)GOTO 2 GOTO 3 2 C=(K+2)**2 FK8=FK2**4+C**4 FK6=FK2**3+C**3 FK4=FK2*FK2+C*C FK2=C+FK2 FLTB=JJ1-K*KK FLTC=JJ1-KK*(K+2) F2=FLTB*FLTC H(I,I+1)=(OK0+FK2*OK2+FK4*OK4+FK6*OK6+FK8*OK8)*DSQRT(F2) H(I+1,I)=H(I,I+1) 3 CONTINUE C C Final touch up: 4 -> O+ and O-, 6 -> E-, 5 -> E+ C c IF(2-N)4,6,5 IF(2-N.ge.0)then ! Fortran2018 if(2-N.gt.0)goto 5 ! Fortran2018 goto 6 ! Fortran2018 endif ! Fortran2018 C C FLTA = +1 for O+, -1 for O- C FLTA=-2*N+7 ! old label 4 H(1,1)=FLTA*(OK0+2.D0*(OK2+OK4+OK6+OK8))*Q+H(1,1) GOTO 6 C 5 H(1,2)=H(1,2)*1.414213562373095D0 H(2,1)=H(1,2) 6 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ASREDS(I,J,N,A,IRED) C C Evaluation of a Wang submatrix for reduction S, where: C C I - on output contains the dimensionality of the 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) PARAMETER (NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) INTEGER(4) K COMMON /BLKH/H * /BLKR/U DIMENSION A(NCONST),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C Q=REAL(J+1)*REAL(J) QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C IF(IRED.GT.0)THEN C=0.5D0*(A(2)+A(3)) DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ELSE C=0.5D0*(A(2)+A(1)) DK2=A(3)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ENDIF DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ 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 IF(IRED.GT.0)THEN OKK2=0.25D0*(A(2)-A(3))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ ELSE OKK2=0.25D0*(A(1)-A(2))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ ENDIF 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 1001 I=1,M ! Fortran2018 DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE 1001 CONTINUE ! Fortran2018 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 FK2=REAL(K)**2 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 c IF(2-N)4,12,5 IF(2-N.ge.0)then ! Fortran2018 if(2-N.gt.0)goto 5 ! Fortran2018 goto 12 ! Fortran2018 endif ! Fortran2018 c c...O+ and O-, FLTA = +1 for O+, -1 for O- FLTA=-2*N+7 ! old label 4 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) 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) 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 diagonalised 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 diagonalised 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) 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 1014 I=1,N ! Fortran2018 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 1014 CONTINUE ! Fortran2018 15 NR=0 IF(N-1.LE.0)GOTO 1000 NMI1=N-1 DO 1030 I=1,NMI1 ! Fortran2018 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 1030 CONTINUE ! Fortran2018 40 DO 70 I=1,NMI1 IF(I-1.LE.0)GOTO 60 c IF(XMAX-X(I))60,70,70 IF(XMAX-X(I).ge.0.0d0)goto 70 ! Fortran2018 60 XMAX=X(I) IPIV=I JPIV=IQ(I) 70 CONTINUE EPS=1.0D-12 c IF(XMAX-EPS)1000,1000,148 IF(XMAX-EPS.le.0.d0)goto 1000 ! Fortran2018 c148 CONTINUE ! no longer necessary NR=NR+1 c130 IF(H(IPIV,IPIV)-H(JPIV,JPIV))140,141,142 IF(H(IPIV,IPIV)-H(JPIV,JPIV).ge.0.d0)then ! Fortran2018 old label 130 IF(H(IPIV,IPIV)-H(JPIV,JPIV).gt.0.d0)goto 142 ! Fortran2018 goto 141 ! Fortran2018 endif ! Fortran2018 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) ! old label 140 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 c IF(H(IPIV,IPIV)-H(JPIV,JPIV))153,153,152 IF(H(IPIV,IPIV)-H(JPIV,JPIV).le.0.d0)goto 153 ! Fortran2018 HTEMP=H(IPIV,IPIV) ! old label 152 H(IPIV,IPIV)=H(JPIV,JPIV) H(JPIV,JPIV)=HTEMP HTEMP=SIGN(1.D0,-SINE)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 CONTINUE c DO 350 I=1,NMI1 c IF(I-IPIV)210,350,200 IF(I-IPIV.le.0)then ! Fortran2018 if(I-IPIV.lt.0)goto 210 ! Fortran2018 goto 350 ! Fortran2018 ENDIF ! Fortran2018 c200 IF(I-JPIV)210,350,210 IF(I-JPIV.eq.0)goto 350 ! Fortran2018 old label 200 210 IF(IQ(I)-IPIV.EQ.0)GOTO 240 c IF(IQ(I)-JPIV)350,240,350 IF(IQ(I)-JPIV.ne.0)goto 350 ! Fortran2018 240 K=IQ(I) HTEMP=H(I,K) ! old label 250 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 c X(IPIV)=0.D0 X(JPIV)=0.D0 DO 530 I=1,N c IF(I-IPIV)370,530,420 IF(I-IPIV.ge.0)then ! Fortran2018 if(i-ipiv.gt.0)goto 420 ! Fortran2018 goto 530 ! Fortran2018 endif ! Fortran2018 HTEMP=H(I,IPIV) ! old label 370 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 c420 IF(I-JPIV)430,530,480 420 IF(I-JPIV.ge.0)then ! Fortran2018 if(i-jpiv.gt.0)goto 480 ! Fortran2018 goto 530 ! Fortran2018 endif ! Fortran2018 HTEMP=H(IPIV,I) ! old label 430 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 X(JPIV)=ABS(H(JPIV,I)) ! old label 510 IQ(JPIV)=I 530 CONTINUE c 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 SUBROUTINE HSORT(M) C C SORTING OF EIGENVALUES AND EIGENVECTORS COMING OUT OF HDIAG C - introduced after it was found that for representation III some C HDIAG eigenvalues did not keep ascending order. C IMPLICIT REAL(8) (A-H,O-Z) PARAMETER (JMAX=250,MAXDIM=JMAX/2+1,NLINES=15000) C COMMON /BLKH/A * /BLKR/U * /SORTCC/WK,IPT REAL(8) A(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),WK(NLINES),D(MAXDIM) INTEGER(4) IPT(NLINES) isort=1 C C...Check whether additional sorting necessary C DO 150 I=1,M-1 IF(A(I+1,I+1).LT.A(I,I))GOTO 151 150 CONTINUE RETURN C C...Transfer eigenvalues from the diagonal of matrix A to D, and eigenvectors C from U to A C 151 DO 115 I=1,M D(I)=A(I,I) DO 1115 J=1,M ! Fortran2018 A(I,J)=U(I,J) 1115 CONTINUE ! Fortran2018 115 CONTINUE C C...Straightforward magnitude sorting C DO 122 I=1,M IPT(I)=I WK(I)=D(I) 122 CONTINUE CALL SORTC(isort,M) ! <--- sortc C C...Return of sorted eigenvectors in U and eigenvalues on the diagonal of A C DO 117 II=1,M ! old label 140 I=IPT(II) DO 1117 J=1,M ! Fortran2018 U(J,II)=A(J,I) 1117 CONTINUE ! Fortran2018 117 CONTINUE DO 118 II=1,M I=IPT(II) A(II,II)=D(I) 118 CONTINUE C RETURN END C C____________________________________________________________________________ C SUBROUTINE SORTC(N,M) PARAMETER (NLINES=15000) C COMMON /SORTCC/WK,IPT INTEGER(4) IPT(NLINES) REAL(8) WK(NLINES),EE C C ... This routine sorts the quantities in vector WK in ascending order C of magnitude and also accordingly rearranges vector IPT of pointers C to original positions of sorted quantities C DO 101 I=N,M-1 J=I 106 J=J+1 c IF(WK(J)-WK(I))103,104,104 IF(WK(J)-WK(I).ge.0.d0)goto 104 ! Fortran2018 EE=WK(I) ! old label 103 WK(I)=WK(J) WK(J)=EE K=IPT(I) IPT(I)=IPT(J) IPT(J)=K 104 IF(J.EQ.M)GOTO 101 GOTO 106 101 CONTINUE C RETURN END C c_____________________________________________________________________________ C C Check whether the string contains any non-FORTRAN input characters C integer function icheck(charg) character charg*(*) c nn=len_trim(charg) do 1 n=1,nn i=ichar(charg(n:n)) if(i.lt.32.or.i.gt.126)goto 2 1 continue c icheck=0 return 2 icheck=i return end C_____________________________________________________________________________ C_____________________________________________________________________________ C C Routines in the order in this listing: C C MP C DATAIN - input of data C MODIF - set up internal arrays and online modifications C GENER - generate blocks of lines for specified bands C DATOUT - output of data C DFBET - calculate DFBETA coefficients C hjjset - calculate HJJ coefficients C lzero - add leading zero to an number in an output string (if necessary) C CONFOR - format fitted constants for printout C CONFIX - format fixed constants for printout C ASREDA - set up A reduced Hamiltonian C ASREDS - set up S reduced Hamiltonian C FPLUS - evaluate the Fplus function for the Hamiltonian C FMINUS - evaluate the Fminus function for the Hamiltonian C HDIAG - diagonalise the Hamiltonian matrix C HSORT - sort eigenvalues and eigenvectors out of HDIAG C SORTC - simple sorting of a vector in ascending order of magnitude C icheck - function to check string for non-FORTRAN characters C C_____________________________________________________________________________ C_____________________________________________________________________________