SUBROUTINE RUSSEL(TEM,PG)
c     =========================
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      COMMON/COMFH1/C(600,5),PPMOL(600),APMLOG(600),
     *              XIP(100),XIP2(100),CCOMP(100),UIIDUI(100),
     *              P(100),FP(100),XKP(100),XK2(100),EPS,SWITER,
     *              NELEM(5,600),NATO(5,600),MMAX(600),
     *              NELEMX(100),NMETAL,NMOLEC,NIMAX
        DIMENSION FX(100),DFX(100),Z(100),PREV(100),WA(100),
     *            UIIDU2(100)
C
C       ECONST=4.342945E-1
        ECONST=4.3426E-1
        XKCON=6.667343E-1
        EPSDIE=5.0E-5
        T=5040.4/TEM
        PGLOG=log10(PG)
        tk=1./(tem*1.38054e-16)
C
C    HEH=helium/hydrogen ratio by number
C
c       HEH=CCOMP(2)/CCOMP(1)
        HEH=YTOT(1)-UN
C
C    evaluation of log XKP(MOL)
C
        DO J=1,NMOLEC
           APLOGJ=C(J,5)
           DO K=1,4
              KM5=5-K
              APLOGJ=APLOGJ*T + C(J,KM5)
           END DO
           APMLOG(J)=APLOGJ
        END DO
        apmlog(1)=-log10(1.0353e-16/tem/sqrt(tem)*tk*exp(8762.9/tem))
        DHH=(((0.1196952E-02*T-0.2125713E-01)*T+0.1545253E+00)*T
     *     -0.5161452E+01)*T+0.1277356E+02
        DHH=EXP(DHH/ECONST)
C
C  evaluation of the ionization constants
C
        TEM25=TEM**2*SQRT(TEM)
        DO I=1,NMETAL
           NELEMI = NELEMX(I)
*
* calculation of the partition functions following Irwin (1981)
C
           call mpartf(nelemi,1,0,tem,g0,dulog)
           call mpartf(nelemi,2,0,tem,g1,dulog)
           call mpartf(nelemi,3,0,tem,g2,dulog)
           uiidui(nelemi)=g1/g0*xkcon
           uiidu2(nelemi)=g2/g1*xkcon
           uiidui(nelemi)=g1/g0*xkcon
           XKP(NELEMI)=UIIDUI(NELEMI)*TEM25*
     *                 EXP(-XIP(NELEMI)*T/ECONST)
           XK2(NELEMI)=UIIDU2(NELEMI)*TEM25*
     *                 EXP(-XIP2(NELEMI)*T/ECONST)
           xk2(nelemi)=max(xk2(nelemi),1.d-70)
        END DO
        XK2(1)=0.
C
C   preliminary value of PH at high temperatures
C
        HKP=XKP(1)
        IF(T.LT.0.6) THEN
           PPH=SQRT(HKP*(PG/(1.0+HEH)+HKP))-HKP
           PH=PPH**2/HKP
         ELSE
           IF(PG/DHH.LE.0.1) THEN
              PH=PG/(1.0+HEH)
            ELSE
              PH=0.5 * (SQRT(DHH*(DHH+4.0 *PG/(1.0+HEH)))-DHH)
           END IF
        END IF
C
C  evaluation of the fictitious pressures of hydrogen
C     PG=PH+PHH+2.0*PPH+HEH*(PH+2.0*PHH+PPH)
C
        U=(1.0+2.0*HEH)/DHH
        Q=1.0+HEH
        R=(2.0+HEH)*SQRT(HKP)
        S=-1.0*PG
        X=SQRT(PH)
C
C       Russell iterations
C
        ITERAT=0
   10   CONTINUE
        F=((U*X**2+Q)*X+R)*X+S
        DF=2.0*(2.0*U*X**2+Q)*X+R
        XR=X-F/DF
