! ! 5/29/98 btf: Original code from Marina Galand, here modified for ! use in tiegcm (tgcm13) ! 9/9/99 btf: added implicit none for tgcm13mt. ! subroutine low_proton(alp,eflux,qia,lat) implicit none ! ! Determination of electron and ion production rates ! - noted QTI and QIA respectively - ! induced by a proton beam of an average energy of ! ALP keV and an incident energy flux of EFLUX erg.cm-2.s-1. ! Restriction : 2 < ALP < 40 keV ! Note : ALP = 2*E0 (for a Maxwellian) ! include 'params.h' include 'blnk.h' include 'cons.h' include 'buff.h' include 'index.h' ! ! Args: real,intent(in) :: eflux(zimxp),alp(zimxp) integer :: lat ! ! qia is output ion production for N2+, O2+, O+, NO+, N+ real,intent(out) :: qia(zimxp,zkmxp,5) ! ! Locals: integer :: i,k,n,m,ntk,npo2k,npo1k,nzk,ii real :: ALPV(zimxp),EL,KL,RANGE,F_E,R_norm(zimxp),LAMBDA,ETA, + W(6),shape(3),Kd(4),E,QIA_TOT real :: + xnn2(zimxp,zkmxp), ! number density of n2 + xno2(zimxp,zkmxp), ! number density of o2 + xno(zimxp,zkmxp) ! number density of o real :: rho(zimxp,zkmxp),zpht(zimxp,zkmxp) real :: pkt_barm(zimxp,zkmxp) ! ! qti = output electron production, not currently used in tiegcm: real :: qti(zimxp,zkmxp) ! ! Parameters real,parameter :: + EminV=100., + Mean_mass=4.23e-23, + tab_shape(9)=(/1.,0.,0., 0.,0.1,1., 0.,0.7,1./), + tab_EL(6)=(/0.87, 0.84, 0.883, 0.9, 0.92, 0.92/), + tab_KL(6)=(/4.5e-18,8.1e-18,7.0e-18,7.0e-18,7.0e-18,8.e-18/), + tab_E0(6)=(/13.5, 6.5, 3.5, 2.5, 1.5, 0./), + P_W1(7)=(/37.90, 2.95e-1, 2.03e-1, -1.57e-2, 3.44e-4, 0., 0./), + P_W2(7)=(/156.68, -63.86, 17.33, -2.32, 1.641e-1, -5.83e-3, + 8.21e-5/), + P_W3(7)=(/8.01e-1, 2.09, -2.91e-2, 1.58e-4, 0., 0., 0./), + P_W4(7)=(/14.26, 2.73, -8.55e-2, 2.93e-3, 0., 0., 0./), + P_W5(7)=(/28.81, 32.69, -1.83, 4.02e-2, 0., 0., 0./), + P_W6(7)=(/18.94, 1.32, -8.25e-2, 3.80e-3, -7.06e-5, + 0., 0./), + P_Kd1(3)=(/3.85e-2, 3.03e-2, -4.74e-4/), + P_Kd2(3)=(/4.56e-2, 1.27e-2, -1.15e-4/), + P_Kd3(3)=(/8.90e-2, 3.95e-2, -4.21e-4/), + P_Kd4(3)=(/3.15e-1, 2.19e-2, -3.92e-4/), + E_Kd(4)=(/30., 50., 40., 30./) ! ! exp(-z)*p0 1. ! pkt_barm = ---------- * ------------------ ! 1.38e-16*tn o2/32.+o/16.+n2/28. ! ! Where: o2, o, n2 are mass mixing ratios. ! Used in number densities calculation below. ! ntk = nj+nt-1 ! tn npo2k = nj+nps-1 ! o2 (mmr) npo1k = nj+nps2-1 ! o (mmr) do k=1,kmax ntk = ntk+1 npo2k = npo2k+1 npo1k = npo1k+1 do i=1,len1 pkt_barm(i,k) = exps(k)*c(81) / (c(84)*f(i,ntk)* + (f(i,npo2k)/rmass(1)+f(i,npo1k)/rmass(2)+ + (1.-f(i,npo2k)-f(i,npo1k))/rmass(3))) xnn2(i,k) = max(.00001,(1.-f(i,npo1k)-f(i,npo2k))) ! n2 mmr enddo enddo ! ! Calculate number densities of O2, O, and N2 from mass mixing ! ratios in f-array. Store in locals xno2, xno, xnn2. ! Also sum total number density rho = o2+o+n2 (gm/cm3) ! npo2k = nj+nps ! o2 (mmr) npo1k = nj+nps2 ! o (mmr) nzk = nj+nz ! geopotential ht (cm) do i=1,len2 rho(i,1) = (f(i,npo2k)+f(i,npo1k)+xnn2(i,1))*pkt_barm(i,1)* + 1.66e-24 xno2(i,1) = f(i,npo2k)*pkt_barm(i,1)/rmass(1) xno(i,1) = f(i,npo1k)*pkt_barm(i,1)/rmass(2) xnn2(i,1) = xnn2(i,1) *pkt_barm(i,1)/rmass(3) zpht(i,1) = f(i,nzk) enddo ! ! Loop down from top: do k=kmax,1,-1 ! ! Longitude loop: do i=1,len1 ! ! RANGE (g.cm-2) ALPV(i)=ALP(i)*1.e3 do ii=6,1,-1 if (ALP(i).gt.2.*tab_E0(ii)) then EL=tab_EL(ii) KL=tab_KL(ii) endif enddo RANGE = Mean_mass/KL/(1.-EL)*(ALPV(i)**(1.-EL)-EminV**(1.-EL)) F_E = ALPV(i)**(1.-EL)/KL/(1.-EL)/RANGE*Mean_mass ! ! R_norm : Normalized atmospheric scattering depth (g.cm-2) if (k.eq.kmax) then R_norm(i)=0. else R_norm(i)=R_norm(i)+(RHO(i,k+1)+RHO(i,k))/2* + (zpht(i,k+1)-zpht(i,k))/RANGE endif ! ! LAMBDA : Normalized energy deposition function if (abs(1.-EL)/(1.-EL)*(F_E-R_norm(i)).gt.0.) then LAMBDA=KL*RANGE/Mean_mass/(ALPV(i)-EminV)* + (((1.-EL)*KL*RANGE/Mean_mass* + (F_E-R_norm(i)))**(EL/(1.-EL))) else LAMBDA=0. endif ! ! ETA : Energy deposition rate (eV.cm-2.s-1) ETA=EFLUX(i)*1.e-7/(1.602e-19)/RANGE*RHO(i,K)*LAMBDA ! do ii=1,6 W(ii)=0. enddo do ii=1,7 W(1)=W(1)+P_W1(ii)*(ALP(i)/2)**(ii-1) W(2)=W(2)+P_W2(ii)*(ALP(i)/2)**(ii-1) W(3)=W(3)+P_W3(ii)*(ALP(i)/2)**(ii-1) W(4)=W(4)+P_W4(ii)*(ALP(i)/2)**(ii-1) W(5)=W(5)+P_W5(ii)*(ALP(i)/2)**(ii-1) W(6)=W(6)+P_W6(ii)*(ALP(i)/2)**(ii-1) enddo ! ! QTI : Electron production rate (cm-3.s-1) QTI(i,K)=ETA/W(6) ! ! QIA(i,K,n) : Ion production rate (cm-3.s-1) ! n=1,2,3,4,5 -> N2+,O2+,O+,NO+,N+ ! shape(1)=tab_shape(1)*xnn2(i,k)/(tab_shape(1)*xnn2(i,k)+ + tab_shape(2)*xno2(i,k)+tab_shape(3)*xno(i,k)) shape(2)=tab_shape(5)*xno2(i,k)/(tab_shape(4)*xnn2(i,k)+ + tab_shape(5)*xno2(i,k)+tab_shape(6)*xno(i,k)) shape(3)=tab_shape(9)*xno(i,k)/(tab_shape(7)*xnn2(i,k)+ + tab_shape(8)*xno2(i,k)+tab_shape(9)*xno(i,k)) ! do n=1,4 Kd(n)=0. if (ALP(i).gt.E_Kd(n)) then E=E_Kd(n) else E=ALP(i) endif select case (n) case (1) do m=1,3 Kd(n)=Kd(n)+P_Kd1(m)*E**(m-1) enddo case (2) do m=1,3 Kd(n)=Kd(n)+P_Kd2(m)*E**(m-1) enddo case (3) do m=1,3 Kd(n)=Kd(n)+P_Kd3(m)*E**(m-1) enddo case (4) do m=1,3 Kd(n)=Kd(n)+P_Kd4(m)*E**(m-1) enddo end select enddo ! QIA(i,K,1)=shape(1)*(ETA/W(1)/(1+Kd(1))+ETA/W(2)/(1+Kd(2))) QIA(i,K,5)=shape(1)*(ETA/W(1)*Kd(1)/(1+Kd(1))+ + ETA/W(2)*Kd(2)/(1+Kd(2))) QIA(i,K,2)=shape(2)*(ETA/W(3)/(1+Kd(3))+ETA/W(4)/(1+Kd(4))) QIA(i,K,3)=shape(2)*(ETA/W(3)*Kd(3)/(1+Kd(3))+ + ETA/W(4)*Kd(4)/(1+Kd(4)))+shape(3)*ETA/W(5) QIA_TOT=QIA(i,k,1)+QIA(i,k,2)+QIA(i,k,3)+QIA(i,k,5) if (QIA_TOT.gt.1.e-20) then do n=1,3 QIA(i,k,n)=QIA(i,k,n)*QTI(i,k)/QIA_TOT enddo QIA(i,k,5)=QIA(i,k,5)*QTI(i,k)/QIA_TOT endif ! write(6,"(1x,i5,2x,6e12.3)") K,i,zpht(i,K),QIA(i,K,1), ! + QIA(i,k,2),QIA(i,k,3),QIA(i,k,5),QTI(i,k) enddo enddo return end