SUBROUTINE BHED(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     ii) 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     iii) the part of matrices B and C corresponding to the
C     z-m (z-distance versus mass-depth coordinate) relation,
C     ie. the (NFREQE+INZD)-th row of matrices B and C, however, the
C     elements of C are treated separately
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'
      COMMON/SURFEX/EXTJ(MFREQ),EXTH(MFREQ)
      COMMON/CMATZD/CZZ,CZN,CZE,CZM
C
      NHE=NFREQE+INHE
      NRE=NFREQE+INRE
      NPC=NFREQE+INPC
      NSE=NFREQE+INSE-1
c
      if(inhe.le.0) go to 100
      IJ1=1
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.
      DO I=1,NLVEXP
         HEX(I)=0.
      END DO
C
      IF(ID.GT.1) GO TO 50
C
C *** Upper boundary condition (ID=1)
C
C     1. possibility - the same as in stellar atmospheres
C     Basically, linearized eq. (7-10) of Mihalas (1978)
C
      IF(IBCHE.LE.0) THEN
      X1=PCK/DENS(ID)
      IF(NFREQE.GT.0) THEN
      DO IJ=IJ1,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*WDEP0(IJ)*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*PSI0(NHE)/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)
C
      GRAV=QGRAV*ZD(1)
      VECL(NHE)=GRAV-BOLK*TEMP(ID)*PSI0(NHE)/DM(ID)-
     *          X1*(GRD+FPRD(ID))-VT0/WMM(ID)*DENS(ID)
      GO TO 100
      ELSE IF(IBCHE.EQ.1) THEN
C
C     2. possibility - specifically disk - Hubeny (1990), Eq. (4.19)
C        newer variant
C
C
      IF(NFREQE.GT.0) THEN
      DO IJ=IJ1,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
         END IF
      END DO
      END IF
C
      CCC=PCK/QGRAV
      HR1=CCC*(GRD+FPRD(1))/DENS(1)
      PG1=BOLK*PSI0(NHE)*TEMP(1)
      HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
      X=(ZD(1)-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
      GGG=DENS(1)*HG1*F1
      RF1=DENS(1)*F1D
      CCD=CCC*F1D
C
      DO IJ=1,NFREQE
         B(NHE,IJ)=-CCD*WDEP0(IJ)*FH(IJFR(IJ))*ABSO0(IJ)
      END DO
C
C     columns corresponding to total particle density and temperature
C
      B(NHE,NHE)=B(NHE,NHE)+(GGG+HR1*RF1)/PSI0(NHE)
      IF(INRE.GT.0) B(NHE,NRE)=
     *   (GGG-RF1*ZD(1)+RF1*HR1)*HALF/TEMP(1)-CCD*(HEXT+HEIT(ID))
      IF(INZD.GT.0) B(NHE,NZD)=RF1
      IF(INPC.GT.0) B(NHE,NPC)=-CCD*(HEXN+HEIN(ID))
      DO II=1,NLVEXP
         B(NHE,NSE+II)=-CCD*(HEX(II)+HEIP(II,ID))
      END DO
C
C     The rhs vector 
C
      VECL(NHE)=DM(1)-GGG
      GO TO 100
      ELSE IF(IBCHE.EQ.2) THEN
C
C     3. possibility - specifically disk - Hubeny (1990), Eq. (4.19)
C        older variant
C
      IF(NFREQE.GT.0) THEN
      DO IJ=IJ1,NFREQE
         IJT=IJFR(IJ)
         IF(.NOT.LSKIP(ID,IJT)) THEN 
            FLUXW=W(IJT)*(FH(IJT)*RAD0(IJ)-HEXTRD(IJT))
            GRD=GRD+FLUXW*ABSO0(IJ)
         END IF
      END DO
      END IF
      CCC=PCK/QGRAV
      PR1=CCC*(GRD+FPRD(1))/DENS(1)
      PG1=BOLK*PSI0(NHE)*TEMP(1)
      HG1=SQRT(TWO*PG1/DENS(1)/QGRAV)
      X=(ZD(1)-PR1)/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
      GGG=HG1*QGRAV*HALF/F1
C
C     columns corresponding to total particle density and temperature
C
      B(NHE,NHE)=BOLK*TEMP(1)
      IF(INRE.GT.0) B(NHE,NFREQE+INRE)=PG1/TEMP(1)
C
C     The rhs vector 
C
      VECL(NHE)=DM(1)*GGG-PG1
      GO TO 100
      END IF
C
C *** Normal depth point (ID > 1)
C
C     Columns (for matrices A and B) corresponding to mean intensities
C
   50 IF(NFREQE.GT.0) THEN
      DO IJ=IJ1,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)
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*PSIM(NHE)+PCK*HEITM(ID)
         B(NHE,NRE)=BOLK*PSI0(NHE)+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     column corresponding to ZD (z-distance)
