      PROGRAM LAZAREV
*
*	Version 2.0	D P Steele, 15 August 1990
*
*	PURPOSE
*
*	This program implements a method of Lazarev, V. I., Geomagnetism 
*	and Aeronomy, v. 7, pp. 219-223, 1967, "for computing the 
*	[energy deposition] rate at various levels in the upper 
*	atmosphere during injection of electron beams of various 
*	energies."
*
*	The user chooses either a measured HiLat/DMSP spectrum or a 
*	model Maxwellian spectrum with specified energy flux and 
*	characteristic energy, either isotropic or field-aligned 
*	electron precipitation, and a date for the model atmosphere the 
*	Maxwellian spectrum will penetrate. The program loads the 
*	measured spectrum (if needed) and model atmosphere data for the
*	specified date, and uses Lazarev's approximation to calculate 
*	profiles of energy deposition rate (erg/cm^3 s), ionization rate 
*	(ion pairs/cm^3 s), N2 Lyman-Birge-Hopfield system excitation 
*	rates (excitations/cm^3 s) and N2+ First Negative system (0,1) 
*	band excitation rates (excitations/cm^3 s) leading to emission 
*	of 4278 A photons. The program selects the altitude range and 
*	step sizes on the basis of the electron spectrum. 
*
*	This program has been almost completely rewritten since the 
*	first version. It is now modular, and accepts DMSP spectra as 
*	well as HiLat data. All of the above data are now output rather 
*	than only one type of profile. The altitude range and step size 
*	choices are shameless imitations of Doug Strickland's SETUP 
*	program to his auroral code. The output format is compatible 
*	with IDL so that plotting of the results should be routine.
*
*	Modified 17 Jan 1991:
*	ReadSpectrum now fills in any gaps in the measured spectrum so 
*	there are no zero fluxes. Zero fluxes are replaced by fluxes no 
*	more than 1% of the adjacent fluxes.
*
*	USAGE
*
*	To compile, link, and run this program, use the following DCL 
*	commands:
*
*	$ FORTRAN LAZAREV
*	$ LINK LAZAREV,USER$:[STEELE.LAZAREV]E_MU_INTEGRATION, -
*	               USER$:[STEELE.LAZAREV]READDMSP, -
*		       USER$:[STEELE.UTILITIES]RTNLIB/LIB
*	$ RUN LAZAREV
*
*	REQUIREMENTS
*
*	The program obtains the center energies of the HiLat electron 
*	spectrometer channels from the file 
*
*	USER$:[STEELE.HILAT]HILDAT.DAT.
*
*	The corresponding data for the DMSP F6 and F7 satellites are 
*	obtained from the files
*
*	USER$:[STEELE.DMSP]DMSP_F6_DATA.DAT and DMSP_F7_DATA.DAT.
*
*	The program reads HiLat electron spectrometer differential 
*	number fluxes from a HiLat particle data file with name of the
*	form 
*
*	HILATXPD_yyyydddhhmmss.DAT
*
*	where	yyyy	are the 4 digits of the Year of the pass,
*		ddd	are the 3 digits of the Day,
*		hh	are the 2 digits of the start hour of the pass, 
*		mm	are the 2 digits of the start minute, and
*		ss	are the 2 digits of the start second.
*	The file of this type containing the data from the desired HiLat 
*	pass must be in WORK$:[STEELE.HILAT] when the program is run. 
*
*	The file name form for DMSP particle data is
*
*	DMSP#_yyyydddhhmmss.DAT
*
*	where # indicates the number of the satellite (e.g., 6 => F6) 
*	and the other symbols have the same meaning as above.
*
*	The program obtains MSIS-86 model atmosphere data for the 
*	selected date from a data file with name of the form
*
*	mmmddyy.MSI
*
*	where	mmm is the first 3 letters of a month name, 
*		dd is the 2 digits of the day of the month, and 
*		yy is the last 2 digits of the year. 
*	The file is read from directory WORK$:[STEELE.MSIS]. The file
*	containing the atmospheric profile for the desired day must be
*	available when the program is run. The preferred file format is 
*	as shown below, excluding the asterisks (*) at left.
*
* LATI/LONG= 58.2/256.3  H= 110.0  F10.7 = 72.6  DAY:353  LT: 1.0  AP=  0.0
* MLAT/MLON= 66.9/311.6            F10.7M= 73.8           UT: 8.0
* --------------------------------------------------------------------------
*    H/KM        NUMBER DENSITIES/CM-3     MASS DENSITY  TEMPERATURE/K
*               O         N2         O2       G*CM-3       TN   TEXO
*    85.0   3.380E+11  8.067E+12  2.025E+12  4.965E-10     196   749
*    90.0   3.172E+11  5.705E+12  1.405E+12  3.515E-10     199   749
*    95.0   2.860E+11  4.024E+12  9.683E+11  2.483E-10     204   749
*   100.0   2.501E+11  2.833E+12  6.640E+11  1.751E-10     212   749
*   105.0   2.132E+11  1.994E+12  4.537E+11  1.234E-10     224   749
*   110.0   1.781E+11  1.404E+12  3.097E+11  8.704E-11     240   749
* --------------------------------------------------------------------------
*    H/KM             NUMBER DENSITIES/CM-3            TEMPERATURE/K
*               HE         AR         H          N         TN
*    85.0   9.943E+07  7.544E+10  2.424E+07  1.206E+05     196
*    90.0   8.229E+07  5.044E+10  2.039E+07  1.365E+05     199
*    95.0   6.949E+07  3.331E+10  1.736E+07  1.545E+05     204
*   100.0   5.971E+07  2.177E+10  1.489E+07  1.747E+05     212
*   105.0   5.206E+07  1.412E+10  1.281E+07  1.970E+05     224
*   110.0   4.589E+07  9.127E+09  1.100E+07  2.210E+05     240
*
*	It is preferred that the file tabulate data at 5-km intervals 
*	from 90 km to 150 km, every 10 km from 160 to 200 km, and every 
*	25 km from 225 to 500 km, for a total of 30 heights of 
*	tabulation. The table MUST NOT contain more than 50 heights of 
*	tabulation or the program may abort.
*
*	The program uses the ratio between the N2 number densities at a 
*	given height from the Jacchia 1977 and MSIS 1986 model
*	atmospheres to normalize the LBH production rate derived from 
*	Strickland et al. The Jacchia 77 model atmosphere is read in
*	from the file WORK$:[STEELE.MSIS]JACCHIAEXT.77.
*
*	The program gets tabulated values of the pitch angle integral 
*	from one of the two binary data files 
*
*	USER$:[STEELE.LAZAREV]PAI_OUT_ISO.DAT
*
*	or
*
*	USER$:[STEELE.LAZAREV]PAI_OUT_FAL.DAT.
*
*	These files must be accessible to the program.

