SUBROUTINE RHSGEN(ID)
C     =====================
C
C     Auxiliary procedure for SOLVE
C     controls evaluation of rhs vector when
C     matrices A,B, and C are kept fixed.
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'
      PARAMETER (XCON=8.0935D-21,YCON=1.68638E-10)
      PARAMETER (SIXTH=UN/6.D0,
     *           THIRD=UN/3.D0)
C     DIMENSION AJ(MLEVEL)
      DIMENSION POPP(MLEVEL)
      COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
C
      ispl=isplin
      if(isplin.ge.5) isplin=isplin-5
      IJ1=1
C
C     in the case of Compton scattering - boundary condition
C     for the highest frequency
C
      if(icompt.gt.0.and.icombc.gt.0.and.ijex(1).gt.0) then
         IJ1=2
       ij=1
       iji=nfreq
         zj1=exp(-hk*freq(ij)/temp(id))
         zj2=exp(-hk*freq(ij+1)/temp(id))
       dlt=delj(iji-1,id)
       if(ichcoo.eq.0) then
          zj0=un/(hk*sqrt(freq(ij)*freq(ij+1))/temp(id))
          zxx=un-3.*zj0+(un-dlt)*zj1+dlt*zj2
          combid=zj0/dlnfr(iji-1)+(un-dlt)*zxx
            comaid=-zj0/dlnfr(iji-1)+dlt*zxx
        else
            e2=ycon*temp(id)
            zxx0=xcon*freq(ij)*(un+zj1)-3.*e2
            zxxm=xcon*freq(ij+1)*(un+zj2)-3.*e2
            zxx=(un-dlt)*zxx0+dlt*zxxm
          combid=e2/dlnfr(iji-1)+(un-dlt)*zxx
          comaid=-e2/dlnfr(iji-1)+dlt*zxx
       end if
       vecl(ij)=-combid*rad(iji,id)-comaid*rad(iji-1,id)
      end if
C
C     evaluation of the opacity, emissivity, scattering, and
C     their derivatives, at the current depth point ID;
C
      DO I=1,NN
         PSI0(I)=PSY0(I,ID)
      END DO
      IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            WDEP0(IJ)=W(IJT)
            RAD0(IJ)=RADEX(IJ,ID)
            FK0(IJ)=FAKEX(IJ,ID)
            ABSO0(IJ)=ABSOEX(IJ,ID)
            EMIS0(IJ)=EMISEX(IJ,ID)
            SCAT0(IJ)=SCATEX(IJ,ID)
         END DO
      END IF
C
      IF(ID.GT.1) THEN
         DO I=1,NN
            PSIM(I)=PSY0(I,ID-1)
         END DO
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            RADM(IJ)=RADEX(IJ,ID-1)
            FKM(IJ)=FAKEX(IJ,ID-1)
            ABSOM(IJ)=ABSOEX(IJ,ID-1)
            EMISM(IJ)=EMISEX(IJ,ID-1)
            SCATM(IJ)=SCATEX(IJ,ID-1)
         END DO
         END IF
      END IF
C
      IF(ID.LT.ND) THEN
         DO I=1,NN
            PSIP(I)=PSY0(I,ID+1)
         END DO
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            RADP(IJ)=RADEX(IJ,ID+1)
            FKP(IJ)=FAKEX(IJ,ID+1)
            ABSOP(IJ)=ABSOEX(IJ,ID+1)
            EMISP(IJ)=EMISEX(IJ,ID+1)
            SCATP(IJ)=SCATEX(IJ,ID+1)
         END DO
         END IF
      END IF
c
C
C     ------------------------------------------------------------
C     Actual evaluation of the rhs vector VECL
C     ------------------------------------------------------------
C
      DO I=1,NN
         VECL(I)=0.
      END DO
C
      IF(NFREQE.LE.0) GO TO 100
