c******************************************************************
c
c                 i_run.for 
c
c****************************************************************
c------------------------------------------------------------------
c       Main progamme to calculate ion 3-D velocity distribution 
c       functions from the Helios 1 and 2 raw data
c------------------------------------------------------------------
c
c-----------------------------------------------------------------
c      Initiate new common blocks needed for raw data transfer
c      and to construct the ion 3-D distribution function
c------------------------------------------------------------------ 
c
	Character*60 ltit1,ltit2
	real ct3(33,9,9),v1(33)
  
      Character*27 file,l1,l2,l3
      Character*10 amod,ael,aaz,aen  
      Character*10  mSec, seq ,  n0mSec, n0seq,  
     &                  BmSec, Bseq, AmSec, Aseq	
		 
      integer Year, DOY, Helios,Mode, BitRate, Format, DisMod, 
     &    shift, AzShift, I1A3,I2AB ,  HelCarRot,   EarCarRot
     
c     B   indicates  HDM2 data

      integer BInit(8,4), BQw(5,4),BMaxEl(4),
     &            BMaxAz(4), BMaxEn(4), BMass(4),
     &            BI1B(32), BI1Aint(32), BI2AB(8,32),
     &            BI1a3(6,7,32), BmSecE(4)

c     A   indicates  HDM1 data

      integer  AInit(8,4), AQw(5,4),    AMaxEl(4),
     &            AMaxAz(4),  AMaxEn(4),   AMass(4),
     &            AI1B(32),   AI1Aint(32), AI2AB(8,32),
     &            AI1a3(7,7,32), AmSecE(4)

       real I1Aint(5), I1B(5),min,min0

	dimension s(4)
cc	double precision s
	real*8 s

C     Common blocks for transfer of data to subroutines


      common/ MODES /  mSec,seq,Year, DOY, Helios,Mode, BitRate, 
     &                  Format, DisMod, shift, AzShift, I1A3,I2AB

      common/ ORBIT / SpinRate, Pitch, Aspect, HelLngAsc, 
     &               HelDisSun,  HelVrad, HelVnorm,  HelCarLng, 
     &               HelCarLat, HelCarRot , EarDisSun, EarCarLng, 
     &               EarCarLat, EarCarRot,  HSEangle 
 

      common/ PARAM / VpI1A, TpI1A, rNpI1A, AZpI1A , ELpI1A ,
     &                 Valpha, Talpha, rNalpha , ZeroRateI1A,
     &                 VpI1B, TpI1B, rNpI1B, ZeroRateI1B,
     &                 ZeroVarI1B

      common/ MAGN / Bx, By, Bz, BxSig, BySig, BzSig, 
     &                I1Aint, I1B

      common/ NDMc/ n0mSec, n0seq, n0Init(8), 
     &                n0Qw(5), n0MaxEl,n0MaxAz, n0MaxEn, n0Mass, 
     &                n0I1B(32),n0I1Aint(32), nd0I2AB(8,16),
     &                n0I1a3(5,5,9)

      common/ HDM1c / AmSec, Aseq, AmSecE, 
     &                  AInit,AQw, AMaxEl, AMaxAz,AMaxEn,
     &                  AMass, AI1B, AI1Aint, AI2AB, AI1a3

      common/ HDM2c / BmSec, Bseq, BmSecE, 
     &                  BInit,BQw, BMaxEl, BMaxAz,BMaxEn,
     &                  BMass, BI1B, BI1Aint, BI2AB, BI1a3

      common/ sise /iel,iaz,ien,lel,laz,len
C
      COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      COMMON/KANAL/VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
      COMMON/KANAL3/VOLT3(18),AZIMU3(18,2),ELEVA3(11)
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/MAXIMUM/MAXINT,RATINT,MAX3D(3),RAT3D,MAXI1B,RATI1B
      COMMON/AA/NEN,NAZ,NEL,LA1,LA2,LA3
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      COMMON/OUTPUT/ OUT(17,8)

      common/time/ IJAHR,ITAG,ISTUN,IMIN,ISEK

C
      DIMENSION EV(33),EEL(9),EAZ(9),FALF(9,9,9),PARA(3),MAXALF(3),R(9)
      DOUBLE PRECISION F(33,9,9),V(33),AZ(9),EL(9),DF
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      LOGICAL HDM,TEST,IPL,ALFA,EXTRA
C
      INTEGER Hhelios, MAX3D, MAXINT, MAXI1B     
      DATA EV/33*0./
      DATA EEL/9*0./
      DATA EAZ/9*0./

	ipl=.false.
c##	ipl=.true.

	write(6,*)'##### ipl=', ipl
C 
      DO 191 I = 1,32
      FINT(I) = 0.0
      FI1B(I) = 0.0
      DO 192 J = 1,7
      DO 193 K = 1,7
  193 FIS1(K,J,I) = 0.0
  192 CONTINUE
  191 CONTINUE
C
      TEST = .FALSE.
      HDM = .TRUE.
      EXTRA = .FALSE.

c      ALFA = .TRUE.            ! put on alfa separation
       alfa = .false.           ! no alfa separation
C
      NALFA = 0 
      IEXTRA = 0
      IALFA = 0
      IFERR = 0  
      INST = 1                                                    
      IGRENZ = 0                           
      NPRO1 = 3
      NPRO2 = 3
      NPRO3 = 3
