program sodcon

!  g77 -g -o sodcon sodcon.f

      implicit none

      integer i,j,k,ni,nj,nk,kk,iq0,nbelow0
      integer nvmx,npmx,ntype,angle_type,nscale
      real scalef,phi2,pmx,xii,dx,xsize,ysize,xcen,ycen
      real xsta,ysta,red,green,blue,xss,yss,xsss,fnum
      real txdiff,tydiff,rtmp9,xx1,xx2,yy1,yy2,height
      real ff1,ff2,txold,tyold,y_crd_sc,y_crd_fl,scx
      integer nsec,nxs,nys,ijk,jkl

c	This program reads in .SOD and .COD files and plot these
c       in Bunge-convention ODF space.  At current stage, sections
c       are made every 5 degrees (phi2 for .COD, and phi1 for .SOD).

c      This program is made available on the understanding that its
c      use will be appropriately acknowledged.  No warranty of any
c      kind is avaiable.  Please inform tony rollett of any errors that
c      you find, rollett@andrew.cmu.edu
c                  ADR   March 2004

!  last edited for HODF input,  Apr 07

!  minor correction to allow for each section going to 60 degrees (HEX), iv 08

	real fint(20,20,20), vodf(20,20,20)
	real scint(0:6),scleg(0:6),scleg2(0:6)
	integer nvint(20,20,20)
        integer afile_len
        parameter (afile_len = 50)
	character afile*50, file*50, out1*50
        character outfile*50, outfile_jpg*50
       
	real array(19,19),cval(0:10),xold,yold,xdiff,ydiff,xorig,yorig
	real colrval(3,0:10)
c  colors are in COLRVAL
	integer nval   !  number of contour values
	character*5 nomen,seclab,label
	real delth,rm,delom,pm
	integer iw,jw,iper(3),iavg,ngr
	logical rescale
	real data_min,data_max,data_delta,auto_val(6)
	logical l_com_line,l_ticks,l_lines,l_hcp
	integer q_contours
	character arg_line*20,ticks*1,lines*1,hcp*1
	real xclip(5),yclip(5)
	integer nclip  !  for clipping off hcp boxes
	integer ispec,ioffpp,spvall,ilegg,ilabb,nhii,ndeccn,nlbll,
     1     mscall,ldsh,hgtlab
	common/conpar/ispec,ioffpp,spvall,ilegg,ilabb,nhii,ndeccn,nlbll,
     1     mscall,ldsh,hgtlab
        integer point
        character inline*130
        real xa,ya
        integer upper

!  CODE::
	print *,' '
	print *,' '
	print *,'     ***********************************************'
	print *,'     *                                             *'
	print *,'     *     .COD and .SOD plotting in ODF space     *'
	print *,'     *            by Paul Seungyong Lee            *'
	print *,'     *            modified by Tony Rollett         *'
	print *,'     *                                             *'
	print *,'     *         Carnegie Mellon University          *'
	print *,'     *                                             *'
	print *,'     ***********************************************'
	print *,' '

        if(iargc().lt.1) then
           print*,'Usage for SODCON:'
           print*,
     &      './sodcon input_data contours ticks lines hcp [1=yes]'
        endif
	l_com_line = .false.
        if ( iargc().ge. 1 ) then
           CALL GetArg(1,afile)
           if(afile.ne.'') l_com_line = .true.
           print*,'reading input from ',afile
        else
	   print *,'Name of input file (.SOD, .SMH, .HODF or .COD) ?'
	   read(*,1011) afile
 1011	   format(a)
        end if
!	print*,'debug: afile = ',afile
!  detect an argument for the input file name; look for user input if not present

        q_contours = -99
        if ( iargc().ge. 2 ) then
           CALL GetArg(2,arg_line)
!	print*,'debug: arg_line = ',arg_line
           read(arg_line,"(i4)") q_contours
           print*,' contours choice = ',q_contours
        end if

        ticks = ''
        l_ticks = .false.
        if ( iargc().ge. 3 ) then
           CALL GetArg(3,ticks)
           if(ticks.eq.'1') l_ticks = .true.
           print*,'ticks on ? ', l_ticks
        else
	   print*,'Tick marks?  1 = YES '
	   read(*,"(a)") ticks
	   if(ticks.eq.'1') l_ticks = .true.
        end if
