SpectraRust/tlusty/extracted/inifrt.f
2026-03-19 14:05:33 +08:00

264 lines
7.9 KiB
Fortran

SUBROUTINE INIFRT
C =================
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
C
DIMENSION FREQCO(MFREQC),WCO(MFREQC),FRLEV(MLEVEL),IENS(MLEVEL)
dimension WCHCO(MFREQC),IJXCO(MFREQC)
common/ijflar/ijfl(mlevel)
C
DFEDG=0.01
NEND=NFTAIL
DIVEND=DFTAIL
NJC=NFREQC/5
DNX=UN-UN/FLOAT(NJC)
CALL INDEXX(NLEVEL,ENION,IENS)
DO IL=1,NLEVEL
ILS=IENS(NLEVEL-IL+1)
FRLEV(IL)=ENION(ILS)/H
IJFL(IL)=0
END DO
C
IF(FRCMAX.LE.0.) THEN
FREQCO(1)=8.E11*TEFF
ELSE
FREQCO(1)=FRCMAX
END IF
C
IL0=1
IF(FREQCO(1).LT.CFRMAX*FRLEV(IL0) .AND. CFRMAX.GT.UN) THEN
FREQCO(1)=CFRMAX*FRLEV(IL0)
GO TO 10
END IF
DO WHILE (FREQCO(1).LT.FRLEV(IL0) .AND. IL0.LT.NLEVEL)
ILS=IENS(NLEVEL-IL0+1)
ILN=NNEXT(IEL(ILS))
ITR0=ITRA(ILS,ILN)
INDEXP(ITR0)=0
IFR0(ITR0)=0
IFR1(ITR0)=0
IF(FRLEV(IL0).GT.0.) WRITE(10,159) IL0,FRLEV(IL0),ILS,ILN,ITR0
159 FORMAT(' Edge at frequency larger than FRCMAX',I5,E12.4,3I7)
IL0=IL0+1
END DO
10 CONTINUE
C
NFREQC=NEND+1
FREQCO(NEND)=(un+dfedg)*FRLEV(il0)
FREQCO(NEND+1)=(un-dfedg)*FRLEV(il0)
NEND1=NEND/2+1
XEND=UN/FLOAT(NEND1-1)
FREQCO(NEND1)=FREQCO(1)-(UN-DIVEND)*(FREQCO(1)-FREQCO(NEND))
IJXCO(NEND+1)=1
C
C high-frequency tail of the spectrum - a two-part Simpson integration
C
C 1st part - from the highest frequency FRCMAX to a division freq.
C
IJXCO(1)=1
IJXCO(NEND1)=1
D121=XEND*(FREQCO(1)-FREQCO(NEND1))
if(icompt.gt.0) d121=(freqco(1)/freqco(nend1))**xend
DO IJ=2,NEND1-1
FREQCO(IJ)=FREQCO(IJ-1)-D121
if(icompt.gt.0) freqco(ij)=freqco(ij-1)/d121
IJXCO(IJ)=2
END DO
D121=THIRD*(FREQCO(1)-FREQCO(2))
DO IJ=2,NEND1-1,2
WCO(IJ)=4.*D121
WCO(IJ-1)=WCO(IJ-1)+D121
WCO(IJ+1)=WCO(IJ+1)+D121
WCHCO(IJ)=0.
END DO
C
C 2nd part - from the division freq to the first discontinuity
C
if(nend1.lt.nend) then
IJXCO(NEND)=1
IJXCO(NEND+1)=1
D121=XEND*(FREQCO(NEND1)-FREQCO(NEND))
if(icompt.gt.0) d121=(freqco(nend1)/freqco(nend))**xend
DO IJ=NEND1+1,NEND-1
FREQCO(IJ)=FREQCO(IJ-1)-D121
if(icompt.gt.0) freqco(ij)=freqco(ij-1)/d121
IJXCO(IJ)=2
END DO
D121=THIRD*(FREQCO(NEND1)-FREQCO(NEND1+1))
DO IJ=NEND1+1,NEND-1,2
WCO(IJ)=4.*D121
WCO(IJ-1)=WCO(IJ-1)+D121
WCO(IJ+1)=WCO(IJ+1)+D121
WCHCO(IJ)=0.
END DO
end if
C
C the 1st discontinuity - the one with the highest frequency
C
HAEND=HALF*(FREQCO(NEND)-FREQCO(NEND+1))
XEND=UN/FLOAT(NEND-1)
WCO(NEND)=WCO(NEND)+HAEND
WCO(NEND+1)=WCO(NEND+1)+HAEND
WCHCO(1)=0.
WCHCO(NEND)=HAEND
ILS=IENS(NLEVEL)
IJFL(ILS)=NEND
IL0=NLEVEL
IF(FRCMIN.LE.0.) FRCMIN=1.D12
FRCLST=FRLEV(IL0)
DO WHILE(FRCLST.LT.FRCMIN)
IF(FRLEV(IL0).GT.0.) WRITE(10,160) IL0,FRLEV(IL0)
160 FORMAT(' Edge at frequency smaller than 1.d12',i5,e12.4)
IL0=IL0-1
FRCLST=FRLEV(IL0)
END DO
IL0=2
C
100 CONTINUE
FRC0=DNX*FREQCO(NFREQC)
IF(FRC0.LT.FRCLST) THEN
NFREQC=NFREQC+2
FREQCO(NFREQC-1)=(un+dfedg)*FRCLST
FREQCO(NFREQC)=(un-dfedg)*FRCLST
IJXCO(NFREQC-1)=1
IJXCO(NFREQC)=1
WCO(NFREQC)=WCO(NFREQC)+HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCO(NFREQC-1)=WCO(NFREQC-1)+
* HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC))
WCO(NFREQC-2)=WCO(NFREQC-2)+
* HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC-1))
WCHCO(NFREQC-1)=HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCHCO(NFREQC-2)=HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC-1))
ILS=IENS(NLEVEL-IL0+1)
IJFL(ILS)=NFREQC-1
IF(IL0.LT.NLEVEL) THEN
DO IL=IL0+1,NLEVEL
ILS=IENS(NLEVEL-IL+1)
IJFL(ILS)=NFREQC-1
IF(FRLEV(IL).LT.FRCMIN) IJFL(ILS)=0
END DO
END IF
D121=XEND*(FREQCO(NFREQC)-FRCMIN)
DO IJ=NFREQC+1,NFREQC+NEND-1
FREQCO(IJ)=FREQCO(IJ-1)-D121
IJXCO(IJ)=2
WCHCO(IJ)=0.
END DO
IJXCO(NFREQC+NEND-1)=1
DO IJ=NFREQC+1,NFREQC+NEND-2,2
WCO(IJ)=FTH*D121
WCO(IJ-1)=WCO(IJ-1)+THIRD*D121
WCO(IJ+1)=WCO(IJ+1)+THIRD*D121
END DO
WCHCO(NFREQC)=THIRD*D121
NFREQC=NFREQC+NEND-1
GO TO 200
END IF
DF0=FRLEV(IL0)+0.1*(FREQCO(NFREQC)-FRC0)
FRTL=(UN+DFEDG)*FRLEV(IL0)
IF(FRC0.GT.DF0) THEN
NFREQC=NFREQC+1
FREQCO(NFREQC)=FRC0
IJXCO(NFREQC)=2
WCO(NFREQC)=WCO(NFREQC)+HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCO(NFREQC-1)=WCO(NFREQC-1)+
* HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCHCO(NFREQC-1)=HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
ELSE IF(FRTL.LT.FREQCO(NFREQC)) THEN
NFREQC=NFREQC+2
FREQCO(NFREQC-1)=FRTL
FREQCO(NFREQC)=(un-dfedg)*FRLEV(IL0)
IJXCO(NFREQC-1)=1
IJXCO(NFREQC)=1
WCO(NFREQC)=WCO(NFREQC)+HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCO(NFREQC-1)=WCO(NFREQC-1)+
* HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC))
WCO(NFREQC-2)=WCO(NFREQC-2)+
* HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC-1))
WCHCO(NFREQC-1)=HALF*(FREQCO(NFREQC-1)-FREQCO(NFREQC))
WCHCO(NFREQC-2)=HALF*(FREQCO(NFREQC-2)-FREQCO(NFREQC-1))
ILS=IENS(NLEVEL-IL0+1)
IJFL(ILS)=NFREQC-1
IL0=IL0+1
ELSE
ILS=IENS(NLEVEL-IL0+1)
IJFL(ILS)=NFREQC-1
IL0=IL0+1
END IF
GO TO 100
200 CONTINUE
C
nend2=nfreqc
IF(IFRSET.EQ.0) THEN
DO IJ=1,NUMFREQ
FREQCO(IJ+NEND2)=FRTAB(IJ)
END DO
NFR=NUMFREQ+NEND2
ELSE IF(IFRSET.GT.0) THEN
NFRTAB=IFRSET
DO IJ=1,NFRTAB
FR=LOG(FRTAB(NUMFREQ))+
* LOG(FRTAB(1)/FRTAB(NUMFREQ))*(IJ-1)/(NFRTAB-1)
FREQCO(NEND2+NFRTAB-IJ+1)=EXP(FR)
NFR=NFRTAB+NEND2
END DO
END IF
C
IF(FRTAB(NUMFREQ).GT.FRCMIN) THEN
DO IJ=1,NEND
FR=LOG(FRCMIN)+(LOG((UN-DFEDG)*FRTAB(NUMFREQ)/FRCMIN)*
* (IJ-1)/(NEND-1))
FREQCO(NFR+NEND-IJ+1)=EXP(FR)
END DO
NFR=NFR+NEND
END IF
NFREQC=NFR
C
WCO(1)=0.5*(FREQCO(1)-FREQCO(2))
WCO(NFR)=0.5*(FREQCO(NFR-1)-FREQCO(NFR))
DO IJ=2,NFR-1
WCO(IJ)=0.5*(FREQCO(IJ-1)-FREQCO(IJ+1))
END DO
C
DO IJ=NFREQ,1,-1
FREQ(IJ+NFREQC)=FREQ(IJ)
W(IJ+NFREQC)=W(IJ)
PROF(IJ+NFREQC)=PROF(IJ)
IJALI(IJ+NFREQC)=IJALI(IJ)
END DO
DO IJ=1,NFREQC
FREQ(IJ)=FREQCO(IJ)
W(IJ)=WCO(IJ)
IJALI(IJ)=1
IJX(IJ)=1
PROF(IJ)=0.
END DO
DO 20 ITR=1,NTRANS
IF(.NOT.LINE(ITR).OR.INDEXP(ITR).EQ.0) GO TO 20
IF(IFR0(ITR).GT.0) IFR0(ITR)=IFR0(ITR)+NFREQC
IF(IFR1(ITR).GT.0) IFR1(ITR)=IFR1(ITR)+NFREQC
20 CONTINUE
NFREQ=NFREQ+NFREQC
C
C determination of the first (IFR0) and the last (IFR1) frequency
C for each explicit continuum - only if they were not already read
C
DO 30 ITR=1,NTRANS
IF(LINE(ITR)) GO TO 30
IF(IFR0(ITR).LE.0) IFR0(ITR)=1
IF(IFR1(ITR).GT.0) GO TO 30
IF1=0
FR01=FR0(ITR)
MODE=INDEXP(ITR)
IF(IABS(MODE).EQ.5.OR.IABS(MODE).EQ.15) FR01=FR0PC(ITR)
DO IJ=1,NFREQC
IF(FREQ(IJ).GE.FR01) IF1=IJ
END DO
IFR1(ITR)=IF1
30 CONTINUE
RETURN
END