C boc**********************************************************

 
        subroutine iplot(ndev)


c            in " i-plot.for " to plot the distribution
c                functions
c             link with ,dislin/lib	
c************************************************************* 

c*********************************************************
c************************************************************
      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 intc     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


C     Common blocks for transfer of data to subroutines




      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/ 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


C
      COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
      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

C
      DIMENSION EV(33),EEL(9),EAZ(9),FALF(9,9,9),PARA(3),MAXALF(3),R(9)
C
      INTEGER Hhelios, MAX3D, MAXINT, MAXI1B     
      DATA EV/33*0./
      DATA EEL/9*0./
      DATA EAZ/9*0./
c*********************************************************
c***********************************************************

      LOGICAL HDM,TEST,IPL,ALFA,EXTRA
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      DOUBLE PRECISION F(33,9,9),V(33),AZ(9),EL(9),DF

	common/channel/v,az,el
	dimension s(4)
	double precision s
cc	real*8 s

          character*60 ltit1,ltit2,lx,ly,lz,lv
	common/time/iyear,iday,ihour,minut,isec
        common/ sise /iel,iaz,ien,lel,laz,len

      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3


      common/ MODES /  mSec,seq,Year, DOY, Helios,Mode, BitRate, 
     &                  Format, DisMod, shift, AzShift, I1A3,I2AB
cc
cc                                                           	
	  REAL ct2(110,52),xv(110),yel(52)       ! for plot s(4) in a fram
	  real ct3(110,20,52)	                 ! of xv and elevation angel
	 real angle(2,16),alp(16)                ! yel for diferent azimuth
                                                 ! angle alp

          real vbx(100),vby(100),fbx(100,100),vbx2(2),vby2(2)
cc                              ! for plot s(4) in Vbx and Vby plane

	  real v1(9),v2(100),v3(100),fbo(100,100,9) 
cc                  ! for plot s(4) in plans perpendicular to Bo

           common/maximum/ vxs,vys,vzs,emax
	

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


	write(4,*)'vxs=', vxs,' vys=',vys,' vzs=', vzs
	write(4,*)'emax=',emax

c                         s(1)    V in cm/km (200-1300),nv
c                         s(2)    alpha Azimut angle (-60,40)degree,nal
c                         s(3)    elevation angle (-25,25)degree,nel   


	write(ltit2,902)'icount=', icount,'ipl=',ipl,'alfa=',alfa
 902    format(a7,i2,2x,a4,l2,2x,a5,l2)
	write(4,*)ltit2
	

	write(4,*) '*** iplut ***'
	write(4,*)'##### ipl=', ipl

c	write(4,*)'s(1)=',s(1),' s(2)=',s(2),' s(3)=',s(3)
c	write(4,*)
c	write(4,*)'s(4)=',s(4)
	
c	ly='H2'
c	lz='AU'
c	write(ltit2,101) ly, iy,idoy,ih,min, Heldissun,lz
c 101    format(a2,i3,i4,i3,f6.2,f7.2,a2)
c	write(lv,102)  'Vp=',  Vpi1a
c 102    format (a4,f7.1)
c	write(6,*) ltit2,lv
	
c******************************************


	nv=110
	nel=51
	write(4,*)'iaz =',iaz, ' lc2 (number od data points)=', lc2

	do k=1,lc2
          kaz=iaz+k-1
           alpha=az(k)
	  alp(k)=alpha
	  s(2)=alpha

c	  write(4,*)'ish=',ish,' kaz=', kaz,' alpha=',alpha
	do i=1,51
	  inel=i-26                                         

	  yel(i)=float(inel)
	    
	     do j=1,110
	
	     xv(j)=  200. +float(j) * 10.

	     s(1)=xv(j)
	     s(1)=s(1)*10.**5

	     s(3)=yel(i)
	
	     call cint(s,ipl,ierr)
	     
             if(ierr.eq.1) then
                 ct2(j,i)= 10.**(-20)
	         ct3(j,k,i)=10.**(-20)
              end if
	     if(ierr.eq.0) then
	       cc=dexp(s(4))
	       ct2(j,i)=cc*10.**20              ! in unit of 10**(-20)

	       ct3(j,k,i)=cc*10.**20

	     end if

           end do
	end do
	end do
c************************************************



c	lx='iaz ='
	write(4,103)iyear,iday,ihour,minut,isec
	write(ltit1,103)iyear,iday,ihour,minut,isec
 103    format(5(i5,x))
	
	write(4,*) ltit1

c** plot contours in the frame, v and elevation engle for 9 azimuth angles   

c@         call multicontour(alp,laz,xv,nv,yel,nel,ct3,ndev,ltit1,ltit2)

cc  plot contours in the frame of elevation and azimuth angles for 9 v's

c@        call interenmm(alp,lc2,xv,nv,yel,nel,ct3,ndev,ltit1,ltit2)

	 kg=2
	if(kg.eq.1)stop
