* Tue Jun 21 14:11:02 EDT 1994 /plas1/h1/mzh/src/test mjsanl.f
C   Mon Oct 18 14:45:37 EDT 1993 /plas5/h1/vgr/src/mjs mjsanl
* Fri Jan 11 15:22:10 EST 1991 /herschel/h1/rlm/src/test mjsanl.f
c
c XFULL has been modified as well as MJSANL to use 
c full analysis on GS5 Neptune data - checkout is 
c still in progress (as well as debugging!)
c 
c 10:30 AM  4/23/90 rlm
c
c
c Begin final modifications and check out to support 
c Neptune encounter operations
c RLM  3:25 PM Saturday July 15, 1989
c
c Modifications to support thermal cycle 
c testing on Voyager 1
c
c January 5, 1989
c
c Begin making modifications to use PMODEL to handle eletron to ion 
c and ion to electron feedthrough in fits 
c 2:12 PM, Sunday, October 2, 1988 - rlm
c
c
      COMMON/PLTSIZ/ W,H
      COMMON/ANSWER/ANS(150),CURRNT(512)
      common/full/curnew(544),ifull,jnei
      common/sedr/isedr,ist(6),isp(6)
      COMMON/SPEC/SBLK(60)
      INTEGER*2 JTB(6),JDATA(512),JNE,JTLMOD,JCLK
      DIMENSION TEMP(3)
      integer i
      integer isedr
      character*4 itype
      call initio( 1, 36,i)
      call initio( 1,-13,i)
      call initio( 1, 29,i)
 3    format(A130)
 4    read (29,3,end=5)
      go to 4
 5    continue
c
c  Initialize numerical data output file  --- asz 10/25/93
c
      call initio( 1, 28,i)
 6    read (28,3,end=7)
      go to 6
 7    continue
c
      CALL NXTMOD (ANS,JTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
     *ITYPE)
c
      call mjsanl(jtb)
      end
      SUBROUTINE MJSANL(jtb)
c
c  ******************************************************************** 
c  *                                                                  *
c  *                           MAIN                                   *
c  *                                                                  *
c  *      GENERAL MAGNETOSPHERIC PLASMA DATA ANALYSIS PROGRAM:        *
c  *                                                                  *
c  *      AN INTERACTIVE SYSTEM FOR PERFORMING NON-LINEAR LEAST       *
c  *      SQUARES FITS TO DATA OBTAINED BY THE "PLASMA SUBSYSTEM"     *
c  *      (PLS) EXPERIMENT ON BOARD THE VOYAGER 1 AND 2 SPACECRAFT    *
c  *      DURING THEIR ENCOUNTERS WITH THE OUTER PLANETS JUPITER,     *
c  *      SATURN, URANUS, AND NEPTUNE.                                *
c  *                                                                  *
c  *                                                                  *
c  *      ALSO CONTAINS COMMON BLOCK "BMODEL"                         *
c  *                                                                  *
c  ********************************************************************
c
      DIMENSION ITB(6),KTB(6)
      INTEGER*2 JTB(6),JDATA(512),JNE,JTLMOD,JCLK
c
c Differentiate between ntype and nntype to decide
c whether to read in initialization parameters from
c answer array (i.e., Mission Tape, or Summary or LONG 
c spectral tape during cruise = use of PHASE3).
c
c June 13, 1990 RLM
c
      integer*4 ntype
      integer*4 nntype
c
c  Input data type common block
c
      common/tptype/jtype,ntype
c
c  JTYPE C4:  The data tape/file type string returned by NXTMOD
c             in the last variable ITYPE. Legal possibilities are:
c             'EDR ','SUM ','SPL ', and 'VGRT' corresponding to
c             EDR, summary, spectral, and mission tapes, respectively.
c  NTYPE I4:  Index number of the tape/file type: 1=EDR; 2=SUM;
c             3=SPL; 4=VGRT. Half way through mjsanl.f, it is reset
c             to 4 for SUM and Long SPL files. NTYPE=4 then indicates
c             data tapes/files with available ANS array.
c
c Add to use SETJTL 3/23/90
c
      integer*2 jtlmo,jtloff
c
c 7/19/89 NET for Neptune
c
      integer*2 njtb(6)
c
      character*4 ifrmt1
      COMMON/HEADER/PROJID,FILEID,ISCID,TAPEID(2),IGENDT,IGENTM,FIPID(2)
     1,IFGEND,IFGENT,DPTJ1(2),IFRMT1(2),DPTJ2(2),IFRMT2(2),DPTJ3(2),IFRM
     2T3(2),DPTJ4(2),IFRMT4(2),HSPARE(18)
c
c  The only element of this common block used is the S/C id
c  stored in ISCID. 0 = Voyager 2; 1 = Voyager 1. It is loaded
c  in phase2 and is needed only in encounter mode.
c
      integer iplt,i5,icur,irec,nrep
      common /cntrl/ iplt,i5,icur,irec,nrep
c
c	iplt	Plotting flag (1=Yes, 0=No)
c	i5	input file number, default 5, alt = 15
c	icur	number for type of analysis, 
c        if negative calculate currents
c	irec	record results on unit # 17
c	nrep	number of iterations befor ploting
c
      COMMON/ANSWER/ANS(150),CURRNT(512)
c
c  The ANS array is primarily filled in NXTMOD. Only the following
c  elements are filled here:
c   ANS(20):  In encounter mode: The maximum response of any cup
c             to rigid corotating plasma.
c             In cruise mode: The proton number density copied from
c             ANS(38).
c   ANS(21):  The cup number which gave the maximum response above.
c   CURRNT(512):  The real currents (fA) in each channel for each cup.
c                 The currents are calculated in KNTCUR and passed back.
c
      common/crt/icrt,ns
c
c  ICRT: Flag to correct the L mode currents for noise. 1=Yes, 0=No
c  NS:   The number of spectra to average. 0=no averaging
c
      common/sedr/isedr,ist(6),isp(6)
c
c  ISEDR: SEDR tape/file type index: 0=No SEDR available;
c         2=SSEDR (via phase2); 3=MSSEDR (via phase3)
c  IST(6): The beginning of the SSEDR file (integer array)
c  ISP(6): The end of the SSEDR file.
c
      common/sel/isel
c
c  ISEL: Index of the main mode of operation:
c        0 = Cruise mode; 1 = Encounter mode
c
      COMMON/BMODEL/IBMODL,IDSRNN,IK(10)
c
c  BMODEL is only initialized in mjsanl.
c
c  Array curnew is for data in full up display (ions and electrons)
c
      common/full/curnew(544),ifull,jnei
c
c  CURNEW(544): real currents (fA) in each channel for each cup
c  IFULL: 0=no joint ion/electron analysis; 1= Joint L,E1,E2 analysis;
c         2= Joint M,E1,E2 analysis
c  JNEI: The # of channels in a cup for the current mode
c        128 for M, 16 for L,E1, and E2 modes
c 
c  Pass down step information for PMODEL 
c 
      common/pstep/rstl,rstm,rste1,rste2,
     +rstle,rstme,rste1i,rste2i
      COMMON/MID/NUMFOR
c
c  NUMFOR: The planet id: 0=Cruise mode; 1=Jupiter; 2=Saturn
c                         3=Uranus; 4=Neptune
c
      COMMON/PLTSIZ/W,H
c
c  W: pagelength in inches
c  H: pagewidth in inches
c
c  Character strings to hold the name of the mission data
c  and short sedr file names
c
      character*64 ftdata,ftsedr
c
c  Character strings to hold the start and stop times of the
c  attached vgr and mssedr tapes
c
      character*24 sedrtime, vgrtime
c
c  Character string to hold type of run from the RUN_TYPE environment
c
      character*3 run_type