C
C     -----------------------------------------------------------
C     1. Radiative transfer components
C     -----------------------------------------------------------
C
C     For  ID = 1  -  upper boundary condition
C
      IF(ID.GT.1) GO TO 50
      DDP=DELDMZ(1) 
      DO IJ=1,NFREQE
         IJT=IJFR(IJ)
         IF(IZSCAL.EQ.0) THEN
            OMEG0=ABSO0(IJ)/DENS(ID)
            OMEGP=ABSOP(IJ)/DENS(ID+1)
          ELSE
            OMEG0=ABSO0(IJ)
            OMEGP=ABSOP(IJ)
         END IF
         DZP=OMEG0+OMEGP
         DTAUP=DZP*DDP
         ALF1=(FK0(IJ)*RAD0(IJ)-FKP(IJ)*RADP(IJ))/DTAUP
         CHIEL0=SCAT0(IJ)
         CHIELP=SCATP(IJ)
         S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
         BS=HALF*DTAUP
         CS=0.
         C2=0.
         SP=0.
         BET2=0.
         GAM2=0.
c
c        additional terms for Compton scattering
c
         if(icompt.gt.0) then
          call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
          s0=s0+cms
       end if
C     
         IF(ISPLIN.GT.0) THEN
C
C      Spline collocation and/or Hermitian method (ISPLIN=1 or 2) -
C      both give the same expression for the boundary conditions
C
            BS=DTAUP*THIRD
            CS=HALF*BS
            SP=(EMISP(IJ)+CHIELP*RADP(IJ))/ABSOP(IJ)
            C2=CS/ABSOP(IJ)
            GAM2=CS*(RADP(IJ)-SP)
         END IF
C
C      auxiliary quantities
C
         ALF2=BS*(RAD0(IJ)-S0)
         BET2=ALF2+GAM2
C
C *** the IJ-th element of the rhs vector
C
         VECL(IJ)=ALF1+BET2+FH(IJT)*RAD0(IJ)-S0*Q0(IJT)
         IF(IWINBL.LT.0) VECL(IJ)=VECL(IJ)-HEXTRD(IJT)
C
      END DO  
      GO TO 100
C
C     For  1 < ID < ND  -  normal depth point
C
   50 DDM=DELDMZ(ID-1)
      IF(ID.EQ.ND) GO TO 80
      DDP=DELDMZ(ID)
      DO IJ=1,NFREQE
         IF(IZSCAL.EQ.0) THEN
            OMEG0=ABSO0(IJ)/DENS(ID)
            OMEGP=ABSOP(IJ)/DENS(ID+1)
            OMEGM=ABSOM(IJ)/DENS(ID-1)
          ELSE
            OMEG0=ABSO0(IJ)
            OMEGP=ABSOP(IJ)
            OMEGM=ABSOM(IJ)
         END IF
         DZP=OMEG0+OMEGP
         DZM=OMEG0+OMEGM
         DTAUP=DZP*DDP
         DTAUM=DZM*DDM
         DTAU0=HALF*(DTAUP+DTAUM)
         FRD=FK0(IJ)*RAD0(IJ)
         ALF1=(FRD-FKP(IJ)*RADP(IJ))/DTAUP/DTAU0
         GAM1=(FRD-FKM(IJ)*RADM(IJ))/DTAUM/DTAU0
         BET1=ALF1+GAM1
         BS=UN
         CHIELM=SCATM(IJ)
         CHIEL0=SCAT0(IJ)
         CHIELP=SCATP(IJ)
         S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
         BET2=0.
c
c        additional terms for Compton scattering
c
         if(icompt.gt.0) then
          call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
          s0=s0+cms
       end if
C     
         IF(ISPLIN.EQ.1) THEN
C
C     spline collocation (ISPLIN=1)
C
            AS=DTAUM/DTAU0*SIXTH
            CS=DTAUP/DTAU0*SIXTH
            BS=0.666666666666667D0
            SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ)
            SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ)
            ALF2=AS*(RADM(IJ)-SM)
            GAM2=CS*(RADP(IJ)-SP)
            BET2=ALF2+GAM2
          ELSE IF(ISPLIN.EQ.2) THEN
C
C     Hermitian method (ISPLIN=2)
C
            AS=DTAUP*DTAUP/DTAUM/DTAU0
            CS=DTAUM*DTAUM/DTAUP/DTAU0
            SM=(EMISM(IJ)+RADM(IJ)*CHIELM)/ABSOM(IJ)
            SP=(EMISP(IJ)+RADP(IJ)*CHIELP)/ABSOP(IJ)
            AS=(UN-HALF*AS)*SIXTH
            CS=(UN-HALF*CS)*SIXTH
            BS=UN-AS-CS
            BET2=AS*(RADM(IJ)-SM)+CS*(RADP(IJ)-SP)
         END IF
