SUBROUTINE MATCON(ID)
C     =====================
C
C     Evaluation of the terms in matrices A and B of complete
C     linearization due to convection.
C     Two rows are modified:
C      NRE  - energy balance (i.e. temperature)
C      NDEL - new row corresponding to a new model parameter DELTA
C
C     Model parameter DELTA is defined as dln(T)/dln(P);
C     FCONV is the convective flux
C
C     The equations and corresponding matrix elements are similar
C     to those considered by Grenfell, Astr.Ap. 20, 293 (1972).
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      INCLUDE 'ARRAY1.FOR'
      COMMON/CUBCON/ACNV,BCNV,DEL,GRDADB,DELMDE,RHO,FLXTOT,GRAVD
C
      IF(HMIX0.LE.0.) RETURN
      NHE=NFREQE+INHE
      NRE=NFREQE+INRE
      NPC=NFREQE+INPC
      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.
         IF(INDL.GT.0) B(NDEL,NDEL)=UN
        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
         if(ilgder.eq.0) then
            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
            TT=T*T-TM*TM
            DDT0= DLT/HALF*TM/TT
            DDTM=-DDT0*T/TM
          else
            T0=SQRT(T*TM)
            P0=SQRT(P*PM)
            PG0=SQRT(PG*PGM)
            PR0=SQRT(PRAD*PRADM)
            AB0=SQRT(ABROSD(ID)*ABROSD(ID-1))
            DLP=UN/LOG(P/PM)
            DLT=LOG(T/TM)*DLP
            DDT0= DLP/T
            DDTM=-DLP/TM
         end if         
         DELTA(ID)=DLT
         IF(INDL.GT.0) THEN
            B(NDEL,NDEL)=-UN
            VECL(NDEL)=DELTA(ID)-DLT
            DDP0=0.
            DDPM=0.
            IF(IPRESS.GT.0) THEN
                if(ilgder.eq.0) then
                   PP0=P*P-PM*PM
                   DDP0=-DLT/HALF*PM/PP0
                   DDPM=-DDP0*P/PM
                 else
                   PP0=LOG(T/TM)*DLP*DLP
                   DDP0=-PP0/P
                   DDPM=PP0/PM
                end if
            END IF
C
C           the row of matrix A corresponding to DELTA
C
            if(inhe.gt.0) A(NDEL,NHE)=BOLK*TM*DDPM
            A(NDEL,NRE)=PGM/TM*DDPM+DDTM
C
C           the row of matrix B corresponding to DELTA
C
            if(inhe.gt.0) B(NDEL,NHE)=BOLK*T*DDP0
            B(NDEL,NRE)=PG/T*DDP0+DDT0
         END IF
C
C        convective flux and its derivatives
C        (derivativs wrt. T and P(gas) - calculated numerically;
C        derivative wrt. DELTA calculated analytically)
C
         if(idisk.eq.1) gravd=zd(id)*qgrav
         CALL CONVEC(ID,T0,P0,PG0,PR0,AB0,DLT,FLXCNV,VCON)
         FLXC(ID)=FLXCNV
         DHCDD=0.
         IF(DELMDE.GT.0.) DHCDD=1.5D0/DELMDE*FLXCNV
         T1=1.001D0*T0
         CALL CONVEC(ID,T1,P0,PG0,PR0,AB0,DLT,FLXC1,VCON)
         DHCDT0=(FLXC1-FLXCNV)*1.D3*HALF
         DHCDT =DHCDT0/T
         DHCDTM=DHCDT0/TM
         DHCDP=0.
         IF(IPRESS.GT.0) THEN
            PG1=1.001*PG0
            CALL CONVEC(ID,T0,P0,PG1,PR0,AB0,DLT,FLXC2,VCON)
            DHCDP=(FLXC2-FLXCNV)*1.D3/PG0*HALF
            IF(IPRESS.GT.1) THEN
               P1=1.001*P0
               CALL CONVEC(ID,T0,P1,PG0,PR0,AB0,DLT,FLXC3,VCON)
               DHCDPT=(FLXC3-FLXCNV)*1.D3/P0*HALF
               DHCDP=DHCDP+DHCDPT
               DHCDT=DHCDT+DHCDPT*4.*PR0/T0
            END IF
         END IF
         IF(INDL.EQ.0) THEN
            DHCDT=DHCDT+DHCDD*DDT0
            DHCDTM=DHCDTM+DHCDD*DDTM
         END IF