C                                          
C                          
C-----------------------------------------------------------------------
C----------------- BLOCK 0  --------------------------------------------
C-----------------------------------------------------------------------
C----------------Read input parameters and raw data --------------------
C-----------------------------------------------------------------------
C
      READ(5,1212) Hhelios,Iyear,Idoy,Ihour,Iminute,NCYC
 1212 Format(i4,i4,i4,i3,i3,i3)
C      READ(5,*) NPRO1,NPRO2,NPRO3,ALFA
      ALFA = .TRUE.                        
      write(6,*) Hhelios,Iyear,Idoy,Ihour,Iminute,NPRO1,NPRO2,NPRO3,ALFA,NCYC
C 1220 Format(3X,4(I4),E10.4/3(I3),6H ALFA=L3,I3)
C      Read time and alfa, and integration precision given by NPRO1,2,3
C
C
C------------------INIT produces calibration data-----------------------------
C
       if ( Hhelios .eq. 1) CALL INIT
       if ( Hhelios .eq. 2) CALL INITB
       if ( Hhelios .eq. 1) CALL INIT31
       if ( Hhelios .eq. 2) CALL INB3
C    Azimuth angle shift in Helios 2 perihelion
       if ( Hhelios .eq. 2) CALL AZIMA 
C                          
C     Open data file
C
      l1 = '$tmp3:[tu]H'
c*      l1 = '$tmp3:[marsch]H'
      l2 = '_'
      l3 = '.aat'
      write(file,123) l1,Hhelios,Iyear,l2,Idoy,l3
  123 format(a15,i1,i2,a1,i3,a4)
C
      write(6,*) file
      open( UNIT = 2, FILE = file, STATUS = 'old')
C
C     Read NCYC number of spectra
C
      mcyc = NCYC
      IF (mcyc .EQ. 0) GO TO 9999
  100 CONTINUE
C
      Rminute = FLOAT(Iminute)

c***************************************************************

      call data(Hhelios,Iyear,Idoy,Ihour,Rminute)

c***************************************************************
C
C    Identification of arrays and other variables;  read orbit data
C
      READ(mSec,*) RSEC
      Write(6,6666) RSEC
 6666 Format(6H RSEC=E10.4) 
C
      if ( shift .eq. 0) ISHIF = 1
      if ( shift .eq. 1) ISHIF = 2
C
      write(6,*) 'shift=',ISHIF,shift,'AzShift=',AzShift
C
      ORA4 = HelDisSun
      ORA8 = HelVrad
      ORA9 = HelVnorm
      ORA15 = Aspect
      ORA17 = SpinRate
      write(6,1221) ORA4,ORA8,ORA9,ORA15,ORA17
 1221 Format('  orbital data',5(2XE10.4))
C
C    Transfer count rate data into general common blocks
C                

C    Define array coordinates and lenghts
C
      if (mode .eq. 0) HDM = .FALSE.
C
      Write(6,1234) HDM
 1234 Format(3X,2L3)
C
C     Identify instrument 
      if (I1A3 .eq. 1) INST = 1
      if (I1A3 .eq. 2) INST = 3
      WRITE(6,4321) INST
 4321 FORMAT(8H   INST=,I2)
C***************************************************
c      JEL = IEL
c      JAZ = IAZ
c      JEN = IEN
C      
c      LLEL = LEL
c      LLAZ = LAZ
c      LLEN = LEN
C
c      MEN = JEN
c      MAZ = JAZ
c      MEL = JEL
C
c      LC1 = LLEN
c      LC2 = LLAZ
c      LC3 = LLEL
c*************************************************
c *                             nov.29, 1995
      mEL = IEL
      mAZ = IAZ
      mEN = IEN
C      
      lc3 = LEL
      lc2 = LAZ
      lc1 = LEN
c*************************************************

C
      Write(6,8080)  MEN,MAZ,MEL,LC1,LC2,LC3
 8080 Format('  counting rate array coordinates',3(I3),2X,3(I3))  
C
      Write(6,8081)  mode, format
 8081 Format(2X,8H   mode=I3,8H format=I3)
C
          if ( mode .eq. 0) then
C      
      MAXINT = 1
      MAXI1B = 1
      RATINT = 0.
      RATI1B = 0.
C
      Do 1 I = 1,32
      FINT(I) = n0I1Aint(I)
      FI1B(I) = n0I1B(I)
C
      if (FINT(I) .lt. MAXINT) go to 51
      MAXINT = I
      RATINT = FINT(I)
   51 if (FI1B(I) .lt. MAXI1B) go to 1
      MAXI1B = I
      RATI1B = FI1B(I)
C
    1 Continue
C *************************************************               
c      MAX3D(3) = JEN
c      MAX3D(2) = JAZ
c      MAX3D(1) = JEL
C****************************************************
      MAX3D(3) = mEN                  ! nov. 29,1995
      MAX3D(2) = mAZ
      MAX3D(1) = mEL

      RAT3D = 0.
C
C$      Do 2 I = JEN,LLEN
c$      Do 3 J = JAZ,LLAZ
c$      Do 4 K = JEL,LLEL
c****************************************      
      Do 2 I = 1,9                      !?????
      Do 3 J = 1,5
      Do 4 K = 1,5
      FIS1(K,J,I) = n0I1A3(K,J,I)
C
      if (FIS1(K,J,I) .lt. RAT3D) go to 4
      MAX3D(1) = K
      MAX3D(2) = J
      MAX3D(3) = I
      RAT3D = FIS1(K,J,I)