C
C *** the IJ-th element of the rhs vector
C
         VECL(IJ)=BET1+BET2+BS*(RAD0(IJ)-S0)
C
      END DO   
      GO TO 100
C
C   For ID=ND  -  lower boundary condition
C
   80 T=TEMP(ID)
      IF(TEMPBD.GT.0.) T=TEMPBD
      HKT=HK/T
C
c     the case of stellar atmospheres
C
      IF(IDISK.EQ.0.OR.IFZ0.LT.0) THEN
C
C     1. the case NDRE.le.ND, ie. radiative equilibrium equation at this
C        depth is treated as differential equation.
C
C     2. the case NDRE > ND, ie. radiative equilibrium equation at this
C        depth is treated as integral equation.
C        Boundary condition of the transfer equation is then complex,
C        because it has to incorporate the constraint of the total flux
C        equal to SIG4P*Teff**4 (integral form of the radiative
C        equation does not fix the total flux)
C
C        Details of this approach are presented in Mihalas, Heasley, Auer
C        NCAR-TN/STR-104 (1975).
C
C        This option is presented here for the sake of continuity with
C        previous work; otherwise the first option (NDRE.le.ND) is
C        recommended due to its better numerical properties
C
C        auxiliary quantities for the second option
C
         IF(ID.LT.NDRE.AND.INRE.NE.0) THEN
         SUMB=0.
         SUMF=0.
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            FR=FREQ(IJT)
            FR15=FR*1.D-15
            WD=WDEP0(IJ)
            X=HKT*FR
            EX=EXP(X)
            PLAN=BN*FR15*FR15*FR15/(EX-UN)*RRDIL
            DPLAN=PLAN*X/T/(UN-UN/EX)*WD
            SUMB=SUMB+(PLAN-RAD0(IJ))*WD
            FI=DPLAN/ABSO0(IJ)
            SUMF=SUMF+FI
         END DO
         FL=SIG4P*TEFF*TEFF*TEFF*TEFF
         ZZ=(FL-HALF*SUMB)/SUMF
         END IF
C
C        auxiliary quantities for both options
C
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            CHIELM=SCATM(IJ)
            CHIEL0=SCAT0(IJ)
            IF(IZSCAL.EQ.0) THEN
               OMEG0=ABSO0(IJ)/DENS(ID)
               OMEGM=ABSOM(IJ)/DENS(ID-1)
             ELSE
               OMEG0=ABSO0(IJ)
               OMEGM=ABSOM(IJ)
            END IF
            DZM=OMEG0+OMEGM
            DTAUM=DZM*DDM
            GAM1=(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))/DTAUM
            BET2=0.
C
C        second-order boundary condition
C
            IF(IBC.GT.0.AND.IBC.LT.4) THEN
               BS=DTAUM*HALF
               S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
c
c              additional terms for Compton scattering
c
               if(icompt.gt.0) then
                call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
                s0=s0+cms
              end if
C     
               BET2=BS*(RAD0(IJ)-S0)
            END IF
C
            FR=FREQ(IJT)
            X=HKT*FR
            EX=EXP(X)
            PLAN=BN*(FR*1.D-15)**3/(EX-UN)*RRDIL
            IF(INRE.EQ.0.OR.ID.GE.NDRE) THEN
               DPLAN=BN*(FR*1.D-15)**3/(EXP(HK*FR/TEMP(ID-1))-UN)
               GAM1=GAM1-(PLAN-DPLAN)/DTAUM*THIRD
            ELSE
               DPLAN=PLAN*X/T/(UN-UN/EX)
               FI=DPLAN/ABSO0(IJ)
               X1=FI*ZZ
               GAM1=GAM1-X1
            END IF
C
C *** the IJ-th element of the rhs vector
C
            IF(IBC.EQ.0.OR.IBC.EQ.4) THEN
               VECL(IJ)=GAM1+BET2-HALF*(PLAN-RAD0(IJ))
             ELSE
               VECL(IJ)=GAM1+BET2-HALF*PLAN+FHD(IJT)*RAD0(IJ)
            END IF
         END DO  
