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