************************************************************************
*                                                                      *
*	Definitions and declarations                                   *
*                                                                      *
************************************************************************

      CHARACTER MonthName*9, Resp*1, DateTime*13, Sat*5, OutString*15
      CHARACTER PADist(2)*9/'Isotropic','B-aligned'/, Root*10
      LOGICAL   Measured, Model
      INTEGER   Year, Day, DoY, PAParm, Month, NSteps, Num
      REAL      E(40), Flux(40)
      REAL      MSHt(50), MSMW(50), MSND(50), MSAD(50)
      REAL      MSO(50), MSN2(50), MSO2(50), MSTN(50)
      REAL      S83J77Ht(24), S83J77N2(24), S83J77O2(24), S83J77O(24) 
      REAL      LogIbyLBH(24), S83J77AD(24), TotFlux, E0
      REAL      TabHt(30), TabAD(30), TabN2(30), TabO2(30), TabO(30)
      REAL      TabND(30), TabMW(30)
      REAL      EDepos(30), Ionize(30), LBHPro(30), I4278A(30)
      REAL      NormAtmDep(0:3200),PAIntgrl(0:3200)
      REAL*8    Start, End

      COMMON    /Depths/ NormAtmDep
      COMMON    /Int_Vals/ PAIntgrl

      EXTERNAL  Message, DayToYMD, HMS_To_Sec, GetChar, QRomb, Lagrange
      EXTERNAL  IntSpc, ReadDMSP, ISec_To_HMS

************************************************************************
*                                                                      *
*	Main program                                                   *
*                                                                      *
************************************************************************

      CALL GetInputs( Measured, Model, Sat, DateTime, Year, Month,
     &                Day, Start, End, TotFlux, E0, PAParm, Root, 
     &                OutString)
      CALL GetPADist( PAParm)
      IF (Measured) THEN
         IF (Sat .NE. 'Table') CALL GetDetParms( Sat, NSteps, E)
         CALL ReadSpectrum( Sat, DateTime, Start, End, PAParm, Root, 
     &                      NSteps, E, Flux, TotFlux, E0)
      ELSE
         CALL SetUpSpectrumGrid( E0, NSteps, E)
         CALL CalculateSpectrum( TotFlux, E0, PAParm, NSteps, E, Flux)
      END IF	! (Measured)
      CALL GetAtmosData( Year, Month, Day, Num, MSHt, MSMW, MSND, 
     &                   MSAD, MSO, MSN2, MSO2, MSTN, S83J77Ht, S83J77O,
     &                   S83J77O2, S83J77N2, LogIbyLBH, S83J77AD) 
      CALL CalculateHeights( NSteps, E, Flux, TotFlux, PAParm, Num, 
     &                       MSHt, MSAD, TabHt)
      CALL CalculateEDepos( NSteps, E, Flux, TotFlux, Num, MSHt, MSAD, 
     &                      MSND, MSMW, TabHt, EDepos, TabAD, TabND, 
     &                      TabMW) 
      CALL CalculateIonization( EDepos, Ionize)
      CALL CalculateLBHProduction( Num, TabHt, S83J77Ht, MSHt, 
     &                             LogIbyLBH, S83J77N2, S83J77O2, 
     &                             S83J77O, MSN2, MSO2, MSO, Ionize, 
     &                             LBHPro, TabN2, TabO2, TabO) 
      CALL Calculate4278Production( TabHt, Ionize, TabN2, TabO2, TabO, 
     &                              Num, I4278A)
      CALL OutputResults( Sat, OutString, PADist( PAParm), Start, End, 
     &                    NSteps, E, Flux, TotFlux, E0, TabHt, TabAD,
     &                    EDepos, Ionize, LBHPro, I4278A, TabN2, TabO2,
     &                    TabO) 
      STOP 'LAZAREV done'
      END

************************************************************************
*                                                                      *
*	SUBROUTINE GetInputs                                           *
*                                                                      *
*	Query user whether a measured or a model electron spectrum is  *
*	to be processed. If a measured spectrum, get the satellite     *
*	name and the data file descriptor. If a model Maxwellian       *
*	spectrum (with low-energy tail), get the total energy flux     *
*	and characteristic energy, as well as the date for the model   *
*	atmosphere. Get pitch-angle distribution to be assumed.        *
*                                                                      *
************************************************************************

      SUBROUTINE GetInputs( Measured, Model, Sat, DateTime, Year, Month, 
     &                      Day, Start, End, TotFlux, E0, PAParm, Root, 
     &                      OutString)

      CHARACTER Resp*1, Sat*5, MonthName*9, OutString*15, DateTime*13
      CHARACTER PADist(2)*9/'Isotropic','B-aligned'/, Root*10
      LOGICAL   Measured, Model, Tabulated
      INTEGER   DoY, Year, Day, HMS1, PAParm, C
      REAL*8    DHMS1, DT, Start, End

      TYPE 114
114   FORMAT (' This program will compute energy deposition rate, ',
     &        'ionization rate, LBH '/' production rate, and 4278A ',
     &        'production rate as a function of altitude for'/' a ',
     &        'specified electron spectrum incident on a specified ',
     &        'model atmosphere.')
      TYPE *,'Use an Actual measured spectrum or Model spectrum? [A/M]'
      CALL GetChar(Resp)
      Measured=(IN_SET(Resp,'Aa'))
      Model=.NOT.Measured
      IF (Measured) THEN
         TYPE 115
115      FORMAT(' As Measured or as Interpolated by Strickland SETUP ',
     &           'code? [M/I]')
         CALL GetChar( Resp)
         Tabulated=(IN_SET(Resp,'Ii'))
         IF (Tabulated) THEN
            Sat = 'Table'
            TYPE 116
116         FORMAT(' Enter the root part of the FMXroot.DAT file ',
     &              'containing the spectrum:')
            ACCEPT 117, Root
117         FORMAT( A)
            DECODE( 7, 118, Root(1:7)) MonthName, Day, Year
118         FORMAT( A3, I2, I2)
            DO I = 1, 3		! convert MonthName to UPPER CASE
               C = ICHAR( MonthName(I:I))
               IF (C .GT. 96) MonthName(I:I) = CHAR( C - 32)
            END DO
            Year = Year + 1900
            CALL DayToYMD( -2, Year, DoY, MonthName, Month, Day)
            ENCODE( 7, 119, DateTime(1:7)) Year, DoY
119         FORMAT( I4, I3, I6)
            IF (DoY .LT. 100) DateTime(5:5) = '0'
            DateTime(8:13) = 'FMX'//Root(8:10)
            IF (DateTime(13:13) .EQ. ' ') DateTime(13:13) = 'A'
            OutString(3:15) = DateTime
            Start = 0.0D0
            End = Start
         ELSE
            TYPE *,'Which satellite (Hilat, DMSP F6, DMSP F7)? [H/6/7]'
            CALL GetChar( Resp)
            IF (IN_SET( Resp, 'Hh')) THEN
               Sat = 'HiLat'
            ELSE IF (IN_SET( Resp, '6')) THEN
               Sat = 'DMSP6'
            ELSE IF (IN_SET( Resp, '7')) THEN
               Sat = 'DMSP7'
            ELSE
               STOP 'Unknown satellite. STOP!'
            ENDIF
            TYPE 111, Sat
111         FORMAT (' Measured ', A5, ' spectrum selected:'/
     &              ' Enter the descriptor for the data file:') 
            ACCEPT 112,DateTime
