SUBROUTINE HESOLV
C     =================
C
C     Solution of the coupled system of hydrostatic equilibrium equation
C     and the z-m relation; with a given (generally depth-dependent) sound
C     speed, defined as total pressure/density.
C     Numerical solution by a Newton-Raphson method
C
C     Input:  P     - initial total pressure
C             VSND2 - sound speed squared
C             HG1   - gas pressure scale height at the surface
C             RR1   - ratio of radiation and gas pressure scale heights
C                     at the surface
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      COMMON/PRSAUX/VSND2(MDEPTH),HG1,HR1,RR1
      DIMENSION P(MDEPTH),B(2,2),C(2,2),VL(2),
     *          D(2,2,MDEPTH),ANU(2,MDEPTH)
C
      DATA ERROR /1.D-4/
C
C     Density for a given total pressure and sound speed
C
      DO ID=1,ND
         P(ID)=PTOTAL(ID)
         vsnd2(id)=p(id)/dens(id)
      END DO
C
C     Consistent z-values
C
c     ZD(ND)=ZND
c     DO IID=1,ND-1
c        ID=ND-IID
c        ZD(ID)=ZD(ID+1)+HALF*(DM(ID+1)-DM(ID))*(UN/DENS(ID)+
c    *                   UN/DENS(ID+1))
c     END DO
C
C     Basic Newton-Raphson iteration loop
C
      ITERH=0
   30 CONTINUE
      ITERH=ITERH+1
C
C     -------------------
C     Forward elimination
C     -------------------
C
C     Upper boundary condition
C
      ID=1
      X=ZD(1)/HG1-RR1
      IF(X.LT.3.) THEN
         IF(X.LT.0.) X=0.
         F1=1.772453851D0*EXP(X*X)*ERFCX(X)
       ELSE
         F1=(UN-HALF/X/X)/X
      END IF
      BET0=HALF/DENS(ID)/P(ID)
      BETP=HALF/DENS(ID+1)/P(ID+1)
      GAMA=UN/(DM(ID+1)-DM(ID))
      B(1,1)=F1
      B(1,2)=TWO*(X*F1-UN)*P(1)/HG1
      B(2,1)=BET0
      B(2,2)=GAMA
      C(1,1)=0.
      C(1,2)=0.
      C(2,1)=-BETP
      C(2,2)=GAMA
      VL(1)=DM(ID)*2.D0*VSND2(ID)/HG1-P(ID)*F1
      VL(2)=BET0*P(ID)+BETP*P(ID+1)-GAMA*(ZD(ID)-ZD(ID+1))
      ANU(1,ID)=0.
      ANU(2,ID)=0.
      CALL MATINV(B,2,2)
      DO I=1,2
         DO J=1,2
            S=0.
            DO K=1,2
               S=S+B(I,K)*C(K,J)
            END DO
            D(I,J,ID)=S
            ANU(I,ID)=ANU(I,ID)+B(I,J)*VL(J)
         END DO
      END DO
C
C     Normal depth points  1 < ID < ND
C
      DO ID=2,ND-1
         BET0=BETP
         BETP=HALF/DENS(ID+1)/P(ID+1)
         GAMA=UN/(DM(ID+1)-DM(ID))
         DMD=HALF*(DM(ID+1)-DM(ID-1))
         AA=UN/(DM(ID)-DM(ID-1))/DMD
         CC=GAMA/DMD
         BB=AA+CC
         BQ=QGRAV/P(ID)/DENS(ID)
         B(1,1)=BB+BQ-AA*D(1,1,ID-1)
         B(1,2)=-AA*D(1,2,ID-1)
         B(2,1)=BET0
         B(2,2)=GAMA
         C(1,1)=CC
         C(1,2)=0.
         C(2,1)=-BETP
         C(2,2)=GAMA
         VL(1)=AA*P(ID-1)+CC*P(ID+1)-(BB-BQ)*P(ID)+AA*ANU(1,ID-1)
         VL(2)=BET0*P(ID)+BETP*P(ID+1)-GAMA*(ZD(ID)-ZD(ID+1))
         CALL MATINV(B,2,2)
         ANU(1,ID)=0.
         ANU(2,ID)=0.
         DO I=1,2
            DO J=1,2
               S=0.
               DO K=1,2
                  S=S+B(I,K)*C(K,J)
               END DO
               D(I,J,ID)=S
               ANU(I,ID)=ANU(I,ID)+B(I,J)*VL(J)
            END DO
         END DO
      END DO  