!	print*,'debug: ticks = ',ticks

        lines = ''
        l_lines = .false.
        if ( iargc().ge. 4 ) then
           CALL GetArg(4,lines)
           if(lines.eq.'1') l_lines = .true.
           print*,'lines on ? ', l_lines
        else
c  ***  found mistake in line above 4 july 05: was l_ticks instead of l_lines
	   print*,'Line contours (vs solid)?  1 = YES '
	   read(*,"(a)") lines
	   if(lines.eq.'1') l_lines = .true.
         end if
!	print*,'debug: lines = ',lines

         hcp = ''
         if ( iargc().ge. 5 ) then
            CALL GetArg(5,hcp)
            if(hcp.eq.'1') l_hcp = .true.
            print*,'plotting as HCP ? ', l_hcp
         else
            l_hcp = .false.
            print*,'Limit phi2 to 60 degrees for HCP?  1 = YES '
            read(*,"(a)") hcp
            if(hcp.eq.'1') l_hcp = .true.
         endif

	data_min=9999999.
	data_max=0.
	ispec=1
	ilegg=0
	ilabb=0
	nhii=-1
c       turns off legend, labels, hi/lo labels

	print *,' '
c	print *,'Choose the scaling method.'
c	print *,'    Linear ------------ 1 '
c	print *,'    Log --------------- 2 '
c	read(*,*)nscale
	nscale = 1  !  stick to linear scaling for now

	kk=50
	do 154, i=50,1,-1
c       find the right-most occurrence of a "."
	   if(afile(i:i).eq.'.') then
	      kk=i
	      goto 155
	   endif
 154	continue
 155	continue
	file=afile(kk+1:kk+3)
	if(file.eq.'sod'.or.file.eq.'SOD'.or.file.eq.'smh') then
	   ntype=1
	   goto 888
	endif
	if(file.eq.'cod'.or.file.eq.'COD'.or.file.eq.'cmh') then
	   ntype=2
	   goto 888
	endif
	if(file.eq.'son'.or.file.eq.'SON'.or.file.eq.'snh') then
	   ntype=3
	   goto 888
	endif
	if(file.eq.'con'.or.file.eq.'CON'.or.file.eq.'cnh') then
	   ntype=4
	   goto 888
	endif
        point = 0
        point = index(afile,'.hodf')
        if(point.eq.0) point = index(afile,'.HODF')
	if(point.gt.0) then
	   ntype=5
	   goto 888
	endif
	write(*,*)
     &       'This program digests only HODF,SOD,COD,SON and CON files'
	goto 999

 888	open(11,file=afile,status='old')

!  now we read the first header line to check for ranges
        read(11,*)
        read(11,1332) nomen,delth,rm,delom,pm
        print*, nomen,delth,rm,delom,pm
        rewind(11)
        if ( pm .eq. 60. ) then
           print*,'azimuth goes to 60 degrees - hexagonal assumed '
           l_hcp = .true.
        elseif ( pm .ne. 90. ) then
           write(*,*) 'The range of phi1 is not 90, but =' , pm
           write(*,*) 'therefore this is not suitable for SODCON'
           if ( pm .eq. 360. ) then
              print*,'The range is 360 degrees, so try SOD4CON'
           endif
           stop
        endif

!       detect the type of Euler angles in use
        angle_type=1
        if(nomen(5:5).eq.'B'.or.nomen(4:4).eq.'B') then
           write(*,*) 'detected Bunge angles'
        endif
c       assume Bunge angles
        if(nomen(5:5).eq.'K'.or.nomen(4:4).eq.'K') then
           write(*,*) 'detected Kocks angles'
           angle_type=2
        endif
c       Kocks angles
       
c	write(*,*) 'Enter a name for the OUTPUT file, e.g. sample.cmh.ps'
c	read(*,1011) afile
	outfile = afile(1:(kk+3))//'.ps'
	outfile_jpg = afile(1:(kk+3))//'.jpg'
	write(*,*) 'Using ',outfile,' as output filename'
	call newdev(outfile,3)
c	call newdev(afile,3)
c       
c       
	if(ntype.eq.2.or.ntype.eq.4) then
	   ff1=1
	   ff2=2
        else
           ff1=2
           ff2=1
	endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
c	read in the input data
c       
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c       see below for dealing with header lines

c 	do 21, i=1,2
c	read(11,*)
c       read(11,"(a5)") nomen
c       
c       21	continue
c       
c       read the first header lines