C
C     the case of accretion disks
C
      ELSE
C
      DO IJ=IJ1,NFREQE
         CHIELM=SCATM(IJ)
         CHIEL0=SCAT0(IJ)
         IF(IZSCAL.EQ.0) THEN
            OMEG0=ABSO0(IJ)/DENS(ID)
            OMEGM=ABSOM(IJ)/DENS(ID-1)
          ELSE
            OMEG0=ABSO0(IJ)
            OMEGM=ABSOM(IJ)
         END IF
         DZM=OMEG0+OMEGM
         DTAUM=DZM*DDM
         FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
         GAM1=FRD/DTAUM
         BS=DTAUM*HALF
         S0=(EMIS0(IJ)+CHIEL0*RAD0(IJ))/ABSO0(IJ)
         GAM2=BS*(RAD0(IJ)-S0)
C
C *** the IJ-th element of the rhs vector
C
         VECL(IJ)=GAM1+GAM2
      END DO
      END IF
C
C     -----------------------------------------------------------
C     2. Hydrostatic equilibrium
C     -----------------------------------------------------------
C
  100 isplin=ispl
      IF(INHE.EQ.0) GO TO 170
      NHE=NFREQE+INHE
      GRD=0.
C
      IF(ID.EQ.1) THEN
C
C     Upper boundary condition (ID=1)
C     Basically, linearized eq. (7-10) of Mihalas (1978)
C
         IF(IDISK.EQ.0.OR.IBCHE.EQ.0) THEN
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IJT=IJFR(IJ)
            IF(.NOT.LSKIP(ID,IJT)) THEN 
            FLUXW=W(IJT)*FH(IJT)*RAD0(IJ)
            GRD=GRD+FLUXW*ABSO0(IJ)
            END IF
         END DO
         END IF
C
         X1=PCK/DENS(ID)
         VT0=HALF*VTURB(ID)*VTURB(ID)/DM(ID)*WMM(ID)
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)*PSI0(NHE)/DM(ID)-
     *            X1*(GRD+FPRD(ID))-VT0/WMM(ID)*DENS(ID)
C
         ELSE
         IF(IBCHE.EQ.1) THEN
C
C        specifically disk - newer variant
C
         CCC=PCK/QGRAV
         HR1=CCC*SIG4P*TEFF**4*ABROSD(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=DENS(1)*HG1*F1
         VECL(NHE)=DM(1)-GGG
C
C        specifically disk - older variant
C
         ELSE IF(IBCHE.EQ.2) THEN
         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)
            GRD=GRD+FLUXW*ABSO0(IJ)
            END IF
         END DO
         END IF
C
         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        The rhs vector 
C
         VECL(NHE)=DM(1)*GGG-PG1
         ELSE
         VECL(NHE)=PGAS0-BOLK*TEMP(1)*PSI0(NHE)
         END IF
         END IF
C
      ELSE
C
C     Normal depth point (ID > 1)
C
C
         IF(NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            IF(.NOT.LSKIP(ID,IJFR(IJ))) 
     *      GRD=GRD+(FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ))*W(IJFR(IJ))
         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     the rhs vector
C     again, which accounts for the total radiation pressure in the
C     fixed-option transitions (array FPRD, generated by FIXLIN)
C
         IF(IDISK.EQ.1) GRAV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
         IF(IZSCAL.EQ.0) THEN
         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-1)*DENS(ID-1)
         ELSE
         GRAVZ=GRAV*(ZD(ID)-ZD(ID-1))
         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-1)*DENS(ID-1)
         END IF
      END IF
C
C     -----------------------------------------------------------
C     2a. z-m relation
C     -----------------------------------------------------------
C
  170 IF(INZD.LE.0.OR.ID.EQ.ND.OR.IDISK.EQ.0) GO TO 200
      NZD=NFREQE+INZD
      DDP=(DM(ID+1)-DM(ID))*HALF
      VECL(NZD)=ZD(ID+1)-ZD(ID)+DDP/DENS(ID)+DDP/DENS(ID+1)
C
C     -----------------------------------------------------------
C     3. Radiative equilibrium
C     -----------------------------------------------------------
C
  200 IF(INRE.EQ.0) GO TO 250
      NRE=NFREQE+INRE
