fmq/SpectraRust
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
SpectraRust/tlusty/extracted/inifrc.f
Asfmq 1e30b7bc63 git push -u origin main
2026-03-19 14:05:33 +08:00

361 lines
11 KiB
Fortran
Raw Permalink Blame History

SUBROUTINE INIFRC(IALIEX)
C =========================
C
C Setup continuum frequencies, including frequencies around
C ionization limits
C IALIEX=0 : setup frequencies, all ALI
C IALIEX=1 : change IJALI for explicit frequencies
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
INCLUDE 'ODFPAR.FOR'
common/ijflar/ijfl(mlevel)
PARAMETER (THIRD=UN/3.,FTH=4./3.,TENLG = 2.302585093)
DIMENSION FRLEV(MLEVEL),IENS(MLEVEL),IJXCO(MFREQC),
* FREQCO(MFREQC),WCO(MFREQC),WCHCO(MFREQC)
C
if(ioptab.lt.0) return
C
dfedg=0.000001
if(icompt.gt.0) dfedg=0.01
DFEDG=0.01
IF(IALIEX.EQ.1) THEN
C
if(nffix.lt.0) then
do ij=1,nfreqc
ijali(ij)=0
end do
return
end if
c
IF(NFFIX.EQ.2) RETURN
DO 10 IT=1,NTRANS
IF(LINE(IT)) GOTO 10
IF(IFC1(IT).EQ.0) GOTO 10
if(iadop(iatm(ilow(it))).gt.0) go to 10
IJFL0=IJFL(ILOW(IT))+1
IF(IFC1(IT).LT.100) THEN
DO IJ=IFC0(IT),IFC1(IT)
IJFLS=IJFL0-IJ
IF(IJFLS.GE.1) THEN
IJALI(IJFLS)=0
IJX(IJFLS)=1
END IF
END DO
ELSE
DO IJ=IJFL0,1,-1
IJALI(IJ)=0
IJX(IJ)=1
END DO
END IF
10 CONTINUE
C
IF(ICOMPT.GT.0.AND.FRLCOM.GT.0) THEN
DO IJ=1,NFREQ
IF(FREQ(IJ).GT.FRLCOM) THEN
IJALI(IJ)=0
IJX(IJ)=1
END IF
END DO
END IF
RETURN
END IF
C
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
IF(FRCMAX.LE.0.) THEN
FREQCO(1)=8.D11*TEFF
ELSE
FREQCO(1)=FRCMAX
END IF
IL0=1
IF(FREQCO(1).LT.CFRMAX*FRLEV(IL0) .AND. CFRMAX.GT.UN) THEN
FREQCO(1)=CFRMAX*FRLEV(IL0)
GO TO 20
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
20 CONTINUE
C
if(nftail.le.0) then
if(nftail.eq.-1.or.dftail.eq.0.) then
NEND1=NFREQC
FREQCO(NEND1)=FRCMIN
IJXCO(1)=1
IJXCO(NEND1)=1
XEND=UN/FLOAT(NEND1-1)
D121=(FREQCO(1)/FREQCO(NEND1))**XEND
DO IJ=2,NEND1-1
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
else
c
NEND1=-NFTAIL
FREQCO(NEND1)=DFTAIL*FREQCO(1)
IJXCO(1)=1
IJXCO(NEND1)=1
XEND=UN/FLOAT(NEND1-1)
D121=(FREQCO(1)/FREQCO(NEND1))**XEND
DO IJ=2,NEND1-1
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
NEND2=NFREQC
IJXCO(NEND1)=1
IJXCO(NEND1+1)=1
XEND=UN/FLOAT(NEND2-NEND1)
d121=(freqco(nend1)/frcmin)**xend
DO IJ=NEND1+1,NEND2
freqco(ij)=freqco(ij-1)/d121
IJXCO(IJ)=2
END DO
D121=THIRD*(FREQCO(NEND1)-FREQCO(NEND1+1))
DO IJ=NEND1+1,NEND2-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
else
C
NEND=NFTAIL
DIVEND=DFTAIL
NJC=NFREQC/5
DNX=UN-UN/FLOAT(NJC)
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
write(*,*) 'ils,ijfl',ils,ijfl(ils),freqco(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 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
end if
C
200 SUMWC=0.
C
IF(NFREAD.GT.0) THEN
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)
WCH(IJ)=WCHCO(IJ)
IJALI(IJ)=1
IJX(IJ)=IJXCO(IJ)
PROF(IJ)=0.
END DO
DO 320 ITR=1,NTRANS
IF(.NOT.LINE(ITR).OR.INDEXP(ITR).EQ.0) GO TO 320
IF(IFR0(ITR).GT.0) IFR0(ITR)=IFR0(ITR)+NFREQC
IF(IFR1(ITR).GT.0) IFR1(ITR)=IFR1(ITR)+NFREQC
320 CONTINUE
NFREQ=NFREQ+NFREQC
END IF
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 340 ITR=1,NTRANS
IF(LINE(ITR)) GO TO 340
IF(IFR0(ITR).LE.0) IFR0(ITR)=1
IF(IFR1(ITR).GT.0) GO TO 340
IF1=0
FR01=FR0(ITR)
MODE=INDEXP(ITR)
IF(IABS(MODE).EQ.5.OR.IABS(MODE).EQ.15) FR01=FR0PC(ITR)
DO 330 IJ=1,NFREQC
IF(FREQ(IJ).GE.FR01) IF1=IJ
330 CONTINUE
IFR1(ITR)=IF1
340 CONTINUE
if(nfrecl.gt.nfreq) nfrecl=nfreq
nfreqe=nfrecl
do ij=1,nfreqe
ijali(ij)=0
end do
c
RETURN
END
Reference in New Issue View Git Blame Copy Permalink