c
      dimension bcupn(4)
c
c  Set up array for averaging currents
c
      dimension curave(512)
      dimension ncount(512)
c
c  Common block for elapsed CPU time. Not used currently.
c
      COMMON/TIMExx/ITT,IAA
      COMMON/RTIM/RCUR(6)
c
c  RCUR(6): The current six element time array of the
c           most recently read spectra.
c
      character*4 TIMESR(5)
      COMMON/CDATES/SPARA(17),TIMESR
c
c  There are many time related stuff in this common block,
c  but mjsanl uses only TIMESR
c
c  TIMESR(5):  When printing side by side the 5 4-character
c              strings they give the time of the current spectra
c              in the format: '19xx xxx xxxx:xx.xxx'
c
      COMMON/SCID/VOY,IME
c
c  VOY: 453=Voyager 1; 372=Voyager 2
c  IME: The number of channels in a cup for current mode.
c       128 for M, 16 for L,E1, and E2 modes
c
      COMMON/SCALE/YB,IDEC,XB,XSCALE
c
c  The coordinate limits of the current vs. channel # for each cup.
c  YB:  The Y axis lower limit (fA)
c  IDEC:  The number of logarithmic decades on the Y axis
c  XB:  The X axis lower limit in channel number
c  XSCALE:  The length of the X axis scale in channel numbers
c
      COMMON/SPEC/SBLK(60)
c
c  SBLK(60): Trajectory and coordinate transformation data
c            filled in xcru(phase3.f) and xspn(phase2.f).
c
      COMMON/PLCONS/SPARE(726),CN(4,3),COSB(3),COS2B(3),SIN2B(3)
     *  ,CALCUR(12),CLOCK(3),CAP(2),GAIN(2),V0,ALPHB
     *  ,NQ(4),MASS(4)
c
c  Out of all of these general plasma analysis parameters,
c  mjsanl uses only V0: 453.=Voy. 1; 372.=Voy.2
c
      dimension cupn(3,4)
      character*4 itype,jtype
      character*4 idtype(4)
c
c  Offset correction currents for L long (GS5) at Neptune
c
      dimension curoff(64)
c
c  SET UP CONSTANTS USED IN XSIMUL AND PMODEL
c
      LOGICAL*4 LD
      COMMON/NLLS/ICPLT,MP,ME,RSTEP,R1,R2,IPRT,LD
c
c  ICPLT: Fitting mode id: 0=no fit; 1=PMODEL; 3=XJC; 4=ABSO;
c         5=ABSOM; 6=BiMax. ABSO; 7=Cold beam; 9=Very cold beam;
c         >10=Terminate.
c  MP: Spectral type: 0=L; 1=M; 2=E1; 3=E2
c  ME: Spectral type: 0=E2; 1=E1; 2=M; 3=L
c  RSTEP,R1,R2: PMODEL integration constants
c  IPRT: Unit number for text out put (normally 6 for screen)
c  LD: Old single precision flag. No longer used.
c
      COMMON/ACCU/EPS,PVSTP,EVSTP,IBG1,IIIBG1,IBG2,IIIBG2,IPHGB,
     * NDOCUP(5),STPSZ,EPSMOD
c
c  The only elements used/set in mjsanl are:
c   IBG1: Flag to turn on the PMODEL detail in CPINT. 6=Yes, 0=No
c   EPSMOD: A value of 0.0 turns on all feedthrough.
c
c  Set up common block with RMARK values for E1, E2, and ion mode
c
      common/lim/rmarkp(4,3)
c
c  RMARK(i,j): The currents in fA at particular digital current levels.
c              where i: 1=L; 2=M; 3=E1; 4=E2 modes
c              and   j: 1=digital current at channel 0; 2=128; 3=255.
c
c  Time of spectra variables containing mSec from 1977
c
      REAL*8 TBEG,TSPEC
      LOGICAL LFIRST, REW
c
c  Currents in the cups
c
      REAL BARR(129,5),BPLT(129,5)
c
c  Cup temperatures in Centigrades and current levels in KNTCUR
c
      DIMENSION TEMP(3),RMARK(3)
c
c  Status word masking variables
c
      integer*2 i2m(2)
      equivalence (kmask,i2m(1))
c
c  Flag for first time around
c
      DATA LFIRST/.TRUE./
c
c  Status word mask corresponding to KMASK=000F0000
c
      data i2m /15,0/
c
c  Possible input data tape/file types
c
      data idtype/'EDR ','SUM ','SPL ','VGRT'/
c
c  Astronomical unit in km
c
      DATA AU/1.495979E8/
c
c  Clear initially all the currents
c
      data BARR/645*0.0/,BPLT/645*0.0/
c
c  Cup Normal Unit Vectors
c
      DATA CUPN/-.2962,-.1710,-.9397,+.2962,-.1710,-.9397,
     *          +0.0,.3420,-.9397,.7310,.6816,-.0349/
c
c Offset currents for L long (GS5) at Neptune
c
c The offsets averaged in the solar wind is the file:
c /data3/asz/ions/noise.ave
c The first four lines are the averages, the second four are the
c standard deviations.  The four lines correspond to the four
c cups in the order A,B,C,D.  The numbers in each line are the
c L mode channels in the order from the lower to higher.
c
c The values used below are selected for offset currents greater than 
c standard deviation currents.  THESE SHOULD BE REGARDED AS PRELIMINARY!!!
c 3/23/90 rlm
c
      data curoff/0.,0.,0.,0.,0.,0.,0.,0., 0., 0.,
     +62.,249.,334.,967.,1105.,1349.,
c End A cup
     +0.,0.,0.,0.,0.,0.,0.,0.,
     +0.,0.,0.,0.,0.,0.,0.,0.,
c End B cup
     +0.,0.,0.,0.,0.,0.,0.,0., 0., 0.,
     +95.,228.,201.,217.,276.,298.,
c End C cup
     +0.,0.,0.,0.,0.,0.,0.,0.,
     +0.,0.,0.,0.,0.,0.,0.,0./
c End D cup
c
c==========================================================
c       B E G I N   E X E C U T I O N
c==========================================================
c
c  Initialize trajectory file flag to no traj data available
c
      isedr=0
c
c  Initialize SC velocity in SC coordinates to 0.0
c
      do 1201 i=1,3
         sblk(39+i)=0.0
 1201 continue
c
c  Initialize rotation matrix from SC to 
c  cylindrical coordinates to the identity matrix
c
      do 21 i=1,9
         sblk(43+i)=0.0
   21 continue
      sblk(44)=1.0
      sblk(48)=1.0
      sblk(52)=1.0
c
c  End matrix initialization
c
c  Determine the Version of Program
c
c  The RUN_TYPE environment variable is set in the RUN command
c  executed from the script file. It value is a 3 character string
c  which is 'vgr' for cruise mode (mission tapes) and 'spl' for
c  encounter mode (spectral tapes). Read this variable to
c  determine what is the appropriate mode of operation.
c
      call getenv('RUN_TYPE',run_type)
      if (run_type .eq. 'vgr') then
         isel=0
         write(6,*)
         write(6,*)'========================================'
         write(6,*)'    Cruise Mode (Mission tape)'
         write(6,*)'========================================'
         write(6,*)
      else if (run_type .eq. 'spl') then
         isel=1
         write(6,*) 
         write(6,*)'========================================' 
         write(6,*)'    Encounter Mode (Spectral tape)'
         write(6,*)'========================================' 
         write(6,*)
      else
         write(6,*)
         write(6,*)'ERROR! Unknown mode of operation.'
         write(6,*)' The environment RUN_TYPE is incorrectly set in RUN'
         write(6,*)
         stop
      endif
   89 if (isel .eq. 1) go to 201