C
    4 Continue
    3 Continue
    2 Continue
C
C      MAX3D(1) = n0MaxEl
C      MAX3D(2) = n0MaxAz
C      MAX3D(3) = n0MaxEn
C      RAT3D = FIS1(n0MaxEl,n0MaxAz,n0MaxEn)
C      IMAX = MAX3D(3)
C      MAXINT = IMAX
C      RATINT = FINT(IMAX) 
C      RATI1B = FI1B(IMAX)
C
             else
C  
      if (format .eq. 2  .or. format .eq. 3 ) then
C        
      MAXINT = 1
      MAXI1B = 1
      RATINT = 0.
      RATI1B = 0.
C     
      Do 5  I = 1,32
      FINT(I) = BI1Aint(I)
      FI1B(I) = BI1B(I)
C
      if (FINT(I) .lt. RATINT) go to 65
      MAXINT = I
      RATINT = FINT(I)
   65 if (FI1B(I) .lt. RATI1B) go to 5
      MAXI1B = I
      RATI1B = FI1B(I)
C
    5 Continue
C************************************************************
c      MAX3D(3) = JEN
c      MAX3D(2) = JAZ
c      MAX3D(1) = JEL
c************************************************************
      MAX3D(3) = mEN
      MAX3D(2) = mAZ
      MAX3D(1) = mEL

      RAT3D = 0.
C************************************************************
c      Do 6 I = JEN,LLEN
c      Do 7 J = JAZ,LLAZ
c      Do 8 K = JEL,LLEL
c**********************************************************
      Do 6 I = 1,32
      Do 7 J = 1,7
      Do 8 K = 1,6                       !????????
      FIS1(K,J,I) = BI1A3(K,J,I)
C
      if ( FIS1(K,J,I) .lt. RAT3D) go to 8
      RAT3D = FIS1(K,J,I)
      MAX3D(1) = K
      MAX3D(2) = J
      MAX3D(3) = I
C     
    8 Continue
    7 Continue
    6 Continue
C                         
C      MAX3D(1) = BMaxEl(1)
C      MAX3D(2) = BMaxAz(1)
C      MAX3D(3) = BMaxEn(1)
C      RAT3D = FIS1(BMaxEl(1),BMaxAz(1),BMaxEn(1))
C      IMAX = MAX3D(3)
C      RATINT = FINT(IMAX) 
C      RATI1B = FI1B(IMAX)
C
             else
C   
      MAXINT = 1
      MAXI1B = 1
      RATINT = 0.
      RATI1B = 0.
C      
      Do 9  I = 1,32
      FINT(I) = AI1Aint(I)
      FI1B(I) = AI1B(I)
C
      if (FINT(I) .lt. RATINT) go to 59
      MAXINT = I
      RATINT = FINT(I)
   59 if ( FI1B(I) .lt. RATI1B) go to 9
      MAXI1B = I
      RATI1B = FI1B(I)
C
    9 Continue
C*********************************************
c      MAX3D(3) = JEN
c      MAX3D(2) = JAZ
c      MAX3D(1) = JEL
c***********************************************      
        MAX3D(3) = mEN
        MAX3D(2) = mAZ
        MAX3D(1) = mEL
 
        RAT3D = 0.
C@@@                                      
                    write(4,*)' lc1=',lc1,' lc2=',lc2,' lc3=',lc3
                             write(4,*)'jen=',jen,'llen=',llen
	                     write(4,*)'jaz=',jaz,'llaz=',llaz
	                     write(4,*)'jel=',jel,'llel=',llel
	 
c      Do 10 I = JEN,LLEN
c      Do 11 J = JAZ,LLAZ
c      Do 12 K = JEL,LLEL

      Do 10 I = 1,32
      Do 11 J = 1,7
      Do 12 K = 1,7
      FIS1(K,J,I) = AI1A3(K,J,I)
C
      if ( FIS1(K,J,I) .lt. RAT3D) go to 12
      RAT3D = FIS1(K,J,I)
      MAX3D(3) = I
      MAX3D(2) = J
      MAX3D(1) = K
C     
   12 Continue
   11 Continue
   10 Continue

	             write(4,*)'#####  FIS1(k,j,i) #### before block 1'
c                             do i=1,32
c	                     write(4,*)'******i= for energy',i,'##'
c	                     do j=1,7 
c	                     write(4,*)'***** j= for az',j
c	                     write(4,*)(fis1(k,j,i),k=1,7)
c	                     end do
c	                     end do

	kg=2
	if(kg.eq.1)stop
C
C      MAX3D(1) = AMaxEl(1)
C      MAX3D(2) = AMaxAz(1)
C      MAX3D(3) = AMaxEn(1)
C      RAT3D = FIS1(AMaxEl(1),AMaxAz(1),AMaxEn(1))
C      IMAX = MAX3D(3)
C      MAXINT = IMAX
C      RATINT = FINT(IMAX) 
C      MAXI1B = IMAX
C      RATI1B = FI1B(IMAX)
C
            end if
C
      end if
C                                                             
C-----------------------------------------------------------------------
C----------------- BLOCK  1 --------------------------------------------
C-----------------------------------------------------------------------
C----------------- Definitions and constants ---------------------------
C----------------------------------------------------------------------- 
      PI = 3.14159265      
      PIH = PI/2.
      UMRF = PI/180.
      BAS = 0.001
C
      WRITE(6,117) (MAX3D(I),I=1,3), RAT3D, MAXINT, RATINT
  117 FORMAT('  max  ',3I3,2XE10.4,2XI3,2XE10.4)
