anchor.diff:  !  Assumes AK135 unless otherwise stated.
   !  Input depth in km.
         implicit none
  -   
  +
         real(rs),intent(in) :: depth
         real(rs)            :: pressure
         character(len=*),optional :: model
         character(len=5)          :: model_in
         real(rs)                  :: rho, r, P, r0, r1, dep
         real(rs),parameter        :: dr=10.   ! integration width
  -            
  +
         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)') "pressure: model '",trim(model_in),"' not recognised."
  -		 write(0,'(a)') '   Available models are: ', available_models
  -		 stop
  -	  endif
  +     if (trim(model_in) /= 'PREM' .and. trim(model_in) /= 'AK135') then
  +       write(0,'(a,a,a)') "pressure: 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)') 'pressure: Depth must be in range 0--6371 km'
  -		 stop
  -	  endif
  -	  
  +     if (depth < 0. .or. depth > 6371.) then
  +       write(0,'(a)') 'pressure: Depth must be in range 0--6371 km'
  +       stop
  +     endif
  +
   !  Get radius required
  -	  r = 6371. - depth
  +     r = 6371. - depth
   
   !  Calculate pressure by integrating downwards
         P = 0.
  
positive.code:          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 ak135
  !-------------------------------------------------------------------------------
  
negative.code:       enddo
  
        pressure = P / 1.e9   ! GPa
  
     end function pressure
  !-------------------------------------------------------------------------------
  
  