C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C ASQ1P - FIRST ORDER PREDICTIVE PROGRAM FOR ASYMMETRIC TOP WITH ONE C QUADRUPOLAR NUCLEUS C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C - ASQ1P carries out first order single-quadrupole calculation for an C asymmetric top and will output frequency shifts and coefficients a and b C of X.aa and X.bb in: C Freq_hyperfine = Freq_hypefinefree + a X_aa + b X_bb C The coefficients can be used in an external least squares program for C fitting X_aa and X_bb from just the splittings, or for deriving C hyperfine free frequencies C C - Frequency shifts are printed out in groups according to the F selection C rule and are also sorted in ascending order. Relative intensities are C calculated C C - P, Q, and R-branch transitions are supported C C - Uppper limits on J and I are JMAX and 9/2 resp. C C This program is (very) loosely based on the program ASQUP by A.Travis, C Dept. of Chemistry, University of Exeter, the calculation is straight out C of Gordy & Cook and Townes & Schawlow C C C ver. 14.X.2005 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C C Modification history: C c 1989: Creation c 29.05.97: Use of centre frequency and annotated input file c 6.01.99: ICONTR=-1 option c 14.10.05: General debugging and updating of code C C_____________________________________________________________________________ C C Data file is to have the following structure: C C LINE.1: COMMENT C LINE.2: kappa: -0.731 C LINE.3: spin*2: 3 C LINE.4: chiaa: -53.397 C LINE.5: chibb: 25.112 C LINE.6: output (0,1,2): 0 C LINE.7:====================------------------------------ C LINE.8: J' K-1' K+1' J'' K-1'' K+1'' Cutoff Fcent (use -ve cutoff for C LINE.9:-------------------------------------------------- freq. range of plot) C LINE.10: J' K-1' K+1' J'' K-1'' K+1'' CUTOFF FCENT C . . . C . . . C LINE.N: -1 C c - numerical values in lines 2-7 are to be reedited as necessary, and C have to start not earlier than the 21st column C C - KAPPA, CHIAA, CHIBB are to be input as REAL, the rest are INTEGER's C C - CUTOFF is the intensity cutoff for autoscaling of frequency scale C of the stick pattern display (a negative value of this parameter C specifies a fixed, symmetric frequency range) C c - FCENT is the central frequency of the multiplet, to which calculated c frequency shifts are added C C - output option ICONTR determines the amount of output c =-1 only the sticks C = 0 sorted components and sticks only C = 1 as above and also unsorted components and dE/dKappa C = 2 as 1 and also contributions of X.aa and X.bb to the energy levels C C - lines 7-9 are ignored by the program but should not be changed since C they serve as header lines to the information that follows C C - line.10 contains transitions for which the quadrupole structure C is to be calculated. Free format is allowed. This line can be C followed by any number of similar lines, until terminated by C line N containing -1 C C NOTE that in the output for P-types quantum numbers are reversed C so watch the directions of the transition arrowa C C C The program is semi-interactive in that it will ask for the input and C output file names and display some intermediate information on the C terminal. C_____________________________________________________________________________ C C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (JMAX=120, MAXDIM=2*JMAX+1) COMMON /HBLK/H/SORTCC/WK,IORDER DIMENSION J2(2),P2(2,3),COEFF(2,2) INTEGER*4 T(2) REAL*8 JQN,FQN,IQN DIMENSION ILAB(2,10),RES(2,2,10),CO(MAXDIM,2),WK(MAXDIM) DIMENSION H(MAXDIM,MAXDIM),E(MAXDIM,2),FREQ(MAXDIM), * RELINT(MAXDIM),Y(10) INTEGER*4 QNO(MAXDIM,2),IORDER(MAXDIM),NF(2) CHARACTER FILNAM*30,COMENT*78,DIGIT(10),BASE(79),STICKS(79,4), * dummy*20,arrow*2,line*80 C DATA DIGIT/'0','1','2','3','4','5','6','7','8','9'/ DELTAK=0.00001D0 C C...INPUT C WRITE(*,5500) 5500 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | ASQ1P - First order prediction for asymmetric rotor ', * 'with one',T79,'|'/ * ' | quadrupolar nucleus',T79,'|'/ * ' |',76(1H_),'|'/' version 14.X.2005',T64,'Zbigniew KISIEL'/) C 5502 WRITE(*,5501)'Name of data file: ' 5501 FORMAT(1X/1X,A,$) READ(*,'(A)',ERR=5502)FILNAM OPEN(2,FILE=FILNAM,status='old',ERR=5502) READ(2,'(A)' ,ERR=3002)COMENT READ(2,'(a,f20.10)',ERR=3002) dummy,AKAPPA READ(2,'(a,i20 )',ERR=3002) dummy,IA read(2,'(a,f20.10)',ERR=3002) dummy,CHIA read(2,'(a,f20.10)',ERR=3002) dummy,CHIB read(2,'(a,i20 )',ERR=3002) dummy,ICONTR read(2,'(a/a/a)', ERR=3002)dummy,dummy,dummy IQN=FLOAT(IA)*0.5D0 C 5503 WRITE(*,5501)'Name of output file: ' READ(*,'(A)',ERR=5503)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=5503) C WRITE(3,'(1X,78(1H_)//1X,A/1X,78(1H_)//)')COMENT WRITE(3,5523) 5523 FORMAT(' HYPERFINE STRUCTURE PREDICTION FOR AN ASYMMETRIC ', 1 'TOP WITH ONE QUADRUPOLAR ATOM'/ 1 ' (First order calculation only)'/) WRITE(3,'(1X,A,I2,''/2''/)')' Spin = ',IA WRITE(3,5522)' Kappa =',AKAPPA WRITE(3,5522)' X.aa =',CHIA WRITE(3,5522)' X.bb =',CHIB CHIC=-(CHIA+CHIB) WRITE(3,5522)' X.cc =',CHIC 5522 FORMAT(1X,A,F12.6) C C C...MAIN LOOP OVER TRANSITIONS C C INITIALISATION OF ARRAYS C WRITE(*,5510) 5510 FORMAT(1X/' NOW WORKING ON TRANSITION:'/) C C if(icontr.eq.0)WRITE(3,17) 17 FORMAT (1X/20X,' C(X.aa)',5X,' C(X.bb)',5X,'Frequency',3X, * 'Relative Intens.') C 5505 READ(2,'(A)',err=3001,end=3001)line READ(LINE,*,ERR=3003,end=3003) J2(1),KM1U,KP1U,J2(2),KM1L,KP1L, * CUTOFF,FCENT goto 3004 c 3003 read(line,*,err=3001)jjj IF(jjj.LT.0)GOTO 5504 goto 3001 c 3004 arrow='<-' pqr=1.d0 c c...set up reversal for P-types if(j2(1).lt.j2(2))then n=J2(2) J2(2)=J2(1) J2(1)=n n=km1u km1u=km1l km1l=n n=kp1u kp1u=kp1l kp1l=n arrow='->' pqr=-1.d0 endif WRITE(3,10)J2(1),KM1U,KP1U,arrow,J2(2),KM1L,KP1L,FCENT 10 FORMAT(1X//1X,I3,',',I3,',',I3,2x,a2,I4,',',I3,',',I3, * ' F.cent = ',F13.6/) WRITE(*,5511) J2(1),KM1U,KP1U,arrow,J2(2),KM1L,KP1L 5511 FORMAT(15X,I3,',',I3,',',I3,2x,a2,I4,',',I3,',',I3) C T(1)=KM1U-KP1U T(2)=KM1L-KP1L JQN=J2(1) TRM1=JQN+IQN+1.0D0 TRM2=JQN+IQN TRM3=JQN-IQN-1.0D0 TRM4=JQN-IQN TRM5=JQN*(JQN+1.0D0) TRM6=IQN*(IQN+1.0D0) C C C...SET UP THE HAMILTONIAN (The incremental method of evaluating the C derivatives wrt. kappa is numerically unstable and is critically c dependent on the value of DELTAK - should really be replaced by c use of eigenvectors) C DO 1001 MN=1,2 J=J2(MN) JJP1=J*(J+1) L=2*J+1 C DO 2001 N=1,2 DO 1000 I=1,L DO 1000 JJ=1,L H(I,JJ)=0.0D0 1000 CONTINUE AN=AKAPPA IF (N.EQ.2)AN=AKAPPA+DELTAK TERM1=0.5D0*(AN-1.D0) TERM2=-0.5D0*(AN+1.D0) DO 2000 KK=1,L K=KK-(J+1) H(KK,KK)=TERM1*(JJP1-K**2)+K**2 2000 CONTINUE DO 3000 KK=1,L-2 K=KK-(J+1) H(KK,KK+2)=TERM2*DSQRT(0.25D0*(JJP1-K*(K+1)) 1 *(JJP1-(K+1)*(K+2))) H(KK+2,KK)=H(KK,KK+2) 3000 CONTINUE C C...DIAGONALIZE THE HAMILTONIAN C CALL JACK(L) DO 2001 I=1,L E(I,N)=H(I,I) 2001 CONTINUE C C...CALCULATE dE(kappa)/dkappa C DO 5000 LM=1,L WK(LM)=E(LM,1) IORDER(LM)=LM E(LM,2)=(E(LM,2)-E(LM,1))/DELTAK 5000 CONTINUE C C...SORT THE EIGENVALUES C CALL SORTC(1,L) C C...SELECT THE REQUIRED EIGENVALUE AND CALCULATE PA2,PB2,PC2 C NUM=T(MN)+J+1 IF(NUM.GT.L)GOTO 308 EKAPPA=WK(NUM) DEKAPA=E(IORDER(NUM),2) C IF(icontr.gt.0)then IF(MN.EQ.1)WRITE(3,5517)J2(1),KM1U,KP1U,EKAPPA,DEKAPA IF(MN.EQ.2)WRITE(3,5517)J2(2),KM1L,KP1L,EKAPPA,DEKAPA endif 5517 FORMAT(1X,I2,',',I2,',',I2,': E(kappa) =',F12.6,' dE/dk =', 1 F12.6) C C Eq.(9.97)-(9.99) G&CIII P2(MN,1)=0.5D0*(JJP1+EKAPPA-(AKAPPA+1.0D0) 5 *DEKAPA) P2(MN,2)=DEKAPA P2(MN,3)=0.5D0*(JJP1-EKAPPA+(AKAPPA-1.0D0) 5 *DEKAPA) 1001 CONTINUE GO TO 310 C 308 WRITE(3,309)J2(1),KM1U,KP1U,J2(2),KM1L,KP1L 309 FORMAT(1H ,'**** ERROR in quantum numbers: ',I3,',',I3, 1 ',',I3,' and ',I3,',',I3,',',I3/) GOTO 5505 C C...OBTAIN COEFFICIENTS OF XAA AND XBB in Eq.(9.100) G&CIII C 310 DO 2 K=1,2 TERM1=J2(K)*(J2(K)+1) IF(J2(K).EQ.0)TERM1=1.D50 TERM1=2.0D0/TERM1 COEFF(K,1)=TERM1*(P2(K,1)-P2(K,3)) COEFF(K,2)=TERM1*(P2(K,2)-P2(K,3)) 2 CONTINUE C C...WRITE OUT RESULTS FOR ENERGY LEVELS C C OBTAIN COEFFICIENTS FOR EACH F C IF(ICONTR.LT.2)GOTO 5506 WRITE(3,5507)J2(1),KM1U,KP1U,J2(2),KM1L,KP1L 5507 FORMAT(1X//' Coefficients of X.aa and X.bb and Y(I,J,F) for ', 1 I2,',',I2,',',I2,' and ',I2,',',I2,',',I2/) C 5506 DO 5 K=1,2 J=J2(K)*2 MAXF=IA+J MINF=IABS(IA-J) NF(K)=(MAXF-MINF)/2+1 N=1 DO 3 IF=MINF,MAXF,2 Y(N)=CASMIR(J,IA,IF) ILAB(K,N)=IF RES(K,1,N)=COEFF(K,1)*Y(N) RES(K,2,N)=COEFF(K,2)*Y(N) N=N+1 3 CONTINUE IF(ICONTR.LT.2)GOTO 5 DO 22 N=1,NF(K) WRITE(3,9)J2(K),ILAB(K,N),(RES(K,L,N),L=1,2),Y(N) 22 CONTINUE 9 FORMAT (4X,I2,3X,I3,'/2',6X,F9.6,4X,F9.6,F12.6) WRITE(3,5521) 5 CONTINUE C C...CALCULATE THE TRANSITIONS, Relative intensities are according to C Eq. 6-6a,6-6b of T&S C IF(ICONTR.GT.2)WRITE(3,17) C C...F <- F-1 C SUMINT=0. LNO=0 NFPLUS=0 JSELEC=J2(1)-J2(2) DO 13 M=1,NF(2) DO 13 N=1,NF(1) IF((ILAB(1,N)-ILAB(2,M)).NE.2)GOTO 13 NFPLUS=NFPLUS+1 LNO=LNO+1 CO(LNO,1)=RES(1,1,N)-RES(2,1,M) CO(LNO,2)=RES(1,2,N)-RES(2,2,M) FREQ(LNO)=pqr*( (CO(LNO,1)*CHIA)+(CO(LNO,2)*CHIB) ) QNO(LNO,1)=ILAB(1,N) QNO(LNO,2)=ILAB(2,M) FQN=FLOAT(ILAB(1,N))*0.5D0 if(fqn.eq.0.d0)then relint(lno)=0.0d0 goto 5542 endif IF(JSELEC.eq.0)THEN RELINT(LNO)=-((TRM1+FQN)*(TRM4+FQN)*(TRM1-FQN)*(TRM4-FQN)) 1 /FQN ELSE RELINT(LNO)=((TRM1+FQN)*(TRM2+FQN)*(TRM4+FQN)*(TRM3+FQN)) 1 /FQN ENDIF 5542 SUMINT=SUMINT+RELINT(LNO) 13 CONTINUE C C...F <- F C NFZERO=0 DO 16 M=1,NF(2) DO 16 N=1,NF(1) IF((ILAB(1,N)-ILAB(2,M)).NE.0)GOTO 16 NFZERO=NFZERO+1 LNO=LNO+1 CO(LNO,1)=RES(1,1,N)-RES(2,1,M) CO(LNO,2)=RES(1,2,N)-RES(2,2,M) FREQ(LNO)=PQR*( (CO(LNO,1)*CHIA)+(CO(LNO,2)*CHIB) ) QNO(LNO,1)=ILAB(1,N) QNO(LNO,2)=ILAB(2,M) FQN=FLOAT(ILAB(1,N))*0.5D0 FFP1=FQN*(FQN+1.D0) if(ffp1.eq.0.d0)then relint(lno)=0.0d0 goto 5545 endif IF(JSELEC.eq.0)then RELINT(LNO)=((TRM5-TRM6+FFP1)**2*(2.D0*FQN+1.D0))/FFP1 else RELINT(LNO)=-((TRM1+FQN)*(TRM4+FQN)*(TRM2-FQN)*(TRM3-FQN) 1 *(2.D0*FQN+1.D0))/FFP1 endif 5545 SUMINT=SUMINT+RELINT(LNO) 16 CONTINUE C C...F <- F+1 C DO 5509 M=1,NF(2) DO 5509 N=1,NF(1) IF((ILAB(1,N)-ILAB(2,M)).NE.-2)GOTO 5509 LNO=LNO+1 CO(LNO,1)=RES(1,1,N)-RES(2,1,M) CO(LNO,2)=RES(1,2,N)-RES(2,2,M) FREQ(LNO)=PQR*( (CO(LNO,1)*CHIA)+(CO(LNO,2)*CHIB) ) QNO(LNO,1)=ILAB(1,N) QNO(LNO,2)=ILAB(2,M) FQN=FLOAT(ILAB(1,N))*0.5D0 if(fqn+1.d0.eq.0.d0)then relint(lno)=0.0d0 goto 5548 endif IF(JSELEC.eq.0)then RELINT(LNO)=-((TRM1+1.D0+FQN)*(TRM4+1.D0+FQN) 1 *(TRM1-1.D0-FQN) 1 *(TRM4-1.D0-FQN))/(FQN+1.D0) else RELINT(LNO)=((TRM2-FQN)*(TRM1-2.D0-FQN)*(TRM3-FQN) 1 *(TRM3-1.D0-FQN))/(FQN+1.D0) endif 5548 SUMINT=SUMINT+RELINT(LNO) 5509 CONTINUE C C...WRITE OUT TRANSITIONS (GROUPED IN F SELECTION RULE AND ALSO SORTED C ACCORDING TO FREQUENCY) C DO 5515 N=1,LNO RELINT(N)=RELINT(N)/SUMINT WK(N)=FREQ(N) IORDER(N)=N IF(ICONTR.GT.0)THEN WRITE(3,104) QNO(N,1),arrow,QNO(N,2),(CO(N,J),J=1,2), * FREQ(N)+FCENT,RELINT(N) 104 FORMAT(3X,I3,'/2 ',a2,I3,'/2',1X,4F13.6) IF(N.EQ.NFPLUS)WRITE(3,5521) IF(N.EQ.NFPLUS+NFZERO)WRITE(3,5521) ENDIF 5515 CONTINUE C C...sorted frequencies IF(ICONTR.gt.0)THEN WRITE(3,5521) WRITE(3,5521) ENDIF 5521 FORMAT(1X) CALL SORTC(1,LNO) DO 5520 N=1,LNO M=IORDER(N) if(icontr.ge.0)WRITE(3,104) * QNO(M,1),arrow,QNO(M,2),(CO(M,J),J=1,2), * FREQ(M)+FCENT,RELINT(M) 5520 CONTINUE C C...baseline for sticks NPLINE=79 IF(CUTOFF.GE.0.0D0)THEN FMIN=100.D0 FMAX=-100.D0 DO 806 N=1,LNO IF(RELINT(N).LT.CUTOFF)GOTO 806 IF(FREQ(N).LT.FMIN)FMIN=FREQ(N) IF(FREQ(N).GT.FMAX)FMAX=FREQ(N) 806 CONTINUE RANGE=(FMAX-FMIN) fmax=fmax+0.05D0*RANGE FMIN=FMIN-0.05D0*RANGE ELSE FMIN=CUTOFF FMAX=-CUTOFF ENDIF C RANGE=(FMAX-FMIN) SCALE=NPLINE/RANGE DO 800 N=1,NPLINE BASE(N)='-' STICKS(N,1)=' ' STICKS(N,2)=' ' STICKS(N,3)=' ' STICKS(N,4)=' ' 800 CONTINUE C MARK=1 801 NMARKS=RANGE/MARK IF(NMARKS.GT.12)THEN MARK=MARK*2 GOTO 801 ENDIF DO 802 N=INT(FMIN),INT(FMAX) IF(N.EQ.MARK*(N/MARK))THEN NNN=(N-FMIN)*SCALE+1 NNNN=IABS(N-10*(N/10)) IF(NNN.LT.1.OR.NNN.GT.NPLINE)GOTO 802 BASE(NNN)=DIGIT(NNNN+1) ENDIF 802 CONTINUE c c...sticks DO 803 N=1,LNO NN=(FREQ(N)-FMIN)*SCALE+1 IF(NN.LT.1.OR.NN.GT.NPLINE)GOTO 803 IF(RELINT(N).GE.0.3D0) STICKS(NN,1)='|' IF(RELINT(N).GE.0.2D0) STICKS(NN,2)='|' IF(RELINT(N).GE.0.1D0) STICKS(NN,3)='|' IF(RELINT(N).GE.0.02D0)STICKS(NN,4)='|' IF(RELINT(N).LT.0.02D0)STICKS(NN,4)='.' 803 CONTINUE C WRITE(3,5521) WRITE(3,804)(STICKS(N,1),N=1,NPLINE) WRITE(3,804)(STICKS(N,2),N=1,NPLINE) WRITE(3,804)(STICKS(N,3),N=1,NPLINE) WRITE(3,804)(STICKS(N,4),N=1,NPLINE) WRITE(3,804)BASE 804 FORMAT(79A1) WRITE(3,805)FMIN,FMAX 805 FORMAT('Baseline between ',F5.1,' and',F5.1) C GOTO 5505 C 5504 CLOSE(2) CLOSE(3) stop C C...Error conditions C 3002 WRITE(*,'(1x//'' **** ERROR on reading top datafile lines''//)') CLOSE(2) STOP C 3001 WRITE(*,'(1x// * '' **** ERROR reading transition or end of datafile''//)') CLOSE(2) CLOSE(3) STOP C STOP END C C_____________________________________________________________________________ C FUNCTION CASMIR(J,I,F) C C FUNCTION TO DETERMINE CASIMIRS FUNCTION C (NOTE: function parameters are twice the values of the quantum numbers) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER*4 F C C=0.25D0*(F*(F+2)-J*(J+2)-I*(I+2)) DENOM=(J-1)*(J+3)*I*(I-1) CASMIR=(0.75D0*C*(C+1.D0)-0.0625D0*(I*(I+2)*J*(J+2)))/DENOM C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SORTC(N,M) PARAMETER (JMAX=120, MAXDIM=2*JMAX+1) COMMON /SORTCC/WK,IPT INTEGER*4 IPT(MAXDIM) REAL*8 WK(MAXDIM),EE C C ... This routine sorts the quantities in vector WK in ascending order C of magnitude and also accordingly rearranges vector IPT of pointers C to original positions of sorted quantities C DO 101 I=N,M-1 J=I 106 J=J+1 IF(WK(J)-WK(I))103,104,104 103 EE=WK(I) WK(I)=WK(J) WK(J)=EE K=IPT(I) IPT(I)=IPT(J) IPT(J)=K 104 IF(J.EQ.M)GOTO 101 GOTO 106 101 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE JACK(N) c c JACOBI type matrix diagonalisation routine for real, symmetric matrices c c H - array to be diagonalised (contains the eigenvalues on the diagonal c on output c S - contains the eigenvectors on output (in order of eigenvalues) c N - dimension of the array to be diagonalised contained in the top c left hand corner of H c c IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (JMAX=120, MAXDIM=2*JMAX+1) COMMON /HBLK/H REAL*8 H(MAXDIM,MAXDIM),S(MAXDIM,MAXDIM) c c ... initialisation and setup of limits c ss=0.7071067811865475d0 n1=n-1 sum=0.d0 do 20 i=1,n1 s(i,i)=1.d0 k=i+1 do 20 m=k,n sum=sum+h(i,m)**2 s(i,m)=0.d0 s(m,i)=0.d0 20 continue s(n,n)=1.d0 unu=dsqrt(2.d0*sum) ene=n c c ... convergence criterion c unufi=1.d-12 c c c ... iteration sequence - determination of SIN and COS for rotation c 54 unu=unu/ene 55 ind=0 do 89 m=2,n l=m-1 do 89 i=1,l if(h(i,m))60,89,60 60 q=dabs(h(i,m)) if(q-unu)89,62,62 62 ind=1 if(h(i,i)-h(m,m))68,64,68 64 si=dsign(ss,h(i,m)) co=ss goto 77 68 x1=-h(i,m)*2.d0/(h(i,i)-h(m,m)) x1a=dabs(x1) if(1.d30-x1a)66,66,71 66 si=dsign(ss,x1) co=ss goto 77 71 if(1.d-30-x1a)72,89,89 72 t=datan(x1)*0.5d0 si=dsin(t) co=dcos(t) c c ... transformation loop c 77 do 79 j=1,n hji=h(j,i) sji=s(j,i) h(j,i)=hji*co-h(j,m)*si h(j,m)=hji*si+h(j,m)*co s(j,i)=sji*co-s(j,m)*si s(j,m)=sji*si+s(j,m)*co if(j.eq.i.or.j.eq.m)goto 79 h(i,j)=h(j,i) h(m,j)=h(j,m) 79 continue c h(i,i)=co*h(i,i)-si*h(m,i) h(m,m)=co*h(m,m)+si*h(i,m) h(i,m)=0.d0 h(m,i)=0.d0 89 continue if(ind)92,91,55 91 if(unu-unufi)92,92,54 c 92 return end C C_____________________________________________________________________________ C_____________________________________________________________________________