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

65 lines
1.6 KiB
Fortran

SUBROUTINE WNSTOR(ID)
C =====================
C
C Stores occupation probabilities for hydrogen levels
C in common WNCOM for further use
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
PARAMETER (SIXTH=UN/6.,CCOR=0.09)
parameter (p1=0.1402,p2=0.1285,p3=un,p4=3.15,p5=4.)
parameter (tkn=3.01,ckn=5.33333333,cb0=8.59d14,f23=-2./3.)
if(ioptab.lt.0) return
C
cb=cb0*bergfc
ANE=ELEC(ID)
A=CCOR*EXP(SIXTH*LOG(ANE))/SQRT(TEMP(ID))
z=un
x=exp(p4*log(un+p3*a))
c1=p1*(x+p5*(z-un)*a*a*a)
c2=p2*x
beta0=cb*z*z*z*exp(f23*log(ane))
DO I=1,NLMX
XN=I
if(xn.le.tkn) then
xkn=un
else
xn1=un/(xn+un)
xkn=ckn*xn*xn1*xn1
end if
beta=beta0*xkn*xi2(i)*xi2(i)
f=(c1*beta*beta*beta)/(un+c2*beta*sqrt(beta))
WNHINT(I,ID)=f/(un+f)
END DO
C
C array WOP - occupation probabilities for explicit levels
C (if ifwop>1, occ. probabilities have been initialized
C for iron-peak elements and are not updated)
C
do ii=1,nlevel
if(ifwop(ii).le.0) then
wop(ii,id)=un
else if(ifwop(ii).eq.1) then
ie=iel(ii)
nq=nquant(ii)
if(iz(ie).eq.1) then
wop(ii,id)=wnhint(nq,id)
else
z=iz(ie)
xn=nq
wop(ii,id)=wn(xn,a,ane,z)
end if
end if
if(ifwop(ii).gt.1.and.lte) wop(ii,id)=un
end do
RETURN
END