145 lines
3.8 KiB
Fortran
145 lines
3.8 KiB
Fortran
SUBROUTINE SABOLF(ID)
|
|
C =====================
|
|
C
|
|
C Saha-Boltzmann factors (SBF)
|
|
C and "upper sums" - sum of Saha-Boltzmann factors for upper, LTE,
|
|
C levels which are not included explicitly (USUM), and derivatives
|
|
C wrt. temperature (T) and electron density (DUSUMN)
|
|
C
|
|
C Input: ID - depth index
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
PARAMETER (UH=1.5)
|
|
PARAMETER (CMAX=2.154D4,CCON=2.0706D-16)
|
|
C
|
|
C DCHI - approximate lowering of ionization potential for neutrals
|
|
C Actual lowering is DCHI*effective charge, and is considered only
|
|
C if IUPSUM(ION).GT.0
|
|
C
|
|
c
|
|
if(ioptab.lt.0) return
|
|
c
|
|
T=TEMP(ID)
|
|
SQT=SQRT(T)
|
|
ANE=ELEC(ID)
|
|
STANE=SQRT(T/ANE)
|
|
XMAX=CMAX*SQRT(STANE)
|
|
TK=BOLK*T
|
|
CON=CCON/T/SQT
|
|
C
|
|
C Saha-Boltzmann factors
|
|
C
|
|
DO ION=1,NION
|
|
QZ=IZ(ION)
|
|
CFN=CON/G(NNEXT(ION))
|
|
DCH=0.
|
|
IUPS=IUPSUM(ION)
|
|
SSBF=0.
|
|
DSSBFT=0.
|
|
USUM(ION)=0.
|
|
DUSUMT(ION)=0.
|
|
DUSUMN(ION)=0.
|
|
nlst=nlast(ion)
|
|
if(ifwop(nlst).ge.0) then
|
|
nl1up=nquant(nlst)+1
|
|
else
|
|
nl1up=nquant(nlst)
|
|
end if
|
|
DO II=NFIRST(ION),NLAST(ION)
|
|
if(ifwop(ii).lt.0) then
|
|
E=EH*QZ*QZ/TK
|
|
SUM=0.
|
|
DO J=nl1up,NLMX
|
|
XJ=J
|
|
XI=J*J
|
|
X=E/XI
|
|
FI=XI*EXP(X)*WNHINT(J,ID)
|
|
SUM=SUM+FI
|
|
END DO
|
|
g(ii)=sum*two
|
|
gmer(imrg(ii),id)=g(ii)
|
|
end if
|
|
X=ENION(II)/TK
|
|
if(x.gt.110.) x=110.
|
|
SB=CFN*G(II)*EXP(X)
|
|
SBF(II)=SB
|
|
SSBF=SSBF+SB
|
|
DSBF(II)=-(UH+ENION(II)/TK)/T
|
|
DSSBFT=DSSBFT+SB*DSBF(II)
|
|
END DO
|
|
C
|
|
C Upper sums
|
|
C
|
|
if(ifwop(nlst).ge.0) then
|
|
IF(ION.EQ.IELHM) THEN
|
|
USUM(ION)=0.
|
|
DUSUMT(ION)=0.
|
|
DUSUMN(ION)=0.
|
|
GO TO 50
|
|
END IF
|
|
IF(IUPS.EQ.0) THEN
|
|
C
|
|
C 1. More exact approach - using (exact) partition functions
|
|
C
|
|
IAT=NUMAT(IATM(NFIRST(ION)))
|
|
XMX=XMAX*SQRT(QZ)
|
|
CALL PARTF(IAT,IZ(ION),T,ANE,XMX,U,DUT,DUN)
|
|
EE=ENION(NFIRST(ION))/TK
|
|
IF(EE.GT.110.) EE=110.
|
|
CFE=CFN*EXP(EE)
|
|
USUM(ION)=CFE*U-SSBF
|
|
DUSUMT(ION)=CFE*(DUT-U/T*(UH+EE))-DSSBFT
|
|
DUSUMN(ION)=CFE*DUN
|
|
xx=(ssbf-sbf(nfirst(ion)))/sbf(nfirst(ion))
|
|
IF(USUM(ION).LT.0.or.ee.ge.109.or.xx.lt.1.e-7) THEN
|
|
USUM(ION)=0.
|
|
DUSUMT(ION)=0.
|
|
DUSUMN(ION)=0.
|
|
END IF
|
|
C
|
|
C 2. Approximate approach - summation over fixed number of upper
|
|
C levels, assumed hydrogenic (ie. their ionization energy and
|
|
C statistical weight hydrogenic)
|
|
C
|
|
ELSE IF(IUPS.GT.0) THEN
|
|
SUM=0.
|
|
DSUM=0.
|
|
E=EH*QZ*QZ/TK
|
|
DO J=NQUANT(NLAST(ION))+1,IUPS
|
|
XI=J*J
|
|
X=E/XI
|
|
FI=XI*EXP(X)
|
|
SUM=SUM+FI
|
|
DSUM=DSUM-FI*(UH+X)/T
|
|
END DO
|
|
USUM(ION)=SUM*CON*TWO
|
|
DUSUMT(ION)=DSUM*CON*TWO
|
|
DUSUMN(ION)=0.
|
|
C
|
|
c 3. occupation probability form
|
|
c
|
|
ELSE
|
|
SUM=0.
|
|
DSUM=0.
|
|
E=EH*QZ*QZ/TK
|
|
DO J=NQUANT(NLAST(ION))+1,NLMX
|
|
XJ=J
|
|
XI=J*J
|
|
X=E/XI
|
|
FI=XI*EXP(X)*WNHINT(J,ID)
|
|
SUM=SUM+FI
|
|
DSUM=DSUM-FI*(UH+X)/T
|
|
END DO
|
|
USUM(ION)=SUM*CON*TWO
|
|
DUSUMT(ION)=DSUM*CON*TWO
|
|
DUSUMN(ION)=0.
|
|
END IF
|
|
end if
|
|
50 CONTINUE
|
|
END DO
|
|
RETURN
|
|
END
|