SUBROUTINE ELDEN
      implicit none
C     ****
C     ****     CALCULATES ELECTRON DENSITY AND ION NUMBER DENSITIES
C     ****
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      include "crates.h"
      real :: rmn4s,rmn2d,rmno,brn2d,cee
      COMMON/MASS/RMN4S,RMN2D,RMNO,BRN2D,CEE
!
! Local:
      integer :: i,k
      integer :: ntk,nmsk,nnpk,nqnpk,nopk,npn2dk,nps1k,nps2k,nqnopk,
     |  nqo2pk,nqn2pk,npnok,npn4sk,nno2k,nn2pk,no2pk,nnopk,nek
C     ****
C     ****     S12 = N*MBAR =PO*EXP(-Z)*MBAR/(K*T)
C     ****
      NTK=NJ+NT-1
      NMSK=NJ+NMS-1
      DO 2 K=1,KMAX
      NTK=NTK+1
      NMSK=NMSK+1
      DO 2 I=1,LEN1
      S12(I,K)=C(81)*EXPS(K)*.5*(F(I,NMSK)+F(I,NMSK+1))/(C(84)*F(I,NTK))
    2 CONTINUE
C     ****
C     ****     CALCULATE N(N+)
C     ****
      NNPK=NNP
      NQNPK=NQNP
      NOPK=NJ+NOP
      NPN2DK=NJ+NPN2D
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      DO 3 I=1,LEN2
      S11(I,1)=(.5*(F(I,NQNPK)+F(I,NQNPK+1))
     1+RK10(I,1)*F(I,NOPK)*F(I,NPN2DK)*S12(I,1)/
     1RMN2D)/(S12(I,1)*((RK6(I,1)+RK7(I,1))*F(I,NPS1K)/RMASS(1)+RK8(I,1)
     2*F(I,NPS2K)/RMASS(2)))
    3 CONTINUE
      DO 4 I=1,LEN2
      F(I,NNPK)=S11(I,1)
    4 CONTINUE
C     ****
C     ****     CALCULATE A, B, C, D, E, (F+G) AND H FOR QUARTIC
C     ****
      NQNOPK=NQNOP
      NQO2PK=NQO2P
      NQN2PK=NQN2P
      NOPK=NJNP+NOP
      NNPK=NNP
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NPNOK=NJ+NPNO
      NPN4SK=NJ+NPN4S
      NNO2K=NNO2
      DO 5 I=1,LEN2
C     ****
C     ****     S10 = A = QI(NO+)+K2*N(O+)*N(N2)+K7*N(N+)*N(02)+B9*N(NO)
C     ****
      S10(I,1)=.5*(F(I,NQNOPK)+F(I,NQNOPK+1))+
     1S12(I,1)*(RK2(I,1)*F(I,NOPK)*(1.-F(I,NPS1K)-
     1F(I,NPS2K))/RMASS(3)+RK7(I,1)*F(I,NNPK)*F(I,NPS1K)/RMASS(1)+
     2.5*(BETA9(I,1)+BETA9(I+LEN1,1))*F(I,NPNOK)/RMNO)
C     ****
C     ****     S9 = B = QI(O2+)+K1*N(O+)*N(O2)+K6*N(N+)*N(02)
C     ****
      S9(I,1)=.5*(F(I,NQO2PK)+F(I,NQO2PK+1))+
     1S12(I,1)*(RK1(I,1)*F(I,NOPK)+RK6(I,1)*
     2F(I,NNPK))*F(I,NPS1K)/RMASS(1)+RK26(I,1)*
     3XIOP2D(I,1)*F(I,NPS1K)/RMASS(1)
C     ****
C     ****     S8 = C = K4*N(N4S)+K5*N(NO)
C     ****
      S8(I,1)=S12(I,1)*(RK4(I,1)*F(I,NPN4SK)/RMN4S+RK5(I,1)*F(I,NPNOK)/
     1RMNO)
C     ****
C     ****     S7 = D = QI(N2+)
C     ****
      S7(I,1)=.5*(F(I,NQN2PK)+F(I,NQN2PK+1))+(RK16(I,1)*XIOP2P(I,1)
     1         +RK23(I,1)*XIOP2D(I,1))*(1.-F(I,NPS1K)-F(I,NPS2K))
     2         /RMASS(3)
C     ****
C     ****     S6 = E =K3*N(O)+K9*N(O2)
C     ****
      S6(I,1)=S12(I,1)*(RK3(I,1)*F(I,NPS2K)/RMASS(2)+RK9(I,1)*F(I,NPS1K)
     1/RMASS(1))
