anchor.diff:  endif
   
   call get_command_argument(1,model); read(model,*) P
  -model = 'AK135'
  +model = 'PREM'
   if (command_argument_count() == 2) call get_command_argument(2,model)
  -
  +if (model == 'PREM') tol = 1.
   !  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
  
positive.code:    endif
     P0 = pressure(depth0,model=trim(model))
     P1 = pressure(depth1,model=trim(model))
     write(0,'(a)',advance='no') '.'
  enddo
  write(0,'(a)') ''
  write(*,'(i4)') int(depth1)
  
  end program GPa2km
  !-------------------------------------------------------------------------------
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
  