c
c  Select trajectory option in Cruise Mode
c
      write (6,*)'Trajectory file selection:'
      write (6,*)'0 ==> Do not mount trajectory file'
      write (6,*)'1 ==> Use short cruise file (via PHASE3)'
      read (5,*,end=89,err=89) isedr
c
c  The flag variable ISEDR should have the value 0 for no trajectory
c  or 3 for the use of PHASE3.
c
      isedr = isedr * 3
      go to 202
c
c  Select trajectory option in Encounter Mode
c
c
  201 continue
      write(6,*)'Trajectory file selection:'
      write(6,*)'0 ==> Do not mount trajectory file'
      write(6,*)'1 ==> Use short encounter file (via PHASE2)'
      read(5,*,end=89,err=89) isedr
c
c  The flag variable ISEDR should have the value 0 for no trajectory 
c  or 2 for the use of PHASE2.
c
      isedr = isedr * 2
c
c  First call to initialize and start reading file
c
  202 if (isedr .eq. 2)  CALL PHASE2(JTB,IEOD)
      if (isedr .eq. 3)  CALL PHASE3(JTB,IEOD)
      if (isedr .gt. 3)  go to 89
      if (isedr .lt. 0)  go to 89
c
c  Setup graphics window
c
      call plot( 0., 1., 3)
      call plot( 0., .1, 2)
      call nupage2 (W,H)
c
c  Zero out the common block /BMODEL/
c
      ibmodl=0
      idsrnn=0
      do 7776 i=1,10
        ik(i)=0
 7776 continue
      ik(1)=1
      ik(2)=1
      ik(3)=1
      ik(4)=1
c
c  Initialize the common block /CNTRL/
c
      icur=0
      i5 = 5
      iplt = 0
      irec = 17
      call initio (1,irec,itemp)
c
c  Define the plot size in inches
c
      W=11.0
      H=8.5
c
c  Setup runtime info
c
      CALL WRTIME(0)
      CALL TIMER
      CALL HOTIME
c
c  Initialize the text output unit number.
c  WARNING! This flag flag is used inconsistently. Changing
c  its value will redirect only some of the output.
c
      IPRT=6
c
c  Decimal time initialization
c
      TSPEC=0.0
c
c  Termination flag: 0=keep current spectra, 1=read next spectra
c                    2=terminate, 3=same data
c
      iread=1
c
c  Differentiate on type of sedr data
c
      if (isedr .eq. 0) go to 5998
      if (isedr .eq. 2) go to 2
c
c  Report mission and sedr tapes connected for cruise mode
c  and the beginning and times of these files.
c
      call getenv('ft32',ftdata)
      call getenv('ft33',ftsedr)
      call getenv('SSEDR_TIME',sedrtime)
      call getenv('VGR_TIME',vgrtime)
      write(6,*)
      write (6,83) ftsedr
   83 format('Cruise SEDR = ',a64)
      write (6,84) sedrtime
   84 format('SEDR file covers the time range:    ',a24)
      write(6,*)
      write (6,85) ftdata
   85 format('Cruise DATA = ',a64)
      write (6,86) vgrtime
   86 format('Mission file covers the time range: ',a24)
      write(6,*)
      goto 111
c
c  Warn that no trajectory data is connected
c
 5998 continue
      write (6,*)'---  No Trajectory Info Available ----'
      go to 111
c
c  Report encounter tape start and stop times
c
  2   CONTINUE
      GO TO (5,10,15,20),NUMFOR
   20 WRITE(6,6)
    6 FORMAT(22x,'***  NEPTUNE DATA DISPLAY  ***')
      WRITE(6,6000)
 6000 FORMAT('       SEDR DATA FILE COVERS:    ',5x,
     +'SPECTRAL DATA FILE STARTS AT:    ')
      WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
 6001 FORMAT('V2 - ',2I4,I3,I2,' TO '2I4,I3,I2,5x,'V2 - '2I4,I3,I2)
 6003 FORMAT('V1 - ',2I4,I3,I2,' TO '2I4,I3,I2,5x,'V1 - '2I4,I3,I2)
      WRITE(6,*)
      GO TO 111 
   15 WRITE(6,7)
    7 FORMAT(22x,'***  URANUS DATA DISPLAY  ***')
      WRITE(6,6000)
      WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
      WRITE(6,*)
      GO TO 111
   10 WRITE (6,4)
    4 FORMAT (22x,'***  SATURN DATA DISPLAY  ***')
      WRITE(6,6000)
      IF (ISCID .EQ. 1) THEN
         WRITE(6,6003)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
      ELSE
         WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
      ENDIF
      WRITE(6,*)
      GO TO 111
    5 WRITE(6,8)
    8 FORMAT(22x,'***  JUPITER DATA DISPLAY  ***')
      WRITE(6,6000)
      IF (ISCID .EQ.1) THEN
         WRITE(6,6003)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
      ELSE 
         WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4) 
      ENDIF
      WRITE(6,*)
  111 continue
      if (iread .ne. 0) go to 1333
      write (6,1334) TIMESR
 1334 format('Time of last spectrum is ',5a4)

 1333 continue
c
c  Option to do joint analysis of E1, L, E2 or E2, M, E1 
c  data sets with PMODEL and all feedthrough code turned on
c
      ifull=0
c
c  Modified to work on GS5 data with L modes; 
c  only long modes called - last spectrum read must 
c  be L-long or M.
c
      write (6,3550) 
 3550 format ('Do joint ion/electron analysis?',
     +' 2=M, E1, E2; 1=L, E1, E2; 0=NO')
      read (5,*,end=111,err=111) ifull
      if (ifull .eq. 0) go to 3551
      if ((ifull .eq. 1) .or. (ifull .eq. 2)) go to 1121
      go to 111
c     
c  Set flags appropriately and skip other inquiries
c
 1121 isptyp=1
      ipltdc=0
      ijtl=0
c
c  CLEAR BARR; SET ALL ELEMENTS TO ZERO
c
      BARR(1,1) = 0.0
      CALL MVC (BARR(1,1),BARR(2,1),2576)
c
c  VOY has already been set in phase2,3 and stored in /SCID/
c  Load this value in V0 in /PLCONS/
c
      V0=VOY
c
c  Continue execution
c
      if (isel .ne. 0) go to 597
      go to 777
c
c  Select type of data to be examined
c
 3551 continue
c
c  CLEAR BARR; SET ALL ELEMENTS TO ZERO
c
      BARR(1,1) = 0.0
      CALL MVC (BARR(1,1),BARR(2,1),2576)
      V0=VOY
      write (6,3196)
 3196 format (1h ,'Select Data Mode:')
      write (6,3195)
 3195 format (1h ,'1 = Plasma, 2 = DC Return,'
     +,' 3 = MODCAL and CURCAL, 4 = All')
            read (5,*,END=111,Err=111) isptyp
      ijtl=0
c
c Option to plot summed plasma current in DC return
c format; skip analysis
c
 6711 continue   
      WRITE(IPRT,770)
  770 FORMAT(1H ,'Integrate Plasma Data? (Analysis',
     +' Skipped)')
      write (iprt, 769)
  769 format(1h ,'2=YES for DC return spectra')
      write (iprt, 768)
  768 format(1h ,'1=YES for plasma spectra')
      write (iprt, 767)
  767 format(1h ,'0=NO')
      ipltdc = 0 
      READ(5,*,END=6711,Err=6711) ipltdc
      if ((isptyp .ne. 1) .and. (isptyp .ne. 4)) go to 777
c
c Select spectral type if plasma data
c
  597 continue
      write (6,3199)
 3199 format (1h ,'Select Spectral Type:')
      write (6,3198)
 3198 format (1h ,'    0 = all, 1 = ions only,' 
     +,' 2 = electrons only')
      write (6,3197)
 3197 format (1h ,'    3 = L only, 4 = M only, 5 = E1 only,'
     +,' 6 = E2 only, 7 = LS, E1S, E2S only')
      read (5,*,END=111,Err=111) ijtl