C
C        additional terms in matrices A and B (the row corresponding to
C        energy balance, i.e. T) due to convection
C
C
C    **  1. differential equation form
C
         if(redif(id).gt.0) then
         IF(ICONV.GT.0) THEN
         IF(INHE.GT.0) A(NRE,NHE)=A(NRE,NHE)-DHCDP*BOLK*TM*redif(id)
         A(NRE,NRE)=A(NRE,NRE)-(DHCDP*PGM/TM+DHCDTM)*redif(id)
         IF(INHE.GT.0) B(NRE,NHE)=B(NRE,NHE)+DHCDP*BOLK*T*redif(id)
         B(NRE,NRE)=B(NRE,NRE)+(DHCDP*PG/T+DHCDT)*redif(id)
         IF(INDL.GT.0) B(NRE,NDEL)=B(NRE,NDEL)+DHCDD*redif(id)
         END IF
         VECL(NRE)=VECL(NRE)-FLXC(ID)*redif(id)
         end if
C
C    **  2. integral equation form
C
         if(reint(id).gt.0.AND.ICONV.LE.2) then
         TP=TEMP(ID+1)
         PTP=PTOTAL(ID+1)
         PGP=PGS(ID+1)
         PRADP=PTP-PGP-HALF*DENS(ID+1)*VTURB(ID+1)**2
         if(ilgder.eq.0) then
            T0=HALF*(T+TP)
            P0=HALF*(P+PTP)
            PG0=HALF*(PG+PGP)
            PR0=HALF*(PRAD+PRADP)
            AB0=HALF*(ABROSD(ID)+ABROSD(ID+1))
            DLT=(TP-T)/(PTP-P)*P0/T0
           else
            T0=SQRT(T*TP)
            P0=SQRT(P*PTP)
            PG0=SQRT(PG*PGP)
            PR0=SQRT(PRAD*PRADP)
            AB0=SQRT(ABROSD(ID)*ABROSD(ID+1))
            DLP=UN/LOG(PTP/P)
            DLT=LOG(TP/T)*DLP
            DDTP0= DLP/TP
            DDTPM=-DLP/T
         end if
         CALL CONVEC(ID,T0,P0,PG0,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,PG0,PR0,AB0,DLT,FLXC1,VCON)
         DHCDT0=(FLXC1-FLXCNV)*1.D3*HALF
         DHCDTP=DHCDT0/TP
         DHCDTU=DHCDT0/T
         DHCDPP=0.
         IF(IPRESS.GT.0) THEN
            PG1=1.001*PG0
            CALL CONVEC(ID,T0,P0,PG1,PR0,AB0,DLT,FLXC2,VCON)
            DHCDPP=(FLXC2-FLXCNV)*1.D3/PG0*HALF
            IF(IPRESS.GT.1) THEN
               P1=1.001*P0
               CALL CONVEC(ID,T0,P1,PG0,PR0,AB0,DLT,FLXC3,VCON)
               DHCDPT=(FLXC3-FLXCNV)*1.D3/P0*HALF
               DHCDPP=DHCDPP+DHCDPT
               DHCDTP=DHCDTP+DHCDPT*4.*PR0/T0
            END IF
         END IF
         IF(INDL.EQ.0) THEN
            DHCDTP=DHCDTP+DHCDDP*DDTP0
            DHCDTU=DHCDTU+DHCDDP*DDTPM
         END IF
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
         DELHC=WMM(ID)/DELM*(FLXCNV-FLXC(ID))
         IF(ICONV.GT.0) THEN
         IF(INHE.GT.0) THEN
         A(NRE,NHE)=A(NRE,NHE)-RDELM*DHCDP*BOLK*TM*reint(id)
         B(NRE,NHE)=B(NRE,NHE)+
     *             (DELHC+RDELM*(DHCDP-DHCDPP)*BOLK*T)*reint(id)
         C(NRE,NHE)=C(NRE,NHE)+RDELM*DHCDPP*BOLK*TP*reint(id)
         END IF
         A(NRE,NRE)=A(NRE,NRE)-RDELM*(DHCDP*PGM/TM+DHCDTM)*reint(id)
         B(NRE,NRE)=B(NRE,NRE)+
     *              RDELM*((DHCDPP-DHCDP)*PG/T+DHCDTP-DHCDTU)*reint(id)
         C(NRE,NRE)=C(NRE,NRE)+RDELM*(DHCDPP*PGP/TP+DHCDTP)*reint(id)
         IF(INPC.GT.0) B(NRE,NPC)=B(NRE,NPC)-DELHC*reint(id)
         IF(INDL.GT.0) THEN
         B(NRE,NDEL)=B(NRE,NDEL)-RDELM*DHCDD*reint(id)
         C(NRE,NDEL)=C(NRE,NDEL)+RDELM*DHCDDP*reint(id)
         END IF
         END IF
         VECL(NRE)=VECL(NRE)-RDELM*(FLXCNV-FLXC(ID))*reint(id)
         end if
      END IF
      RETURN
      END