c**********************************************************************

c to plot contours in the frame of Vx and Vby = (ex X eB X ex). V
cc	real vbx(110),vby(110),fbx(110,110),vbx2(2),vby2(2)

cc plotbx
c                                    V cos(e) cos(alp) =Vx    
c                                    V cos(e) sin(alp) =Vy
c                                    V sin(e)          =Ve



          Bper=-sqrt(by**2 +bz**2)     ! + or - ??????

	 vbx2(1)=-300.
	 vby2(1)=vbx2(1)*(Bper/Bx)
		
	 vbx2(2)=300.
	 vby2(2)=vbx2(2)*(Bper/Bx)

	do i=1,100
	  vbx(i)= -300. + float(i)*6.

	     do j=1,100
	
	     vby(j)= -300. +float(j) * 6.

	     vx=           vbx(i) + vxs*10.**(-5)
	     vy= (By/Bper)*vby(j) + vys*10.**(-5)
	     vz= (Bz/Bper)*vby(j) + vzs*10.**(-5)
             vm=sqrt(vx**2 + vy**2 + vz**2)		  
	     
             alpha= atand (vy/vx)
	     eps = asind (vz/vm)
c@	 write(4,*)' Vbx Vby plot alpha=',alpha,' epsilon=',eps,' vm=',vm

	     s(1)=vm
	     s(1)=s(1)*10.**5
	     s(2)=alpha
	     s(3)=eps
	
	     call cint(s,ipl,ierr)
	     
             if(ierr.eq.1) then
                 fbxij= 10.**(-30)
            end if
	     if(ierr.eq.0) then
	       cc=dexp(s(4)) 
cc	       fbxij=cc*10.**20              ! in unit of 10**(-20)
	       fbxij=cc
	     end if
	       fbx(i,j)=fbxij/emax           ! in unit of maximum
	       fbx(i,j)= alog10(fbx(i,j))
           end do
	end do
	n=100
	


c@        call bxcontour(vbx2,vby2,vbx,n,vby,n,fbx,ndev,ltit1,ltit2)


	write(4,*)'@@@@@ iplot bxcontour @@@@'
	kg=2
	if(kg.eq.1) go to 100
c**************************************end of plot in plane Vx and Vby
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c********************************** plot in plans perpendicular to Bo.

C                          Bo=Bo.e1
c                          e2=ez x e1
c                          e3=e1 x e2
c     
c  vx = v1 cos aplhab cos phib + v2 (-sin phib) + v3 (-sin alphab cos phib) 
c                          
c  vy = v1 cos alphab sin phib + v2  cos phib   + v3 (-sin alphab sin phib) 
c
c  vz = v1 sin alphab                           + V3 ( cos alphab )

c
cplotbo

	alpb=atand(bz/sqrt(bx**2+by**2))
	phib=atand(by/bx)
c************************************************************
c V1,v2,v3 > Vm,alpha,phi

	vxm=vxs*10.**(-5)
	vym=vys*10.**(-5)
	vzm=vzs*10.**(-5)
             vmm=sqrt(vxm**2 + vym**2 + vzm**2)		  
	write(4,*)'vmm=',vmm,'alpb=',alpb,'phib=',phib

	v1c=0.
	v2c=0.
c	v3c=-125.
        v3c=0.
         vxc = v1c*cosd(alpb)*cosd(phib) + v2c*(-sind(phib)) 
     &                        + v3c*(-sind(alpb)*cosd(phib)) 
                          
         vyc = v1c*cosd(alpb)*sind(phib) + v2c*cosd(phib)   
     &                        + v3c*(-sind(alpb)*sind(phib)) 

         vzc = v1c*sind(alpb)                           
     &                        + V3c* cosd(alpb)


	vxc=vxc +vxs*10.**(-5)
	vyc=vyc +vys*10.**(-5)
	vzc=vzc +vzs*10.**(-5)

             vmc=sqrt(vxc**2 + vyc**2 + vzc**2)		  
             alphac= atand (vyc/vxc)
	     epsc = asind (vzc/vmc)
	write(4,*)'v3c=0, vmc=',vmc,' alphac=',alphac,' epsc=',epsc

