344 lines
9.5 KiB
Fortran
344 lines
9.5 KiB
Fortran
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
|