      subroutine GEOD(r,z,fi,h)
c Program to transform Cartesian to geodetic coordinates
c based on the exact solution (Borkowski,1989)
c Input :  r, z = equatorial [m] and polar [m] components
c Output: fi, h = geodetic coord's (latitude [rad], height [m])
      implicit real*8(a-h,o-z)
c IAU ellipsoid: semimajor axis (a) and inverse flattening (fr)
      data a,fr /6378137.d0,298.257222101d0/
      b = dsign(a - a/fr,z)
      E = ((z + b)*b/a - a)/r
      F = ((z - b)*b/a + a)/r
c Find solution to: t**4 + 2*E*t**3 + 2*F*t - 1 = 0
      P = (E*F + 1.)*4d0/3.d0
      Q = (E*E - F*F)*2.d0
      D = P*P*P + Q*Q
            if(D.ge.0d0) then
       s = dsqrt(D) + Q
       s = dsign(dexp(dlog(dabs(s))/3d0),s)
       v = P/s - s
c Improve the accuracy of numeric values of v
       v = -(Q + Q + v*v*v)/(3*P)
            else
       v = 2.*dsqrt(-P)*dcos(dacos(Q/P/dsqrt(-P))/3.)
            endif
      G = .5*(E + dsqrt(E*E + v))
      t = dsqrt(G*G + (F - v*G)/(G + G - E)) - G
      fi = datan((1.d0 - t*t)*a/(2*b*t))
      h = (r - a*t)*dcos(fi) + (z - b)*dsin(fi)
      return
      end