c*******************************************************************
	do k=1,7
	
	v1(k)=50.*float(k-4)+0.1	
	do i=1,100
	  v2(i)= -300.1 + float(i)*6.

	     do j=1,100
	
	     v3(j)= -300.1 +float(j) * 6.

         vx = v1(k)*cosd(alpb)*cosd(phib) + v2(i)*(-sind(phib)) 
     &                        + v3(j)*(-sind(alpb)*cosd(phib)) 
                          
         vy = v1(k)*cosd(alpb)*sind(phib) + v2(i)*cosd(phib)   
     &                        + v3(j)*(-sind(alpb)*sind(phib)) 

         vz = v1(k)*sind(alpb)                           
     &                        + V3(j)* cosd(alpb)


	vx=vx +vxs*10.**(-5)
	vy=vy +vys*10.**(-5)
	vz=vz +vzs*10.**(-5)

             vm=sqrt(vx**2 + vy**2 + vz**2)		  
             alpha= atand (vy/vx)
	     eps = asind (vz/vm)

	     s(1)=vm
	     s(1)=s(1)*10.**5
	     s(2)=alpha
	     s(3)=eps
	
	     call cint(s,ipl,ierr)
	     
             if(ierr.eq.1) then
                 fboij= 10.**(-30)
            end if
	     if(ierr.eq.0) then
	       cc=dexp(s(4)) 
cc	       fboij=cc*10.**20              ! in unit of 10**(-20)
	       fboij=cc
	     end if
	       fbo(i,j,k)=fboij/emax           ! in unit of maximum
	       fbo(i,j,k)= alog10(fbo(i,j,k))
           end do
	end do
	end do

	n=100
	

         call bomultic(v1,7,v2,n,v3,n,fbo,ndev,ltit1,ltit2)


	write(4,*)'@@@@@ iplot bxcontour @@@@'
	kg=12
	if(kg.eq.1) stop
C***************************************
c*****************************end of plot in plans perpendicular to Bo.


c        call oneddf(k,iaz,mel,men,icount,ndev,ltit1,ltit2)

c        call surface(k,iaz,mel,men,icount,ndev,ltit1,ltit2)



cc         call multisurface(k,iaz,mel,men,icount,ndev,ltit1,ltit2)

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

	
  100   continue 

	return
	end


c******************************************************************
C########################################################
        subroutine interenmm(alp,laz,xv,nv,yel,nel,ct3,ndev
     &                          ,ltit1,ltit2)
c      plote after interpretation, in the frame of alpha and epsilon
c                   for diferent velocity values

c********************************************************************
        character*60 ltit1,ltit2,liv

	  real alp(16),xv(110),yel(52)
	 real ct3(110,20,52),ct2(9,51)

	 dimension iv(9)

cc	kv=5                ! high speed
	kv=1                ! low speed
	if(kv.eq.5)then	
	  iv(1)=33	
	  iv(2)=37	
	  iv(3)=41	
	  iv(4)=46	
	  iv(5)=51	
	  iv(6)=56
	  iv(7)=61
	  iv(8)=68
	  iv(9)=75
	else
	  iv(1)=10	
	  iv(2)=13	
	  iv(3)=15	
	  iv(4)=18	
	  iv(5)=21	
	  iv(6)=24
	  iv(7)=27
	  iv(8)=31
	  iv(9)=35
        end if


	write(4,*)' e mu con ct3(20,3,10', ct3(20,3,10) 
	write(4,*)'sub int ***** alp,laz=',laz
cc	write(4,756)(alp(j),j=1,9)

	write(4,*)'***** xv,nv=',nv
cc	write(4,756)(xv(i),i=1,9)

	write(4,*)'sub int ***** yel,nel=',nel
cc	write(4,756)(yel(k),k=1,9)
 756    format(9(2xe10.4))	
	
        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
        call pagera
        call complx


	call titlin(ltit1,3)
	call titlin(ltit2,2)

        call name(' azimuth','x')
        call name('elvation','y')

        call intax

	call height(30)


c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c##	nya=2650                        ! for 7 sub-figures
	nya=2800                        ! for 9 sub-figures

	do igrap=1,9
	nyai=nya-(igrap-1)*290
	call axspos(300,nyai)
	call axslen(290,290)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if
        call graf(-30.,20.,-30.,5.,-25.,25.,-25.,5.)
     

c@@@
	k=igrap

	ivv=int(xv(iv(k)))
	write(liv,101)  ivv,'km/s,after cinti'
 101	format(i4,x,a17)

	write(4,*)ivv
	nym=nyposn(6.)
	nxm=nxposn(30.)
	call messag(liv,nxm,nym)

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

	do i=1,nel
	    
	     do j=1,laz
	
	       ct2(j,i)=ct3(iv(k),j,i)              ! in unit of 10**(-20)
c                                            ! f(en,asimuth,elvation)
	end do
	end do

  
       do i=5,15
          a=(float(i)/3.-5.)
	if(i.gt.11) call solid
	if(i.gt.8.and.i.le.11) call dot
	if(i.ge.5.and.i.le.8) call dash
c	do i=1,9
c	a=(float(i)-9.)
          zlev = 10.**a
c@	call conmat(ct2,9,9,zlev)
         call contur(alp,laz,yel,nel,ct2,zlev)          

       end do
	call solid
       call height(30)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	    
	     if (igrap.eq.9) then

	      CALL TITle
	     end if

        	CALL ENDGRf

	end do
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c                  end of enmmlti after cinti (interpretation)


