anchor.diff:  
   real(8) :: P,P0,P1,P2,depth,depth0,depth1,depth2,dPdz
   character(10) :: model
  -real(8),parameter :: tol=5.  !  Accuracy of the solution / GPa
  +real(8) :: tol=5.  !  Accuracy of the solution / GPa
   real(8),parameter :: dz=100.
   
   if (command_argument_count() /= 1 .and. command_argument_count() /= 2) then
  
positive.code: endif
  
  call get_command_argument(1,model); read(model,*) P
  model = 'AK135'
  if (command_argument_count() == 2) call get_command_argument(2,model)
  
  !  Search for the depth which gives the correct pressure using Newton-Raphson
  !  f(x)=0: x_n+1 = x_n - f(x_n) / f'(x_n)
  !  f'(x) approx. by (f(x_n)-f(x_n+dx))/dx
  
negative.code: !enddo
  
  !  Actually, use the method of halves to get the correct depth:
  !  This is fine for PREM which has an analytical answer, but layered models
  !  will required some thought.
  if (model /= 'PREM') write(0,'(a)',advance='no') 'GPa2km is thinking'
  !  starting guesses
  depth0 = 0.      ;  P0 = 0.    ! Top
  depth2 = 6371.   ;  P2 = 365.  ! Bottom
  