C
      IF(INZD.GT.0) THEN
         A(NHE,NFREQE+INZD)=-QGRAV*(DM(ID)-DM(ID-1))*HALF
         B(NHE,NFREQE+INZD)=-QGRAV*(DM(ID)-DM(ID-1))*HALF
      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
C     Since ZD(ID) is in fact z-distance corresponding to a midpoint
C     between depth points ID and ID+1 (as follows from the numerical
C     representation of the relation between DM and ZD), and ZD(ID-1)
C     corresponds to a midpoint between ID and ID-1, z-distance for
C     the point ID is better approximated by the mean value of ZD(ID)
C     and ZD(ID-1)
C
      GRAV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
      VECL(NHE)=GRAV*(DM(ID)-DM(ID-1))-
     *          BOLK*(TEMP(ID)*PSI0(NHE)-TEMP(ID-1)*PSIM(NHE))-
     *          PCK*(GRD+FPRD(ID))-
     *          VT0/WMM(ID)*DENS(ID)+VTM/WMM(ID)*DENS(ID-1)
C
C ***********  Linearized z-m (z-distance vers. mass-depth) relation
C
C     Note: since there are only at most four non-zero elements of
C     matrix C, they are stored separately in CZZ,CZN,CZE,CZM;
C     when multiplying any matrix by matrix C, these terms must be
C     treated separately - see SOLVE
C
  100 IF(INZD.LE.0) RETURN
      NZD=NFREQE+INZD
C
C *** lower boundary condition [ie. z(ND)=0 ]
C
      B(NZD,NZD)=UN
      IF(ID.EQ.ND) RETURN
C
C *** normal depth point
C
      DDP=(DM(ID+1)-DM(ID))*HALF
C
C     column corresponding to ZD
C
      B(NZD,NZD)=UN
      CZZ=-UN
C
C     column corresponding to total particle density
C
      X1=GN*WMM(ID)*DDP
      IF(INHE.GT.0) THEN
         B(NZD,NFREQE+INHE)=X1/DENS(ID)/DENS(ID)
         CZN=X1/DENS(ID+1)/DENS(ID+1)
      END IF
C
C     column corresponding to electron density
C
      IF(INPC.GT.0) THEN
         B(NZD,NFREQE+INPC)=-X1/DENS(ID)/DENS(ID)
         CZE=-X1/DENS(ID+1)/DENS(ID+1)
      END IF
C
C     column corresponding to NMP
C
      IF(INMP.GT.0) THEN
         B(NZD,NFREQE+INMP)=DDP/DENS(ID)/PSI0(NFREQE+INMP)
         CZM=DDP/DENS(ID+1)/PSIP(NFREQE+INMP)
      END IF
C
C     the element of the rhs vector
C
      VECL(NZD)=ZD(ID+1)-ZD(ID)+DDP/DENS(ID)+DDP/DENS(ID+1)
      RETURN
      END