cBo
c************************************************************
         subroutine bomultic(v1,n1,v2,n2,v3,n3,fbo,ndev,ltit1,ltit2)

c       plot in a plane perpendicular to Bo

c********************************************************************
        character*60 ltit1,ltit2,ltit3,liaz

         REAL v1(9),v2(100),v3(100),x1(2),y1(2),x2(2),y2(2)
	  real fbo(100,100,9),ct2(100,100)	
c	write(4,*)'***** bomulti, n1=',n1
	ltit3='ions.c;i_plot.f;bomul;cut per Bo'

	x1(1)=-200.
        y1(1)=0.

	x1(2)=200.
	y1(2)=0.

	y2(1)=-200.
        x2(1)=0.

	y2(2)=200.
	x2(2)=0.

	
        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
        call pagera
        call complx

	call titlin(ltit1,3)
	call titlin(ltit2,2)
	call titlin(ltit3,1)

	
        call name(' VLx','x')
        call name('VLy','y')

        call intax
c        call axspos(450,2670)
C        call graf(7.,31.,7.,2.,1.,7.,1.,1.)



c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c	nya=2800                              ! for 9 pannels
c######################### do loop for 9 figures ###################
c	do igrap=1,n1
c	nyai=nya-(igrap-1)*290                 ! for 9 pannels
c	call axspos(300,nyai)
c	call axslen(600,290)

c*

	
	nya=2750                              ! for 7 pannels

	do igrap=1,n1
	nyai=nya-(igrap-1)*350
	call axspos(350,nyai)
	call axslen(350,350)

	call height(20)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if

        call graf(-200.,200.,-200.,100.,-200.,200.,-200.,100.)


c@@@
	ii=int(v1(igrap))
	write(liaz,101)  ii
 101	format(i4)
	nym=nyposn(20.)
	nxm=nxposn(300.)
	call messag(liaz,nxm,nym)


	k=igrap 


cc         REAL v1(9),v2(100),v3(100)
cc	  real fbo(100,100,9),ct2(110,100)	

c***************************************************************
	do i=1,n2
	     do j=1,n3
	
	       ct2(i,j)=fbo(i,j,k)     ! exponet -20.
c         write(4,*)'v1=',v1(k),' v2=',v2(i),' v3=',v3	
c	 write(4,*)'i=',i,' j=', j,' ct2=',ct2(i,j)	     
	     end do
	end do

  

  
       do i=1,9

	if(i.lt.6)then
	call dash
          a=(float(i)-1.)/2.-3.
	end if

	if(i.gt.5) then
	call solid
       a=alog10(float(i-5)*0.2)
	end if

	zlev=a
c	  call labels('float','contur')

          call contur(v2,n2,v3,n3,ct2,zlev)          
	call dot
	call curve(x1,y1,2)
	call curve(x2,y2,2)
	call solid
       end do

 
       call height(30)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	    	write(4,*)'n1=',n1,'n2=',n2,'n3=',n3
	write(4,*)'igrap=',igrap,' n1=',n1
	     if (igrap.eq.n1) then

	      CALL TITle
	     end if
        	CALL ENDGRf
c######################### end of do loop for 9 figures ###################
	end do

c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       

c	CALL ENDGRf
c	CALL DISFIN
c        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
c       end



c******************************************************************
 
        subroutine multisurface(k,iaz,mel,men,iz,ndev,ltit1,ltit2)
                                         
	
c*******************************************************************
        character*60 ltit1,ltit2,liaz


	real x(7),y(32),zmat(7,32)
	integer iz(7,7,32)
		  REAL ct2(110,50),xv(110),yel(52)
	  real ct3(110,50,20)	
	 real angle(2,16) 


C ***       ***   THE NAMBLE IN LPOL SHOULD BE CHANGE ACCORDING
C ***       ***           TO THE CALCULATION          *******


 	n=25
 

!	iel=(ielr -2)+1     ielr 2-8    ielr=iel +1   iel 1-7                   
!	iaz=(iazr -5)+1     iazr 5-11   iazr=iaz +4   iaz 1-7               

	write(4,*)'i_plot, multi ndev=', ndev

        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')   !p vietical, L horizontal 
  	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

	CALL DISINI
	CALL SCMPLX
        call complx

cc	call chncrv('line')

c	call pagera
c        call axspos(200,2400)
c        call axslen(1800,1800)
	
	call intax
c 	call axslen(970,1070)
c 	call axslen(1000,1500)

	call name('y-energy channel','y')
	call name('x-elevation','x')
	call name('counts','z')
	call titlin(ltit1,2)
	call titlin(ltit2,1)


C	call ydezim(-1)
C	call origin(0,-1800)

	call labels('log','z')
	
	call scale('log','z')

c        CALL GRAF(-6.,-2.,-6.,1.,3.,8.,3.,1.)
C  ***           XLEF XRIT      YMIN YMAX 
c	call ydezim(1)
c	call xdezim(-1)
c	call digits(1,'y')
c	call digits(-1,'x')


