anchor.diff:           endif
         endif
      enddo
  -   
  +
      if (twodepths) then
   !  Extrapolate linearly to find density
         if (present(rho)) rho = den1 + (depth-dep1)*(den2-den1)/(dep2-dep1)
         if (present(vp))   vp = vp1  + (depth-dep1)*(vp2-vp1) / (dep2-dep1)
         if (present(vs))   vs = vs1  + (depth-dep1)*(vs2-vs1) / (dep2-dep1)
      endif
  -   
  +
      return
      end subroutine prem
   !-------------------------------------------------------------------------------
  
positive.code: !  Input depth in km, converted to m for calculation within the function.
  !  g = GM/r**2
  !  M = 4π∫ rho(r) r**2 dr
  	  implicit none
  	  
  	  real(rs),intent(in) :: depth
  	  character(len=*),optional :: model
  	  character(len=5) :: model_in
  	  real(rs) :: gravity,r,r0,r1,M,rho,z
  	  integer :: ir
  	  real(rs),parameter :: dr=1.
  	  	  
  !  Check if we've specified PREM instead, and stop for bad input
  	  model_in = 'AK135'
  	  if (present(model)) model_in = trim(model)
  	  if (trim(model_in) /= 'PREM' .and. trim(model_in) /= 'AK135') then
  		 write(0,'(a,a,a)') "gravity: model '",trim(model_in),"' not recognised."
  		 write(0,'(a)') '   Available models are: ', available_models
  		 stop
  	  endif
  !  Check our depth isn't silly
  	  if (depth < 0. .or. depth > 6371.) then
  		 write(0,'(a)') 'gravity: Depth must be in range 0--6371 km'
  		 stop
  	  endif
  	  
  !  Get radius required
  	  r = 6371. - depth
  	  
  	  M = 0.
  	  r0 = 0.
  !  Calculate mass up to radius required
  	  do ir=1,int(r),int(dr)
  		 r1 = real(ir)                ! Radius to top of shell
  		 z = 6371. - (r0 + r1)/2.     ! Depth to middle of shell
  		 if (trim(model_in) == 'AK135') call ak135(z,rho=rho)
  		 if (trim(model_in) == 'PREM' ) call prem(z,rho=rho)
  		 rho = rho * 1.e3             ! Convert to kg/m^3
  		 M = M + rho*(4.)*pi*((r1 * 1.e3)**2 - (r0 - 1.e3)**2)*dr*1.e3   ! Add mass within shell
  		 r0 = r1                      ! Next layer upwards
  	  enddo
  	  
  !  Calculate g
  	  gravity = bigG *  M / (r * 1.e3)**2
  	  
     end function gravity
  !-------------------------------------------------------------------------------
  
  
negative.code:       real(rs), parameter, private :: to_km = 111.194926644559 ;      
  
        real(rs), parameter, private :: big_number = 10.e36 ;     
  
        real(rs), parameter, private :: bigG = 6.67428e-11  ! Newton's gravitational constant
  
  !  ** Available models:
        character(len=32),parameter,private :: available_models(2) = (/ 'AK135', &
                                                                        'PREM ' /)
  