112         FORMAT (A13)
            DECODE (7,104,DateTime(1:7)) Year,DoY
104         FORMAT (I4,I3)
            CALL DayToYMD(1,Year,DoY,MonthName,Month,Day)
            OutString(3:15) = DateTime	! copy descriptor to OutString
            TYPE 113
113         FORMAT (' Enter start time, interval length for spectra:')
            ACCEPT *,DHMS1,DT
            HMS1 = INT(DHMS1)
            Frac = SNGL(DHMS1-DBLE(HMS1))
            Start = HMS_To_Sec(HMS1) + Frac
            End = Start + DT
            Root = 'Not needed'
         END IF
         TotFlux = 0.0
         E0 = 0.0
      ELSE
         TYPE 102
102      FORMAT (' Model Maxwellian spectrum selected:'/
     &           ' Enter the total energy flux (erg/cm2 s) and ',
     &              'characteristic energy (keV;<10):')
         ACCEPT *,TotFlux,E0
         ENCODE( 3, 108, OutString(10:12)) E0	! load Maxwellian
108      FORMAT( F3.1)				! characteristic energy
         OutString(11:11) = 'P'			! into identifier for
         OutString(13:15) = 'KEV'		! output file name
         TYPE 103
103      FORMAT (' Enter the date (yy mm dd) of the model atmosphere ',
     &            'to be used:')
         ACCEPT *,Year,Month,Day
         CALL DayToYMD(-1,Year,DoY,MonthName,Month,Day)
         ENCODE( 7, 106, OutString(3:9)) 1000*(1900 + Year) + DoY