c
c Select time of first data to be viewed
c
  777 continue
      WRITE (6,100)
  100 FORMAT(1H ,'INPUT START TIME: IYR, IDAY, IHR, IMIN, ISEC, IMSEC')
      itb(6)=-1
      READ (5,*,END=111,Err=111)ITB
      if (itb(6) .eq. -1) go to 111
c
c  Calculate the number of msec since the beginning of 1977
c  to the time specified by the array ITB
c  WARNING: This routine does NOT correct for leap years.
c
      CALL xTYME(ITB,TBEG)
c
c  IF THE DESIRED SPECTRUM IS EARLIER IN TIME THAN THE
c  CURRENT ONE, A CALL TO VGRREW IS MADE.  THIS WILL REWIND THE TAPE.
c
c  phase2 and phase3 takes care of rewinding the traj tapes
c
c  First time around TSPEC is zero
c
      if (TBEG .ge. TSPEC) go to 774
c
c  The routines VGRREW are located in the nxtvgr.f and nxtspl.f
c  modules. They have a slightly different format.
c
      IF((TBEG.LT.TSPEC) .and. (ntype .ne. 4)) CALL VGRREW
      if (ntype .eq. 4) then
         REW = .true.
         CALL VGRREW(REW)
      endif
  774 continue
c
c  Set read flag to read next spectra
c
      IREAD=1
c
c  For combination modes, do not allow averaging
c
      if ((ijtl .lt. 3) .or. (ijtl .gt. 6)) go to 11
c
c  Option for averaging over NS spectra before doing analysis
c  (Similar, but not the same, as implementation via CAVE in KNTCUR -
c  not implemented in this code).  rlm  3/7/90
c
      write (6,*) "Input number of spectra to average (0= no ave.)"
      read (5,*) ns
      nsc=0
      do 5003 k=1,512
         curave(k)=0.0
         ncount(k)=0
 5003 continue
      GO TO 11
    1 CONTINUE
  199 continue
      WRITE (6,200)
  200 FORMAT(1H ,'READ NEXT SPECTRUM? 1=YES, 0=NO, 2=TERMINATE',
     +' 3=SAME DATA'/)
      write (6,2200)
 2200 format (1h ,'"0=NO" allows reset of spectral type selection')
      iread = -1
      READ(5,*,END=199,Err=199) IREAD
      if (iread .eq. -1)  go to 199
      IF (IREAD .EQ. 1) GO TO 11
      IF (IREAD .EQ. 0) GO TO 2
      IF (IREAD .EQ. 2) GO TO 1000
      IF (IREAD .EQ. 3) GO TO 1101
      GO TO 199
   11 CONTINUE
      CALL NXTMOD (ANS,nJTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
     *ITYPE)
c
c  Check for short versus long data (GS-5)
c
      kstat=lstat
      call SETJTL( JTLMO, JTLOFF, KSTAT, nJTB, JTLMOD, JCLK)
c
c  Returned:
c   JTLMO: =JTLMOD unless short spectra then =JTLMOD+4
c   JTLOFF:  4 for short, 0 for long modes
c   KSTAT: only the low 2 bytes are kept, rest zeroed out.
c
c  Copy over time as read by NXTMOD
c
      jtb(1)=njtb(1)
      jtb(2)=njtb(2)
      jtb(3)=njtb(3)
      jtb(4)=njtb(4)
      jtb(5)=njtb(5)
      jtb(6)=njtb(6)
c
c  CHECK FOR END OF FILE OR ERROR
c
      IF (IEOD .EQ. 0) GO TO 1001
      WRITE (6,198)
  198 FORMAT ('END OF DATA OR ERROR ENCOUNTERED IN SPECTRAL FILE')
      WRITE (6,197)
  197 FORMAT ('DO YOU WISH TO CONTINUE WITH CURRENT S/C AND PLANET?')
  195 continue
      WRITE (6,196)
  196 FORMAT ('TYPE 1 FOR YES OR 0 TO TERMINATE')
      iend = -1
      READ (5,*,END=195,Err=195) IEND
      if (iend .eq. -1)  go to 195
      IF (IEND .EQ. 1) GO TO 111
      IF (IEND .EQ. 0) GO TO 1000
      GO TO 195
c
c  INITIALIZE
c
 1001 IF ( .NOT. LFIRST) GO TO 1101
      LFIRST=.FALSE.
      CALL RUNBEG
      CALL PLSBEG
      CALL CONRD
c
c  Initialize values in the common block /SCALE/
c
      YB=.1
      IDEC=6
      XB=0.
c
c  Put data tape/file type into /TPTYPE/
c
      jtype=itype
c
c  Find corresponding flag number and put it into /TPTYPE/
c   NTYPE=1 - EDR; 2 - SUM; 3 - SPL; 4 - VGRT
c
      do 505 ii=1,4
         if (itype .eq. idtype(ii)) nntype=ii
  505 continue
      ntype=nntype
c
c  End initialization
c
 1101 continue
c
c  Reset JTLMO, JTLOFF and KSTAT
c
      call GETJTL( JTLMO, JTLOFF, KSTAT)
c
c  Check type for full up analysis and reject unselected modes
c
      if ((ifull .eq. 1) .and. (jtlmo .eq. 2)) go to 11
      if ((ifull .eq. 2) .and. (jtlmo .eq. 1)) go to 11
c
c  DECODE STATUS TO DETERMINE INSTRUMENT MODE
c
c  Keep only one byte of the status word
c
      IPLS=IAND(LSTAT,KMASK)
c
c  Move that one byte to the lowest position
c
      IPLS=IPLS/2**16
c
c  The meaning of the lowest byte of the status word:
c  Status byte = XX Y Z
c                84 2 1
c
c   Z = 1/0 = Modulator On/Off
c   Y = 1/0 = Plasma / DC Return mode
c   XX = 00 = PLS reverse grid
c        01 = Modulator voltage monitor
c        10 = Current calibration
c        11 = PLS normal grid
c
c  In case all modes were selected, do not check the status byte
c
      if (isptyp .eq. 4) go to 900
c
c IDC - MODULATOR ON (Normal .OR. Reversed)
c
      IF ((IPLS .EQ. 13 .OR. IPLS .EQ. 1) .and. isptyp .eq. 2) GO TO 112
c
c IDC - MODULATOR OFF (Normal .OR. Reversed)
c
      IF ((IPLS .EQ. 12 .OR. IPLS .EQ. 0) .and. isptyp .eq. 2) go to 112
c
c PLASMA DATA  - MODULATOR ON (Normal .OR. Reversed)
c
      IF ((IPLS .EQ. 15 .OR. IPLS .EQ. 3) .and. isptyp .eq. 1) GO TO 900
c
c PLASMA DATA  - MODULATOR OFF (Normal .OR. Reversed)
c
      IF ((IPLS .EQ. 14 .OR. IPLS .EQ. 2) .and. isptyp .eq. 1) GO TO 900
c
c Modulator Voltage Monitor (MVM / MODCAL)
c (Voltage Monitor On + Modulator On)
c
      IF (IPLS .EQ. 7 .and. isptyp .eq. 3) go to 112
c
c Current Calibration (CURCAL)
c  Calibrate Measurement Chain (Modulator Off) = A
c
      IF (IPLS .EQ. 10 .and. isptyp .eq. 3) go to 112
c
c  Background Correction (Modulator On) = B
c
      IF (IPLS .EQ. 11 .and. isptyp .eq. 3) go to 112
c
c Illegal Status: IPLS = 4, 5, 6, 8, 9
c
      GO TO 11
  900 CONTINUE