!  @@@@@@@@@@@@@@@@@@
        upper = 19
        if ( ( ntype .eq. 2 .or. ntype .eq. 4 ) .and. l_hcp ) upper = 13
        if(ntype.le.4) then
!  we have popLA input
           
	do 725 k = 1 , upper

	   read(11,*)
	   read(11,1332) nomen,delth,rm,delom,pm,iw,jw
     1	,(iper(i),i=1,3),iavg,ngr,seclab,label
 1332	   format(a5,4f5.1,5i2,2i5,2a5)

	   scalef=100./float(iavg)
	   rescale=(abs(scalef-1.0).gt.0.05)
	   if(rescale) write(*,*) 'Rescaling intensities for k=',k,scalef
c       re-scale all the intensities in the section by the ratio
c       of the nominal average (100) divided by IAVG in the file
c       but require more than a 5% change to do it
c       
	   do 724, j=1,19
	      read(11,31,err=932) (nvint(i,j,k),i=1,19)
 31	      format(1x,19i4)
	      if(rescale) then
		 do 723, i=1,19
		    nvint(i,j,k) = int(float(nvint(i,j,k))*scalef)
 723		 continue
	      endif
c       
 724	   continue
	   read(11,*,end=725)
c       blank line after each section
c       
 725	continue  !  end of input

!  now check for max, min
 	nvmx=nvint(1,1,1)
	do k = 1 , upper
	   do  j = 1,19
	      do i = 1,19
		 if(nvint(i,j,k).le.0.and.nscale.eq.2) nvint(i,j,k)=1
c       reset intensity to 1 if we are using a log scale
		 if(nvint(i,j,k).gt.nvmx) nvmx=nvint(i,j,k)
 	      enddo

	      do i=1,19
		 if(nvint(i,j,k).lt.data_min) data_min = nvint(i,j,k)
		 if(nvint(i,j,k).gt.data_max) data_max = nvint(i,j,k)
	      enddo
	   enddo
 	enddo

c  changed the scaling to have 1.0 = 1 MRD
	do nk = 1 , upper
	   do nj = 1,19
	      do 149, ni = 1,19
!		 if(nvint(ni,nj,nk).gt.npmx) then
!		    nvint(ni,nj,nk)=npmx
!		 endif
		 if(nscale.eq.2) then
		    fint(ni,nj,nk) = log(float(nvint(ni,nj,nk))*0.01)
c       write(*,*)fint(ni,nj,nk)
		 else
                    fint(ni,nj,nk) = float(nvint(ni,nj,nk))*0.01
		 endif
c       write(*,*) fint(i,j,k)
 149	      continue
	   enddo
	enddo

c$$$	do ni = 1,19
c$$$           print*,'section for k = ',ni
c$$$           do nj = 1,19
c$$$              print"(20f5.1)",(fint(ni,nj,nk),nk = 1,19)
c$$$           enddo
c$$$        enddo

        endif  !  if(ntype.le.4) then  [popLA]
!  @@@@@@@@@@@@@@@@@@

! ######################
        if(ntype.eq.5) then
!  we have RESMAT .HODF input; for now, cubic only
           angle_type=1  !  assume Bunge
           rescale = .false.
           nvmx = 0
           npmx = 0
           print*,'Reading header lines '
           do i = 1,18
              read(11,"(a)") inline
              print*,inline
           enddo
           do i = 1,19
              read(11,*) phi2
              print*,'reading section at phi2 = ',phi2,' degrees'
              do j = 1,19
                 read(11,*) (fint(k,j,i),k=1,73)
                 do k = 1,73
                    if(fint(k,j,i).gt.float(nvmx)/100.) then
                       nvmx = int(100.*fint(k,j,i))
                    endif
                    if(fint(k,j,i).gt.float(npmx)/100.) then
                       npmx = int(100.*fint(k,j,i))
                    endif
                    if(fint(k,j,i).lt.data_min/100.) then
                       data_min = fint(k,j,i)*100.
                    endif
                    if(fint(k,j,i).gt.data_max/100.) then
                       data_max = fint(k,j,i)*100.
                    endif
                 enddo   !  k = 
              enddo   !   j = 
           enddo  !  i = 
        endif