c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



	call view3d(-5.,-5.,4.,'abs')
        



c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
CC	nya=2050
c*	nya=2650                              ! for 7 pannels
	nya=2900                              ! for 9 pannels

	do igrap=1,7
c*	nyai=nya-(igrap-1)*327                 ! for 7 pannels
c*	call axspos(350,nyai)
c*	call axslen(600,400)

	nyai=nya-(igrap-1)*290                 ! for 9 pannels
	call axspos(300,nyai)
	call axslen(600,290)
c@@@
	write(liaz,101) 'iaz=', igrap
 101	format(a4,i2)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if
	call solid

        CALL GRAF3d(1.,7.,1.,1.,31.,7.,31.,-2.,0.,2.,0.,1.)

	nym=nyai+70
	nxm=250

	call messag(liaz,nxm,nym)

	iaz=igrap
	    

	do iel = 1,mel

	  x(iel)=float(iel)
	 do i=1,n
	   j=i+6

	   y(i)=float(j)
	   ien=j
           zmat(iel,i)=float(iz(iel,iaz,ien))
	   if(zmat(iel,i).lt.1.e-2) zmat(iel,i)=1.e-2
	 end do
	end do
c	write(6,*)((zmat(i,j),j=1,25),i=1,7)

	  call thkcrv(2)
        	call shlsur

         CALL surfce(x,7,y,25,zmat)  ! page 100 dislin 5.0

	  call dot
c         grid in the xy plane
          call grfini(-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.)
          call nograf
          call graf(1.,7.,1.,1.,31.,7.,31.,-2.)
	  call grid(1,1)
          call grffin

c         grid in the yz plane
           call grfini(1.,-1.,-1.,1.,1.,-1.,1.,1.,1.)
          call nograf
          call graf(31.,7.,31.,-2.,0.,4.,0.,1.)

	  call grid(1,1)
          call grffin

c         grid in the xz plane
          call grfini(-1.,1.,-1.,1.,1.,-1.,1.,1.,1.)
          call nograf
          call graf(1.,7.,1.,1.,0.,4.,0.,1.)
	  call grid(1,1)
          call grffin
	
c	call grid3d(1,1,'all')
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!




	     if (igrap.eq.7) then

c	      call messag(ltit1,1000,1000)
c	      call messag(ltit2,1000,700)
	      CALL TITle
	     end if

        	CALL ENDGRf

	end do

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c#############################################################
c************************************************************
        subroutine multicontour(alp,laz,xv,nv,yel,nel,ct3,ndev
     &                          ,ltit1,ltit2)


c********************************************************************
        character*60 ltit1,ltit2,liaz

         REAL alp(16),xv(110),yel(52)
	  real ct3(110,20,52),ct2(110,52)	
	write(4,*)' e mu con ct3(20,3,10', ct3(20,3,10) 
c	write(4,*)'***** alp,laz=',laz
c	write(4,756)(alp(j),j=1,9)

c	write(4,*)'***** xv,nv=',nv
c	write(4,756)(xv(i),i=1,9)

c	write(4,*)'***** yel,nel=',nel
c	write(4,756)(yel(k),k=1,9)
 756    format(9(2xe10.4))	

        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
        call pagera
        call complx


	call titlin(ltit1,3)
	call titlin(ltit2,2)

        call name(' V','x')
        call name('elvation','y')

        call intax
c        call axspos(450,2670)
C        call graf(7.,31.,7.,2.,1.,7.,1.,1.)

	call height(15)


c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c*
c*	nya=2650                              ! for 7 pannels
c*	do igrap=1,laz
c*	nyai=nya-(igrap-1)*327
c*	call axspos(350,nyai)
c*	call axslen(600,400)

	nya=2800                              ! for 9 pannels

	do igrap=1,laz

	nyai=nya-(igrap-1)*290                 ! for 9 pannels
	call axspos(300,nyai)
	call axslen(600,290)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if
        call graf(150.,1300.,200.,200.,-25.,25.,-25.,10.)


c@@@
	ialp=int(alp(igrap))
	write(liaz,101)  ialp
 101	format(i3)
	nym=nyposn(20.)
	nxm=nxposn(220.)
	call messag(liaz,nxm,nym)


	iaz=igrap 


c***************************************************************
	do i=1,nel
	    
	     do j=1,nv
	
c@	       ct2(j,i)=ct3(j,iaz,i)              ! in unit of 10**(-20)
	       ct2(j,i)=alog10(ct3(j,iaz,i))     ! exponet -20.
		     
	end do
	end do

  
       do i=1,9