c
c CONTINUE ONLY for selected modes if plasma data
c
      if (jtlmo .gt. 4) go to 3193
      if (ijtl .eq. 0) go to 112
      idiff=ijtl*1.-jtlmod*1.
      if (idiff .eq. 2) go to 112
      if (ijtl .eq. 1 .and. jtlmod .lt. 3) go to 112
      if (ijtl .eq. 2 .and. ((jtlmod .eq. 3) .or. (jtlmod .eq. 4))) 
     +go to 112
      go to 11
c
c  Only accept short modes if it was specificaly requested
c
 3193 continue
      if ((ijtl .eq. 7) .and. (jtlmo .ge. 5)) go to 112
      go to 11
c
c  SEARCH FOR FIRST SPECTRUM AFTER START TIME
c
  112 continue
      jtlold=jtlmod
      mp=jtlmod-1
      me=4-jtlmod
c
c  Store current spectra time in /RTIM/
c
      DO 75 N=1,6
         RCUR(N)=JTB(N)
         KTB(N)=JTB(N)
   75 CONTINUE
c
c  Convert current time to msec since 1977
c
      CALL xTYME(KTB,TSPEC)
c
c  If current time is before the requested time, read next spectra
c
      IF (TSPEC .LT. TBEG) GO TO 11
c
c  DATA NUMBER TO ENGINEERING UNIT (FEMTOAMPERE) CONVERSION (NO FILTER)
c  (KNTCUR does no decoding for MVM; uses stored lookup table for 
c  DC returns - apparently from Anton's calibration data for -10 de C;
c  PLASMA and CURCAL currents converted based on current status word
c  for amplifier gain and capacitor and on JCLK for current integration 
c  time). Negative JNE turns off vgranl filtering in KNTCUR.
c
      CALL KNTCUR(ANS,JTB,JDATA,-JNE,LSTAT,JTLMOD,JCLK,TEMP,CURRNT,KSTAT
     *,IXCEL,CNOISE,RMARK)
c
c  INPUT: The standard calling sequence described in NXTMOD
c  OUTPUT:
c   CURRNT(512): The measured currents in each channel (fA)
c                if -1 - suspect data; -2 - data not received
c   KSTAT: Best guess status word
c   IXCEL: Quality flag: 1 = Ok. If larger positive, it is the
c          number of suspect channels. If negative, some channels
c          are saturated.
c   CNOISE: The noise level for the current instrument settings
c   RMARK(3): The current levels for digital currents 0, 128, 255.
c
c  Skip averaging if it was not asked for
c
      if (ns .eq. 0) go to 5005
c
c  Collect spectra for averaging. Accept only currents above
c  the minimum threshold.
c
      nsc=nsc+1
      do 5000 k=1,4*jne
c        if (currnt(k) .gt. rmark(1)) ncount(k)=ncount(k)+1
         if (currnt(k) .ge. rmark(1)) ncount(k)=ncount(k)+1
c        if (currnt(k) .gt. rmark(1)) curave(k)=curave(k)+currnt(k)
         if (currnt(k) .ge. rmark(1)) curave(k)=curave(k)+currnt(k)
 5000 continue
      jtlold=jtlmod
      if (nsc .lt. ns) go to 11
c
c  Once collected the right number of spectra, divide by count
c
      do 5001 k=1,4*jne
         if (ncount(k) .ne. 0) currnt(k)=curave(k)/(ncount(k)*1.)
 5001 continue
c
c  Reset variables used for averaging
c
      nsc=0
      do 5002 k=1,4*jne
         ncount(k)=0
         curave(k)=0.0
 5002 continue
 5005 continue
c
c  Transfer threshold information to common block /LIM/
c
      do 700 j=1,3
         rmarkp(jtlmod,j)=rmark(j)
  700 continue 
c
c For Cruise data and SUM or long SPL tapes, load
c initial values from answer array; key on 
c ntype = 4 (passed through common block TPTYPE)
c
c Check if Summary tape
c
      if (nntype .eq. 2) ntype=4
c
c Check if long spectral tape
c
      if ((nntype.eq.3) .and. (ans(38).lt.1000.0) .and.
     & (ans(38) .ne. 0.0)) ntype=4 
  775 continue
c
c  Get trajectory information for the requested time
c
      if (isedr .eq. 3)  CALL PHASE3(JTB,IEOD)
      if (isedr .eq. 2)  CALL PHASE2(JTB,IEOD)
c
c  Do not clear barr each time at this point
c
c  TRANSFER CURRENTS TO BARR
c  Put E1 or E2 currents into D-cup
c
      if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 1450
c
c  Place the # of channels in a cup into /FULL/
c
      jnei=jne
      DO 495 I=1,4
         DO 595 K=1,JNE
            BARR(K,I)=CURRNT(K+JNE*(I-1))
c
c  Also put ion currents into CURNEW in /FULL/
c
            curnew(K+JNE*(I-1))=CURRNT(K+JNE*(I-1))
  595    CONTINUE
  495 CONTINUE
      go to 1350
c
c  Move electron currents into D cup slot in BARR if IFULL=0
c
 1450 continue
      if ((ifull .eq. 1) .or. (ifull .eq. 2)) go to 3553
      DO 3595 K=1,JNE
         BARR(K,4)=CURRNT(K)
         BPLT(K,4)=BARR(K,4)
 3595 CONTINUE
      go to 3554
c 
c  Put E1 electrons into Barr(1-16,5) for full up analysis
c  Put E2 electrons into Barr(17-32,5)
c
 3553 continue
      do 3555 k=1,16
         j=(jtlmod-3)*16+k
         barr(j,5)=currnt(k)
 3555 continue
c
c  Move electrons into D cup slot in CURRNT 
c  for both full and not full up analysis
c
 3554 continue
      DO 3495 k=1,JNE
         CURRNT(K+JNE*3)=BARR(K,4)
         j=(jtlmod-3)*16+k
         curnew(j+jnei*4)=barr(j,5)
 3495 continue
 1350 continue
c
c  Place the # of channels per cup into /SCID/ and /SCALE/
c
      IME=JNE
      XSCALE=JNE
c
c  WRITE OUT TIME OF CURRENT SPECTRUM STORED in KTB
c
      CALL YDOYHM(KTB(1),KTB(2),KTB(3),KTB(4),KTB(5)+KTB(6)/1000.)
c
c  YDOYHM fills the /CDATES/ common block including the 20 character
c  string TIMESR containing the current spectrum time.
c
      WRITE(IPRT,38) TIMESR
   38 FORMAT('TIME OF CURRENT SPECTRUM IS ',5A4)
      IF (JNE .EQ. 16 .AND. JTLMOD  .EQ. 1) WRITE (IPRT,747)
  747 FORMAT('CURRENT SPECTRUM IS AN  L-MODE')
      IF (JNE .EQ. 16 .AND. JTLMOD  .EQ. 4) WRITE (IPRT,47)
   47 FORMAT('CURRENT SPECTRUM IS AN E2-MODE')
      IF (JNE .EQ. 16 .AND. JTLMOD  .EQ. 3) WRITE (IPRT,746)
  746 FORMAT('CURRENT SPECTRUM IS AN E1-MODE')
      IF (JNE .EQ. 128) WRITE (IPRT,48)
   48 FORMAT('CURRENT SPECTRUM IS AN M-MODE')
c
c  For Titan Proposal, if not at Saturn, skip this section
c
      if (numfor .ne. 2) go to 8889
c
c  Assume Titan radius = 2575 km and calculate the
c  distance of the S/C from the cloud tops.
c
      titalt=sblk(55)-2575
      write (6,8888) sblk(55),titalt
 8888 format ('  Distance to Titan center =',f9.2,
     +' = altitude of ',f9.2)