C
C      LLC1 = MAX3D(3) + 3
C      MMEN = MAX3D(3) - 3
C      DO 111 I=MMEN,LLC1
C      DO 112 K=MEL,LC3
C      WRITE(6,113) ( FIS1(K,J,I),J=MAZ,LC2 )
C 113  FORMAT( 9(3XE10.4))
C 112  CONTINUE
C 111  CONTINUE
C 
C      DO 114 I=MEN,LLC1
C      WRITE(6,115) FINT(I),FI1B(I)
C  115 FORMAT( 2(3XE10.2))
C  114 CONTINUE
C                         
C-----------------------------------------------------------------------
C---------------- BLOCK  2 ---------------------------------------------
C-----------------------------------------------------------------------
C----------------------Test count rates?--------------------------------
C-----------------------------------------------------------------------
      IF ( ALFA ) IALFA = 1
      IF ( EXTRA) IEXTRA = 1
      IF ( .NOT.   TEST  ) GO TO 300
 200  BAS = 0.0
      TN = 50.
      VXT = 5.E7
      VYT = 0.
      VZT = 0.
      TXT = 1.0E4
      TYT = 4.0E4
      TZT = 4.0E4
      ITAG = 348
      ISTUN = 20
      IMIN = 33
      ISEC = 5
      B = 5.
      PHIB = 0.0
      EPSIB = 0.0
      BX = B*COS(PHIB*UMRF)*COS(EPSIB*UMRF)
      BY = B*SIN(PHIB*UMRF)*COS(EPSIB*UMRF)
      BZ = B*SIN(EPSIB*UMRF)
      MEN = 1
      MAZ = 5
      MEL = 2
      HDM = .TRUE.
      ISHIF = 1
      LC1 = 32
      LC2 = 7
      LC3 = 7
      GJAHR = IJAHR
      GTAG = ITAG
      INST = 1
      GSTUN = ISTUN
      GMIN  = IMIN
      GSEK  = ISEC
      MAX3D(1) = 4
      MAX3D(2) = 5
      MAX3D(3) = 6
      MENDIF = 0
      MAZDIF = 0
      MELDIF = 0                         
C
      CALL  ZTEST(TN,VXT,VYT,VZT,TXT,TYT,TZT,MEN,MAZ,MEL,HDM,
     1LC1,LC2,LC3,ISHIF,Hhelios)
      GO TO 700
C-----------------------------------------------------------------------
C----------------- BLOCK  3 --------------------------------------------
C-----------------------------------------------------------------------
C------- Determine grid points; correct count rates; -------------------
C------- eventually extrapolate count rates on the array edges ---------
C-----------------------------------------------------------------------
  300 CONTINUE
  301 NALFA = 0
      IF ( IALFA .EQ. 1 ) ALFA = .TRUE.
      IF ( IEXTRA .EQ. 1) EXTRA = .TRUE.
C
  304 CONTINUE
      MENDIF = 0
      MAZDIF = 0
      MELDIF = 0
C                          
C       The subroutine "Rand" evaluates the number of counts contained
C       beyond the limb of the distribution function and determines the logical
C       variable "EXTRA"; if  EXTRA = .TRUE. ->  call EXTRAP
C
      XFE = 0.
      CALL RAND(XFE,HDM,EXTRA,IQAL,INST)
      WRITE(6,354) EXTRA, IQUAL
  354 FORMAT(8H  EXTRA=L3,8H  IQUAL=,I2)
C
      IF ( .NOT. EXTRA  ) GO TO 400
  309 CALL EXTRAP(HDM,EXTRA,F,MENDIF,MAZDIF,MELDIF)              
C
      WRITE(6,3317) EXTRA
 3317 FORMAT('  Randextrapolation; extra=',L3)
C      DO 3311 I=5,7
C      DO 3312 K=5,7
C      WRITE(6,3313) ( F(K,J,I),J=7,13 )
C 3313 FORMAT( 9(1XE10.4))
C 3312 CONTINUE
C 3311 CONTINUE 
C 
C
C-----------------------------------------------------------------------
C---------------- BLOCK  4 ---------------- Time -----------------------
C-----------------------------------------------------------------------
  400 CONTINUE
      IJAHR = Iyear
      ITAG  = Idoy
C
      RSTUN = RSEC/3600000
      ISTUN = INT(RSTUN)
      ISEC   = 3600000*ISTUN
      RMIN  = (RSEC - ISEC)/60000
      IMIN  = INT(RMIN)
      LSEC = IMIN*60000 
      SEC  = (RSEC - ISEC - LSEC)/1000
      ISEK  = INT(SEC)
C
      OUT(13,8) = IJAHR
      OUT(14,8) = ITAG
      OUT(15,8) = ISTUN
      OUT(16,8) = IMIN
      OUT(17,8) = ISEK
C
      LP1 = LC1
      LP2 = LC2
      LP3 = LC3
      MPEN = MEN
      MPAZ = MAZ
      MPEL = MEL
C                                                        
      WRITE(6,450 ) IJAHR,ITAG,ISTUN,IMIN,ISEK
  450 FORMAT(6H JAHR=I5,6H  TAG=I5,7H  STUN=I5,6H  MIN=I5,6H  SEK=I5)
