!
      subroutine elden(tn,barm,op,o2,o1,ox,n2,n2d,no,n4s,co2,xiop2p,
     |  xiop2d,z,nplus,n2p,nop,o2p,electrons,lev0,lev1,lon0,lon1,lat)
!
! Solve for electron density.
! Also calculate ions N2+ (n2p), NO+ (nop), and O2+ (o2p).
!
      use cons_module,only: p0,expz,boltz,rmassinv_o2,rmassinv_n2d,
     |  rmassinv_o1,rmassinv_n2,rmassinv_no,rmassinv_n4s,rmassinv_co2
      use chemrates_module,only: rk1,rk2,rk3,rk4,rk5,rk6,rk7,rk8,
     |  rk9,rk13,rk16,rk23,rk24,rk27,ra1,ra2,ra3
      use qrj_module,only: qnp,qnop,qo2p,qn2p
      use addfld_module,only: addfld
      use diags_module,only: mkdiag_TEC,mkdiag_HNMF2
      implicit none
!
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
!
! Input args: press vs longitude input fields (2d (k,i)):
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: 
     |  tn,    ! neutral temperature (deg K)
     |  barm,  ! mean molecular weight
     |  op,    ! O+ ion (current time-step)
     |  o2,    ! molecular oxygen (mmr)
     |  o1,    ! atomic oxygen (mmr)
     |  ox,    ! atomic oxygen (mmr)
     |  n2,    ! molecular nitrogen (mmr)
     |  n2d,   ! n2d
     |  no,    ! nitric oxide
     |  n4s,   ! n4s
     |  co2,   ! co2
     |  xiop2p,! O+(2P) from oplus
     |  xiop2d,! O+(2D) from oplus
     |  nplus, ! N+ from oplus
     |  z      ! geopotential (needed by mkdiag_TEC)
!
! Output args (particles/cm3):
      real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: 
     |  n2p,       ! N2+ output
     |  nop,       ! NO+ output
     |  o2p,       ! O2+ output
     |  electrons  ! electron density (output to f4d(ne))
!
! Local:
      integer :: k,i
      integer ::i0,i1
      real,dimension(lev0:lev1,lon0:lon1) :: 
     |  xnmbarm,        ! for conversion from mmr to cm3
     |  a0,a1,a2,a3,a4, ! coefficients for quartic solver
     |  a,b,c,d,e,fg,h, ! terms for quartic coefficients
     |  root            ! output from quartic solver
      real,dimension(lon0:lon1) ::
     |  tec             ! diagnostic total electron content
!
! Convenience ints for addfld calls:
      i0 = lon0
      i1 = lon1