c
 8889 continue
c
c  Find B field projections along cup normals
c
      if (ntype .eq. 4) WRITE(6,5666)
 5666 FORMAT('PROJECTION OF Magnetic Field ALONG CUP NORMALS',
     *       ' (Gamma)')
      if (ntype .eq. 4) WRITE(6,5667)
 5667 FORMAT(1H ,5X,5X,'BA',9X,'BB',9X,'BC',9X,'BD',9x,'B mag')
      DO 5671 I=1,4
         BCUPN(I)=-(ans(5)*CUPN(1,I)+ans(6)*CUPN(2,I)+
     *                     ans(7)*CUPN(3,I))
 5671 CONTINUE 
      if (ntype .eq. 4) WRITE(6,5672)(BCUPN(I),I=1,4),ans(8)
 5672 FORMAT(3X,5(2x,F9.3)/)
c
c  End insert for B field projection
c
c  WRITE OUT TYPE OF DATA
c
      IF (IPLS .EQ. 3) WRITE (6,910)
  910 FORMAT ('MODULATOR ON, PLASMA, ** REVERSED SUPPRESSOR **')
      IF (IPLS .EQ. 15) WRITE (6,911)
  911 FORMAT ('MODULATOR ON, PLASMA, NORMAL SUPPRESSOR')
      IF (IPLS .EQ. 2) WRITE (6,3910)
 3910 FORMAT ('MODULATOR *OFF*, PLASMA, ** REVERSED SUPPRESSOR **')
      IF (IPLS .EQ. 14) WRITE (6,3911)
 3911 FORMAT ('MODULATOR *OFF*, PLASMA, NORMAL SUPPRESSOR')
      IF (IPLS .EQ. 13) WRITE (6,912)
  912 FORMAT ('MODULATOR ON, DC RET, NORMAL SUPPRESSOR')
      IF (IPLS .EQ. 1) WRITE (6,913)
  913 FORMAT ('MODULATOR ON, DC RET, ** REVERSED SUPPRESSOR **')
      IF (IPLS .EQ. 12) WRITE (6,914)
  914 FORMAT ('MODULATOR *OFF*, DC RET, NORMAL SUPPRESSOR')
      IF (IPLS .EQ. 0) WRITE (6,915)
  915 FORMAT ('MODULATOR *OFF*, DC RET, ** REVERSED SUPPRESSOR **')
      If (ipls .eq. 7) write (6,916)
  916 format ('Modulator Voltage Monitor Test')
      if (ipls .eq. 10) write (6,917)
  917 format ('Current Calibration Test - Modulator Off')
      if (ipls .eq. 11) write (6,918)
  918 format ('Current Calibration - Mod. On - Bckgnd. Cal.')
  101 CONTINUE
c
c  Voltage or current calibration
c
      if ((ipls .eq. 7) .or. (ipls .eq. 10) .or. 
     +(ipls .eq. 11)) go to 6695
c
c  Anything but plasma mode with normal or reverse grid with
c  either the modulator on or off.
c
      if ((ipls .ne. 3) .and. (ipls .ne. 15) .and. 
     +(ipls .ne. 2) .and. (ipls .ne. 14)) go to 995
c
c  For electron modes, skip the next max response section
c
      if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 507
c ----------------------------------------------------------------------
c
c  FIND CUP WITH MAXIMUM RESPONSE TO COROTATING, COLD BEAM OF PROTONS
c
      AMAXH=-1.
      IPH=0
      DO 500 K=1,4
         IF (AMAXH .GE. SBLK(8+K)) GO TO 500
         AMAXH=SBLK(8+K)
         IPH=K
  500 CONTINUE
      ANS(20)=0.
      IF (AMAXH .LE. 0. .OR. SBLK(IPH+12) .LE. 0.) GO TO 507
      SUM=0.
      DO 501 I=1,JNE
         SUM=SUM+BARR(I,IPH)
  501 CONTINUE
      RCUP=9.264E-4
      IF (IPH .EQ. 4) RCUP=1.110E-3
      ANS(20)=RCUP*SUM/(AMAXH*SBLK(IPH+12))
      ANS(21)=IPH
c
c Get density from mission tape if it is not the filll value of `preset`
c
      if ((ntype .eq. 4) .and. (ans(38) .lt. 1.e6)) ans(20)=ans(38)
c --------------------------------------------------------------------
 507  CONTINUE
c
c  ASK IF TO PLOT REAL CURRENTS
c
  670 continue
      WRITE(IPRT,70)
   70 FORMAT('PLOT DATA? 1=YES, 0=NO')
      iplt = -1
      READ(5,*,END=670,Err=670) IPLT
      if (iplt .eq. -1)  go to 670
      IF (IPLT .EQ. 1) GO TO 71
      if (iplt .eq. 0 .and. ipltdc .eq. 1) go to 6695
      IF (IPLT .EQ. 0) GO TO 60
      GO TO 670
c
c  Find thresholds and limits
c
   71 CONTINUE
      thresh=rmark(1)
      WRITE (IPRT,53) RMARK(1)
   53 FORMAT('THRESHOLD  IS ',1PE12.4,' FEMTOAMPS.')
      WRITE (IPRT,153) RMARK(3)
  153 FORMAT('SATURATION IS ',1PE12.4,' FEMTOAMPS.'/)
      WRITE (6,154)
  154 FORMAT ('MINIMUM AND MAXIMUM CURRENTS (IN FEMTOAMPS) ARE:')
      if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 3154
      DO 95 J=1,4
         IS=JNE*(J-1)+1
         CALL LIMITS (CURRNT(IS),IME,AMIN,AMAX,IMIN,IMAX)
         WRITE (IPRT,50) AMIN,AMAX,J
  50     FORMAT(2(1X,1PE12.3),' IN CUP ',I5)
  95  continue
      go to 995
 3154 continue
      CALL LIMITS (BPLT(1,4),IME,AMIN,AMAX,IMIN,IMAX)
      WRITE (IPRT,350) AMIN,AMAX
  350 FORMAT(2(1X,1PE12.3))
  995 continue
      write (iprt,*)
      WRITE (IPRT,430)
  430 FORMAT('WRITE OUT ALL CURRENTS? 1=YES, 0=NO')
      irite = -1
      READ (5,*,END=995,Err=995) IRITE
      if (irite .eq. -1)  go to 995
      IF (IRITE .le. 0) GO TO 6697
      IF (IRITE .EQ. 1) GO TO 1695
      GO TO 995
 1695 continue
      if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 4695
c
c  Write Out Plasma Currents for L and M Modes
c
 1696 continue
c
c  For M mode only give a way out
c
      IF (MP .EQ. 1) THEN
         WRITE (6,940)
         READ (5,*,END=1696,Err=1696) IAN
         IF (IAN .EQ. 0) GO TO 6697
      ENDIF
      DO 1941 J=1,4
         WRITE(IPRT,440) J
 440     FORMAT(' VOYAGER DATA  -   CURRENTS (FA)  CUP=',I3/)
         IM = 16
         IF (MP.EQ.1 .AND. J.NE.5) IM=128
         WRITE (IPRT,450) (BARR(I,J),I=1,IM)
  450    FORMAT(17(1X,1P8E10.3/))
 1940    continue
         if ( (j .ne. 2) .or. (mp .ne. 1)) go to 1941
         WRITE (6,940)
  940    FORMAT ('TYPE 1 TO CONTINUE ')
         READ (5,*,END=1940,Err=1940) IAN
 1941 continue
 1942 continue
      WRITE (6,940)
      READ (5,*,END=1942,Err=1942) IAN
c
c  Finished option for writing out currents
c  Now do plots; use PLTVGR for normal ions, PLTELE for normal 
c  electrons; and, PLTFULL for full analysis
c
 6697 continue
      if (ifull .ne. 0) go to 3556
      if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 6698
