379 lines
10 KiB
Fortran
379 lines
10 KiB
Fortran
SUBROUTINE HESOL6
|
|
C =================
|
|
C
|
|
C Solution of the coupled system of
|
|
C 1) hydrostatic equilibrium equation;
|
|
C 2) definitions of Ptotal, Pgas and rho;
|
|
C 3) state equation
|
|
C 4) z-m relation
|
|
C with a given temperature and radiation pressure
|
|
C Numerical solution by a Newton-Raphson method
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
PARAMETER (MP=6,NP=6,IP=1,IG=2,IR=3,IN=4,IE=5,IZ=6)
|
|
PARAMETER (NITERH=15)
|
|
COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
|
|
DIMENSION P(MDEPTH),ANTT(MDEPTH),CHNG(MDEPTH),
|
|
* VEC(MP,MDEPTH),VEC1(MP,MDEPTH),
|
|
* VEC2(MP,MDEPTH),VEC3(MP,MDEPTH),
|
|
* SOL(MP),ERR(MP),CHN(MP),
|
|
* A(MP,MP),B(MP,MP),C(MP,MP),VL(MP),
|
|
* D(MP,MP,MDEPTH),ANU(MP,MDEPTH),
|
|
* anerl(mdepth)
|
|
C
|
|
DATA ERROR /1.D-4/
|
|
C
|
|
C radiation pressure scale feight at the surface
|
|
C
|
|
ID=1
|
|
IF(ITER.EQ.0) THEN
|
|
HR1=SIG4P*TEFF**4*PCK*ABROSD(ID)/QGRAV
|
|
ELSE
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
IF(.NOT.LSKIP(ID,IJT)) THEN
|
|
FLUXW=W(IJT)*(FH(IJT)*RADEX(IJ,ID)-HEXTRD(IJT))
|
|
GRD=GRD+FLUXW*ABSOE1(IJ)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
HR1=PCK/QGRAV*(GRD+FPRD(ID))/DENS(ID)
|
|
END IF
|
|
do id=1,nd
|
|
ANTT(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
PGS(ID)=ANTT(ID)*BOLK*TEMP(ID)
|
|
PTOTAL(ID)=PRADT(ID)+PGS(ID)
|
|
end do
|
|
C
|
|
C Basic Newton-Raphson iteration loop
|
|
C
|
|
ITERH=0
|
|
LAC2H=.FALSE.
|
|
IACH=6
|
|
IACDH=4
|
|
IACH0=IACH-3
|
|
10 CONTINUE
|
|
ITERH=ITERH+1
|
|
C
|
|
C -------------------
|
|
C Forward elimination
|
|
C -------------------
|
|
C
|
|
DO ID=1,ND
|
|
ANTT(ID)=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
PGS(ID)=ANTT(ID)*BOLK*TEMP(ID)
|
|
PTOTAL(ID)=PRADT(ID)+PGS(ID)
|
|
P(ID)=PTOTAL(ID)
|
|
VEC(IP,ID)=PTOTAL(ID)
|
|
VEC(IG,ID)=PGS(ID)
|
|
VEC(IR,ID)=DENS(ID)
|
|
VEC(IN,ID)=ANTT(ID)
|
|
VEC(IE,ID)=ELEC(ID)
|
|
VEC(IZ,ID)=ZD(ID)
|
|
C
|
|
DO I=1,NP
|
|
VL(I)=0.
|
|
DO J=1,NP
|
|
B(I,J)=0.
|
|
A(I,J)=0.
|
|
C(I,J)=0.
|
|
END DO
|
|
END DO
|
|
C
|
|
IF(ID.EQ.1) THEN
|
|
HG1=SQRT(TWO*PGS(ID)/DENS(ID)/QGRAV)
|
|
X=(ZD(ID)-HR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
X1=X*1.01
|
|
F1D=0.
|
|
IF(X1.LT.3.) THEN
|
|
F1D=8.86226925D-1*EXP(X1*X1)*ERFCX(X1)
|
|
ELSE
|
|
F1D=HALF*(UN-HALF/X1/X1)/X1
|
|
END IF
|
|
IF(X.GT.0.) F1D=(F1D-F1)*100./X
|
|
B(IG,IG)=DENS(ID)*HG1*HALF/PGS(ID)*(F1-F1D*X)
|
|
B(IG,IR)=HG1*HALF*(F1+F1D*X)
|
|
B(IG,IZ)=DENS(ID)*F1D
|
|
VL(IG)=DM(ID)-DENS(ID)*HG1*F1
|
|
B(IP,IP)=UN
|
|
B(IP,IG)=-UN
|
|
VL(IP)=PRADT(ID)+PGS(ID)-P(ID)
|
|
ELSE
|
|
DMM=UN/(DM(ID)-DM(ID-1))
|
|
QG=HALF*QGRAV
|
|
B(IP,IP)=DMM
|
|
B(IP,IZ)=-QG
|
|
A(IP,IP)=DMM
|
|
A(IP,IZ)=QG
|
|
VL(IP)=QG*(ZD(ID)+ZD(ID-1))-(P(ID)-P(ID-1))*DMM
|
|
B(IG,IP)=UN
|
|
B(IG,IG)=-UN
|
|
VL(IG)=PRADT(ID)+PGS(ID)-P(ID)
|
|
END IF
|
|
C
|
|
B(IR,IR)=UN
|
|
B(IR,IN)=-WMM(ID)
|
|
B(IR,IE)=WMM(ID)
|
|
VL(IR)=wmm(id)*(antt(id)-elec(id))-dens(id)
|
|
C
|
|
B(IN,IG)=UN
|
|
B(IN,IN)=-BOLK*TEMP(ID)
|
|
VL(IN)=antt(id)*bolk*temp(id)-pgs(id)
|
|
C
|
|
T=TEMP(ID)
|
|
AN=ANTT(ID)
|
|
ANERL(ID)=ELEC(ID)/AN
|
|
B(IE,IE)=UN
|
|
b(ie,in)=-anerl(id)
|
|
vl(ie)=anerl(id)*antt(id)-elec(id)
|
|
C
|
|
IF(ID.LT.ND) THEN
|
|
DMP=(DM(ID+1)-DM(ID))*HALF
|
|
B(IZ,IR)=DMP/DENS(ID)**2
|
|
B(IZ,IZ)=UN
|
|
C(IZ,IR)=-DMP/DENS(ID+1)**2
|
|
C(IZ,IZ)=UN
|
|
VL(IZ)=ZD(ID+1)-ZD(ID)+DMP/DENS(ID)+DMP/DENS(ID+1)
|
|
ELSE
|
|
B(IZ,IZ)=UN
|
|
VL(IZ)=0.
|
|
END IF
|
|
C
|
|
IF(ID.GT.1) THEN
|
|
B(IP,IR)=B(IP,IR)-
|
|
* A(IP,IP)*D(IP,IR,ID-1)-
|
|
* A(IP,IZ)*D(IZ,IR,ID-1)
|
|
B(IP,IZ)=B(IP,IZ)-
|
|
* A(IP,IP)*D(IP,IZ,ID-1)-
|
|
* A(IP,IZ)*D(IZ,IZ,ID-1)
|
|
VL(IP)=VL(IP)+A(IP,IP)*ANU(IP,ID-1)+
|
|
* A(IP,IZ)*ANU(IZ,ID-1)
|
|
END IF
|
|
C
|
|
CALL MATINV(B,NP,MP)
|
|
C
|
|
DO I=1,NP
|
|
SUM=0.
|
|
DO J=1,NP
|
|
SUM=SUM+B(I,J)*VL(J)
|
|
END DO
|
|
ANU(I,ID)=SUM
|
|
END DO
|
|
IF(ID.LT.ND) THEN
|
|
DO I=1,NP
|
|
D(I,IR,ID)=B(I,IZ)*C(IZ,IR)
|
|
D(I,IZ,ID)=B(I,IZ)*C(IZ,IZ)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
C -------------------
|
|
C backsubstitution
|
|
C -------------------
|
|
C
|
|
CHANM=0
|
|
DO ID=ND,1,-1
|
|
CHNG(ID)=0.
|
|
IF(ID.EQ.ND) THEN
|
|
DO I=1,NP
|
|
SOL(I)=ANU(I,ID)
|
|
END DO
|
|
ELSE
|
|
DO I=1,NP
|
|
ANU(I,ID)=ANU(I,ID)+D(I,IR,ID)*ANU(IR,ID+1)+
|
|
* D(I,IZ,ID)*ANU(IZ,ID+1)
|
|
SOL(I)=ANU(I,ID)
|
|
END DO
|
|
END IF
|
|
DO I=1,NP
|
|
CHAN=0.
|
|
IF(VEC(I,ID).NE.0.) CHAN=SOL(I)/VEC(I,ID)
|
|
CHN(I)=CHAN
|
|
IF(ABS(CHAN).GT.CHANM) CHANM=ABS(CHAN)
|
|
IF(ABS(CHAN).GT.CHNG(ID)) CHNG(ID)=ABS(CHAN)
|
|
IF(CHAN.LT.-0.99) CHAN=-0.99
|
|
IF(CHAN.GT.99.00) CHAN=99.00
|
|
VEC(I,ID)=VEC(I,ID)*(UN+CHAN)
|
|
END DO
|
|
C
|
|
PTOTAL(ID)=VEC(IP,ID)
|
|
P(ID)=PTOTAL(ID)
|
|
PGS(ID)=VEC(IG,ID)
|
|
DENS(ID)=VEC(IR,ID)
|
|
ANTT(ID)=VEC(IN,ID)
|
|
ELEC(ID)=VEC(IE,ID)
|
|
ZD(ID)=VEC(IZ,ID)
|
|
C
|
|
IF(ID.EQ.1) THEN
|
|
HG1=SQRT(TWO*PGS(ID)/DENS(ID)/QGRAV)
|
|
X=(ZD(ID)-HR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
ERR(IP)=(DENS(ID)*HG1*F1-DM(ID))/DM(ID)
|
|
ELSE IF(ID.LT.ND) THEN
|
|
ERR(IP)=(P(ID+1)-P(ID))/(DM(ID+1)-DM(ID))*
|
|
* TWO/QGRAV/(ZD(ID+1)+ZD(ID))-UN
|
|
ELSE
|
|
ERR(IP)=0.
|
|
END IF
|
|
IF(P(ID).NE.0.) ERR(IG)=
|
|
* (PRADT(ID)+PGS(ID)-P(ID))/P(ID)
|
|
IF(DENS(ID).NE.0.) ERR(IR)=
|
|
* (DENS(ID)-(ANTT(ID)-ELEC(ID))*WMM(ID))/DENS(ID)
|
|
IF(PGS(ID).NE.0.) ERR(IN)=
|
|
* (PGS(ID)-ANTT(ID)*BOLK*TEMP(ID))/PGS(ID)
|
|
IF(ELEC(ID).NE.0.)
|
|
* ERR(IE)=(ELEC(ID)-ANerl(id)*antt(id))/ELEC(ID)
|
|
IF(ID.LT.ND) THEN
|
|
ERR(IZ)=(ZD(ID)-ZD(ID+1))*TWO/(DM(ID+1)-DM(ID))/
|
|
* (UN/DENS(ID)+UN/DENS(ID+1))-UN
|
|
ELSE
|
|
ERR(IZ)=ZD(ID)
|
|
END IF
|
|
END DO
|
|
C
|
|
C acceleration
|
|
C
|
|
IF(NITERH.LT.IACH .OR. ITERH.LT.IACH0) GO TO 100
|
|
ipngH=1
|
|
if(iacdH.gt.0) ipngH=mod((iterH-iacH),iacdH)
|
|
if(.not.lac2H) then
|
|
IPT=MOD(ITERH,3)
|
|
IPT0=MOD(IACH,3)
|
|
IPT1=MOD((IACH+1),3)
|
|
IPT2=MOD((IACH+2),3)
|
|
IF(ITERH.EQ.IACH0) THEN
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC3(IX,ID)=VEC(IX,ID)
|
|
END DO
|
|
END DO
|
|
ELSE IF(IPT.EQ.IPT1) THEN
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC2(IX,ID)=VEC(IX,ID)
|
|
END DO
|
|
END DO
|
|
ELSE IF(IPT.EQ.IPT2) THEN
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC1(IX,ID)=VEC(IX,ID)
|
|
END DO
|
|
END DO
|
|
ENDIF
|
|
else if (ipngH.ne.0) then
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC3(IX,ID)=VEC2(IX,ID)
|
|
END DO
|
|
END DO
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC2(IX,ID)=VEC1(IX,ID)
|
|
END DO
|
|
END DO
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC1(IX,ID)=VEC(IX,ID)
|
|
END DO
|
|
END DO
|
|
GO TO 100
|
|
end if
|
|
|
|
IF(ITERH.LT.IACH) GO TO 100
|
|
C
|
|
A1=0.
|
|
B1=0.
|
|
B2=0.
|
|
C1=0.
|
|
C2=0.
|
|
DO IX=1,NP
|
|
DO ID=1,ND
|
|
WT=0.
|
|
IF(VEC(IX,ID).NE.0.) WT=1./ABS(VEC(IX,ID))
|
|
D0=VEC(IX,ID)-VEC1(IX,ID)
|
|
D1=D0-VEC1(IX,ID)+VEC2(IX,ID)
|
|
D2=D0-VEC2(IX,ID)+VEC3(IX,ID)
|
|
A1=A1+WT*D1*D1
|
|
B1=B1+WT*D1*D2
|
|
B2=B2+WT*D2*D2
|
|
C1=C1+WT*D0*D1
|
|
C2=C2+WT*D0*D2
|
|
END DO
|
|
END DO
|
|
AB=B2*A1-B1*B1
|
|
IF(AB.EQ.0.) THEN
|
|
IACH=IACH+IACDH
|
|
IACH0=IACH-3
|
|
GO TO 100
|
|
END IF
|
|
AA=(B2*C1-B1*C2)/AB
|
|
BB=(A1*C2-B1*C1)/AB
|
|
DO ID=1,ND
|
|
DO IX=1,NP
|
|
VEC(IX,ID)=(UN-AA-BB)*VEC(IX,ID)+AA*VEC1(IX,ID)+
|
|
* BB*VEC2(IX,ID)
|
|
END DO
|
|
END DO
|
|
LAC2H=.TRUE.
|
|
100 CONTINUE
|
|
C
|
|
DO ID=1,ND
|
|
PTOTAL(ID)=VEC(IP,ID)
|
|
P(ID)=PTOTAL(ID)
|
|
PGS(ID)=VEC(IG,ID)
|
|
DENS(ID)=VEC(IR,ID)
|
|
ANTT(ID)=VEC(IN,ID)
|
|
ELEC(ID)=VEC(IE,ID)
|
|
ZD(ID)=VEC(IZ,ID)
|
|
C
|
|
IF(ID.EQ.1) THEN
|
|
HG1=SQRT(TWO*PGS(ID)/DENS(ID)/QGRAV)
|
|
X=(ZD(ID)-HR1)/HG1
|
|
IF(X.LT.3.) THEN
|
|
IF(X.LT.0.) X=0.
|
|
F1=8.86226925D-1*EXP(X*X)*ERFCX(X)
|
|
ELSE
|
|
F1=HALF*(UN-HALF/X/X)/X
|
|
END IF
|
|
ERR(IP)=(DENS(ID)*HG1*F1-DM(ID))/DM(ID)
|
|
ELSE IF(ID.LT.ND) THEN
|
|
ERR(IP)=(P(ID)-P(ID-1))/(DM(ID)-DM(ID-1))*
|
|
* TWO/QGRAV/(ZD(ID)+ZD(ID-1))-UN
|
|
END IF
|
|
IF(P(ID).NE.0.) ERR(IG)=
|
|
* (PRADT(ID)+PGS(ID)-P(ID))/P(ID)
|
|
IF(DENS(ID).NE.0.) ERR(IR)=
|
|
* (DENS(ID)-(ANTT(ID)-ELEC(ID))*WMM(ID))/DENS(ID)
|
|
IF(PGS(ID).NE.0.) ERR(IN)=
|
|
* (PGS(ID)-ANTT(ID)*BOLK*TEMP(ID))/PGS(ID)
|
|
IF(ELEC(ID).NE.0.)
|
|
* ERR(IE)=(ELEC(ID)-ANerl(id)*antt(id))/ELEC(ID)
|
|
IF(ID.GT.1) THEN
|
|
ERR(IZ)=(ZD(ID-1)-ZD(ID))*TWO/(DM(ID)-DM(ID-1))/
|
|
* (UN/DENS(ID-1)+UN/DENS(ID))-UN
|
|
ELSE
|
|
ERR(IZ)=0.
|
|
END IF
|
|
END DO
|
|
C
|
|
C Convergence criterion for the Newton-Raphson method
|
|
C
|
|
IF(CHANM.GT.ERROR.AND.ITERH.LT.NITERH) GO TO 10
|
|
C
|
|
RETURN
|
|
END
|