c#          a=(float(i)-9.)
c#          zlev = 10.**a
	a=(float(i)-1.)-4.
	zlev=a
	  call labels('float','contur')
          call contur(xv,nv,yel,nel,ct2,zlev)          

       end do

       call height(30)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	    
	     if (igrap.eq.laz) then

	      CALL TITle
	     end if

        	CALL ENDGRf

	end do
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c***********************************************************

C########################################################
        subroutine enmmltic(az,laz,v,nv,el,nel,f,ndev
     &                          ,ltit1,ltit2)


c********************************************************************
        character*60 ltit1,ltit2,ltit3,liv

         double precision az(9),el(9),v(33),f(33,9,9)
	  real ct2(9,9),ct3(33,9,9),xv(33),alp(9),yel(9)

	ltit3='enmmltic, befor cinti'
	do i=1,9
	alp(i)=az(i)
	yel(i)=el(i)
	end do
	
	do 707 i=1,nv
c	write(4,*)'###### i= for V(i))', i
	xv(i)=v(i)*10.**(-5)
	do 707 j=1,laz
c	write(4,*)'###### j= for az(j)',j
c	write(4,754)(ct3(i,j,k),k=1,nel)
	do k=1,nel
	ct3(i,j,k)=f(i,j,k)*10.**20
	end do
 707   continue
 754   format(9(2xe10.4))



	write(4,*)' e mu con ct3(20,3,10', ct3(20,3,10) 
	write(4,*)'***** alp,laz=',laz
	write(4,756)(alp(j),j=1,9)

	write(4,*)'***** xv,nv=',nv
	write(4,756)(xv(i),i=1,9)

	write(4,*)'***** yel,nel=',nel
	write(4,756)(yel(k),k=1,9)
 756    format(9(2xe10.4))	
	
        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
c        call pagera
        call complx


	call titlin(ltit1,3)
	call titlin(ltit2,2)
	call titlin(ltit3,1)

        call name(' azimuth','x')
        call name('elvation','y')

        call intax

	call height(30)


c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c##	nya=2650                        ! for 7 sub-figures
	nya=2800                        ! for 9 sub-figures

	do igrap=1,9
	nyai=nya-(igrap-1)*290
	call axspos(300,nyai)
	call axslen(290,290)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if
        call graf(-30.,20.,-30.,5.,-25.,25.,-25.,5.)
c@        call graf(1.,9.,1.,1.,1.,9.,1.,1.)
     

c@@@
	if(nv.eq.33)ko=16
	if(nv.eq.13)ko=1
	if(nv.eq.15)ko=1
	if(nv.eq.9)ko=0
	
	if(nv.ne.33 .and.nv.ne.13.and.nv.ne.15.and.nv.ne.9)then
	write(4,*)'** STOP *for nv not eq to 13, 33, 15, 9 enmmltic'
	stop
	end if
 
	k=igrap+ko
	ivv=int(xv(k))
	write(liv,101) k, ivv
 101	format(i2,2x,i4)

	write(4,*)ivv, k
	nym=nyposn(6.)
	nxm=nxposn(30.)
	call messag(liv,nxm,nym)

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

	do i=1,nel
	    
	     do j=1,laz
	
	       ct2(j,i)=ct3(k,j,i)              ! in unit of 10**(-20)
c                                            ! f(en,asimuth,elvation)
	end do
	end do

  
       do i=5,18
          a=(float(i)/3.-5.)
	if(i.gt.15) call dashl
	if(i.gt.11) call solid
	if(i.gt.8.and.i.le.11) call dot
	if(i.ge.5.and.i.le.8) call dash
c	do i=1,9
c	a=(float(i)-9.)
          zlev = 10.**a
c@	call conmat(ct2,9,9,zlev)
         call contur(alp,laz,yel,nel,ct2,zlev)          

       end do
	call solid
       call height(30)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	    
	     if (igrap.eq.9) then

	      CALL TITle
	     end if

        	CALL ENDGRf

	end do
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c                  end of enmmlti
C***************************************************************

C########################################################
        subroutine mmlticontour(ct2,az,laz,v,nv,el,nel,f,ndev
     &                          ,ltit1,ltit2)


c********************************************************************
        character*60 ltit1,ltit2,ltit3,liaz

         double precision az(9),el(9),v(33),f(33,9,9)

	  dimension ct2(nv,nel)
           real ct3(33,9,9),xv(33),alp(9),yel(9)

	ltit3='mmiticontour'
	do i=1,9
	alp(i)=az(i)
	yel(i)=el(i)
	end do
	
	do 707 i=1,nv
c	write(4,*)'###### i= for V(i))', i
	xv(i)=v(i)*10.**(-5)
	do 707 j=1,laz
