      SUBROUTINE ORORA
C     ****
C     ****     CALCULATES AURORAL ADDITIONS TO IONIZATION RATES
C     ****     ON ENTRY:
C     ****     T3=DRIZL                 (PROFILE FOR AURORAL DRIZZLE)
C     ****     T5=CUSP                  (PROFILE FOR AURORAL CUSP)
C     ****     T6,T7=ALFA(1), ALFA(2)   (CHARACTERISTIC PARTICLE ENERGIE
C     ****     T8,T9=F(1),F(2)          (PARTICLE FLUXES)
C     ****
      include "params.h"
      include "blnk.h"
      include "vscr.h"
      include "cons.h"
      include "index.h"
      include "buff.h"
      include "phys.h"
      COMMON /OVALR/ RRAD(2),H0,RH,RROTH,E0,REE,RROTE,FC,ALFAC,FD,ALFAD
     1, ALFK,ALF6,ALF21,RROT6,RROT21,RD6,RD6V,RD21,RH6,RH21,RT6,RT21
     2, ALFA0,RALFA,ALFA20,RALFA2,E20,RE2
c
c Low energy auroral proton input (see ALFALP and EFLUXLP from input):
      common/aurlp/ alfa_lp(zimxp),eflux_lp(zimxp)
C     ****
C     ****
C     ****     ARRAYS AUREFF AND BDRIZ DEFINE VERTICAL PROFILES FOR
C     ****     AURORAL HEATING EFFICIENCY AND BACKGROUND DRIZZLE
C     ****     RESPECTIVELY
C     ****
      DIMENSION AUREFF(ZKMXP),BDRIZ(ZKMXP)
      DATA AUREFF
     1/0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,0.55,
     20.50,0.45,0.37,0.30,0.20,0.16,0.13,0.10,0.08,0.07,0.06,0.05,
     30.05,0.05,0.05,0.05/
      DATA BDRIZ
     1/2.17E-4,4.48E-3,3.78E-2,1.55E-1,3.26E-1,3.85E-1,3.50E-1,2.81E-1,
     22.04E-1,1.38E-1,8.74E-2,5.23E-2,3.01E-2,1.68E-2,9.27E-3,5.08E-3,
     32.79E-3,1.55E-3,8.62E-4,4.85E-4,2.76E-4,1.58E-4,9.11E-5,5.28E-5,
     43.08E-5,1.79E-5,1.04E-5,6.01E-6,3.49E-6/
!
! Ion production from low energy proton source (see low_proton)
      real :: qia(zimxp,zkmxp,5)
!     real,parameter :: eflux=0.4, alp=10.
C     ****
C     ****     NO CALCULATION IF ISHUNK=0
C     ****
      IF(ISHUNK.EQ.0) then
        qia = 1.e-20
!       call addfsech('QIA1_N2+',qia(:,:,1),zimxp,zkmxp,zkmx,j)
!       call addfsech('QIA2_O2+',qia(:,:,2),zimxp,zkmxp,zkmx,j)
!       call addfsech('QIA3_O+ ',qia(:,:,3),zimxp,zkmxp,zkmx,j)
!!      call addfsech('QIA4_NO+',qia(:,:,4),zimxp,zkmxp,zkmx,j)
!       call addfsech('QIA5_N+ ',qia(:,:,5),zimxp,zkmxp,zkmx,j)
        RETURN
      endif
C     ****
C     ****     CALCULATE TOTAL IONIZATION, QI
C     ****
      DO 1 K=1,KMAX
      DO 1 I=1,LEN1
C     ****
C     ****     S7,S8=X(1),X(2)  (AURORA)
C     ****     S9=X  (CUSP)
C     ****     S10=X  (DRIZZLE)
C     ****     WHERE X=(N/4.E-6)**0.606/ALFA
C     ****     AND   N=P0*EXP(-Z)/G
C     ****
      S6(I,K)=(C(81)*EXPS(K)/(C(54)*4.E-6))**0.606
      S7(I,K)=S6(I,K)/T6(I)
      S8(I,K)=S6(I,K)/T7(I)
      S9(I,K)=S6(I,K)/ALFAC
      S10(I,K)=S6(I,K)/ALFAD
    1 CONTINUE
C     ****
C     ****     S3,S4 = F(X(1)),F(X(2)) (AURORA)
C     ****     S5 = F(X) (CUSP)
C     ****     S6 = F(X) (DRIZZLE)
C     ****
      CALL AION(S7,S3,LEN2)
      CALL AION(S8,S4,LEN2)
      CALL AION(S9,S5,LEN2)
      CALL AION(S10,S6,LEN2)
!
! Contribution of low energy protons:
! qia(zimxp,zkmxp,5) is output ion production for N2+, O2+, O+, NO+, N+
! (top level zkmxp not defined)
!
      call low_proton(alfa_lp,eflux_lp,qia,j)
!
! Put qia on secondary histories:
!
!     call addfsech('QIA1_N2+',qia(:,:,1),zimxp,zkmxp,zkmx,j)
!     call addfsech('QIA2_O2+',qia(:,:,2),zimxp,zkmxp,zkmx,j)
!     call addfsech('QIA3_O+ ',qia(:,:,3),zimxp,zkmxp,zkmx,j)
!!    call addfsech('QIA4_NO+',qia(:,:,4),zimxp,zkmxp,zkmx,j)
!     call addfsech('QIA5_N+ ',qia(:,:,5),zimxp,zkmxp,zkmx,j)
!
      DO 4 K=1,KMAX
      DO 4 I=1,LEN1
C     ****
C     ****     S7,S8=ALFA(1)*F(1),ALFA(2)*F(2) (AURORA)
C     ****     S9=ALFA*F (CUSP)
C     ****     S10=ALFA*F (DRIZZLE)
C     ****
      S7(I,K)=T6(I)*T8(I)
      S8(I,K)=T7(I)*T9(I)
      S9(I,K)=T5(I)*ALFAC*FC
      S10(I,K)=T3(I)*ALFAD*FD
    4 CONTINUE
      NMSK=NJ+NMS
      NPS1K=NJ+NPS
      NPS2K=NJ+NPS2
      NTK=NJ+NT
      DO 2 I=1,LEN2
C     ****
C     ****     S1=SUM(QI)=QT
C     ****       =SUM(F*ALFA*FX(X)/(35.E-3*H))
C     ****
      S11(I,1)=C(54)*.5*(F(I,NMSK)+F(I,NMSK+1))/(35.E-3*C(57)*F(I,NTK))
      S1(I,1)=S7(I,1)*S3(I,1)+S8(I,1)*S4(I,1)
      S1(I,1)=S1(I,1)+S9(I,1)*S5(I,1)
      S1(I,1)=(S1(I,1)+S10(I,1)*S6(I,1))*S11(I,1)
C     ****
C     ****     S3=DENOM=0.92*PS3/RMASS(3)+1.5*PS1/RMASS(1)+0.56*PS2/
C     ****     RMASS(2)
C     ****
      S3(I,1)=0.92*(1.-F(I,NPS1K)-F(I,NPS2K))/RMASS(3)+1.5*F(I,NPS1K)/
     1RMASS(1)+0.56*F(I,NPS2K)/RMASS(2)
    2 CONTINUE
C     ****
C     ****      ADD IONIZATION CONTRIBUTIONS TO NQO2P, NQOP, NQN2P, NQNP
C     ****
C     ****      S13 = QO2P, S14 = QOP, S15 = QN2P         (K+1/2)
c
c 6/98: Include production due to low energy protons.
c qia(zimxp,zkmxp,5) is output ion production from low_proton
c   (called above) for N2+, O2+, O+, NO+, N+ respectively 
c   (currently, NO+ is not defined)
c
      NPS1K=NJ+NPS
      NPS2K=NPS1K+KMAXP1
      DO 40 I=1,LEN2
        S13(I,1)=S1(I,1)*F(I,NPS1K)/(RMASS(1)*S3(I,1))+qia(i,1,2)
        S14(I,1)=S1(I,1)*(0.5*F(I,NPS1K)/RMASS(1)+
     1    0.56*F(I,NPS2K)/RMASS(2))/S3(I,1)+qia(i,1,3)
        S15(I,1)=S1(I,1)*0.7*merge(1.-F(I,NPS1K)-F(I,NPS2K),0.,1.-
     1    F(I,NPS1K)-F(I,NPS2K)>=0.)/(RMASS(3)*S3(I,1))+qia(i,1,1)
   40 CONTINUE
C     ****      SET NQO2P, NQOP, NQN2P
      DO 35 M=1,3
      K=(3-M)*KMAXP1+1
      NQPK=NQO2P+(M-1)*KMAXP1
C     ****      LEVELS 2 THRU KMAX
      DO 36 I=LEN1+1,LEN2
      F(I,NQPK)=F(I,NQPK)+SQRT(S15(I,K)*S15(I,K-1))
   36 CONTINUE
C     ****      LEVELS 1 AND KMAXP1
      DO 37 I=1,LEN1
      F(I,NQPK) = F(I,NQPK)+1.5*S15(I,K)-0.5*S15(I,K+1)
      F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+1.5*S15(I,K+KMAXM1)-0.5*
     1            S15(I,K+KMAX-2)
   37 CONTINUE
   35 CONTINUE
C     ****      SET NQNP
      NQPK=NQNP
C     ****      LEVELS 2 THRU KMAX
      DO 38 I=LEN1+1, LEN2
      F(I,NQPK)=F(I,NQPK)+.22/.7*SQRT(S15(I,1)*S15(I,0))
   38 CONTINUE
C     ****      LEVELS 1 AND KMAXP1
      DO 39 I=1,LEN1
      F(I,NQPK) = F(I,NQPK)+.22/.7*(1.5*S15(I,1)-0.5*S15(I,2))
      F(I,NQPK+KMAX) = F(I,NQPK+KMAX)+.22/.7*(1.5*S15(I,KMAX)
     1            -0.5*S15(I,KMAXM1))
   39 CONTINUE
      RETURN
      END
C