! ######################

        print*,' minimum value = ',data_min/100.
     &       ,' max value = ',data_max/100.
	 data_delta=(data_max-data_min)/6.
	 print*,'difference/6. = ',data_delta/100.
	 print*,' automatic levels = '
	 do i=1,6
	    auto_val(i)=(data_min+(float(i-1)*data_delta)
     1   	 +(data_delta*0.5))*0.01
	    print*,' auto contour level ',i,' = ',auto_val(i)
	 enddo

      if(.not.l_com_line) then
	 print*,'Now we deal with contour levels'
	 print*,'NB: values given here are multiples of random (uniform)'
	 print*,'  commonly written as MRD units'
	 print*,'NB: we assume popLA format files, such that the'
	 print*,'intensities are stored as i4 integers, scaled such'
	 print*,'100 == 1MRD, unless IAVG is not equal to 100'
	 print*,'in which case, inspect the code for meaning of'
	 print*,'the re-scaling operation'
	 print*
      endif

	write(*,492) nvmx
 492	format('Maximum intensity =',i7)
c     1  //'Enter the maximum intensity to be plotted.')
	npmx=nvmx+1  !  seems unnecessary now to read this in
c 492	format('Maximum intensity =',i7,/,
c     &         'Put the maximum intensity to be plotted.')
c	read(*,*)npmx
	write(*,*)
	pmx=float(npmx)

c       making scale bar

        nval = 6
!	if(l_com_line) then
        if ( q_contours.gt. 0 ) then
	   iq0 = q_contours
	else
	   write(*,*) 'Automatic choice of levels? 0=yes'
	   write(*,*) 'Enter 1 for 0.5 / 1 / 2 / 4 / 8 / 16 '
	   write(*,*) 'Enter 2 for -0.5 / -0.25 / 0 / 0.25 / 0.50 / 1.0'
	   write(*,*) 'Enter 3 for 0.04 / 0.1 / 0.3 / 0.5 / 0.7 / 0.9 '
	   write(*,*) 'Enter 4 for 0.009 /0.013 / 0.04 / 0.2 / 0.6 /0.9'
	   print*,    'Enter 5 for auto levels above, based on min-max'
	   read(*,*) iq0
	endif

	if(iq0.eq.0) then
	   do 385,i=0,5
	      xii=float(i)/5.
	      scint(i)=xii
	      scleg(i)=(1.-xii)*pmx
	      scleg2(i)=exp((1.-xii)*log(pmx))
c	write(*,*)scint(i),scleg(i),scleg2(i)
	      if(nscale.eq.1) then 
		 cval(i+1)=scleg(i)
	      else
		 cval(i+1)=scleg2(i)
	      endif
 385	   continue
	elseif(iq0.eq.1) then
	   scleg(0) = 0.5
	   cval(1) = 0.5
	   scleg(1) = 1.
	   cval(2) = 1.
	   scleg(2) = 2.
	   cval(3) = 2.
	   scleg(3) = 4.
	   cval(4) = 4.
	   scleg(4) = 8.
	   cval(5) = 8.
	   scleg(5) = 16.
	   cval(6) = 16.
	   nbelow0=1
	   nscale=1
	elseif(iq0.eq.2) then
	   scleg(0)=-0.50
	   cval(1)=-0.50
	   scleg(1)=-0.25
	   cval(2)=-0.25
	   scleg(2)=0.
	   cval(3)=0.
	   scleg(3)=0.25
	   cval(4)=0.25
	   scleg(4)=0.50
	   cval(5)=0.5
	   scleg(5)=1.
	   cval(6)=1.
	   nbelow0=2
	   nscale=1
	elseif(iq0.eq.3) then
	   scleg(0)=0.0400
	   cval(1)=0.0400
	   scleg(1)=0.1000
	   cval(2)=0.1000
	   scleg(2)=0.3000
	   cval(3)=0.3000
	   scleg(3)=0.5000
	   cval(4)=0.5000
	   scleg(4)=0.7000
	   cval(5)=0.7000
	   scleg(5)=0.9000
	   cval(6)=0.9000
	   nbelow0=1
	   nscale=1
        elseif(iq0.eq.4) then
	   scleg(0)=0.0090
	   cval(1)=0.0090
	   scleg(1)=0.0130
	   cval(2)=0.0130
	   scleg(2)=0.0400
	   cval(3)=0.0400
	   scleg(3)=0.2000
	   cval(4)=0.2000
	   scleg(4)=0.6000
	   cval(5)=0.6000
	   scleg(5)=0.9000
	   cval(6)=0.9000
	   nbelow0=1
	   nscale=1
	elseif(iq0.eq.5) then
	   do i=1,6
	      scleg(i-1)=auto_val(i)
	      cval(i)=auto_val(i)
	   enddo
	else
	   nbelow0=0
	   nscale=1
	   write(*,*) 'Enter 6 contour values.... '
	   do 387, i=0,5
	      write(*,*) 'Enter the value for the ',i,'th value'
	      read(*,*) scleg(i)
	      cval(i+1)=scleg(i)
	      if(scleg(i).le.1.0) nbelow0 = nbelow0+1
