SUBROUTINE RYBENE
C     =================
C
C     complementing partial matrices of complete linearization
c     corresponding to the contribution of ALI frequencies and the
c     convection to the radiative/convective equilibrium equation in
c     the Rybicki formalism
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ALIPAR.FOR'
      INCLUDE 'ARRAY1.FOR'
      COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
      COMMON/RYBMTX/RA(MDEPTH),RB(MDEPTH),RC(MDEPTH),VR(MDEPTH),
     *              UA(MDEPTH),UB(MDEPTH),UC(MDEPTH),    
     *              VA(MDEPTH),VB(MDEPTH),VC(MDEPTH),WR(MDEPTH), 
     *              WM(MDEPTH,MDEPTH)
      common/deridt/dert
C
C     contribution from the radiative equilibrium part of the energy eq.
c
      flxto0=sig4p*teff**4
      DO ID=1,ND
         WR(ID)=WR(ID)+FCOOL(ID)
         IF(IDISK.EQ.1) WR(ID)=WR(ID)-reint(id)*TVISC(ID)
         if(reint(id).gt.0.) then
            IF(ID.GT.1) WM(ID,ID-1)
     *                 =WM(ID,ID-1)+AREIT(ID)*dens(id)*reint(id)
            WM(ID,ID)  =WM(ID,ID)  + REIT(ID)*dens(id)*reint(id)
            IF(IDISK.EQ.1) WM(ID,ID)=WM(ID,ID)+
     *                     DTVIST(ID)*reint(id)
            IF(ID.LT.ND) WM(ID,ID+1)
     *                 =WM(ID,ID+1)+CREIT(ID)*dens(id)*reint(id)
         end if
c
         if(idisk.eq.0) then
            flxtot=flxto0
            gravd=grav
          else
            flxtot=flxto0*(1.d0-thetav(id))
            gravd=qgrav*zd(id)
         end if
c
         if(redif(id).gt.0) then
            WR(ID)=WR(ID)+flxtot*redif(id)
            IF(ID.GT.1) WM(ID,ID-1)=WM(ID,ID-1)+REDTM(ID)*REDIF(ID)
            WM(ID,ID)=  WM(ID,ID)+  REDT (ID)*REDIF(ID)
            IF(ID.LT.ND) WM(ID,ID+1)=WM(ID,ID+1)+REDTP(ID)*REDIF(ID)
         end if
      END DO
