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

277 lines
8.4 KiB
Fortran

SUBROUTINE REFLEV(ID,IMODE)
C ===========================
C
C determination of the global reference level, and
C determination of the LTE reference levels and corresponding
C quantities
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
INCLUDE 'ITERAT.FOR'
DIMENSION XSBF(MLEVEL)
if(ioptab.lt.0) return
C
C determination of the reference level (if required)
C
LREFP=.TRUE.
LREFP=.false.
IF(MODREF.GE.1) THEN
IF(IMODE.EQ.1.AND.(ITER.LE.1.OR.KANT(ITER).EQ.0)) THEN
DO IAT=1,NATOM
PMAX=0.
IREF=0
DO I=N0A(IAT),NKA(IAT)
IF(POPUL(I,ID).GE.PMAX) THEN
PMAX=POPUL(I,ID)
IREF=I
END IF
END DO
IF(IREF.NE.NKA(IAT).AND.IREF.NE.NFIRST(IEL(IREF)))
* IREF=NFIRST(IEL(IREF))
IF(MODREF.EQ.2.AND.IREF.LT.NFIRST(IEL(NKA(IAT))))
* IREF=NFIRST(IEL(NKA(IAT)))
NREF(IAT)=IREF
LREFP=LREFP .AND. NREFS(IAT,ID).EQ.NREF(IAT)
NREFS(IAT,ID)=NREF(IAT)
END DO
ELSE
DO IAT=1,NATOM
NREF(IAT)=NREFS(IAT,ID)
END DO
END IF
ELSE
DO IAT=1,NATOM
LREFP=LREFP .AND. NREFS(IAT,ID).EQ.NREF(IAT)
NREFS(IAT,ID)=NREF(IAT)
END DO
END IF
C
C set to the zeroing mode levels with SBF*ANE large;
C check whether reference level is not zerod, and is so, re-establish
C the reference level.
C
IF(ITER.LE.NITZER) THEN
DO IAT=1,NATOM
N1=N0A(IAT)
NK=NKA(IAT)
ISBMX=0
XSBMX=0.
IREF=NREFS(IAT,ID)
DO II=N1,NK
IF(ITER.LE.NITZER) IPZERO(II,ID)=0
XSBF(II)=SBF(II)*ELEC(ID)
IF(.NOT.LTE) THEN
ITR=ITRA(II,NNEXT(IEL(II)))
IF(ITR.GT.0) THEN
IF(RRU(ITR,ID).GT.1.E-30)
* XSBF(II)=SBF(II)*ELEC(ID)*RRD(ITR,ID)/RRU(ITR,ID)
END IF
END IF
END DO
IF(LTE) THEN
IREF=N1
DO II=N1+1,NK
IF(ILK(II).GT.0) THEN
IREF=II
IF(XSBF(II).GT.1.) GO TO 10
END IF
END DO
10 CONTINUE
NREFS(IAT,ID)=IREF
NREF(IAT)=IREF
END IF
IF(IREF.GT.N0A(IAT)) THEN
X=1.
DO II=IREF-1,N1,-1
IF(ILK(II).GT.0.OR.II.EQ.N0A(IAT)) THEN
X=X*XSBF(II)
IF(X.LT.POPZR2) THEN
DO III=N1,NLAST(IEL(II))
IPZERO(III,ID)=1
END DO
GO TO 20
END IF
END IF
END DO
END IF
20 CONTINUE
IF(IREF.LT.NK) THEN
X=1.
DO II=IREF+1,NK
IF(ILK(II).GT.0) THEN
X=X*XSBF(NFIRST(ILK(II)))
IF(X.GT.1./POPZR2) THEN
NFIR=NFIRST(IEL(II))
IF(II.EQ.NK) NFIR=NK
DO III=NFIR,NK
IPZERO(III,ID)=1
END DO
GO TO 30
END IF
END IF
END DO
END IF
30 CONTINUE
END DO
END IF
C
C determination of the LTE reference levels and corresponding
C quantities
C
ELEC1(ID)=UN/ELEC(ID)
DO IAT=1,NATOM
IFSUP=0
IREF=NREF(IAT)
DO I=N0A(IAT),NKA(IAT)
C
C generalized LTE reference level formalism
C
IF(IABS(IMODL(I)).EQ.1.OR.IABS(IMODL(I)).EQ.2) THEN
IN=NNEXT(IEL(I))
IF(I.LT.IREF.OR.POPUL(I,ID).LT.POPUL(IN,ID)) THEN
ILTREF(I,ID)=IN
SBPSI(I,ID)=SBF(I)*ELEC(ID)*WOP(I,ID)*BFAC(I,ID)
DSBPST(I,ID)=DSBF(I)
DSBPSN(I,ID)=ELEC1(ID)
ELSE IF(I.GT.IREF) THEN
I1=NFIRST(IEL(I))
ILTREF(I,ID)=I1
SBPSI(I,ID)=SBF(I)*WOP(I,ID)*BFAC(I,ID)/
* (SBF(I1)*WOP(I1,ID)*BFAC(I1,ID))
DSBPST(I,ID)=DSBF(I)-DSBF(I1)
DSBPSN(I,ID)=0.
END IF
C
C original scheme
C
ELSE IF(IABS(IMODL(I)).EQ.3.OR.IMODL(I).EQ.0) THEN
ILTREF(I,ID)=NNEXT(IEL(I))
SBPSI(I,ID)=SBF(I)*ELEC(ID)*WOP(I,ID)*BFAC(I,ID)
DSBPST(I,ID)=DSBF(I)
DSBPSN(I,ID)=ELEC1(ID)
ELSE IF(IABS(IMODL(I)).EQ.5) THEN
IFSUP=1
END IF
END DO
C
C super-reference level formalism
C
LRF=LTE.AND.IREF.NE.NKA(IAT)
IF(IFSUP.EQ.1.OR.LRF) THEN
XA=UN
XT=0.
XN=0.
ILTREF(IREF,ID)=IREF
IF(IREF.GT.N0A(IAT)) THEN
DO I=IREF-1,N0A(IAT),-1
ILTREF(I,ID)=IREF
SBPSI(I,ID)=XA*SBF(I)*ELEC(ID)*BFAC(I,ID)*WOP(I,ID)
if(sbpsi(i,id).lt.popzr2) then
sbpsi(i,id)=0.
else
DSBPST(I,ID)=XT+DSBF(I)
DSBPSN(I,ID)=XN+ELEC1(ID)
end if
IF(I.EQ.NFIRST(IEL(I))) THEN
XA=SBPSI(I,ID)
XT=DSBPST(I,ID)
XN=DSBPSN(I,ID)
if(xn.eq.0.) ilk(i)=0
END IF
END DO
END IF
C
XA=UN
XT=0.
XN=0.
IF(IREF.LT.NKA(IAT)) THEN
IF(IREF.EQ.NLAST(IEL(IREF))) THEN
XA=UN/(SBF(IREF)*BFAC(IREF,ID)*WOP(IREF,ID)*ELEC(ID))
XT=-DSBF(IREF)
XN=-ELEC1(ID)
END IF
DO I=IREF+1,NKA(IAT)-1
ILTREF(I,ID)=IREF
I1=NFIRST(IEL(I))
SBB1=UN/(SBF(I1)*BFAC(I1,ID)*WOP(I1,ID))
SBPSI(I,ID)=XA*SBF(I)*BFAC(I,ID)*WOP(I,ID)*SBB1
if(sbpsi(i,id).lt.popzr2) then
sbpsi(i,id)=0.
else
DSBPST(I,ID)=XT+DSBF(I)-DSBF(I1)
DSBPSN(I,ID)=XN
end if
IF(I.EQ.NLAST(IEL(I))) THEN
XA=XA*SBB1*ELEC1(ID)
XT=XT-DSBF(I1)
XN=XN-ELEC1(ID)
END IF
END DO
I=NKA(IAT)
ILTREF(I,ID)=IREF
SBPSI(I,ID)=XA
DSBPST(I,ID)=XT
DSBPSN(I,ID)=XN
END IF
END IF
C
C different meaning of SBPSI for IMODL=2, i.e. for levels
C with fixed ratio of population to that of reference level
C
DO I=N0A(IAT),NKA(IAT)
IF(IABS(IMODL(I)).EQ.2) THEN
IF(POPUL(I,ID).EQ.0..OR.POPUL(ILTREF(I,ID),ID).EQ.0.)
* THEN
IN=NNEXT(IEL(I))
IF(I.LT.IREF.OR.POPUL(I,ID).LT.POPUL(IN,ID)) THEN
ILTREF(I,ID)=IN
SBPSI(I,ID)=SBF(I)*ELEC(ID)*WOP(I,ID)*BFAC(I,ID)
ELSE IF(I.GT.IREF) THEN
I1=NFIRST(IEL(I))
ILTREF(I,ID)=I1
SBPSI(I,ID)=SBF(I)*WOP(I,ID)*BFAC(I,ID)/
* (SBF(I1)*WOP(I1,ID)*BFAC(I1,ID))
END IF
ELSE
SBPSI(I,ID)=POPUL(I,ID)/POPUL(ILTREF(I,ID),ID)
END IF
DSBPST(I,ID)=0.
DSBPSN(I,ID)=0.
END IF
END DO
C
C different meaning of SBPSI for IMODL=6 or 7, i.e. for levels
C with fixed ratio of population to that of a "guiding" level
C
if(imode.eq.1) then
DO I=N0A(IAT),NKA(IAT)
IF(IABS(IMODL(I)).EQ.6) THEN
ILTREF(I,ID)=IGUIDE(I)
if(POPUL(ILTREF(I,ID),ID).gt.0.) then
SBPSI(I,ID)=POPUL(I,ID)/POPUL(ILTREF(I,ID),ID)
end if
DSBPST(I,ID)=0.
DSBPSN(I,ID)=0.
END IF
END DO
end if
END DO
C
C true LTE reference level
C
DO IAT=1,NATOM
DO I=N0A(IAT),NKA(IAT)
IF(IABS(IMODL(I)).LT.6) THEN
ILTERF(I,ID)=ILTREF(I,ID)
SBLPSI(I,ID)=SBPSI(I,ID)
ELSE
ILTERF(I,ID)=NNEXT(IEL(I))
SBLPSI(I,ID)=SBF(I)*ELEC(ID)*WOP(I,ID)*BFAC(I,ID)
END IF
END DO
END DO
RETURN
END