C
C-----------------------------------------------------------------------
C------------- BLOCK  5 ------------------------------------------------
C-------- Magnetic field data ------------------------------------------
C
  500 CONTINUE           
      Bx = Bx/100.
      By = By/100.
      Bz = Bz/100.
      BxSig = BxSig/100.
      BySig = BySig/100.
      BzSig = BzSig/100.
      write(6,*) ' magnetic field',Bx,By,Bz
C       
      B = SQRT(Bx*Bx + By*By + Bz*Bz)
      PHIB = ATAN2(By,Bx)/UMRF
      EPSIB = ATAN2(Bz,SQRT(Bx*Bx + By*By))/UMRF
      SBX = BxSig
      SBY = BySig
      SBZ = BzSig
      SB = SQRT(SBX*SBX + SBY*SBY + SBZ*SBZ)
      OUT(2,2) = B
      OUT(3,2) = EPSIB
      OUT(4,2) = PHIB
      OUT(5,2) = SBX
      OUT(6,2) = SBY
      OUT(7,2) = SBZ
      OUT(17,7)= SB
  502 CONTINUE
      WRITE(6,552) B,PHIB,EPSIB,SBX,SBY,SBZ,SB
  552 FORMAT(4H  B=E10.4,7H  PHIB=E10.4,8H  EPSIB=E10.4/5H SBX=E10.4,
     15H SBY=E10.4,5H SBZ=E10.4,4H SB=E10.4)
C     Calculate the angle between radial and magnetic field direction
      RPHI = ACOS(Bx/B)/UMRF
      OUT(1,2) = RPHI
      WRITE (6,553) RPHI
  553 FORMAT(7H  RPHI=E10.4)
C     
      IF ( INST .EQ. 3 ) GO TO 700
C
C-----------------------------------------------------------------------
C--------------- BLOCK  6 ----------------------------------------------
C-----------------------------------------------------------------------
  600 IF ( .NOT. ALFA ) GO TO 700
C
C         Separate the alpha particle counts from proton counts
C
  601 CONTINUE
      LC1 = LP1
      LC2 = LP2
      LC3 = LP3
      MEN = MPEN
      MAZ = MPAZ
      MEL = MPEL
      write(6,*) 'ALFA',ALFA
C
      KZAHL = 0
      MZAHL = 3
      CALL RALFA(FALF,ISHIF,HDM,ALFA,EXTRA,F,MENDIF,MAZDIF,MELDIF,VER,
     1MAXALF,MZAHL,KZAHL)
C
C      do 623 i=1,9
C      do 623 j=1,9
C      write(6,*) (FALF(i,j,k), k=1,5)
C      write(6,*) (FALF(i,j,k), k=6,9)
C  623 continue
C
  603 CONTINUE
C
C-----------------------------------------------------------------------
C----------------- BLOCK  7 --------------------------------------------
C-----------------------------------------------------------------------
C------ Construct the distribution function from count rates and -------
C-------calibration data -----------------------------------------------
C-----------------------------------------------------------------------
C
  700 CONTINUE
C
      CALL VERT(F,V,AZ,EL,MV,MA,ME,VXS,VYS,VZS,HDM,BAS,EXTRA,INST,MENDIF
     1,MAZDIF,MELDIF,ALFA,TEST,ISHIF,Hhelios)
  701 CONTINUE
C
      ICOUNT = 3
c**           icount = 0                 !  The lower limit counting rat

      IF ( NALFA .EQ. 1 ) ICOUNT = 2
C
      CALL GRENZ(F,V,AZ,EL,MV,MA,ME,INST,ICOUNT,IGREN)
C
      LH1 = (LC1 - 1)/2
      LH2 = (LC2 - 1)/2
      LH3 = (LC3 - 1)/2
      WRITE(6,751) LC1,LC2,LC3,LH1,LH2,LH3
  751 FORMAT(3X,6(I3))
      write(6,*) '  IGRENZ =',IGRENZ
      IF ( IGRENZ .EQ. 1) GO TO 301
C
c                                         to check the distribution             
       DO 707 I=1,LC1
c	write(4,*)'###### I= (for v(I))',I
	 v1(i)=v(i)*10.**(-5)                ! in unit km/s  
       DO 707 J=1,LC2
c	write(4,*)'###### J= (FOR az(j))',J
c       WRITE(4,754) ( F(I,J,K),K=1,LC3)
  754   FORMAT(9(2XE10.4))
	 do k=1,lc3
	    ct3(i,j,k)=f(i,j,k)*10.**20
	 end do
  707   CONTINUE
c	write(4,*)'Lc1=,',lc1,'*****  V(i) i=1-9 *****'
c         WRITE(4,756) ( V(I),I=1,9)
c	write(4,*) '***************   V(i) i= 10-19 ***'
c          WRITE(4,756) ( V(I),I=10,19)
c	write(4,*) 'lc2=', lc2, '***** az(j), j=1,9 ****'
c          WRITE(4,756) ( AZ(J),J=1,9)
c	write(4,*) 'lc3=', lc3, '***** el(k), j=1,9 ****'
c          WRITE(4,756) ( EL(K),K=1,9)
  756 FORMAT(9(2XE10.4))
C
C     Transform into LOGARITHM
c
c***********************************************
c                     plot F before intepretation ( call cinti)
 
	  ndev = 1
	ltit1='F'
	write(ltit2,903) IJAHR,ITAG,ISTUN,IMIN,ISEK


 903    format(5(i5,x))

	write(4,*)'###,l-mm',ltit1,ltit2

          write(4,*)'### l-mm lc1=',lc1,' lc2=',lc2,' lc3=',lc3

         call mmlticontour(AZ,lc2,V,LC1,EL,LC3,f,ndev,ltit1,ltit2)
	
	kg=1
	if(kg.eq.1)stop