C
C     Lower boundary condition
C
      ID=ND
      AA=TWO/(DM(ID)-DM(ID-1))**2
      BQ=QGRAV/P(ID)/DENS(ID)
      B(1,1)=AA+BQ-AA*D(1,1,ID-1)
      B(1,2)=-AA*D(1,2,ID-1)
      B(2,1)=0.
      B(2,2)=UN
      VL(1)=QGRAV/DENS(ID)+AA*(P(ID-1)-P(ID)+ANU(1,ID-1))
      VL(2)=0.
      CALL MATINV(B,2,2)
      ANU(1,ID)=0.
      ANU(2,ID)=0.
      DO I=1,2
         DO J=1,2
            S=0.
            DO K=1,2
               S=S+B(I,K)*C(K,J)
            END DO
            D(I,J,ID)=S
            ANU(I,ID)=ANU(I,ID)+B(I,J)*VL(J)
         END DO
      END DO
C
C     ------------
C     Backsolution
C     ------------
C
      P(ID)=P(ID)+ANU(1,ID)
      ZD(ID)=ZND
      CHMAXX=ABS(ANU(1,ID)/P(ID))
      DO IID=1,ND-1
         ID=ND-IID
         DO I=1,2
            DO J=1,2
               ANU(I,ID)=ANU(I,ID)+D(I,J,ID)*ANU(J,ID+1)
            END DO
         END DO
         CH1=ANU(1,ID)/P(ID)
         CH2=ANU(2,ID)/ZD(ID)
         IF(ABS(CH1).GT.CHMAXX) CHMAXX=ABS(CH1)
         IF(ABS(CH2).GT.CHMAXX) CHMAXX=ABS(CH2)
         IF(CH1.LT.-0.9D0) CH1=-0.9D0
         IF(CH1.GT.9.D0) CH1=9.D0
         P(ID)=P(ID)*(UN+CH1)
      END DO
C
C     Recalculate density for the new total pressure
C
      DO ID=1,ND
         DENS(ID)=P(ID)/VSND2(ID)
      END DO
C
C     New  z-values
C
      ZD(ND)=ZND
      DO IID=1,ND-1
         ID=ND-IID
         ZD(ID)=ZD(ID+1)+HALF*(DM(ID+1)-DM(ID))*(UN/DENS(ID)+
     *                   UN/DENS(ID+1))
      END DO
C
C     Convergence criterion for the Newton-Raphson method
C
      IF(IPRING.GE.1) WRITE(6,601) ITERH,CHMAXX
  601 FORMAT(/' solution of hydrostatic eq. + z-m relation:',
     *       'iter = ',I3,' max.rel.chan. =',1PD10.2)
      IF(CHMAXX.GT.ERROR.AND.ITERH.LT.10) GO TO 30
C
      DO ID=1,ND
         X=PGS(ID)/PTOTAL(ID)
         PTOTAL(ID)=P(ID)
         PGS(ID)=X*P(ID)
      END DO
C
C     Recalculation of the populations
C
      if(ih2p.ge.0) then
         ID=1
         ANEREL=ELEC(ID)/(DENS(ID)/WMM(ID)+ELEC(ID))
         DO ID=1,ND
            CALL RHONEN(ID,TEMP(ID),DENS(ID),AN,ANE)
            ELEC(ID)=ANE
            CALL WNSTOR(ID)
            CALL STEQEQ(ID,POP,1)
         END DO
      end if
      RETURN
      END