c
c  Go to subroutine PLTVGR to actually do plotting
c
      if ((ipls .ne. 3) .and. (ipls .ne. 15) .and. 
     +(ipls .ne. 2) .and. (ipls .ne. 14)) go to 6695
c
c  Plot only for plasma mode
c
      CALL PLTVGR(BARR,ANS,RMARK)
c
c  Case of DC style plots and no analysis for plasma data
c
      if (ipltdc .eq. 1) go to 6695
      go to 142
c
c  Write Out Plasma Currents for E1 and E2 Modes
c
 4695 continue
      IF (JTLMOD .EQ. 3) WRITE(IPRT,6440)
 6440 FORMAT(' VOYAGER DATA  -   CURRENTS (FA)  E1 MODE'/)
      IF (JTLMOD .EQ. 4) WRITE(IPRT,4440)
 4440 FORMAT(' VOYAGER DATA  -   CURRENTS (FA)  E2 MODE'/)
      WRITE (IPRT,6450) (BARR(I,4),I=1,16)
 6450 FORMAT(17(1X,1P8E10.3/))
c
c  Go to subroutine PLTELE to actually do plotting
c
 6698 continue
      if ((ipls .ne. 3) .and. (ipls .ne. 15) .and.
     +(ipls .ne. 2) .and. (ipls .ne. 14)) go to 6695
c
c  Plot only for plasma mode
c
      CALL PLTELE(BPLT,ANS,JTLMOD)
c
c  Case of DC style plots and no analysis for plasma data
c
      if (ipltdc .eq. 1) go to 6695
      go to 142
 6695 CONTINUE
c
c  Plot out and/or analyze calibration, DC, return, or PLASMA MOD OFF
c  in subroutine CALANL
c
      call calanl(ANS,JTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
     *ITYPE,ipls,ipltdc,RMARK,thresh)
c
c  Return to query "READ NEXT SPECTRUM?"
c
       go to 1
c
c  OPTION TO DO NONLINEAR LEAST SQUARES ANALYSIS on plasma data
c  (data may be from either normal or reversed suppressor modes)
c
c  N.B. If DC plots for plasma data are selected (ipltdc=1), then
c  the logic is routed through calanl and no analysis is 
c  performed on this spectrum/spectral set.
c
c
c  Plot call for full analysis scheme
c
 3556 continue
      call PLTFULL(barr, ans) 
  142 continue
      i5 = 5
   60 CONTINUE
      WRITE (IPRT,42)
   42 FORMAT('DO NONLINEAR LEAST SQUARES ANALYSIS WITH CURRENT SPECT',
     +'RA? INPUT NUMBER of ITERATIONS.')
      write (iprt,1142)
 1142 format ('A negative number uses unit 15 as a remote input',
     +' file')
      if (ipls .eq. 3) write (6,242)
  242 format (' WARNING - REVERSED SUPPRESSOR CONFIGURATION')
      ifit = -999
      read ( i5,*,END=142,Err=142) IFIT
      if (ifit .eq. -999) go to 142
      nrep = iabs(ifit)
      IF (IFIT .EQ. 0) GO TO 1
      IF (IFIT .gt. 0) GO TO 8999
      i5 = 15
      call initioco( 1, i5, itemp)
C
C  SELECT PRECISION
C
 8999 continue
      if (ifull .eq. 0) go to 3558
      write (6,3559)
 3559 format ('Do you really want to do joint ion/electron analysis '
     +'with the current spectra? 1=YES, 0=NO')
      read (i5,*, end=8999,err=8999) if2
      if (if2 .eq. 0) go to 111
      if (if2 .ne. 1) go to 8999
c
c Really do want to continue, so set up all defaults 
c for PMODEL integration steps
c
c Include unity response XJC option in as check out
c of code; no cross talk between ion and electron modes
c
 3761 continue
      write(6,*)'Do you want to use PMODEL with ion/electron ',
     & 'interaction (1)'
      write(6,*)'or unity response with no ion/electron interaction ',
     & '(2)'
      read(i5,*,end=3761,err=3761) icplt
      if (icplt .eq. 2) icplt=3
      if ((icplt .ne. 1) .and. (icplt .ne. 3)) go to 3761
c
c Option for intermediate output
c
 6659 write (6,6660)
 6660 format ('Turn on PMODEL detail in CPINT? 1=Yes, 0=No')
      read (i5,*,end=6659,err=6659) ibg1
      if (ibg1 .eq. 1) ibg1=6
      if ((ibg1 .ne. 0) .and. (ibg1 .ne. 6)) go to 6659
c
c Set old single precision flag
c
      ld=.false.
c
c Set nominal integration parameters for PMODEL
c
c  EPSMOD=0.0 turns on all feedthrough
c  R1=1.4 sets break from linear to logarithmic integration at 1.4 
c    thermal widths
c  R2=3.0 sets break from logarithmic to inverse integration at 3.0
c    thermal widths
c  Empirically, these values seem to do the integration rapidly and
c    accurately
c
c RSTEP of 7 does adequate job for E2 feedthrough into M
c RSTEP > 4 is probably necessary for a reasonable job
c
c     RSTEP=7.0
c     rstep=3.5
c
c Species into nominal mode
c
c The following values were used for fit shown at 3rd Brice 
c meeting in Lindau, FRG, Oct. 10-13, 1988.  Fit was to 
c two O+ and three e- Maxwellians.  Both n and w were fit for 
c each (but NOT the velocity vector).  Four iteractions gave 
c convergence; 16 cpu minutes per iteration on PLASMA using the 
c floating point accelerator board.
c
cc    rstl=5.0
cc    rstm=0.1
cc    rste1=4.0
cc    rste2=0.7
cc    rstle=2.0    
cc    rstme=0.5 
cc    rste1i=0.2 
cc    rste2i=0.08
c
c New defaults (test with V1 Saturn data 10/26/88)
c
      rstl=6.0
c
c M mode parameters neither checked nor verified
c 10/26/88 rlm
c
      rstm=0.5
      rste1=1.0
      rste2=0.5
      rstle=2.0    
      rstme=0.5 
      rste1i=0.5 
      rste2i=0.1
 3662 continue
      write (6,3661) rstl,rstm,rste1,rste2
 3661 format ('Default for RSTEPs (direct):',4(1x,f6.3),
     +' for L, M,',' E1, E2, respectively.')
 3665 write (6,3664)
 3664 format ('Do you wish to override? 0=NO')
      read (i5,*,end=3662, err=3662) ior
      if (ior .eq. 0) go to 3672
c
c Change defaults
c
 3666 continue
      write (6,3668)
 3668 format ('Input new values (direct)')
      read (i5,*,end=3666,err=3666) rstl,rstm,rste1,rste2
      write (6,3663) rstl,rstm,rste1,rste2
 3663 format ('RSTEPs = ',4(1x,f6.3))
c 
c Species into feedthrough mode
c 
 3672 continue
      write (6,3671) rstle,rstme,rste1i,rste2i
 3671 format ('Default for RSTEPs (fdthr):',4(1x,f6.3),
     +' for L, M,',' E1, E2, respectively.')
 3675 write (6,3674)
 3674 format ('Do you wish to override? 0=NO')
      read (i5,*,end=3672, err=3672) ior   
      if (ior .eq. 0) go to 3682
c 
c Change defaults 
c 
 3676 continue
      write (6,3678) 
 3678 format ('Input new values (fdthr)') 
      read (i5,*,end=3676,err=3676) rstle,rstme,rste1i,rste2i 
      write (6,3673) rstle,rstme,rste1i,rste2i 
 3673 format ('RSTEPs = ',4(1x,f6.3))
