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