Added subroutine atmospheric_refraction() to compute accurate atmospheric refraction

* Added subroutine atmospheric_refraction() to compute accurate atmospheric refraction:
  * uses function aref() that computes the refraction for a given observed position
    by integrating through the troposphere and stratosphere
    * taken from Hohenkerk & Sinclair, HMNAO technical note 63 (1985)
  * change this to computing the refraction for a given true (theoretical) position
    * using an iterative method: expensive
    * ideally, adapt the integral in aref() for a direct method
Save topocentric altitude corrected for refraction in planpos(31)

* Save topocentric altitude corrected for refraction in planpos(31):
  * + do compute refraction in h<0, since it may bring the object above the horizon