C============================================================================== C C X I A M C C Rotational, Centrifugal Distortion, Internal Rotation Calculation (V2.5e) C Holger Hartwig 08-Nov-96 (hartwig@phc.uni-kiel.de) C C Please cite: H.Hartwig and H.Dreizler, Z.Naturforsch, 51a (1996) 923. C C============================================================================== C C This file collects the constituent modules for XIAM in the order: C C iam iamio iamint iamadj iamm iamv iamv2 iamfit mgetx iamsys iamlib C C Three remaining INCLUDEd modules: C C iam_.for, iamdata_.for, mgetx_.for C C are in separate files and have to be available on compilation. C Simplest compilation command with Compaq Visual Fortran is: C C df -fast -static xiamall C zk C------------------------------------------------------------------------------ C C module IAM.FOR C C------------------------------------------------------------------------------ C program xiam C Fitting of Rotational Constants and Centrifugal Distortion Constants C and Internal-rotation parameters of a 3 fold top C Syntax of Input-File: C C Phase Convention: =0.5*dsqrt(J*(J+1)-K*(K+1)) C implicit none include 'iam_.for' C arrays for fit-proc C DIMFIT: Dimension for paramaters to be fitted simultaneouly real*8 covar(DIMFIT,DIMFIT), alpha(DIMFIT,DIMFIT), $ evec(DIMFIT,DIMFIT), w(DIMFIT), $ beta(DIMFIT), freed(DIMFIT) C One fit paramater is a linearcombination of the parameters a C PArameter Linear Combinations as fit parameters C DIMPLC maximum number parameters in one fitparameters real*8 palc(DIMFIT,-1:DIMPLC) integer pali(DIMFIT, 0:DIMPLC,2) C dfit controls the differences quotient calculation integer dfit(DIMFIT) C a is the parameter array for the Hamiltonian real*8 a(DIMPAR,DIMVB), anew(DIMPAR,DIMVB), da(DIMPAR,DIMVB) C ifit controls if analytic derivatives are to be calculated integer ifit(DIMPAR,DIMVB) C real*8 ab(DIMPAR) real*8 chi2(DIMVB,DIMVV,2),wght(DIMVB,DIMVV,2) integer ndat(DIMVB,DIMVV,2) real*8 sigsq,sigsqold,stepw,cmax integer npar,i,izyk,j,iconv,iq,nsvfit,fstat,ift,nfit,itop,deriv integer myints,oib,maxf,maxb,maxg,ifunpr,ming integer ij,ib,iv,is,jc,ic,i2f,isf real*8 maxchn character*15 astr character*12 dstr character*12 binfname logical sortend real*8 pi,indeg,inkj,inkc,incm,fold,condno integer myand,myor external myand,myor integer sig_stat common/sig_com/sig_stat include 'iamdata_.for' call mysignal() pi=dacos(-1.0d0) inkj=3.9903132D-04 inkc=3.9903132D-04/4.186903D0 indeg=180.0d0/pi incm=1.0d0/29.97925d0 write(*,'(A,A)') $ ' Rotational, Centrifugal Distortion,' $ ,' Internal Rotation Calculation (V2.5e)' write(*,'(23X,A)') $ 'Holger Hartwig 08-Nov-96 (hartwig@phc.uni-kiel.de)' write(*,'(/,2A,/)') ' Please cite:', $ ' H.Hartwig and H.Dreizler, Z.Naturforsch, 51a (1996) 923.' write(*,'(A,$)') ' Calculation date and time: ' call mydate() write(*,'(A)') ' Type help now for the list of parameters : ' call parinp(a,palc,pali,ifit,dfit,npar,nfit) call mysignal() C the array todo is the list of (good) quantum no.s which C have to be calculated to assign all transitions C the good quantum no.s are j, f, b and sigma, no symmetry adapted k used size(S_MAXK)=0 ncalc=0 maxf=0 maxg=0 maxb=0 ming=100 C qlin is the list of transitions do ic=1, DIMTDO todo(ic,Q_STAT)=0 end do C copy qlin into todo do i=1, ctlint(C_NDATA) do iq=Q_UP, Q_LO if (maxg.lt.qlin(i,Q_S,iq)) maxg=qlin(i,Q_S,iq) if (ming.gt.qlin(i,Q_S,iq)) ming=qlin(i,Q_S,iq) if (maxb.lt.qlin(i,Q_B,iq)) maxb=qlin(i,Q_B,iq) if (maxf.lt.qlin(i,Q_F,iq)) maxf=qlin(i,Q_F,iq) c vb(qlin(i,Q_B,iq))=qlin(i,Q_V1,iq) if (size(S_MAXK).lt.qlin(i,Q_J,iq)) $ size(S_MAXK)=qlin(i,Q_J,iq) do ic=1, ncalc if ((qlin(i,Q_J,iq).eq.todo(ic,Q_J)) $ .and.(qlin(i,Q_S,iq).eq.todo(ic,Q_S)) $ .and.(qlin(i,Q_B,iq).eq.todo(ic,Q_B)) $ .and.(qlin(i,Q_F,iq).eq.todo(ic,Q_F))) then goto 50 end if end do ncalc=ncalc+1 if (ncalc.gt.DIMTDO) stop 'ERROR: todo >DIMTDO' todo(ncalc,Q_J)=qlin(i,Q_J,iq) todo(ncalc,Q_S)=qlin(i,Q_S,iq) todo(ncalc,Q_B)=qlin(i,Q_B,iq) todo(ncalc,Q_F)=qlin(i,Q_F,iq) ic=ncalc 50 continue if ((dln(i,LN_ERR).ne.NOFIT).or.(ctlint(C_DFRQ).gt.0)) $ todo(ic,Q_STAT)=myor(todo(ic,Q_STAT),2) end do end do C for intensities all j,b,f,gam are needed if (ctlint(C_INTS).ge.2) then write(*,'(/,A)') ' Intensity Calculation Mode ' if (ctlint(C_NZYK).ne.1) then ctlint(C_NZYK)=1 write(*,'(A)') ' Number of iteration cylcles reset to one !' end if c if (ctlint(C_NTOP).gt.0) ming=min(1,ming) do ib=1, size(S_NB) ! maxb do ij=0,size(S_MAXK) C because nq-hfs intensities are not implemented: ignore f C do i2f= i2f=-1 isf=0 if (ctlint(C_INTS).eq.2) isf=1 do is=isf, size(S_G) do ic=1, ncalc if ( (ij.eq.todo(ic,Q_J)) $ .and.(is.eq.todo(ic,Q_S)) $ .and.(ib.eq.todo(ic,Q_B)) $ .and.(i2f.eq.todo(ic,Q_F))) goto 51 end do ncalc=ncalc+1 if (ncalc.gt.DIMTDO) stop 'ERROR: todo >DIMTDO' todo(ncalc,Q_J)=ij todo(ncalc,Q_S)=is todo(ncalc,Q_B)=ib todo(ncalc,Q_F)=i2f 51 continue end do C end do end do end do end if C sort b values in first place, followed by ascending J values do iq=1, ncalc sortend=.true. do i=2, ncalc if (todo(i-1,Q_B).gt.todo(i,Q_B)) then call swptdo(i,i-1) sortend=.false. end if end do if (sortend) goto 13 end do 13 continue do iq=1, ncalc sortend=.true. do i=2, ncalc if ((todo(i-1,Q_J).gt.todo(i,Q_J)) $ .and.(todo(i-1,Q_B).eq.todo(i,Q_B))) then call swptdo(i,i-1) sortend=.false. end if end do if (sortend) goto 12 end do 12 continue do iq=1, ncalc sortend=.true. do i=2, ncalc if ((todo(i-1,Q_S).gt.todo(i,Q_S)) $ .and.(todo(i-1,Q_J).eq.todo(i,Q_J)) $ .and.(todo(i-1,Q_B).eq.todo(i,Q_B))) then call swptdo(i,i-1) sortend=.false. end if end do if (sortend) goto 11 end do 11 continue write(*,'(/,2X,A,I4)') '\\ Maximal K = J =',size(S_MAXK) if (myand(ctlint(C_PRI),AP_ST).ne.0) then write(*,'(/,A)') ' ToDo: list of quantum no.' write(*,'(3X,A)') ' J Sym B F Stat' do i=1,ncalc write(*,'(3X,5I4)') todo(i,Q_J),todo(i,Q_S) $ ,todo(i,Q_B),todo(i,Q_F),todo(i,Q_STAT) end do end if C -------------------------------------------------------------- oib =-1 do ic=1, ncalc ib =todo(ic,Q_B) if (ib.ne.oib) then oib=ib ctlint(C_WOODS)=ctlnb(CB_WDS,ib) ctlint(C_ADJF) =ctlnb(CB_ADJ,ib) write(*,'(2X,2(A,I3))') '\\ B=',ib,' adj=',ctlint(C_ADJF) if (ctlint(C_NTOP).le.0) goto 33 if (myand(ctlint(C_WOODS),1).ne.0) write(*,'(2X,A)') $ '\\ (1) calculate torsional integrals ' if (myand(ctlint(C_WOODS),2).ne.0) write(*,'(2X,A)') $ '\\ (2) use scaled torsional integrals ' if (myand(ctlint(C_WOODS),4).ne.0) write(*,'(2X,A)') $ '\\ (4) use torsional integrals in the rotation matrix ' if (myand(ctlint(C_WOODS),8).ne.0) write(*,'(2X,A)') $ '\\ (8) use torsional integrals in H_ir for other top(s)' if (myand(ctlint(C_WOODS),16).ne.0) write(*,'(2X,A)') $ '\\ (16)use torsional integrals in twotop terms' if (myand(ctlint(C_WOODS),32).ne.0) write(*,'(2X,A)') $ '\\ (32)use torsional integrals in rigid rotor H_rr ' if (myand(ctlint(C_WOODS),64).ne.0) write(*,'(2X,A)') $ '\\ (64)multiply torsional integrals like Demaison' 33 continue end if end do do i=1, DIMPAR do ib=1, size(S_NB) anew(i,ib)=a(i,ib) end do end do sigsqold=1.0d30 iconv=1 deriv=1 fstat=0 stepw=1.0d0 myints=ctlint(C_INTS) ctlint(C_INTS)=0 do izyk=1, ctlint(C_NZYK) ifunpr=0 if (ctlint(C_EVAL).ne.0) write(20,'(/,A,i5)') ' Iteration',izyk if (izyk.ne.1) write(*,'(/,A,/,A,i3)') $ ' -------------------------------------------------------', $ ' Iteration :', izyk-1 xcy=0 if (myand(ctlint(C_XPR),XP_EC).ne.0) xcy=1 if ((izyk.eq.1).and.(myand(ctlint(C_XPR),XP_FI).ne.0)) xcy=1 if ((iconv.ge.2).and.(myand(ctlint(C_XPR),XP_LA).ne.0)) xcy=1 if ((myand(ctlint(C_PRI),AP_PC).ne.0).and.(xcy.ge.1)) $ call prpar(a) C -- calc. the spectrum if (ctlint(C_NZYK).eq.1) ctlint(C_INTS)=myints if (iconv.eq.2) ctlint(C_INTS)=myints call calc1(anew,ifit,dfit,npar,nfit,palc,pali,0) if (ctlint(C_NZYK).eq.1) goto 40 C -- get sigsq fstat=0 call lmfit(ctlint(C_NDATA),npar,size(S_NB) $ ,DIMFIT,DIMPAR,DIMVB,DIMPLC $ ,ctlint(C_PRI),nsvfit,nfit,ifit,dfit $ ,alpha,covar,evec,beta,w,a,anew,da,freed,palc,pali $ ,ctlpar(C_ROFIT),sigsq,sigsqold,ctlpar(C_EPS) $ ,stepw,fstat,ctlint(C_ORGER),0 $ ,ctlpar(C_LMBDA),ctlint(C_FITSC),ctlint(C_SVDER)) C -- test if sigsq was improved if ((sigsq/sigsqold).le.ctlpar(C_CNVG)) iconv=1 if (((sigsq/sigsqold).gt.ctlpar(C_CNVG)) $ .and.((sigsq/sigsqold).le.(2.0d0-ctlpar(C_CNVG)))) $ iconv=iconv+1 C -- switch to better derivatives if the fit is not converging if ((sigsq/sigsqold).gt.(2.0d0-ctlpar(C_CNVG))) then if (iconv.gt.0) iconv=0 if (iconv.le.0) iconv=iconv-1 end if C-- set deriv to 3 at the the first time; deriv is set to 2 later on. if ((iconv.eq.-3).and.(deriv.eq.1)) then deriv=3 write(*,*)'Switching to better derivatives (takes more time)' end if if ((sigsq.lt.sigsqold).and.(myand(ctlint(C_XPR),XP_CC).ne.0)) $ xcy=1 write(*,'(A,D12.6,A,F8.6,A,I2)') $ ' Sigma:',sqrt(sigsq) $ ,' Sigma/OldSigma:',sqrt(sigsq/sigsqold) $ ,' conv:',iconv C -- print the parameters if anew is better than before if ((izyk.ne.1).and.(xcy.ge.1)) then write(*,'(/,15X,A,5X,A)') 'Parameters','Change' call pra(anew,da,ifit,myand(ctlint(C_PRI),AP_PL),1,0) end if C -- print the parameter with the max change if (izyk.ne.1) then maxchn=0.0 ic=0 do i=1, npar do ib=1, size(S_NB) if ((anew(i,ib).ne.0.0).and.(ifit(i,ib).ne.0)) then if (abs(da(i,ib)/anew(i,ib)).gt.maxchn) then maxchn=abs(da(i,ib)/anew(i,ib)) ic=i ij=ib end if end if end do end do if (ic.gt.0) then call pra_f(anew(ic,ij),da(ic,ij),astr,dstr) write(*,'(1X,2A,I1,4A,F8.3,A)') parstr(ic),'(',ij,')' $ ,astr,' ',dstr,maxchn*100.0,'% Max. Change' end if end if C -- print the transitions if (((myand(ctlint(C_PRI),AP_TL).ne.0).and.(xcy.ge.1)).or. $ ((myand(ctlint(C_PRI),AP_TF).ne.0).and.(izyk.eq.1))) then write(*,*) ifunpr=1 call funpr(ctlpar(C_ROFIT),ndat,chi2,wght) end if C -- if sigsq is better then oldsigsq calc the dqu derivatives if ((sigsq.lt.sigsqold).or.(ctlint(C_DFRQ).gt.1) $ .or.(deriv.eq.3)) then call calc1(anew,ifit,dfit,npar,nfit,palc,pali,deriv) if (ctlint(C_DFRQ).ne.0) then write(21,'(/,A,i5)') ' Iteration',izyk call prderiv(dfit,nfit) end if end if if (deriv.eq.3) deriv=2 C -- if convergence is achived: fstat=-1 to calc. the errors fstat=1 if (iconv.gt.2) fstat=-1 if (izyk.eq.ctlint(C_NZYK)) fstat=-1001 if (sig_stat.eq.1) fstat=-1002 if (ctlpar(C_LMBDA).gt.1.0d3) fstat=-1003 call lmfit(ctlint(C_NDATA),npar,size(S_NB) $ ,DIMFIT,DIMPAR,DIMVB,DIMPLC $ ,ctlint(C_PRI) $ ,nsvfit,nfit,ifit,dfit $ ,alpha,covar,evec,beta,w,a,anew,da,freed,palc,pali $ ,ctlpar(C_ROFIT),sigsq,sigsqold,ctlpar(C_EPS) $ ,stepw,fstat,ctlint(C_ORGER),0 $ ,ctlpar(C_LMBDA),ctlint(C_FITSC),ctlint(C_SVDER)) if (fstat.ge.0) then do ib=1, size(S_NB) call adjusta(anew(1,ib),npar,ctlnb(CB_ADJ,ib)) end do end if condno=0.0d0 do i = 1,nfit if (w(i).gt.0.0d0) then condno=w(nfit)/w(i) goto 60 end if end do 60 continue write(*,'(/,A,I3,A,F6.4,A,D9.3,A,D9.3)') $ ' indep.par:',nsvfit, $ ' stepw:',stepw, $ ' lambda:',ctlpar(C_LMBDA), $ ' cond.no:',condno if (((myand(ctlint(C_PRI),AP_SV).ne.0).and.(xcy.ge.1)) $ .and.(fstat.gt.-1)) then write(*,'(/,A,A,/)') ' Eigenvalues and ' $ ,'Eigenvector Matrix of SVD-FIT ' do i = 1, nfit write(*,'(1D12.6,2X,30F6.3)') w(i), $ (evec(j,i),j=1,nfit) end do end if if (sigsq.lt.sigsqold) sigsqold=sigsq C -- stop here if convergence if (fstat.lt.0) goto 20 end do 30 continue write(*,'(/,A)') ' No Convergence, fit aborted ' 20 continue if ((iconv.eq.0).or.(ctlint(C_INTS).ne.myints) $ .or.(sigsq.gt.sigsqold)) then ctlint(C_INTS)=myints ifunpr=0 write(*,'(/,A)') ' Recalculation of the spectrum' call calc2(a,ifit,dfit,npar,nfit,palc,pali,0) end if izyk=izyk-1 write(*,'(A,/,A,I3,/)') $ '##########################################################', $ ' End at Cycle ',izyk 40 continue if (ifunpr.eq.0) then call funpr(ctlpar(C_ROFIT),ndat,chi2,wght) end if if (nfit.gt.0) then write(*,'(/,A)')' RMS deviations (MHz), B and V sorted' do j=1, 2 if (j.eq.1) write(*,'(2A3,A4,A)') 'B','V','n',' splittings MHz' if (j.eq.2) write(*,'(2A3,A4,A)') 'B','V','n',' abs. freq. MHz' do ib=1, DIMVB do iv=1, DIMVV if (ndat(ib,iv,j).ne.0) then write(*,'(2I3,I4,2F18.6)') ib,iv,ndat(ib,iv,j) $ ,sqrt(chi2(ib,iv,j)/wght(ib,iv,j))*1.D3 $ ,sqrt(chi2(ib,iv,j)/wght(ib,iv,j))*1.D3 $ *(ndat(ib,iv,j)+nfit)/ndat(ib,iv,j) end if end do end do end do end if if (ctlint(C_NZYK).eq.1) goto 100 if (nsvfit.ne.nfit) then write(*,'(/,A,I3,A,I3,A)')' **** Warning **** : only',nsvfit $ ,' of',nfit,' Parameters are linear independent!' if (ctlint(C_SVDER).ne.0) write(*,'(A)') $ ' Errors may be wrong! set ''svderr'' to 0 recommended' end if write(*,'(/,12X,A)') 'Parameters and Errors' call pra(a,da,ifit,min(myand(ctlint(C_PRI),AP_PL)+1,5),1,0) write(*,'(/,A,F15.6,A,/)') ' Standard Deviation' $ ,dsqrt(sigsqold)*1000.0d0,' MHz' do ib=1, size(S_NB) write(*,'(A,I3)')'------------------------------------- B =',ib call adjusta(a(1,ib),npar,ctlnb(CB_ADJ,ib)) call prrrp(a(1,ib),da(1,ib),covar,palc,pali,nfit,ib) if (ctlint(C_NTOP).gt.0) then call prirp(a(1,ib),da(1,ib),indeg) call prpot(a(1,ib),da(1,ib),inkj,inkc,incm) if (ctlint(C_NTOP).gt.1) then fold=a(P1_F,ib) call adjusta(a(1,ib),npar,myor(ctlnb(CB_ADJ,ib),1)) if (a(P1_F,ib).ne.fold) write(*,'(A,/,20X,A)') $ '**** Warning **** : F values not consistent ' $ ,'Use adjf = 1 to keep F values right.' do itop=1, ctlint(C_NTOP) ift=DIMPIR*(itop-1) write(*,'(A,1F18.9)') 'F(calc) ',a(P1_F+ift,ib) end do if (ctlint(C_NTOP).gt.2) write(*,'(A,1F18.9)') $ 'F12(calc)',a(P_FF,ib) end if end if end do write(*,'(/,A)') ' Errors of fitted linear combinations' write(*,'(5F15.9)') $ (anew(mod(i-1,DIMPAR)+1,int((i-1)/DIMPAR)+1),i=1, nfit) write(*,'(/,A)') $ ' Correlation Matrix of fitted linear combinations ' cmax=0.0 do i=1, nfit do j=1, i-1 if (abs(covar(j,i)).gt.cmax) then ic=i jc=j cmax=abs(covar(j,i)) end if end do write(*,'(1X,A6,$)') parstr(pali(i,1,1)) write(*,'(20F7.3)') (covar(j,i), j=1,i-1),1.0d0 end do write(*,'(A,I3,A,I3,A,F7.4,A)') ' strongest correlation between' $ ,ic,' and',jc,' (',covar(jc,ic),')' write(*,'(/,A)') ' Freedom Cofreedom Matrix of linear comb.' cmax=2.0 do i=1, nfit do j=1, i-1 if (covar(i,j).lt.cmax) then ic=i jc=j cmax=covar(i,j) end if end do write(*,'(1X,A6,$)') parstr(pali(i,1,1)) write(*,'(20F7.3)') (covar(i,j), j=1,i) end do write(*,'(A,I3,A,I3,A,F7.4,A)') ' minimum cofreedom between' $ ,ic,' and',jc,' (',cmax,')' write(*,'(/,A,A,/)') ' Eigenvalues and ' $ ,'Eigenvector Matrix of SVD-FIT ' do i = 1, nfit write(*,'(1D12.6,2X,30F6.3)') w(i), $ (evec(j,i),j=1,nfit) end do 100 continue write(*,*) if (myand(ctlint(C_PRI),AP_LT).ne.0) call ltxpr(ctlpar(C_ROFIT)) if (myand(ctlint(C_INTS),16).eq.0) then do i=1, ncalc call binnam(todo(i,Q_J),todo(i,Q_S),todo(i,Q_F),todo(i,Q_B) $ ,binfname) c write(*,*) binfname open(99,file=binfname) close(99,status='delete') end do end if stop end C---------------------------------------------------------------------- subroutine sig_func(sig_no) C called with signal(2)=SIGINT=control C C see procedure mysignal() in iamsys.f implicit none integer sig_stat,sig_no common/sig_com/sig_stat sig_stat=1 write(0,*)'Signal no.',sig_no write(0,*)'Control-C: premature termiation of xiam, finishing...' write(*,*)'Control-C: premature termiation of xiam, finishing...' return end C---------------------------------------------------------------------- subroutine funcs(ix,df,dfda,a,sig,nfit,ifit,dfit,idfrq) C interface subroutine between LM Fit and the dnv matrix implicit none include 'iam_.for' integer ix, nfit, idfrq integer ifit(DIMPAR,DIMVB),dfit(DIMFIT) real*8 df, sig, dfda(DIMFIT),a(DIMPAR,DIMVB) C local.. real*8 dfdao(DIMFIT) real*8 dfobs,dfcal,dupda(DIMFIT),dloda(DIMFIT),dfsig,dfor,dfcr integer avg,i,ref logical firstt save firstt external myand,myor integer myand,myor data firstt/.true./ dfor=0.0d0 dfcr=0.0d0 do i=1, nfit dfdao(i)=0.0d0 dfda(i)=0.0d0 end do dfobs=dln(ix,LN_FREQ) dfsig=dln(ix,LN_ERR) C average of the data, normaly avg=ix: this has no effect avg=qlin(ix,Q_AVG,Q_UP) if (avg.eq.0) avg = ix dfcal=((dnv(avg,NV_ENG,Q_UP)-dnv(avg,NV_ENG,Q_LO)) $ +(dnv(ix,NV_ENG,Q_UP)-dnv(ix,NV_ENG,Q_LO)))/2.0d0 do i=1, nfit dupda(i)=(dnv(avg,i+NV_DEF,Q_UP)+dnv(ix,i+NV_DEF,Q_UP))/2.0d0 dloda(i)=(dnv(avg,i+NV_DEF,Q_LO)+dnv(ix,i+NV_DEF,Q_LO))/2.0d0 end do C check if a difference is fitted ref=qlin(ix,Q_REF,Q_UP) if ((ref.ne.0).and.(ref.ne.ix)) then dfor=dln(ref,LN_FREQ) dfcr=dnv(ref,NV_ENG,Q_UP)-dnv(ref,NV_ENG,Q_LO) do i=1, nfit dfdao(i)=-(dnv(ref,i+NV_DEF,Q_UP)-dnv(ref,i+NV_DEF,Q_LO)) end do end if sig=dfsig if ((myand(qlin(ix,Q_STAT,Q_LO),1).ne.1).or. $ (myand(qlin(ix,Q_STAT,Q_UP),1).ne.1)) sig=NOFIT df=dfobs-dfor-dfcal+dfcr if ((ref.ne.0).and.(ref.ne.ix).and.(sig.lt.1d3*ctlpar(C_DEFER))) $ then if (abs(dfor-dfobs-dfcal+dfcr).lt.abs(df)) then if (firstt) write(*,'(/,A,A,/)') ' *Warning*:' $ ,' possible assignment error in line(s) marked with ''x''' if (firstt) firstt=.false. qlin(ix,Q_STAT,Q_UP)=myor(qlin(ix,Q_STAT,Q_UP),8) end if end if do i=1, nfit dfda(i)=dfdao(i)+dupda(i)-dloda(i) end do return end C---------------------------------------------------------------------- subroutine calc1(a,ifit,dfit,npar,nfit,palc,pali,istat) C simple calculation of the spectrum if istat .le. 0 C calculation of derivatives if istat.gt.0 C if istat.gt.1 better derivatives are used implicit none include 'iam_.for' real*8 a(DIMPAR,DIMVB) integer ifit(DIMPAR,DIMVB),dfit(DIMFIT),npar,nfit,istat real*8 palc(DIMFIT,-1:DIMPLC) integer pali(DIMFIT, 0:DIMPLC,2) C local .. integer ifitmp(DIMPAR,DIMVB) real*8 ad(DIMPAR,DIMVB) real*8 dnvtmp(DIMLIN,Q_UP:Q_LO),dnvsav(DIMLIN,Q_UP:Q_LO) integer i,j,pri,ib real*8 diffup,difflo,delta,dsum,devar(DIMPAR) integer myand real*8 myrand external myand,myrand save devar, dnvsav data devar /DIMPAR*1.0/ if (istat.le.0) then C clear the old eigenvalues und derivativs do i=1, ctlint(C_NDATA) do j=1,DIMDNV dnv(i,j,Q_UP)=0.0d0 dnv(i,j,Q_LO)=0.0d0 end do end do xde=0 if (xcy.ge.1) xde=1 ! For output control C calculate with original values, obtain derivatives if dfit(i).gt.0 call calc2(a,ifit,dfit,npar,nfit,palc,pali,0) C save the original frequencies do i=1, ctlint(C_NDATA) dnvsav(i,Q_UP)=dnv(i,NV_ENG,Q_UP) dnvsav(i,Q_LO)=dnv(i,NV_ENG,Q_LO) end do else C calculate derivatives with difference quotient if dfit(i).lt.0 C but do not calc. analtic deriv. here do ib=1, DIMVB do i=1, DIMPAR ifitmp(i,ib)=0 end do end do xde=0 if ((myand(ctlint(C_PRI),XP_DE).ne.0).and.(xcy.ge.1)) xde=1 pri=ctlint(C_PRI) ctlint(C_PRI)=myand(pri,not(AP_EH)) ctlint(C_PRI)=myand(pri,not(AP_MH)) C calculate the ad = parameters + delta do i=1, nfit if (dfit(i).lt.0) then devar(i)=-devar(i) dsum=0.0d0 do j=1, pali(i,0,1) dsum=dsum+dabs(a(pali(i,j,1),pali(i,j,2))*palc(i,j)) end do delta=(palc(i,0)/100.0d0)*dsum*devar(i)/dble(pali(i,0,1)) do ib=1, size(S_NB) do j=1, npar ad(j,ib)=a(j,ib) end do end do do j=1, pali(i,0,1) ad(pali(i,j,1),pali(i,j,2)) $ =ad(pali(i,j,1),pali(i,j,2))+delta*palc(i,j) end do call calc2(ad,ifitmp,dfit,npar,nfit,palc,pali,i) C do an opposite variation for precise derivatives if ((dfit(i).le.-2).or.(istat.gt.1)) then do ib=1, size(S_NB) do j=1, npar ad(j,ib)=a(j,ib) end do end do do j=1, pali(i,0,1) ad(pali(i,j,1),pali(i,j,2)) $ =ad(pali(i,j,1),pali(i,j,2))-delta*palc(i,j) end do do j=1, ctlint(C_NDATA) dnvtmp(j,Q_UP)=dnv(j,NV_ENG,Q_UP) dnvtmp(j,Q_LO)=dnv(j,NV_ENG,Q_LO) end do call calc2(ad,ifitmp,dfit,npar,nfit,palc,pali,i) C calculate differential quotient do j=1,ctlint(C_NDATA) diffup=dnvtmp(j,Q_UP)-dnv(j,NV_ENG,Q_UP) difflo=dnvtmp(j,Q_LO)-dnv(j,NV_ENG,Q_LO) dnv(j,i+NV_DEF,Q_UP)=diffup/(2.0d0*delta) dnv(j,i+NV_DEF,Q_LO)=difflo/(2.0d0*delta) end do else C calculate differential quotient do j=1,ctlint(C_NDATA) diffup=dnv(j,NV_ENG,Q_UP)-dnvsav(j,Q_UP) difflo=dnv(j,NV_ENG,Q_LO)-dnvsav(j,Q_LO) dnv(j,i+NV_DEF,Q_UP)=diffup/delta dnv(j,i+NV_DEF,Q_LO)=difflo/delta end do end if end if end do ctlint(C_PRI)=pri C restore the original frequencies do i=1, ctlint(C_NDATA) dnv(i,NV_ENG,Q_UP)=dnvsav(i,Q_UP) dnv(i,NV_ENG,Q_LO)=dnvsav(i,Q_LO) end do end if return end C---------------------------------------------------------------------- subroutine calc2(a,ifit,dfit,npar,nfit,palc,pali,fistat) C calculation of the eigenvalues C the evalues are put in the field of dnv(1..ndata,NV_ENG,Q_UP/LO) C the deviations DE/DPi in dnv(1..ndata,2-DIMPAR,Q_UP/LO(i)) implicit none include 'iam_.for' real*8 a(DIMPAR,DIMVB) integer ifit(DIMPAR,DIMVB),dfit(DIMFIT),npar,nfit,fistat real*8 palc(DIMFIT,-1:DIMPLC) integer pali(DIMFIT, 0:DIMPLC,2) C work real*8 h(DIMTOT,DIMTOT),evh(DIMTOT) real*8 evalv(DIMV,-DIMSIG:DIMSIG,-DIMJ:DIMJ,DIMTOP) real*8 ovv(DIMV,DIMV,DIMOVV,-DIMSIG:DIMSIG,-DIMJ:DIMJ,DIMTOP) real*8 rotm(-DIMJ:DIMJ,-DIMJ:DIMJ,1:2,DIMTOP) real*8 rott(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV,DIMTOP) real*8 tori(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV, $ -DIMSIG:DIMSIG,DIMTOP) integer qmv(DIMV),oldj(DIMTOP) real*8 ints integer j,gam,qf,ib,oib,ic,i,itop,maxi,mini,iv,ntop integer myand,not external myand save evalv,ovv,rotm,rott,tori,qmv,oib do i=1,ctlint(C_NDATA) qlin(i,Q_STAT,Q_UP)=myand(qlin(i,Q_STAT,Q_UP),(not(1))) qlin(i,Q_STAT,Q_LO)=myand(qlin(i,Q_STAT,Q_LO),(not(1))) qlin(i,Q_STAT,Q_UP)=myand(qlin(i,Q_STAT,Q_UP),(not(8))) qlin(i,Q_STAT,Q_LO)=myand(qlin(i,Q_STAT,Q_LO),(not(8))) qlin(i,Q_GK,Q_UP)=-1000 qlin(i,Q_GK,Q_LO)=-1000 end do do i=1, size(S_G) gamma(i,0)=0 end do oib =-1 do ic=1, ncalc j =todo(ic,Q_J) gam=todo(ic,Q_S) qf =todo(ic,Q_F) ib =todo(ic,Q_B) if ((myand(todo(ic,Q_STAT),2).eq.0).and. $ (fistat.ne.0)) goto 99 ctlint(C_WOODS)=ctlnb(CB_WDS,ib) ctlint(C_ADJF) =ctlnb(CB_ADJ,ib) size(S_VV)=1 do iv=1, DIMVV if (qvv(iv,1,ib).eq.-1) goto 11 size(S_VV)=iv end do 11 continue if (size(S_VV).gt.DIMVV) stop ' max VV > DIM VV' do itop=1, ctlint(C_NTOP) size(S_V+itop)=0 end do do itop=1, ctlint(C_NTOP) mini=99 maxi=-1 do iv=1, size(S_VV) if (qvv(iv,itop,ib).lt.mini) mini=qvv(iv,itop,ib) if (qvv(iv,itop,ib).gt.maxi) maxi=qvv(iv,itop,ib) end do size(S_V+itop)=maxi-mini+1 size(S_MINV+itop)=mini if (size(S_V+itop).lt.0) stop 'error in calc2' end do if (ib.ne.oib) then oib=ib do itop=1,ctlint(C_NTOP) oldj(itop)=0 end do call adjusta(a(1,ib),npar,ctlint(C_ADJF)) if ((myand(ctlint(C_PRI),AP_PC).ne.0) $ .and.(xde.ne.0)) then write(*,'(2(A,I3))') 'Parameter for B=',ib $ ,' deriv.=',fistat write(*,'(A,3F18.9)') 'BJ BK B- ' $ ,a(P_BJ,ib),a(P_BK,ib),a(P_BD,ib) do itop=1, ctlint(C_NTOP) write(*,'(A,3F18.9)') 'rho gamma beta' $ ,a(P1_RHO +DIMPIR*(itop-1),ib) $ ,a(P1_GAMA+DIMPIR*(itop-1),ib) $ ,a(P1_BETA+DIMPIR*(itop-1),ib) end do end if C calc the |m> and |K> part in the rho-system call calmk(ib,h,evalv,ovv,rotm,rott,tori $ ,a(1,ib),qmv,ifit(1,ib),npar,fistat,0) C and save the torsional integrals for the intensities if ((ctlint(C_INTS).gt.1).and.(fistat.eq.0)) $ call wrtori(tori,ib) end if C set up the rotation matrix do itop=1,ctlint(C_NTOP) call rotate(rotm(-DIMJ,-DIMJ,1,itop) $ ,a(P1_BETA+DIMPIR*(itop-1),ib),j,oldj(itop)) end do C !!!!!!!!!!!!!!!!!! ntop !!!!!!!!!!!!!!!!!!! ntop=ctlint(C_NTOP) if (gam.eq.0) ctlint(C_NTOP)=0 C this is the main routine to calculate the spectrum call calvjk(j,gam,qf,ib,h,evalv,ovv,rotm,rott,tori $ ,a(1,ib),qmv,ifit(1,ib),dfit,palc,pali,npar,fistat,evh) C calc the intensities for this ib of all eigenvalues are done if (((ic.eq.ncalc).or.(ib.ne.todo(ic+1,Q_B))) $ .and.(fistat.eq.0).and.(ctlint(C_INTS).ge.1)) then do i=1,ctlint(C_NDATA) if ((qlin(i,Q_B,Q_UP).eq.ib).and.(qlin(i,Q_B,Q_LO).eq.ib)) $ then call intens(i,ints,a(P_MUX,ib),a(P_MUY,ib),a(P_MUZ,ib) $ ,tori) dln(i,LN_INT)=ints end if end do end if ctlint(C_NTOP)=ntop C !!!!!!!!!!!!!!!!!! ntop !!!!!!!!!!!!!!!!!!! 99 continue end do C test if all eigenvalues are calculated and assigned do i=1,ctlint(C_NDATA) if ((myand(qlin(i,Q_STAT,Q_UP),16).eq.16) $ .or.(myand(qlin(i,Q_STAT,Q_LO),16).eq.16) $ .or.(fistat.eq.0)) then if ((myand(qlin(i,Q_STAT,Q_UP),1).ne.1) $ .or.(myand(qlin(i,Q_STAT,Q_LO),1).ne.1)) $ write(*,'(A,I4,A,2I4)') $ ' Assign error at line',i,' reason', $ qlin(i,Q_STAT,Q_UP),qlin(i,Q_STAT,Q_LO) end if end do C call the routines for intensity calculation if ((ctlint(C_INTS).gt.1).and.(fistat.eq.0)) then if (ctlint(C_INTS).eq.2) $ call intall(a(P_MUX,ib),a(P_MUY,ib),a(P_MUZ,ib) $ ,ctlpar(C_TEMP)) if (ctlint(C_INTS).eq.3) $ call intal2(a(P_MUX,ib),a(P_MUY,ib),a(P_MUZ,ib) $ ,ctlpar(C_TEMP)) end if return end C---------------------------------------------------------------------- subroutine calmk(ib,h,evalv,ovv,rotm,rott,tori $ ,a,qmv,ifit,npar,fistat,imaxm) C calculation of the eigenvalues and matrixelements C of the internal rotation part implicit none include 'iam_.for' integer ib,imaxm integer ifit(DIMPAR),npar,fistat integer qmv(DIMV) real*8 a(DIMPAR) real*8 h(DIMTOT,DIMTOT) real*8 evalv(DIMV,-DIMSIG:DIMSIG,-DIMJ:DIMJ,DIMTOP) real*8 ovv(DIMV,DIMV,DIMOVV,-DIMSIG:DIMSIG,-DIMJ:DIMJ,DIMTOP) real*8 rotm(-DIMJ:DIMJ,-DIMJ:DIMJ,1:2,DIMTOP) real*8 rott(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV,DIMTOP) real*8 tori(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV, $ -DIMSIG:DIMSIG,DIMTOP) C work real*8 am(DIMPM) real*8 mvec(DIMM,DIMV,-DIMJ:DIMJ) integer imfit(DIMOVV),sigma,maxm integer i,itop,k,ip,ivc,ivr,isig,im integer sdone(-20:20,DIMTOP) integer myand external myand do itop=1,ctlint(C_NTOP) do isig=-20,20 sdone(isig,itop)=0 end do end do do isig=1,size(S_G) do itop=1,ctlint(C_NTOP) if (myand(ctlint(C_WOODS),4).eq.0) then size(S_FIRV+itop)=size(S_MINV+itop) size(S_MAXV+itop)=size(S_V+itop) else size(S_FIRV+itop)=1 if (size(S_MAXV+itop).eq.0) $ size(S_MAXV+itop)=size(S_V+itop) end if if (size(S_FIRV+itop).le.0) then write(*,'(A,I3)') 'FIRST V: spezify V lines for B=',ib stop 'FIRST V < 0 error' end if sigma=gamma(isig,itop) if (sdone(sigma,itop).eq.0) then sdone(sigma,itop)=1 if (imaxm.eq.0) then maxm=size(S_MAXM+itop) end if if (imaxm.lt.0) then maxm=size(S_MINV+itop)-1+size(S_V+itop) end if if (imaxm.gt.0) then maxm=imaxm end if do i=1, DIMPM am(i)=a(DIMPRR+(itop-1)*DIMPIR+i) imfit(i)=ifit(DIMPRR+(itop-1)*DIMPIR+i) end do imfit(PM_PI)=ifit(P_FF) if (a(P_FF).ne.0.0) imfit(PM_PI)=1 imfit(PM_COS)=ifit(P_VCC) if (a(P_VCC).ne.0.0) imfit(PM_COS)=1 imfit(PM_SIN)=ifit(P_VSS) if (a(P_VSS).ne.0.0) imfit(PM_SIN)=1 do k=-size(S_MAXK),size(S_MAXK) call calcm(sigma,h $ ,evalv(1,sigma,k,itop) $ ,ovv(1,1,1,sigma,k,itop) $ ,mvec(1,1,k) $ ,am $ ,qmv $ ,imfit $ ,k $ ,maxm $ ,size(S_FIRV+itop) $ ,size(S_MAXV+itop)) end do if ((myand(ctlint(C_PRI),AP_MO).ne.0).and.(xde.ge.1)) then write(*,'(A,5I3)') ' itop,sigma,k,B,fistat' $ ,itop,sigma,k,ib,fistat do ip=1,DIMOVV do ivr= 1, size(S_MAXV+itop) write(*,'(I3,A,30F11.6)') ip,' ovv', $ ((ovv(ivr,ivc,ip,sigma,k,itop) $ ,k=-size(S_MAXK),size(S_MAXK)) $ ,ivc=1, size(S_MAXV+itop)) end do end do write(*,*) 'mvec' do im=1,2*size(S_MAXM+itop)+1 write(*,'(40F10.6)') ((mvec(im,ivc,k) $ ,k=-size(S_MAXK),size(S_MAXK)) $ ,ivc=1, size(S_MAXV+itop)) end do end if call caltori(sigma,itop,mvec,tori,fistat $ ,maxm,size(S_FIRV+itop), size(S_MAXV+itop)) end if end do ! itop end do ! isig if ((myand(ctlint(C_PRI),AP_EO).ne.0).and.(xde.ge.1)) then do itop=1, ctlint(C_NTOP) write(*,'(/,A,I3)')' Eigenvalues of one Top. B=',ib do k=-size(S_MAXK),size(S_MAXK) do i=1, size(S_MAXV+itop) do isig=1, size(S_G) if (gamma(isig,itop).eq.-NaQN) goto 10 write(*,'(F15.6,$)') $ evalv(i,gamma(isig,itop),k,itop) end do 10 continue end do write(*,*) end do end do end if return end C---------------------------------------------------------------------- subroutine caltori(sigma,itop,mvec,tori,fistat,maxm,minv,sizv) C calculation of the eigenvalues of one matrix with specified j,k,m C the evalues are put in the field of dnv(1..ndata,NV_ENG,Q_UP/LO) C the deviations DE/DPi in dnv(1..ndata,2-DIMPAR,Q_UP/LO(i)) implicit none include 'iam_.for' integer sigma,itop,fistat,maxm,minv,sizv real*8 tori(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV, $ -DIMSIG:DIMSIG,DIMTOP) real*8 mvec(DIMM,DIMV,-DIMJ:DIMJ) integer kr,kc,vr,vc,im real*8 sum1,sum2, scale integer myand external myand if (myand(ctlint(C_WOODS),1).ne.0) then C calc the torsional integrals do kr=-size(S_MAXK), size(S_MAXK) do kc=-size(S_MAXK), size(S_MAXK) do vr=1, sizv do vc=1, sizv tori(kr,kc,vr,vc,sigma,itop)=0.0d0 do im=1, 2*maxm+1 tori(kr,kc,vr,vc,sigma,itop) = $ tori(kr,kc,vr,vc,sigma,itop) $ +mvec(im,vr,kr)*mvec(im,vc,kc) end do end do end do end do end do do kr=-size(S_MAXK), size(S_MAXK) do vr=1,sizv do vc=vr+1,sizv tori(kr,kr,vr,vc,sigma,itop)=0.0d0 tori(kr,kr,vc,vr,sigma,itop)=0.0d0 end do end do do vr=1,sizv tori(kr,kr,vr,vr,sigma,itop)=1.0d0 end do end do if ((myand(ctlint(C_PRI),AP_TI).ne.0).and.(xde.ge.1)) then write(*,'(/,2(A,I3,3X))') ' torsional integrals sigma=',sigma $ ,'deriv.=',fistat do vr=1,sizv do kr=-size(S_MAXK), size(S_MAXK) write(*,'(20F8.5)') $ ((tori(kr,kc,vr,vc,sigma,itop) $ ,kc=-size(S_MAXK), size(S_MAXK)) $ ,vc=1,sizv) end do end do end if end if if (myand(ctlint(C_WOODS),1).eq.0) then do kr=-size(S_MAXK), size(S_MAXK) do kc=-size(S_MAXK), size(S_MAXK) do vr=1,sizv do vc=1,sizv tori(kr,kc,vr,vc,sigma,itop)=0.0d0 end do end do end do end do do kr=-size(S_MAXK), size(S_MAXK) do kc=-size(S_MAXK), size(S_MAXK) do vr=1,sizv tori(kr,kc,vr,vr,sigma,itop)=1.0d0 end do end do end do end if if (myand(ctlint(C_WOODS),2).ne.0) then do kr=-size(S_MAXK), size(S_MAXK) do kc= -size(S_MAXK), size(S_MAXK) do vr=1,sizv sum1=0.0d0 sum2=0.0d0 do vc=1,vr-1 sum1=sum1+tori(kr,kc,vr,vc,sigma,itop)**2 end do do vc=vr,sizv sum2=sum2+tori(kr,kc,vr,vc,sigma,itop)**2 end do do vc=vr,sizv scale=(dsqrt(1.0d0-sum1))/dsqrt(sum2) tori(kr,kc,vr,vc,sigma,itop)= $ tori(kr,kc,vr,vc,sigma,itop)*scale if (vr.ne.vc) $ tori(kc,kr,vc,vr,sigma,itop)= $ tori(kr,kc,vr,vc,sigma,itop) end do end do end do end do if ((myand(ctlint(C_PRI),AP_TI).ne.0).and.(xde.ge.1)) then write(*,'(/,A)') '----- scaled torsional integrals -----' do vr=1,sizv do kr=-size(S_MAXK), size(S_MAXK) write(*,'(20F8.5)') $ ((tori(kr,kc,vr,vc,sigma,itop) $ ,kc=-size(S_MAXK), size(S_MAXK)) $ ,vc=1,sizv) end do end do end if end if return end C --------------------------------------------------------------------- subroutine maxof(mat, dimrow, dimcol, nrow, fcol, ncol, $ besti, scndi, bestv, scndv) C find the biggest and second abs values in column=fcol..fcol+ncol-1 implicit none integer dimrow,dimcol,nrow,fcol,ncol integer besti(ncol), scndi(ncol) real*8 mat(dimrow,dimcol), bestv(ncol), scndv(ncol) C work integer ic,ir,i if ((ncol+fcol-1).gt.dimcol) stop ' ERROR: DIMCOL exceeded' if (nrow.gt.dimrow) stop ' ERROR: DIMROW exceeded' do ic=1, ncol besti(ic)=1 scndi(ic)=0 bestv(ic)=0.0d0 scndv(ic)=0.0d0 end do do ir=1, nrow do i=1, ncol ic=i+fcol-1 if (scndv(i).le.dabs(mat(ir,ic))) then if (bestv(i).le.dabs(mat(ir,ic))) then scndi(i)=besti(i) scndv(i)=bestv(i) besti(i)=ir bestv(i)=dabs(mat(ir,ic)) else scndi(i)=ir scndv(i)=dabs(mat(ir,ic)) end if end if end do end do do i=1, ncol ic=i+fcol-1 bestv(i)=mat(besti(i),ic) scndv(i)=mat(scndi(i),ic) end do return end C---------------------------------------------------------------------- subroutine savepsi(i,psi) implicit none real*8 psi integer i include 'iam_.for' dln(i,LN_PSI)=psi return end C-------------------------------------------------------------- subroutine swptdo(i,j) implicit none include 'iam_.for' integer i,j,itmp itmp=todo(i,Q_J) todo(i,Q_J)=todo(j,Q_J) todo(j,Q_J)=itmp itmp=todo(i,Q_S) todo(i,Q_S)=todo(j,Q_S) todo(j,Q_S)=itmp itmp=todo(i,Q_B) todo(i,Q_B)=todo(j,Q_B) todo(j,Q_B)=itmp itmp=todo(i,Q_F) todo(i,Q_F)=todo(j,Q_F) todo(j,Q_F)=itmp itmp=todo(i,Q_STAT) todo(i,Q_STAT)=todo(j,Q_STAT) todo(j,Q_STAT)=itmp return end C===================================================================== block data bdata C --------------------------------------------------------------------- include 'iam_.for' include 'iamdata_.for' integer idat data gamstr(1)/'G'/ data gamstr(2)/'S'/ data gamstr(3)/'G'/ data gamstr(4)/'S'/ data gamstr(5)/'G'/ data gamstr(6)/'S'/ C data gamstr /DIMTOP*'G'/ data parstr(P_BJ ) /'BJ '/, parfit(P_BJ ) /0/ data parstr(P_BK ) /'BK '/, parfit(P_BK ) /0/ data parstr(P_BD ) /'B- '/, parfit(P_BD ) /0/ data parstr(P_DJ ) /'DJ '/, parfit(P_DJ ) /0/ data parstr(P_DJK ) /'DJK '/, parfit(P_DJK ) /0/ data parstr(P_DK ) /'DK '/, parfit(P_DK ) /0/ data parstr(P_DJD ) /'dj '/, parfit(P_DJD ) /0/ data parstr(P_DKD ) /'dk '/, parfit(P_DKD ) /0/ C data parstr(P_R6 ) /'R6 '/, parfit(P_R6 ) /0/ data parstr(P_HJ ) /'H_J '/, parfit(P_HJ ) /0/ data parstr(P_HJK ) /'HJK '/, parfit(P_HJK ) /0/ data parstr(P_HKJ ) /'HKJ '/, parfit(P_HKJ ) /0/ data parstr(P_HK ) /'H_K '/, parfit(P_HK ) /0/ data parstr(P_HJD ) /'h_j '/, parfit(P_HJD ) /0/ data parstr(P_HJKD ) /'hjk '/, parfit(P_HJKD ) /0/ data parstr(P_HKD ) /'h_k '/, parfit(P_HKD ) /0/ data parstr(P_QZ ) /'chi_z '/, parfit(P_QZ ) /0/ data parstr(P_QD ) /'chi_- '/, parfit(P_QD ) /0/ data parstr(P_QXY ) /'chi_xy '/, parfit(P_QXY ) /0/ data parstr(P_QXZ ) /'chi_xz '/, parfit(P_QXZ ) /0/ data parstr(P_QYZ ) /'chi_yz '/, parfit(P_QYZ ) /0/ data parstr(P_FF ) /'F12 '/, parfit(P_FF ) /1/ data parstr(P_VSS ) /'Vss '/, parfit(P_VSS ) /1/ data parstr(P_VCC ) /'Vcc '/, parfit(P_VCC ) /1/ data parstr(P_CP ) /'C+ '/, parfit(P_CP ) /0/ data parstr(P_CZ ) /'C_z '/, parfit(P_CZ ) /0/ data parstr(P_CD ) /'C- '/, parfit(P_CD ) /0/ data parstr(P_MUX ) /'mu_x '/, parfit(P_MUX ) /9/ data parstr(P_MUY ) /'mu_y '/, parfit(P_MUY ) /9/ data parstr(P_MUZ ) /'mu_z '/, parfit(P_MUZ ) /9/ data parstr(P_PX ) /'P_x '/, parfit(P_PX ) /0/ data parstr(P_PY ) /'P_y '/, parfit(P_PY ) /0/ data parstr(P_PZ ) /'P_z '/, parfit(P_PZ ) /0/ data parstr(P1_VN1 ) /'V1n_1 '/, parfit(P1_VN1 ) /1/ data parstr(P1_VN2 ) /'V2n_1 '/, parfit(P1_VN2 ) /1/ data parstr(P1_F ) /'F_1 '/, parfit(P1_F ) /1/ data parstr(P1_RHO ) /'rho_1 '/, parfit(P1_RHO ) /1/ data parstr(P1_BETA ) /'beta_1 '/, parfit(P1_BETA ) /1/ data parstr(P1_GAMA ) /'gamma_1'/, parfit(P1_GAMA ) /1/ data parstr(P1_DPIJ ) /'Dpi2J_1'/, parfit(P1_DPIJ ) /1/ data parstr(P1_DPIK ) /'Dpi2K_1'/, parfit(P1_DPIK ) /1/ data parstr(P1_DPID ) /'Dpi2-_1'/, parfit(P1_DPID ) /1/ c data parstr(P1_DPI4 ) /'Dpi4_1 '/, parfit(P1_DPI4 ) /1/ c data parstr(P1_DPIC ) /'Dpi2c_1'/, parfit(P1_DPIC ) /1/ data parstr(P1_DC3J ) /'Dc3J_1 '/, parfit(P1_DC3J ) /1/ data parstr(P1_F0 ) /'F0_1 '/, parfit(P1_F0 ) /1/ data parstr(P1_ANGX ) /'epsil_1'/, parfit(P1_ANGX ) /1/ data parstr(P1_ANGZ ) /'delta_1'/, parfit(P1_ANGZ ) /1/ data parstr(P2_VN1 ) /'V1n_2 '/, parfit(P2_VN1 ) /1/ data parstr(P2_VN2 ) /'V2n_2 '/, parfit(P2_VN2 ) /1/ data parstr(P2_F ) /'F_2 '/, parfit(P2_F ) /1/ data parstr(P2_RHO ) /'rho_2 '/, parfit(P2_RHO ) /1/ data parstr(P2_BETA ) /'beta_2 '/, parfit(P2_BETA ) /1/ data parstr(P2_GAMA ) /'gamma_2'/, parfit(P2_GAMA ) /1/ data parstr(P2_DPIJ ) /'Dpi2J_2'/, parfit(P2_DPIJ ) /1/ data parstr(P2_DPIK ) /'Dpi2K_2'/, parfit(P2_DPIK ) /1/ data parstr(P2_DPID ) /'Dpi2-_2'/, parfit(P2_DPID ) /1/ c data parstr(P2_DPI4 ) /'Dpi4_2 '/, parfit(P2_DPI4 ) /1/ c data parstr(P2_DPIC ) /'Dpi2c_2'/, parfit(P2_DPIC ) /1/ data parstr(P2_DC3J ) /'Dc3J_2 '/, parfit(P2_DC3J ) /1/ data parstr(P2_F0 ) /'F0_2 '/, parfit(P2_F0 ) /1/ data parstr(P2_ANGX ) /'epsil_2'/, parfit(P2_ANGX ) /1/ data parstr(P2_ANGZ ) /'delta_2'/, parfit(P2_ANGZ ) /1/ data parstr(P3_VN1 ) /'V1n_3 '/, parfit(P3_VN1 ) /1/ data parstr(P3_VN2 ) /'V2n_3 '/, parfit(P3_VN2 ) /1/ data parstr(P3_F ) /'F_3 '/, parfit(P3_F ) /1/ data parstr(P3_RHO ) /'rho_3 '/, parfit(P3_RHO ) /1/ data parstr(P3_BETA ) /'beta_3 '/, parfit(P3_BETA ) /1/ data parstr(P3_GAMA ) /'gamma_3'/, parfit(P3_GAMA ) /1/ data parstr(P3_DPIJ ) /'Dpi2J_3'/, parfit(P3_DPIJ ) /1/ data parstr(P3_DPIK ) /'Dpi2K_3'/, parfit(P3_DPIK ) /1/ data parstr(P3_DPID ) /'Dpi2-_3'/, parfit(P3_DPID ) /1/ c data parstr(P3_DPI4 ) /'Dpi4_3 '/, parfit(P3_DPI4 ) /1/ c data parstr(P3_DPIC ) /'Dpi2c_3'/, parfit(P3_DPIC ) /1/ data parstr(P3_DC3J ) /'Dc3J_3 '/, parfit(P3_DC3J ) /1/ data parstr(P3_F0 ) /'F0_3 '/, parfit(P3_F0 ) /1/ data parstr(P3_ANGX ) /'epsil_3'/, parfit(P3_ANGX ) /1/ data parstr(P3_ANGZ ) /'delta_3'/, parfit(P3_ANGZ ) /1/ data ctlstr(C_NZYK ) /'nzyk '/ data ctlstr(C_NCYCL) /'ncycl '/ data ctlstr(C_PRI) /'aprint'/ data ctlstr(C_PRINT) /'print '/ data ctlstr(C_XPR) /'xprint'/ data ctlstr(C_INTS ) /'ints '/ data ctlstr(C_ORGER) /'orger '/ data ctlstr(C_EVAL ) /'eval '/ data ctlstr(C_DFRQ ) /'dfreq '/ data ctlstr(C_MAXM ) /'maxm '/ c data ctlstr(C_MAXM1) /'maxm1 '/ c data ctlstr(C_MAXM2) /'maxm2 '/ c data ctlstr(C_MAXM3) /'maxm3 '/ data ctlstr(C_MAXV) /'maxvm '/ c data ctlstr(C_MAXV1) /'maxvm1'/ c data ctlstr(C_MAXV2) /'maxvm2'/ c data ctlstr(C_MAXV3) /'maxvm3'/ data ctlstr(C_WOODS) /'woods '/ c data ctlstr(C_WOOD1) /'woods1'/ c data ctlstr(C_WOOD2) /'woods2'/ c data ctlstr(C_WOOD3) /'woods3'/ c data ctlstr(C_WOOD4) /'woods4'/ data ctlstr(C_NDATA) /'ndata '/ data ctlstr(C_NFOLD) /'nfold '/ data ctlstr(C_SPIN ) /'spin '/ data ctlstr(C_NTOP ) /'ntop '/ data ctlstr(C_ADJF ) /'adjf '/ c data ctlstr(C_ADJ1 ) /'adjf1 '/ c data ctlstr(C_ADJ2 ) /'adjf2 '/ c data ctlstr(C_ADJ3 ) /'adjf3 '/ c data ctlstr(C_ADJ4 ) /'adjf4 '/ data ctlstr(C_ROFIT) /'rofit '/ data ctlstr(C_DEFER) /'defer '/ data ctlstr(C_EPS ) /'eps '/ data ctlstr(C_WEIGF) /'weigf '/ data ctlstr(C_CNVG ) /'convg '/ data ctlstr(C_LMBDA) /'lambda'/ data ctlstr(C_FITSC) /'fitscl'/ data ctlstr(C_FRQLO) /'freq_l'/ data ctlstr(C_FRQUP) /'freq_h'/ data ctlstr(C_INTLM) /'limit '/ data ctlstr(C_SVDER) /'svderr'/ data ctlstr(C_TEMP ) /'temp '/ data ctlstr(C_RED ) /'reduct'/ data qnostr( 1 ) /'Jup '/ data qnostr( 2 ) /'K-up '/ data qnostr( 3 ) /'K+up '/ data qnostr( 4 ) /'Jlo '/ data qnostr( 5 ) /'K-lo '/ data qnostr( 6 ) /'K+lo '/ data qnostr( 7 ) /'= '/ data qnostr( 8 ) /'Err '/ data qnostr( 9 ) /'Sup '/ data qnostr(10 ) /'Slo '/ data qnostr(11 ) /'V1up '/ data qnostr(12 ) /'V1lo '/ data qnostr(13 ) /'V2up '/ data qnostr(14 ) /'V2lo '/ data qnostr(15 ) /'Bup '/ data qnostr(16 ) /'Blo '/ data qnostr(17 ) /'Fup '/ data qnostr(18 ) /'Flo '/ data qnostr(19 ) /'Tup '/ data qnostr(20 ) /'Tlo '/ data qnostr(21 ) /'tup '/ data qnostr(22 ) /'tlo '/ data qnostr(23 ) /'# '/ data qnostr(24 ) /'& '/ data qnostr(25 ) /'Kup '/ data qnostr(26 ) /'Klo '/ data qnostr(27 ) /'diff '/ data qnostr(28 ) /'GHz '/ data qnostr(29 ) /'MHz '/ data qnostr(30 ) /'cm-1 '/ data (qnostr(MAXQC+idat),idat=1,DIMGAM) /DIMGAM*' '/ data (q_q(idat),idat=1,MAXQC) $ /Q_J ,0 ,0 ,Q_J ,0 ,0 $ ,0 ,0 ,Q_S, Q_S, Q_V1,Q_V1,Q_V2,Q_V2 $ ,Q_B, Q_B, Q_F, Q_F, Q_T ,Q_T ,Q_TJ,Q_TJ $ ,Q_REF,Q_AVG, Q_K, Q_K $ ,0,0,0,0/ data (q_q(MAXQC+idat),idat=1,DIMGAM) /DIMGAM*0/ data (qul(idat),idat=1,MAXQC) /Q_UP,0 ,0 ,Q_LO,0 ,0 $ ,0 ,0 ,Q_UP,Q_LO,Q_UP,Q_LO,Q_UP,Q_LO $ ,Q_UP,Q_LO,Q_UP,Q_LO,Q_UP,Q_LO,Q_UP,Q_LO $ ,0 ,0 , Q_UP,Q_LO $ ,0,0,0,0/ data (qul(MAXQC+idat),idat=1,DIMGAM) /DIMGAM*0/ c data q_q /Q_J, Q_K, Q_J, Q_K c $ ,Q_V1,Q_V2,Q_B, Q_S, Q_V1,Q_V2,Q_B, Q_S c $ ,Q_F, Q_F, Q_T, Q_T, Q_TJ, Q_TJ, Q_REF, Q_AVG c $ ,0,0,0,0,0,0,0,0,0,0,0/ c data qul /Q_UP,Q_UP,Q_LO,Q_LO c $ ,Q_UP,Q_UP,Q_UP, Q_UP,Q_LO,Q_LO,Q_LO, Q_LO c $ ,Q_UP,Q_LO,Q_UP,Q_LO,Q_UP, Q_LO, 0, 0 c $ ,0,0,0,0,0,0,0,0,0,0,0 / C data rpstr(1) /'Ir (a,b,c <-> z,x,y) '/ C data rpstr(2) /'IIr (a,b,c <-> y,z,x) '/ C data rpstr(3) /'IIIr (a,b,c <-> x,y,z) '/ end C C------------------------------------------------------------------------------ C C module IAMIO.FOR C C------------------------------------------------------------------------------ C C ---------------------------------------------------------- subroutine parinp(a,palc,pali,ifit,dfit,npar,nfit) implicit none include 'iam_.for' real*8 a(DIMPAR,DIMVB) real*8 palc(DIMFIT,-1:DIMPLC) integer pali(DIMFIT, 0:DIMPLC,2) integer ifit(DIMPAR,DIMVB), dfit(DIMFIT),npar, nfit integer i,j,iq,iv,ib,ndata,nsplit,ift,apri,xpri,sizeb,ffree real*8 indeg,inkj,pi real*8 msplit,efdata,m2split C local variables: header character*10 helpstr C local variables: reading list of control variables real*8 dctl(DIMCPAR+DIMCINT) integer cdone(DIMCPAR+DIMCINT) integer adj C local variables: reading list of parameters integer padone(DIMPAR,DIMVB) C local variables: reading definition of gammas real*8 dgam(2*DIMTOP) integer gdone(2*DIMTOP),is,it,itop character*10 symstr C local variables: reading linear comb. of parameters to fit real*8 afit(DIMPAR,DIMVB) real*8 sum integer ix,gl,gm,oldix,df,no character*10 fitstr,spaz C local variables: reading list of transitions real*8 dqno(DIMQC) integer qdone(DIMQC),ilr,ila,ilx,iqq,il character*40 fmtstr logical stpflg integer getd,geti,getc,getbuf,myand,myor,s_mark,len_c logical getend external geti,getd,getbuf,getc,myand,myor,s_mark,len_c external getend C include 'iamdata_.for' pi=dacos(-1.0d0) inkj=3.9903132D-04 indeg=180.0d0/pi call fillsp(spaz) do i=1, DIMFIT do j=1, DIMPLC palc(i,j)=0.0d0 pali(i,j,1)=0 pali(i,j,2)=0 end do dfit(i)=0 end do do i=1, DIMPAR do j=1, DIMVB ifit(i,j)=0 end do end do C Header of Input File do while (.true.) if (getbuf(gu,ui).le.0) goto 4 call writebuf(6) call fillsp(helpstr) i=getc(gu,helpstr) if ((helpstr.eq."help").or. $ (helpstr.eq."Help").or. $ (helpstr.eq."HELP")) then write(*,'(A)') ' possible parameters: ' do i=1,DIMPAR write(*,'(X,A)') parstr(i) end do write(*,'(/,A,I3)') ' DIMTOP ',DIMTOP write(*,'(A,I3)') ' DIMJ ',DIMJ write(*,'(/,A)') $ ' See help file "xiam_v25.txt" for more information' stop end if end do 4 continue do i=1, DIMCINT+DIMCPAR cdone(i)=0 dctl(i)=0.0d0 end do do while (.true.) if (getbuf(gu,ui).le.0) goto 6 call getln(gu,ctlstr,dctl,cdone,DIMCINT+DIMCPAR) end do 6 continue do i=1, DIMCINT ctlint(i)=0 end do do i=DIMCINT+1, DIMCINT+DIMCPAR ctlpar(i)=0.0d0 end do ctlint(C_MAXM)=8 ctlint(C_NZYK )=1 ctlint(C_PRI)=0 ctlint(C_XPR)=0 ctlint(C_PRINT)=3 ctlint(C_NFOLD)=3 ctlint(C_FITSC)=0 ctlint(C_WOODS)=33 ctlint(C_ADJF)=0 ctlpar(C_EPS)=1.0d-12 ctlpar(C_DEFER)=1.0d-5 ctlpar(C_CNVG)=0.999d0 ctlpar(C_LMBDA)=0.00001d0 ctlpar(C_FRQLO)=6.0d0 ctlpar(C_FRQUP)=40.0d0 ctlpar(C_INTLM)=0.1d0 ctlpar(C_TEMP)=273.0d0 do i=1, DIMCINT if (cdone(i).ne.0) then ctlint(i)=int(dctl(i)) end if end do do i= DIMCINT+1, DIMCINT+DIMCPAR if (cdone(i).ne.0) then ctlpar(i)=dctl(i) end if end do apri=0 xpri=0 if (ctlint(C_PRINT).eq.1) apri=1 if (ctlint(C_PRINT).eq.2) apri=1 + AP_TL if (ctlint(C_PRINT).ge.2) xpri=XP_CC if (ctlint(C_PRINT).eq.3) apri=1 + AP_TF if (ctlint(C_PRINT).eq.4) apri=2 + AP_TF if (ctlint(C_PRINT).eq.5) apri=2 + AP_TF + AP_TL if (ctlint(C_PRINT).eq.6) apri=2 + AP_TL + AP_LT ctlint(C_PRI)=myor(ctlint(C_PRI),apri) ctlint(C_XPR)=myor(ctlint(C_XPR),xpri) ctlint(C_NZYK)=max(ctlint(C_NZYK),ctlint(C_NCYCL)) do i=1, DIMCINT if (mod(i,4).eq.1) write(*,*) if (ctlstr(i)(1:1).ne.'_') $ write(*,'(3X,A,3X,I5,$)') ctlstr(i),ctlint(i) end do C write(*,*) j=0 do i= DIMCINT+1, DIMCINT+DIMCPAR j=j+1 if (mod(j,3).eq.1) write(*,*) if (len(ctlstr(i)).ne.0) write(*,'(3X,A,3X,D12.7,$)') $ ctlstr(i),ctlpar(i) end do write(*,*) npar=DIMPRR+ctlint(C_NTOP)*DIMPIR if (ctlint(C_NTOP).eq.3) then if (myand(ctlint(C_ADJF),4).ne.0) write(*,*) $ ' I recommend to set adj to 4 for a three top molecule!' end if if (ctlint(C_NTOP).gt.DIMTOP) stop 'ERROR: ntop > DIMTOP' do itop=1, DIMTOP if ((2*ctlint(C_MAXM)+1).gt.DIMM) stop 'ERROR: 2maxm+1 > DIMM' end do if (ctlint(C_EVAL).ne.0) then write(*,'(A)') '\\ writing eigenvalues in file eval.out ' open(20,file='eval.out',status='unknown') end if if (ctlint(C_DFRQ).ne.0) then write(*,'(A,A)') '\\ writing deviation of frequencies in', $ ' file dfreq.out' open(21,file='dfreq.out',status='unknown') end if C set the parameter name acoording to the reduction A or S if (ctlint(C_RED).eq.1) then c parstr(P_DKD)='R6 ' write(*,*) 'Using Watson S Reduction ' end if if (ctlint(C_RED).eq.2) then parstr(P_DKD)='R6 ' write(*,*) 'Using van Eijck-Typke Reduction ' end if if (ctlint(C_RED).eq.0) then write(*,*) 'Using Watson A Reduction ' end if do itop=1, ctlint(C_NTOP) size(S_MAXM+itop)=ctlint(C_MAXM) size(S_MAXV+itop)=ctlint(C_MAXV) end do C -------------- write(*,*) do j=1, DIMVB do i=1, DIMPAR padone(i,j)=0 a(i,j)=0.0d0 end do end do do i=1,DIMPAR+1 if (getbuf(gu,ui).lt.0) stop 'reading a' if (getend(gu)) goto 5 call getxln(gu,parstr,a,padone,DIMPAR,DIMVB) end do 5 continue C fill parameters of not used tops with zeros do ib=1, DIMVB do i=DIMPRR+DIMPIR*ctlint(C_NTOP)+1,DIMPAR a(i,ib)=0.0d0 ifit(i,ib)=0 end do if (ctlint(C_NTOP).le.1) then a(P_FF,ib)=0.0d0 a(P_VSS,ib)=0.0d0 a(P_VCC,ib)=0.0d0 end if end do C get a preliminary value of sizeb sizeb=1 do i=1, DIMPRR+DIMPIR*ctlint(C_NTOP) do ib= DIMVB, 2, -1 if ((a(i,ib).ne.a(i,(ib-1))).and.(a(i,ib).ne.0.0d0)) $ sizeb=max(sizeb,ib) end do end do size(S_NB)=sizeb write(*,*) 'assumed sizeb',sizeb do ib=1, size(S_NB) adj=ctlint(C_ADJF) if ((ctlint(C_INTS).ne.0) .and. (a(P_MUX,ib).eq.0.0d0) $ .and. (a(P_MUY,ib).eq.0.0d0).and. (a(P_MUZ,ib).eq.0.0d0)) $ stop 'ERROR: need mu_x mu_y or mu_z for intensities' if (a(P_BJ,ib).eq.0.0d0) stop 'ERROR: BJ can not be zero' if (ctlint(C_NTOP).ge.1) then C if (a(P_FF,j).eq.a(P1_F,j)) a(P_FF,j)=0.0d0 C if (a(P1_F0,j).eq.a(P1_F,j)) a(P1_F0,j)=0.0d0 if (padone(P_FF,ib).ge.DIMTOP) a(P_FF,ib)=0.0d0 if (padone(P1_F0,ib).ge.2) a(P1_F0,ib)=0.0d0 if (padone(P2_F0,ib).ge.3) a(P2_F0,ib)=0.0d0 if (((a(P1_ANGZ,ib).ne.0.0d0).or.(a(P1_ANGX,ib).ne.0.0d0)) $ .and. $ (a(P1_BETA,ib).eq.0.0d0).and.(a(P1_GAMA,ib).eq.0.0d0) $ .and.(myand(adj,16).eq.0)) then adj=myor(adj,16) write(*,'(A)') ' \\ set (adj or 16)' end if if ((a(P1_F0,ib).ne.0.0d0).and.(a(P1_RHO,ib).eq.0.0d0) $ .and.(myand(adj,8).eq.0)) then adj=myor(adj,8) write(*,'(A)') ' \\ set (adj or 8)' end if if ((a(P1_F,ib).eq.0.0d0).and.(myand(adj,1).eq.0)) then adj=myor(adj,1) write(*,'(A)') ' \\ set (adj or 1)' end if end if if (ctlint(C_NTOP).ge.2) then if ((a(P_FF,ib).eq.0.0d0).and.(myand(adj,2).eq.0)) then adj=myor(adj,2) write(*,'(A)') ' \\ set (adj or 2)' end if end if c if ((a(P1_F0).eq.0.0d0).and.(a(P1_F).eq.0.0d0) c $ .and.(myand(ctlint(C_ADJF),1).eq.0)) then c ctlint(C_ADJF)=myor(ctlint(C_ADJF),1) c write(*,'(A)') ' \\ set (adj or 1)' c end if ctlnb(CB_ADJ,ib)=adj ctlnb(CB_WDS,ib)=ctlint(C_WOODS) call adjusta(a(1,ib),npar,adj) if (myand(adj,1).ne.0) write(*,*) $ '\\ adj 1: adjust F according to rho, beta and gamma' if (myand(adj,2).ne.0) write(*,*) $ '\\ adj 2: adjust F12 according to rho, beta and gamma' if (myand(adj,4).ne.0) write(*,*) $ '\\ adj 4: adjust F (one top case) and ignore F'' ' if (myand(adj,8).ne.0) write(*,*) $ '\\ adj 8: adjust rho according to F0 = 1/(2 I_alpha)' if (myand(adj,16).ne.0) write(*,*) $ '\\ adj 16: adjust beta and gamma according delta + epsil' write(*,'(A,I4)') ' new adj :',adj end do do ib=sizeb+1, DIMVB ctlnb(CB_ADJ,ib)=adj ctlnb(CB_WDS,ib)=ctlint(C_WOODS) end do call pra(a,a,ifit,3,0,0) write(*,*) C ------ do i=1,DIMFIT+1 do is=1, DIMPAR do ib=1, DIMVB padone(is,ib)=0 afit(is,ib)=0.0d0 end do end do if (getbuf(gu,ui).lt.0) stop ' Error reading fit parameters ' if (getend(gu)) goto 8 dfit(i)=0 df=0 gl=getc(gu,fitstr) if (fitstr.eq.'dqx') df=-2 if (fitstr.eq.'dqu') df=-1 if (fitstr.eq.'fit') df=1 if ((i.gt.DIMFIT).and.(df.ne.0)) $ stop ' maximum number (DIMFIT) of fit variables exceeded !' dfit(i)=df if (dfit(i).eq.0) goto 8 C palc(i,0): stepwidth in % for differential quotient calc no=getd(gu,palc(i,0)) if (no.le.0) palc(i,0)=0.1d0 C palc(i,-1): factor to scale the new parameters in fit no=getd(gu,palc(i,-1)) if (no.le.0) palc(i,-1)=1.0d0 if ((dfit(i).gt.0)) $ write(*,'(A,2D8.1,3X,$)') ' fit ',palc(i,0),palc(i,-1) if ((dfit(i).eq.-1)) $ write(*,'(A,2D8.1,3X,$)') ' dqu ',palc(i,0),palc(i,-1) if ((dfit(i).eq.-2)) $ write(*,'(A,2D8.1,3X,$)') ' dqx ',palc(i,0),palc(i,-1) call getxln(gu,parstr,afit,padone,DIMPAR,DIMVB) sum=0.0d0 oldix=0 ix=0 do is=1, DIMPAR oldix=ix do ib=1, DIMVB if (padone(is,ib).eq.1) then ifit(is,ib)=1 if (afit(is,ib).eq.0.0d0) afit(is,ib)=1.0d0 ix=ix+1 if (ix.gt.DIMPLC) stop ' Dimension Error: DIMPLC' pali(i,ix,1)=is pali(i,ix,2)=ib palc(i,ix)=afit(is,ib) sum=sum+dabs(a(is,ib)) end if end do if ((ix-oldix).eq.DIMVB) then write(*,'(3X,1A,3X,1F6.2,$)') $ parstr(is),palc(i,ix) else do j=oldix+1, ix write(*,'(3X,2A,I1,2A,1F5.2,$)') $ parstr(is)(1:len_c(parstr(is))) $ ,'(',pali(i,j,2),')' $ ,spaz(1:(len(parstr(is))-len_c(parstr(is))+1)) $ ,palc(i,j) end do end if end do write(*,*) pali(i,0,1)=ix do j=1, pali(i,0,1) if ((parfit(pali(i,j,1)).ne.0).and.(dfit(i).eq.1)) then c write(*,'(3A)')' Warning: fit changed to dqu!', c $ ' No analytic derivatives for ',parstr(j) dfit(i)=-1 end if end do if ((dfit(i).eq.-1).and.(sum.eq.0.0d0)) $ stop 'ERROR: dqu/x parameter can not be zero !' end do 8 continue nfit=i-1 write(*,*) do ib=1, DIMVB adj=ctlnb(CB_ADJ,ib) if (myand(adj,2).gt.0) then if (ifit(P_FF,ib).ne.0) $ stop ' Fit ERROR: can not fit F12: adjf = 2' end if do itop=1, ctlint(C_NTOP) ift=(itop-1)*DIMPIR if (myand(adj,1).gt.0) then if (ifit(P1_F +ift,ib).ne.0) $ stop ' Fit ERROR: can not fit F: adjf = 1' else if (ifit(P1_F0+ift,ib).ne.0) $ stop ' Fit ERROR: can not fit F0: adjf <> 1' end if if (myand(adj,8).gt.0) then if (ifit(P1_RHO+ift,ib).ne.0) $ stop ' Fit ERROR: can not fit rho: adjf = 8' else if (ifit(P1_F0+ift,ib).ne.0) $ stop ' Fit ERROR: can not fit F0: adjf <> 8' end if if (myand(adj,16).gt.0) then if ((ifit(P1_GAMA+ift,ib).ne.0) $ .or.(ifit(P1_BETA+ift,ib).ne.0)) $ stop ' Fit ERROR: can not fit beta/gamma: adjf = 16' else if ((ifit(P1_ANGZ+ift,ib).ne.0) $ .or.(ifit(P1_ANGX+ift,ib).ne.0)) $ stop ' Fit ERROR: can not fit delta/epsil: adjf <> 16' end if end do ! itop end do ! ib C ------ size(S_G)=0 do i=1, DIMGAM do itop=1, DIMTOP gamma(i,itop)=NaQN end do gamma(i,0)=0 end do do i=1, 2*DIMTOP dgam(i)=0 end do if (ctlint(C_NTOP).ne.0) then do i=1, DIMGAM+1 do itop=1, 2*DIMTOP gdone(itop)=0 end do if (getbuf(gu,ui).lt.0) stop ' Error: reading gamma' if (getend(gu)) goto 7 if (i.gt.DIMGAM) stop ' to many lines reading S! ' C symmetry-species-name beginning with '/' gm=s_mark(gu) gl=getc(gu,symstr) call fillsp(qnostr(MAXQC+i)) if (symstr(1:1).eq.'/') then qnostr(MAXQC+i)=symstr(1:len(qnostr(MAXQC+i))+1) else call g_mark(gu,gm) end if call getln(gu,gamstr,dgam,gdone,2*DIMTOP) do itop=1, ctlint(C_NTOP) if ((int(dgam(2*itop-1)).eq.0).and.(int(dgam(2*itop)).eq.0)) $ gamma(i,itop)=0 if ((int(dgam(2*itop-1)).ne.0).and.(int(dgam(2*itop)).eq.0)) $ gamma(i,itop)=int(dgam(2*itop-1)) if ((int(dgam(2*itop-1)).eq.0).and.(int(dgam(2*itop)).ne.0)) $ gamma(i,itop)=int(dgam(2*itop)) if ((int(dgam(2*itop-1)).ne.0).and.(int(dgam(2*itop)).ne.0)) $ stop 'Error: use G or S as Keyword' end do size(S_G)=i end do 7 continue do is=1, size(S_G) if (qnostr(MAXQC+is)(1:1).ne.' ') $ write(*,'(3X,A,$)') qnostr(MAXQC+is) write(*,'(2X,A,4I4)') 'S ',(gamma(is,it),it=1,ctlint(C_NTOP)) do it=1,ctlint(C_NTOP) if (gamma(is,it).gt.DIMSIG) stop 'ERROR: sigma > DIMSIG' end do end do else gamma(1,1)=-999 end if C ---- do i=1, DIMVV do ib=1, DIMVB do itop=1, DIMTOP qvv(i,itop,ib)=-1 end do end do end do size(S_NB)=1 if (ctlint(C_NTOP).ne.0) then write(*,*) do ib=1, DIMVB+1 if (getbuf(gu,ui).lt.0) stop ' Error: reading qvv' if (getend(gu)) goto 20 do iv=1, DIMVV call fillsp(fitstr) gl=getc(gu,fitstr) if (fitstr(1:1).ne.'V') goto 19 do itop=1, ctlint(C_NTOP) no=geti(gu,j) if (no.le.0) stop ' Error: V no. for each top necessary!' qvv(iv,itop,ib)=j+1 end do end do 19 continue size(S_NB)=ib end do 20 continue do ib=size(S_NB)+1, DIMVB do i=1, DIMVV do itop=1, DIMTOP qvv(i,itop,ib)=qvv(i,itop,size(S_NB)) end do end do end do do ib=1, size(S_NB) do iv=1, DIMVV if (qvv(iv,1,ib).eq.-1) goto 21 write(*,'(A,$)') ' V' do itop=1,ctlint(C_NTOP) write(*,'(I3,$)') qvv(iv,itop,ib)-1 end do end do 21 continue write(*,*) end do end if C ---- write(*,*) if (ctlint(C_NDATA).le.0) ctlint(C_NDATA)=DIMLIN do iq=1, DIMQC dqno(iq)=0.0d0 end do C initial values dqno( 8)=NOFIT ! Err dqno(15)=1 ! Bup dqno(16)=1 ! Blo dqno(17)=-1 ! Fup dqno(18)=-1 ! Flo dqno(19)=0 ! Tup dqno(20)=0 ! Tlo dqno(11)=1 ! V1up dqno(12)=1 ! V1lo ndata =0 msplit =0.0d0 m2split =0.0d0 nsplit =0 efdata =0.0d0 stpflg =.false. do il=1,ctlint(C_NDATA) if (getbuf(gu,ui).lt.0) goto 10 if (getend(gu)) goto 10 dqno(23)=0 ! # dqno(24)=0 ! & dqno(27)=0 ! diff do iq=1, DIMQC qdone(iq)=0 end do call getln(gu,qnostr,dqno,qdone,DIMQC) if ((qdone(25).ne.0).and.(qdone(21).ne.0)) $ write(*,*) 'UP Tau overrides K',il if ((qdone(26).ne.0).and.(qdone(22).ne.0)) $ write(*,*) 'LO Tau overrides K',il if (((qdone(2).ne.0).or.(qdone(3).ne.0)).and.(qdone(21).ne.0)) $ write(*,*) 'UP Tau (t) overrides K- K+',il if (((qdone(5).ne.0).or.(qdone(6).ne.0)).and.(qdone(22).ne.0)) $ write(*,*) 'LO Tau (t) overrides K- K+',il if ((int(dqno(1)).gt.DIMJ).or.(int(dqno(4)).gt.DIMJ)) then write(*,'(A,I4,A)') $ ' Warning Max. J exceeded. Line',il,' error' stpflg=.true. end if if (qdone(19).ne.0) then ! Tup dqno(25)=0.0d0 ! Kup dqno(21)=0.0d0 ! tup qdone(25)=0 dqno( 2)=0.0d0 ! K- up dqno( 3)=0.0d0 ! K+ up end if if (qdone(20).ne.0) then ! Tlo dqno(26)=0.0d0 ! Klo dqno(22)=0.0d0 ! tlo qdone(26)=0 dqno( 5)=0.0d0 ! K- lo dqno( 6)=0.0d0 ! K+ lo end if if (qdone(21).ne.0) then ! tup dqno(25)=0.0d0 ! Kup dqno(19)=0.0d0 ! Tup qdone(25)=0 dqno( 2)=0.0d0 ! K- up dqno( 3)=0.0d0 ! K+ up end if if (qdone(22).ne.0) then ! tlo dqno(20)=0.0d0 ! Tlo dqno(26)=0.0d0 ! Klo qdone(26)=0 dqno( 5)=0.0d0 ! K- lo dqno( 6)=0.0d0 ! K+ lo end if if (qdone(25).ne.0) then dqno( 2)=-10.0d0 ! K- up dqno( 3)=-10.0d0 ! K+ up dqno(21)=0.0d0 ! tup,Q_TJ qdone(21)=0 end if if (qdone(26).ne.0) then dqno( 5)=-20.0d0 ! K- lo dqno( 6)=-20.0d0 ! K+ lo dqno(22)=0.0d0 ! tlo,Q_TJ qdone(22)=0 end if C copy the vector dqno to the transition list qlin do iq=1, DIMQC if ((q_q(iq).gt.0).and.(qul(iq).gt.0)) then qlin(il,q_q(iq),qul(iq))=int(dqno(iq)) end if end do if (qdone(7).ne.0) then ! = (frequency) dln(il,LN_FREQ)=dqno(7) dln(il,LN_ERR)=ctlpar(C_DEFER) if (qdone(8).ne.0) then ! Err dln(il,LN_ERR)=dqno(8) end if ndata=ndata+1 else dln(il,LN_FREQ)=0.0d0 dln(il,LN_ERR)=NOFIT end if if ((qdone(2).ne.0).or.(qdone(3).ne.0)) then qlin(il,Q_TJ,Q_UP)=0 end if if ((qdone(5).ne.0).or.(qdone(6).ne.0)) then qlin(il,Q_TJ,Q_LO)=0 end if C Check for symmetry-labels do i=1,DIMGAM if (qdone(MAXQC+i).ne.0) then if ((qdone(9).eq.0).and.(qdone(10).eq.0)) then qlin(il,Q_S,Q_UP)=i qlin(il,Q_S,Q_LO)=i else write(*,*)' Warning: S and Symmetry Label given (using S)' end if end if end do if ((qlin(il,Q_S,Q_UP).gt.size(S_G)).or. $ (qlin(il,Q_S,Q_LO).gt.size(S_G))) then write(*,*)'ERROR: Symmetry spezies not defined in line',il stpflg=.true. end if C Check for units (GHz default, MHZ and cm) if ((qdone(29).ne.0).and.(qdone(30).eq.0)) then dln(il,LN_FREQ)=dln(il,LN_FREQ)/1000.0d0 if (qdone(8).ne.0) $ dln(il,LN_ERR)=dln(il,LN_ERR)/1000.0d0 end if if ((qdone(29).eq.0).and.(qdone(30).ne.0)) then dln(il,LN_FREQ)=dln(il,LN_FREQ)*29.97925d0 if (qdone(8).ne.0) $ dln(il,LN_ERR)=dln(il,LN_ERR)*29.97925d0 end if C # references (23): ref lines ilr=int(dqno(23)) ilr=ilr+il qlin(il,Q_REF,Q_UP)=ilr qlin(il,Q_REF,Q_LO)=ilr C splitting frequency as input (27) if (qdone(27).ne.0) then ! diff ilx=int(dqno(27)) if (ilx.eq.0) then ilx=ilr else ilx=ilx+il end if if ((dln(il,LN_FREQ).ne.0.0d0).and. $ (dln(ilx,LN_FREQ).ne.0.0d0)) then dln(il,LN_FREQ)=dln(ilx,LN_FREQ)+dln(il,LN_FREQ) else write(0,'(A,2I4)') 'WARNING: diff frequency zero at',il,ilx end if end if if ((ilr.ne.il).and.(dln(il,LN_ERR).ne.NOFIT)) then msplit=msplit+abs(dln(il,LN_FREQ)-dln(ilr,LN_FREQ)) m2split=m2split+(dln(il,LN_FREQ)-dln(ilr,LN_FREQ))**2 nsplit=nsplit+1 endif efdata=efdata+ctlpar(C_DEFER)**2/dln(il,LN_ERR)**2 C Check for an average line (dqno(24)) ila=int(dqno(24)) ila=ila+il qlin(il,Q_AVG,Q_UP)=ila qlin(il,Q_AVG,Q_LO)=ila if (ila.ne.il) then if (dln(ila,LN_ERR).ne.NOFIT) then write(*,'(A,I4,A,I4,A)') $ 'Warning: average freq. in',il,': line',ila,' not used' ilr=qlin(ila,Q_REF,Q_UP) if (ilr.ne.ila) then msplit=msplit-abs(dln(ila,LN_FREQ)-dln(ilr,LN_FREQ)) m2split=m2split-(dln(ila,LN_FREQ)-dln(ilr,LN_FREQ))**2 nsplit=nsplit-1 endif efdata=efdata-ctlpar(C_DEFER)**2/dln(ila,LN_ERR)**2 dln(ila,LN_ERR)=NOFIT efdata=efdata+ctlpar(C_DEFER)**2/dln(ila,LN_ERR)**2 end if end if C calc tau (Q_TJ)(dqno(21/22) from K- dqno(2/5) and K+ dqno(3/6) if (((dqno(21).eq.0).and.(dqno(19).eq.0)).and. $ ((dqno(2).ge.0.0d0).or.(dqno(3).ge.0.0d0))) $ qlin(il,Q_TJ,Q_UP)=int(dqno(2))-int(dqno(3))+1 $ +int(dqno(1)) if (((dqno(22).eq.0).and.(dqno(20).eq.0)).and. $ ((dqno(5).ge.0.0d0).or.(dqno(6).ge.0.0d0))) $ qlin(il,Q_TJ,Q_LO)=int(dqno(5))-int(dqno(6))+1 $ +int(dqno(4)) C Q_TJ numbering of eigenvalues of one V if (qlin(il,Q_TJ,Q_UP).ne.0) then qlin(il,Q_STAT,Q_UP)=myor(qlin(il,Q_STAT,Q_UP),2) end if if (qlin(il,Q_TJ,Q_LO).ne.0) then qlin(il,Q_STAT,Q_LO)=myor(qlin(il,Q_STAT,Q_LO),2) end if C Q_T numbering of eigenvalues of the whole matrix (with several V''s) if (qlin(il,Q_T,Q_UP).ne.0) then qlin(il,Q_STAT,Q_UP)=myor(qlin(il,Q_STAT,Q_UP),4) end if if (qlin(il,Q_T,Q_LO).ne.0) then qlin(il,Q_STAT,Q_LO)=myor(qlin(il,Q_STAT,Q_LO),4) end if if (dln(il,LN_ERR).ne.NOFIT) then qlin(il,Q_STAT,Q_UP)=myor(qlin(il,Q_STAT,Q_UP),16) qlin(il,Q_STAT,Q_LO)=myor(qlin(il,Q_STAT,Q_LO),16) end if write(fmtstr,'(A,I2,A,I2,A)') $ '(',DIMQLP,'I3,3X,',DIMQLP,'I3,F18.7,D12.2)' if (myand(ctlint(C_PRI),AP_IO).ne.0) write(*,fmtstr) $ (qlin(il,iqq,Q_UP),iqq=1,DIMQLP) $ ,(qlin(il,iqq,Q_LO),iqq=1,DIMQLP) $ ,dln(il,LN_FREQ),dln(il,LN_ERR) if (myand(ctlint(C_PRI),AP_IO).ne.0) then if (qlin(il,Q_S,Q_LO).eq.1) then write(0,'(2(3I3,2X))') $ qlin(il,Q_J,Q_UP) $ ,qlin(il,Q_TJ,Q_UP)/2 $ ,(2*qlin(il,Q_J,Q_UP)-qlin(il,Q_TJ,Q_UP)+2)/2 $ ,qlin(il,Q_J,Q_LO) $ ,qlin(il,Q_TJ,Q_LO)/2 $ ,(2*qlin(il,Q_J,Q_LO)-qlin(il,Q_TJ,Q_LO)+2)/2 write(0,'(I4,F15.4)') $ qlin(il,Q_S,Q_LO),dln(il,LN_FREQ)*1000.0d0 end if if (qlin(il,Q_S,Q_LO).gt.1) then write(0,'(I4,F15.4)') $ qlin(il,Q_S,Q_LO),dln(il,LN_FREQ)*1000.0d0 end if end if if ((qlin(il,Q_TJ,Q_UP).le.0).and.(qlin(il,Q_T,Q_UP).le.0) $ .and.(qdone(25).eq.0)) then write(0,'(A,i4)') 'No tau or K (up) in line',il stpflg=.true. end if if ((qlin(il,Q_TJ,Q_LO).le.0).and.(qlin(il,Q_T,Q_LO).le.0) $ .and.(qdone(26).eq.0)) then write(0,'(A,i4)') 'No tau or K (lo) in line',il stpflg=.true. end if if ((qlin(il,Q_TJ,Q_LO).gt.(2*qlin(il,Q_J,Q_LO)+1)).or. $ (qlin(il,Q_TJ,Q_UP).gt.(2*qlin(il,Q_J,Q_UP)+1))) then write(0,'(A,i4)') 'tau > 2J+1 at line',il stpflg=.true. end if if ((qlin(il,Q_V1,Q_LO).le.0).or.(qlin(il,Q_V1,Q_UP).le.0)) then write(0,'(A,i4)') 'V < 1 at line',il stpflg=.true. end if if ((qlin(il,Q_B,Q_LO).le.0).or.(qlin(il,Q_B,Q_UP).le.0)) then write(0,'(A,i4)') 'Error: B < 1 at line',il stpflg=.true. end if if ((qlin(il,Q_F,Q_LO).lt.-1).or.(qlin(il,Q_F,Q_UP).lt.-1))then write(0,'(A,i4)') 'Error: F < -1 at transition no.',il stpflg=.true. end if if (qlin(il,Q_F,Q_LO).ge.0) then if (abs(2*qlin(il,Q_J,Q_LO)-qlin(il,Q_F,Q_LO)) $ .gt.ctlint(C_SPIN)) then write(0,'(A,i4)') 'Error: |J-F| > ! at transition no.',il stpflg=.true. end if end if if (qlin(il,Q_F,Q_UP).ge.0) then if (abs(2*qlin(il,Q_J,Q_UP)-qlin(il,Q_F,Q_UP)) $ .gt.ctlint(C_SPIN)) then write(0,'(A,i4)') 'Error: |J-F| > I at transition no.',il stpflg=.true. end if end if if ((qlin(il,Q_B,Q_LO).gt.DIMVB) $ .or.(qlin(il,Q_B,Q_UP).gt.DIMVB)) then write(0,'(A,i4)') 'B > DIM B' stpflg=.true. end if if ((qlin(il,Q_B,Q_LO).gt.size(S_NB)) $ .or.(qlin(il,Q_B,Q_UP).gt.size(S_NB))) then size(S_NB)=max(qlin(il,Q_B,Q_LO),qlin(il,Q_B,Q_UP)) end if end do 10 continue if (stpflg) then write(*,*) 'INPUT ERROR(S)' write(0,*) 'INPUT ERROR(S)' stop end if ctlint(C_NDATA)=il-1 write(*,'(3(3X,A,I4),/,3X,A,F6.1)') $ ctlstr(C_NDATA),ctlint(C_NDATA) $ ,'Data Points',ndata,'Splittings',nsplit $ ,'Effective Data Points',efdata if (nsplit.ne.0) then msplit=msplit/nsplit m2split=sqrt(m2split/nsplit) write(*,'(3X,2(A,F12.6))') 'Mean Experimental Splitting:',msplit $ ,' squared:',m2split end if C clean up pali and palc for not used ib values do i=1, nfit do j=1, pali(i,0,1) if (pali(i,j,2).gt.size(S_NB)) then pali(i,j,2)=0 pali(i,j,1)=0 palc(i,j)=0.0d0 end if end do end do do i=1, nfit do j=1, pali(i,0,1) if (pali(i,j,2).eq.0) goto 100 end do 100 continue ffree=j do j=ffree, pali(i,0,1) if (pali(i,j,2).gt.0) then pali(i,ffree,1)=pali(i,j,1) pali(i,ffree,2)=pali(i,j,2) palc(i,ffree)=palc(i,j) pali(i,j,1)=0 pali(i,j,2)=0 palc(i,j)=0.0d0 ffree=ffree+1 end if end do pali(i,0,1)=ffree-1 end do C- delete here c do ib=1, size(S_NB) cc do iv=1, DIMVV c if (qvv(iv,1,ib).eq.-1) goto 22 c write(*,'(A,$)') ' V' cc do itop=1,ctlint(C_NTOP) c write(*,'(I3,$)') qvv(iv,itop,ib)-1 c end do c end do c 22 continue c write(*,*) c end do return end C---------------------------------------------------------------------- subroutine prrrp(a,da,covar,palc,pali,nfit,ib) C print rotational constants and errors implicit none include 'iam_.for' real*8 a(DIMPAR),da(DIMPAR),covar(DIMFIT,DIMFIT) real*8 palc(DIMFIT,-1:DIMPLC) integer pali(DIMFIT, 0:DIMPLC,2) integer nfit,ib integer ibj,ibk,ibd,i,j integer ixj,ixk,ixd real*8 cjk,cjd,ckd,kappa real*8 b(3),d(3),bb,dd character*3 c(3),cc logical noc C find correlation coeff between BJ BK and B- C works only if no linear comb. of BJ BK B- is fitted ibj=0 ibk=0 ibd=0 ixj=0 ixk=0 ixd=0 do i=1, nfit do j=1, pali(i,0,1) if ((pali(i,j,1).eq.P_BJ).and.(pali(i,j,2).eq.ib)) then ibj=i ixj=ixj+1 end if end do do j=1, pali(i,0,1) if ((pali(i,j,1).eq.P_BK).and.(pali(i,j,2).eq.ib)) then ibk=i ixk=ixk+1 end if end do do j=1, pali(i,0,1) if ((pali(i,j,1).eq.P_BD).and.(pali(i,j,2).eq.ib)) then ibd=i ixd=ixd+1 end if end do end do C covar(j,i) with j jup' ntop=ctlint(C_NTOP) inti=0.0 intr=0.0 dj=dble(jlo) C muz do ik=1, sizelo(S_K) do iv=1, sizelo(S_VV) i=ik+(iv-1)*sizelo(S_K) qk=qvklo(i,Q_K) dq=dble(qk) intr=intr+2.0*muz*dq*(vrlo(i)*vrup(i)+vilo(i)*viup(i)) inti=inti+2.0*muz*dq*(vrlo(i)*viup(i)-vilo(i)*vrup(i)) end do end do C mux/y do ik2=1, sizelo(S_K)-1 qk2=qvkup(ik2,Q_K) ik1=ik2+1 qk1=qvklo(ik1,Q_K) dq=sqrt((dj-dble(qk2))*(dj+dble(qk2)+1.0)) do iv1=1, sizelo(S_VV) i1=ik1+(iv1-1)*sizelo(S_K) do iv2=1, sizeup(S_VV) i2=ik2+(iv2-1)*sizeup(S_K) tt=1.0 if (gam.ne.0) then do itop=1, ntop tt=tt*tori(qk1,qk2 $ ,qvklo(i1,Q_V+itop)-size(S_MINV+itop)+1 $ ,qvkup(i2,Q_V+itop)-size(S_MINV+itop)+1 $ ,gamma(gam,itop),itop) end do else if (iv1.ne.iv2) tt=0.0 end if intr=intr $ -tt*mux*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) inti=inti $ -tt*mux*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) inti=inti $ -tt*muy*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) intr=intr $ +tt*muy*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) end do end do end do do ik2=2, sizelo(S_K) qk2=qvkup(ik2,Q_K) ik1=ik2-1 qk1=qvklo(ik1,Q_K) dq=sqrt((dj+dble(qk2))*(dj-dble(qk2)+1.0)) do iv1=1, sizelo(S_VV) i1=ik1+(iv1-1)*sizelo(S_K) do iv2=1, sizeup(S_VV) i2=ik2+(iv2-1)*sizeup(S_K) tt=1.0 if (gam.ne.0) then do itop=1, ntop tt=tt*tori(qk1,qk2 $ ,qvklo(i1,Q_V+itop)-size(S_MINV+itop)+1 $ ,qvkup(i2,Q_V+itop)-size(S_MINV+itop)+1 $ ,gamma(gam,itop),itop) end do else if (iv1.ne.iv2) tt=0.0 end if intr=intr $ -tt*mux*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) inti=inti $ -tt*mux*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) inti=inti $ +tt*muy*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) intr=intr $ -tt*muy*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) end do end do end do ints=(intr**2+inti**2)*(2.0*dj+1.0)/(2.0*dj*(2.0*dj+2.0)) return end C ---------------------------------------------------------------------- C Delta J = 1 ---------------------------------------------------------------- subroutine calir(ints,tori,jup,jlo,gam,muz, mux, muy $ ,sizeup,qvkup,vrup,viup $ ,sizelo,qvklo,vrlo,vilo) implicit none include 'iam_.for' integer jup,jlo,gam real*8 ints,muz, mux, muy real*8 tori(-DIMJ:DIMJ,-DIMJ:DIMJ,DIMV,DIMV, $ -DIMSIG:DIMSIG,DIMTOP) integer sizeup(DIMSIZ),sizelo(DIMSIZ) integer qvkup(DIMTOT,Q_K:Q_V+DIMTOP),qvklo(DIMTOT,Q_K:Q_V+DIMTOP) real*8 vrup(DIMTOT),vrlo(DIMTOT),viup(DIMTOT),vilo(DIMTOT) C real*8 intr,inti,dq,dj,tt integer ik,ik1,ik2,qk,qk1,qk2,iv,iv1,iv2,i,i1,i2 integer itop,ntop if ((jlo+1).ne.jup) pause 'CALIR: intens error: jlo+1 <> jup' ntop=ctlint(C_NTOP) inti=0.0 intr=0.0 dj=dble(jlo) C mu_z do ik2=2, sizelo(S_K)+1 qk2=qvkup(ik2,Q_K) ik1=ik2-1 qk1=qvklo(ik1,Q_K) do iv=1, sizelo(S_VV) i1=ik1+(iv-1)*sizelo(S_K) i2=ik2+(iv-1)*sizeup(S_K) c qk1=qvklo(i1,Q_K) c qk2=qvkup(i2,Q_K) if (qk1.ne.qk2) stop 'Error in CALIR : K numbers' dq=2.0*sqrt((dj+dble(qk1)+1.0)*(dj-dble(qk1)+1.0)) intr=intr-muz*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) inti=inti-muz*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) end do end do C mux/y do ik2=1, sizelo(S_K) qk2=qvkup(ik2,Q_K) ik1=ik2 qk1=qvklo(ik1,Q_K) dq=sqrt((dj-dble(qk2))*(dj-dble(qk2)+1.0)) do iv1=1, sizelo(S_VV) i1=ik1+(iv1-1)*sizelo(S_K) do iv2=1, sizeup(S_VV) i2=ik2+(iv2-1)*sizeup(S_K) tt=1.0 if (gam.ne.0) then do itop=1, ntop tt=tt*tori(qk1,qk2 $ ,qvklo(i1,Q_V+itop)-size(S_MINV+itop)+1 $ ,qvkup(i2,Q_V+itop)-size(S_MINV+itop)+1 $ ,gamma(gam,itop),itop) end do else if (iv1.ne.iv2) tt=0.0 end if intr=intr $ +tt*mux*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) inti=inti $ +tt*mux*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) inti=inti $ +tt*muy*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) intr=intr $ -tt*muy*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) end do end do end do do ik2=3, sizelo(S_K)+2 qk2=qvkup(ik2,Q_K) ik1=ik2-2 qk1=qvklo(ik1,Q_K) dq=sqrt((dj+dble(qk2))*(dj+dble(qk2)+1.0)) do iv1=1, sizelo(S_VV) i1=ik1+(iv1-1)*sizelo(S_K) do iv2=1, sizeup(S_VV) i2=ik2+(iv2-1)*sizeup(S_K) tt=1.0 if (gam.ne.0) then do itop=1, ntop tt=tt*tori(qk1,qk2 $ ,qvklo(i1,Q_V+itop)-size(S_MINV+itop)+1 $ ,qvkup(i2,Q_V+itop)-size(S_MINV+itop)+1 $ ,gamma(gam,itop),itop) end do else if (iv1.ne.iv2) tt=0.0 end if intr=intr $ -tt*mux*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) inti=inti $ -tt*mux*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) inti=inti $ +tt*muy*dq*(vrlo(i1)*vrup(i2)+vilo(i1)*viup(i2)) intr=intr $ -tt*muy*dq*(vrlo(i1)*viup(i2)-vilo(i1)*vrup(i2)) end do end do end do ints=(intr**2+inti**2)/(2.0*(2.0*dj+2.0)) return end C ---------------------------------------------------------------------- subroutine binnam(ij,gam,if,ib,binfname) implicit none include 'iam_.for' integer ij,if,gam,ib character*(*) binfname if (if.ge.0) then write(binfname,'(A,2I1,A,2I1,A,I1,A,2I1)') $ 'j',int(ij/10),mod(ij,10), $ 'f',int(if/10),mod(if,10), $ 'b',ib, $ '.s',int(gam/10),mod(gam,10) else write(binfname,'(A,2I1,A,I1,A,2I1)') $ 'j',int(ij/10),mod(ij,10), $ 'b',ib, $ '.s',int(gam/10),mod(gam,10) end if return end C ---------------------------------------------------------------------- subroutine wrvec(zr,zi,d,ij,gam,if,ib,qvk,qmk) implicit none include 'iam_.for' integer ij,if,gam,ib real*8 zr(DIMTOT,DIMTOT),zi(DIMTOT,DIMTOT) real*8 d(DIMTOT) integer qvk(DIMTOT,Q_K:Q_V+DIMTOP),qmk(DIMTOT,DIMQLP) integer ie,id,itop,bu,ios character*30 binfname integer myand external myand call binnam(ij,gam,if,ib,binfname) call getfu(bu) if (myand(ctlint(C_PRI),AP_ST).ne.0) $ write(*,*)' Writing ',binfname open(bu,file=binfname,status='unknown',err=99,iostat=ios $ ,form='unformatted') write(bu,err=99,iostat=ios) (size(ie),ie=1,DIMSIZ) write(bu,err=99,iostat=ios) (qvk(ie,Q_K), $ (qvk(ie,Q_V+itop),itop=1,DIMTOP) $ ,ie=1,size(S_H)) do id=1, size(S_H) c write(bu,'(F18.9,$)',err=99,iostat=ios) d(id) c write(bu,err=99,iostat=ios) d(id) write(bu,err=99,iostat=ios) d(id), (qmk(id,ie),ie=1,DIMQLP), $ (zr(ie,id),zi(ie,id),ie=1,size(S_H)) end do close (bu) return 99 continue write(*,*) 'wrvec Error iostat', ios,ij,gam,if,ib,binfname return end C---------------------------------------------------------------------- subroutine getfu(myunit) implicit none integer myunit,funit save funit data funit/10/ if (funit.gt.50) funit=10 funit=funit+1 myunit=funit c write(0,*)'funit ',myunit return end C---------------------------------------------------------------------- subroutine rdvec(ij,it,is,if,ib,mysize,qvk,vr,vi,e,qmk) implicit none include 'iam_.for' integer qvk(DIMTOT,Q_K:Q_V+DIMTOP) integer mysize(DIMSIZ) integer ij,it,is,if,ib,i,ie,itop,bu integer qmk(DIMQLP) real*8 vr(DIMTOT),vi(DIMTOT),e character*30 binfname call getfu(bu) c write(0,*)'rdvec',bu call binnam(ij,is,if,ib,binfname) open(bu,file=binfname,status='unknown',form='unformatted') read