c       was counting contours below 100, now 1.0;  ADR  dec 01
 387	   continue
	endif

c       end of scale bar information

	colrval(1,0)=0.0
	colrval(2,0)=0.0
	colrval(3,0)=0.0	! black
	colrval(1,1)=0.0
	colrval(2,1)=0.0
	colrval(3,1)=1.0	! blue
	colrval(1,2)=0.0
	colrval(2,2)=1.0
	colrval(3,2)=1.0	! cyan
	colrval(1,3)=0.0
	colrval(2,3)=1.0
	colrval(3,3)=0.0	!  green
	colrval(1,4)=1.0
	colrval(2,4)=1.0
	colrval(3,4)=0.0	! yellow
	colrval(1,5)=1.0
	colrval(2,5)=0.35
	colrval(3,5)=0.0	!  orange
	colrval(1,6)=1.0
	colrval(2,6)=0.0
	colrval(3,6)=0.0	!  red
	colrval(1,7)=1.0
	colrval(2,7)=0.0
	colrval(3,7)=1.0	! magenta
c  add more values if you use more than 6 contour values!

c       
	write(*,32) afile
 32	format('Status : Data read in from ',a)
c       
c       end of data input
c       
	call psinit(.true.)
c       
c       initiating post script plotting
c       
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
c	start of 3 orientation loops (phi1,phi,phi2)
c       phi2 is multiples of 5.
c       phi1 and phi are defined by the resolution of VODF space.
c       
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
	dx=2.5/90.*1.5
c       
c       
	write(*,39)
 39	format('Status : Plugging in data ')
c       
        xsize = 1.5
        ysize = 1.5
	xold=0.
	yold=0.
	xdiff=0.
	ydiff=0.

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
c       vodfmx is the maximum vodf.  vodf will be normalized by vodfmx
c       because gray scale should be in the range of 0 and 1.
c       
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

	do 26 , k = 1 , upper
	   nsec = k-1

c       nsec is the number of ODF section from top left to bottom right.
c       It is from 0 to 18.

	   nys = 5 - nsec/4
c	   nxs=nsec+1
!           nxs = mod(k,4) + 1
           nxs = mod((k-1),4) + 1
c$$$	   if(nxs.gt.16) nxs=nxs-16
c$$$	   if(nxs.gt.12) nxs=nxs-12
c$$$	   if(nxs.gt.8) nxs=nxs-8
c$$$	   if(nxs.gt.4) nxs=nxs-4
!	   write(*,*)'nsec, nxs, nys = ',nsec,nxs,nys

c       nxs is the x-section number from left.
c       nys is the y-section number from top.

c	   xcen = 0.2 + nxs*1.5
	   xcen = -0.5 + nxs*xsize
	   ycen = (nys)*1.9
	   xsta = xcen- (xsize/2.)
	   ysta = ycen+ (xsize/2.)
	   xdiff = xsta-xold
	   ydiff = (ysta-xsize)-yold
	   xold = xsta
	   yold = ysta - xsize
c	   write(*,*)xcen,ycen,xsta,ysta
	   call plot(xdiff,ydiff,-3)


c       xcen, and ycen are the coordinates of centers of ODF sections.
c       These are the same locations as the centers in ODF frames.
c       xsta, and ysta are the coordinates of starting point of each
c       ODF section.
c       
	   do 1020, ijk = 1,19
	      do  jkl = 1,19
		 array(ijk,jkl)=fint(ijk,(20-jkl),k)
c       have to go from bottom to top to obtain the correct order
	      enddo
 1020	   continue
c       copy section into array for contouring
c       
	   if(ntype.eq.1.and.l_hcp) then
	      call defclip_box(xclip,yclip,nclip,(xsize*2./3.),ysize)
	      call GSAV
	      call clipbox(xclip,yclip,nclip)
	   endif
c  set up the restricted box for HCP

	   if(l_lines) then