c
      ittc=abs(nretc)/100
      if(iter.gt.ittc) then
      if(id.le.mod(abs(nretc),100)) then
         if(nretc.lt.0) vecl(nre)=temp(id+1)-temp(id)
         go to 250
      end if
      end if
C
C     integral equation part of the radiative
C     equilibrium equation
C
C     the rhs vector accounts for total net cooling in ALI transitions
C
      VECL(NRE)=FCOOL(ID)
      IF(IDISK.EQ.1) VECL(NRE)=FCOOL(ID)-TVISC(ID)*reint(id)
C
      IF(REINT(ID).GT.0.AND.NFREQE.GT.0) THEN
         DO IJ=1,NFREQE
            HEAT=ABSO0(IJ)-SCAT0(IJ)
            VECL(NRE)=VECL(NRE)-
     *                (HEAT*RAD0(IJ)-EMIS0(IJ))*WDEP0(IJ)*reint(id)
c
c           additional terms for Compton scattering
c
            if(icompt.gt.5) then
               ijt=ijfr(ij)
               call compt0(ijt,id,abso0(ij),cma,cmb,cmc,cme,cms,cmd)
               vecl(nre)=vecl(nre)+abso0(ij)*cms*wdep0(ij)*reint(id)
            end if
c
         END DO
      END IF
C
C
C     differential equation part of the
C     radiative equilibrium equation
C
      if(redif(id).gt.0) then
         TEFFD=TEFF**4
         IF(IDISK.EQ.1) TEFFD=TEFFD*(UN-THETAV(ID))
         VECL(NRE)=vecl(nre)+SIG4P*TEFFD*redif(id)
         IF(ID.GT.1) THEN
            DDM=(DM(ID)-DM(ID-1))*HALF 
            DDM=DELDMZ(ID-1)
            IF(NFREQE.GT.0) THEN
               DO IJ=1,NFREQE
                  OMEG0=ABSO0(IJ)*DENSI(ID)
                  OMEGM=ABSOM(IJ)*DENSI(ID-1)
                  DTAUM=(OMEG0+OMEGM)*DDM
                  FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
                  GAMR=FRD/DTAUM
                  VECL(NRE)=VECL(NRE)-WDEP0(IJ)*GAMR*redif(id)
               END DO
            END IF
          ELSE
            IF(NFREQE.GT.0) THEN
               DO IJ=1,NFREQE
                  IJT=IJFR(IJ)
                  VECL(NRE)=VECL(NRE)-
     *                      WDEP0(IJ)*FH(IJT)*RAD0(IJ)*REDIF(ID)
               END DO
            END IF
         END IF
      end if
C
C     -----------------------------------------------------------
C     4. Statistical equilibrium
C     -----------------------------------------------------------
C
  250 IF(IABS(IFPOPR).GE.3.and.ifpopr.le.5) THEN
         if(inse.gt.0) then
         NSE=NFREQE+INSE-1
         CALL SABOLF(ID)
         CALL LEVGRP(ID,IIEXP,0,POPP)
         CALL RATMAT(ID,IIEXP,0,ESEMAT,BESE)
         IF(IFPOPR.LE.3) THEN
            CALL MATINV(ESEMAT,NLVEXP,MLEVEL)
            DO I=1,NLVEXP
               SUM=0.
               DO J=1,NLVEXP
                  SUM=SUM+ESEMAT(I,J)*BESE(J)
               END DO
               VECL(NSE+I)=POPGRP(I)-SUM
               IF(IGZERO(I,ID).GT.0) VECL(NFREQE+INSE-1+I)=0.
            END DO
          ELSE 
            DO I=1,NLVEXP
               SUM=0.
               DO J=1,NLVEXP
                  SUM=SUM+ESEMAT(I,J)*POPGRP(J)
               END DO
               VECL(NSE+I)=BESE(I)-SUM
              IF(IGZERO(I,ID).GT.0) VECL(NFREQE+INSE-1+I)=0.
           END DO
         END IF
         end if
      END IF
C
C     -----------------------------------------------------------
C     5. charge conservation
C     -----------------------------------------------------------
C
      IF(INPC.EQ.0) GO TO 400
      NPC=NFREQE+INPC
