94 lines
2.4 KiB
Fortran
94 lines
2.4 KiB
Fortran
SUBROUTINE ELCOR(ID)
|
|
C ====================
|
|
C
|
|
C Procedure for a reevaluation of the electron number density
|
|
C from the charge conservation equation in the formal solution
|
|
C step (RESOLV)
|
|
C This procedure is called only if LCHC=false, ie. if the charge
|
|
C conservation equation is not part of the rate matrix
|
|
C
|
|
C Input: ID - depth index
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
COMMON/ADCHAR/QADD(MDEPTH)
|
|
c
|
|
if(ioptab.lt.0.or.ioptab.gt.0) return
|
|
C
|
|
T=TEMP(ID)
|
|
ANE=ELEC(ID)
|
|
AN0=DENS(ID)/WMM(ID)+ANE
|
|
C
|
|
C basic iteration loop for solving simultaneously a non-linear set
|
|
C of statistical equilibrium, charge conservation, and particle
|
|
C conservation equations
|
|
C
|
|
KKK=0
|
|
1 KKK=KKK+1
|
|
IF(IFIXDE.GT.0) THEN
|
|
AN=DENS(ID)/WMM(ID)+ANE
|
|
ELSE
|
|
AN=AN0
|
|
DENS(ID)=(AN-ANE)*WMM(ID)
|
|
END IF
|
|
C
|
|
C determine QQ, the total charge due to non-explicit atoms
|
|
C
|
|
QQ=0.
|
|
ANMNE1=WMM(ID)/DENS(ID)
|
|
if(ifmol.eq.0.or.t.ge.tmolim) then
|
|
CALL STATE(2,ID,T,ANE)
|
|
QQ=Q*ABUND(IATREF,ID)/YTOT(ID)*DENS(ID)/WMM(ID)
|
|
if(ioptab.gt.0) QQ=DENS(ID)/YTOT(ID)/WMM(ID)
|
|
else
|
|
aein=ane
|
|
call moleq(id,t,an,aein,ane,enrg,entt,wm,0)
|
|
qq=qadd(id)
|
|
end if
|
|
C
|
|
RHS=QFIX(ID)+QQ
|
|
DO IAT=1,NATOM
|
|
IF(IIFIX(IAT).NE.1) THEN
|
|
DO I=N0A(IAT),NKA(IAT)
|
|
IL=ILK(I)
|
|
CH=IZ(IEL(I))-1
|
|
IF(IL.GT.0) CH=IZ(IL)+(IZ(IL)-1)*USUM(IL)*ANE
|
|
IF(IMODL(I).GE.0) RHS=RHS+CH*POPUL(I,ID)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
C new electron density
|
|
C
|
|
RHS=HALF*(ANE+RHS)
|
|
ELEC(ID)=RHS
|
|
IF(IFIXDE.EQ.0) DENS(ID)=WMM(ID)*(AN-ELEC(ID))
|
|
ANMA(ID)=DENS(ID)/WMM(ID)
|
|
ANTO(ID)=ANMA(ID)+ELEC(ID)
|
|
RELANE=(RHS-ANE)/ANE
|
|
ANE=RHS
|
|
C
|
|
C second part of the iteration loop - recalculation of all
|
|
C populations with new electron density
|
|
C
|
|
CALL STEQEQ(ID,POP,1)
|
|
C
|
|
C convergence criterion for electron density
|
|
C
|
|
IF(ABS(RELANE).LE.1.D-6) THEN
|
|
CALL WNSTOR(ID)
|
|
RETURN
|
|
ENDIF
|
|
C
|
|
C if convergence is not achieved
|
|
C
|
|
IF(KKK.LT.10) GO TO 1
|
|
WRITE(6,601) ID,RELANE
|
|
WRITE(10,601) ID,RELANE
|
|
601 FORMAT('0 SLOW CONVERGENCE OF ELCOR ID =',I4,' REL =',1PD10.3/)
|
|
CALL WNSTOR(ID)
|
|
RETURN
|
|
END
|