!
subroutine elden(tn,barm,op,op_upd,o2,o1,n2d,no,n4s,xiop2p,xiop2d,
| z,nplus,n2p,nop,o2p,electrons,lev0,lev1,lon0,lon1,lat)
!
! This software is part of the NCAR TIE-GCM. Use is governed by the
! Open Source Academic Research License Agreement contained in the file
! tiegcmlicense.txt.
!
! Solve for electron density.
! Also calculate ions N+ (nplus), 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
use chemrates_module,only: rk1,rk2,rk3,rk4,rk5,rk6,rk7,rk8,
| rk9,rk10,rk16,rk23,rk26,beta9,ra1,ra2,ra3
use qrj_module,only: qnp,qnop,qo2p,qn2p
use addfld_module
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)
| op_upd,! O+ ion (updated, from sub oplus)
| o2, ! molecular oxygen (mmr)
| o1, ! atomic oxygen (mmr)
| n2d, ! n2d
| no, ! nitric oxide
| n4s, ! n4s
| xiop2p,! from oplus
| xiop2d,! from oplus
| z ! geopotential
!
! Output args (particles/cm3):
real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) ::
| nplus, ! N+ output
| n2p, ! N2+ output
| nop, ! NO+ output
| o2p, ! O2+ output
| electrons ! electron density (output to f4d(ne))
!
! Local:
integer :: k,i,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
| xn2, ! n2 (mmr) (1-o2-o)
| root ! output from quartic solver
real,dimension(lon0:lon1) ::
| tec ! diagnostic total electron content
! write(6,"('enter elden: lat=',i2)") lat
i0 = lon0 ; i1 = lon1
!
! call addfld('TN_ELD' ,' ',' ',tn (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('BAR_ELD',' ',' ',barm (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('OP_ELD' ,' ',' ',op (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('OP1_ELD',' ',' ',op_upd(:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('O2_ELD' ,' ',' ',o2 (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('O1_ELD' ,' ',' ',o1 (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('N2D_ELD',' ',' ',n2d (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('NO_ELD' ,' ',' ',no (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('N4S_ELD',' ',' ',n4s (:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('XIOP2P' ,' ',' ',xiop2p(:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,lat)
! call addfld('XIOP2D' ,' ',' ',xiop2d(:,i0:i1),
! | 'lev',lev0,lev1,'lon',lon0,lon1,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))
xn2(k,i) = (1.-o2(k,i)-o1(k,i)) ! n2 (mmr)
!
! N+ ion (output):
nplus(k,i) = (0.5*(qnp(k,i,lat)+qnp(k+1,i,lat))+
| rk10*op(k,i)*n2d(k,i)*xnmbarm(k,i)*rmassinv_n2d) /
| (xnmbarm(k,i)*((rk6+rk7)*o2(k,i)*rmassinv_o2+
| rk8*o1(k,i)*rmassinv_o1))
!
! 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_upd(k,i)*xn2(k,i)*rmassinv_n2+
| rk7*nplus(k,i)*o2(k,i)*rmassinv_o2+
| .5*(beta9(k,i,lat)+beta9(k+1,i,lat))*no(k,i)*rmassinv_no)
!
! B = QI(O2+)+K1*N(O+)*N(O2)+K6*N(N+)*N(02) (s9)
! (very small "diamond diffs" with tgcm15 due to op_upd)
!
b(k,i) = .5*(qo2p(k,i,lat)+qo2p(k+1,i,lat))+xnmbarm(k,i)*
| (rk1(k,i,lat)*op_upd(k,i)+rk6*nplus(k,i))*o2(k,i)*
| rmassinv_o2+rk26*xiop2d(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)
!
d(k,i) = .5*(qn2p(k,i,lat)+qn2p(k+1,i,lat))+(rk16*xiop2p(k,i)+
| rk23*xiop2d(k,i))*xn2(k,i)*rmassinv_n2
!
! 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)
! (very small "diamond diffs" with tgcm15 due to op_upd)
!
fg(k,i) = op_upd(k,i)+nplus(k,i)
!
! H = K9*N(02)
!
h(k,i) = xnmbarm(k,i)*rk9*o2(k,i)*rmassinv_o2
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! 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,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('B_COEF' ,' ',' ',b,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('C_COEF' ,' ',' ',c,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('D_COEF' ,' ',' ',d,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('E_COEF' ,' ',' ',e,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('FG_COEF' ,' ',' ',fg,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('H_COEF' ,' ',' ',h,'lev',lev0,lev1,'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,'lev',lev0,lev1,'lon',i0,i1,lat)
!
! 1/24/08 btf, maute: Minimum Ne is replaced by new values for flux
! parameter al in qinite.F
!
! 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)
!
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))
enddo ! k=lev0,lev1-1
enddo ! i=lon0,lon1
!
! call addfld('NPLUSb' ,' ',' ',nplus(:,i0:i1),
! | 'lev',lev0,lev1,'lon',i0,i1,lat)
! 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)
!
! 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)
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,i0,i1
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
i0 = lon0 ; i1 = lon1
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,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('W2' ,' ',' ',w2,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('W3' ,' ',' ',w3,'lev',lev0,lev1,'lon',i0,i1,lat)
! call addfld('VQROOT',' ',' ',root,'lev',lev0,lev1,'lon',i0,i1,lat)
end subroutine vquart