am 5/29/07
rerun the 30 cases, because the way how the ions in
subroutine elden is calculated is not consistent

the minimum electron density is set before the ions 
(n2p,o2p,nop) are calculated. The sum of all ions is used
to determine the conductivity. Since the Ne is modified
after solving for Ne, and used in the equations for the ions,
the sum of the ions is not equal to the Ne any more.
This effected the vertical drift- peak in the early morning
at 4 LT.

A quick fix: calculate a factor if the Ne < Ne_min
calculate the ions with the original Ne, but then muliply
it by the factor. Also Op is multiplied by the factor
although it's not calculated there. 
Basically we multilpy all the equations by the factor.
 
We should change the production rate in th long run- revisit the problem 
SOON! We don't know how to distribute the production rate to the
different ions.

source code now: line 189 elden.F

! Calculate N2+, O2+, NO+ (cm3)
!
      do i=lon0,lon1      
        do k=lev0,lev1-1
          if (root(k,i) < 1.) root(k,i) = 1.0  ! insure positive Ne from solver
	                                       ! in case there is a problem
          n2p(k,i) = d(k,i)/(e(k,i)+ra3(k,i,lat)*root(k,i))    
          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))
!
! nop = (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))
          if (root(k,i) < elden_min) then
	      fac = elden_min/root(k,i)
	      root(k,i)= elden_min
              n2p(k,i) = fac*n2p(k,i)
              o2p(k,i) = fac*o2p(k,i)
              nop(k,i) = fac*nop(k,i)
	      op(k,i)  = fac*op(k,i)  ! modify the O+ at the current timestep
	  endif    
        enddo ! k=lev0,lev1-1
      enddo ! i=lon0,lon1      


before: line 177 elden.F
!
! Insure positive Ne (at least elden_min):
      do i=lon0,lon1
        do k=lev0,lev1-1
          if (root(k,i) < elden_min) root(k,i) = elden_min
        enddo 
      enddo
!
! Calculate N2+, O2+, NO+ (cm3)
!
      do i=lon0,lon1      
        do k=lev0,lev1-1
          n2p(k,i) = d(k,i)/(e(k,i)+ra3(k,i,lat)*root(k,i))    
          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))
!
! nop = (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