C     ****
C     ****     S5 = F+G = N(O+)+N(N+)
C     ****
      S5(I,1)=F(I,NOPK)+F(I,NNPK)
C     ****
C     ****     S4 = H = K9*N(02)
C     ****
      S4(I,1)=S12(I,1)*RK9(I,1)*F(I,NPS1K)/RMASS(1)
C     ****
C     ****     CALCULATE QUARTIC COEFFICIENTS
C     ****
C     ****     S15 = A0
C     ****
      S15(I,1)=-S6(I,1)*S8(I,1)*(S10(I,1)+S9(I,1)+S7(I,1))
C     ****
C     ****     S14 = A1/4
C     ****
      S14(I,1)=-(RA1(I,1)*(S6(I,1)*(S8(I,1)*S5(I,1)+S9(I,1))+S7(I,1)*
     1(S8(I,1)+S4(I,1)))+RA2(I,1)*(S6(I,1)*(S10(I,1)+S7(I,1))-S4(I,1)*
     2S7(I,1))+RA3(I,1)*S8(I,1)*(S10(I,1)+S9(I,1)))/4.
C     ****
C     ****     S13 = A2/6
C     ****
      S13(I,1)=(RA1(I,1)*(S6(I,1)*S8(I,1)-(RA2(I,1)*S6(I,1)+RA3(I,1)*
     1S8(I,1))*S5(I,1)-RA2(I,1)*S7(I,1)-RA3(I,1)*S9(I,1))-RA2(I,1)*
     2RA3(I,1)*S10(I,1))/6.
C     ****
C     ****     S12 = A3/4
C     ****
      S12(I,1)=(RA1(I,1)*(RA2(I,1)*S6(I,1)+RA3(I,1)*S8(I,1)-RA2(I,1)*
     1RA3(I,1)*S5(I,1)))/4.
C     ****
C     ****     S11 = A4
C     ****
      S11(I,1)=RA1(I,1)*RA2(I,1)*RA3(I,1)
    5 CONTINUE
C     ****
C     ****     SOLVE QUARTIC    ( S5 =  N(E) )
C     ****
      CALL VQUART(S15,S5,S3,S2,S1,LEN3,1,LEN2)
C     ****
C     ****     ENSURE NON NEGATIVE N(E)
C     ****
C     ****
C     ****     CALCULATE N(N2+), N(02+) AND N(NO+)
C     ****
      NN2PK=NN2P
      NO2PK=NJNP+NO2P
      NNOPK=NNOP
      DO 7 I=1,LEN2
c 5/7/98: ensure non-negative Ne:
C     S5(I,1) = merge(S5(I,1),3.3E+3,S5(I,1)-3.3E+3>=0.)
      S5(I,1) = merge(S5(I,1),3.1E+3,S5(I,1)-3.1E+3>=0.)
C     S5(I,1) = merge(S5(I,1),2.8E+3,S5(I,1)-2.8E+3>=0.)
C     S5(I,1) = merge(S5(I,1),2.5E+3,S5(I,1)-2.5E+3>=0.)
      F(I,NN2PK)=S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1))
      F(I,NO2PK)=(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)*S5(I,1)))/
     1(S8(I,1)+RA2(I,1)*S5(I,1))
      F(I,NNOPK)=(S10(I,1)+S7(I,1)*(S6(I,1)-S4(I,1))/(S6(I,1)+RA3(I,1)*
     1S5(I,1))+S8(I,1)*(S9(I,1)+S4(I,1)*S7(I,1)/(S6(I,1)+RA3(I,1)*
     2S5(I,1)))/(S8(I,1)+RA2(I,1)*S5(I,1)))/(RA1(I,1)*S5(I,1))
    7 CONTINUE
C     ****
C     ****
C     ****     PLACE N(E) IN SLOT IN F-ARRAY
C     ****
      NEK=NJNP+NE
      DO 8 I=1,LEN2-LEN1
      F(I,NEK+1)=SQRT(S5(I,1)*S5(I,2))
    8 CONTINUE
      DO 9 I=1,LEN1
      F(I,NEK)=SQRT(S5(I,1)**3/S5(I,2))
      F(I,NEK+KMAX)=SQRT(S5(I,KMAX)**3/S5(I,KMAXM1))
C     F(I,NEK)= 1.5*S5(I,1)-0.5*S5(I,2)
C     F(I,NEK+KMAX)=1.5*S5(I,KMAX)-0.5*S5(I,KMAXM1)
    9 CONTINUE
      RETURN
      END
C