SUBROUTINE BHEZ(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
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)
C
      NHE=NFREQE+INHE
      NRE=NFREQE+INRE
      NPC=NFREQE+INPC
      NZD=NFREQE+INZD
      NSE=NFREQE+INSE-1
c
      if(inhe.le.0) return
      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.EQ.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)
C
      RETURN
      ELSE IF(IBCHE.EQ.1) THEN
C
C     2. possibility - specifically disk - Hubeny (1990), Eq. (4.19)
C        newer 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)
         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
      RETURN
      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
      RETURN
      ELSE
C
C     4. a simple form of the bounary condition P_gas(ID=1)=PGAS0,
C        where PGAS0 is an input parameter
C
      B(NHE,NHE)=BOLK*TEMP(1)
      IF(INRE.GT.0) B(NHE,NRE)=BOLK*PSI0(NHE) 
      VECL(NHE)=PGAS0-BOLK*TEMP(1)*PSI0(NHE)
      RETURN
      END IF
C
C *** Normal depth point (ID > 1)
C
C     Columns (for matrices A and B) corresponding to mean intensities
C
   50 CONTINUE            
      GRAV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
      GRAVZ=GRAV*(ZD(ID)-ZD(ID-1))
      DGRV=GRAVZ*HALF*WMM(ID)
      GRD=0.
      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+DGRV)
      B(NHE,NHE)=BOLK*TEMP(ID)+GN*(VT0+DGRV)
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+DGRV)+PCK*HEINM(ID)
         B(NHE,NPC)=-GN*(VT0+DGRV)+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+DGRV)
         B(NHE,NFREQE+INMP)=GP*(VT0+DGRV)
      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
      VECL(NHE)=-GRAVZ*(DENS(ID)+DENS(ID-1))*HALF-
     *          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
      RETURN
      END