106      FORMAT( I7)
         Sat = 'Model'
         DateTime = 'NoneSpecified'
         Start = 0.0D0
         End = Start
         Root = 'Not needed'
      END IF	! (Measured)
      TYPE *,'Isotropic (0) or field-aligned (1) electron influx?'
      ACCEPT *,PAParm
      PAParm=PAParm+1
      OutString(1:1) = PADist(PAParm)(1:1)	! e.g. IONIZI, IONIZB, ...
      OutString(2:2) = '_'
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE GetPADist                                           *
*                                                                      *
*	This subroutine opens the file containing the energy           *
*	deposition function (tabulated for values of the normalized    *
*	atmospheric depth from 0 to 3.2 by 0.001) corresponding to the *
*	selected pitch angle distribution (isotropic or field-         *
*	aligned), and reads the normalized atmospheric depth and       *
*	energy deposition function values into named COMMON blocks     *
*	Depths and Int_Vals.                                           *
*                                                                      *
************************************************************************

      SUBROUTINE GetPADist( PAParm)
      INTEGER   PAParm
      REAL      NormAtmDep(0:3200), PAIntgrl(0:3200)

      COMMON /Depths/ NormAtmDep
      COMMON /Int_Vals/ PAIntgrl

      IF (PAParm.EQ.1) THEN	! read isotropic integral table
         OPEN (UNIT=30,FILE='USER$:[STEELE.LAZAREV]PAI_OUT_ISO.DAT',
     &         FORM='UNFORMATTED',SHARED,STATUS='OLD',IOSTAT=IOS)
      ELSE			! read field-aligned integral table
         OPEN (UNIT=30,FILE='USER$:[STEELE.LAZAREV]PAI_OUT_FAL.DAT',
     &         FORM='UNFORMATTED',SHARED,STATUS='OLD',IOSTAT=IOS)
      END IF	! (PAParm...
      CALL MESSAGE(IOS)
      READ (30) Zero,PAIMax
      DO I=0,PAIMax	! read in appropriate table of pitch angle 
         READ (30) NormAtmDep(I),PAIntgrl(I)	! integral values
      END DO	! I=...
      CLOSE (30)
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE GetDetParms                                         *
*                                                                      *
*	This subroutine opens a file containing calibration data for   *
*	the electron detector which obtained the spectrum to be used   *
*	as input to this program. The important parameter is the       *
*	sequence of electron energies sampled by the detector.         *
*                                                                      *
************************************************************************

      SUBROUTINE GetDetParms( Sat, NSteps, E)

      CHARACTER FileName*40, Sat*5
      REAL NCAL(20), ECAL(20), E(40)

      IF (Sat.EQ.'HiLat') THEN
         FileName = 'USER$:[STEELE.HILAT]HILDAT.DAT'
         NSteps = 16
      ELSE IF (Sat.EQ.'DMSP6') THEN
         FileName = 'USER$:[STEELE.DMSP]DMSP_F6_DATA.DAT'
         NSteps = 20
      ELSE IF (Sat.EQ.'DMSP7') THEN
         FileName = 'USER$:[STEELE.DMSP]DMSP_F7_DATA.DAT'
         NSteps = 20
      ENDIF      
      OPEN (UNIT=31,FILE=FileName,SHARED,STATUS='OLD',IOSTAT=IOS)
      CALL MESSAGE(IOS)
      IF (IOS.EQ.29) THEN
         TYPE *,'Restore ', FileName, '!!', Char(7)
         STOP 
      END IF
      DO I=1,NSteps	! read center energies and other stuff
         READ (31,*,IOSTAT=IOS) E(I),NCAL(I),ECAL(I)
	 CALL MESSAGE (IOS)
      END DO	! I=...
      CLOSE (31)
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE ReadSpectrum                                        *
*                                                                      *
*	This subroutine accepts parameters identifying the measured    *
*	electron spectrum to be used as input to the energy deposition *
*	code, and returns it in array Flux. All spectra between times  *
*	Start and End are read in and averaged, and the total energy   *
*	flux and "characteristic energy" (one half the average energy) *
*	are calculated from the average dspectrum. For HiLat data in   *
*	Mode 3, only zenith spectra are used, and no account is taken  *
*	of the 40 degree spectra.                                      *
*                                                                      *
*	Modified 17 Jan 1991:                                          *
*	Zero fluxes in the measured spectrum are replaced by fluxes no *
*	more than 1% of the adjacent fluxes. This is to avoid crashing *
*	CalculateEDepos.                                               *
*                                                                      *
************************************************************************

      SUBROUTINE ReadSpectrum( Sat, DateTime, Start, End, PAParm, Root, 
     &                         NSteps, Epts, Flux, TotFlux, EChar)

      CHARACTER DateTime*13, FileName*50, Sat*5, Root*10
      LOGICAL   InWork, StartFound, EndFound, WhoKnows
      INTEGER   Year, Day, DoY, Var, FramCt, Frame, Time1(6), PAParm
      INTEGER   Pos, SatNum, Date, Time
      REAL      PrtFlx(96), TNFlux(12), TEFlux(12), AvgEn(12), OptDat(2)
      REAL      SatPos(11), TFrame(20), Spectrum(0:16,1:120), Epts(40)
      REAL      Flux(40)
      REAL*8    Start, End
*
      INTEGER   NETR, NZBR, MZBR, IMOD, IYD, MASS
      REAL      ALPER, EMAXR, Q, ZLB, ZUB, E0, EAVG, SCALE(3), TINF
      REAL      EDEN350, SCALHGT, SEC, GLAT, GLONG, STL, F107A, F107, AP
      REAL      E(100), FX(100), PITCH(20), ANGZB(40)
      REAL*8    Pi/3.1415926536/
*
      NAMELIST /PARAM/ NETR,NZBR,MZBR,ALPER,EMAXR,Q,ZLB,ZUB,E0,EAVG,
     &                 IMOD,SCALE,TINF,EDEN350,SCALHGT
      NAMELIST /INMSIS/ IYD,SEC,GLAT,GLONG,STL,F107A,F107,AP,MASS
      NAMELIST /FLAG/ WHOKNOWS
      NAMELIST /EGRID/ E,FX,PITCH,ANGZB
*
      IF (Sat.EQ.'HiLat') THEN
         FileName = 'WORK$:[STEELE.HILAT]HILATXPD_'//DateTime//'.DAT'
         INQUIRE (FILE=FileName,EXIST=InWork)
         IF (InWork) THEN
            OPEN (UNIT=32,FILE=FileName,FORM='UNFORMATTED',
     &            RECORDTYPE='VARIABLE',STATUS='OLD',IOSTAT=IOS)
            CALL MESSAGE(IOS)
         ELSE
            STOP 'HiLat file not found - check WORK$:[STEELE.HILAT]!'
         END IF	! (InWork)
         T1 = 0.0
         DO WHILE (T1.LT.Start-0.5)	! skip data until start time
            READ (32,IOSTAT=IOS) FramCt,Frame,Time1,T1,SatPos,Mode,
     &                           PrtFlx,TNFlux,TEFlux,AvgEn,OptDat
            CALL MESSAGE(IOS)
            IF (IOS.LT.0) STOP 'STOP - specified time not in this file' 
         END DO	! WHILE (T1...
         EndFound = .FALSE.
         IX=0			! initialize number of frames and 
         ISpec=0		! total number of spectra read
         IZSpec=0		! number of zenith spectra read
         NSpec=12-(Mode/2)*6	! spectra per 0.5 s in Mode Mode
         RNSpec=REAL(NSpec)	! REAL value for division
         NSteps=96/NSpec	! points per spectrum in Mode Mode
         DO WHILE((.NOT.EndFound).AND.(IOS.NE.-1).AND.(T1.LT.End+0.5))
*				! save spectral data and start times
            IX=IX+1		! increment Frame counter
            TFrame(IX)=T1	! save start time of 0.5 s data Frame
            DO Var=1,NSpec	! for each spectrum read,
               ISpec=ISpec+1	! increment total spectrum counter
               IF (Mode.LE.2) THEN	! only zenith detector operating;
                  IZSpec=IZSpec+1	! increment spectrum counter
                  Spectrum(0,IZSpec)=((REAL(Var)-1.)/RNSpec)*0.5+T1
*				! calculate t(sec) at end of spectrum
                  DO Pos=1,NSteps	! extract all spectra from PrtFlx array
                     Spectrum(Pos,IZSpec)=PrtFlx(Var+(Pos-1)*NSpec)
                  END DO	! Pos=...
               ELSE
                  IF (MOD(ISpec,3).EQ.1) THEN	! for zenith spectra,
                     IZSpec=IZSpec+1	! increment zenith counter,
                     Spectrum(0,IZSpec)=(((Var-1)/3)*3./RNSpec)*0.5+T1
                     DO Pos=1,NSteps	! and extract points from PrtFlx
                        Spectrum(Pos,IZSpec)=PrtFlx(Var+(Pos-1)*NSpec)
                     END DO	! Pos=...
                  END IF	! (MOD(ISpec,3)...
               END IF	! (Mode...
               IF ((.NOT.StartFound).AND.
     &             (Spectrum(0,IZSpec).GT.Start)) THEN
                  IZS=IZSpec	! index for first spectrum to process
                  StartFound=.TRUE.
               ELSE IF ((.NOT.EndFound).AND.
     &            (Spectrum(0,IZSpec).GT.End)) THEN
                  IZE=IZSpec-1	! index for last spectrum to process
                  EndFound=.TRUE.
               END IF	! ((.NOT.StartFound...
            END DO	! Var=...
            IF (.NOT.EndFound) THEN
               READ (32,IOSTAT=IOS) FramCt,Frame,Time1,T1,SatPos,Mode,
     &                              PrtFlx,TNFlux,TEFlux,AvgEn,OptDat
               CALL MESSAGE(IOS)
            END IF	! (.NOT.EndFound)
         END DO	! WHILE ((IOS...
         NF=IX			! total number of data frames
         NS=IZSpec		! total number of zenith spectra read in
         CLOSE (32)
         NSpec=IZE-IZS+1
         DO ISpec=IZS,IZE
            DO IE=1,NSteps
               Flux(IE)=Flux(IE)+Spectrum(IE,ISpec)/NSpec
            END DO	! IE=...
         END DO	! ISpec=...
         Flux(8)=(Flux(8)+Flux(9))/2.0
         Flux(9)=Flux(8)
      ELSE IF (Sat(1:4) .EQ. 'DMSP') THEN
         FileName = 'WORK$:[STEELE.DMSP]'//Sat//'_'//DateTime//'.DAT'
         INQUIRE (FILE=FileName,EXIST=InWork)
         IF (InWork) THEN
            DECODE( 1, 320, Sat(5:5)) SatNum
320         FORMAT( I1)
            DECODE( 13, 330, DateTime) Date, Time
330         FORMAT( I7, I6)
            CALL ReadDMSP( SatNum, Date, Time, Start, End, Flux)
            Flux(10)=(Flux(10)+Flux(11))/2.0
            Flux(11)=Flux(10)
         ELSE
            STOP 'DMSP file not found - check WORK$:[STEELE.DMSP]!'
         END IF	! (InWork)
      ELSE	! Sat = 'Table'
         L = INDEX( Root, ' ')
         IF (L .EQ. 0) L = 11
         FileName = 'WORK$:[STEELE.STRICK]FMX'//Root(1:L-1)//'.DAT'
         INQUIRE (FILE=FileName,EXIST=InWork)
         IF (InWork) THEN
            OPEN( UNIT=32, FILE=FileName, STATUS='OLD', IOSTAT=IOS)
         ELSE
           STOP 'FMX.DAT file not found - check WORK$:[STEELE.STRICK]'
         END IF
         READ( 32, NML=PARAM)
         READ( 32, NML=INMSIS)
         READ( 32, NML=FLAG)
         READ( 32, NML=EGRID)
         CLOSE( 32)
         I=0
         DO WHILE (FX(I+1) .GT. 0.0)
            I = I + 1
            Epts(I) = E(I)/1000.	! convert to keV
            Flux(I) = FX(I)*1000./Pi	! convert to el/cm2 s sr keV
         END DO	! WHILE (FX(I+1)...
         NSteps = I
      END IF	! (Sat .EQ. ...
      IFF = 1
      DO WHILE (Flux(IFF) .LE. 1.0)
         IFF = IFF + 1
      END DO
      IF (IFF .NE. 1) THEN
         DO I = 1, IFF - 1
            Flux(I) = Flux(IFF)*(Epts(IFF)/Epts(I))**2	! E^-2 dependence
         END DO
      END IF
C	The following lines fill in zero fluxes within the spectrum.
      DO IZF = IFF, NSteps
         IF (Flux(IZF) .EQ. 0.) THEN
            Flux(IZF) = 0.01*Flux(IZF-1)
            IF ((IZF .LT. NSteps) .AND. (Flux(IZF+1) .LT. Flux(IZF-1))) 
     &         Flux(IZF) = 0.01*Flux(IZF+1)
         END IF
      END DO
      CALL IntSpc( NSteps, Epts, Flux, 1, 3, FxN, TotFlux, EChar)
      EChar = EChar/2.0	! equivalent to characteristic energy
      IF (PAParm.EQ.2) TotFlux = TotFlux/4.0
      TYPE 340, TotFlux, EChar
340   FORMAT( ' Total energy flux = ', F6.2, ' erg/cm^2 s'/
     &        ' "Characteristic energy" = ', F4.2, ' keV')
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE SetUpSpectrumGrid                                   *
*                                                                      *
*	This subroutine creates an array of energy steps for use in    *
*	modelling the effects of an incident Maxwellian spectrum (with *
*	low-energy tail). The lowest energy is 20 eV and logarithmic   *
*	spacing in energy is used, with a 25% increase in energy from  *
*	one step to the next. The maximum energy integrated is 10*E0,  *
*	where E0 is the characteristic energy of the Maxwellian.       *
*                                                                      *
************************************************************************

      SUBROUTINE SetUpSpectrumGrid( E0, NSteps, E)
      REAL      E(40)

      NSteps = INT( ALOG( 10.0*E0/0.020)/ALOG( 1.25)) + 2
      E(NSteps) = 0.020
      DO I = NSteps - 1, 1, -1
         E(I) = E(I + 1)*1.25
      END DO
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE CalculateSpectrum                                   *
*                                                                      *
*	This subroutine takes the desired total energy flux and        *
*	characteristic energy of a Maxwellian spectrum, and the number *
*	and location of the energy points where it is to be tabulated, *
*	and returns an array containing the number fluxes at the       *
*	desired energies. A low-energy tail is added to the spectrum   *
*	to provide a more realistic energy deposition profile at high  *
*	energies (as described by Meier et al., JGR, 94, 13???, 1989). *
*	The subroutine also performs a numerical integration of the    *
*	tabulated spectrum and returns the calculated values of the    *
*	the total E flux and characteristic energy (one half the mean  *
*	energy).                                                       *
*                                                                      *
************************************************************************

      SUBROUTINE CalculateSpectrum( TotFlux, E0, PAParm, NSteps, E, 
     &                              Flux)
      INTEGER   PAParm
      REAL      E(40), Flux(40), IntFlux, K
      REAL*8    Pi/3.1415926536/
      DATA      ErgPerkeV/1.6021892E-09/

      K = TotFlux/(2.0*Pi*E0*E0*E0*ErgPerkeV)
      IF (PAParm.EQ.2) K = K*4.0
      IF (E0.LT.0.5) THEN
         B = 0.8*E0
      ELSE
         B = 0.1*E0 + 0.35
      END IF
      PhiMax = K*E0*EXP( -1.0)
      DO I = 1, NSteps
         Flux(I) = K*E(I)*EXP(-E(I)/E0)		! Maxwellian
         Flux(I) = Flux(I) + 0.4*PhiMax*(E0/E(I))*EXP( -E(I)/B)	! lo-E tail
      END DO
      CALL IntSpc( NSteps, E, Flux, 1, 3, FxN, IntFLux, E1)
      IF (PAParm.EQ.2) IntFlux = IntFlux/4.0
      E1 = E1/2.0
      TYPE 200, IntFlux, TotFlux
200   FORMAT( ' Flux with low-energy tail = ', F6.2, ' erg/cm2 s (', 
     &         F6.2, ')')
      TYPE 210, E1, E0
210   FORMAT( ' "Characteristic energy" with tail = ', F4.2, ' keV (', 
     &         F4.2, ')')
      TotFlux = IntFlux		! redefine TotFlux to check E conserved?
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE GetAtmosData                                        *
*                                                                      *
*	This subroutine opens a file containing model atmosphere data  *
*	for the specified date, and reads in the atmospheric           *
*	composition as a function of altitude. Various quantities are  *
*	then calculated, including the mean molecular weight, scale    *
*	height, total number density, and atmospheric depth. The       *
*	subroutine can handle arbitrarily-spaced data as long as they  *
*	are monotonic in altitude.                                     *
*                                                                      *
************************************************************************

      SUBROUTINE GetAtmosData( Year, Month, Day, Num, Ht, MW, ND, 
     &                         AD, O, N2, O2, TN, S83J77Ht, S83J77O, 
     &                         S83J77O2, S83J77N2, LogIbyLBH, S83J77AD)

      CHARACTER AtmosModelFile*40, MonthName*9, String*80, Chrs(80)*1
      INTEGER   DoY, Day, Year
      REAL      S83J77Ht(24), S83J77T(24), S83J77O(24), S83J77O2(24)
      REAL      S83J77N2(24), S83J77Num(24), S83J77Den(24)
      REAL      S83J77Pres(24), S83J77AD(24), LogIbyLBH(24)
      REAL      Ht(50), O(50), N2(50), O2(50), He(50), Ar(50), H(50)
      REAL      N(50), MassDen(50), TN(50), Texo, KBoltz
      REAL      MD(50), ND(50), MW(50), SH(50), P(50), AD(50), G(50)

      DATA      G0/9.80665/,Re/6.3712E+06/,KBoltz/1.38066E-23/
      DATA      AvogNum/6.0221367E+26/	! kg-mole

      EQUIVALENCE (String,Chrs)

      CALL DayToYMD( -1, Year, DoY, MonthName, Month, Day)
      AtmosModelFile(1:3)=MonthName(1:3)	! build up AtmosModel-
      ENCODE (2,202,AtmosModelFile(4:5)) Day	! File piece by piece
202   FORMAT (I2)
      IF (AtmosModelFile(4:4).EQ.' ') AtmosModelFile(4:4)='0'
      ENCODE (2,202,AtmosModelFile(6:7)) MOD(Year,100)
      AtmosModelFile(8:11) = '.MSI'
      AtmosModelFile = 'WORK$:[STEELE.MSIS]'//AtmosModelFile

      OPEN (UNIT=33, FILE=AtmosModelFile, FORM='FORMATTED', SHARED,
     &      STATUS='OLD', IOSTAT=IOS)
      CALL MESSAGE(IOS)
      IF (IOS.EQ.29) STOP 
     & 'Atmospheric model file not found - check WORK$:[STEELE.MSIS]!'

      DO I=1,5		! read lines of characters preceding data
         READ (33,31,IOSTAT=IOS) NChrs,(Chrs(J),J=1,NChrs)
         CALL MESSAGE(IOS)
      END DO	! I=...
31    FORMAT (Q, 80A1)
      I=0
      DO WHILE (IOS.NE.64)	! until after non-numeric line found
         I=I+1
         READ (33, *, IOSTAT=IOS) Ht(I), O(I), N2(I), O2(I), MassDen(I),
     &                            TN(I), Texo
      END DO	! WHILE (IOS...
      Num = I-1
      READ (33,31,IOSTAT=IOS) NChrs,(Chrs(J),J=1,NChrs)
      CALL MESSAGE(IOS)
      DO I=1,Num
         READ (33, *, IOSTAT=IOS) HH, He(I), Ar(I), H(I), N(I), TT
         CALL MESSAGE(IOS)
         IF ((HH.NE.Ht(I)).OR.(TT.NE.TN(I))) 
     &      STOP 'Height or temperature error in MSIS data - STOP'
      END DO	! I=...
      CLOSE (33)

*	Calculate the SI values of the tabulated quantities

      DO I=1,Num
         Ht(I)=1.0E+03*Ht(I)				! m
         O(I)=1.0E+06*O(I)				! m**-3
         N2(I)=1.0E+06*N2(I)				! "
         O2(I)=1.0E+06*O2(I)				! "
         He(I)=1.0E+06*He(I)				! "
         Ar(I)=1.0E+06*Ar(I)				! "
         H(I)=1.0E+06*H(I)				! "
         N(I)=1.0E+06*N(I)				! "
         MassDen(I)=1.0E+03*MassDen(I)			! kg m**-3

*	Calculate atmospheric depth and other quantities

         ND(I)=O(I)+N2(I)+O2(I)+He(I)+Ar(I)+H(I)+N(I)	! num dens, m**-3
         MW(I)=MassDen(I)/ND(I)				! mol wt, kg
         G(I)=G0*(Re/(Re+Ht(I)))**2			! g at Ht(I), m s**-2
         SH(I)=KBoltz*TN(I)/(MW(I)*G(I))		! scale ht, m
         Geff=G0*(Re/(Re+Ht(I)+0.5*SH(I)))**2		! g a half scale ht up
         P(I)=ND(I)*KBoltz*TN(I)			! pressure, Pa
         AD(I)=P(I)/Geff				! kg m**-2

*	Convert quantities to familiar units

         Ht(I)=1.0E-03*Ht(I)				! km
         O(I)=1.0E-06*O(I)				! cm**-3
         N2(I)=1.0E-06*N2(I)				! "
         O2(I)=1.0E-06*O2(I)				! "
         ND(I)=1.0E-06*ND(I)				! "
         MW(I)=AvogNum*MW(I)				! amu
         SH(I)=1.0E-03*SH(I)				! km
         AD(I)=0.1*AD(I)				! g cm**-2
      END DO	! I=...

      OPEN (UNIT=12, FILE='WORK$:[STEELE.MSIS]JACCHIAEXT.77',
     &      FORM='FORMATTED', SHARED, STATUS='OLD', IOSTAT=IOS)
      CALL MESSAGE(IOS)
      DO I = 1, 24
         READ (12,*) S83J77Ht(I), S83J77T(I), S83J77O(I), S83J77O2(I),
     &               S83J77N2(I), S83J77Num(I), S83J77Den(I), 
     &               S83J77Pres(I), S83J77AD(I)
      END DO
      CLOSE (12)

      OPEN( UNIT=13, FILE='USER$:[STEELE.LAZAREV]LOGIBYLBH.DAT',
     &      FORM='FORMATTED', STATUS='OLD')
      READ(13,599) LOGIBYLBH
599   FORMAT( 8F7.4)
      CLOSE( 13)

      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE CalculateHeights                                    *
*                                                                      *
*	This subroutine takes the input electron spectrum and esti-    *
*	mates the lowest altitude at which significant ionization is   *
*	produced. This altitude is found in the following way. That    *
*	electron energy is found above which only 5% of the total      *
*	energy flux lies. The atmospheric depth corresponding to the   *
*	effective range of this electron energy is determined.  The    *
*	altitude is found at which this atmospheric depth is reached   *
*	for the specified model atmosphere.  The range between this    *
*	altitude and 500 km is partitioned using altitude steps which  *
*	decrease like exp {-(i/14.5)^2} for i = 0..28.  Thus there are *
*	29 steps and 30 altitudes in the range.  If the lowest         *
*	altitude be Hm, the first step is (500 - Hm)/13.28060, and     *
*	subsequent steps decrease as indicated.  The number 13.28060   *
*	is the sum for i = 0 to 28 of exp {-(i/14.5)^2}.  This method  *
*	of partitioning the height range puts the greatest density of  *
*	steps where most of the ionization takes place.  It is         *
*	borrowed from the SETUP program of Doug Strickland's AURORA    *
*	code.                                                          *
*                                                                      *
************************************************************************

      SUBROUTINE CalculateHeights( NSteps, E, Flux, TotFlux, PAParm, 
     &                             Num, Ht, AD, TabHt)
      INTEGER   PAParm
      REAL      E(40), Flux(40), Ht(50), AD(50), TabHt(30)
      REAL*8    Pi/3.1415926536/
      DATA      ErgPerkeV/1.6021892E-09/

      I = 1
      SumFlux = Flux(1)*E(1)*(E(1)-E(2))
      DO WHILE (SumFlux.LT.0.05*(TotFlux/ErgPerkeV)*PAParm**2/Pi)
         I = I + 1
         SumFlux = SumFlux + Flux(I)*E(I)*(E(I-1)-E(I+1))/2.0
      END DO
      E05 = E(I)
      AD05 = 4.6E-6*E05**1.65	! Lazarev's effective range approximation 
      I = Num
      DO WHILE (AD(I).LT.AD05)
         I = I - 1
      END DO
      H05 = Ht(MAX( I - 1, 1))	! don't fall out the bottom
      Top = 500.0
      DH = (Top - H05)/13.28060
      TabHt(1) = Top
      DO I = 0, 28
         TabHt(I + 2) = TabHt(I + 1) - DH*EXP( -(I/14.5)**2)
      END DO
      TYPE 310, H05, Top
310   FORMAT( ' Height range covered = ', F5.1, ' - ', F5.1, ' km')
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE EstimateTotal                                       *
*                                                                      *
*	This subroutine uses the upper five values in an altitude      *
*	profile of some quantity (e.g., energy deposition rate, column *
*	excitation rate) to estimate the contribution from points      *
*	above the highest calculated value. A decaying exponential is  *
*	fit to the top of the profile and integrated. If the result is *
*	positive (i.e., if the rate is decreasing with increasing      *
*	altitude) it is added to the sum over the tabulated range to   *
*	give the total, which is returned to the calling routine.      *
*                                                                      *
************************************************************************

      SUBROUTINE EstimateTotal( TabHt, Profile, Total)
      REAL      TabHt(30), Profile(30), Total

      Total=Profile(1)*(TabHt(1)-TabHt(2))
      DO IHt=2,29
         Total=Total+Profile(IHt)*(TabHt(IHt-1)-TabHt(IHt+1))/2.0
      END DO	! IHt=...
      Total=Total+Profile(30)*(TabHt(29)-TabHt(30))/2.0
      Total=Total*1.0E5	! multiply quant/cm3 s by cm to get quant/cm2 s
      SX=0.0
      SY=0.0
      SXX=0.0
      SXY=0.0
      DO JJ = 1, 5		! Approximate the profile as a 
         X=TabHt(JJ)*1.0E5	! decaying exponential with height, and 
         Y=ALOG( Profile(JJ))	! integrate the rest of the profile.
         SX=SX+X		! A straight line fit to ln( rate)
         SY=SY+Y		! over the last five calculated points
         SXX=SXX+X*X		! gives parameters needed to evaluate 
         SXY=SXY+X*Y		! contribution from rest of profile.
      END DO	! JJ=...
      PM=5.0*SXX-SX*SX
      PP=5.0*SXY-SX*SY
      Slope=PP/PM
      HIntercept=(PM*SY-PP*SX)/(5.0*PM)
      Top=EXP( HIntercept+Slope*TabHt(1)*1.0E5)/(-Slope)
      IF (Top.GT.0.0) Total=Total+Top	! approximation breaks down if 
					! rate increases with height
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE CalculateEDepos                                     *
*                                                                      *
*	This subroutine accepts the incident electron spectrum and the *
*	altitude profile of atmospheric depth and calculates the       *
*	altitude profile of energy deposition, in the approximation of *
*	Lazarev, at the 30 height points previously tabulated for this *
*	spectrum. The normalized energy deposition curve appropriate   *
*	to the specified pitch angle distribution is used.             *
*                                                                      *
************************************************************************
 
      SUBROUTINE CalculateEDepos( NSteps, E, Flux, TotFlux, Num, Ht, AD,
     &                            ND, MW, TabHt, EDepos, TabAD, TabND, 
     &                            TabMW) 
      REAL      E(40), Flux(40), Ht(50), AD(50), ND(50), MW(50)
      REAL      TabHt(30), TabAD(30), EDepos(30), LnAD(50), Lagrange
      REAL      NormAtmDep(0:3200), PAIntgrl(0:3200), TabND(30) 
      REAL      TabMW(30)
      DATA      GperAMU, Pi/1.66054E-24, 3.1415926536/

      COMMON /Depths/ NormAtmDep
      COMMON /Int_Vals/ PAIntgrl

      EHigh = SQRT( E(1)**3/E(2))
      ELow = SQRT( E(NSteps)**3/E(NSteps-1))
      DO I = 1, Num
         LnAD(I) = LOG( AD(I))
      END DO
      I1 = Num	! begin interpolation at top of atmosphere
      Last = 1
      Inc = -1
      I2 = I1 - Inc
      DO IHt=1,30
         DO WHILE (((Last-I1)/Inc.GT.0).AND.(Ht(I1).GE.TabHt(IHt)))
            I1 = I1 + Inc
         END DO	! WHILE ((I1...
         I2 = I1 - Inc
         TabAD( IHt) = EXP( Lagrange( Ht, LnAD, TabHt(IHt), I1, I2))
         TabND( IHt) = Lagrange( AD, ND, TabAD( IHt), I1, I2)
         TabMW( IHt) = Lagrange( Ht, MW, TabHt( IHt), I1, I2)
         CALL QRomb( ELow, EHigh, EMuIntegral, TabAD(IHt), E, Flux)
         EDepos(IHt)=2.0*Pi*EMuIntegral*TabND(IHt)*TabMW(IHt)*GperAMU
      END DO	! IHt=...	! EDepos(*) contains average E dep profile
      CALL EstimateTotal( TabHt, EDepos, TotEDepos)
      TYPE 410, TotFlux - TotEDepos
410   FORMAT( ' Energy flux - deposited energy = ', F6.3, ' erg/cm^2 s')
      TYPE 420, 100*(TotFlux - TotEDepos)/TotFlux
420   FORMAT( ' Relative error = ', F6.3, '%')
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE CalculateIonization                                 *
*                                                                      *
*	This subroutine performs the trivial calculation of the local  *
*	ionization rate from the local energy deposition rate,         *
*	assuming that 35 eV per ion pair is a good approximation.      *
*                                                                      *
************************************************************************

      SUBROUTINE CalculateIonization( EDepos, Ionize)
      REAL      EDepos(30), Ionize(30)
      DATA      ErgPerkeV/1.6021892E-09/

      DO IHt=1,30	! for each grid height, get the ionization rate
         Ionize(IHt)=EDepos(IHt)/(0.035*ErgPerkeV)
      END DO	! IHt=...
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE CalculateLBHProduction                              *
*                                                                      *
*	This subroutine accepts the calculated ionization rate profile *
*	and the composition profiles for both the model atmosphere     *
*	requested by the user and the Jacchia 77 model atmosphere used *
*	by Strickland et al. (1983) in calculating LBH production rate *
*	profiles, and estimates an LBH production rate profile, using  *
*	height as the independent variable to relate the Strickland    *
*	LBH/ionization ratio to that in the present model atmosphere.  *
*	Modified 18 February 1991 by D P Steele:                       *
*	The LBH production rate was formerly scaled by the ratio of N2 *
*	densities in the user's model atmosphere and in J77, but this  *
*	seemed to break down at high exospheric temperatures. It is    *
*	scaled by the ratio of the weighted N2 fractions in the two    *
*	atmospheres, with the weights given by the peak excitation     *
*	cross-sections of the three major species (kN2, kO2, kO). The  *
*	weights were scaled from graphs in D. Lummerzheim's thesis.    *
*                                                                      *
************************************************************************

      SUBROUTINE CalculateLBHProduction( Num, TabHt, S83J77Ht, Ht, 
     &                                   LogIbyLBH, S83J77N2, S83J77O2, 
     &                                   S83J77O, N2, O2, O, Ionize, 
     &                                   LBHPro, TabN2, TabO2, TabO) 
      REAL      TabHt(30), TabN2(30), TabO2(30), TabO(30), Ionize(30)
      REAL      S83J77Ht(24), S83J77N2(24), S83J77O2(24), S83J77O(24)
      REAL      LnSJN2(24), LnSJO2(24), LnSJO(24), LogIbyLBH(24)
      REAL      Lagrange, kN2/1.00/, kO2/0.79/, kO/0.45/, LBHPro(30)
      REAL      Ht(50), N2(50), O2(50),O(50), LnMSN2(50), LnMSO2(50)
      REAL      LnMSO(50)

      DO I = 1, 24
         LnSJN2(I) = LOG( S83J77N2(I))
         LnSJO2(I) = LOG( S83J77O2(I))
         LnSJO(I) = LOG( S83J77O(I))
      END DO
      DO I = 1, Num
         LnMSN2(I) = LOG( N2(I))
         LnMSO2(I) = LOG( O2(I))
         LnMSO(I) = LOG( O(I))
      END DO
      IJ1 = 23
      IJ3 = 23
      IM1 = Num - 1
      DO I = 1, 30
         DO WHILE ((TabHt(I).LE.S83J77Ht(IJ3)).AND.(IJ3.GT.1))
            IJ3 = IJ3 - 1	! Use height as indep. variable to 
         END DO	! WHILE ((TabHt.. interpolate both ionization/LBH ratio
         IJ4 = IJ3 + 1		! and ln([N2]) in MSIS and J77 models 
         ThisLog = Lagrange( S83J77Ht, LogIbyLBH, TabHt(I), IJ3, IJ4)
         SJN2 = EXP( Lagrange( S83J77Ht, LnSJN2, TabHt(I), IJ3, IJ4))
         SJO2 = EXP( Lagrange( S83J77Ht, LnSJO2, TabHt(I), IJ3, IJ4))
         SJO = EXP( Lagrange( S83J77Ht, LnSJO, TabHt(I), IJ3, IJ4))
         DO WHILE ((TabHt(I).LE.Ht(IM1)).AND.(IM1.GT.1))
            IM1 = IM1 - 1
         END DO	! WHILE ((TabHt(...
         IM2 = IM1 + 1
         TabN2(I) = EXP( Lagrange( Ht, LnMSN2, TabHt(I), IM1, IM2))
         TabO2(I) = EXP( Lagrange( Ht, LnMSO2, TabHt(I), IM1, IM2))
         TabO(I) = EXP( Lagrange( Ht, LnMSO, TabHt(I), IM1, IM2))
         fSJN2 = kN2*SJN2/(kN2*SJN2 + kO2*SJO2 + kO*SJO)
         fMSN2 = kN2*TabN2(I)/(kN2*TabN2(I) + kO2*TabO2(I) + kO*TabO(I))
         LBHPro(I) = 0.7*(fMSN2/fSJN2)*Ionize(I)/(10**ThisLog)
      END DO	! I = ...
      CALL EstimateTotal( TabHt, LBHPro, TotLBHPro)
      TYPE 510, TotLBHPro/1E9
510   FORMAT( ' Total LBH column production rate = ', F6.2, ' kR')
      RETURN
      END

************************************************************************
*                                                                      *
*	SUBROUTINE Calculate4278Production                             *
*                                                                      *
*	This subroutine uses the relative peak values of the total     *
*	ionization cross-sections of the major atmospheric species     *
*	(N2, O2, O) to estimate the volume emission rate profile in    *
*	the N2+ First Negative system (0,1) band at 4278A. The         *
*	equation used comes from AURORA by Vallance Jones.             *
*                                                                      *
************************************************************************

      SUBROUTINE Calculate4278Production( TabHt, Ionize, TabN2, TabO2, 
     &                                    TabO, Num, I4278A)
      REAL      TabHt(30), Ionize(30)
      REAL      TabN2(30), TabO2(30), TabO(30), Lagrange, I4278A(30)
      REAL      O2Wt/1.00/, N2Wt/0.94/, OWt/0.55/, Branch/0.0215574/

      DO I = 1, 30
         WtSum = N2Wt*TabN2(I)+O2Wt*TabO2(I)+OWt*TabO(I)
         F4278 = Branch*N2Wt*TabN2(I)/WtSum
         I4278A(I) = Ionize(I)*F4278
      END DO	! I = ...
      CALL EstimateTotal( TabHt, I4278A, Tot4278A)
      TYPE 610, Tot4278A/1E9
610   FORMAT( ' Total 4278A column intensity = ', F6.2, ' kR')
      RETURN       
      END

************************************************************************
*                                                                      *
*	SUBROUTINE OutputResults                                       *
*                                                                      *
*	This subroutine outputs to a formatted disk file all the       *
*	important results of the program in a format easily readable   *
*	by IDL as well as humans. The same basic output structure is   *
*	repeated several times: a label appears on a single line,      *
*	preceding several lines of numeric data giving the profile of  *
*	the parameter named in the label.                              *
*                                                                      *
************************************************************************

      SUBROUTINE OutputResults( Sat, OutString, PAString, Start, End, 
     &                          NSteps, E, Flux, TotFlux, E0, TabHt, 
     &                          TabAD, EDepos, Ionize, LBHPro, I4278A, 
     &                          TabN2, TabO2, TabO)
      CHARACTER Sat*5, OutString*15, FileName*50, PAString*9, HMS*6
      INTEGER   NSteps
      REAL      TotFlux, E0
      REAL      E(40), Flux(40)
      REAL      TabHt(30), TabAD(30), EDepos(30), Ionize(30), LBHPro(30) 
      REAL      I4278A(30), TabN2(30), TabO2(30), TabO(30)
      REAL*8    Start, End

      IF ((Sat .EQ. 'Model') .OR. (Sat .EQ. 'Table')) THEN
         FileName = 'WORK$:[STEELE.LAZAREV]'//Sat//OutString//'.LAZ'
      ELSE
         ENCODE( 6, 910, HMS) ISec_To_HMS( SNGL( Start))
910      FORMAT( I6)
         DO I = 1, 5
            IF (HMS(I:I).EQ.' ') HMS(I:I) = '0'
         END DO
         FileName = 'WORK$:[STEELE.LAZAREV]'//Sat//OutString(1:9)//
     &              HMS//'.LAZ'
      END IF
      OPEN( 34, FILE = FileName, FORM = 'FORMATTED', STATUS = 'NEW',
     &      IOSTAT = IOS)
      WRITE( 34, 1000, IOSTAT = IOS) SNGL( End - Start), NSteps, 
     &                               TotFlux, E0, PAString
1000  FORMAT( ' Spectrum averaged over ', F5.2, ' seconds'/
     &        1X, I6, ' steps in input spectrum'/
     &        1X, F6.2, ' erg/cm!u2!n s energy flux'/
     &        1X, F6.2, ' keV characteristic energy'/
     &        1X, A9, ' pitch angle distribution assumed')
      CALL Message( IOS)
      WRITE( 34, 1010) 'Energy (keV)'
1010  FORMAT( 1X, A)
      WRITE( 34, 1030) E
1030  FORMAT( 5E13.5)
      WRITE( 34, 1010) 'Flux (electrons/cm!u2!n s sr keV)'
      WRITE( 34, 1030) Flux
      WRITE( 34, 1010) 'Altitude (km)'
      WRITE( 34, 1030) TabHt
      WRITE( 34, 1010) 'Energy deposition rate (erg/cm!u3!n s)'
      WRITE( 34, 1030) EDepos
      WRITE( 34, 1010) 'Ionization rate (ion pairs/cm!u3!n s)'
      WRITE( 34, 1030) Ionize
      WRITE( 34, 1010) 'LBH excitation rate (cm!u-3!n s!u-1!n)'
      WRITE( 34, 1030) LBHPro
      WRITE( 34, 1010) '4278 emission rate (cm!u-3!n s!u-1!n)'
      WRITE( 34, 1030) I4278A
      WRITE( 34, 1010) '[N2] (cm!u-3!n)'
      WRITE( 34, 1030) TabN2
      WRITE( 34, 1010) '[O2] (cm!u-3!n)'
      WRITE( 34, 1030) TabO2
      WRITE( 34, 1010) '[O] (cm!u-3!n)'
      WRITE( 34, 1030) TabO
      WRITE( 34, 1010) 'Atmospheric depth (g/cm!u2!n)'
      WRITE( 34, 1030) TabAD
      CLOSE( 34)
      RETURN
      END