c *********************************************************************
C
      DO 703 I=1,LC1
      DO 703 J=1,LC2
      DO 703 K=1,LC3
      DF = F(I,J,K)
      IF ( DF .LE. 0.)  DF = 1.D-70
  703 F(I,J,K) = DLOG(DF)
      MVZ = MV
      MAZZ= MA
      MEZ = ME
      write(6,*) 'call cinti'
C
      CALL CINTI(F,V,AZ,EL,MVZ,MAZZ,MEZ)
C 
C     write(6,*) 'check interpolation cinti'
C      DO 777 I=1,LC1                  
C      DO 777 J=1,LC2
C      WRITE(6,784) ( F(I,J,K),K=1,LC3)
C  784 FORMAT(9(2XE10.4))
C  777 CONTINUE
C


c**************** to plot the iso-contours of the distribution s(4) *****
c                              after intepretation
c                                   
c	ipl =.false.
c	write(4,*) 'ipl=',ipl
c	s(1)=700.
c	s(1)=s(1)*10.**5
c	s(2)=-15.                     ! alpha can not be zoro
c	s(3)=0.1                     ! epsilon can not be sero
c	s(4)=10.**(-10)
c	call cint(s,ipl,ierr)
cc                              ipl = .False. (or ipl = .true.)
cc                              if ierr = 1, s(4)=0.
cc                                 ierr = 0, s(4)= dexp(s(4))

c	write(4,*) 's(1)=',s(1),'s(2)=',s(2),'s(3)=',s(3)
c	write(4,*)
c	write(4,*) 'ierr=',ierr, 's(4)=',s(4)



	call iplot

	write(4,*)'***i_run.for***'
c********** finish ploting the iso-contours of the distribution ****



      FMAX = SNGL(F(MV,MA,ME))

      VL = 20.E5
  704 CALL MAXSU(VXS,VYS,VZS,FMAX,VL,IPL)
      EMAX = EXP(FMAX)
      WRITE(6,753) VXS,VYS,VZS,EMAX
  753 FORMAT(4(2X,E10.4))
      write(6,*) 'Integration limits: velocity,elevation,azimuth'
      WRITE(6,755) DVO,DVU,DELO,DELU,DAZO,DAZU
      write(6,*) '.......................... *********************'
  755 FORMAT(6(2X,D10.4))
  705 CONTINUE
C-----------------------------------------------------------------------
C----------------- BLOCK  8 --------------------------------------------
C-----------------------------------------------------------------------
C------- Calculate the moments-----------------------------------------
C-----------------------------------------------------------------------
C
  800 CONTINUE
      CALL MOMENT(MV,MA,ME,RHO,BV,BAZ,BEL,T1,T2,T3,T1AZ,T1EL,T2AZ,T2EL,
     1T3AZ,T3EL,NPRO1,NPRO2,NPRO3,VX,VY,VZ,QX,QY,QZ,ALFA,NALFA,
     2R,PHIQ,EPSQ,MOFE)
C      Write(6,*) 'Moment integration completed'






      WRITE(6,*) 'Moment integration',MOFE
      IF ( MOFE .EQ. 1) GO TO 304
      Write(6,*) nalfa,alfa
      IF ( NALFA .EQ. 0 .AND. RHO .LT. 1.)  GO TO 301
      SKAL = ( Bx*QX + By*QY + Bz*QZ )/B
      QPARX = Bx/B*SKAL
      QPARY = By/B*SKAL
      QPARZ = Bz/B*SKAL
      QSENX = QX - QPARX
      QSENY = QY - QPARY
      QSENZ = QZ - QPARZ
      QPAR = SKAL
      QSEN = SQRT(QSENX**2 + QSENY**2 + QSENZ**2)
      SKA1 = (R(1)*Bx + R(2)*By + R(3)*Bz)/B
      SKA2 = (R(4)*Bx + R(5)*By + R(6)*Bz)/B
      SKA3 = (R(7)*Bx + R(8)*By + R(9)*Bz)/B
      TPAR = T1*SKA1**2 + T2*SKA2**2 + T3*SKA3**2
      TSEN = ( T1+T2+T3-TPAR)/2.
      Write(6,*)  'TPAR          TSEN            QPAR           QSEN'
      Write(6,*) TPAR,TSEN,QPAR,QSEN
C
C     BERECHNUNG NORMIERTER WAERMESTROEME
      XMASS = 1.673E-24
      IF ( NALFA .EQ. 1) XMASS = XMASS*4.
      VTHERM = SQRT( (T1+T2+T3)/XMASS*1.381E-16)
      XNOR = RHO*XMASS/2.*(VTHERM**3)
      QPANOR = QPAR/XNOR
      QSENOR = QSEN/XNOR
      Q = SQRT( QX*QX + QY*QY + QZ*QZ )
      QNOR = Q/XNOR
      VTHERM = VTHERM/1.E5
C     Calculate the radial temperature
      TR = T1*R(1)**2 + T2*R(4)**2 + T3*R(7)**2
      ANISO = TPAR/TSEN
      V12 = T1/T2
      V13 = T1/T3
      V23 = T2/T3
