SUBROUTINE BHE(ID)
C     ==================
C
C     The part of matrices A and B corresponding to the hydrostatic
C     equilibrium equation,
C     i.e. the (NFREQE+INHE)-th row;
C     and, if desired (INMP > 0), the part corresponding to the
C     definition equation for the fictitious massive particle density,
C     ie. the (NFREQE+INMP)-th row.
C
C     Input:  ID - depth index
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'ATOMIC.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ARRAY1.FOR'
      INCLUDE 'ALIPAR.FOR'
C
      NHE=NFREQE+INHE
      NRE=NFREQE+INRE
      NPC=NFREQE+INPC
      NSE=NFREQE+INSE-1
c
c     the case of fixed mass density
c
      if(ifixde.gt.0) then
         b(nhe,nhe)=un
         b(nhe,npc)=-un
         vecl(nhe)=dens(id)/wmm(id)+elec(id)-totn(id) 
         return
      end if
C
C ***********  Linearized equation for the fictitious massive particle
C              density
C
      IF(INMP.GT.0) THEN
         NMP=NFREQE+INMP
         B(NMP,NMP)=-UN
         B(NMP,NHE)=UN
         IF(INPC.GT.0) B(NMP,NPC)=-UN
      END IF
C
C ***********  Linearized hydrostatic equilibrium
C
      HEXT=0.
      HEXN=0.
      GRD=0.
      FLUXW=0.
      IF(ID.GT.1) GO TO 50
C
C *** Upper boundary condition (ID=1)
C     Basically, linearized eq. (7-10) of Mihalas (1978)
C
      DO I=1,NLVEXP
         HEX(I)=0.
      END DO
      x1=0.
      IF(NFREQE.GT.0.AND.IFPRAD.GT.0) THEN
         X1=PCK/DENS(ID)
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            IF(.NOT.LSKIP(ID,IJT)) THEN
               FLUXW=W(IJT)*(FH(IJT)*RAD0(IJ)-HEXTRD(IJT))
               GRD=GRD+FLUXW*ABSO0(IJ)
               HEXN=HEXN+FLUXW*DABN0(IJ)
               HEXT=HEXT+FLUXW*DABT0(IJ)
               DO I=1,NLVEXP
                  HEX(I)=HEX(I)+FLUXW*DRCH0(I,IJ)
               END DO
C
C              Columns corresponding to mean intensities
C
               B(NHE,IJ)=X1*W(IJT)*FH(IJT)*ABSO0(IJ)
            END IF
         END DO
      END IF
C
      RTN=X1*WMM(ID)/DENS(ID)*(GRD+FPRD(ID))
      VT0=HALF*VTURB(ID)*VTURB(ID)/DM(ID)*WMM(ID)
C
C     columns corresponding to total particle density, fictitious
C     massive particle density, temperature, and electron density,
C     respectively
C
      B(NHE,NHE)=BOLK*TEMP(ID)/DM(ID)-GN*(RTN-VT0)
      IF(INMP.GT.0) B(NHE,NFREQE+INMP)=GP*(VT0-RTN)
      IF(INRE.GT.0) THEN
         B(NHE,NRE)=BOLK*TOTN(ID)/DM(1)+X1*(HEXT+HEIT(ID))
         C(NHE,NRE)=X1*HEITP(ID)
      END IF
      IF(INPC.GT.0) THEN
         B(NHE,NPC)=X1*(HEXN+HEIN(ID))+GN*(RTN-VT0)
         C(NHE,NPC)=X1*HEINP(ID)
      END IF
C
C     Columns corresponding to populations
C
      DO II=1,NLVEXP
          B(NHE,NSE+II)=B(NHE,NSE+II)+X1*(HEX(II)+HEIP(II,ID))
          C(NHE,NSE+II)=C(NHE,NSE+II)+X1*HEIPP(II,ID)
      END DO
C
C     The rhs vector also accounts for the total radiation pressure in
C     the fixed-option transitions (array FPRD, generated by FIXLIN)
C
      VECL(NHE)=GRAV-BOLK*TEMP(ID)*TOTN(ID)/DM(ID)-
     *          X1*(GRD+FPRD(ID))-VT0/WMM(ID)*DENS(ID)
      RETURN
C
C *** Normal depth point (ID > 1)
C
C     Columns (for matrices A and B) corresponding to mean intensities
C
   50 CONTINUE
      IF(NFREQE.GT.0.and.ifprad.gt.0) THEN
      DO IJ=1,NFREQE
         IF(.NOT.LSKIP(ID,IJFR(IJ))) THEN
         GRD=GRD+(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))*W(IJFR(IJ))
         A(NHE,IJ)=-PCK*W(IJFR(IJ))*FKM(IJ)
         B(NHE,IJ)=PCK*W(IJFR(IJ))*FK0(IJ)
         END IF
      END DO
      END IF
C
      VT0=HALF*VTURB(ID)*VTURB(ID)*WMM(ID)
      VTM=HALF*VTURB(ID-1)*VTURB(ID-1)*WMM(ID-1)
C
C     columns corresponding to total particle density
C
      A(NHE,NHE)=-BOLK*TEMP(ID-1)-GN*VTM
      B(NHE,NHE)=BOLK*TEMP(ID)+GN*VT0
C
C     columns corresponding to temperature
C
      IF(INRE.GT.0) THEN
         A(NHE,NRE)=-BOLK*TOTN(ID-1)+PCK*HEITM(ID)
         B(NHE,NRE)=BOLK*TOTN(ID)+PCK*HEIT(ID)
         C(NHE,NRE)=PCK*HEITP(ID)
      END IF
C
C     columns corresponding to electron density
C
      IF(INPC.GT.0) THEN
         A(NHE,NPC)=GN*VTM+PCK*HEINM(ID)
         B(NHE,NPC)=-GN*VT0+PCK*HEIN(ID)
         C(NHE,NPC)=PCK*HEINP(ID)
      END IF
C
C     columns corresponding to NMP
C
      IF(INMP.GT.0) THEN
         A(NHE,NFREQE+INMP)=-GP*VTM
         B(NHE,NFREQE+INMP)=GP*VT0
      END IF
C
C     columns corresponding to populations
C
      DO II=1,NLVEXP
         A(NHE,NSE+II)=A(NHE,NSE+II)+PCK*HEIPM(II,ID)
         B(NHE,NSE+II)=B(NHE,NSE+II)+PCK*HEIP(II,ID)
         C(NHE,NSE+II)=C(NHE,NSE+II)+PCK*HEIPP(II,ID)
      END DO
C
C     the rhs vector
C     again, which accounts for the total radiation pressure in the
C     fixed-option transitions (array FPRD)
C
      VECL(NHE)=GRAV*(DM(ID)-DM(ID-1))-
     *          BOLK*(TEMP(ID)*TOTN(ID)-TEMP(ID-1)*TOTN(ID-1))-
     *          PCK*(GRD+FPRD(ID))-
     *          VT0/WMM(ID)*DENS(ID)+VTM/WMM(ID-1)*DENS(ID-1)
      RETURN
      END