C
C     This part is very similar to procedure ELCOR (obviously);
C     array AJ has the meaning of coefficients of the charge conserv.
C     ie. charge conservation is written
C           AJ * (vector of populations) = electron density
C
      T=TEMP(ID)
      ANE=ELEC(ID)
      CALL STATE(2,ID,T,ANE)
C
      VPC=QFIX(ID)+Q*DENS(ID)/WMM(ID)/YTOT(ID)
      IF(IOPTAB.EQ.0) VPC=VPC*ABUND(IATREF,ID)
      DO IAT=1,NATOM
         IF(IIFIX(IAT).NE.1) THEN
         DO I=N0A(IAT),NKA(IAT)
            IL=ILK(I)
            CH=IZ(IEL(I))-1
            IF(IL.GT.0) CH=IZ(IL)+(IZ(IL)-1)*ANE*USUMS(IL,ID)
            IF(IMODL(I).GE.0) VPC=VPC+CH*POPUL(I,ID)
         END DO
         END IF
      END DO
C
      VECL(NPC)=ANE-VPC
C
C     -----------------------------------------------------------
C     6. Convection
C     -----------------------------------------------------------
C
  400 IF(HMIX0.gt.0.) then
      NRE=NFREQE+INRE
      NDEL=NFREQE+INDL
C
C     Upper boundary condition (ID=1)
C
      ANEREL=ELEC(1)/(DENS(1)/WMM(1)+ELEC(1))
      IF(ID.EQ.1) THEN
         DELTA(ID)=0.
         FLXC(ID)=0.
      ELSE
C
C     Normal depth point 1 < ID < ND
C
         T=TEMP(ID)
         P=PTOTAL(ID)
         PG=PGS(ID)
         PRAD=P-PG-HALF*DENS(ID)*VTURB(ID)**2
         TM=TEMP(ID-1)
         PM=PTOTAL(ID-1)
         PGM=PGS(ID-1)
         PRADM=PM-PGM-HALF*DENS(ID-1)*VTURB(ID-1)**2
         T0=HALF*(T+TM)
         P0=HALF*(P+PM)
         PG0=HALF*(PG+PGM)
         PR0=HALF*(PRAD+PRADM)
         AB0=HALF*(ABROSD(ID)+ABROSD(ID-1))
         DLT=(T-TM)/(P-PM)*P0/T0
         DELTA(ID)=DLT
         IF(INDL.GT.0) VECL(NDEL)=DELTA(ID)-DLT
         IF(IDISK.EQ.1) GRAVD=ZD(ID)*QGRAV
C
C        convective flux
C
         CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
         FLXC(ID)=FLXCNV
         if(redif(id).gt.0) then
            VECL(NRE)=VECL(NRE)-FLXC(ID)*redif(id)
         end if
C
         if(reint(id).gt.0.AND.ICONV.LE.2) then
         TP=TEMP(ID+1)
         PM=PTOTAL(ID+1)
         PGP=PGS(ID+1)
         PRADP=PM-PGP-HALF*DENS(ID+1)*VTURB(ID+1)**2
         T0=HALF*(T+TP)
         P0=HALF*(P+PM)
         PG0=HALF*(PG+PGP)
         PR0=HALF*(PRAD+PRADP)
         AB0=HALF*(ABROSD(ID)+ABROSD(ID+1))
         DLT=(TP-T)/(PM-P)*(PM+P)/(TP+T)
         CALL CONVEC(ID+1,T0,P0,PG0,PR0,AB0,DLT,FLXCP,VCON)
         DELM=(DM(ID+1)-DM(ID))*HALF
         RDELM=DENS(ID)/DELM
         VECL(NRE)=VECL(NRE)-RDELM*(FLXCP-FLXCNV)*reint(id)
         end if
      END IF
      END IF
C
C     skip rows corresponding to fully-zeroed populations
C
      NSE=NFREQE+INSE
      INONZ=NSE
      DO II=NSE,NN0
         IF(IGZERT(II-NSE+1).EQ.0) THEN
            IF(INONZ.NE.II) VECL(INONZ)=VECL(II)
            INONZ=INONZ+1
         END IF
      END DO
C
      RETURN
      END