program SHtoGravPot
      use SHTOOLS
      use PlanetsConstants

      implicit none
      integer, parameter ::   degmax = 3600
      real*8 ::               lat, lon, cilm(2,degmax+1, degmax+1),interval
      real*8 ::               r_grav,omega,gm,param(2),a,f
      real*8 ::               rad(3601,7201),theta(3601,7201),phi(3601,7201),total(3601,7201)
      integer ::              nmax, lmax, l, i,nlat,nlong,lmaxfile
      character*80 ::         shfile
      character*1 ::          GRAVPOT
      
      write(*,*) 'Input file name of Spherical Harmonics'
      read(*,*) shfile
      call SHRead(shfile, cilm, lmaxfile, header=param(1:2))
      gm = param(1)
      r_grav = param(2)
      lmax=lmaxfile
      write(*,*) 'Input Semi-Major Axis of the Flattened Ellipsoid [m]'
      read(*,*) a
      write(*,*) 'Imput Flattening of Refference Ellipsoid'
      read(*,*) f
      omega=Omega_Moon
      write(*,*) ' Rotation Rate of the Moon [rad/sec]', omega
      write(*,*) 'Input Grid Spacing in Degree'
      read(*,*) interval
      write(*,*) 'Make Potential Grid (U) or Gravity Grid (G)'
      read(*,*) GRAVPOT
      call MakeGravGrid2D (rad,cilm,lmax,r_grav,a,f,gm,GRAVPOT,interval,nlat,nlong, &
    &      theta=theta,phi=phi,total=total,omega=omega) 
      open (10,file='radgrd.out',status='unknown')
      open (20,file='thetagrd.out',status='unknown')
      open (30,file='phigrd.out',status='unknown')
      open (40,file='totalgrd.out',status='unknown')
      do l=1,nlat
      do i=1,nlong
      lat=90.0d0-interval*dble(l-1)
      lon=0.0d0+interval*dble(i-1)
      if (GRAVPOT.eq."G") then
      write (10,*) lon, lat, rad(l,i)*(-100000.d0)
      write (20,*) lon, lat, theta(l,i)*100000.d0
      write (30,*) lon, lat, phi(l,i)*100000.d0
      write (40,*) lon, lat, total(l,i)*100000.d0
      else 
      write (10,*) lon, lat, rad(l,i)
      write (20,*) lon, lat, theta(l,i)
      write (30,*) lon, lat, phi(l,i)
      write (40,*) lon, lat, total(l,i)
      endif
      enddo
      enddo
      stop
end program  SHtoGravPot
