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

116 lines
3.4 KiB
Fortran

SUBROUTINE LYMLIN(IJ)
C =====================
C
C opacity and emissibvity in first 30 Lyman lines
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
C
PARAMETER (SIXTH=1./6.,TTW=2./3.,
* OS0=0.02654,CPP=4.1412E-16,CPJ=157803.,
* C00=1.25E-9,C18=2.997925E18)
PARAMETER (MLEVL=30)
C
DIMENSION ABLYM(MDEPTH),EMLYM(MDEPTH),
* XKIJL(MLEVL),FIJL(MLEVL),WL0L(MLEVL),FR0L(MLEVL),
* F00(MDEPTH),DOP0(MDEPTH),PJ(MLEVL,MDEPTH),
* ABTR(MLEVL,MDEPTH),EMTR(MLEVL,MDEPTH),FID(MLEVL,MDEPTH),
* AD0(MLEVL,MDEPTH),DIV0(MLEVL,MDEPTH),
* DBET0(MLEVL,MDEPTH),BETAD0(MLEVL,MDEPTH)
C
DATA INIP/1/
C
IF(IELH.EQ.0) RETURN
C
FR=FREQ(IJ)
IF(FR.GT.3.28805E15.OR.FR.LT.1.5E15) RETURN
C
IF(INIP.EQ.1) THEN
N0H=NFIRST(IELH)
N1H=NLAST(IELH)
NKH=NNEXT(IELH)
NLH=N1H-N0H
XII=1.
DO J=2,30
CALL STARK0(1,J,1,XKIJ0,WL00,FIJ0)
XKIJL(J)=XKIJ0
FIJL(J)=FIJ0
WL0L(J)=WL00
FR0L(J)=C18/WL0L(J)
END DO
DO ID=1,ND
T=TEMP(ID)
T1=1./T
SQT=SQRT(T)
ANE=ELEC(ID)
ANP=POPUL(NKH,ID)
F00(ID)=C00*ANE**TTW
DOP0(ID)=1.E8*SQRT(1.65E8*T+VTURB(ID))
C
P0=CPP*ANE*ANP*T1/SQT
P1=POPUL(N0H,ID)
DO J=2,30
X=J*J
JJ=N0H+J-1
XJJ=1./JJ
IF(J.LE.NLH) THEN
PJ(J,ID)=POPUL(JJ,ID)
ELSE
PJ(J,ID)=P0*EXP(CPJ/(X*T))*X*WNHINT(J,ID)
END IF
ABTR(J,ID)=P1*WNHINT(J,ID)
EMTR(J,ID)=PJ(J,ID)*XJJ*EXP(CPJ*(XII-XJJ)*T1)
C
FXK=F00(ID)*XKIJL(J)
DBETA=WL0L(J)*WL0L(J)/(C18*FXK)
FID(J,ID)=OS0*FIJL(J)*DBETA
DOP=DOP0(ID)/WL0L(J)
BETAD=DOP*DBETA
CALL DIVSTR(1)
AD0(J,ID)=ADH
DIV0(J,ID)=DIVH
DBET0(J,ID)=DBETA
BETAD0(J,ID)=BETAD
END DO
END DO
INIP=0
END IF
C
FR=FREQ(IJ)
WL=C18/FR
F15=FR*1.E-15
DO ID=1,ND
ABLYM(ID)=0.
EMLYM(ID)=0.
DO J=2,30
BETA=DBET0(J,ID)*ABS(FR-FR0L(J))
BETAD=BETAD0(J,ID)
ADH=AD0(J,ID)
DIVH=DIV0(J,ID)
SG=STARKA(BETA,TWO)*FID(J,ID)
ABLYM(ID)=ABLYM(ID)+SG*ABTR(J,ID)
EMLYM(ID)=EMLYM(ID)+SG*EMTR(J,ID)
c if(wl.gt.1120.0.and.wl.lt.1120.3.and.id.eq.50)
c * write(6,600) j,fr,fr0l(j),abs(fr-fr0l(j)),
c * dbet0(j,id),beta,fid(j,id),sg,abtr(j,id),
c * emtr(j,id),sg*abtr(j,id),ablym(id)
c if(wl.gt.1120.0.and.wl.lt.1120.3.and.id.eq.50)
c * write(6,600) j,fid(j,id),dbet0(j,id),wl0l(j),xkijl(j),
c * dop0(id),fijl(j)
c 600 format('lymsig',i4,1p11e11.3)
END DO
XKT=EXP(-4.79928e-11*FR/TEMP(ID))
XKB=XKT*1.4743E-2*F15*F15*F15
ABLYM(ID)=ABLYM(ID)-XKT*EMLYM(ID)
EMLYM(ID)=XKB*EMLYM(ID)
ABSO1(ID)=ABSO1(ID)+ABLYM(ID)
EMIS1(ID)=EMIS1(ID)+EMLYM(ID)
c if(wl.gt.1120.0.and.wl.lt.1120.3)
c * write(6,601) ij,wl,ablym(50),emlym(50),xkt,xkb
c 601 format('lymlin',i6,f10.3,1p4e11.3)
END DO
RETURN
END