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

118 lines
3.5 KiB
Fortran

SUBROUTINE INIBL1(IGRD)
C =======================
C
C AUXILIARY INITIALIZATION PROCEDURE
C
INCLUDE 'PARAMS.FOR'
INCLUDE 'MODELP.FOR'
INCLUDE 'LINDAT.FOR'
INCLUDE 'SYNTHP.FOR'
INCLUDE 'WINCOM.FOR'
COMMON/LIMPAR/ALAM0,ALAM1,FRMIN,FRLAST,FRLI0,FRLIM
COMMON/BLAPAR/RELOP,SPACE0,CUTOF0,TSTD,DSTD,ALAMC
common/alsave/ALAM0s,ALASTs,CUTOF0s,CUTOFSs,RELOPs,SPACEs
common/plaopa/plalin,plcint,chcint
common/conabs/absoc(mfreqc),emisc(mfreqc),scatc(mfreqc),
* plac(mfreqc)
parameter (un=1.,bnc=1.4743e-2,hkc=4.79928e4,
* clc=2.997925e17)
DIMENSION CROSS(MCROSS,MFRQ),
* ABSO(MFREQ),EMIS(MFREQ),SCAT(MFREQ)
C
C auxiliary quantities for dissolved fractions
C
DO ID=1,ND
CALL DWNFR0(ID)
CALL WNSTOR(ID)
anh2(id)=0.
anhm(id)=0.
anch(id)=0.
anoh(id)=0.
END DO
CALL TINT
c
c reset wavelengths in case of opacity grid calculations
c
if(igrd.ge.0) then
alam0=alam0s
if(alam0s.eq.0.) alam0=5.e7/temp(1)/10.
if(alam0s.lt.0.) alam0=-5.e7/temp(1)/alam0s
alast=alasts
if(alasts.eq.0.) alast=5.e7/temp(1)*20.
if(alasts.lt.0.) alast=-5.e7/temp(1)*alasts
c if(alast.gt.1.e5) alast=1.e5
cutof0=cutof0s
cutofs=cutofss
relop=relops
if(relops.eq.0) then
relop=1.e-15
if(temp(1).lt.2.e6) relop=1.e-6
if(temp(1).lt.1.e6) relop=1.e-5
if(temp(1).lt.1.e5) relop=1.e-4
end if
space=spaces
ALAMC=(ALAM0+ALAST)*0.5
if(space.eq.0.) space=4.3e-8*sqrt(temp(idstd))*alamc
if(space.lt.0.) space=-5.72e-8*sqrt(temp(idstd))*alamc*space
SPACF=2.997925E18/ALAMC/ALAMC*SPACE
CUTOF0=0.1*CUTOF0
SPACE0=SPACE*0.1
ALAM0=1.D-1*ALAM0
ALAST=1.D-1*ALAST
ALAMC=ALAMC*0.1
ALST00=ALAST
FRLAST=CLC/ALAST
c
nfreqc=ifix(real(cutofs,4))
if(nfreqc.eq.0) nfreqc=mfreq
all0=log(alam0)
all1=log(alast)
dlc=(all1-all0)/(nfreqc-1)
xcc0=hkc/temp(1)
do ijc=1,nfreqc
wlamc(ijc)=exp(all0+(ijc-1)*dlc)
freqc(ijc)=clc/wlamc(ijc)
c frc=freqc(ijc)*1.e-15
c plac(ijc)=bnc*frc**3/(exp(xcc0*frc)-un)
end do
id=1
CALL CROSEW(CROSS)
CALL OPACON(ID,CROSS,ABSOC,EMISC,SCATC)
wc0=(freqc(1)-freqc(2))*0.5
wc1=(freqc(nfreqc-1)-freqc(nfreqc))*0.5
do ijc=2,nfreqc-1
absoc(ijc)=min(absoc(ijc),1.e30)
write(26,642) wlamc(ijc)*10.,log(absoc(ijc)/dens(1))
end do
642 format(f11.3,1p5e13.5)
c
do ijc=1,nfreqc
abstdw(ijc,id)=absoc(ijc)
end do
c
end if
c
c calculate the characteristic standard opacity
c
IF(IMODE.LE.2.and.imode.ge.-2) THEN
if(ifwin.le.0) then
CALL CROSET(CROSS)
DO ID=1,ND
CALL OPAC(ID,CROSS,ABSO,EMIS,SCAT)
ABSTD(ID)=MIN(ABSO(1)+SCAT(1),ABSO(2)+SCAT(2))
END DO
else
CALL CROSEW(CROSS)
DO ID=1,ND
CALL OPACW(ID,CROSS,ABSO,EMIS,ABSOC,EMISC,SCATC,0)
DO IJ=1,NFREQC
denscon(id)=1.
ABSTDW(IJ,ID)=ABSOC(IJ)/DENSCON(ID)
END DO
END DO
end if
END IF
C
RETURN
END