program plot

c  This program converts 3D and 2D model in typical Harvard format and 
c  converts into a format that can be used to plot the model in GMT.

c  Output file contains
c   longitude  latitude  value
c  where the value is in % for 3D models and in km for 2D models.

cccccccccc
c  Subroutines:
c	kernel is in sub_kernel.f
c	readModel is in sub_readModel.f 
c	switchLatLon in sub_switchLatLon
c	chebyshev_mantle in sub_chebyshev.f 
cccccccccc

c  July 29, 1999.
c  Last Modified:  September 07, 1999.


      real block_size
      parameter (block_size = 3.0)
c  block_size = size of the blocks in degrees for plotting in GMT

      integer lmax, kmax
      parameter (lmax = 6, kmax = 13)
c  lmax = maximum angular degree.
c  kmax = maximum radial order. 

      real Xlm(0:lmax, 0:lmax)
c  Legendre function and factors for a given l,m at a given point.

      integer nj, ni
      parameter (nj = Int(180.0*360.0/(block_size*block_size)) )
      parameter (ni = (lmax+1)*(lmax+1) )
c  nj = number of data points
c  ni = number of spherical coefficients

      real Gkernel(nj, ni)
c  Gkernel is the matrix containing the kernel.
c  Gkernel elements correpond to spherical harmonics at a given point.

      real A(0:lmax,0:lmax,0:kmax), B(0:lmax,0:lmax,0:kmax)
c  model coefficients in real spherical harmonics

      real depth, radius
c-- depth = depth we are interested in
c-- radius = radius of the Earth corresponding to the above depth

      real Tr(0:kmax) 
c  single precisioned Chebyshev polynomials.

      integer i
c  counter to keep track of model position

      real model(ni)
c  model vector containing spherical harmonic coefficients.

      real data(nj), data_GMT(nj)
c  data = vector containing data at various points (fwd modelling result)
c  data_GMT = data vector containing values in GMT order

      logical there
c  this is .true. if kernel is already stored.

      character*200 filename
c  filename to which resulting data in GMT format is saved.

      character*10 ans_depth

      character*45 file_kernel

      real lat, lon
c  lat = latitude
c  lon = longitude

      character*5, ans_changeModel
c  this starts with 'y' if another model is going to be used.

      integer choiceD
c  this is 1 for 2D model and 2 for 3D model.

      real tmp_block_size
c  this is a temporary value which makes sure that the kernel
c  matrix has been calculated using the same block_size.

      integer k, l, m, j
c  variables used in do loops

c==================PROGRAM============================

      do l = 0, lmax
	 do m = 0, lmax
	    do k = 0, kmax
	       A(l,m,k) = 0.0
	       B(l,m,k) = 0.0  ! initial cleaning up
	    end do  ! end k loop
	 end do  ! end m loop
      end do  ! end l loop
      do k = 1, nj
	 data(k) = 0.0
	 data_GMT(k) = 0.0  ! initial cleaning up
      end do  ! end k loop
      tmp_block_size = 0.0

c  First construct the kernel.
      inquire(file = 'kernel.dat', exist = there)
      if (.not.there) then
         call kernel(block_size, lmax, nj, ni, Gkernel, Xlm)
      elseif (there) then  ! need to read kernel
	 file_kernel = 'kernel.dat'

         open(110, file = file_kernel)
         rewind(110)
         read(110,*), tmp_block_size

	 if (tmp_block_size.eq.block_size) then

            do j = 1, nj
	       do i = 1, ni
	          Gkernel(j,i) = 0.0  ! initial cleaning
	       end do  ! end i loop
	       read(110,*), (Gkernel(j,i), i = 1, ni)  ! reading elements
            end do  ! end j loop

            close(110)

	 elseif (tmp_block_size.ne.block_size) then
	   
	    close(110)

	    ! make the kernel
	    call kernel(block_size, lmax, nj, ni, Gkernel, Xlm)

	 end if  ! if (tmp_block_size.eq.block_size)...

      end if  ! if (.not.there)

c  Read model
 35   print*, 'Choose one of the following:'
      print*, '	1.  expand a 2D (topography) model'
      print*, '	2.  expand a 3D (S, P, or R) model'
      read*, choiceD

      if ((choiceD.ne.1).and.(choiceD.ne.2)) then
         print*, choiceD, ' is not allowed.'
	 go to 35
      end if  ! if ((choiceD.ne.1).and.(choiceD.ne.2))

      call readModel(A, B, lmax, kmax, choiceD)

