c--gsmain.f
c-----------------------------------------------------------
c--this line is optional
      program GSMAIN
c--use data to set some variables
      data yp,zp,xq,yq/4*0.0/,a/0.012/,rho/5000/,zq/0.2/
c--set up a do loop for various values of xp
c--usually we have integers for the do loop variable
      do i1=1,11
c--every 100 m starting at 0
	  xp=0.1*(i1-1)
c--this is where we call the subroutine
        call sphere(xq,yq,zq,a,rho,xp,yp,zp,gx,gy,gz)
c--write out the values, using free format
        write (*,*) xp,gz
      end do
c--stop execution
      stop
      end
c--gspherec.f
c---------------------------------------------------------------------------
c--FORTRAN subroutines have a statement and a name and arguments (variables)
c--subroutines cannot run alone
c--this subroutine can be called by
c-- call sphere(....)
      subroutine sphere(xq,yq,zq,a,rho,xp,yp,zp,gx,gy,gz)
c
c  SPHERE calculates the three components of gravitational attraction
c  at a single point due to a uniform sphere. 
c  
c  Input parameters:
c    Observation point is (xp,yp,zp), and center of sphere is at 
c    (xq,yq,zq).  Radius of sphere is a and density is rho.  Density
c    in units of kg/(m**3). All distance parameters in units of km.
c
c  Ouput parameters:
c    Gravitational components (gx,gy,gz) in units of mGal.
c
c--FORTRAN assumes types, i.e. all variables are real except integers (i,j,k,l,m,n)
c--so here we override the default and declare km2m as real
c--statements begin in column 7
      real km2m
c--now we have a data statement to inialize variables within /.../
c--1.e5 is 10^5	
      data gamma/6.67e-11/,si2mg/1.e5/,pi/3.14159265/,km2m/1.e3/
c--here we set some values
      ierror=0
      rx=xp-xq
      ry=yp-yq
      rz=zp-zq
c--compute distance from P to center of sphere
      r=sqrt(rx**2+ry**2+rz**2)
c--test to see if r=0, use pause to exit subroutine
      if(r.eq.0.)pause 'SPHERE:  Bad argument detected.'
c--need r^3
      r3=r**3
c--compute total mass of sphere
      tmass=4.*pi*rho*(a**3)/3.
c--now compute components of g
      gx=-gamma*tmass*rx/r3
      gy=-gamma*tmass*ry/r3
      gz=-gamma*tmass*rz/r3
      gx=gx*si2mg*km2m
      gy=gy*si2mg*km2m
      gz=gz*si2mg*km2m
c--this tells subroutine to exit
      return
c--required to have an end statement
      end