C
      IF (NALFA .EQ. 1 ) GO TO 802
      OUT(6,4)  = QPAR
      OUT(7,4)  = QSEN
      OUT(8,4)  = TPAR
      OUT(9,4)  = TSEN
      OUT(10,4) = QPANOR
      OUT(11,4) = QSENOR
      OUT(12,4) = QNOR
      OUT(13,4) = VTHERM
      OUT(4,8)  = ANISO
      OUT(5,8)  = V12
      OUT(6,8)  = V13
      OUT(7,8)  = V23
      OUT(12,8) = VX
      OUT(13,8) = VY
      OUT(14,8) = VZ
      OUT(11,6) = VXS - VX
      OUT(12,6) = VYS - VY
      OUT(13,6) = VZS - VZ
C      OUT(10,7) = SQRT( (VX-VFIT(1))**2 + (VY-VFIT(2))**2 + (VZ-VFIT(3))
C     1**2 )/1.E5
      IF (IFERR .EQ. 1) OUT(10,7) = 0.
C     WRITE(6,1301) OUT(10,7)
 1301 FORMAT(6H DIFV=E10.4)
      IF ( .NOT. ALFA ) GO TO 803
      OUT(1,7)  = RHO
      OUT(2,7)  = T1+T2+T3
      GO TO 803
  802 OUT(6,6)  = QPAR
      OUT(7,6)  = QSEN
      OUT(8,6)  = TPAR
      OUT(9,6)  = TSEN
      OUT(14,4) = QPANOR
      OUT(15,4) = QSENOR
      OUT(16,4) = QNOR
      OUT(17,4) = VTHERM
      OUT(8,8)  = ANISO
      OUT(9,8)  = V12
      OUT(10,8) = V13
      OUT(11,8) = V23
      OUT(15,8) = VX
      OUT(16,8) = VY
      OUT(17,8) = VZ
      OUT(13,7) = VXS - VX
      OUT(14,7) = VYS - VY
      OUT(15,7) = VZS - VZ
C      OUT(12,7) = SQRT( (VX-VFIT(1))**2 + (VY-VFIT(2))**2 + (VZ-VFIT(3))
C     1**2 )/1.E5
      IF ( IFERR .EQ. 1) OUT(12,7) = 0.
C     WRITE(6,1302) OUT(12,7)
 1302 FORMAT(6H VDIF=E10.4)
      OUT(1,7)  = RHO/OUT(1,7)
      OUT(2,7)  =(T1+T2+T3)/OUT(2,7)
  803 CONTINUE
      IF ( TEST ) GO TO 807
C
C     Correct density and velocity by means of orbit parameters
      PARA(1) = BV
      PARA(2) = 0.
      PARA(3) = RHO
      WEL = BEL
      WAZ = BAZ
      AWAZ = WAZ
      AWEL = WEL
      ORA15 = 90.
      ORA17 = 60.
      CALL KORORB(PARA,WEL,WAZ,AWEL,AWAZ,Hhelios,
     & ORA4,ORA8,ORA9,ORA15,ORA17)
C
      IF ( NALFA .EQ. 1 .AND. ALFA ) GO TO 805
      OUT(1,3) = PARA(3)
      OUT(2,3) = PARA(1)
      OUT(3,3) = AWAZ
      OUT(4,3) = AWEL
      GO TO 806
  805 OUT(1,5) = PARA(3)
      OUT(2,5) = PARA(1)
      OUT(3,5) = AWAZ
      OUT(4,5) = AWEL
  806 CONTINUE
C
  807 CONTINUE
      write(6,*) ' sw parameter',para(3),awaz,awel,para(1)
C
  819 IF( TEST ) GO TO 910
C     WRITE(6,8091) OUT(15,2),OUT(14,2),OUT(16,2),OUT(17,2),OUT(3,7),OUT
C    1(4,7),OUT(1,8),OUT(2,8),OUT(3,8)
C     WRITE(6,8091) OUT(5,7),OUT(6,7),OUT(7,7),OUT(8,7)
 8091 FORMAT(1X,10(1X,E10.4))
  810 CONTINUE
C                    
C
C-----------------------------------------------------------------------
C----------------- BLOCK 9 ---------------------------------------------
C-----------------------------------------------------------------------
C--------- Evaluation of the alpha particles ---------------------------
C-----------------------------------------------------------------------
C
  910 CONTINUE
      IF ( NALFA .EQ. 1) GO TO 1103
      IF ( .NOT. ALFA  ) GO TO 1103
C
      write(6,*) 'alpha particle evaluation'
      DO 1000 I=1,33
      V(I) = 0.D0
      DO 1000 J=1,9
      AZ(J) = 0.D0
      DO 1000 K=1,9
      EL(K) = 0.D0
 1000 F(I,J,K) = 0.D0
      CALL ALVER(F,V,AZ,EL,FALF,MVALF,MAALF,MEALF,AVXS,AVYS,AVZS,BAS,
     1 ISHIF,Hhelios)
      LC1 = LA1
      LC2 = LA2
      LC3 = LA3
      MV = MVALF
      ME = MEALF
      MA = MAALF
      VXS = AVXS
      VYS = AVYS
      VZS = AVZS
      NALFA = 1
      write(6,*) 'NALFA =',NALFA
      GO TO 701
C
C-----------------------------------------------------------------------
C----------------- BLOCK 11 --------------------------------------------
C----------- write out moments and solar wind parameters ---------------
C-----------------------------------------------------------------------
C
C     IF ( NALFA .EQ. 0 ) GO TO 601
 1103 CONTINUE