c
C     contribution from convection
C
      IF(HMIX0.LE.0.OR.ICONV.LE.0) RETURN
      if(dert.eq.0.) dert=0.01
      DO 10 ID=idconz,ND
         T=TEMP(ID)
         P=PTOTAL(ID)
         TM=TEMP(ID-1)
         PM=PTOTAL(ID-1)
         pr0=0.
         if(icentr.eq.0.or.id.eq.nd) then
            T0=SQRT(T*TM)
            P0=SQRT(P*PM)
            AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
            pr0=sqrt(pradt(id)*pradt(id-1))*pck
            DLP=UN/LOG(P/PM)
            DLT=LOG(T/TM)*DLP
            DELTA(ID)=DLT
            DDT0= DLP/T
            DDTM=-DLP/TM
            CALL CONVEC(ID,T0,P0,P0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
            FLXC(ID)=FLXCNV
            DHCDD=0.
            IF(DELMDE.GT.0.) THEN
               bb=bcnv*half
               cord=un-bb/sqrt(bb*bb+(dlt-grdadb))
               DHCDD=1.5D0/DELMDE*FLXCNV*cord
            END IF
            if(id.lt.icbegp-2) then
               flxc(id)=0.
               go to 10
            end if
            T1=(un+dert)*T0
            CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
            DHCDT0=(FLXC1-FLXCNV)*HALF/dert
            DHCDT =DHCDT0/T+DHCDD*DDT0
            DHCDTM=DHCDT0/TM+DHCDD*DDTM
          else
            dlm=log(p/pm)
            dlp=log(ptotal(id+1)/p)
            dl0=dlm+dlp
            palf=dlp/dlm/dl0
            pgam=dlm/dlp/dl0
            dlt=palf*log(t/tm)+pgam*log(temp(id+1)/t)
            delta(id)=dlt
            ddtm=-palf/tm
            ddtp=pgam/temp(id+1)
            ddt0=(palf-pgam)/t
            pr0=sqrt(pradt(id)*pradt(id-1))*pck
            call convec(id,t,p,p-pr0,pr0,abrosd(id),dlt,flxcnv,vcon)
            flxc(id)=flxcnv
            dhcdd=0.
            IF(DELMDE.GT.0.) then
               bb=bcnv*half
               cord=un-bb/sqrt(bb*bb+(dlt-grdadb))
               DHCDD=1.5D0/DELMDE*FLXCNV*cord
            end if
            if(id.lt.icbegp-2) then
               flxc(id)=0.
               go to 10
            end if
            T1=(un+dert)*T
            CALL CONVEC(ID,T1,P,p-pr0,pr0,abrosd(id),DLT,FLXC1,VCON)
            DHCDT0=(FLXC1-FLXCNV)*HALF/dert
            DHCDT =DHCDT0/T+DHCDD*DDT0
            DHCDTM=DHCDD*DDTM
            DHCDTP=DHCDD*DDTP
         end if
    
C
C        additional terms in matrix WM
C
C    **  1. differential equation form
C
         if(redif(id).gt.0) then
            WM(ID,ID-1)=WM(ID,ID-1)+DHCDTM*redif(id)
            WM(ID,ID)=  WM(ID,ID)+DHCDT*redif(id)
            WR(ID)=WR(ID)-FLXC(ID)*redif(id)
            if(icentr.ne.0.and.id.lt.nd) then
               WM(ID,ID+1)=WM(ID,ID+1)+DHCDTP*redif(id)
            end if
        END IF
C
C    **  2. integral equation form
C
         if(reint(id).gt.0.) then
         if(id.lt.nd) then
            TP=TEMP(ID+1)
            PTP=PTOTAL(ID+1)
            T0=SQRT(T*TP)
            P0=SQRT(P*PTP)
            pr0=sqrt(pradt(id)*pradt(id-1))*pck
            AB0=SQRT(ABROSD(ID)*ABROSD(ID+1))
            DLP=UN/LOG(PTP/P)
            DLT=LOG(TP/T)*DLP
            DDTP0= DLP/TP
            DDTPM=-DLP/T
            CALL CONVEC(ID,T0,P0,p0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
            DHCDDP=0.
            IF(DELMDE.GT.0.) DHCDDP=1.5D0/DELMDE*FLXCNV
            T1=1.001D0*T0
            CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
            DHCDT0=(FLXC1-FLXCNV)*500.
            DHCDTP=DHCDT0/TP+DHCDDP*DDTP0
            DHCDTU=DHCDT0/T+DHCDDP*DDTPM
C
C           additional terms in matrices A and B (the row corresponding to
C           energy balance, i.e. T) due to convection
C
            DELM=(DM(ID+1)-DM(ID-1))*HALF
            RDELM=DENS(ID)/DELM*reint(id)
            IF(ICONV.GT.0) THEN
               WM(ID,ID-1)=WM(ID,ID-1)-RDELM*DHCDTM
               WM(ID,ID)=WM(ID,ID)+RDELM*(DHCDTP-DHCDTU)
               WM(ID,ID+1)=WM(ID,ID+1)+RDELM*DHCDTP
            END IF
            WR(ID)=WR(ID)-RDELM*(FLXCNV-FLXC(ID))
          ELSE
            TP=TEMP(ID-2)
            PTP=PTOTAL(ID-2)
            T0=SQRT(T*TP)
            P0=SQRT(P*PTP)
            pr0=sqrt(pradt(id)*pradt(id-1))*pck
            AB0=SQRT(ABROSD(ID)*ABROSD(ID-2))
            DLP=UN/LOG(PTP/P)
            DLT=LOG(TP/T)*DLP
            DDTP0= DLP/TP
            DDTPM=-DLP/T
            CALL CONVEC(ID,T0,P0,p0-pr0,pr0,AB0,DLT,FLXCNV,VCON)
            if(flxcnv.le.0.) go to 10
               DHCDDP=0.
               IF(DELMDE.GT.0.) DHCDDP=1.5D0/DELMDE*FLXCNV
               T1=1.001D0*T0
               CALL CONVEC(ID,T1,P0,p0-pr0,pr0,AB0,DLT,FLXC1,VCON)
               DHCDTP=(FLXC1-FLXCNV)*1.D3/T0*HALF+DHCDDP*DDTP0
C
C           additional terms in matrices A and B (the row corresponding to
C           energy balance, i.e. T) due to convection
C
               DELM=(DM(ID)-DM(ID-2))*HALF
               RDELM=DENS(ID)/DELM*reint(id)
               DELHC=WMM(ID)/DELM*(FLXCNV-FLXC(ID))
               WM(ID,ID-1)=WM(ID,ID-1)+RDELM*(DHCDT-DHCDTP)
               WM(ID,ID)=WM(ID,ID)+RDELM*DHCDT
               WR(ID)=WR(ID)-RDELM*(FLXC(ID)-FLXCNV)
            END IF
         END IF
   10 CONTINUE
c
      RETURN
      END