c	if(nbelow0.gt.0) then
	      ldsh=6
	      red=0.
	      green=0.
	      blue=1.
	      call setcolr(red,green,blue)
	      CALL SETLW (.01)
c	      if(ntype.eq.1.and.l_hcp) then
c		 call conrec(array,13,19,19,xsize,ysize,cval(1),2)
c	      else
		 call conrec(array,19,19,19,xsize,ysize,cval(1),2)
c	      endif
c	endif
	      ldsh=0
	      red=0.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.02)
c	      if(ntype.eq.1.and.l_hcp) then
c		 call conrec(array,13,19,19,xsize,ysize,cval(3),2)
c	      else
		 call conrec(array,19,19,19,xsize,ysize,cval(3),2)
c	      endif

	      ldsh=0
	      red=1.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.03)
c	      if(ntype.eq.1.and.l_hcp) then
c		 call conrec(array,13,19,19,xsize,ysize,cval(5),2)
c	      else
		 call conrec(array,19,19,19,xsize,ysize,cval(5),2)
c	      endif
	   else
c	      if(ntype.eq.1.and.l_hcp) then
c		 call concolr(array,13,19,19,xsize,ysize,cval(0),colrval,
c     +7,0,0.)
c	      else
		 call concolr(array,19,19,19,xsize,ysize,cval(0),colrval,
     +7,0,0.)
c	      endif
	      red=0.
	      green=0.
	      blue=0.
	      call setcolr(red,green,blue)
	      CALL SETLW (.01)
c	      if(ntype.eq.1.and.l_hcp) then
c		 call conrec(array,13,19,19,xsize,ysize,cval(0),6)
c	      else
		 call conrec(array,19,19,19,xsize,ysize,cval(0),6)
c	      endif
	   endif

	   CALL SETLW (.005)
	   if(l_ticks) then
	      if(ntype.eq.1.and.l_hcp) then
		 call border((xsize*2./3.),ysize,1111,1111,1,9,1,9)
	      else
		 call border(xsize,ysize,1111,1111,1,9,1,9)
	      endif
	   else
	      if(ntype.eq.1.and.l_hcp) then
		 call border((xsize*2./3.),ysize,0000,1111,1,9,1,9)
	      else
		 call border(xsize,ysize,0000,1111,1,9,1,9)
	      endif
	   endif

	   if(ntype.eq.1.and.l_hcp) then
	      call GREST
	   endif

c  now to label the subplot
	   ldsh=0
	   red=0.
	   green=0.
	   blue=0.
	   call setcolr(red,green,blue)

	   xss = 0.33
	   if(ntype.eq.1.and.l_hcp) then
	      xss = 0.15
	   endif
	   yss = -0.21
	   if(angle_type.eq.1) then 
	      call grksym(xss,yss,.15,45,0.,1,1)
	      call keknum(xss+.13,yss-.07,.09,ff2,0.,-1,1)
	   endif
	   if(angle_type.eq.2) then
	      if(ntype.eq.2) call grksym(xss,yss,.15,45,0.,1,1)
	      if(ntype.eq.1) call grksym(xss,yss,.15,23,0.,1,1)
	   endif
	   xsss = xss+0.23
	   call keksym(xsss,yss,.15,'  =',0.,3,1)
	   fnum = float(k-1)*5.
c       fnum=float((i+(ky-1)*4)*5-5)
	   call keknum(xsss+.4,yss,.15,fnum,0.,-1,1)

 26     continue                !  end of subplot loop (each box)

        xorig = -1.*xsta - xsize
        yorig = -1.*(ysta - xsize)
c        xorig = xsize
c        yorig = 0.
        call plot(xorig,yorig,-3)
c	   write(*,*) 'resetting origin to ',xorig,yorig
c       
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
c       frame setting
c       
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        write(*,239)
 239    format('Status : completed contour plots')