C
      IF ( TEST ) GO TO 1161
C     WRITE(6,1159)  ( OUT(I,6) ,I=14,17)
 1159 FORMAT(3H   3(E10.4),1H F6.2)
C     WRITE(6,1180) OUT(5,3),OUT(6,3),OUT(7,3),OUT(9,
C    17),OUT(5,5),OUT(6,5),OUT(7,5),OUT(11,7)
 1180 FORMAT(8(2X,E10.4))
      DR = SQRT( OUT(15,6)**2 +OUT(14,6)**2 )
      IF ( .NOT. ALFA )  GO TO 1161
      DALF = ATAN2(OUT(15,6),OUT(14,6))/UMRF
      DEPS = ATAN2(OUT(16,6),DR)/UMRF
C     WRITE(6,1160) DALF,DEPS,OUT(17,6)
 1160 FORMAT(3H   2(F7.2),1H F6.2)
C     Alfvenvelocity 
C
 1161 VA = B/SQRT(1.673*4.*PI*( OUT(1,3) + OUT(1,5)*4.))*100.
      OUT(10,6 ) = VA
C
C      WRITE(6,1153) ITAG,ISTUN,IMIN,ISEK,OUT(8,1),INST,HDM,EXTRA,IFE,
C     & ALFA,VER,OUT(16,2),OUT(14,2),B,PHIB,EPSIB,RPHI,VA
 1153 FORMAT('      ZEIT      R(AU)  INST    HDM   EXTRA   XFE   ALFA
     &VER    ST1   B(GAMMA)  PHIB    EPSIB   RPHI    VA',
     &/2H  4I3,2H  F5.3,2H  I3,2(L7),3H   I5,2H  L4,3H   F5.3,
     & 3H   F6.1,2H  F5.1,1H F7.2,1H F7.2,1H F7.2,1H F7.2)
C      WRITE(6,1154) 
 1154 FORMAT('     N    V    VAZ   VEL   T1    T2    T3     T1AZ   T1EL
     &  T2AZ   T2EL   T3AZ   T3EL   Q     PHIQ  EPSQ  TPAR  TSEN   QPAR
     & QSEN'/'  CM**-3 KM/S  GRAD  GRAD 1000K 1000K 1000K   GRAD   GRAD
     &   GRAD   GRAD   GRAD   GRAD 10**-4  GRAD  GRAD     1000K        1
     &0*-4   '/'
     &                            ERG/CM**2/S                        ERG
     &/CM**2/S')
C      WRITE(6,*) 'PROTONEN'
C      WRITE(6,*) 'ALPHA-PARTICLES  (Q IN 10**-5 ERG/CM**2/SEC) '
      T1 = OUT(11,3)/1000.
      T2 = OUT(12,3)/1000.
      T3 = OUT(13,3)/1000.
      Q =  OUT(3,4)*1.E4
      TPAR = OUT(8,4)/1000.
      TSEN = OUT(9,4)/1000.
      QPAR = OUT(6,4)*1.E4
      QSEN = OUT(7,4)*1.E4
C       WRITE(6,1156) ( OUT(I,3),I=1,4),T1,T2,T3,(OUT(I,3),I=14,17),
C     & OUT(1,4),OUT(2,4),Q,OUT(4,4),OUT(5,4),TPAR,TSEN,QPAR,QSEN
C 1156  FORMAT(F7.2,F6.1,2(F6.2),3(F6.0)/6(F7.1),F7.2,2(F6.1),
C     & 2(F6.0),2(F7.2)/)
       WRITE(6,*) T1,T2,T3,TPAR,TSEN,Q,QPAR,QSEN
C 1156 FORMAT(5(F9.2)/3(F9.2)/)
C
      T1 = OUT(11,5)/1000.
      T2 = OUT(12,5)/1000.
      T3 = OUT(13,5)/1000.
      Q =  OUT(3,6)*1.E5
      TPAR = OUT(8,6)/1000.
      TSEN = OUT(9,6)/1000.
      QPAR = OUT(6,6)*1.E5
      QSEN = OUT(7,6)*1.E5
C       WRITE(6,1158) ( OUT(I,5),I=1,4),T1,T2,T3,(OUT(I,5),I=14,17),
C     & OUT(1,6),OUT(2,6),Q,OUT(4,6),OUT(5,6),TPAR,TSEN,QPAR,QSEN
C 1158  FORMAT(F7.2,F6.1,2(F6.2),3(F6.0)/6(F7.1),F7.2,2(F6.1),
C     & 2(F6.0),2(F7.2)/)
       WRITE(6,*) T1,T2,T3,TPAR,TSEN,Q,QPAR,QSEN
C 1158 FORMAT(5(F9.2)/3(F9.2)/)

C
C      OUT(16,7) = NSPEK
C
C     CALL PAROUT
C
C      DO 1176 I=1,17
C      WRITE(6,1179)  ( OUT(I,M) , M=1,8 )
C 1179 FORMAT(8(2X,E10.4))
C 1176 CONTINUE
C      DO 1177 II=1,17
C      DO 1177 KK=1,8
C 1177 OUT(II,KK) = 0.
C
      mcyc = mcyc - 1
      WRITE(6,9898) mcyc
 9898 Format(' Next spectrum; mcyc=',I4//)     
      IF ( mcyc .GT. 0 ) GO TO 100
C
 9999 CLOSE(2)
      STOP
      END
C