c  Convert a 3D model to 2D.
      if (choiceD.eq.2) then  ! if 3D model

 25      print*, 'Enter depth in km.'
         read*, depth

	 if ((depth.lt.24.4).or.(depth.gt.2891.0)) then  ! above MOHO
	    print*, 'Enter depth within the mantle (between 24.4 and 2891.0).'
	    go to 25
         end if  ! if ((depth.lt.24.4).or.(depth.gt.2891.0))

         radius = 6371.0 - depth

         do k = 0, kmax
	    Tr(k) = 0.0  ! cleaning up
         end do  ! end k loop
     
         call chebyshev_mantle (Tr, kmax, radius)
	 ! calculating radial functions at a given depth

         do i = 1, ni  ! cleaning up
	    model(i) = 0.0
         end do  ! end i loop

         i = 1  ! initial value
         do l = 0, lmax
	    do m = 0, l

	       do k = 0, kmax  ! summing over radial terms
	          model(i) = model(i) + A(l,m,k)*Tr(k)
	          if (m.ne.0) then
		     model(i+1) = model(i+1) + B(l,m,k)*Tr(k)
	          end if  ! if (m.ne.0)
	       end do  ! end k loop

	       i = i + 1
	       if (m.ne.0) then
	          i = i + 1
	       end if  ! if (m.ne.0)

	    end do  ! end m loop
         end do  ! end l loop

      elseif (choiceD.eq.1) then  ! if 2D model

         do i = 1, ni  ! cleaning up
	    model(i) = 0.0
         end do  ! end i loop

         i = 1  ! initial value
         do l = 0, lmax
	    do m = 0, l

	       model(i) = A(l,m,0)
	       if (m.ne.0) then
		  model(i+1) = B(l,m,0)
	       end if  ! if (m.ne.0)

	       i = i + 1
	       if (m.ne.0) then
	          i = i + 1
	       end if  ! if (m.ne.0)

	    end do  ! end m loop
         end do  ! end l loop

      end if  ! if (choiceD.eq.2)...

      if (i-1.ne.ni) then  ! checking length
	 print*, '!!Error(plot.f): i-1 = ', i-1, ' but ni = ', ni
         print*, 'Quitting...'
	 go to 15
      end if  ! if (i-1.ne.ni)

      do k = 1, nj
	 data(k) = 0.0
	 data_GMT(k) = 0.0  ! cleaning up
      end do  ! end k loop

c  now doing foward calculations.
      do j = 1, nj
	 do i = 1, ni
	    data(j) = data(j) + Gkernel(j,i)*model(i)
	 end do  ! end i loop
      end do  ! end j loop

c  saving the data in GMT format
      print*, 'Enter the filename to which GMT format result is saved.'
      read*, filename

	! re-arranging the data vector for GMT format.
      call switchLatLon(block_size, nj, data, data_GMT)

      open(120, file = filename)
      rewind(120)

      j = 1  ! initial position

      do lon = block_size/2.0, 360.0, block_size
         do lat = block_size/2.0 - 90.0, 90.0, block_size

	    if (choiceD.eq.2) then  ! if 3D model
	       write(120,*), lon, lat, data_GMT(j)*100.0
		   ! factor of 100.0 for converting into %
	    elseif (choiceD.eq.1) then  ! if 2D model
	       write(120,*), lon, lat, data_GMT(j)*6371.0
		   ! factor of 6371.0 for converting into km
	    end if  ! if (choiceD.eq.2)...

	    j = j + 1

	 end do  ! end lat/lon loop
      end do  ! end lon/lat loop

      close(120)

c  checking for dimension
      if (j-1.ne.nj) then
	 print*, '!!Error(plot.f-save): j-1 = ', j-1, ' but nj = ', nj
	 print*, 'Quitting plot.f ...'
	 go to 15
      end if  ! if (j-1.ne.nj)

      if (choiceD.eq.2) then  ! if 3D model
         print*, 'Try another depth?'
         read*, ans_depth
         if (ans_depth(1:1).eq.'y') then
	    go to 25
         end if  ! if (ans_depth(1:1).eq.'y')
      end if  ! if (choiceD.eq.2)

      print*, 'Try another model?'
      read*, ans_changeModel
      if (ans_changeModel(1:1).eq.'y') then
	 go to 35
      end if  ! if (ans_changeModel(1:1).eq.'y')

 15   END