!        xa = 5.7
        xa = 6.5
        ya = 2.4
        call arrow(xa,ya,xa+0.7,ya,.08,20.,1)
        call arrow(xa,ya,xa,ya-0.7,.08,20.,1)
        if(angle_type.eq.1) then 
           call grksym(xa+0.15,1.85,.15,21,0.,1,1)
           call grksym(xa+0.4,ya-0.25,.15,45,0.,1,1)
           call keknum(xa+0.5,ya-0.37,.09,ff1,0.,-1,1)
        endif
        if(angle_type.eq.2) then
           call grksym(xa+0.15,ya-0.55,.15,8,0.,1,1)
           if(ntype.eq.2) call grksym(xa+0.55,ya-0.25,.15,23,0.,1,1)
           if(ntype.eq.1) call grksym(xa+0.55,ya-0.25,.15,45,0.,1,1)
        endif
        call keknum(xa,ya-1.0,.15,90.,0.,-1,1)
        call keknum(xa+1.0,ya-0.07,.15,90.,0.,-1,1)
        
        y_crd_fl = 0.15
        call setfnt(31)
        call keksym(1.3,y_crd_fl,.14,10hInput file,0.,10,0)
!        call keksym(1.7,y_crd_fl,.14,1h:,0.,1,0)
        call keksym(2.5,y_crd_fl,.14,afile,0.,afile_len,0)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     
c     making scale list
c     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


        if(l_lines) then
           print*,'Status : making scale list'
	   y_crd_sc = 0.45
	   call keksym(1.3,y_crd_sc,.14,12hContours at ,0.,12,0)
	   red=0.
	   green=0.
	   blue=1.
	   call setcolr(red,green,blue)
	   do 386, i = 1,6
	      if(i.gt.2) then
		 red=0.
		 green=0.
		 blue=0.
		 call setcolr(red,green,blue)
	      endif
	      if(i.gt.4) then
		 red=1.
		 green=0.
		 blue=0.
		 call setcolr(red,green,blue)
	      endif
	      scx = 2.5+float(i)/10.*6.5
c     call rectfilg(scx,.31,scx+.7,.31,.3,scint(i))
	      if(nscale.eq.1) then
!		 call keknum(scx+.1,y_crd_sc,.14,scleg(i),0.,3,0)
		 call keknum(scx+.1,y_crd_sc,.14,cval(i),0.,3,0)
	      endif
	      if(nscale.eq.2) then
!		 call keknum(scx+.1,y_crd_sc,.14,scleg2(i),0.,3,0)
		 call keknum(scx+.1,y_crd_sc,.14,cval(i),0.,3,0)
	      endif
 386	   continue
        end if

        red=0.
        green=0.
        blue=0.
        call setcolr(red,green,blue)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       
c       Making frames of ODF space using post script language.
c       Sections were made with constant phi2 with every 5 degrees.
c       
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

!  now to make a legend
           
        if(.not.l_lines ) then
           print*,'Status : making legend'
           xa = 8.0
           ya = 0.9

           call setlw(0.025)
           rtmp9 = 0.25
           call keksym ( xa , ya
     $          , .15 , 'MRD' , 0. , 3, 0 )
           do i = 1 , nval
              xx1 = xa + (xsize/2.)
              xx2 = xx1 + 2.*rtmp9
              yy1 = ya + (ysize*0.01) + float(i)*rtmp9
              yy2 = yy1
              height = rtmp9
!                 call keknum(-0.8+xx1, yy1-(height/3.), 0.15
!     &                , cval(i), 0., 2, 0)
              call keknum(-0.8+xx1, yy1+(height*0.7) , 0.15
     &             , cval(i), 0., 2, 0)
              CALL RECTFILC(XX1,YY1,XX2,YY2,HEIGHT,colrval(1,i)
     &             ,colrval(2,i),colrval(3,i))
              CALL SLDLIN (xX1, yY1, xx1-0.10, yY1, 0.015)
           end do
           yy1 = yy1 + rtmp9
           CALL SLDLIN (xX1, yY1, xx1-0.10, yY1, 0.015)
        end if

        call plotnd

        call system('convert '//outfile//' '//outfile_jpg)
        call system('rm '//outfile)

        goto 999
 932    write(*,*)'Input file is not in the right format.'
 999    stop
        end

c
c ________________________________________
c
	subroutine defclip_box(xclip,yclip,nclip,xsize,ysize)

c  defines the clipping region
c  in this case a box

	implicit none

	real xclip(5),yclip(5)
	integer nclip
	real xsize,ysize

c  CODE::

	nclip = 5

	xclip(1) = 0.
	yclip(1) = 0.
	xclip(2) = xsize
	yclip(2) = 0.
	xclip(3) = xsize
	yclip(3) = ysize
	xclip(4) = 0.
	yclip(4) = ysize
	xclip(5) = 0.
	yclip(5) = 0.

	return
	end

c____________

	include 'psplot.txt'