SUBROUTINE LINSEL C ================= C C Exclude weakest lines C (i.e. set them to detailed radiative balance, and C exclude their frequencies) C C Selection of lines based on ratio of line core flux to C continuum flux C C Non-standard parameters: STRL1 (default value 0.001) C STRL2 ( 0.02) C C STRL2 allows to reduce the number of frequency points in C "intermediate-strength" lines C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ODFPAR.FOR' INCLUDE 'ALIPAR.FOR' DIMENSION PRFTMP(MDEPTH) C parameter (SIXTH=UN/6.,FTH=4./3.) C NLSW=0 NLSI=0 NLSS=0 NLSTO=0 C IF(ISPODF.GE.1) GO TO 190 C CALL OPAINI(1) C C Normal lines C DO 10 ITR=1,NTRANS IF(LINEXP(ITR)) GO TO 10 IF(LEXP(ITR) .OR. IEL(ILOW(ITR)).EQ.IELH) THEN NLSTO=NLSTO+1 NLSS=NLSS+1 GO TO 10 ENDIF MODE=IABS(INDEXP(ITR)) IF(MODE.GE.2 .AND. MODE.LE.4) GO TO 10 IKA=IFR0(ITR) IF(IKA.EQ.0) GO TO 10 IKB=IFR1(ITR) CALL OPACF1(IKA) CALL RTEFR1(IKA) FLUXA=FH(IKA)*RAD1(1) CALL OPACF1(IKB) CALL RTEFR1(IKB) FLUXB=FH(IKB)*RAD1(1) IK0=(IKA+IKB)/2 IF(MODE.EQ.2) IK0=IFR1(ITR)-1 CALL OPACF1(IK0) CALL RTEFR1(IK0) FLUX0=FH(IK0)*RAD1(1) RHAB=UN-TWO*FLUX0/(FLUXA+FLUXB) NFK0=IKB-IKA+1 NLSTO=NLSTO+1 IF(ABS(RHAB).LT.STRL1) THEN NLSW=NLSW+1 LINEXP(ITR)=.TRUE. INDEXP(ITR)=0 DO IJ=IFR0(ITR),IFR1(ITR) IJX(IJ)=-1 IJLIN(IJ)=0 END DO IFR0(ITR)=0 IFR1(ITR)=0 KFR0(ITR)=0 KFR1(ITR)=0 LALI(ITR)=.FALSE. ELSE IF(ABS(RHAB).LT.STRL2) THEN NLSI=NLSI+1 IFR0(ITR)=IKA+3 IFR1(ITR)=IKB-3 KFR0(ITR)=KIJ(IFR0(ITR)) KFR1(ITR)=KIJ(IFR1(ITR)) DO IJ=IKA,IFR0(ITR) IJX(IJ)=-1 IJLIN(IJ)=0 DO ID=1,ND PRFLIN(ID,IJ)=0. END DO END DO IJX(IFR0(ITR))=1 IJLIN(IFR0(ITR))=ITR DO IJ=IFR1(ITR),IKB IJX(IJ)=-1 IJLIN(IJ)=0 DO ID=1,ND PRFLIN(ID,IJ)=0. END DO END DO IJX(IFR1(ITR))=1 IJLIN(IFR1(ITR))=ITR ELSE NLSS=NLSS+1 ENDIF NFK1=IFR1(ITR)-IFR0(ITR) IF(NFK1.GT.0) NFK1=NFK1+1 WRITE(82,601) ITR,ILOW(ITR),IUP(ITR), * NFK0,NFK1,RHAB 601 FORMAT(I6,2I5,2I7,1PE12.4) 10 CONTINUE C C superlines C DO 20 ITR=1,NTRANS IF(LINEXP(ITR)) GO TO 20 MODE=IABS(INDEXP(ITR)) IF(MODE.NE.3 .AND. MODE.NE.4) GO TO 20 IF(LEXP(ITR)) THEN NLSTO=NLSTO+1 NLSS=NLSS+1 GO TO 20 ENDIF IKA=IFR0(ITR) IF(IKA.EQ.0) GO TO 20 IKB=IFR1(ITR) NFK0=IFR1(ITR)-IFR0(ITR)+1 NFK1=NFK0 RHAB=0. RHABMX=0. PRFA=PRFLIN(1,IKA+1) PRFB=PRFLIN(1,IKB-1) IF(PRFA.GT.PRFB) THEN IK2=IKB-1 DO WHILE (IK2.GT.IKA .AND. RHAB.LT.STRL1) CALL OPACF1(IK2) CALL RTEFR1(IK2) FLUX2=FH(IK2)*RAD1(1) DO ID=1,ND PRFTMP(ID)=PRFLIN(ID,IK2) PRFLIN(ID,IK2)=0. END DO CALL OPACF1(IK2) CALL RTEFR1(IK2) FLUX1=FH(IK2)*RAD1(1) RHAB=ABS(UN-FLUX2/FLUX1) IF(RHAB.GT.RHABMX) RHABMX=RHAB IK20=IK2 IF(RHAB.LT.STRL1) THEN IJX(IK2+1)=-1 IJLIN(IK2+1)=0 IK2=IK2-1 ELSE IK2=IK2+1 IFR1(ITR)=IK2 KFR1(ITR)=KIJ(IK2) DO ID=1,ND PRFLIN(ID,IK20)=real(PRFTMP(ID)) END DO END IF END DO NFK1=IFR1(ITR)-IFR0(ITR)+1 IF(IK2.EQ.IKA) NFK1=0 ELSE IK2=IKA+1 DO WHILE (IK2.LT.IKB .AND. RHAB.LT.STRL1) CALL OPACF1(IK2) CALL RTEFR1(IK2) FLUX2=FH(IK2)*RAD1(1) DO ID=1,ND PRFTMP(ID)=PRFLIN(ID,IK2) PRFLIN(ID,IK2)=0. END DO CALL OPACF1(IK2) CALL RTEFR1(IK2) FLUX1=FH(IK2)*RAD1(1) RHAB=ABS(UN-FLUX2/FLUX1) IF(RHAB.GT.RHABMX) RHABMX=RHAB IK20=IK2 IF(RHAB.LT.STRL1) THEN IJX(IK2-1)=-1 IJLIN(IK2-1)=0 IK2=IK2+1 ELSE IK2=IK2-1 IFR0(ITR)=IK2 KFR0(ITR)=KIJ(IK2) DO ID=1,ND PRFLIN(ID,IK20)=real(PRFTMP(ID)) END DO END IF END DO NFK1=IFR1(ITR)-IFR0(ITR)+1 IF(IK2.EQ.IKB) NFK1=0 END IF IF(NFK1.EQ.0) THEN NLSW=NLSW+1 LINEXP(ITR)=.TRUE. INDEXP(ITR)=0 IFR0(ITR)=0 IFR1(ITR)=0 KFR0(ITR)=0 KFR1(ITR)=0 DO IJ=IKA,IKB IJX(IJ)=-1 IJLIN(IJ)=0 END DO LALI(ITR)=.FALSE. ELSE IF(NFK1.EQ.NFK0) THEN NLSS=NLSS+1 ELSE NLSI=NLSI+1 END IF NLSTO=NLSTO+1 WRITE(82,601) ITR,ILOW(ITR),IUP(ITR), * NFK0,NFK1,RHABMX 20 CONTINUE C WRITE(6,602) NLSTO,NLSW,NLSI,NLSS 602 FORMAT(' Total number of lines :',i8,/, * ' Number of weak lines :',i8,/, * ' Intermediate lines :',i8,/, * ' Number of strong lines:',i8/) C C lines or ODFs associated with each frequency C NLIMAX=0 DO IJ=1,NFREQ NLINES(IJ)=0 DO 50 IT=1,NTRANS IF(LINEXP(IT)) GO TO 50 IF(KIJ(IJ).LT.KFR0(IT)) GO TO 50 IF(KIJ(IJ).GT.KFR1(IT)) GO TO 50 IF(IJLIN(IJ).EQ.IT) GO TO 50 NLINES(IJ)=NLINES(IJ)+1 IF(NLINES(IJ).GT.MITJ) * CALL QUIT('Too many overlappins-nlines(ij).gt.mitj', * nlines(ij),mitj) ITRLIN(NLINES(IJ),IJ)=int2(IT) 50 CONTINUE IF(NLINES(IJ).GT.NLIMAX) NLIMAX=NLINES(IJ) END DO WRITE(6,603) NLIMAX 603 FORMAT(' MAXIMUM NUMBER OF OVERLAPPING TRANSITIONS: ',I3/) C C recalculate weights for frequency integration C after the exclusion of some frequencies C NPPX=0 DO 100 IJ=1,NFREQ IF(IJX(IJ).GT.0) NPPX=NPPX+1 W(IJ)=0. KJ0=KIJ(IJ) IF(IJX(JIK(KJ0)).EQ.-1) GO TO 100 IF(KJ0.GE.2 .AND. KJ0.LT.NFREQ) THEN IK1=KJ0-1 DO WHILE (IJX(JIK(IK1)).EQ.-1) IK1=IK1-1 END DO IK2=KJ0+1 DO WHILE (IJX(JIK(IK2)).EQ.-1) IK2=IK2+1 END DO W(IJ)=HALF*ABS(FREQ(JIK(IK1))-FREQ(JIK(IK2))) ELSE IF(KJ0.EQ.1) THEN W(IJ)=HALF*ABS(FREQ(JIK(KJ0))-FREQ(JIK(KJ0+1))) ELSE IF(KJ0.EQ.NFREQ) THEN W(IJ)=HALF*ABS(FREQ(JIK(KJ0-1))-FREQ(JIK(KJ0))) END IF 100 CONTINUE C C Correction for Simpson weights C JK1=JIK(1) DO IJ=2,NFREQ,2 JK2=JIK(IJ) JK3=JIK(IJ+1) IF(IJLIN(JK2).NE.0 .OR. IJLIN(JK3).NE.0) GO TO 130 IF(WCH(JK2).NE.0.) GO TO 130 W(JK1)=W(JK1)-SIXTH*W(JK2) W(JK3)=W(JK3)-SIXTH*W(JK2) W(JK2)=W(JK2)*FTH JK1=JK3 END DO 130 JK1=JIK(NFREQ) DO IJ=NFREQ-1,3,-2 JK2=JIK(IJ) JK3=JIK(IJ-1) IF(IJLIN(JK2).NE.0 .OR. IJLIN(JK3).NE.0) GO TO 150 IF(WCH(JK2).NE.0.) GO TO 150 W(JK1)=W(JK1)-SIXTH*W(JK2) W(JK3)=W(JK3)-SIXTH*W(JK2) W(JK2)=W(JK2)*FTH JK1=JK3 END DO 150 CONTINUE C DO IJ=1,NFREQ W0E(IJ)=W(IJ)*PI4H/FREQ(IJ) IF(IJALI(IJ).GT.0) WC(IJ)=W(IJ) END DO C C check accuracy of weights for integration C 190 Z0=0. Z1=0. Z2=0. ZH=0. T1=TEFF T2=TWO*TEFF T3=HALF*TEFF X1=HK/T1 X2=HK/T2 X3=HK/T3 DO 200 IJ=1,NFREQ Z0=Z0+W(IJ) X15=FREQ(IJ)*1.D-15 BNZ=BN*X15*X15*X15 FX1=FREQ(IJ)*X1 IF(FX1.GT.100.) GO TO 200 Z1=Z1+W(IJ)*BNZ/(EXP(FREQ(IJ)*X1)-1) Z2=Z2+W(IJ)*BNZ/(EXP(FREQ(IJ)*X2)-1) ZH=ZH+W(IJ)*BNZ/(EXP(FREQ(IJ)*X3)-1) 200 CONTINUE T1S=SQRT(SQRT(0.25*Z1/SIG4P)) T1ER=T1S/T1-UN T2S=SQRT(SQRT(0.25*Z2/SIG4P)) T2ER=T2S/T2-UN T3S=SQRT(SQRT(0.25*ZH/SIG4P)) T3ER=T3S/T3-UN JK1=JIK(1) JK2=JIK(NFREQ) Z00=FREQ(JK1)-FREQ(JK2) WRITE(6,701) FREQ(JK1),FREQ(JK2),Z00,Z0,T3,T3ER,T1,T1ER,T2,T2ER 701 FORMAT(/' ACCURACY OF INTEGRATIONS:',/, * ' Interval:',1p4e16.8,/, * 15x,' Planck functions:',9x,0pf12.0,4x,1pe12.4,/, * 42x,0pf12.0,4x,1pe12.4,/,42x,0pf12.0,4x,1pe12.4,/) IF(ISPODF.GT.0) NPPX=NFREQ WRITE(6,702) NFREQ,NPPX 702 FORMAT(' TOTAL NUMBER OF FREQUENCIES:',I8,/, * ' SELECTED FREQUENCIES: ',I8) C RETURN END