!
! Save inputs to secondary histories:
!     call addfld('TN_ELD' ,' ',' ',tn    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('BAR_ELD',' ',' ',barm  (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('OP_ELD' ,' ',' ',op    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O2_ELD' ,' ',' ',o2    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O1_ELD' ,' ',' ',o1    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('N2D_ELD',' ',' ',n2d   (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NO_ELD' ,' ',' ',no    (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('N4S_ELD',' ',' ',n4s   (:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('XIOP2P' ,' ',' ',xiop2p(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('XIOP2D' ,' ',' ',xiop2d(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
      do i=lon0,lon1
        do k=lev0,lev1-1
          xnmbarm(k,i) = p0*expz(k)*.5*(barm(k,i)+barm(k+1,i))/
     |      (boltz*tn(k,i))
!
! Set up terms for quartic coefficients:
!
! A = QI(NO+)+K2*N(O+)*N(N2)+K7*N(N+)*N(02)+B9*N(NO) (s10)
!
          a(k,i) = .5*(qnop(k,i,lat)+qnop(k+1,i,lat))+xnmbarm(k,i)*
     |      (rk2(k,i,lat)*op(k,i)*n2(k,i)*rmassinv_n2+
     |       rk7*nplus(k,i)*o2(k,i)*rmassinv_o2)
!
! B = QI(O2+)+K1*N(O+)*N(O2)+K6*N(N+)*N(02) (s9)
!
          b(k,i) = .5*(qo2p(k,i,lat)+qo2p(k+1,i,lat))+xnmbarm(k,i)*
     |      ((rk1(k,i,lat)*op(k,i)+rk6*nplus(k,i))*o2(k,i)*
     |       rmassinv_o2+rk27*xiop2d(k,i)*o2(k,i)*rmassinv_o2+
     |       rk13*op(k,i)*co2(k,i)*rmassinv_co2+rk23*xiop2p(k,i)*
     |       o2(k,i)*rmassinv_o2)
!
! C = K4*N(N4S)+K5*N(NO) (s8)
!
          c(k,i) = xnmbarm(k,i)*(rk4*n4s(k,i)*rmassinv_n4s+
     |       rk5*no(k,i)*rmassinv_no)
!
! D = QI(N2+) (s7)
! (Note: xnmbar(k,i) at end is missing in tgcm15/tiegcm1)
!
          d(k,i) = .5*(qn2p(k,i,lat)+qn2p(k+1,i,lat))+(rk16*xiop2p(k,i)+
     |      rk24*xiop2d(k,i))*n2(k,i)*rmassinv_n2*xnmbarm(k,i)
!
! E =K3*N(O)+K9*N(O2) (s6)
!
          e(k,i) = xnmbarm(k,i)*(rk3(k,i,lat)*o1(k,i)*rmassinv_o1+
     |       rk9*o2(k,i)*rmassinv_o2)
!
! F+G = N(O+)+N(N+) (s5)
!
          fg(k,i) = op(k,i)+nplus(k,i)
!
! H = K9*N(02) (s4)
!
          h(k,i) = xnmbarm(k,i)*rk9*o2(k,i)*rmassinv_o2
!
          if (a(k,i) < 1.e-20) a(k,i) = 1.e-20    ! s10
          if (b(k,i) < 1.e-20) b(k,i) = 1.e-20    ! s9
          if (c(k,i) < 1.e-20) c(k,i) = 1.e-20    ! s8
          if (d(k,i) < 1.e-20) d(k,i) = 1.e-20    ! s7
          if (e(k,i) < 1.e-20) e(k,i) = 1.e-20    ! s6
          if (fg(k,i) < 1.e-20) fg(k,i) = 1.e-20  ! s5
          if (h(k,i) < 1.e-20) h(k,i) = 1.e-20    ! s4
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Coefficients for quartic solver: a0,a1,a2,a3,a4
!
! A0 = -e*c*(a+b+d) (s15)
!
          a0(k,i) = -e(k,i) * c(k,i) * (a(k,i) + b(k,i) + d(k,i))
!
! A1 = -(ra1*(e*(c*fg+b)+d*(c+h))+ra2*(e*(a+d)-h*d)+ra3*c*(a+b))/4. (s14)
!
          a1(k,i) = -(ra1(k,i,lat)*(e(k,i)*(c(k,i)*fg(k,i)+b(k,i))+
     |      d(k,i)*(c(k,i)+h(k,i)))+ra2(k,i,lat)*(e(k,i)*(a(k,i)+
     |      d(k,i))-h(k,i)*d(k,i))+ra3(k,i,lat)*c(k,i)*(a(k,i)+b(k,i)))/
     |      4.
!
! A2 = (ra1*(e*c-(ra2*e+ra3*c)*fg-ra2*d-ra3*b)-ra2*ra3*a)/6. (s13)
!
          a2(k,i) = (ra1(k,i,lat)*(e(k,i)*c(k,i)-(ra2(k,i,lat)*e(k,i)+
     |      ra3(k,i,lat)*c(k,i))*fg(k,i)-ra2(k,i,lat)*d(k,i)-
     |      ra3(k,i,lat)*b(k,i))-ra2(k,i,lat)*ra3(k,i,lat)*a(k,i))/6.
!
! A3 = (ra1*(ra2*e+ra3*c-ra2*ra3*fg))/4. (s12)
!
          a3(k,i) = (ra1(k,i,lat)*(ra2(k,i,lat)*e(k,i)+ra3(k,i,lat)*
     |      c(k,i)-ra2(k,i,lat)*ra3(k,i,lat)*fg(k,i)))/4.
!
! A4 = ra1*ra2*ra3 (s11)
!
          a4(k,i) = ra1(k,i,lat)*ra2(k,i,lat)*ra3(k,i,lat)
!
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1
!
!     call addfld('XNMBARM' ,' ',' ',xnmbarm,
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NPLUS'   ,' ',' ',nplus(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('A_COEF'  ,' ',' ',a(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('B_COEF'  ,' ',' ',b(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('C_COEF'  ,' ',' ',c(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('D_COEF'  ,' ',' ',d(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('E_COEF'  ,' ',' ',e(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('FG_COEF' ,' ',' ',fg(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('H_COEF'  ,' ',' ',h(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('A0'      ,' ',' ',a0(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('A1'      ,' ',' ',a1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('A2'      ,' ',' ',a2(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('A3'      ,' ',' ',a3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!     call addfld('A4'      ,' ',' ',a4(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! Solve quartic. Vquart returns electron density Ne in root:
!
      call vquart(a0,a1,a2,a3,a4,root,lev0,lev1,lon0,lon1,lat)
!
!     call addfld('ROOT',' ',' ',root(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',i0,i1,lat)
!
! Insure positive Ne (at least 3100):
!     where(root < 3.1e3) root = 3.1e3
!     do i=lon0,lon1
!       do k=lev0,lev1-1
!         if (root(k,i) < 3.1e3) root(k,i) = 3.1e3
!       enddo
!     enddo
!
! Calculate N2+, O2+, NO+ (cm3)
! (tgcm24: a=s10, b=s9, c=s8, d=s7, e=s6, fg=s5, h=s4)
!
      do i=lon0,lon1      
        do k=lev0,lev1-1
!
! N2+ = d/(e+ra3*root)
          n2p(k,i) = d(k,i)/(e(k,i)+ra3(k,i,lat)*root(k,i))    
!
! O2+ = (b+h*d/(e+ra3*root))/(c+ra2*root)
          o2p(k,i) = (b(k,i)+h(k,i)*d(k,i)/(e(k,i)+ra3(k,i,lat)*
     |      root(k,i)))/(c(k,i)+ra2(k,i,lat)*root(k,i))
!
! NO+ = (a+d*(e-h)/(e+ra3*root)+c*(b+h*d/(e+ra3*root))/(c+ra2*root))/(ra1*root)
          nop(k,i)=(a(k,i)+d(k,i)*(e(k,i)-h(k,i))/(e(k,i)+ra3(k,i,lat)*
     |     root(k,i))+c(k,i)*(b(k,i)+h(k,i)*d(k,i)/(e(k,i)+ra3(k,i,lat)*
     |     root(k,i)))/(c(k,i)+ra2(k,i,lat)*root(k,i)))/(ra1(k,i,lat)*
     |     root(k,i))
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1      
!
!     call addfld('N2P_ELD',' ',' ',n2p(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('O2P_ELD',' ',' ',o2p(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!     call addfld('NOP_ELD',' ',' ',nop(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
! Top level of O2+
      do i=lon0,lon1
        o2p(lev1,i) = o2p(lev1-1,i)
      enddo
!
! Transfer root to electrons output array:
      do i=lon0,lon1
        do k=lev0,lev1-2
          electrons(k+1,i) = sqrt(root(k,i)*root(k+1,i))
        enddo ! k=lev0,lev1-2
!
! Lower and upper boundaries:
        electrons(lev0,i) = sqrt(root(lev0  ,i)**3/root(lev0+1,i))
        electrons(lev1,i) = sqrt(root(lev1-1,i)**3/root(lev1-2,i))
      enddo ! i=lon0,lon1
!
!     call addfld('NE_ELDEN',' ',' ',electrons(:,i0:i1),
!    |  'lev',lev0,lev1,'lon',i0,i1,lat)
!
!
! Save diagnostic tec (total electron content):
      call mkdiag_TEC('TEC',tec,z,electrons,lev0,lev1,lon0,lon1,lat)
!
! Save diagnostics hmf2, nmf2:
      call mkdiag_HNMF2('HMF2',z,electrons,lev0,lev1,lon0,lon1,lat)
      call mkdiag_HNMF2('NMF2',z,electrons,lev0,lev1,lon0,lon1,lat)
      call mkdiag_HNMF2('FOF2',z,electrons,lev0,lev1,lon0,lon1,lat)

      end subroutine elden
!-----------------------------------------------------------------------
      subroutine vquart(a0,a1,a2,a3,a4,root,lev0,lev1,lon0,lon1,lat)
      use addfld_module,only: addfld
      implicit none
!
! Determines five roots of the equation:
!   a4*x**4 + 4.*a3*x**3 + 6.*a2*x**2 + 4.*a1*x + a0 = 0.
!
! Procedure is specificlly designed for real quartics with real roots
! only one of which is positive.
!
! This is called by elden for electron density.
! 
! Args:
      integer,intent(in) :: lev0,lev1,lon0,lon1,lat
      real,dimension(lev0:lev1,lon0:lon1),intent(in)  :: a0,a1,a2,a3,a4
      real,dimension(lev0:lev1,lon0:lon1),intent(out) :: root
!
! Local:
      integer :: k,i,nlevs
      real,dimension(lev0:lev1,lon0:lon1) :: w1,w2,w3 ! work arrays
      real,parameter :: e=1.e-300 ! largest exponent on ieee is about 307
!
      nlevs = lev1-lev0+1
      do i=lon0,lon1
        do k=lev0,lev1-1
!
! w1 = ch
          w1(k,i) = -(a4(k,i)*a0(k,i)-4.*a3(k,i)*a1(k,i)+3.*a2(k,i)**2)/
     |                12.
!
! w2 = cg
          w2(k,i) = (a4(k,i)*(a2(k,i)*a0(k,i)-a1(k,i)**2)-a3(k,i)*
     |      (a3(k,i)*a0(k,i)-a1(k,i)*a2(k,i))+a2(k,i)*(a3(k,i)*a1(k,i)-
     |      a2(k,i)**2))/4.
!
! root=rlam=-2.*real((.5*(cmplx(cg,0.)+csqrt(cmplx(cg**2+4.
!      *ch**3+e,0.)))+cmplx(e,0.))**(1./3.))
! 
          root(k,i) = -2.*real((.5*(cmplx(w2(k,i),0.)+
     |      csqrt(cmplx(w2(k,i)**2+4.*w1(k,i)**3+e,0.)))+
     |      cmplx(e,0.))**(1./3.))
!
! W1=P=SQRT(A(5)*RLAM+A(4)**2-A(5)*A(3)+E)
!
          w1(k,i) = a4(k,i)*root(k,i)+a3(k,i)**2-a4(k,i)*a2(k,i)+e
          if (w1(k,i) < 0.) w1(k,i) = 0.
          w1(k,i) = sqrt(w1(k,i))
!
! W2=Q=SQRT((2.*RLAM+A(3))**2-A(5)*A(1)+E)
!
          w2(k,i) = sqrt((2.*root(k,i)+a2(k,i))**2-a4(k,i)*a0(k,i)+e)
!
! W3=PQ=2.*A(4)*RLAM+A(4)*A(3)-A(5)*A(2)+E
!
          w3(k,i) = 2.*a3(k,i)*root(k,i)+a3(k,i)*a2(k,i)-a4(k,i)*a1(k,i)
     |      +e
!
!  W1=P=SIGN(P,Q*PQ)
!
          w1(k,i) = sign(w1(k,i),w2(k,i)*w3(k,i))
!
! W3=P-A4
!
          w3(k,i) = w1(k,i)-a3(k,i)
!
! Final evaluation of root:
!
          root(k,i) = (w3(k,i)+sqrt(w3(k,i)**2-a4(k,i)*(a2(k,i)+2.*
     |      root(k,i)-w2(k,i))))/a4(k,i)
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1

!     call addfld('W1'    ,' ',' ',w1(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('W2'    ,' ',' ',w2(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('W3'    ,' ',' ',w3(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)
!     call addfld('VQROOT',' ',' ',root(lev0:lev1-1,:),
!    |  'lev',lev0,lev1-1,'lon',lon0,lon1,lat)

      end subroutine vquart
