207 lines
6.1 KiB
Fortran
207 lines
6.1 KiB
Fortran
SUBROUTINE RATMAT(ID,IICAL,IMODE,A,B)
|
|
C =====================================
|
|
C
|
|
C Evaluation of the rate matrix
|
|
C
|
|
C More specifically, the set of statistical equilibrium + constraint
|
|
C equations (see further) is represented as
|
|
C
|
|
C A * vector of populations = B
|
|
C
|
|
C (see Mihalas, 1978) pp.138-139 for details
|
|
C
|
|
C Constraint equations:
|
|
C the rows corresponding to the last explicit level of each atom
|
|
C (ie. a one-level ion representing the highest ionization degree
|
|
C of the given atom) are the abundance definition equations
|
|
C
|
|
C Input: ID - depth index
|
|
C Input transmitted by COMMON blocks:
|
|
C RRU - array of upward radiative rates in all transitions
|
|
C RRD - array of downward radiative rates in all transitions
|
|
C COLRAT - array of collisional rates (only upward) in all
|
|
C transitions
|
|
C SBF - array of Saha-Boltzmann factors
|
|
C USUM - array of upper sums
|
|
C
|
|
C Output: A - rate matrix
|
|
C B - the right-hand-side vector
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
DIMENSION A(MLEVEL,MLEVEL),B(MLEVEL)
|
|
DIMENSION AIJ(MTRANS),AJI(MTRANS)
|
|
DIMENSION SBW(MLEVEL),IICAL(MLEVEL)
|
|
DIMENSION LLTE(MLEVEL)
|
|
c DATA ICOMP /0/
|
|
c
|
|
if(ioptab.lt.0) return
|
|
c
|
|
if(ipslte.ne.0) then
|
|
do itr=1,ntrans
|
|
rru(itr,id)=0.
|
|
rrd(itr,id)=0.
|
|
drdt(itr,id)=0.
|
|
end do
|
|
end if
|
|
c
|
|
T=TEMP(ID)
|
|
ANE=ELEC(ID)
|
|
HKT=HK/T
|
|
TK=HKT/H
|
|
DO I=1,NLEVEL
|
|
B(I)=0.
|
|
SBW(I)=ANE*SBF(I)*WOP(I,ID)
|
|
ILT=ILTION(IEL(I))
|
|
LLT=ILT.EQ.1.AND.IMODE.EQ.0
|
|
LLTE(I)=LLT.OR.LTE.OR.ILTLEV(I).GE.1.OR.ILT.GE.2
|
|
LLTE(I)=LLTE(I).OR.ID.GE.IDLTE
|
|
DO J=1,NLEVEL
|
|
A(J,I)=0.
|
|
END DO
|
|
END DO
|
|
C
|
|
C determine reference levels
|
|
C
|
|
CALL REFLEV(ID,IABS(IMODE))
|
|
C
|
|
C ******* First part - the rows corresponding to true statistical
|
|
C equilibrium equations
|
|
C i.e. NLTE rate equations
|
|
C
|
|
C 1) simple expression in the case of LTE
|
|
C
|
|
DO IAT=1,NATOM
|
|
IF(IIFIX(IAT).NE.1.OR.IMODE.LT.0) THEN
|
|
NREFI=NREFS(IAT,ID)
|
|
DO I=N0A(IAT),NKA(IAT)
|
|
II=IABS(IICAL(I))
|
|
IF(I.NE.NREFI.AND.II.NE.0.AND.LLTE(I)) THEN
|
|
A(II,II)=UN
|
|
N=IABS(IICAL(ILTERF(I,ID)))
|
|
A(II,N)=A(II,N)-SBLPSI(I,ID)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
C 2) NLTE rate equation
|
|
C
|
|
IF(.NOT.LTE) THEN
|
|
DO 40 ITR=1,NTRANS
|
|
I=ILOW(ITR)
|
|
IF(IIFIX(IATM(I)).EQ.1) GO TO 40
|
|
J=IUP(ITR)
|
|
NKE=NNEXT(IEL(I))
|
|
C
|
|
C upward total rate
|
|
C
|
|
AIJ(ITR)=(COLRAT(ITR,ID)+RRU(ITR,ID))*WOP(J,ID)
|
|
C
|
|
C downward total rate
|
|
C
|
|
IF(LINE(ITR)) THEN
|
|
C
|
|
C bound-bound
|
|
C
|
|
AJI(ITR)=(COLTAR(ITR,ID)+RRD(ITR,ID)*
|
|
* G(I)/G(J)*EXP(HKT*FR0(ITR)))*WOP(I,ID)
|
|
ELSE
|
|
C
|
|
C Bound-free
|
|
C Quantity CORR is a factor which allows the bound-free
|
|
C transition to end at an excited state
|
|
C
|
|
CORR=UN
|
|
IF(NKE.NE.J) CORR=G(NKE)/G(J)*
|
|
* EXP((ENION(NKE)-ENION(J))*TK)
|
|
AJI(ITR)=COLTAR(ITR,ID)*WOP(I,ID)+
|
|
* RRD(ITR,ID)*SBW(I)*CORR
|
|
END IF
|
|
IF(IICAL(I).LT.0.AND.IICAL(J).LT.0) THEN
|
|
AIJ(ITR)=AIJ(ITR)*SBPSI(I,ID)
|
|
AJI(ITR)=AJI(ITR)*SBPSI(J,ID)
|
|
END IF
|
|
40 CONTINUE
|
|
C
|
|
C Elements of the rate matrix
|
|
C
|
|
DO 41 ITR=1,NTRANS
|
|
I=ILOW(ITR)
|
|
IF(IIFIX(IATM(I)).EQ.1) GO TO 41
|
|
NREFI=NREFS(IATM(I),ID)
|
|
J=IUP(ITR)
|
|
II=IABS(IICAL(I))
|
|
JJ=IABS(IICAL(J))
|
|
IF(IPZERO(I,ID).NE.0.OR.IPZERO(J,ID).NE.0) GO TO 41
|
|
IF(I.NE.NREFI.AND.II.GT.0.AND..NOT.LLTE(I)) THEN
|
|
A(II,II)=A(II,II)+AIJ(ITR)
|
|
IF(JJ.GT.0) THEN
|
|
A(II,JJ)=A(II,JJ)-AJI(ITR)
|
|
ELSE
|
|
JJJ=IICAL(ILTREF(J,ID))
|
|
A(II,JJJ)=A(II,JJJ)-AJI(ITR)*SBPSI(J,ID)
|
|
END IF
|
|
END IF
|
|
IF(J.NE.NREFI.AND.JJ.GT.0.AND..NOT.LLTE(J)) THEN
|
|
A(JJ,JJ)=A(JJ,JJ)+AJI(ITR)
|
|
IF(II.GT.0) THEN
|
|
A(JJ,II)=A(JJ,II)-AIJ(ITR)
|
|
ELSE
|
|
III=IICAL(ILTREF(I,ID))
|
|
A(JJ,III)=A(JJ,III)-AIJ(ITR)*SBPSI(I,ID)
|
|
END IF
|
|
END IF
|
|
41 CONTINUE
|
|
END IF
|
|
C
|
|
C reset the rate matrix elements for "small" populations
|
|
C
|
|
DO I=1,NLEVEL
|
|
II=IICAL(I)
|
|
IF(II.GT.0) THEN
|
|
IF(IPZERO(I,ID).GT.0) THEN
|
|
DO J=1,NLEVEL
|
|
A(II,J)=0.
|
|
END DO
|
|
A(II,II)=1.
|
|
END IF
|
|
ELSE IF(II.LT.0) THEN
|
|
IF(IGZERO(-II,ID).GT.0) THEN
|
|
DO J=1,NLEVEL
|
|
A(-II,J)=0.
|
|
END DO
|
|
A(-II,-II)=1.
|
|
END IF
|
|
END IF
|
|
END DO
|
|
C
|
|
C ******** Second part - abundance definition equations
|
|
C
|
|
DO 100 IAT=1,NATOM
|
|
IF(IIFIX(IAT).EQ.1.AND.IMODE.GE.0) GO TO 100
|
|
NREFII=IABS(IICAL(NREFS(IAT,ID)))
|
|
DO I=N0A(IAT),NKA(IAT)
|
|
IL=ILK(I)
|
|
II=IICAL(I)
|
|
IF(II.GT.0) THEN
|
|
A(NREFII,II)=A(NREFII,II)+UN
|
|
IF(IL.NE.0) A(NREFII,II)=A(NREFII,II)+ANE*USUM(IL)
|
|
ELSE IF(II.LT.0) THEN
|
|
A(NREFII,-II)=A(NREFII,-II)+SBPSI(I,ID)
|
|
IF(IL.NE.0) A(NREFII,-II)=A(NREFII,-II)+ANE*USUM(IL)*
|
|
* SBPSI(I,ID)
|
|
ELSE
|
|
III=IICAL(ILTREF(I,ID))
|
|
A(NREFII,III)=A(NREFII,III)+SBPSI(I,ID)
|
|
END IF
|
|
END DO
|
|
B(NREFII)=B(NREFII)+
|
|
* DENS(ID)/WMM(ID)/YTOT(ID)*ABUND(IAT,ID)
|
|
100 CONTINUE
|
|
C
|
|
RETURN
|
|
END
|