c	write(4,*)'###### j= for az(j)',j
c	write(4,754)(ct3(i,j,k),k=1,nel)
	do k=1,nel
	ct3(i,j,k)=f(i,j,k)*10.**20
	end do
 707   continue
 754   format(9(2xe10.4))



	write(4,*)' e mu con ct3(20,3,10', ct3(20,3,10) 
	write(4,*)'***** alp,laz=',laz
	write(4,756)(alp(j),j=1,9)

	write(4,*)'***** xv,nv=',nv
	write(4,756)(xv(i),i=1,9)

	write(4,*)'***** yel,nel=',nel
	write(4,756)(yel(k),k=1,9)
 756    format(9(2xe10.4))	
	
        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
c        call pagera
        call complx


	call titlin(ltit1,3)
	call titlin(ltit2,2)
	call titlin(ltit3,1)

        call name(' V','x')
        call name('elvation','y')

        call intax
c        call axspos(450,2670)
C        call graf(7.,31.,7.,2.,1.,7.,1.,1.)

	call height(30)


c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

c##	nya=2650                        ! for 7 sub-figures
	nya=2800                        ! for 9 sub-figures

	do igrap=1,laz
	nyai=nya-(igrap-1)*290
	call axspos(300,nyai)
	call axslen(600,290)

	call setgrf('ticks','name','ticks','ticks')

	if(igrap .eq. 1) then
	call setgrf('name','name','ticks','ticks')
	end if
        call graf(150.,1300.,200.,200.,-25.,25.,-25.,10.)


c@@@
	ialp=int(alp(igrap))
	write(liaz,101)  ialp
 101	format(i3)

	write(4,*)ialp, liaz
	nym=nyposn(20.)
	nxm=nxposn(220.)
	call messag(liaz,nxm,nym)


	iaz=igrap 


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


	do i=1,nel
	    
	     do j=1,nv

            ct2(j,i)=ct3(j,iaz,i)              ! in unit of 10**(-20)
   
	end do
	end do

  
       do i=5,18
          a=(float(i)/3.-5.)
	if(i.gt.15) call dashl
	if(i.gt.11) call solid
	if(i.gt.8.and.i.le.11) call dot
	if(i.ge.5.and.i.le.8) call dash
          zlev = 10.**a

        call contur(xv,nv,yel,nel,ct2,zlev)

       end do
	call solid
       call height(30)

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	    
	     if (igrap.eq.laz) then

	      CALL TITle
	     end if

        	CALL ENDGRf

	end do
c@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       

c	CALL ENDGRf
c	CALL DISFIN
c        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
c       end




c******************************************************************
 
        subroutine surface(k,iaz,mel,men,iz,ndev,ltit1,ltit2)
                                         
	
c*******************************************************************
        character*60 ltit1,ltit2

	REAL  X(7),Y(32),zmat(7,32)
        integer iz(7,7,32)
	
C ***       ***   THE NAMBLE IN LPOL SHOULD BE CHANGE ACCORDING
C ***       ***           TO THE CALCULATION          *******


 	n=25
 

!	iel=(ielr -2)+1     ielr 2-8    ielr=iel +1   iel 1-7                   
!	iaz=(iazr -5)+1     iazr 5-11   iazr=iaz +4   iaz 1-7               



        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4p')
	call metafL('cons')
        ELSE
	CALL SETpag('da4p')
	call metafL('gksl')
        END IF
	CALL DISINI
	CALL SCMPLX
        call complx

cc	call chncrv('line')

c	call pagera
        call axspos(200,2400)
        call axslen(1800,1800)
	
	call intax
c 	call axslen(970,1070)
c 	call axslen(1000,1500)

	call name('y-energy channel','y')
	call name('x-elevation','x')
	call name('counts','z')
	call titlin(ltit1,3)
	call titlin(ltit2,2)

C	call ydezim(-1)
C	call origin(0,-1800)

	call labels('log','z')
	
	call scale('log','z')

c        CALL GRAF(-6.,-2.,-6.,1.,3.,8.,3.,1.)
C  ***           XLEF XRIT      YMIN YMAX 
c	call ydezim(1)
c	call xdezim(-1)
c	call digits(1,'y')
c	call digits(-1,'x')


	call view3d(-5.,-5.,4.,'abs')
        
        CALL GRAF3d(1.,7.,1.,1.,31.,7.,31.,-2.,0.,4.,0.,1.)

	call height(50)

	write (6,*) 'mel',mel
	do iel = 1,mel

	  x(iel)=float(iel)
	 do i=1,n
	   j=i+6

	   y(i)=float(j)
	   ien=j
           zmat(iel,i)=float(iz(iel,iaz,ien))
	   if(zmat(iel,i).lt.1.e-2) zmat(iel,i)=1.e-2
	 end do
	end do
c	write(6,*)((zmat(i,j),j=1,25),i=1,7)

	  call thkcrv(2)
	call shlsur

          CALL surfce(x,7,y,25,zmat)  ! page 100 dislin 5.0

	  call dot
c         grid in the xy plane
          call grfini(-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.)
          call nograf
          call graf(1.,7.,1.,1.,31.,7.,31.,-2.)
	  call grid(1,1)
          call grffin