C
        IF(ABS((X-XR)/XR).GT.EPSDIE) THEN
           ITERAT=ITERAT+1
           IF(ITERAT.GT.50) THEN
              WRITE(6,710) TEM,PG,X,XR,PH
  710    FORMAT(1H1, ' NOT CONVERGE IN RUSSEL '/// 'TEM=',F9.2,5X,'PG=',
     *        E12.5,5X,'X1=',E12.5,5X,'X2=',E12.5,5X,'PH=',E12.5/////)
            ELSE
              X=XR
              GO TO 10
           END IF 
        END IF
        PH=XR**2
        PHH=PH**2/DHH
        PPH=SQRT(HKP*PH)
        FPH=PH+2.0*PHH+PPH
        P(100)=PPH
C
C   evaluation of the fictitious pressure of each element
C
        DO I=1,NMETAL
           NELEMI=NELEMX(I)
           FP(NELEMI)=CCOMP(NELEMI)*FPH
        END DO
C
C   check of initialization
C
        PE=P(99)
C
C    Russell equations
C
        NITERR = 0
   20   CONTINUE
        DO I=1,NMETAL
           NELEMI=NELEMX(I)
c          FX(NELEMI)=-FP(NELEMI)+P(NELEMI)*(1.0+XKP(NELEMI)/PE)
c          DFX(NELEMI)=1.0+XKP(NELEMI)/PE
           DFX(NELEMI)=1.0+XKP(NELEMI)/PE*(1.0+XK2(NELEMI)/PE)
           FX(NELEMI)=-FP(NELEMI)+P(NELEMI)*DFX(NELEMI)
        END DO
C
        SPNION=0.0
        spnplu=0.
        DO J=1,NMOLEC
           MMAXJ=MMAX(J)
           PMOLJL=-APMLOG(J)
           DO M=1,MMAXJ
              NELEMJ=NELEM(M,J)
              NATOMJ=NATO(M,J)
              PMOLJL=PMOLJL+DFLOAT(NATOMJ)*log10(P(NELEMJ))
           END DO
C
           PMOLJ=EXP(PMOLJL/ECONST)
           DO M=1,MMAXJ
              NELEMJ=NELEM(M,J)
              NATOMJ=NATO(M,J)
              ATOMJ=DFLOAT(NATOMJ)
              IF(NELEMJ.EQ.99) then
                 if(natomj.ge.0) then
                    SPNION=SPNION+PMOLJ*NATOMJ
                  else
                    SPNPLU=SPNPLU-PMOLJ*NATOMJ
                 end if
              end if
              DO I=1,NMETAL
                 NELEMI=NELEMX(I)
                 IF(NELEMJ.EQ.NELEMI) THEN
                    FX(NELEMI)=FX(NELEMI)+ATOMJ*PMOLJ
                    DFX(NELEMI)=DFX(NELEMI)+ATOMJ**2*
     *                          PMOLJ/P(NELEMI)
                 END IF
              END DO
           END DO
           PPMOL(J)=PMOLJ
        END DO
C
C   solution of the Russell equations by Newton-Raphson method
C
        DO I=1,NMETAL
           NELEMI=NELEMX(I)
           WA(I)=log10(P(NELEMI)+1.0D-70)
        END DO
        IMAXP1=NMETAL+1
        WA(IMAXP1)=log10(PE+1.0D-70)                  
        DELTRS = 0.0 
        DO I=1,NMETAL
           NELEMI=NELEMX(I)
           PREV(NELEMI)=P(NELEMI)-FX(NELEMI)/DFX(NELEMI)
           PREV(NELEMI)=ABS(PREV(NELEMI))
           IF(PREV(NELEMI).LT.1.0D-70) PREV(NELEMI)=1.0D-70
           Z(NELEMI)=PREV(NELEMI)/P(NELEMI)
           DELTRS=DELTRS+ABS(Z(NELEMI)-1.0)
           IF(SWITER.GT.0.0) THEN
              P(NELEMI)=(PREV(NELEMI)+P(NELEMI))*0.5
            ELSE
              P(NELEMI)=PREV(NELEMI)
           END IF
        END DO
C
C   ionization equilibrium
C
        PEREV =0.0
        DO I=1,NMETAL
           NELEMI = NELEMX(I)
           PEREV=PEREV+XKP(NELEMI)*P(NELEMI)*(1.+xk2(nelemi)/pe)
        END DO
C
        PEREV=SQRT(PEREV/(1.0+SPNION/PE))
        DELTRS=DELTRS+ABS((PE-PEREV)/PE)
        PE=(PEREV+PE)*0.5
        P(99)=PE
        IF(DELTRS.GT.EPS) THEN
           NITERR=NITERR+1
           IF(NITERR.LE.NIMAX) THEN
              GO TO 20
            ELSE
              WRITE(6,605) NIMAX
           END IF
        END IF
  605   FORMAT(1H0,'*DOES NOT CONVERGE AFTER ',I4,' ITERATIONS')
C
        RETURN
        END