c
c End RSTEP specification
c
 3682 continue
      R1=0.7
      R2=3.0
      EPSMOD=0.0
 3045 continue
      WRITE (6,3040) R1,R2
 3040 FORMAT(' CURRENT PMODEL INTEGRATION CONSTANT SETTINGS ARE:'/,
     1'  R1 = ',F5.2,'  R2 = ',F5.2/,
     2'  DO YOU WISH TO OVERRIDE? 1=YES, 0=NO')
	iover = -1
      read ( i5,*,END=3045,Err=3045) IOVER
      if (iover .eq. -1) go to 3045
      IF (IOVER .EQ. 0) GO TO 3020
      IF (IOVER .EQ. 1) GO TO 3142
      GO TO 3045
 3142 continue
      WRITE (6,3042)
 3042 FORMAT('R1 = ')
      read ( i5,*,END=3142,Err=3142) R1
 3143 continue
      WRITE (6,3043)
 3043 FORMAT('R2 = ')
      read ( i5,*,END=3143,Err=3143) R2
      GO TO 3045
 3020 continue
      write(6,3560)
 3560 format ('Start joint ion/electron analysis sequence')
c
c Call XFULL to start analysis
c
      go to 4021
c
c Normal analysis/options package follows
c
 3558 if (jtlmod .lt. 3) go to 8998
      write (iprt,243)
  243 format ('Analysis code not implemented for E1 or E2')
      go to 1
 8998 continue
c
c No differentiation between single and double precision 
c on SUNs - this is a leftover from O.I.T./GSFC IBM
c 3/10/89 - rlm
c
      ld=.false.
c
c  SELECT SIMULATION OPTION
c
      go to 8996
 8997 continue
      i5 = 5
 8996 continue
      WRITE (IPRT,41)
   41 FORMAT (" SELECT FITTING OPTION:",
     |  " if negative calculate current " )
      WRITE (6,43)
   43 FORMAT (8X,'0 ==> DO NOT DO FIT')
      WRITE (6,44)
   44 FORMAT (8X,'1 ==> FIT USING PMODEL')
c      WRITE (6,45)
c   45 FORMAT (1H ,8X,'2 ==> FIT USING XJC')
      WRITE (6,52)
   52 FORMAT (8X,'3 ==> FIT USING XJC (UNITY RESPONCE)')
      WRITE (6,46)
   46 FORMAT (8X,'4 ==> FIT USING ABSO (L OR M MODES - LABCUR,LDCUR)')
      WRITE (6,55)
   55 FORMAT (8X,'5 ==> FIT USING ABSOM (M-MODE ONLY - ABCCUR,DCUR)')
      WRITE (6,123)
 123  FORMAT (8X,'6 ==> FIT USING ABSO (BIMAXWELLIAN VERSION)')
      WRITE (6,51)
   51 FORMAT (8X,'7 ==> FIT USING COLD BEAM (CUPINT,DCPINT)')
c     WRITE (6,49)
c  49 FORMAT (8X,'8 ==> FIT USING RECTANGULAR APPROX (XCUR)')
      WRITE (6,149)
  149 FORMAT (8X,'9 ==> FIT USING VERY COLD BEAM (CUPEXP)')
c      WRITE (6,163)
c  163 FORMAT (1H ,8X,'10 ==> FIT USING FAST VERSION OF XCUR (XTEST)')
      WRITE (6,263)
  263 FORMAT (5X,'GT 10 ==> TERMINATE EXECUTION ',$)
      icplt = -999
      read ( i5,*,END=8997,Err=8997) ICPLT
      if (icplt .eq. -999)  go to 8997
      IF (ICPLT .EQ. 0) GO TO 1
      icur = icplt
      icplt = iabs(icplt)
      IF (ICPLT .GT. 10) GO TO 1000
c
c  CONSISTENCY CHECK
c
      if ( (icplt .eq. 2) .or. (icplt .eq. 8) .or. (icplt .eq. 10))
     |   go to 8997
      IF (ICPLT .GE. 7) GO TO 4020
      IF (ICPLT .EQ. 1 .AND. JTLMOD .EQ. 2) GO TO 4000
      IF (ICPLT .EQ. 5 .AND. JTLMOD .EQ. 1) GO TO 4010
      IF (ICPLT .EQ. 1) GO TO 4030
      GO TO 4020
 4000 continue
      WRITE (6,4001)
 4001 FORMAT (1H /,' ** WARNING **  -- "PMODEL" SELECTED FOR M-MODE ANAL
     1YSIS. LARGE AMOUNTS OF CPU'/' TIME MAY BE USED. DO YOU WISH TO CON
     2TINUE? 1=YES, 0=NO'/)
      icont = -1
      read ( i5,*,END=4000,Err=4000) ICONT
      if ( icont .eq. -1) go to 4000
      IF (ICONT .EQ. 0) GO TO 60
      IF (ICONT .EQ. 1) GO TO 4030
      GO TO 4000
 4010 WRITE (6,4011)
 4011 FORMAT (1H /,' "ABSOM" CAN SIMULATE M-MODE DATA ONLY'//)
      GO TO 60
c
c  SET PMODEL INTEGRATION CONSTANTS
c
 4030 CONTINUE
      RSTEP=7.0
      R1=1.4
      R2=3.0
      EPSMOD=0.0
 4045 continue
      WRITE (6,4040) RSTEP,R1,R2,EPSMOD
 4040 FORMAT(1H ,' CURRENT PMODEL INTEGRATION CONSTANT SETTINGS ARE:'//
     1'  RSTEP = ',F5.2/
     1'  R1 = ',F5.2/'  R2 = ',F5.2/'  EPSMOD = ',F5.2//
     2'  DO YOU WISH TO OVERRIDE? 1=YES, 0=NO')
      iover = -1
      read ( i5,*,END=4045,Err=4045) IOVER
      if (iover .eq. -1) go to 4045
      IF (IOVER .EQ. 0) GO TO 4020
      IF (IOVER .EQ. 1) GO TO 4121
      GO TO 4045
 4121 continue
      WRITE (6,4041)
 4041 FORMAT('RSTEP = ')
      read ( i5,*,END=4121,Err=4121) RSTEP
 4142 continue
      WRITE (6,4042)
 4042 FORMAT('R1 = ')
      read ( i5,*,END=4142,Err=4142) R1
 4143 continue
      WRITE (6,4043)
 4043 FORMAT('R2 = ')
      read ( i5,*,END=4143,Err=4143) R2
 4144 continue
      WRITE (6,4044)
 4044 FORMAT('EPSMOD = ')
      read ( i5,*,END=4144,Err=4144) EPSMOD
      GO TO 4045
c
c Do special full up analysis
c
 4021 continue
      if (jtlmod .gt. 2) go to 61
      call xfull(ANS,JNE,JTLMOD,CURRNT,IXCEL,CNOISE,RMARK)
      go to 60
   61 continue
      write (6,62)
   62 format ('  Last spectrum read must be L or M')
C
C DO NONLINEAR LEAST SQUARES ANALYSIS
C
 4020 CONTINUE
      icrt=0
      if (jtlmod .ne. 1) go to 8112
      if (iread .eq. 3) go to 8112
      write (6,8111)
 8111 format ('  Correct L mode currents for noise? 1=Yes, 0=No')
      read (i5,*) icrt
      if (icrt .ne. 1) go to 8112
      do 8113 ikj=1,64
      currnt(ikj)=currnt(ikj)-curoff(ikj)
      if (currnt(ikj) .le. rmark(1)) currnt(ikj)=rmark(1)
 8113 continue
 8112 continue
      CALL XANAL(ANS,JNE,JTLMOD,CURRNT,IXCEL,CNOISE,RMARK)
      GO TO 60
 1000 CONTINUE
      RETURN
      END