module assig integer eseign(0:2,0:20,-20:20),ipoint integer aseign(0:2,0:20,0:1,0:20),ipoint1 end module assig PROGRAM ROTFIT c my insertion USE MSFLIB USE MSIMSL c also it is necessary to change name of the array from deriv to deriv1 c in order to avoid redeclaration (deriv is declared in MSFLIB) IMPLICIT REAL*8 ( A-H,O-Z ) * PROGRAM ROTFIT (INPUT,OUTPUT,TAPE31,TAPE32,TAPE33,TAPE34,TAPE35, * 1 TAPE36,TAPE5=INPUT,TAPE6=OUTPUT) C******************************************************************* C I.KLEINER AND M.GODEFROID (FREE UNIVERSITY OF BRUSSELS, C 50, AV. F-D. ROOSEVELT, 1050 BRUSSELS, BELGIUM, 1987-91) ********************************************************************* C 20 aout 1992 C C THIS PROGRAM IS A LEAST-SQUARES PROCEDURE ABLE TO FIT C DATA FROM A MOLECULE WITH TORSION MOOVING.IT IS BASED C ON THE LEAST-SQUARES PROGRAM WRITTEN BY M.GODEFROID C "RTFIT" AND ON THE ROTATOR-TORSION PROGRAM"ROTTOR". C THE V MATRIX ELEMENTS ARE CALCULATED AS THE DERIVATIVES C OF THE ENERGY TO THE PARAMETER TO BE FITTED C THIS VERSION CAN TRAIT MICROWVE OR INFRARED DATA C THIS VERSION IS THE COPY OF TORFIT BUT WAS MODIFIED IN C ORDER TO RUN IN THE CDC NBS COMPUTER AND IS A WEIGHT LEAST-SQUARES C THIS PROGRAM INCLUDES V9 TERM,DABJ*J*(J+1),DABK*K**2,RHOJ AND C RHOK C FITV9B CAN TRAIT HIGH J VALUES THANKS A SPECIAL PROCEDURE OF C PICKING UP EIGENVALUES THAT INVOLVES A SUM OF THE SQUARE OF C COEFFIECIENTS C(K) ET C(-K) (SEE FUNCTION IPOSA AND IPOSE) C IN THIS VERSION ONLY THE WEIGHTS OF THE MW HAS CHANGED C (EACH MW LINE CAN RECEIVE A DIFFERENT WEIGHT). C IN THIS VERSION THE OUTPUT HAS CHANGED (SEE LINECAL) C NEW TERMS HAVE BEEN ADDED:DC1J,C2J,AK2J,GVJ,AMV,ANV,BK2 C C DAC,DBC,C4,ODELTA,AK3J,AK3K C RHO HAS BEEN REPLACED BY RHORHO AND IS ONLY FITTED FOR HTORS. C FITWGT9 IS A VERSION WHERE THE DIAGONAL MATRIX ELEMENTS IN K C AND OFF-DIAG. IN VT HAVE BEEN ADDED,CONTRARLY TO HERBST. C NEW PARAMETER IN FITWG10:D DELTA C FTJON10:D JON HAS MODIFIED THE PROGRAM... C AND I HAVE ADD AK1J,AK1K,AK2K,AK5K C FASTFIT:VERSION THAT HAS BEEN MODIFIED IN ORDER TO GAIN TIME C (SEE SUBROUTINES SETUP AND VSET). C THE ENERGIES AND EIGENVECTORS ARE CALCULATED AT THE VERY BEGINNI C OF EACH CYCLE FOR J=0 TO 14 FOR EVERY K AND VT STATES AND THEN C THOSE ENERGIES AND EIGENVECTORS ARE STORED IN 2 BIG MATRICES C (DERIV AND ENER) C THIS IS THE LAST (IMPROVED) VERSION OF GIROU (CIRCE):IKSRC18 C FOR VT=2+IPOSE,IPOSA MODIFIED (IKSRC20) C belgi28: changes by Ortigoso, AOUT 1992 C belgi30 : changes june 1993; the IPOSA has now a parity calc C form (23)* operator and the basis set has been extended from 7 to 9 C CAN TREAT UP TP VT=4 C belgi32 : has deltab and odeltb C belgi32b : has c4J,C4K,C1K,C2K and C12 C belgi37 : trying to rotate the eigenvectors to fix the labellingstuff C belgi38 ;labelling scheme for vt=0,1,2 same as in belgi32c C but for vt=3 call WIGNER and energy surface results C PERFORMS DIAGONALISATION(CRAY SUBROUTINES) C belgi42.f : DACJ,DBCJ,ODABJ,RHOB,RHOC c c 23 sep 99 JTH has made some minor changes. See cjth lines. IBLK is c no longer necessary in the input file. C************************************************************************ CHARACTER ABC*702 CHARACTER*6 PARA,REF CHARACTER*1 PR,PROBS,PRINF,PRSUP PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/EIGEN1/EGVL(2*NDMX),EGVC(NDMX,2*NDMX) COMMON/TOR/A(9,NTORMX,2*KDMX),ETOR(2*KDMX,9) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/PAR/IND(2,80),NP(2),X(160),ITO COMMON/BIG2/DERIV1(10,961,80),ENER(10,961) COMMON/EXPDAT / ETRANS(80000),IVOBS(2,80000), * NOBS(2,80000),KOBS(2,80000),IBLK(2,80000),IW(80000),W(80000), * NDATA,NADAT,EPSI,NOFIT,NFITMW,NFITIR,NFIMW0 ,NFIMW1 ,NITT, * IAST(80000),NFI004 ,NFI080 ,NFI100 ,NFI020 ,NFI200 ,NFI045 , & NGST0 ,NFI010 ,SW070,NFI070 ,NFIMW2,IKAKC,JMAX,NFI1M0, & NFIR10,NFIR21,NFIR32,NFIMW3,NFIMW4,SWA,SWE,NFITMWA,NFITMWE COMMON/CNTL/NIT,IEND,DELTAS,S,TME(80000),TERMST,NP12,SMW,SIR, 1 SMW0,SMW1,SIG004,SIG080,SIG10,SIG100,SIG020,SIG200,SIG045, & SGST0,SMW2,SIG1M0,SMW3,SIG32,SMW4 COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117),VTLDV(160,160),VTLDTM (160),DELTAX (160) DIMENSION UNIT(160,160),V(160),V1(160),V2(160) DIMENSION WKSPCE(160),DELX(160),DELXST(160),IBLSTO(2) DIMENSION IPIV(160) INTEGER N, LDA, LWORK, IFAIL, INFO c my insertion integer*2 ihr, imin, isec, i100th DIMENSION UL(160,160),VR(160,160),SVD(160),COV(160,160) DIMENSION COR(160,160) c my insertion C C THIS PROGRAM MAY TRAIT (80000) OBSERVED DATA AND FIT (160) C PARAMETERS (UPPER AND/OR LOWER). C C CALL THE SUBROUTINE WHICH READ THE DATA' C C SEC=TCPU() C WRITE(6,*) ' BEGIN PROGRAM ',SEC c my insertion c open(unit=10,file='etype.txt',access='sequential',status='unknown') c open(unit=11,file='Atype.txt',access='sequential',status='unknown') WRITE(6,*) ' BEGIN PROGRAM ' c my insertion time CALL GETTIM (ihr, imin, isec, i100th) write(6,1012)ihr, imin, isec, i100th 1012 format(1x,'time',1x,i2,':',i2,':',i2,':',i2) c my insertion CALL DATA1(CST) C TERM= EVIB(2)-EVIB(1) TERMST=TERM C C BEGINNING OF THE FITTING PROCEDURE C NP12=NP(1)+NP(2) NPP=NP12 C PRINT*, 'NPP= ',NPP,' NP12= ',NP12 IF(ITO.NE.0) NPP=NP12+1 C C INITIALISATIONS C 23 NIT=1 SST=0. DO 30 I=1,NPP DELXST(I)=0.0 30 CONTINUE 25 DO 2 I=1,NPP V(I)=0.0 VTLDTM (I)=0.0 DO 3 J=1,NPP VTLDV(I,J)=0.0 3 CONTINUE 2 CONTINUE IBLSTO(1)=0.0 IBLSTO(2)=0.0 S=0.0 SMW=0.0 SWA=0.0 SWE=0. SIR=0.0 SIR10=0.0 SIR21=0.0 SIR32=0.0 SMW0=0.0 SMW1=0. SMW2=0. SMW3=0. SMW4=0. SIG004=0.0 SIG080=0.0 SIG010=0.0 SIG100=0.0 SIG020=0.0 SIG200=0.0 SIG045=0.0 SIG1M0=0.0 SGST0=0.0 SW070=0.0 NOFIT=0 NFITMW=0 NFITIR=0 NFIR10=0 NFIR21=0 NFIR32=0 NFIMW0 =0 NFIMW1 =0 NFIMW2=0 NFI004 =0 NFI080 =0 NFI010 =0 NFI100 =0 NFI020 =0 NFI200 =0 NFI045 =0 NFI1M0 =0 NGST0=0 NFI070 =0 NFIMW3=0 NFIMW4=0 NFITMWA=0 NFITMWE=0 IF(NP12.GT.(160)) THEN WRITE(6,*) 'TOO MANY PARAMETERS TO FIT' STOP ENDIF C C BEGINNING OF THE DO LOOP ON DATA,IF ONE DOESN'T WANT TO C FIT A DATA,IAST CAN'T BE EQUAL TO 1.NDATA IS THE TOTAL NUMBER C OF DATA INCLUDINGSIGMA=0 AND 1. C C SEC=TCPU() C WRITE(6,*) ' BEGIN DO 4 ',SEC WRITE(6,*) ' BEGIN DO 4 ' c my insertion time CALL GETTIM (ihr, imin, isec, i100th) write(6,1012)ihr, imin, isec, i100th c my insertion c my insertion c VTLDV - A matrix c VTLDTV - b vector c npp - number of parameters c K - index for measured transitions c I - index for fitted parameters c V(I) - derivative for given parameter c W(K) - weight of given measurement c CST(II,K) - old value c DELX(I) - variation of Ith parameter c XX - new value c DELTAX(I) - standard deviation c R - variation/standard deviation c RJTH - constant/standard deviation c my insertion C WE HAVE MODIFIED THE PROCEDURE HERE:DIN THE MATRIX DERIV, ONEA C IS CALCULATING THE ENERGIES AND DERIVATIVES FOR ALL THE J,K LEVELS C AT ONCE AND THEN ONE WILL PICK UP THE "GOOD" VALUES. CALL SETUP(CST) DO 4 K=1,NDATA IF (NOBS(1,K).GT.JMAX.OR.NOBS(2,K).GT.JMAX) GO TO 4 IF (IAST(K).EQ.0) GO TO 4 NOFIT=NOFIT+1 IF(IW(K).GE.1000) NFITMW=NFITMW+1 cjth added two cards here. isig=0 if(k.gt.nadat) isig=1 IF(IW(K).GE.1000.AND.ISIG.EQ.0) NFITMWA=NFITMWA+1 IF(IW(K).GE.1000.AND.ISIG.EQ.1) NFITMWE=NFITMWE+1 IF(IW(K).EQ.1.OR.IW(K).EQ.2.OR.IW(K).EQ.3) NFITIR=NFITIR+1 IF(IW(K).EQ.1) NFIR10=NFIR10+1 IF(IW(K).EQ.2) NFIR21=NFIR21+1 IF(IW(K).EQ.3) NFIR32=NFIR32+1 IF(IW(K).GE.1000.AND.IVOBS(1,K).EQ.0) NFIMW0 =NFIMW0 +1 IF(IW(K).GE.1000.AND.IVOBS(1,K).EQ.1) NFIMW1 =NFIMW1 +1 IF(IW(K).GE.1000.AND.IVOBS(1,K).EQ.3) NFIMW3 =NFIMW3 +1 IF(IW(K).GE.1000.AND.IVOBS(1,K).EQ.4) NFIMW4 =NFIMW4 +1 IF(IW(K).GE.1000.AND.IVOBS(1,K).EQ.2) NFIMW2=NFIMW2+1 C IF(IW(K).EQ.2000) NFITGD=NFITGD+1 IF(IW(K).EQ.10) NGST0=NGST0+1 IF(IW(K).EQ.1004.OR.IW(K).EQ.1005.OR.IW(K).EQ.1003) 1 NFI004 =NFI004 +1 IF(IW(K).EQ.1080) NFI080 =NFI080 +1 IF(IW(K).EQ.1010.OR.IW(K).EQ.1008) * NFI010 =NFI010 +1 IF(IW(K).EQ.1070) NFI070 =NFI070 +1 IF(IW(K).EQ.1100) NFI100 =NFI100 +1 IF(IW(K).EQ.1020) NFI020 =NFI020 +1 IF(IW(K).EQ.1200.OR.IW(K).EQ.1180.OR.IW(K).EQ.1160.OR. * IW(K).EQ.1150.OR.IW(K).EQ.1130) NFI200 =NFI200 +1 IF(IW(K).EQ.1030.OR.IW(K).EQ.1040.OR.IW(K).EQ.1050) * NFI045 =NFI045 +1 IF(IW(K).EQ.1990) NFI1M0 =NFI1M0 +1 cjth I need to get the derivatives, even if NITT=0, in order to generate a new cjth table similar to the one that Juan gave me. cjuan This command must be out to get Juan's table. IF(NITT.EQ.0) GO TO 4 C C SETUP THE MATRIX V,FOR CURRENT DATA DO 3000 I=1,NP12 V1(I)=0.0 V2(I)=0.0 3000 CONTINUE C C HERE WE ARE CALCULATING THE CORRECT INDICES IN ORDER TO FETCHA C THE ENERGIES AND DERIVATIVES CORRESPONDING TO THE OBSERVED TRANSIT C C II=1 C CALL VSET(II,V1,ELOWER,K,CST) C II=2 C CALL VSET(II,V2,EUPPER,K,CST) C DO 26 I=1,NP12 C PRINT*,' V1= ',V1(I),' V2= ',V2(I),' I= ',I DO 331 LO=1,2 N=NOBS(LO,K) KK=KOBS(LO,K) IV=IVOBS(LO,K) PR=PROBS(LO,K) IF(PR.EQ.'+') IPR=1 IF(PR.EQ.'-') IPR=0 IF(K.GT.NADAT) THEN ISIG=1 IVV=IV+6 INK=N**2+N+KK+1 ELSE ISIG=0 IVV=IV+1 INK=N**2+N+KK+1-IPR*N ENDIF IF(LO.EQ.1) THEN ELOWER=ENER(IVV,INK) ELSE EUPPER=ENER(IVV,INK) ENDIF DO 26 I=1,NP12 IF(LO.EQ.1) THEN V1(I)= DERIV1(IVV,INK,I) ELSE V2(I)=DERIV1(IVV,INK,I) ENDIF V(I)=V2(I)-V1(I) 26 CONTINUE 331 CONTINUE cjth Here I attempt to construct a derivative table for dehyperfining cjth similar to what Juan gave me. Be sure to comment out below cjuan's! itable=0 c itable=1 if(itable.eq.0) go to 332 if(nitt.ne.0) go to 332 ejuan = (eupper+term-elower)*29979.2458 write(6,1332) nobs(2,k),kobs(2,k),nobs(1,k),kobs(1,k), 1 ejuan,v2(1),v2(2),v2(3),v2(4),(v1(i),i=1,4) 1332 format(4i3,f15.2,8f9.4) 332 continue C C IF THE VIB.ENERGY HAS TO BE FIT,THE DIAG ELEMENT OF THE C V MATRIX WILL CONTAIN 1. C HERE,IT SEEMS TO BE AN ERROR ON THE PREVIOUS WAY OF CALCULATION! C IF(ITO.NE.0) THEN V(NPP)=1. C V(NPP)=0.0 C DO 8 I=1,NROTOR C V(NPP)=V(NPP)+EGVC(I,NBR(2))**2 C PRINT*,' V(NPP)= ',V(NPP) ENDIF C C SET UP THE VTILDE V MATRIX(NP*NP) C DO 10 I=1,NPP DO 810 J=1,NPP VTLDV(I,J)=VTLDV(I,J)+V(I)*V(J)*W(K) 810 CONTINUE 10 CONTINUE C C OBSERVED MINUS CALCULATED C TME(K)=ETRANS(K)-TERM-EUPPER +ELOWER C STANDARD DEVIATION FOR EACH LINE,WEIGHTED S=S+TME(K)**2*W(K) if(iw(K).GT.1000) then etrans2=etrans(K)*29979.2458 else etrans2=etrans(K) endif c print*,'s=',s,'etrans=',etrans2,'K=',k cjth added two if's here. IF(IW(K).GE.1000.and.k.le.nadat) then SWA=SWA+TME(K)**2*W(K) ELSEIF(IW(K).GE.1000.and.k.gt.nadat) then SWE=SWE+TME(K)**2*W(K) end if IF(IW(K).GE.1000) THEN SMW=SMW+TME(K)**2*W(K) IF(IVOBS(1,K).EQ.0.AND.IW(K).GE.1000) THEN SMW0=SMW0+TME(K)**2*W(K) ELSEIF(IVOBS(1,K).EQ.1.AND.IW(K).GE.1000) THEN SMW1=SMW1+TME(K)**2*W(K) ELSEIF(IVOBS(1,K).EQ.2.AND.IW(K).GE.1000) THEN SMW2=SMW2+TME(K)**2*W(K) ELSEIF(IVOBS(1,K).EQ.3.AND.IW(K).GE.1000) THEN SMW3=SMW3+TME(K)**2*W(K) ELSEIF(IVOBS(1,K).EQ.4.AND.IW(K).GE.1000) THEN SMW4=SMW4+TME(K)**2*W(K) ENDIF C ELSEIF (IW(K).EQ.2000) THEN C SWGOOD=SWGOOD+TME(K)**2*W(K) IF(IW(K).EQ.1004.OR.IW(K).EQ.1005.OR.IW(K).EQ.1003) THEN SIG004=SIG004+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1080) THEN SIG080=SIG080+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1010.OR.IW(K).EQ.1008) THEN SIG010=SIG010+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1100) THEN SIG100=SIG100+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1020) THEN SIG020=SIG020+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1070) THEN SW070=SW070+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1200.OR.IW(K).EQ.1180.OR.IW(K).EQ. 1 1160.OR.IW(K).EQ.1150.OR.IW(K).EQ.1130) THEN SIG200=SIG200+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1030.OR.IW(K).EQ.1040.OR.IW(K).EQ. 1 1050) THEN SIG045=SIG045+TME(K)**2*W(K) ELSEIF(IW(K).EQ.1990) THEN SIG1M0=SIG1M0+TME(K)**2*W(K) ENDIF ELSEIF(IW(K).EQ.1) THEN SIR10=SIR10+TME(K)**2*W(K) SIR=SIR+TME(K)**2*W(K) ELSEIF(IW(K).EQ.2) THEN SIR21=SIR21+TME(K)**2*W(K) SIR=SIR+TME(K)**2*W(K) ELSEIF(IW(K).EQ.3) THEN SIR32=SIR32+TME(K)**2*W(K) SIR=SIR+TME(K)**2*W(K) ELSEIF(IW(K).EQ.10) THEN SGST0=SGST0+TME(K)**2*W(K) ENDIF C IF(IW(K).GE.1000) TME(K)=TME(K)*29979.2458 IF(IBGTME.EQ.1) WRITE(6,300) K,TME(K),EUPPER,ELOWER,S 300 FORMAT(' K= ',I5,2X,' TME= ',F17.9,2X,' EUPPER= ',F16.8, 1 2X,' ELOWER= ',F16.8,' S= ',F20.8) DO 11 I=1,NPP VTLDTM (I)=VTLDTM (I)+V(I)*TME(K)*W(K) 11 CONTINUE C C END OF THE LOOP ON DATA C 4 CONTINUE C SEC=TCPU() C WRITE(6,*) ' END DO 4 ',SEC WRITE(6,*) ' END DO 4 ' c my insertion time CALL GETTIM (ihr, imin, isec, i100th) write(6,1012)ihr, imin, isec, i100th c my insertion IF(NITT.EQ.0) GO TO 244 PRINT*,'S(BEFORE FITTING)= ',S S=SQRT(S/(NOFIT-NPP)) IF(NFITMW.EQ.0) THEN SMW=0.0 ELSE SMW=SQRT(SMW/NFITMW) ENDIF IF(NFITMWA.EQ.0) THEN SWA=0.0 ELSE SWA=SQRT(SWA/NFITMWA) ENDIF IF(NFITMWE.EQ.0) THEN SWE=0.0 ELSE SWE=SQRT(SWE/NFITMWE) ENDIF IF (NFITIR.EQ.0) THEN SIR=0.0 ELSE SIR=SQRT(SIR/NFITIR) ENDIF IF (NFIR10.EQ.0) THEN SIR10=0.0 ELSE SIR10=SQRT(SIR10/NFIR10) ENDIF IF (NFIR21.EQ.0) THEN SIR21=0.0 ELSE SIR21=SQRT(SIR21/NFIR21) ENDIF IF (NFIR32.EQ.0) THEN SIR32=0.0 ELSE SIR32=SQRT(SIR32/NFIR32) ENDIF IF(NFIMW0 .EQ.0) THEN SMW0=0. ELSE SMW0=SQRT(SMW0/NFIMW0 ) ENDIF IF(NFIMW1 .EQ.0) THEN SMW1=0.0 ELSE SMW1=SQRT(SMW1/NFIMW1 ) ENDIF IF(NFIMW2.EQ.0) THEN SMW2=0. ELSE SMW2=SQRT(SMW2/NFIMW2) ENDIF IF(NFIMW3.EQ.0) THEN SMW3=0. ELSE SMW3=SQRT(SMW3/NFIMW3) ENDIF IF(NFIMW4.EQ.0) THEN SMW4=0. ELSE SMW4=SQRT(SMW4/NFIMW4) ENDIF IF(NFI004 .EQ.0) THEN SIG004=0. ELSE SIG004=SQRT(SIG004/NFI004 ) ENDIF IF(NFI070 .EQ.0) THEN SW070=0.0 ELSE SW070=SQRT(SW070/NFI070 ) ENDIF IF(NFI080 .EQ.0) THEN SIG080=0.0 ELSE SIG080=SQRT(SIG080/NFI080 ) ENDIF IF(NFI010 .EQ.0) THEN SIG010=0.0 ELSE SIG010=SQRT(SIG010/NFI010 ) ENDIF IF(NFI100 .EQ.0) THEN SIG100=0.0 ELSE SIG100=SQRT(SIG100/NFI100 ) ENDIF IF(NFI020 .EQ.0) THEN SIG020=0.0 ELSE SIG020=SQRT(SIG020/NFI020 ) ENDIF IF(NFI200 .EQ.0) THEN SIG200=0.0 ELSE SIG200=SQRT(SIG200/NFI200 ) ENDIF IF(NFI045 .EQ.0) THEN SIG045=0.0 ELSE SIG045=SQRT(SIG045/NFI045 ) ENDIF IF(NFI1M0 .EQ.0) THEN SIG1M0=0.0 ELSE SIG1M0=SQRT(SIG1M0/NFI1M0 ) ENDIF IF(NGST0.EQ.0) THEN SGST0=0. ELSE SGST0=SQRT(SGST0/NGST0) ENDIF PRINT*,' S= ',S,' NOFIT= ',NOFIT,'NPP= ',NPP,' SMW= ',SMW, & ' SIR= ',SIR,' NFITMW= ',NFITMW,' NFITIR= ',NFITIR,NFIMW1 , & ' SMW0= ',SMW0,' SMW1= ',SMW1,' NFITMW0= ',NFIMW0 , & ' NFITMW1= ',NFIMW1 ,' SIG004= ',SIG004,' NFIT004= ', & NFI004 ,' SIG080= ',SIG080,' NFIT080= ',NFI080 , & ' SIG010= ',SIG010,' NFIT010= ',NFI010 ,' SIG100= ',SIG100, & ' NFIT100= ',NFI100 ,' SIG020= ',SIG020,' NFIT020= ', & NFI020 ,' SIG200= ',SIG200,' NFIT200= ',NFI200 , & ' SIG045= ',SIG045,' NFIT045= ',NFI045 ,' SGST0= ',SGST0, & ' NGST0= ',NGST0,' SW070= ',SW070,' NFIT070= ',NFI070, & ' SMW2= ',SMW2,' NFIMW2= ',NFIMW2,' SIG1M0= ',SIG1M0, & ' NFIT1M0= ',NFI1M0,' SIG10= ',SIR10,' NFIT10= ',NFIR10, & ' SIG21= ',SIR21,' NFIT21= ',NFIR21,' SMW3= ',SMW3, & ' NFIMW3= ',NFIMW3,' SIR32= ',SIR32,' NFIR32=',NFIR32, & ' SMW4= ',SMW4,' NFIMW4= ', NFIMW4, & ' SMWA=',SWA,' NFITMWA=',NFITMWA,' SWME=',SWE, & ' NFITMWE=',NFITMWE WRITE(6,901) NOFIT 901 FORMAT(/,' NUMBER OF DATA INCLUDED IN THE FIT= ',I6) IF(IBGVTD .NE.0) THEN WRITE(6,51) CALL MATOUT(VTLDV,160,160,NPP,NPP) ENDIF 51 FORMAT(/,' VTLDV MATRIX BEFORE INVERSION',/) C C SET UP THE VTLDV MATRIX WITH ONE ON THE DIAGONAL. C DELX IS THE SQUARE ROOT OF THE DIAGONAL C DO 12 I=1,NPP DELX(I)= SQRT(VTLDV(I,I)) 12 CONTINUE DO 13 I=1,NPP DO 14 J=1,NPP VTLDV(I,J)=VTLDV(I,J)/(DELX(I)*DELX(J)) 14 CONTINUE 13 CONTINUE IF(IBGVTD .NE.0) THEN WRITE(6,51) CALL MATOUT(VTLDV,160,160,NPP,NPP) ENDIF c my insertion c singular value decomposition of the VTLDV IPATH=1 TOL=1.D-50 svd=0. vr=0. ul=0. call DLSVRR(NPP,NPP,VTLDV,160,IPATH,TOL,IRANK,SVD,UL,160,VR,160) write(6,*)'rank of the matrix',irank c my insertion C C INVERSION OF THE VTLDV MATRIX.FROM THIS POINT,VTLDV C BECOMES THE INVERSE MATRIX C cjth The next 8 lines or so needed fixing to run on the amur. c CALL F01AAF(VTLDV,160,NPP,UNIT,160,WKSPCE,IFAIL) c my insertion call DLINRG(NPP,VTLDV,160,VTLDV,160) c CALL F07ADF(NPP,NPP,VTLDV,160,IPIV,INFO) c CALL F07AJF(NPP,VTLDV,160,IPIV,WKSPCE,160,IFAIL) do 144 i=1,npp do 145 j=1,npp UNIT(i,j)=VTLDV(i,j) 145 continue 144 continue C CALL SGETRI(NPP,UNIT,160,WKSPCE,IFAIL) PRINT*,' IFAIL= ',IFAIL IF(IBGVTD .NE.0) THEN PRINT*,' UNIT AFTER INVERSION...' CALL MATOUT(UNIT,160,160,NPP,NPP) ENDIF DO 18 I=1,NPP DO 19 J=1,NPP VTLDV(I,J)=UNIT(I,J)/(DELX(I)*DELX(J)) 19 CONTINUE 18 CONTINUE IF(IBGVTD .NE.0) THEN PRINT*,' VTLDV AFTER INVERSION' CALL MATOUT(VTLDV,160,160,NPP,NPP) ENDIF C C FROM THIS POINT,DELX BECOMES THE VARIATION OF C THE PARAMETERS C DO 21 I=1,NPP DELX(I)=0.0 DO 800 K=1,NPP DELX(I)= VTLDV(I,K)*VTLDTM (K)+DELX(I) 800 CONTINUE 21 CONTINUE ACC=1. WRITE(6,*) 'ACC= ',ACC WRITE(6,302) C ***** Format changes by JTH 302 FORMAT(/,'STATE',2X,'CONSTANT',4X,'OLD VALUE',13X, 1 'VARIATION',13X,'NEW VALUE',13X,'STD.DEV.',12X,'VARIA/STV', 2 5X,'CONST/STV') DO 20 I=1,NPP DELTAX (I)=SQRT(VTLDV(I,I))*S IF(I.LE.NP12) THEN IF(I.GT.NP(1)) THEN II=2 J=I-NP(1) ELSE II=1 J=I ENDIF K=IND(II,J) XX=CST(II,K)+DELX(I)*ACC R= DELX(I)/DELTAX (I) RJTH=XX/DELTAX(I) WRITE(6,301) II,ABC(6*(K-1)+1:6*K),CST(II,K),DELX(I), & XX,DELTAX (I),R,RJTH CST(II,K)=XX 301 FORMAT(I4,2X,A6,4F22.15,F14.7,E14.2) ELSE II=2 VIBEN=TERM+DELX(I)*ACC R= DELX(I)/DELTAX (I) WRITE(6,301) II,' TERM ',TERM,DELX(I),VIBEN,DELTAX (I),R TERM=VIBEN ENDIF DELXST(I)=DELX(I) 20 CONTINUE C******NEW MATERIAL BELOW DO 322 I=1,NPP IF(I.LE.NP12) THEN IF(I.GT.NP(1)) THEN II=2 J=I-NP(1) ELSE II=1 J=I ENDIF K=IND(II,J) WRITE(6,323) ABC(6*(K-1)+1:6*K),CST(II,K) 323 FORMAT(2X,A6,' = ',F24.17,',') ENDIF 322 CONTINUE C******NEWMATERIAL ABOVE WRITE(6,900) S 900 FORMAT(/,' STANDARD DEVIATION BEFORE= ',F20.6) DELTAS=S-SST SST=S CALL CNTRL(NPP,DELX,DELTAX ,CST) c my insertion time CALL GETTIM (ihr, imin, isec, i100th) write(6,1012)ihr, imin, isec, i100th if(iend.eq.0)then c output of singular values and sigular right vectors write(6,*)' singular values and singular right vectors' icc=npp/8 if(icc*8.eq.npp)icc=icc-1 do icn=0,icc nppp=8+8*icn if(8+8*icn.gt.npp)nppp=npp write(6,*) ' ' write(6,314)(svd(ivr),ivr=1+8*icn,nppp) write(6,*) ' ' do isv=1,npp write(6,314)(VR(isv,ivr),ivr=1+8*icn,nppp) enddo enddo 314 format(1x,8E10.3) c covariation matrix do isv=1,npp do ivr=1,isv cov(isv,ivr)=0. do ikr=1,npp cov(isv,ivr)=cov(isv,ivr)+VR(isv,ikr)*VR(ivr,ikr)/(SVD(ikr)**2) enddo enddo enddo c correlation matrix do isv=1,npp do ivr=1,isv cor(isv,ivr)=cov(isv,ivr)/dsqrt(dabs(cov(isv,isv)*cov(ivr,ivr))) enddo enddo c output of the correlation matrix write(6,*)' correlation matrix ' do icn=0,icc write(6,*) ' ' do isv=1+8*icn,npp c decoding of parameter name IF(Isv.LE.NP12) THEN IF(Isv.GT.NP(1)) THEN II=2 J=Isv-NP(1) ELSE II=1 J=Isv ENDIF Ksv=IND(II,J) endif c end of decoding nppp=8+8*icn if(8+8*icn.gt.isv)nppp=isv write(6,315)isv,ABC(6*(Ksv-1)+1:6*Ksv), 1 (cor(isv,ivr),ivr=1+8*icn,nppp) enddo write(6,316)(ivr,ivr=1+8*icn,nppp) enddo 315 format(1x,I3,1x,A6,1x,8F9.4) 316 format(12x,8(6x,I3)) endif c my insertion C C CALL CHEKPTX(IFILES,1) C IF(IEND) 23,24,25 C24 SEC=TCPU() C WRITE(6,*) ' BEGIN SETUP LINECAL ',SEC 24 CALL SETUP(CST) C******THENEXT STATEMENT USED TO BE NUMBER 24. 244 CALL LINECA (CST,NPP) C SEC=TCPU() C WRITE(6,*) ' END LINECAL ',SEC WRITE(6,*) ' END LINECAL ' c my insertion time CALL GETTIM (ihr, imin, isec, i100th) write(6,1012)ihr, imin, isec, i100th c my insertion END C************************************************************** SUBROUTINE ASET(K,EGVL,EGVC,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************ C STORAGE OF THE COEFICIENTS OF THE TORSION EIGENVECTORS C CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION EGVL(NDMX),EGVC(NDMX,NDMX),ETOR(KDMX,9) DIMENSION A(9,NTORMX,KDMX) IK=K+N+1 DO 1 I=1,(09) DO 2 J=1,NDIMTO A(I,J,IK)=EGVC(J,I) 2 CONTINUE ETOR(IK,I)=EGVL(I) C PRINT*,'ETOR=',ETOR(IK,I),'IK=',IK,'I=',I 1 CONTINUE IF(K.EQ.-N) GO TO 16 DO 3 I=1,(09) IF((A(I,10,IK)/A(I,10,IK-1)).LT.0.0) THEN DO 4 J=1,NDIMTO A(I,J,IK)=-A(I,J,IK) 4 CONTINUE ENDIF 3 CONTINUE 16 RETURN END C************************************************************ SUBROUTINE CNTRL(NPP,DELX,DELTAX ,CST) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************ C MONITORS THE FITTING(CONCONVERGENCE,REJECTION OF TRANSITIONS, C DIVERGENCE) C CHARACTER*6 PARA,REF CHARACTER*1 PR,PROBS INTEGER REJ COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/PAR/IND(2,80),NP(2),X(160),ITO COMMON/EXPDAT / ETRANS(80000),IVOBS(2,80000), * NOBS(2,80000),KOBS(2,80000),IBLK(2,80000),IW(80000),W(80000), * NDATA,NADAT,EPSI,NOFIT,NFITMW,NFITIR,NFIMW0 ,NFIMW1 ,NITT, * IAST(80000),NFI004 ,NFI080 ,NFI100 ,NFI020 ,NFI200 ,NFI045 , & NGST0 ,NFI010 ,SW070,NFI070 ,NFIMW2,IKAKC,JMAX,NFI1M0, & NFIR10,NFIR21,NFIR32,NFIMW3,NFIMW4,SWA,SWE,NFITMWA,NFITMWE COMMON/CNTL/NIT,IEND,DELTAS,S,TME(80000),TERMST,NP12,SMW,SIR, 1 SMW0,SMW1,SIG004,SIG080,SIG10,SIG100,SIG020,SIG200,SIG045, & SGST0,SMW2,SIG1M0,SMW3,SIG32,SMW4 COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION DELX(160),DELTAX (160),CST(2,117) C THE CONVERGENCE CRITERION HAS BEEN CHANGED (NOV.1990). CC I = 1 IF(NIT.EQ.1) THEN IEND=1 WRITE(6,307) NIT NIT=NIT+1 RETURN ENDIF IF(NIT.EQ.NITT) THEN IEND=0 C WRITE(6,'("FINISHED",I4," PRESET ITERATIONS")') NIT RETURN ENDIF RS=DELTAS/S DO 6 I=1,NPP RI=ABS(DELX(I))/DELTAX (I) IF(RI.GT.1.E-3) THEN IF(RS.GT.1.E-1) THEN CC6 IF (ABS(DELX(I)) .GE. DELTAXL(I)) THEN C IF (DELTAS .GT. 0 .AND. NIT .GE. 4) THEN c WRITE(6,*) ' PROGRAM IS DIVERGING ...' c IF(NITT.LE.8) GO TO 54 c STOP c54 CONTINUE END IF CC IF(NIT.GT.NITT) STOP IEND = 1 WRITE(6,307) NIT NIT = NIT + 1 RETURN END IF CC I = I + 1 CC IF (I .LE. NPP) GO TO 6 C 6 CONTINUE REJ = 0 C CALL REJEC(TME,S,REJ) IF (REJ .EQ. 0) THEN WRITE(6,*) ' PROGRAM HAS CONVERGED AFTER ',NIT, 1 ' ITERATIONS' IEND = 0 RETURN ELSE DO 94 I = 1,NP12 IF (I .GT. NP(1)) THEN II = 2 J = I - NP(1) ELSE II = 1 J = I END IF CST(II,IND(II,J)) = X(I) 94 CONTINUE TERM = TERMST IEND = -1 WRITE(6,807) RETURN END IF 307 FORMAT(/,' *** END OF CYCLE NO. ',I3,' *** ') 807 FORMAT(//,' *** NEW LEAST SQUARE FITTING PROCEDURE ***', 1 /,' ==========================================') END SUBROUTINE DATA1(CST) use assig IMPLICIT REAL*8 ( A-H,O-Z ) C********************************************************** C THIS SUBROUTINE READS THE EXPERIMENTAL VALUES AND CALLS 2 C SUBROUTINES:DDATA11 WHICH READS THE VALUES OF THE CONSTANTS C AND PARAM WHICH READS WHICH PARAMETERS MUST BE FITTED C*********************************************************** C READ THE CONSTANT VALUES AND DECIDE WHICH PARAMETERS TO FIT C I=1 LOWER STATE;I=2 UPPER STATE C CHARACTER ABC*702 CHARACTER*1 PR,PROBS,PRINF,PRSUP CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/PAR/IND(2,80),NP(2),X(160),ITO COMMON/EXPDAT / ETRANS(80000),IVOBS(2,80000), * NOBS(2,80000),KOBS(2,80000),IBLK(2,80000),IW(80000),W(80000), * NDATA,NADAT,EPSI,NOFIT,NFITMW,NFITIR,NFIMW0 ,NFIMW1 ,NITT, * IAST(80000),NFI004 ,NFI080 ,NFI100 ,NFI020 ,NFI200 ,NFI045 & ,NGST0 ,NFI010 ,SW070,NFI070 ,NFIMW2,IKAKC,JMAX,NFI1M0, & NFIR10,NFIR21,NFIR32,NFIMW3,NFIMW4,SWA,SWE,NFITMWA,NFITMWE COMMON/SYMB/ABC COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117) c my insertion open(unit=5,file='input.txt',access='sequential',status='old') c the following cards are for Ilyushyn label scheme c my insertion c open(unit=22,file='Evt01.dat',access='sequential',status='old') c iklflag=0 c ikl=0 c do while (iklflag.ne.1) c ikl=ikl+1 c read(22,*,end=1001)ivt1,j1,kkk1,neig1 c eseign(ivt1,j1,kkk1)=neig1 c cycle c1001 iklflag=1 c exit c enddo c ipoint=ikl-1 c close(22) c open(unit=23,file='Avt01.dat',access='sequential',status='old') c iklflag=0 c ikl=0 c do while (iklflag.ne.1) c ikl=ikl+1 c read(23,*,end=1002)ivt1,j1,kkk1,ipar1,neig1 c if(ipar1.eq.2)then c aseign(ivt1,j1,0,kkk1)=neig1 c else c aseign(ivt1,j1,1,kkk1)=neig1 c endif c cycle c1002 iklflag=1 c exit c enddo c ipoint1=ikl-1 c close(23) WRITE(6,*) ' THIS OUTPUT IS FROM BELGI38 ' READ(5,*) EPSI,IBUG1,IBUG2,IBUG3,IBGTME,IBGVTD ,IBUG4 WRITE(6,*) ' EPSI= ',EPSI,' IBUG1= ',IBUG1,' IBUG2= ', * IBUG2,' IBUG3= ',IBUG3,' IBGTME= ',IBGTME,' IBGVTDV= ', * IBGVTD ,'IBUG4= ',IBUG4 READ(5,20) KTRONC READ(5,20) IRMW 20 FORMAT(I3) WRITE(6,21)KTRONC,IRMW 21 FORMAT(' KTRONC = ',I3,' IRMW= ',I3) IF(KTRONC.GT.(10)) THEN WRITE(6,*) ' DIMENSION KTRONC TOO LOW ' STOP ENDIF C THE NEXT CARD WAS ADDED FOR NBS COMPUTER. EVIB(1)=0. EVIB(2)=0. IF(IRMW.EQ.1) THEN DO 2 I=1,1 CALL DATA11(I,CST) READ (5,*) EVIB(I) CALL PARAM(I,CST) C PRINT*,' NP APRES PARAM= ',NP(I) 2 CONTINUE ELSEIF(IRMW.EQ.2) THEN DO 3 I=1,2 CALL DATA11(I,CST) C PRINT*, 'I= ',I,' IRMW= ',IRMW READ (5,*) EVIB(I) WRITE(6,*) ' EVIB= ',EVIB(I),'I= ',I PRINT*,' EVIB(1)= ',EVIB(1),' EVIB(2)= ',EVIB(2) CALL PARAM(I,CST) C PRINT*,' NP APRES PARAM= ',NP(I) 3 CONTINUE ENDIF READ(5,*) ITO,NITT,IKAKC,JMAX WRITE(6,*) ' VIB.ENERGY TO BE FITTED? ITO= ',ITO WRITE(6,*)' NUMBER OF ITERATIONS= ',NITT WRITE(6,*)' MAXIMUM VALUE OF J= ',JMAX C I=1 C C IBLK:DNUMBER OF THE N BLOCK;NBR:DPLACE IN THIS BLOCK C 1 IF(IRMW.EQ.1) THEN cjth I have changed these reads to follow des conventions internationales! c READ(5,23,END=10) ETRANS(I),IVOBS(1,I),NOBS(1,I), READ(5,23,END=10) ETRANS(I),IVOBS(2,I),NOBS(2,I), c A KOBS(1,I),PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I), A KOBS(2,I),PROBS(2,I),IBLK(2,I),IVOBS(1,I),NOBS(1,I), c B KOBS(2,I),PROBS(2,I),IBLK(2,I),IAST(I),IW(I),REF(I) B KOBS(1,I),PROBS(1,I),IBLK(1,I),IAST(I),IW(I),REF(I) cjuan This command must be in to generate Juan's table. c iast(i)=1 C cjth I have added the following card in several places to eliminate vt=1. c if(ivobs(1,i).eq.1.and.ivobs(2,i).eq.1) go to 1 cjth I added these cards to make it unnecessary to input the IBLK values. c IBLK are not necessary anymore! iblk(1,I)=nobs(1,I)+1 iblk(2,I)=nobs(2,I)+1 c IJK debug 9 octobre 2006, change of sign of Pg-rho.Pa if(probs(1,I).eq.' ') then KOBS(1,I)=-KOBS(1,I) KOBS(2,I)=-KOBS(2,I) endif c cjth I added these cards to calculate an entire vt=1 spectrum c They must be taken out immediately after the calculation!!!!!! c ivobs(1,I) = 1 c ivobs(2,I) = 1 c IF(ETRANS(I).EQ.10000.OR.ETRANS(I).EQ.10000.000) THEN REF(I)=' ' ENDIF c write(6,23) ETRANS(I),IVOBS(1,I),NOBS(1,I), c A KOBS(1,I),PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I), c B KOBS(2,I),PROBS(2,I),IBLK(2,I),IAST(I),IW(I),REF(I) c those lines are only for CH3COOH! cjth I have turned off this vt=1 quantum number exchange (which were used c for acetic acid)! IF(IVOBS(1,I).EQ.9.AND.IAST(I).EQ.1) THEN c IF(IVOBS(1,I).EQ.1.AND.IAST(I).EQ.1) THEN NSUP=NOBS(2,I) NINF=NOBS(1,I) KSUP=KOBS(2,I) KINF=KOBS(1,I) PRSUP=PROBS(2,I) PRINF=PROBS(1,I) NOBS(2,I)=NINF NOBS(1,I)=NSUP KOBS(2,I)=KINF KOBS(1,I)=KSUP PROBS(2,I)=PRINF PROBS(1,I)=PRSUP ENDIF c IF(IVOBS(1,I).EQ.1.AND.NOBS(1,I).GT.17) IAST(I)=0 c cjth***THE NEXT WRITE IS USEFUL FOR DEBUGGING TYPOS IN THE DATA FILE. JTH c WRITE(6,24) ETRANS(I),IVOBS(1,I),NOBS(1,I),KOBS(1,I) c A ,PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I),KOBS(2,I) c B ,PROBS(2,I),IBLK(2,I),IAST(I),I,IW(I),REF(I) C THE WEIGHT (W) IS EQUAL TO THE INVERSE OF THE SQUARE OF C EXPERIMENTAL UNCERTAINTY IF(IW(I).EQ.1005) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=3.59502E+13 ELSEIF(IW(I).EQ.1000) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=8.98755E+12 ELSEIF(IW(I).EQ.1010) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=8.98755E+12 ELSEIF(IW(I).EQ.1040) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=5.61722E+11 ELSEIF(IW(I).EQ.1004) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=5.61722E+13 ELSEIF(IW(I).EQ.1003) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=9.986168E+13 ELSEIF(IW(I).EQ.1020) THEN ETRANS(I)=ETRANS(I)/29979.2458 W (I)=2.24688E+12 ELSEIF(IW(I).EQ.1050) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=3.59502E+11 ELSEIF(IW(I).EQ.1030) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=9.98617E+11 ELSEIF(IW(I).EQ.1008) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=1.4043E+13 ELSEIF(IW(I).EQ.1100) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=8.98755E+10 ELSEIF(IW(I).EQ.1990) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=8.98755E+8 ELSEIF(IW(I).EQ.1180) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=2.77393E+10 ELSEIF(IW(I).EQ.1130) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=5.31807E+10 ELSEIF(IW(I).EQ.1150) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=3.51076E+10 ELSEIF(IW(I).EQ.1160) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=3.99446E+10 ELSEIF(IW(I).EQ.1200) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=2.24688E+10 ELSEIF(IW(I).EQ.1080) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=1.4043E+11 ELSEIF(IW(I).EQ.1070) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=1.83419E+11 ELSEIF(IW(I).EQ.2000) THEN ETRANS(I)=ETRANS(I)/29979.2458 W(I)=2.24688E+8 ELSEIF(IW(I).EQ.1000.AND.IVOBS(1,I).EQ.1) THEN W(I)=5.61722E+11 ETRANS(I)=ETRANS(I)/29979.2458 ELSE IF(IW(I).EQ.10) THEN W(I)=4000000.0 ELSE IF(IW(I).EQ.2.OR.IW(I).EQ.3) THEN W(I)=4000000.0 ELSE W(I)=8E+6 ENDIF IF(ETRANS(I).EQ.-900.0000) THEN NADAT=I-1 GO TO 1 ENDIF I=I+1 GO TO 1 10 NDATA=I-1 PRINT*,NDATA,NADAT ELSE W(I)=1. C those format from fitk5jb| cjth I have changed these reads to follow standard conventions! c READ(5,225,END=11) ETRANS(I),IVOBS(1,I),NOBS(1,I), READ(5,225,END=11) ETRANS(I),IVOBS(2,I),NOBS(2,I), c A KOBS(1,I),PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I), A KOBS(2,I),PROBS(2,I),IBLK(2,I),IVOBS(1,I),NOBS(1,I), c B KOBS(2,I),PROBS(2,I),IBLK(2,I),IAST(I),IW(I) B KOBS(1,I),PROBS(1,I),IBLK(1,I),IAST(I),IW(I) cjuan This command must be in to generate Juan's table. c iast(i)=1 C cjth I have added the following card in several places to eliminate vt=1. c if(ivobs(1,i).eq.1.and.ivobs(2,i).eq.1) go to 1 cjth I added these cards to make it unnecessary to input the IBLK values. iblk(1,I)=nobs(1,I)+1 iblk(2,I)=nobs(2,I)+1 c IJK debug 9 octobre 2006, change of sign of Pg-rho.Pa if(probs(1,I).eq.' ') then KOBS(1,I)=-KOBS(1,I) KOBS(2,I)=-KOBS(2,I) endif c C IF(IVOBS(2,I).EQ.3) IVOBS(2,I)=4 C IF(IVOBS(1,I).EQ.3) IVOBS(1,I)=4 WRITE(6,226) ETRANS(I),IVOBS(1,I),NOBS(1,I),KOBS(1,I) A,PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I),KOBS(2,I) B,PROBS(2,I),IBLK(2,I),IAST(I),I,IW(I) 225 FORMAT (F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1I3,1X,A1,1X,I3,1X,I3,1X,I4) 226 FORMAT (1X,F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1I3,1X,A1,1X,I3,1X,I3,1X,I4,1X,I4) C the following cards to be turned in| C READ(5,23,END=11) ETRANS(I),IVOBS(1,I),NOBS(1,I), C A KOBS(1,I),PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I), C B KOBS(2,I),PROBS(2,I),IBLK(2,I),IAST(I),IW(I),REF(I) C WRITE(6,26) ETRANS(I),IVOBS(1,I),NOBS(1,I),KOBS(1,I) C A ,PROBS(1,I),IBLK(1,I),IVOBS(2,I),NOBS(2,I),KOBS(2,I) C B ,PROBS(2,I),IBLK(2,I),IAST(I),I,IW(I),REF(I) C23 FORMAT (F12.5,I2,2X,I2,2X,I3,2X,A1,1X,I2,2X,I1,2X,I2,2X, C & I3,2X,A1,1X,I2,2X,I1,3X,I4,3X,A1) C25 FORMAT (F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, C 1 I3,1X,A1,1X,I3,1X,I3,1X,I4,3X,A1) IF(ETRANS(I).EQ.-900.)THEN NADAT=I-1 GO TO 1 ENDIF I=I+1 GO TO 1 11 NDATA=I-1 PRINT*,NDATA,NADAT ENDIF C C NADAT:DNBREOF DATA A TYPE(SIGMA=0);NDATA:DNBRE OF TOTAL DATA C IF(NDATA.GT.(80000))THEN WRITE(6,*) ' TOO MANY DATA...' STOP ENDIF C3 FORMAT (F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, C 1I3,1X,A1,1X,I3,1X,I3,1X,I4) 23 FORMAT (F12.5,I2,2X,I2,2X,I3,2X,A1,1X,I2,2X,I1,2X,I2,2X, & I3,2X,A1,1X,I2,2X,I1,3X,I4,3X,A6) 25 FORMAT (F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1 I3,1X,A1,1X,I3,1X,I3,1X,I4,3X,A6) 26 FORMAT (1X,F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1 I3,1X,A1,1X,I3,1X,I3,1X,I4,1X,I4,3X,A6) 24 FORMAT (1X,F12.5,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1 I3,1X,A1,1X,I3,1X,I3,1X,I4,1X,I4,3X,A6) IF(IRMW.EQ.1) NP(2)=0 cjth I turned on this debugging statement. I turned it back off. c PRINT*, 'NP(1) A LA FIN DE DATA ',NP(1),'NP(2)= ',NP(2) RETURN END C************************************************************** SUBROUTINE DATA11(LO,CST) IMPLICIT REAL*8 ( A-H,O-Z ) C C READING OF THE CONSTANTS FOR EACH STATE AND PUTTING THEM IN C CST VECTOR C CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*4 ANSWER CHARACTER*6 PARA,REF COMMON/SYMB/ABC COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM DIMENSION CST(2,117) NAMELIST/CSTE/OA,B,C,DJ,DJK,DK,ODELN,ODELK,HJ,HJK 1 ,HKJ,HK,OHJ,OHJK,OHK,FV,OFV,GV,ALV,AK1,AK2,AK3, 2 AK4,AK5,AK6,AK7,DAB,ODAB,DABJ,DABK,OLV,C1,C2,C3,HLKJ,AK5J, & C2J, 3 GVJ,AK2J,C1J,ANV,BK2,AMV,DAC,DBC,C4,ODELTA,DELTA,AK1J,AK5K, 4 AK2K,AK1K,AK3J,AK3K,BK1,AK7J,C11,ODAB6,AK6K,AK7K,AK66,AK6J, 5 AK3B,AK4B,DELTAB,ODELTB,C4J,C4K,C1K,C2K,C12, 6 ANVJ,AMVJ,AK3JJ,BK2J,BK1J,AK3KJ,BK1K,AK3KK,C3J,C12J,AK2JJ, 7 AK1JJ,AK2JK,AK1JK,AK3BJ,AK4BJ,AK3BB,AK4BB,ESPOIR,AK9,V12J, 8 V12K,DAC12,DBC12,DACJ,DBCJ,ODABJ,RHOB,RHOC, 9 F,RHO,V3,V6,V9,V12,RHORHO, ODABK,DACK,V9J,V9K,C11J,C11K, ! DABJJ,DABJK,ODABJJ,ODABJK C * Modification Denis Girou (CIRCE) * ABC = 'OA B C DJ DJK DK ODELN ODELK HJ HJK * 1 HKJ HK OHJ OHJK OHK FV OFV GV ALV AK1 AK2 * 2 AK3 AK4 AK5 AK6 AK7 DAB ODAB DABJ DABK OLV C1 * 3 C2 C3 HLKJ AK5J C2J GVJ AK2J C1J ANV BK2 AMV * 4 DAC DBC C4 ODELTADELTA AK1J AK5K AK2K AK1K AK3J AK3K * 5 F RHO V3 V6 V9 RHORHO' ABC(1:210)= & 'OA B C DJ DJK DK ODELN ODELK HJ HJK 1HKJ HK OHJ OHJK OHK FV OFV GV ALV AK1 AK2 2AK3 AK4 AK5 AK6 AK7 DAB ODAB DABJ DABK OLV C1 3C2 C3 HLKJ ' ABC(211:702)= 'AK5J C2J GVJ AK2J C1J ANV BK2 AMV 4DAC DBC C4 ODELTADELTA AK1J AK5K AK2K AK1K AK3J AK3K 5BK1 AK7J C11 ODAB6 AK6K AK7K AK66 AK6J AK3B AK4B DELTAB 6ODELTBC4J C4K C1K C2K C12 ANVJ AMVJ AK3JJ BK2J BK1J 7AK3KJ BK1K AK3KK C3J C12J AK2JJ AK1JJ AK2JK AK1JK AK3BJ AK4BJ 8AK3BB AK4BB ESPOIRAK9 V12J V12K DAC12 DBC12 DACJ DBCJ ODABJ 9RHOB RHOC F RHO V3 V6 V9 V12 RHORHOODABK DACK &V9J V9K C11J C11K DABJJ DABJK ODABJJODABJK' C INITIALISATION C OA=0.0 B=0.0 C=0.0 DJ=0.0 DJK=0.0 DK=0.0 ODELN=0.0 ODELK=0.0 HJ=0.0 HJK=0.0 HKJ=0.0 HK=0.0 OHJ=0.0 OHJK=0.0 OHK=0.0 FV=0.0 OFV=0.0 GV=0.0 ALV=0.0 AK1=0.0 AK2=0.0 AK3=0.0 AK4=0.0 AK5=0.0 AK6=0.0 AK7=0.0 DAB=0.0 DBC=0.0 C4=0.0 ODELTA=0.0 AK1J=0.0 AK5K=0.0 AK2K=0.0 AK1K=0.0 AK3J=0.0 AK3K=0.0 DELTA=0.0 ODAB=0.0 DABJ=0.0 DABK=0.0 OLV=0.0 C1=0.0 C2=0.0 C3=0.0 HLKJ=0.0 AK5J=0.0 C2J=0.0 GVJ=0.0 AK2J=0.0 C1J=0.0 ANV=0.0 BK2=0.0 AMV=0.0 DAC=0.0 DBC12=0.0 DAC12=0.0 BK1=0.0 AK7J=0.0 C11=0.0 ODAB6=0.0 AK6K=0.0 AK7K=0.0 AK66=0.0 AK6J=0.0 AK3B=0.0 AK4B=0.0 DELTAB=0.0 ODELTB=0.0 C4J=0.0 C4K=0.0 C1K=0.0 C2K=0.0 C12=0.0 ANVJ=0.0 AMVJ=0.0 AK3JJ=0.0 BK2J=0.0 BK1J=0.0 AK3KJ=0.0 F=0.0 RHO=0.0 V3=0.0 V6=0.0 V9=0.0 V12=0.0 V12J=0.0 V12K=0.0 BK1K=0.0 AK3KK=0.0 C3J=0.0 C12J=0.0 AK2JJ=0.0 AK1JJ=0.0 AK2JK=0.0 AK1JK=0.0 AK3BJ=0.0 AK4BJ=0.0 AK3BB=0.0 AK4BB=0.0 ESPOIR=0.0 DACJ=0.0 DBCJ=0.0 ODABJ=0.0 RHOB=0.0 RHOC=0.0 AK9=0.0 C RHOJ=0.0 C RHOK=0.0 RHORHO=0.0 ODABK=0.0 DACK=0.0 V9J=0.0 V9K=0.0 C11J=0. C11K=0. DABJJ=0. DABJK=0. ODABJJ=0. ODABJK=0. C C INITIALISATION C DO 3 II=1,(117) CST(LO,II)=0.0 3 CONTINUE C ASK THE USER FOR THE CONSTANT VALUES WRITE(6,132) 132 FORMAT(' CONSTANTS IN CM-1 OR IN MHZ? ',':DA') READ(5,'(A4)') ANSWER WRITE(6,'(A4)') ANSWER CONV=1. IF(ANSWER.EQ.'MHZ') CONV=29979.2458 C READ FROM NAMELIST READ(5,CSTE) WRITE(6,CSTE) C PUT THE CONSTANT VALUES IN CST VECTOR c my insertion c OPEN (35, FILE = 'dat', ACCESS = 'SEQUENTIAL', STATUS = 'NEW') c my insertion OPEN (35, FILE = 'dat', ACCESS = 'SEQUENTIAL', STATUS = 'unknown', ! form='unformatted') WRITE(35) OA,B,C,DJ,DJK,DK,ODELN,ODELK,HJ,HJK,HKJ,HK, 1 OHJ,OHJK,OHK,FV,OFV,GV,ALV,AK1,AK2,AK3,AK4,AK5,AK6,AK7, 2 DAB,ODAB,DABJ,DABK,OLV,C1,C2,C3,HLKJ,AK5J,C2J,GVJ,AK2J,C1J 3 ,ANV,BK2,AMV,DAC,DBC,C4,ODELTA,DELTA,AK1J,AK5K,AK2K,AK1K,AK3J 4 ,AK3K,BK1,AK7J,C11,ODAB6,AK6K,AK7K,AK66,AK6J,AK3B,AK4B,DELTAB 5 ,ODELTB,C4J,C4K,C1K,C2K,C12,ANVJ,AMVJ,AK3JJ,BK2J,BK1J,AK3KJ 6 ,BK1K,AK3KK,C3J,C12J,AK2JJ,AK1JJ,AK2JK,AK1JK,AK3BJ,AK4BJ 7 ,AK3BB,AK4BB,ESPOIR,AK9,V12J,V12K,DAC12,DBC12,DACJ,DBCJ 8 ,ODABJ,RHOB,RHOC 9 ,F,RHO,V3,V6,V9,V12,RHORHO,ODABK,DACK,V9J,V9K,C11J,C11K ! ,DABJJ,DABJK,ODABJJ,ODABJK REWIND 35 READ(35) (CST(LO,II),II=1,(117)) REWIND 35 DO 18 III=1,100 CST(LO,III)=CST(LO,III)/CONV 18 CONTINUE CST(LO,108)=CST(LO,108)/CONV CST(LO,109)=CST(LO,109)/CONV CST(LO,110)=CST(LO,110)/CONV CST(LO,111)=CST(LO,111)/CONV CST(LO,112)=CST(LO,112)/CONV CST(LO,113)=CST(LO,113)/CONV CST(LO,114)=CST(LO,114)/CONV CST(LO,115)=CST(LO,115)/CONV CST(LO,116)=CST(LO,116)/CONV CST(LO,117)=CST(LO,117)/CONV C PRINT*,CST(LO,38) C PRINT*,RHO,CST(LO,1) RETURN END C*********************************************************** SUBROUTINE ENCAL(E,LO,EGVL,EGVC,A) use assig IMPLICIT REAL*8 ( A-H,O-Z ) C*********************************************************** CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/ROTOR/NROTOR,NDIMTO COMMON/QUANT/ISIG,IV,N,KK COMMON/QUANT2/KKK(NDMX),NVT(NDMX),IPRITY(NDMX) COMMON/C2C/PERC2L,PER23L(NDMX) COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION EGVL(NDMX),EGVC(NDMX,NDMX),A(9,NTORMX,KDMX) DIMENSION EGVC4(KDMX,NTORMX,NDMX) C TAKE THE GOOD EIGENVALUE(E) AND ASSOCIATE A NUMBER(NBR) C TO THE EIGENVECTOR CORRESPONDING TO THE K DATA c call IPOSA or IPOSE or call eseign or aseign nbr(LO)=1 C IF(ISIG.EQ.0) THEN DO 5 J=1,NROTOR IVTORP=1 DO 4 I=1,NROTOR IU=(2*N+1)*IVTORP+1 IF(I.EQ.IU) IVTORP=IVTORP+1 II=I-(2*N+1)*(IVTORP-1) EGVC4(II,IVTORP,J)=EGVC(I,J) C IF(N.EQ.1.AND.J.LE.12) THEN C PRINT*,'EGVC4=',EGVC4(II,IVTORP,J),'II=',II,'J=',J C 1,'IVTORP=',IVTORP,'EGVC=',EGVC(I,J),'I=',I C ENDIF 4 CONTINUE 5 CONTINUE C HERE WE ADDED ANOTHER FUNCTION : IN IPOSA ALL THE EIGENVALUES C ARE LABELLED, IN IPOSA2 WE ONLY HAVE TO GET THE RIGHT ONE. C THE WIGNER SUBROUTINE IS DOING THE ROTATION OF THE EIGENVECT. C SO WE PASS FROM THE RAM TO THE PAM AXIS SYTEM. IF (PR.EQ.'+') IPR=1 IF (PR.EQ.'-') IPR=0 c NBR(LO)=aseign(iv,n,ipr,kk) IF(KK.EQ.0.AND.IV.EQ.0)THEN NBR(LO)=IPOSA(EGVC,EGVL,A,EGVC4) ELSE NBR(LO)=IPOSA2(EGVC,EGVL,A,EGVC4,PR) ENDIF ELSE c NBR(LO)=eseign(iv,n,kk) IF(KK.EQ.(-N).AND.IV.EQ.0) THEN NBR(LO)=IPOSE(EGVC,EGVL) ELSE NBR(LO)=IPOSE2(EGVC,EGVL) ENDIF ENDIF if(nbr(LO).gt.ndmx)nbr(LO)=ndmx if(nbr(LO).lt.1)nbr(LO)=1 C PRINT*,' NBR= ',NBR(LO),' LO= ',LO E=EGVL(NBR(LO)) IVTORP=1 DO 3 I=1,NROTOR IU=(2*N+1)*IVTORP+1 IF(I.EQ.IU) IVTORP=IVTORP+1 II=I-(2*N+1)*(IVTORP-1) EGVC3(II,IVTORP)=EGVC(I,NBR(LO)) IF(IVTORP.GE.9) EGVC3(II,IVTORP)=0.0 C PRINT*,' EGVC(I,NBR(LO))= ',EGVC(I,NBR(LO)) C PRINT*,' EGVC3= ',EGVC3(II,IVTORP),' IVTORP= ',IVTORP 3 CONTINUE RETURN END CC********************************************************** FUNCTION H0 (K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C*********************************************************** C SETS UP THE DIAGONAL ELEMENTS OF THE HAMILTONIAN CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK DIMENSION CST(2,117),A(9,NTORMX,KDMX),ETOR(KDMX,9) OA=CST(LO,1) B=CST(LO,2) C=CST(LO,3) DJ=CST(LO,4) DJK=CST(LO,5) DK=CST(LO,6) ODELN=CST(LO,7) ODELK=CST(LO,8) HJ=CST(LO,9) HJK=CST(LO,10) HKJ=CST(LO,11) HK=CST(LO,12) OHJ=CST(LO,13) OHJK=CST(LO,14) OHK=CST(LO,15) FV=CST(LO,16) OFV=CST(LO,17) GV=CST(LO,18) ALV=CST(LO,19) AK1=CST(LO,20) AK2=CST(LO,21) AK3=CST(LO,22) AK4=CST(LO,23) AK5=CST(LO,24) AK6=CST(LO,25) AK7=CST(LO,26) DAB=CST(LO,27) ODAB=CST(LO,28) ODABJ=CST(LO,98) OLV=CST(LO,31) C1=CST(LO,32) C2=CST(LO,33) C3=CST(LO,34) HLKJ=CST(LO,35) AK5J=CST(LO,36) C2J=CST(LO,37) GVJ=CST(LO,38) AK2J=CST(LO,39) C1J=CST(LO,40) ANV=CST(LO,41) BK2=CST(LO,42) AMV=CST(LO,43) AK1J=CST(LO,49) AK5K=CST(LO,50) AK2K=CST(LO,51) AK1K=CST(LO,52) AK3J=CST(LO,53) AK3K=CST(LO,54) BK1=CST(LO,55) AK7J=CST(LO,56) AK6K=CST(LO,59) AK7K=CST(LO,60) AK66=CST(LO,61) AK6J=CST(LO,62) AK3B=CST(LO,63) AK4B=CST(LO,64) ANVJ=CST(LO,72) AMVJ=CST(LO,73) AK3JJ=CST(LO,74) BK2J=CST(LO,75) BK1J=CST(LO,76) AK3KJ=CST(LO,77) BK1K=CST(LO,78) AK3KK=CST(LO,79) C3J=CST(LO,80) C12J=CST(LO,81) AK2JJ=CST(LO,82) AK1JJ=CST(LO,83) AK2JK=CST(LO,84) AK1JK=CST(LO,85) AK3BJ=CST(LO,86) AK4BJ=CST(LO,87) AK3BB=CST(LO,88) AK4BB=CST(LO,89) AK9=CST(LO,91) V12J=CST(LO,92) V12K=CST(LO,93) V9J=CST(LO,110) V9K=CST(LO,111) RHO=CST(LO,102) RHORHO=CST(LO,107) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) IK=K+N +1 C INITIALISATION A1=0.0 A2=0.0 A3=0.0 A4=0.0 A4B=0.0 A4BB=0.0 A5BB=0.0 A5=0.0 A5B=0.0 A6=0.0 A7=0.0 A7=0.0 A8=0.0 A8B=0.0 A9B=0.0 A14B=0.0 A15B=0.0 A9=0.0 A10=0.0 A11=0.0 A12=0.0 A13=0.0 A14=0.0 A15=0.0 A16=0.0 A17=0.0 A10B=0.0 A11B=0.0 A16B=0.0 A17B=0.0 A18=0.0 A19=0.0 A20=0.0 A21=0.0 DO 9 II=1,NDIMTO KL=II-KTRON1 COEFA= A(IVTOR,II,IK) COFAAA =A(IVTORP,II,IK) IF (II.EQ.NDIMTO ) THEN COEFAA=0.0 ELSE COEFAA=A(IVTORP,II+1,IK) ENDIF IF(II.EQ.1) THEN COEFA4=0.0 ELSE COEFA4=A(IVTORP,II-1,IK) ENDIF IF(II.GT.(NDIMTO -2)) THEN COEFA5=0.0 ELSE COEFA5=A(IVTORP,II+2,IK) ENDIF IF(II.LE.2) THEN COEFA6=0.0 ELSE COEFA6=A(IVTORP,II-2,IK) ENDIF IF(II.GT.(NDIMTO -4)) THEN COEFA7=0.0 ELSE COEFA7=A(IVTORP,II+4,IK) ENDIF IF(II.LE.4) THEN COEFA8=0.0 ELSE COEFA8=A(IVTORP,II-4,IK) ENDIF IF(II.GT.(NDIMTO -3)) THEN COEFA9=0.0 ELSE COEFA9=A(IVTORP,II+3,IK) ENDIF IF(II.LE.3) THEN COEFA10=0.0 ELSE COEFA10=A(IVTORP,II-3,IK) ENDIF C RHOEFF=RHO+RHOJ*N*(N+1)+RHOK*K**2 RHOEFF=RHO RHOKSG =3.*KL+ISIG+RHOEFF*K RHOKK=3*(KL+1)+ISIG+RHOEFF*K RHOKK2=3*(KL+2)+ISIG+RHOEFF*K RHOKM1=3.*(KL-1)+ISIG+RHOEFF*K RHOKM2=3.*(KL-2)+ISIG+RHOEFF*K A1=A1+COEFA*COEFAA A2=A2+(COEFA*COFAAA )*(RHOKSG **2) A3=A3+COEFA*COFAAA *RHOKSG A4=A4+COEFA*COFAAA *RHOKSG **3 A4B=A4B+COEFA*COFAAA *RHOKSG **5 A4BB=A4BB+COEFA*COFAAA *RHOKSG **7 A5=A5+COEFA*COFAAA *RHOKSG **4 A5B=A5B+COEFA*COFAAA *RHOKSG **6 A5BB=A5BB+COEFA*COFAAA *RHOKSG **8 A6=A6+COEFA*COFAAA A7=A7+COEFA*COEFA4 A8=A8+COEFA*COEFA4*RHOKSG A8B=A8B+COEFA*COEFA6*RHOKSG A9=A9+COEFA*COEFAA*RHOKSG A9B=A9B+COEFA*COEFA5*RHOKSG A10=A10+COEFA*COEFA4*RHOKSG **2 A10B=A10B+COEFA*COEFA6*RHOKSG **2 A11=A11+COEFA*COEFAA*RHOKSG **2 A11B=A11B+COEFA*COEFA5*RHOKSG **2 A12=A12+COEFA*COEFA5 A13=A13+COEFA*COEFA6 A18=A18+COEFA*COEFA7 A19=A19+COEFA*COEFA8 A20=A20+COEFA*COEFA9 A21=A21+COEFA*COEFA10 A14=A14+COEFA*COEFAA*RHOKK A14B=A14B+COEFA*COEFA5*RHOKK2 A15=A15+COEFA*COEFA4*RHOKM1 A15B=A15B+COEFA*COEFA6*RHOKM2 A16=A16+COEFA*COEFAA*RHOKK**2 A16B=A16B+COEFA*COEFA5*RHOKK2**2 A17=A17+COEFA*COEFA4*RHOKM1**2 A17B=A17B+COEFA*COEFA6*RHOKM2**2 9 CONTINUE C PRINT*,' A1= ',A1,' A2= ',A2,' A3= ',A3,' A4= ',A4,' A5= ', C 1A5 H0=(0.5*(B+C)*(N*(N+1)-K**2)+OA*K**2)*A6 H0=H0+(FV+OFV*N*(N+1))*(A6-0.5*A7-0.5*A1)*N*(N+1) H0=H0+(ANV*N*(N+1))*(A6-0.5*A13-0.5*A12) H0=H0+(V12J*N*(N+1))*(A6-0.5*A18-0.5*A19) H0=H0+(V12K*K**2)*(A6-0.5*A18-0.5*A19) H0=H0+(V9J*N*(N+1))*(A6-0.5*A20-0.5*A21) H0=H0+(V9K*K**2)*(A6-0.5*A20-0.5*A21) H0=H0+(ANVJ*N*(N+1)*N*(N+1))*(A6-0.5*A13-0.5*A12) H0=H0+(A6-0.5*A12-0.5*A13)*(K**2)*BK2 H0=H0+(A6-0.5*A12-0.5*A13)*(K**2)*BK2J*N*(N+1) H0=H0+GV*A2*N*(N+1) H0=H0+GVJ*A2*N**2*(N+1)**2 H0=H0+ALV*K*A3*N*(N+1) H0=H0-K**2*DJK*N*(N+1)*A6 H0=H0-N*(N+1)*DJ*N*(N+1)*A6 H0=H0+AK1*K**3*A3 H0=H0+AK1K*K**5*A3 H0=H0+AK1JK*K**5*A3*N*(N+1) H0=H0+AK1J*K**3*A3*N*(N+1) H0=H0+AK1JJ*K**3*A3*N*(N+1)*N*(N+1) H0=H0+AK2*K**2*A2 H0=H0+AK2K*K**4*A2 H0=H0+AK2JK*K**4*A2*N*(N+1) H0=H0+AK2J*K**2*A2*N*(N+1) H0=H0+AK2JJ*K**2*A2*N*(N+1)*N*(N+1) C PRINT*,' H0= ',H0 H0=H0+AK3*K*A4 H0=H0+AK3B*K*A4B H0=H0+AK3BB*K*A4BB H0=H0+AK3BJ*K*A4B*N*(N+1) H0=H0+AK3J*K*A4*N*(N+1) H0=H0+AK3JJ*K*A4*N*(N+1)*N*(N+1) H0=H0+AK3K*K*A4*(K**2) H0=H0+AK3KK*K*A4*(K**4) H0=H0+AK3KJ*K*A4*(K**2)*N*(N+1) H0=H0+AK4*A5 H0=H0+AK4B*A5B H0=H0+AK4BB*A5BB H0=H0+AK4BJ*A5B*N*(N+1) H0=H0+AMV*A5*N*(N+1) H0=H0+AMVJ*A5*N*(N+1)*N*(N+1) H0=H0+BK1*A5*K**2 H0=H0+BK1K*A5*K**4 H0=H0+BK1J*A5*K**2*N*(N+1) H0=H0+AK5*K**2*(A6-0.5*A7-0.5*A1) H0=H0+AK5K*K**4*(A6-0.5*A7-0.5*A1) H0=H0+AK5J*K**2*(A6-0.5*A7-0.5*A1)*N*(N+1) H0=H0+AK6*K*(A3-0.5*A8-0.5*A9+A3-0.5*A14-0.5*A15) H0=H0+AK66*K*(A3-0.5*A8B-0.5*A9B+A3-0.5*A14B-0.5*A15B) H0=H0+AK6K*K*(A3-0.5*A8-0.5*A9+A3-0.5*A14-0.5*A15)*K**2 H0=H0+AK6J*K*(A3-0.5*A8-0.5*A9+A3-0.5*A14-0.5*A15)*N*(N+1) H0=H0+AK7*(A2-0.5*A10-0.5*A11+A2-0.5*A16-0.5*A17) H0=H0+AK9*(A2-0.5*A10B-0.5*A11B+A2-0.5*A16B-0.5*A17B) H0=H0+AK7K*(A2-0.5*A10-0.5*A11+A2-0.5*A16-0.5*A17)*K**2 H0=H0+AK7J*(A2-0.5*A10-0.5*A11+A2-0.5*A16-0.5*A17)*N*(N+1) C PRINT*,' AK7=',AK7,' A1=',A1,'RHOKSIG=',RHOKSIG,'RHOKK=', C 1RHOKK,'N=',N,'K=',K,'IVTOR=',IVTOR H0=H0-DK*K**4*A6 H0=H0+HJ*(N*(N+1))**3*A6 H0=H0+HJK*(N*(N+1)*K)**2*A6 H0=H0+HKJ*N*(N+1)*K**4*A6 H0=H0+HLKJ*K**6*N*(N+1)*A6 H0=H0+HK*K**6*A6 C PRINT*,' CONV = ',CONV C H0 = H0/CONV C PRINT*,' H0= ',H0 IF(IVTOR.EQ.IVTORP) H0=H0+ETOR(IK,IVTOR) C PRINT*,ETOR(IK,IVTOR) RETURN END C********************************************************** FUNCTION H1(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C********************************************************* C SETS UP THE ELEMENT K/K+OR-1,VT/VT' OF THE HAMILTONIAN C (ODAB) CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK DIMENSION CST(2,117),ETOR(KDMX,9),A(9,NTORMX,KDMX) RHO=CST(LO,102) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) ODAB=CST(LO,28) ODABJ=CST(LO,98) ODABJJ=CST(LO,116) ODABJK=CST(LO,117) ODABK=CST(LO,108) DAC=CST(LO,44) DACJ=CST(LO,96) DACK=CST(LO,109) DAC12=CST(LO,94) ODAB6=CST(LO,58) IK=K+N +1 ISIGN=1 IF (KP.EQ.(K-1)) ISIGN=-1 A1=0.0 A2=0.0 A3=0.0 A4=0.0 A5=0.0 A6=0.0 A3B=0.0 A4B=0.0 DO 10 II=1,NDIMTO COEFA=A(IVTOR,II,IK) COEFA1=A(IVTORP,II,IK+ISIGN) IF (II.EQ.NDIMTO ) THEN COEFAA=0.0 COFAA1 =0.0 ELSE COEFAA=A(IVTOR,II+1,IK) COFAA1 =A(IVTORP,II+1,IK+ISIGN) ENDIF IF(II.GT.(NDIMTO-2)) THEN COEFA3=0.0 COEFA4=0.0 ELSE COEFA3=A(IVTOR,II+2,IK) COEFA4=A(IVTORP,II+2,IK+ISIGN) ENDIF IF(II.EQ.1) THEN COEFA2 =0.0 ELSE COEFA2 =A(IVTORP,II-1,IK+ISIGN) ENDIF C FOR SIN 12 alpha :ndimto-4 C FOR SIN 9 alpha :ndimto-3 C IF(II.GT.(NDIMTO-4)) THEN IF(II.GT.(NDIMTO-3)) THEN COFAA12=0.0 ELSE C COFAA12=A(IVTORP,II+4,IK+ISIGN) COFAA12=A(IVTORP,II+3,IK+ISIGN) ENDIF C IF(II.LE.4) THEN IF(II.LE.3) THEN COEFA12=0.0 ELSE C COEFA12=A(IVTORP,II-4,IK+ISIGN) COEFA12=A(IVTORP,II-3,IK+ISIGN) ENDIF A1=A1+COEFA*COEFA1 A2=A2+COEFAA*COEFA1 A3=A3+COEFA*COFAA1 A3B=A3B+COEFA*COFAA12 A4=A4+COEFA*COEFA2 A4B=A4B+COEFA*COEFA12 A5=A5+COEFA1*COEFA3 A6=A6+COEFA*COEFA4 10 CONTINUE C PRINT*,' A1= ',A1,' A2= ',A2,' A3= ',A3 NN=N*(N+1)-K*(K+ISIGN) ANN=FLOAT(NN) H1=ODAB*SQRT(ANN)*(K+0.5*ISIGN) H1=H1+ODABJ*SQRT(ANN)*(K+0.5*ISIGN)*(N*(N+1)) H1=H1+ODABJJ*SQRT(ANN)*(K+0.5*ISIGN)*(N*(N+1)*N*(N+1)) H1=H1+ODABK*SQRT(ANN)*0.5*(K**3+(K+1.0*ISIGN)**3) H1=H1+ODABJK*SQRT(ANN)*N*(N+1)*0.5*(K**3+(K+1.0*ISIGN)**3) H1=H1*(A1-0.5*(A2+A3)) H1=H1+ODAB6*SQRT(ANN)*(K+0.5*ISIGN)*(A1-0.5*(A5+A6)) H1=H1+(-ISIGN)*(DAC/2.0)*SQRT(ANN)*(K+0.5*ISIGN)*(A4-A3) H1=H1+(-ISIGN)*(DACJ/2.0)*SQRT(ANN)*(K+0.5*ISIGN)*(A4-A3) 1*(N*(N+1)) H1=H1+(-ISIGN)*(DACK/2.0)*SQRT(ANN) &*0.5*(K**3+(K+1.0*ISIGN)**3) &*(A4-A3) H1=H1+(-ISIGN)*(DAC12/2.0)*SQRT(ANN)*(K+0.5*ISIGN)* 1(A4B-A3B) C H1 = H1/CONV RETURN END C*********************************************************** FUNCTION H2(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C C SETS UP THE K/K+OR-1 ELEMENTS(DAB) CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK DIMENSION CST(2,117),ETOR(KDMX,9),A(9,NTORMX,KDMX) RHO=CST(LO,102) DAB=CST(LO,27) DABJ=CST(LO,29) DABK=CST(LO,30) ODABK=CST(LO,108) DACK=CST(LO,109) ODELTA=CST(LO,47) DELTA=CST(LO,48) DELTAB=CST(LO,65) ODELTB=CST(LO,66) ESPOIR=CST(LO,90) DABJJ=CST(LO,114) DABJK=CST(LO,115) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) IK=K+N +1 ISIGN=1 IF (KP.EQ.(K-1)) ISIGN=-1 A1=0.0 A2=0.0 A3=0.0 A2B=0.0 A3B=0.0 A4=0.0 DO 11 II=1,NDIMTO KL=II-KTRON1 COEFA=A(IVTOR,II,IK) COEFA1=A(IVTORP,II,IK+ISIGN) C RHOEFF=RHO+RHOJ*N*(N+1)+RHOK*K**2 RHOEFF=RHO RHOKSG =(3.*KL+ISIG+RHOEFF*(K+ISIGN))*((K+ISIGN)**2) RHOKSH =(3.*KL+ISIG+RHOEFF*(K+ISIGN))**3*((K+ISIGN)**2) RHOKS2 =(3.*KL+ISIG+RHOEFF*K)*(K**2) RHOKSI =(3.*KL+ISIG+RHOEFF*K)**3*(K**2) RHOKS3 =3.*KL+ISIG+RHOEFF*K RHOKS4 =3.*KL+ISIG+RHOEFF*(K+ISIGN) A1=A1+COEFA*COEFA1 A2=A2+COEFA*COEFA1*(RHOKSG +RHOKS2 ) A2B=A2B+COEFA*COEFA1*(RHOKSH +RHOKSI) A3=A3+COEFA*COEFA1*(RHOKS3 **2*K+RHOKS4 **2*(K+ISIGN)) A3B=A3B+COEFA*COEFA1*(RHOKS3 **4*K+RHOKS4 **4*(K+ISIGN)) A4=A4+COEFA*COEFA1*(RHOKS3)**3 11 CONTINUE C PRINT*,' A1= ',A1 C PRINT*,' A1= ',A1 NN=N*(N+1)-K*(K+ISIGN) ANN=FLOAT(NN) C H2=DAB*SQRT(ANN)*(K+0.5*ISIGN)*A1+DABJ*SQRT(ANN)*(K+0.5*ISIGN) C * *A1*N*(N+1)+DABK*SQRT(ANN)*(K+0.5*ISIGN)*A1*K**2 C******DABK TERM MADE HERMITIAN BY JTH 18 APR 92. C H2=DAB*SQRT(ANN)*(K+0.5*ISIGN)*A1+DABJ*SQRT(ANN)*(K+0.5*ISIGN) C * *A1*N*(N+1)+DABK*SQRT(ANN)*(K+0.5*ISIGN)*A1*K**2 H2=DAB*SQRT(ANN)*(K+0.5*ISIGN)*A1 * +DABJ*SQRT(ANN)*(K+0.5*ISIGN)*A1*N*(N+1) * +DABK*SQRT(ANN)*0.5*(K**3+(K+1.0*ISIGN)**3)*A1 H2=H2+DABJJ*SQRT(ANN)*(K+0.5*ISIGN)*A1*N*(N+1)*N*(N+1) H2=H2+DABJK*SQRT(ANN)*0.5*(K**3+(K+1.0*ISIGN)**3)*A1*N*(N+1) H2=H2+ODELTA*A2*SQRT(ANN)*0.5 H2=H2+ODELTB*A2B*SQRT(ANN)*0.5 H2=H2+DELTA*A3*SQRT(ANN)*0.5 H2=H2+DELTAB*A3B*SQRT(ANN)*0.5 H2=H2+ESPOIR*A4*SQRT(ANN)*0.5 C H2 = H2/CONV C PRINT*,' H2= ',H2 RETURN END C************************************************************ FUNCTION H3(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C C SETS UP THE K/K+-2 ELEMENTS CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK DIMENSION CST(2,117),A(9,NTORMX,KDMX),ETOR(KDMX,9) OA=CST(LO,1) B=CST(LO,2) C=CST(LO,3) DJ=CST(LO,4) DJK=CST(LO,5) DK=CST(LO,6) ODELN=CST(LO,7) ODELK=CST(LO,8) HJ=CST(LO,9) HJK=CST(LO,10) HKJ=CST(LO,11) HK=CST(LO,12) OHJ=CST(LO,13) OHJK=CST(LO,14) OHK=CST(LO,15) FV=CST(LO,16) OFV=CST(LO,17) GV=CST(LO,18) ALV=CST(LO,19) AK1=CST(LO,20) AK2=CST(LO,21) AK3=CST(LO,22) AK4=CST(LO,23) AK5=CST(LO,24) AK6=CST(LO,25) AK7=CST(LO,26) DAB=CST(LO,27) ODAB=CST(LO,28) OLV=CST(LO,31) C1=CST(LO,32) C2=CST(LO,33) C3=CST(LO,34) C11=CST(LO,57) C11J=CST(LO,112) C11K=CST(LO,113) HLKJ=CST(LO,35) RHO=CST(LO,102) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) C2J=CST(LO,37) GVJ=CST(LO,38) AK2J=CST(LO,39) C1J=CST(LO,40) DBC=CST(LO,45) DBCJ=CST(LO,97) DBC12=CST(LO,95) C4=CST(LO,46) C4J=CST(LO,67) C4K=CST(LO,68) C1K=CST(LO,69) C2K=CST(LO,70) C12=CST(LO,71) C3J=CST(LO,80) C12J=CST(LO,81) IK=K+N +1 ISIGN=2 C PRINT*,' C1= ',C1,' C2= ',C2,' C3= ',C3,' OHJK= ', C 1 OHJK,' ODELN= ',ODELN,' ODELK= ',ODELK,' OHJ= ',OHJ,' OHK= ' C 2 ,OHK C PRINT*,'C4J=',C4J IF (KP.EQ.(K-2)) ISIGN=-2 A1=0.0 A2=0.0 A3=0.0 A4=0.0 A5=0.0 A6=0.0 A7=0.0 A8=0.0 A9=0.0 A10=0.0 A11=0.0 A4B=0.0 A6B=0.0 DO 12 II=1,NDIMTO KL=II-KTRON1 COEFA=A(IVTOR,II,IK) COEFA2=A(IVTORP,II,IK+ISIGN) IF (II.EQ.NDIMTO ) THEN COEFA1=0.0 COFAA2 =0.0 ELSE COEFA1=A(IVTOR,II+1,IK) COFAA2 =A(IVTORP,II+1,IK+ISIGN) ENDIF IF(II.EQ.1) THEN COEFA3 =0.0 ELSE COEFA3 =A(IVTORP,II-1,IK+ISIGN) ENDIF CC IF(II.LE.4) THEN IF(II.LE.3) THEN COEFA12=0.0 ELSE C COEFA12=A(IVTORP,II-4,IK+ISIGN) COEFA12=A(IVTORP,II-3,IK+ISIGN) ENDIF IF(II.GT.(NDIMTO-2)) THEN COEFA4=0.0 COEFAA4=0.0 ELSE COEFA4=A(IVTOR,II+2,IK) COEFAA4=A(IVTORP,II+2,IK+ISIGN) ENDIF C IF(II.GT.(NDIMTO-4)) THEN IF(II.GT.(NDIMTO-3)) THEN COFAA12=0.0 ELSE C COFAA12=A(IVTORP,II+4,IK+ISIGN) COFAA12=A(IVTORP,II+3,IK+ISIGN) ENDIF C RHOEFF=RHO+RHOJ*N*(N+1)+RHOK*K**2 RHOEFF=RHO RHOKSG =3.*KL+ISIG+RHOEFF*K RHOKS2 =(3.*KL+ISIG+RHOEFF*(K+ISIGN))*(K+ISIGN) RHOKS3 =(3.*KL+ISIG+RHOEFF*K)*(K) RHOK2=3.*KL+ISIG+RHOEFF*(K+ISIGN) RHOKS4=(3.*KL+ISIG+RHOEFF*K)**3 RHOKS4=RHOKS4*K RHOKS5=(3.*KL+ISIG+RHOEFF*(K+ISIGN))**3 RHOKS5=RHOKS5*(K+ISIGN) RHOKS6=(3.*KL+ISIG+RHOEFF*K)*(K) RHOKS6=RHOKS6*(K+ISIGN)**2 RHOKS7=(3.*KL+ISIG+RHOEFF*(K+ISIGN))*(K+ISIGN) RHOKS7=RHOKS7*(K)**2 A1=A1+COEFA*COEFA2 C PRINT*,' COEFA= ',COEFA,' COEFA1= ',COEFA1,' COEFA2= ', C 1COEFA2,' COEFAA2= ',COEFAA2 A2=A2+COEFA*COEFA2*(RHOKSG **2+RHOK2**2) A3=A3+COEFA1*COEFA2 A4=A4+COEFA*COFAA2 A4B=A4B+COEFA*COFAA12 A5=A5+COEFA*COEFA2*(RHOKSG **4+RHOK2**4) A6=A6+COEFA*COEFA3 A6B=A6B+COEFA*COEFA12 A7=A7+COEFA*COEFA2*(RHOKS2 +RHOKS3 ) A8=A8+COEFA2*COEFA4 A9=A9+COEFA*COEFAA4 A10=A10+COEFA*COEFA2*(RHOKS4+RHOKS5) A11=A11+COEFA*COEFA2*(RHOKS6+RHOKS7) 12 CONTINUE C PRINT*,' A1= ',A1,' A2= ',A2,' A3= ',A3,' A4= ',A4 C PRINT*,' K= ',K,' KP= ',KP,' IVTOR= ',IVTOR,' IVTORP= ', C 1IVTORP,' ISIGN= ',ISIGN H3=(B-C)/4.-ODELN*N*(N+1)-0.5*ODELK*(K**2+(K+ISIGN)**2)+ * OHJ*(N*(N+1))**2+0.5*OHJK*N*(N+1)*(K**2+(K+ISIGN)**2) * +0.5*OHK*(K**4+(K+ISIGN)**4) C PRINT*,' H3= ',H3 H3=A1*H3 H3=H3+0.5*C1*A2 H3=H3+0.5*C1J*A2*N*(N+1) C PRINT*,' H3= ',H3 H3=H3+0.5*C2*(A1-0.5*(A3+A4)) H3=H3+0.5*C2J*(A1-0.5*(A3+A4))*(N*(N+1)) H3=H3+0.5*C3*A5 H3=H3+0.5*C3J*A5*N*(N+1) H3=H3+0.5*C4*A7 H3=H3+0.5*C12*A10 H3=H3+0.5*C12J*A10*N*(N+1) H3=H3+0.5*C4J*A7*N*(N+1) C IF(N.EQ.1) PRINT*,'H3=',H3 H3=H3+0.5*C1K*(K**2+(K+ISIGN)**2)*A2 H3=H3+0.5*C2K*(A1-0.5*(A3+A4))*(K**2+(K+ISIGN)**2) H3=H3+0.5*C4K*A11 H3=H3+0.5*C11*(A1-0.5*(A8+A9)) H3=H3+0.5*C11J*(A1-0.5*(A8+A9))*N*(N+1) H3=H3+0.5*C11K*(A1-0.5*(A8+A9))*(K**2+(K+ISIGN)**2) ANN=N*(N+1)-K*(K+0.5*ISIGN) ALL=N*(N+1)-(K+0.5*ISIGN)*(K+ISIGN) H3=H3*SQRT(ANN*ALL) C PRINT*,' H3= ',H3 H3=H3+((A6-A4)*(-ISIGN)*0.5*DBC*(1./4.))*(SQRT(ANN*ALL)) H3=H3+((A6-A4)*(-ISIGN)*0.5*DBCJ*(1./4.))*(SQRT(ANN*ALL)) 1*(N*(N+1)) H3=H3+((A6B-A4B)*(-ISIGN)*0.5*DBC12*(1./4.))*(SQRT(ANN*ALL)) RETURN END C************************************************************** FUNCTION H4(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************** C SETS UP THE ELEMENT K/K,VT,VT'(OLV) CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK DIMENSION CST(2,117),ETOR(KDMX,9),A(9,NTORMX,KDMX) RHO=CST(LO,102) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) OLV=CST(LO,31) IK=K+N +1 A1=0.0 DO 25 II=1,NDIMTO KL=II-KTRON1 COEFA=A(IVTOR,II,IK) COEFAA=A(IVTORP,II,IK) C RHOEFF=RHO+RHOJ*N*(N+1)+RHOK*K**2 RHOEFF=RHO RHOKSG =3.*KL+ISIG+RHOEFF*K A1=A1+COEFA*COEFAA*RHOKSG 25 CONTINUE C H4=OLV*N*(N+1)*K**3*A1 C OLV SHOULD BE A J(J+1) DEPENDANCE OF ALV H4=OLV*N*(N+1)*K*A1*N*(N+1) C H4=H4/CONV RETURN END C************************************************************ SUBROUTINE HRTSET(CST,LO,ETOR,A) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************ C SETS UP THE ROTATION-TORSION HAMILTONIAN CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117),A(9,NTORMX,KDMX),ETOR(KDMX,9) DO 3 I=1,NROTOR DO 4 J=1,I C INITIALISATION H(I,J)=0.0 4 CONTINUE 3 CONTINUE DO 7 I=1,NROTOR DO 8 J=1,I C LEFT-SIDE MATRIX ELEMENT IVTOR=(I-1)/(2*N+1)+1 K=I+N*(1 -2*IVTOR)-IVTOR C RIGHT-SIDE MATRIX ELEMENT IVTORP=(J-1)/(2*N+1)+1 KP=J+N*(1 -2*IVTORP)-IVTORP C C K=KP C PRINT*,' K= ',K,' KP= ',KP,' IVTOR= ',IVTOR,' IVTORP= ', C 1IVTORP IF (K.EQ.KP) THEN C IF(IVTOR.EQ.IVTORP) THEN H(I,J)= H0(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) 1 + H4(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) C PRINT*,' (I,J)= ',H(I,J),' I= ',I,' J= ',J C PRINT*,'H0=',H0(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A), C 1'K=',K,'IVTOR=',IVTOR,'KP=',KP,'IVTORP=',IVTORP, C 2'LO=',LO C PRINT*,'H4=',H4(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) C ENDIF ELSE IF (IABS(K-KP).EQ.1) THEN IF (IABS(IVTOR-IVTORP).LT.9) THEN H(I,J)= H2(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) 1 + H1(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) C PRINT*,' H(I,J)= ',H(I,J),' I= ',I,' J= ',J C PRINT*,'H2=',H2(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) C PRINT*,'H1=',H1(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) ENDIF ELSE IF(IABS(K-KP).EQ.2) THEN IF (IABS(IVTOR-IVTORP).LT.9) THEN H(I,J)= H3(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) C PRINT*,' H(I,J)= ',H(I,J),' I= ',I,' J= ',J C PRINT*,'H3=',H3(K,IVTOR,KP,IVTORP,CST,LO,ETOR,A) ENDIF ENDIF 8 CONTINUE 7 CONTINUE DO 20 I=2,NROTOR DO 21 J=1,I-1 H(J,I)=H(I,J) 21 CONTINUE 20 CONTINUE RETURN END C******************************************************** SUBROUTINE HTORSE (K,CST,LO) IMPLICIT REAL*8 ( A-H,O-Z ) C C SETS UP THE TORSION HAMILTANIAN C C THIS SUBROUTINE SETS UP THE TORSION HAMILTONIAN C CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK COMMON/ROTOR/NROTOR,NDIMTO COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117) C C LOOP ON THE BASIS FUNCTIONS F=CST(LO,101) RHO=CST(LO,102) C RHOJ=CST(LO,60) C RHOK=CST(LO,61) RHORHO=CST(LO,107) V3=CST(LO,103) V6=CST(LO,104) V9=CST(LO,105) V12=CST(LO,106) C INITIALISATION DO 3 I=1,NDIMTO DO 4 J=1,I H(I,J)=0.0 4 CONTINUE 3 CONTINUE C C KTOR:DLEFT SIDE OF THE MATRIX ELEMENT C KTORP:DRIGHT SIDE OF THE MATRIX ELEMENT C DO 10 I=1,NDIMTO DO 11 J=1,I KTOR=I - KTRON1 KTORP=J - KTRON1 C RHOEFF=RHORHO+RHOJ*N*(N+1)+RHOK*K**2 RHOEFF=RHORHO C C KTOR=KTORP C C PRINT*,' KTOR= ',KTOR,' KTORP= ',KTORP IF (KTOR.EQ.KTORP) THEN c H(I,J) = F*(3.*KTOR+ISIG+RHOEFF*K)**2+V3/2. c & +V6/2.+V9/2.+V12/2. cdebugik, 2Frho.pa.pg sign changed! H(I,J) = F*(3.*KTOR+ISIG-RHOEFF*K)**2+V3/2. & +V6/2.+V9/2.+V12/2. C PRINT*,' H= ',H(I,J) ELSE IF(IABS(KTOR-KTORP).EQ.1) THEN C C KTOR=KTORP+1 OR KTOR=KTORP-1 C H(I,J) = -V3/4. ELSE IF (IABS(KTOR-KTORP).EQ.2) THEN C C KTOR=KTORP+2 OR KTOR=KTORP-2 C H(I,J) = -V6/4. ELSE IF(IABS(KTOR-KTORP).EQ.3) THEN H(I,J)=-V9/4. ELSE IF(IABS(KTOR-KTORP).EQ.4) THEN H(I,J)=-V12/4. ENDIF 11 CONTINUE 10 CONTINUE C SYMETRISES THE MATRIX DO 5 I=2,NDIMTO DO 6 J=1,I-1 H(J,I) = H(I,J) 6 CONTINUE 5 CONTINUE RETURN END c******************************************************************** SUBROUTINE PARITY(A,EGVC4,II) c my insertion IMPLICIT REAL*8 ( A-H,O-Z ) C THIS SUBROUTINE CALCULATES THE TRUE PARITY OF THE CURRENT STATE C (+or-).(-1)J+VT = "true" parity C**************************************************************** CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/ROTOR/NROTOR,NDIMTO COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/QUANT2/KKK(NDMX),NVT(NDMX),IPRITY(NDMX) COMMON/C2C/PERC2L,PER23L(NDMX) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION A(9,NTORMX,KDMX) DIMENSION EGVC4(KDMX,NTORMX,NDMX) C NLOW=N NLOW1 = NLOW+1 C P OR PRIME STANDS FOR LEFT OR LOWER KNP2=2*NLOW KLOW=KK IVVLO=IV PERC2L=0.0 PER23L(II)=0. C INITIALIZATION OF THE SUM ON TORSIONAL EIGENFUNCTIONS AN00A=0.0 AN00B=0.0 AN00C=0.0 AN01A=0.0 AN01B=0.0 AN01C=0.0 AN0M1A=0.0 AN0M1B=0.0 AN0M1C=0.0 AN10A=0.0 AN10B=0.0 AN10C=0.0 AN11A=0.0 AN11B=0.0 AN11C=0.0 AN1M1A=0.0 AN1M1B=0.0 AN1M1C=0.0 ANM10A=0.0 ANM10B=0.0 ANM10C=0.0 ANM11A=0.0 ANM11B=0.0 ANM11C=0.0 ANM1MA=0.0 ANM1MB=0.0 ANM1MC=0.0 C C EXPECTATION VALUE OF C2(C) AND OF (23)* FOR CURRENT STATE C DO 304 IVTORP=1,9 DO 302 IVTOR1=1,9 DO 300 ILOW=1,2*NLOW+1 KBASLO=ILOW-NLOW-1 IOPPLO=-KBASLO+NLOW+1 NPKLO=NLOW+KBASLO EXPPLO=1. IF(MOD(NPKLO,2).NE.0) EXPPLO=-1. A2LO=0.0 A3LO=0.0 DO 298 III=1,NDIMTO IIIOPP=2*KTRONC+1-(III-1) COEFAL=A(IVTORP,III,ILOW) COEFAAL=A(IVTOR1,III,IOPPLO) IF(ISIG.EQ.0) THEN COEF4L=A(IVTOR1,IIIOPP,IOPPLO) ELSE COEF4L=0.0 ENDIF A2LO=A2LO+COEFAL*COEFAAL A3LO=A3LO+COEFAL*COEF4L 298 CONTINUE PERC2L=PERC2L+A2LO*EXPPLO* 1 EGVC4(ILOW,IVTORP,II)*EGVC4(IOPPLO,IVTOR1,II) PER23L(II)=PER23L(II)+A3LO*EXPPLO* 1 EGVC4(ILOW,IVTORP,II)*EGVC4(IOPPLO,IVTOR1,II) C IF(N.EQ.1) PRINT*,'EGVC4(K)=',EGVC4(ILOW,IVTORP,II), C 1 'EGVC4(-K)=',EGVC4(IOPPLO,IVTOR1,II),'ILOW=',ILOW, C 2 'IVTORP=',IVTORP,'II=',II C PRINT*,'PER23L=',PER23L,'A3LO=',A3LO,'NLOW=',NLOW,'KK=',KK,'IV=',IV C 1,'EGVC3=',EGVC3(ILOW,IVTORP),'EGVC3(-K)=',EGVC3(IOPPLO,IVTOR1) 300 CONTINUE 302 CONTINUE 304 CONTINUE C PRINT*,'PER23L=',PER23L,'N=',N,'IV=',IV,'KK=',KK RETURN END C*********************************************************** FUNCTION IREP(SYMBOL) IMPLICIT REAL*8 ( A-H,O-Z ) C THIS FUNCTION RECOGNIZE THE SYMBOLIC CONSTANT ENCOUNTERED C ON TAPE 5 AND RETURNS ITS POSITION IN THE ABC STRING CHARACTER SYMBOL*6,ABC*702 COMMON/SYMB/ABC DO 1 K=1,(117) IF(SYMBOL.EQ.ABC(6*(K-1)+1:6*K)) THEN IREP=K RETURN END IF 1 CONTINUE WRITE(6,100) SYMBOL 100 FORMAT(//,' NO MATCHING CONSTANT SYMBOL FOR ',A6,/) STOP END SUBROUTINE LINECA (CST,NPP) IMPLICIT REAL*8 ( A-H,O-Z ) C THIS SUBROUTINE CALCULATES AND WRITES C 1) THE OBS-CALC FOR ALL THE EXPERIMENTAL DATA C 2) THE ENERGY OF THE LINES GIVEN C C********************************************************* CHARACTER ABC*702 CHARACTER*6 PARA,REF CHARACTER*1 AST CHARACTER*1 PR,PROBS,PRSUP,PRINF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/EIGEN1/EGVL(2*NDMX),EGVC(NDMX,2*NDMX) COMMON/TOR/A(9,NTORMX,2*KDMX),ETOR(2*KDMX,9) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/BIG2/DERIV1(10,961,80),ENER(10,961) COMMON/EXPDAT / ETRANS(80000),IVOBS(2,80000), * NOBS(2,80000),KOBS(2,80000),IBLK(2,80000),IW(80000),W(80000), * NDATA,NADAT,EPSI,NOFIT,NFITMW,NFITIR,NFIMW0 ,NFIMW1 ,NITT, * IAST(80000),NFI004 ,NFI080 ,NFI100 ,NFI020 ,NFI200 ,NFI045 , & NGST0 ,NFI010 ,SW070,NFI070 ,NFIMW2,IKAKC,JMAX,NFI1M0, & NFIR10,NFIR21,NFIR32,NFIMW3,NFIMW4,SWA,SWE,NFITMWA,NFITMWE COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117),IBLSTO(2),EE(2),KC(2,80000) S=0. SMW=0. SWA=0.0 SWE=0. SMW0=0.0 SMW1=0.0 SMW2=0.0 SMW3=0.0 SMW4=0.0 SIR=0. SIR10=0. SIR21=0. SIR32=0. SIG004=0. SIG080=0.0 SIG010=0.0 SIG100=0.0 SIG020=0.0 SIG200=0.0 SIG045=0.0 SIG1M0=0.0 SW070=0.0 SIG32=0.0 WMW=0. WMWA=0. WMWE=0. WMW0=0.0 WMW1=0.0 WMW2=0.0 WMW3=0.0 WMW4=0.0 WIR=0. WIR10=0. WIR21=0. WIR32=0. WIG004=0. WIG080=0.0 WIG010=0.0 WIG100=0.0 WIG020=0.0 WIG200=0.0 WIG045=0.0 WIG1M0=0.0 WIG070=0.0 SGST0=0. IBLSTO(1)=0 IBLSTO(2)=0 WRITE(6,602) 602 FORMAT(/,4X,' UPPER STATE',3X,' LOWER STATE',5X, 1 ' OBSERVATION ',8X,' CALCULATION',6X,' TME ',4X,' WEIGHT', 2 4X,' FITTED OR NO' 3 /,2(2X,' VT J KA KC PARITY')) DO 22 K=1,NDATA IF(NOBS(1,K).GT.JMAX.OR.NOBS(2,K).GT.JMAX) GO TO 22 C IF(K.GT.NADAT) THEN C ISIG=1 C ELSE C ISIG=0 C ENDIF C DO 24 II=1,2 C IF(II.EQ.1) THEN C NI=1 C M=1 C ELSE C NI=NDMX+1 C M=KDMX+1 C ENDIF C N=NOBS(II,K) C IV=IVOBS(II,K) C KK=KOBS(II,K) C PR=PROBS(II,K) C LOO=II C III=II C IF(IRMW.EQ.1.AND.II.EQ.2) III=1 C IF(IBLK(LOO,K).NE.IBLSTO(LOO)) THEN C IBLSTO(LOO)=IBLK(LOO,K) C CALL ROTTOR(CST,III,EGVL(NI),EGVC(1,NI),ETOR(M,1) C 1 ,A(1,1,M)) C ENDIF C NROTOR=(2*N+1)*9 C CALL ENCAL(E,LOO,EGVL(NI),EGVC(1,NI)) C EE(LOO)=E C24 CONTINUE DO 31 LO=1,2 N=NOBS(LO,K) KK=KOBS(LO,K) IV=IVOBS(LO,K) PR=PROBS(LO,K) IF(PR.EQ.'+') IPR=1 IF(PR.EQ.'-') IPR=0 IF(K.GT.NADAT) THEN ISIG=1 IVV=IV+6 INK=N**2+N+KK+1 ELSE ISIG=0 IVV=IV+1 INK=N**2+N+KK+1-IPR*N ENDIF EE(LO)=ENER(IVV,INK) 31 CONTINUE ECALC=TERM+EE(2)-EE(1) C IF(ECALC.LT.0.) THEN c NSUP=NOBS(2,K) c NINF=NOBS(1,K) c KSUP=KOBS(2,K) c KINF=KOBS(1,K) c PRSUP=PROBS(2,K) c PRINF=PROBS(1,K) c NOBS(2,K)=NINF c NOBS(1,K)=NSUP c KOBS(2,K)=KINF c KOBS(1,K)=KSUP c PROBS(2,K)=PRINF c PROBS(1,K)=PRSUP c ESUP=EE(2) c EINF=EE(1) c EE(2)=EINF c EE(1)=ESUP c ECALC=TERM+EE(2)-EE(1) c ENDIF C HERE, THE USER HAS THE POSSIBILITY OF HAVING A PRINT OUT C INCLUDING THE JKAKC NOTATION AND THE INTERNAL ROTOR NOTATION C BUT, FOR THE E STATES,SINCE THERE IS NO SYMMETRY RULES ALLOW C THE CONVERSION, IT IS BASED ON ENERGY LEVELS ORDER,M WHICH cjth IS MOLECULE DEPENDANT:DHERE IT IS FOR CH3COOH ONLY! IF(IKAKC.NE.1) GO TO 60 IF(ISIG.EQ.0) THEN IF(MOD(KOBS(1,K),2).EQ.0) THEN IF(PROBS(1,K).EQ.'+') THEN KC(1,K)=NOBS(1,K)-KOBS(1,K) ELSE KC(1,K)=NOBS(1,K)-KOBS(1,K)+1 ENDIF ELSEIF(MOD(KOBS(1,K),2).NE.0) THEN IF(PROBS(1,K).EQ.'+') THEN KC(1,K)=NOBS(1,K)-KOBS(1,K)+1 ELSE KC(1,K)=NOBS(1,K)-KOBS(1,K) ENDIF ENDIF IF(MOD(KOBS(2,K),2).EQ.0) THEN IF(PROBS(2,K).EQ.'+') THEN KC(2,K)=NOBS(2,K)-KOBS(2,K) ELSE KC(2,K)=NOBS(2,K)-KOBS(2,K)+1 ENDIF ELSEIF(MOD(KOBS(2,K),2).NE.0) THEN IF(PROBS(2,K).EQ.'+') THEN KC(2,K)=NOBS(2,K)-KOBS(2,K)+1 ELSE KC(2,K)=NOBS(2,K)-KOBS(2,K) ENDIF ENDIF ELSE ccik debug change of rho sign KT1K=KOBS(1,K)*(-1)**IVOBS(1,K) c IF(KT1K.EQ.0.OR.KT1K.EQ.-1) THEN IF(KT1K.EQ.0.OR.KT1K.EQ.1) THEN KC(1,K)=NOBS(1,K) c ELSEIF(KT1K.EQ.1.OR.KT1K.EQ.-2) THEN ELSEIF(KT1K.EQ.-1.OR.KT1K.EQ.2) THEN KC(1,K)=NOBS(1,K)-1 c ELSEIF(KT1K.EQ.2.OR.KT1K.EQ.-3) THEN ELSEIF(KT1K.EQ.-2.OR.KT1K.EQ.3) THEN KC(1,K)=NOBS(1,K)-2 c ELSEIF(KT1K.EQ.3.OR.KT1K.EQ.-4) THEN ELSEIF(KT1K.EQ.-3.OR.KT1K.EQ.4) THEN KC(1,K)=NOBS(1,K)-3 cjth I changed the Kc labels below to be correct for CH3COOH c ELSEIF(KT1K.EQ.4.OR.KT1K.EQ.-5) THEN ELSEIF(KT1K.EQ.-4.OR.KT1K.EQ.5) THEN KC(1,K)=NOBS(1,K)-4 ELSEIF(KT1K.EQ.-5.OR.KT1K.EQ.6) THEN KC(1,K)=NOBS(1,K)-5 ELSEIF(KT1K.EQ.-6.OR.KT1K.EQ.7) THEN KC(1,K)=NOBS(1,K)-6 ELSEIF(KT1K.EQ.-7.OR.KT1K.EQ.8) THEN KC(1,K)=NOBS(1,K)-7 ELSEIF(KT1K.EQ.-8.OR.KT1K.EQ.9) THEN KC(1,K)=NOBS(1,K)-8 ELSEIF(KT1K.EQ.-9.OR.KT1K.EQ.10) THEN KC(1,K)=NOBS(1,K)-9 ELSEIF(KT1K.EQ.-10.OR.KT1K.EQ.11) THEN KC(1,K)=NOBS(1,K)-10 ELSEIF(KT1K.EQ.-11.OR.KT1K.EQ.12) THEN KC(1,K)=NOBS(1,K)-11 ELSEIF(KT1K.EQ.-12.OR.KT1K.EQ.13) THEN KC(1,K)=NOBS(1,K)-12 ELSEIF(KT1K.EQ.-13.OR.KT1K.EQ.14) THEN KC(1,K)=NOBS(1,K)-13 ELSEIF(KT1K.EQ.-14.OR.KT1K.EQ.15) THEN KC(1,K)=NOBS(1,K)-14 ELSEIF(KT1K.EQ.-15.OR.KT1K.EQ.16) THEN KC(1,K)=NOBS(1,K)-15 ENDIF KT2K=KOBS(2,K)*(-1)**IVOBS(2,K) c IF(KT2K.EQ.0.OR.KT2K.EQ.-1) THEN IF(KT2K.EQ.0.OR.KT2K.EQ.1) THEN KC(2,K)=NOBS(2,K) c ELSEIF(KT2K.EQ.1.OR.KT2K.EQ.-2) THEN ELSEIF(KT2K.EQ.-1.OR.KT2K.EQ.2) THEN KC(2,K)=NOBS(2,K)-1 ELSEIF(KT2K.EQ.-2.OR.KT2K.EQ.3) THEN KC(2,K)=NOBS(2,K)-2 ELSEIF(KT2K.EQ.-3.OR.KT2K.EQ.4) THEN KC(2,K)=NOBS(2,K)-3 cjth I changed the Kc labeling below to be correct for CH3COOH. ELSEIF(KT2K.EQ.-4.OR.KT2K.EQ.5) THEN KC(2,K)=NOBS(2,K)-4 ELSEIF(KT2K.EQ.-5.OR.KT2K.EQ.6) THEN KC(2,K)=NOBS(2,K)-5 ELSEIF(KT2K.EQ.-6.OR.KT2K.EQ.7) THEN KC(2,K)=NOBS(2,K)-6 ELSEIF(KT2K.EQ.-7.OR.KT2K.EQ.8) THEN KC(2,K)=NOBS(2,K)-7 ELSEIF(KT2K.EQ.-8.OR.KT2K.EQ.9) THEN KC(2,K)=NOBS(2,K)-8 ELSEIF(KT2K.EQ.-9.OR.KT2K.EQ.10) THEN KC(2,K)=NOBS(2,K)-9 ELSEIF(KT2K.EQ.-10.OR.KT2K.EQ.11) THEN KC(2,K)=NOBS(2,K)-10 ELSEIF(KT2K.EQ.-11.OR.KT2K.EQ.12) THEN KC(2,K)=NOBS(2,K)-11 ELSEIF(KT2K.EQ.-12.OR.KT2K.EQ.13) THEN KC(2,K)=NOBS(2,K)-12 ELSEIF(KT2K.EQ.-13.OR.KT2K.EQ.14) THEN KC(2,K)=NOBS(2,K)-13 ELSEIF(KT2K.EQ.-14.OR.KT2K.EQ.15) THEN KC(2,K)=NOBS(2,K)-14 ELSEIF(KT2K.EQ.-15.OR.KT2K.EQ.16) THEN KC(2,K)=NOBS(2,K)-15 ENDIF ENDIF 60 IF(IAST(K).EQ.1) THEN EOBS=ETRANS(K) TME2=EOBS-ECALC S=S+TME2**2*W(K) IF(IW(K).GE.1000) THEN WMW=WMW+TME2**2*W(K) SMW=SMW+TME2**2 cjth added some if's here to separate A and E stuff. if(k.le.nadat) then SWA=SWA+TME2**2 WMWA=WMWA+TME2**2*W(K) elseif(k.gt.nadat) then SWE=SWE+TME2**2 WMWE=WMWE+TME2**2*W(K) endif IF(IVOBS(1,K).EQ.0.AND.IW(K).GE.1000) THEN WMW0=WMW0+TME2**2*W(K) SMW0=SMW0+TME2**2 ELSEIF(IVOBS(1,K).EQ.1.AND.IW(K).GE.1000) THEN WMW1=WMW1+TME2**2*W(K) SMW1=SMW1+TME2**2 ELSEIF(IVOBS(1,K).EQ.2.AND.IW(K).GE.1000) THEN WMW2=WMW2+TME2**2*W(K) SMW2=SMW2+TME2**2 ELSEIF(IVOBS(1,K).EQ.3.AND.IW(K).GE.1000) THEN WMW3=WMW3+TME2**2*W(K) SMW3=SMW3+TME2**2 ELSEIF(IVOBS(1,K).EQ.4.AND.IW(K).GE.1000) THEN WMW4=WMW4+TME2**2*W(K) SMW4=SMW4+TME2**2 C ELSEIF(IW(K).EQ.2000) THEN C SWGOOD=SWGOOD+TME2**2*W(K) ENDIF IF(IW(K).EQ.1004.OR.IW(K).EQ.1005.OR.IW(K).EQ.1003) THEN WIG004=WIG004+TME2**2*W(K) SIG004=SIG004+TME2**2 ELSEIF(IW(K).EQ.1080) THEN WIG080=WIG080+TME2**2*W(K) SIG080=SIG080+TME2**2 ELSEIF(IW(K).EQ.1070) THEN WIG070=WIG070+TME2**2*W(K) SW070=SW070+TME2**2 ELSEIF(IW(K).EQ.1010.OR.IW(K).EQ.1008) THEN WIG010=WIG010+TME2**2*W(K) SIG010=SIG010+TME2**2 ELSEIF(IW(K).EQ.1100) THEN WIG100=WIG100+TME2**2*W(K) SIG100=SIG100+TME2**2 ELSEIF(IW(K).EQ.1020) THEN WIG020=WIG020+TME2**2*W(K) SIG020=SIG020+TME2**2 ELSEIF(IW(K).EQ.1200.OR.IW(K).EQ.1180.OR.IW(K).EQ.1160 & .OR.IW(K).EQ.1150.OR.IW(K).EQ.1130) THEN WIG200=WIG200+TME2**2*W(K) SIG200=SIG200+TME2**2 ELSEIF(IW(K).EQ.1030.OR.IW(K).EQ.1040.OR.IW(K).EQ.1050) & THEN WIG045=WIG045+TME2**2*W(K) SIG045=SIG045+TME2**2 ELSEIF(IW(K).EQ.1990) THEN WIG1M0=WIG1M0+TME2**2*W(K) SIG1M0=SIG1M0+TME2**2 ENDIF EOBS=EOBS*29979.2458 ECALC=ECALC*29979.2458 TME2=TME2*29979.2458 ELSEIF(IW(K).EQ.1) THEN WIR10=WIR10+TME2**2*W(K) SIR10=SIR10+TME2**2 WIR=WIR+TME2**2*W(K) SIR=SIR+TME2**2 ELSEIF(IW(K).EQ.2) THEN WIR21=WIR21+TME2**2*W(K) SIR21=SIR21+TME2**2 WIR=WIR+TME2**2*W(K) SIR=SIR+TME2**2 ELSEIF(IW(K).EQ.3) THEN WIR32=WIR32+TME2**2*W(K) SIR32=SIR32+TME2**2 WIR=WIR+TME2**2*W(K) SIR=SIR+TME2**2 ELSEIF(IW(K).EQ.10) THEN WGST0=WGST0+TME2**2*W(K) SGST0=SGST0+TME2**2 ENDIF C IF(IW(K).GE.1000) IW(K)=IW(K)-1000 C IF(IW(K).EQ.10) IW(K)=10000 C WRITE(6,604) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) C 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), C 2 PROBS(1,K),EOBS,ECALC,TME2,IW(K),IAST(K),REF(K) C604 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, C 1 1X,I3,1X,A1,11X,F12.4,1X,F16.4,1X,F16.4,1X,I5,1X,I3,1X,A1) IF(IW(K).GE.1000) THEN IF(IW(K).EQ.2000) IW(K)=2 IF(IW(K).LT.2000) IW(K)=IW(K)-1000 C C HERE THE USER HAS THE CHOICE BETWEEN A COMPACT FORMAT (THE ONE C WHICH IS NOW PRECEDED BY "C" OR A LESS DENSE FORMAT C c IF(REF(K).NE.'NIST '.OR.REF(K).NE.'KHARK ') then ccjth I modified the write statement to print out E' and E". WRITE(6,604) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), cc 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,REF(K) 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,REF(K),EE(2),EE(1) if(tme2.ge.4.) print*,'attention!' c endif 604 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, 1 1X,I3,1X,A1,1X,F12.3,'(',I3,')',1X,F12.3,2X,F7.3, cc 2 1X,A6) 2 1X,A6,3X,2f13.8) C WRITE(6,604) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) C 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), C 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,IAST(K),REF(K) C604 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, C 1 1X,I3,1X,A1,11X,F12.3,'(',I3,')',1X,F16.3,1X,F16.3,1X,I5, C 21X,A1) ELSEIF(IW(K).EQ.10) THEN IW(K)=5 c WRITE(6,634) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) c 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), c 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,IAST(K),REF(K) c 3 ,EE(2),EE(1) 634 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, 1 1X,I3,1X,A1,11X,F12.4,'(',I1,')',1X,F16.4,1X,F16.4,1X,I5, 2 1X,A6,1X,F9.4,1X,F9.4) ELSEIF(IW(K).EQ.1) THEN IW(K)=4 C C HERE AGAIN THE USER CAN CHOICE A FORMAT C C WRITE(6,635) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) C 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), C 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,IAST(K),REF(K) C635 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, C 1 1X,I3,1X,A1,11X,F12.4,'(',I1,')',1X,F16.4,1X,F16.4,1X,I5, C 21X,A1) cc WRITE(6,635) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) c 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), c 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,REF(K),EE(2),EE(1) 635 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3, 1 1X,I3,1X,A1,1X,F12.4,'(',I1,')',1X,F12.4,1X,F7.4, 2 1X,A6,1X,F9.4,1X,F9.4) ELSEIF(IW(K).EQ.2.OR.IW(K).EQ.3)THEN IW(K)=5 C C HERE AGAIN THE USER CAN CHOICE A FORMAT C C WRITE(6,635) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) C 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), C 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,IAST(K),REF(K) cc WRITE(6,635) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K) cc 1 ,PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), cc 2 PROBS(1,K),EOBS,IW(K),ECALC,TME2,REF(K),EE(2),EE(1) ENDIF ELSEIF(IAST(K).NE.1) THEN EOBS=ETRANS(K) TME3=EOBS-ECALC IF(IW(K).GE.1000) THEN IF(IW(K).LT.2000) IW(K)=IW(K)-1000 EOBS=EOBS*29979.2458 ECALC=ECALC*29979.2458 TME3=TME3*29979.2458 AST='*' ccjth I modified the write statement to get E' and E" printed out. WRITE(6,605) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K), 1 PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), cc 2 PROBS(1,K),EOBS,IW(K),ECALC,TME3,AST,REF(K) 2 PROBS(1,K),EOBS,IW(K),ECALC,TME3,AST,REF(K),EE(2),EE(1) ELSE TME2=EOBS-ECALC AST='*' c WRITE(6,608) IVOBS(2,K),NOBS(2,K),KOBS(2,K),KC(2,K), c 1 PROBS(2,K),IVOBS(1,K),NOBS(1,K),KOBS(1,K),KC(1,K), c 2 PROBS(1,K),EOBS,ECALC,TME2,AST,REF(K),EE(2),EE(1) ENDIF C605 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, C 1 I3,1X,A1,11X,F12.4,1X,F16.4,1X,F16.4,1X,A1,1X,A1) 605 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, C*******FORMAT CHANGE JTH C 1 I3,1X,A1,1X,F12.3,1X,F12.3,1X,F7.3,1X,A1,1X,A1) 1 I3,1X,A1,1X,F12.3,'(',I3,')',1X,F12.3,2X,F7.3 cc 2 ,1X,A1,1X,A6) 2 ,1X,A1,1X,A6,1x,2f13.8) C605 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, C 1 I3,1X,A1,11X,F12.3,1X,F16.3,1X,F16.3,1X,A1,1X,A1) 608 FORMAT(1X,I3,1X,I3,1X,I3,1X,I3,1X,A1,1X,I3,1X,I3,1X,I3,1X, 1 I3,1X,A1,11X,F12.4,1X,F16.4,1X,F16.4,1X,A1,1X,A6, 2 1X,F9.4,1X,F9.4) ENDIF 22 CONTINUE IF((NOFIT-NPP).LE.0) GO TO 23 S=SQRT(S/(NOFIT-NPP)) 23 IF(NFITMW.EQ.0) THEN WMW=0.0 SMW=0.0 ELSE WMW=SQRT(WMW/NFITMW) SMW=SQRT(SMW/NFITMW) ENDIF IF(NFITMWA.EQ.0) THEN WMWA=0.0 SWA=0.0 ELSE WMWA=SQRT(WMWA/NFITMWA) SWA=SQRT(SWA/NFITMWA) ENDIF IF(NFITMWE.EQ.0) THEN WMWE=0.0 SWE=0.0 ELSE WMWE=SQRT(WMWE/NFITMWE) SWE=SQRT(SWE/NFITMWE) ENDIF IF(NFITIR.EQ.0) THEN WIR=0.0 SIR=0.0 ELSE WIR=SQRT(WIR/NFITIR) SIR=SQRT(SIR/NFITIR) ENDIF IF(NFIR10.EQ.0) THEN WIR10=0.0 SIR10=0.0 ELSE WIR10=SQRT(WIR10/NFIR10) SIR10=SQRT(SIR10/NFIR10) ENDIF IF(NFIR21.EQ.0) THEN WIR21=0.0 SIR21=0.0 ELSE WIR21=SQRT(WIR21/NFIR21) SIR21=SQRT(SIR21/NFIR21) ENDIF IF(NFIR32.EQ.0) THEN WIR32=0.0 SIR32=0.0 ELSE WIR32=SQRT(WIR32/NFIR32) SIR32=SQRT(SIR32/NFIR32) ENDIF IF(NFIMW0 .EQ.0) THEN WMW0=0.0 SMW0=0.0 ELSE WMW0=SQRT(WMW0/NFIMW0 ) SMW0=SQRT(SMW0/NFIMW0 ) ENDIF IF(NFI070 .EQ.0) THEN SW070=0.0 WIG070=0.0 ELSE SW070=SQRT(SW070/NFI070 ) WIG070=SQRT(WIG070/NFI070 ) ENDIF IF(NFIMW1 .EQ.0) THEN WMW1=0.0 SMW1=0.0 ELSE WMW1=SQRT(WMW1/NFIMW1 ) SMW1=SQRT(SMW1/NFIMW1 ) ENDIF IF(NFIMW2.EQ.0) THEN WMW2=0. SMW2=0. ELSE WMW2=SQRT(WMW2/NFIMW2) SMW2=SQRT(SMW2/NFIMW2) ENDIF IF(NFIMW3.EQ.0) THEN WMW3=0. SMW3=0. ELSE WMW3=SQRT(WMW3/NFIMW3) SMW3=SQRT(SMW3/NFIMW3) ENDIF IF(NFIMW4.EQ.0) THEN WMW4=0. SMW4=0. ELSE WMW4=SQRT(WMW4/NFIMW4) SMW4=SQRT(SMW4/NFIMW4) ENDIF IF(NFI004 .EQ.0) THEN WIG004=0. SIG004=0. ELSE WIG004=SQRT(WIG004/NFI004 ) SIG004=SQRT(SIG004/NFI004 ) ENDIF IF(NFI080 .EQ.0) THEN WIG080=0.0 SIG080=0.0 ELSE WIG080=SQRT(WIG080/NFI080 ) SIG080=SQRT(SIG080/NFI080 ) ENDIF IF(NFI010 .EQ.0) THEN WIG010=0.0 SIG010=0.0 ELSE WIG010=SQRT(WIG010/NFI010 ) SIG010=SQRT(SIG010/NFI010 ) ENDIF IF(NFI100 .EQ.0) THEN WIG100=0.0 SIG100=0.0 ELSE WIG100=SQRT(WIG100/NFI100 ) SIG100=SQRT(SIG100/NFI100 ) ENDIF IF(NFI020 .EQ.0) THEN WIG020=0.0 SIG020=0.0 ELSE WIG020=SQRT(WIG020/NFI020 ) SIG020=SQRT(SIG020/NFI020 ) ENDIF IF(NFI200 .EQ.0) THEN WIG200=0.0 SIG200=0.0 ELSE WIG200=SQRT(WIG200/NFI200 ) SIG200=SQRT(SIG200/NFI200 ) ENDIF IF(NFI045 .EQ.0) THEN WIG045=0.0 SIG045=0.0 ELSE WIG045=SQRT(WIG045/NFI045 ) SIG045=SQRT(SIG045/NFI045 ) ENDIF IF(NFI1M0 .EQ.0) THEN WIG1M0=0.0 SIG1M0=0.0 ELSE WIG1M0=SQRT(WIG1M0/NFI1M0 ) SIG1M0=SQRT(SIG1M0/NFI1M0 ) ENDIF IF(NGST0.EQ.0) THEN WGST0=0. SGST0=0. ELSE WGST0=SQRT(WGST0/NGST0) SGST0=SQRT(SGST0/NGST0) ENDIF WRITE(6,900) S,NOFIT WRITE(6,901) WMW,NFITMW WRITE(6,801) WMWA,NFITMWA WRITE(6,802) WMWE,NFITMWE SMW0=SMW0*29979.2458 SMW1=SMW1*29979.2458 SMW2=SMW2*29979.2458 SMW3=SMW3*29979.2458 SMW4=SMW4*29979.2458 SWA=SWA*29979.2458 SWE=SWE*29979.2458 WRITE(6,903) WMW0,NFIMW0 WRITE(6,904) WMW1,NFIMW1 WRITE(6,700) WMW2,NFIMW2 WRITE(6,600) WMW3,NFIMW3 WRITE(6,500) WMW4,NFIMW4 WRITE(6,902) WIR,NFITIR WRITE(6,912) WIR10,NFIR10 WRITE(6,913) WIR21,NFIR21 WRITE(6,914) WIR32,NFIR32 WRITE(6,899) WGST0,NGST0 WRITE(6,920) WIG004,NFI004 WRITE(6,931) WIG010,NFI010 WRITE(6,933) WIG020,NFI020 WRITE(6,935) WIG045,NFI045 WRITE(6,936) WIG070,NFI070 WRITE(6,930) WIG080,NFI080 WRITE(6,932) WIG100,NFI100 WRITE(6,934) WIG200,NFI200 WRITE(6,937) WIG1M0,NFI1M0 C SIR=SIR*0.000353553 C SMW0=SMW0*0.078 C SMW1=SMW1*0.04 CCC=29979.2458 SIG004=SIG004*CCC SIG080=SIG080*CCC SIG010=SIG010*CCC SIG100=SIG100*CCC SIG020=SIG020*CCC SIG200=SIG200*CCC SIG045=SIG045*CCC SIG1M0=SIG1M0*CCC SW070=SW070*CCC C SGST0=SGST0*0.0005 WRITE(6,906) SMW0 WRITE(6,907) SMW1 WRITE(6,701) SMW2 WRITE(6,601) SMW3 WRITE(6,501) SMW4 WRITE(6,803) SWA WRITE(6,804) SWE WRITE(6,905) SIR WRITE(6,915) SIR10 WRITE(6,916) SIR21 WRITE(6,917) SIR32 WRITE(6,888) SGST0 WRITE(6,890) SIG004 WRITE(6,892) SIG010 WRITE(6,894) SIG020 WRITE(6,896) SIG045 WRITE(6,897) SW070 WRITE(6,891) SIG080 WRITE(6,893) SIG100 WRITE(6,895) SIG200 WRITE(6,898) SIG1M0 900 FORMAT(/,' STD.DEV.(UNITLESS) = ',F20.8,' NB.DATA= ',I4) 901 FORMAT(/,' RMS.DEV.MW(UNITLESS) = ',F20.4,2X,I4) 801 FORMAT(/,' RMS.DEV.MW TYPE A(UNITLESS)= ',F20.4,2X,I4) 802 FORMAT(/,' RMS.DEV.MW TYPE E(UNITLESS)= ',F20.4,2X,I4) 902 FORMAT(/,' RMS.DEV.IR(UNITLESS) = ',F20.6,2X,I4) 912 FORMAT(/,' RMS.DEV.IR(UNITLESS) V=1-0 = ',F20.6,2X,I4) 913 FORMAT(/,' RMS.DEV.IR(UNITLESS) V=2-1 = ',F20.6,2X,I4) 914 FORMAT(/,' RMS.DEV.IR(UNITLESS) V=3-2 = ',F20.6,2X,I4) 903 FORMAT(/,' RMS.DEV.MW(UNITLESS)VT=0 = ',F20.4,2X,I4) 904 FORMAT(/,' RMS.DEV.MW(UNITLESS)VT=1 = ',F20.4,2X,I4) 700 FORMAT(/,' RMS.DEV.MW(UNITLESS)VT=2 = ',F20.4,2X,I4) 600 FORMAT(/,' RMS.DEV.MW(UNITLESS)VT=3 = ',F20.4,2X,I4) 500 FORMAT(/,' RMS.DEV.MW(UNITLESS)VT=4 = ',F20.4,2X,I4) 899 FORMAT(/,' RMS.DEV.GROUND STATE DIFF. VT=0) = ',F20.4,2X,I4) 920 FORMAT(/,' STD.MW. W=0.003MHZ',F20.6,2X,I4) 930 FORMAT(/,' STD.MW. W=0.080MHZ',F20.6,2X,I4) 931 FORMAT(/,' STD.MW. W=0.010MHZ',F20.6,2X,I4) 932 FORMAT(/,' RMS.MW. W=0.100MHZ',F20.6,2X,I4) 933 FORMAT(/,' RMS.MW. W=0.020MHZ',F20.6,2X,I4) 934 FORMAT(/,' RMS.MW. W=0.200MHZ',F20.6,2X,I4) 935 FORMAT(/,' RMS.MW. W=0.045MHZ',F20.6,2X,I4) 936 FORMAT(/,' RMS.MW. W=0.070MHZ',F20.6,2X,I4) 937 FORMAT(/,' RMS.MW. W=1.000MHZ',F20.6,2X,I4) 905 FORMAT(/,' RMS.DEV.IR(CM-1) = ',F20.6) 915 FORMAT(/,' RMS.DEV.IR(CM-1) V=1-0 = ',F20.6) 916 FORMAT(/,' RMS.DEV.IR(CM-1) V=2-1 = ',F20.6) 917 FORMAT(/,' RMS.DEV.IR(CM-1) V=3-2 = ',F20.6) 906 FORMAT(/,' RMS.DEV.MICRONDE VT=0(MHZ)= ',F20.4) 907 FORMAT(/,' RMS.DEV.MICRONDE VT=1(MHZ)= ',F20.4) 701 FORMAT(/,' RMS.DEV.MICRONDE VT=2(MHZ)= ',F20.4) 601 FORMAT(/,' RMS.DEV.MICRONDE VT=3(MHZ)= ',F20.4) 501 FORMAT(/,' RMS.DEV.MICRONDE VT=4(MHZ)= ',F20.4) 803 FORMAT(/,' RMS.DEV.MICRONDE "A" (MHZ)= ',F20.4) 804 FORMAT(/,' RMS.DEV.MICRONDE "E" (MHZ)= ',F20.4) 888 FORMAT(/,' RMS.DEV.GROUND STATE DIFF.VT=0 (CM-1)) = ',F20.6) 890 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.003MHZ',F20.4,' MHZ') 891 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.080MHZ',F20.4,' MHZ') 892 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.010MHZ',F20.4,' MHZ') 893 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.100MHZ',F20.4,' MHZ') 894 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.020MHZ',F20.4,' MHZ') 895 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.200MHZ',F20.4,' MHZ') 896 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.045MHZ',F20.4,' MHZ') 897 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=0.070MHZ',F20.4,' MHZ') 898 FORMAT(/,' RMS.DEV.MICROWAVE,WEIGHT=1.000MHZ',F20.4,' MHZ') RETURN END C************************************************************* SUBROUTINE MATOU2(X,Z,M,N,MM,NN) IMPLICIT REAL*8 ( A-H,O-Z ) C C ****************************************************************** DIMENSION X(M,N),Z(M) 1004 FORMAT(/) 1003 FORMAT(4X,9F14.7) 1000 FORMAT(9(11X,I3)) 1002 FORMAT(I4,9F14.6) C IFLG=1 ILOWER=1 1 IUPPER=ILOWER+8 IF(IUPPER-NN)3,2,2 2 IUPPER=NN IFLG=0 3 WRITE(6,1000)(J,J=ILOWER,IUPPER) WRITE(6,1004) WRITE(6,1003) (Z(J),J=ILOWER,IUPPER) WRITE(6,1004) DO 4 I=1,MM 4 WRITE(6,1002)I,(X(I,J),J=ILOWER,IUPPER) WRITE(6,1004) IF(IFLG)5,6,5 5 ILOWER=ILOWER+9 GO TO 1 6 RETURN END C************************************************************ C************************************************************* SUBROUTINE MATOUT(X,M,N,MM,NN) IMPLICIT REAL*8 ( A-H,O-Z ) C C ****************************************************************** DIMENSION X(M,N) 1004 FORMAT(/) 1000 FORMAT(9(11X,I3)) 1002 FORMAT(I4,9F14.6) C IFLG=1 ILOWER=1 1 IUPPER=ILOWER+8 IF(IUPPER-NN)3,2,2 2 IUPPER=NN IFLG=0 3 WRITE(6,1000)(J,J=ILOWER,IUPPER) WRITE(6,1004) DO 4 I=1,MM WRITE(6,1002)I,(X(I,J),J=ILOWER,IUPPER) 4 CONTINUE WRITE(6,1004) IF(IFLG)5,6,5 5 ILOWER=ILOWER+9 GO TO 1 6 RETURN END C************************************************************ SUBROUTINE PARAM(I,CST) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************ C READING OF THE PARAMETERS TO BE FITTED AND STORING THEM IN C X VECTOR(ONLY ONE VECTOR FOR THE 2 STATES) PARA(LO,L) C RETAINS WHICH PARAMETER C CHARACTER ABC*702 CHARACTER*6 PARA,REF CHARACTER*1 PR,PROBS COMMON/PAR/IND(2,80),NP(2),X(160),ITO COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/SYMB/ABC DIMENSION CST(2,117) READ(5,20) NP(I) 20 FORMAT(I3) WRITE(6,*) ' NBRE OF PARAMETERS TO FIT ',NP(I) IF(NP(I).NE.0) THEN K=0 IF(I.EQ.2) K=NP(1) WRITE(6,*) ' PARAMETER TO BE FIT ' DO 2 J=1,NP(I) IF(I.EQ.2.AND.NP(1).NE.0) THEN JJ=J+NP(1) ELSE JJ=J ENDIF READ(5,'(A6)') PARA(I,JJ) WRITE(6,*) ' PARA = ',PARA(I,JJ) IK=IREP(PARA(I,JJ)) IND(I,J)=IK X(K+J)=CST(I,IK) 2 CONTINUE ENDIF C PRINT*,' NP(1)',NP(1),' NP(2)= ',NP(2) RETURN END SUBROUTINE REJEC(REJ) IMPLICIT REAL*8 ( A-H,O-Z ) C C THIS SUBROUTINE REJECTS AN EXPERIMENTAL VALUE IF THE C CORRESPONDING TME IS GREATER THAN 3*SIGMA (S). C not used in the present version. C ********************************************************* CHARACTER*6 PARA,REF CHARACTER*1 PR,PROBS INTEGER REJ COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/PAR/IND(2,80),NP(2),X(160),ITO COMMON/EXPDAT / ETRANS(80000),IVOBS(2,80000), * NOBS(2,80000),KOBS(2,80000),IBLK(2,80000),IW(80000),W(80000), * NDATA,NADAT,EPSI,NOFIT,NFITMW,NFITIR,NFIMW0 ,NFIMW1 ,NITT, * IAST(80000),NFI004 ,NFI080 ,NFI100 ,NFI020 ,NFI200 ,NFI045 , & NGST0 ,NFI010 ,SW070,NFI070 ,NFIMW2,IKAKC,JMAX,NFI1M0, & NFIR10,NFIR21,NFIR32,NFIMW3,NFIMW4,SWA,SWE,NFITMWA,NFITMWE COMMON/CNTL/NIT,IEND,DELTAS,S,TME(80000),TERMST,NP12,SMW,SIR, 1 SMW0,SMW1,SIG004,SIG080,SIG10,SIG100,SIG020,SIG200,SIG045, & SGST0,SMW2,SIG1M0,SMW3,SIG32,SMW4 COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) C REJLIM = 3*S WRITE(6,803) REJLIM IREJL = 10000 * REJLIM DO 88 I = 1, NDATA IF (IAST(I) .EQ. 0) THEN ITME = 10000 * ABS(TME(I)) IF (ITME .GT. IREJL) THEN IAST(I) = 2 WRITE(6,800) I,ETRANS(I) REJ = 1 END IF END IF 88 CONTINUE 800 FORMAT(' LINE NO ',I5,' AT',F14.4,' HAS BEEN REJECTED') 803 FORMAT(' THE REJECTION LIMIT IS ',F10.4) RETURN END C*********************************************************** SUBROUTINE RESUM1(LO,EGVL) IMPLICIT REAL*8 ( A-H,O-Z ) C********************************************************** C WRITES THE ENERGIES OF TORSION LEVELS C C CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) COMMON/ROTOR/NROTOR,NDIMTO DIMENSION CST(2,117),EGVL(NDMX) DO 1 I=1,NDIMTO WRITE(6,53) I-1,EGVL(I) 53 FORMAT(/,' VT= ',I3,' E= ',F16.10,' CM-1') 1 CONTINUE RETURN END C*********************************************************** SUBROUTINE RESUM2(LO,EGVL,EGVC) IMPLICIT REAL*8 ( A-H,O-Z ) C*********************************************************** C WRITES THE ENERGIES OF THE ROTATION-TORSION LEVELS(E STATES) C CHARACTER*1 PUR(3*(2*(30)+1)) CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION IPST(3*(2*(30)+1)),EGVL(NDMX),EGVC(NDMX,NDMX) WRITE(6,*)' N= ',N,' ISIG= ',ISIG WRITE(6,50) 50 FORMAT(/,' TYPE E ',/) NN=2*N+1 NNT=NN*4 DO 1 J=1,NNT CST=0.0 DO 2 I=1,NNT C=EGVC(I,J)**2 IF(C.GT.CST) THEN CST=C IPST(J)=I ENDIF 2 CONTINUE C TEST THE PURITY PUR(J)=' ' DO 4 K=1,NNT IF(K.EQ.IPST(J)) GO TO 4 CC=EGVC(K,J)**2 IF(CC.GE.(CST*0.7D0)) THEN PUR(J)='*' GO TO 5 ENDIF 4 CONTINUE 5 IVT=0 DO 3 KKK=1,(04) IN=NN*KKK IF (IPST(J).GT.IN) IVT=IVT+1 3 CONTINUE K=IPST(J)-(N+1)-IVT*NN EGVL(J)=EGVL(J)+EVIB(LO) WRITE(6,55) IVT,K,EGVL(J),N,PUR(J) WRITE(LO+32,55) IVT,K,EGVL(J),N,PUR(J) 55 FORMAT (' VT= ',I3,' K= ',I3,' E= ',F16.10,' CM-1 ',' N= ',I3, & ' PURITY= ',A1) 1 CONTINUE RETURN END C************************************************************* SUBROUTINE RESUM3(LO,EGVL,EGVC) IMPLICIT REAL*8 ( A-H,O-Z ) C************************************************************** C WRITES THE ENERGIES OF THE ROT-TORSION LEVELS (A STATES) C CHARACTER PAR*1 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/QUANT/ISIG,IV,N,KK COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION EGVL(NDMX),EGVC(NDMX,NDMX) WRITE(6,*) ' N= ',N,' ISIG= ',ISIG WRITE(6,50) 50 FORMAT(/,' TYPE A ',/) NN=2*N+1 NNT=NN*4 DO 1 J=1,NNT IVT=0 DO 3 KKK=1,(04) IN=NN*KKK IF (J.GT.IN) IVT=IVT+1 3 CONTINUE K=(N*J+N-IVT*N*NN)/NN IF (K.EQ.0) THEN PAR= '+' ELSE IF (MOD(K,2).NE.0) THEN IF(MOD(IVT,2).EQ.0) THEN IF(MOD(J,2).EQ.0) THEN PAR='+' ELSE PAR='-' ENDIF ELSE IF(MOD(J,2).EQ.0) THEN PAR='-' ELSE PAR='+' ENDIF ENDIF ELSE IF(MOD(K,2).EQ.0) THEN IF(MOD(IVT,2).EQ.0) THEN IF(MOD(J,2).EQ.0) THEN PAR='-' ELSE PAR='+' ENDIF ELSE IF(MOD(J,2).EQ.0) THEN PAR='+' ELSE PAR='-' ENDIF ENDIF ENDIF EGVL(J)=EGVL(J)+EVIB(LO) WRITE(6,56) IVT,K,EGVL(J),N,PAR WRITE(LO+30,56) IVT,K,EGVL(J),N,PAR 56 FORMAT (' VT= ',I3,' K= ',I3,' E= ',F16.10,' CM-1 ',' N= ',I3, & ' PAR= ',A1) 1 CONTINUE RETURN END SUBROUTINE ROTTOR(CST,LO,EGVL,EGVC,ETOR,A,nitt) IMPLICIT REAL*8 ( A-H,O-Z ) C*********************************************************** C CALCULATES THE ENERGIES OF LEVELS AND TRANSITIONS FOR EACH C N,LO BLOCK.SAME SUBROUTINE AS THE ROTOR PROGRAM C C CHARACTER ABC*702 CHARACTER*1 PR,PROBS CHARACTER*6 PARA,REF PARAMETER(NDMX=(2*(30)+1)*9,KDMX=2*(30)+1,NTORMX=2*(10)+1) COMMON/EIGEN/H(NDMX,NDMX),WORK(NDMX),EGVC3(KDMX,9),NBR(2) COMMON/DAT/CONV,IBUG1,IBUG2,IBUG3,IBUG4,KTRONC,KTRON1 , * IBGTME,IBGVTD ,EVIB(2),IRMW,TERM COMMON/ROTOR/NROTOR,NDIMTO COMMON/SYMB/ABC COMMON/QUANT/ISIG,IV,N,KK COMMON/QUANT2/KKK(NDMX),NVT(NDMX),IPRITY(NDMX) COMMON/C2C/PERC2L,PER23L(NDMX) COMMON/CHARA/PARA(2,160),PR,PROBS(2,80000),REF(80000) DIMENSION CST(2,117),UNIT(160,160),WKSPCE(160),EGVL(NDMX) DIMENSION EGVC(NDMX,NDMX),ETOR(KDMX,9),A(9,NTORMX,KDMX) DIMENSION EGVC4(KDMX,NTORMX,NDMX) dimension oldegv(ndmx,ndmx),oldegv1(ndmx,ndmx) dimension oldegv2(ndmx,ndmx),oldegv3(ndmx,ndmx) real*8 sum123(0:8),evect(ndmx),sum1(9) integer irespos(1) C KTRON1 =KTRONC+1 NDIMTO =2*KTRONC+1 if(nitt.eq.0)then jmax=20 c overlaps for torsional eigen functions if(n.eq.0.and.isig.eq.0)then c open(unit=13,file='ators.txt',access='sequential',status='unknown') c overlaps for A-type levels DO K1=0,jmax k=-k1 CALL HTORSE (K,CST,LO) CALL TRED2(NDIMTO,NDMX,H,EGVL,WORK,EGVC) CALL TQL2(NDIMTO,NDMX,EGVL,WORK,EGVC,IERR) c cycle over eigenvalues for a given k value do iji=1,ndimto evect=0.d0 c torsional energy c write(13,1924)k,egvl(iji),iji,isig do iki=1,ndimto evect(iki)=egvc(iki,iji)**2 enddo c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c mixing ratio with the same J c for A-type calculate mixing ratio with the same pr23l evect=0.d0 c cycle over all eigen vectors for previous N do i=1,ndimto if(i.eq.iji)cycle c calculating of overlap integral do iki=1,nrotor evect(i)=evect(i)+dabs(egvc(iki,iji)*egvc(iki,i)) enddo enddo irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) c overlap integrals with previous K evect=0.d0 if(k.lt.0)then c cycle over all eigen vectors for previous N do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output two largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-2 evect=0.d0 if(k.lt.-1)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv1(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-3 evect=0.d0 if(k.lt.-2)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv2(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-4 evect=0.d0 if(k.lt.-3)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv3(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo write(13,1917) enddo if(k.lt.-2)then c clean the store for n-4 eigen vectors and put there new ones oldegv3=0.d0 do j=1,ndimto do iji=1,ndimto oldegv3(iji,j)=oldegv2(iji,j) enddo enddo endif if(k.lt.-1)then c clean the store for n-3 eigen vectors and put there new ones oldegv2=0.d0 do j=1,ndimto do iji=1,ndimto oldegv2(iji,j)=oldegv1(iji,j) enddo enddo endif if(k.lt.0)then c clean the store for n-2 eigen vectors and put there new ones oldegv1=0.d0 do j=1,ndimto do iji=1,ndimto oldegv1(iji,j)=oldegv(iji,j) enddo enddo endif c clean the store for n-1 eigen vectors and put there new ones oldegv=0.d0 do j=1,ndimto do iji=1,ndimto oldegv(iji,j)=egvc(iji,j) enddo enddo 1924 format(1x,I6,F15.8,I4,I4,$) enddo DO K=0,jmax CALL HTORSE (K,CST,LO) CALL TRED2(NDIMTO,NDMX,H,EGVL,WORK,EGVC) CALL TQL2(NDIMTO,NDMX,EGVL,WORK,EGVC,IERR) c cycle over eigenvalues for a given k value do iji=1,ndimto evect=0.d0 c torsional energy write(13,1924)k,egvl(iji),iji,isig do iki=1,ndimto evect(iki)=egvc(iki,iji)**2 enddo c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c mixing ratio with the same J c for A-type calculate mixing ratio with the same pr23l evect=0.d0 c cycle over all eigen vectors for previous N do i=1,ndimto if(i.eq.iji)cycle c calculating of overlap integral do iki=1,nrotor evect(i)=evect(i)+dabs(egvc(iki,iji)*egvc(iki,i)) enddo enddo irespos = MAXLOC(evect) write(13,1916)irespos(1),evect(irespos(1)) c overlap integrals with previous K evect=0.d0 if(k.gt.0)then c cycle over all eigen vectors for previous N do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output two largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-2 evect=0.d0 if(k.gt.1)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv1(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-3 evect=0.d0 if(k.gt.2)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv2(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-4 evect=0.d0 if(k.gt.3)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv3(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(13,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c write(13,1917) enddo if(k.gt.2)then c clean the store for n-4 eigen vectors and put there new ones oldegv3=0.d0 do j=1,ndimto do iji=1,ndimto oldegv3(iji,j)=oldegv2(iji,j) enddo enddo endif if(k.gt.1)then c clean the store for n-3 eigen vectors and put there new ones oldegv2=0.d0 do j=1,ndimto do iji=1,ndimto oldegv2(iji,j)=oldegv1(iji,j) enddo enddo endif if(k.gt.0)then c clean the store for n-2 eigen vectors and put there new ones oldegv1=0.d0 do j=1,ndimto do iji=1,ndimto oldegv1(iji,j)=oldegv(iji,j) enddo enddo endif c clean the store for n-1 eigen vectors and put there new ones oldegv=0.d0 do j=1,ndimto do iji=1,ndimto oldegv(iji,j)=egvc(iji,j) enddo enddo enddo c close(13) endif if(n.eq.0.and.isig.eq.1)then c open(unit=14,file='etors.txt',access='sequential',status='unknown') c overlaps for E-type levels DO K1=0,jmax k=-k1 CALL HTORSE (K,CST,LO) CALL TRED2(NDIMTO,NDMX,H,EGVL,WORK,EGVC) CALL TQL2(NDIMTO,NDMX,EGVL,WORK,EGVC,IERR) c cycle over eigenvalues for a given k value do iji=1,ndimto evect=0.d0 c torsional energy c write(14,1924)k,egvl(iji),iji,isig do iki=1,ndimto evect(iki)=egvc(iki,iji)**2 enddo c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c mixing ratio with the same J c for A-type calculate mixing ratio with the same pr23l evect=0.d0 c cycle over all eigen vectors for previous N do i=1,ndimto if(i.eq.iji)cycle c calculating of overlap integral do iki=1,nrotor evect(i)=evect(i)+dabs(egvc(iki,iji)*egvc(iki,i)) enddo enddo irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) c overlap integrals with previous K evect=0.d0 if(k.lt.0)then c cycle over all eigen vectors for previous N do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output two largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-2 evect=0.d0 if(k.lt.-1)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv1(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-3 evect=0.d0 if(k.lt.-2)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv2(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-4 evect=0.d0 if(k.lt.-3)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv3(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c write(14,1917) enddo if(k.lt.-2)then c clean the store for n-4 eigen vectors and put there new ones oldegv3=0.d0 do j=1,ndimto do iji=1,ndimto oldegv3(iji,j)=oldegv2(iji,j) enddo enddo endif if(k.lt.-1)then c clean the store for n-3 eigen vectors and put there new ones oldegv2=0.d0 do j=1,ndimto do iji=1,ndimto oldegv2(iji,j)=oldegv1(iji,j) enddo enddo endif if(k.lt.0)then c clean the store for n-2 eigen vectors and put there new ones oldegv1=0.d0 do j=1,ndimto do iji=1,ndimto oldegv1(iji,j)=oldegv(iji,j) enddo enddo endif c clean the store for n-1 eigen vectors and put there new ones oldegv=0.d0 do j=1,ndimto do iji=1,ndimto oldegv(iji,j)=egvc(iji,j) enddo enddo enddo DO K=0,jmax CALL HTORSE (K,CST,LO) CALL TRED2(NDIMTO,NDMX,H,EGVL,WORK,EGVC) CALL TQL2(NDIMTO,NDMX,EGVL,WORK,EGVC,IERR) c cycle over eigenvalues for a given k value do iji=1,ndimto evect=0.d0 c torsional energy c write(14,1924)k,egvl(iji),iji,isig do iki=1,ndimto evect(iki)=egvc(iki,iji)**2 enddo c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c mixing ratio with the same J c for A-type calculate mixing ratio with the same pr23l evect=0.d0 c cycle over all eigen vectors for previous N do i=1,ndimto if(i.eq.iji)cycle c calculating of overlap integral do iki=1,nrotor evect(i)=evect(i)+dabs(egvc(iki,iji)*egvc(iki,i)) enddo enddo irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) c overlap integrals with previous K evect=0.d0 if(k.gt.0)then c cycle over all eigen vectors for previous N do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output two largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-2 evect=0.d0 if(k.gt.1)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv1(iki,i) enddo evect(i)=dabs(evect(i)) enddo endif c output largest overlap integrals c cycle over 2 maximum values of components do iki=1,2 irespos = MAXLOC(evect) c write(14,1916)irespos(1),evect(irespos(1)) evect(irespos(1))=0.d0 enddo c overlap integrals with J-3 evect=0.d0 if(k.gt.2)then c cycle over all eigen vectors for N-2 do i=1,ndimto c calculating of overlap integral do iki=1,ndimto evect(i)=evect(i)+egvc(iki,iji)*oldegv2(iki,i) enddo