c         grid in the yz plane
           call grfini(1.,-1.,-1.,1.,1.,-1.,1.,1.,1.)
          call nograf
          call graf(31.,7.,31.,-2.,0.,4.,0.,1.)

	  call grid(1,1)
          call grffin

c         grid in the xz plane
          call grfini(-1.,1.,-1.,1.,1.,-1.,1.,1.,1.)
          call nograf
          call graf(1.,7.,1.,1.,0.,4.,0.,1.)
	  call grid(1,1)
          call grffin
	
c	call grid3d(1,1,'all')

	CALL TITle


	CALL ENDGRf
	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

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

c******************************************************************
 
        subroutine oneddf(k,iaz,mel,men,iz,ndev,ltit1,ltit2)
                                         
	
c*******************************************************************
        character*60 ltit1,ltit2

	REAL  X(32),Y(32)
        integer iz(7,7,32)
	
C ***       ***   THE NAMBLE IN LPOL SHOULD BE CHANGE ACCORDING
C ***       ***           TO THE CALCULATION          *******


 	n=25
 

!	iel=(ielr -2)+1     ielr 2-8    ielr=iel +1   iel 1-7                   
!	iaz=(iazr -5)+1     iazr 5-11   iazr=iaz +4   iaz 1-7               



        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4L')
	call metafL('cons')
        ELSE
	CALL SETpag('da4l')
	call metafL('gksl')
        END IF
	CALL DISINI
	CALL SCMPLX

	call chncrv('line')
	
	call intax
 	call axslen(970,1070)
c 	call axslen(1000,1500)

	call name('energy channel','x')


	call name('Counts','y')

c	call labels('log','x')
	call labels('log','y')
	
C	call ydezim(-1)
C	call origin(0,-1800)

	call titlin(ltit1,3)
	call titlin(ltit2,2)

c	call scale('log','x')
	call scale('log','y')

c        CALL GRAF(-6.,-2.,-6.,1.,3.,8.,3.,1.)
C  ***           XLEF XRIT      YMIN YMAX 
c	call ydezim(1)
c	call xdezim(-1)
c	call digits(1,'y')
c	call digits(-1,'x')

        CALL GRAF(8.,32.,8.,2.,0.,4.,0.,1.)

	write (6,*) 'mel',mel

	do iel = 1,mel

	 do i=1,n
	   j=i+7
	   x(i)=float(j)
c	     write(6,*)iz(iel,iaz,ien)
	   ien=j
           y(i)=float(iz(iel,iaz,ien))
	   if(y(i).lt.1.e-2) y(i)=1.e-2
	 end do

	  call thkcrv(2)
          CALL CURVE(X,Y,N)
	
	end do

	maz=7
	 do i=1,n
	j=i+7
	   x(i)=float(j)
	c=0.
	do iaz1 =1,maz
	do iel = 1,mel
	   	   
	   ien=j
           c= c + float(iz(iel,iaz1,ien))
	 end do
	end do

           y(i)= c
	   if(y(i).lt.1.e-2) y(i)=1.e-2

	end do

	  
	  call incmrk(1)
	  call marker(15)
	  call hsymbl(25)

	call dot
          CALL CURVE(X,Y,N)

	CALL TITle


	CALL ENDGRf
	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')
	RETURN
	END

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

c********************************************************************	   
 
       subroutine bxcontour(bx,by,x,n,y,n1,z,ndev,ltit1,ltit2)

c********************************************************************
          character*60 ltit1,ltit2

	  REAL z(100,100),x(100),y(100),bx(2),by(2)
	

c	 program ex14_1
c         parameter (n=100)


        IF(NDEV.EQ.1) THEN
	CALL SETpag('da4L')
	call metafL('cons')
        ELSE
	CALL SETpag('da4l')
	call metafL('gksl')
        END IF

c	call setpag('da4p')
        call disini
        call pagera
        call complx


	call titlin(ltit1,3)
	call titlin(ltit2,1)

        call name(' Vx km/s','x')
        call name('Vby','y')

        call intax
c        call axspos(450,2670)
	call axslen(1000,1000)
        call graf(-300.,300.,-300.,100.,-300.,300.,-300.,100.)

	call height(30)
  



  
       do i=1,9

c@          a=(float(i)-9.)
c@          zlev = 10.**a
c@	  call labels('log','contur')
	if(i.lt.6)then
	call dash
          a=(float(i)-1.)/2.-3.
	end if

	if(i.gt.5) then
	call solid
       a=alog10(float(i-5)*0.2)
	end if
	zlev=a
	  call labels('float','contur')

          call contur(x,n,y,n,z,zlev)          
       end do

	call curve(bx,by,2)
       call height(50)
       call title
       

	CALL ENDGRf
	CALL DISFIN
        IF(NDEV.NE.1) CALL SYMFIL('psc2','DELETE')

	return
       end