subroutine ctorh(x, y, r, theta)
! This routine converts from Cartesion coordinates
! to radial coordinates, where theta is in
! degrees measured clockwise from
! the +y-axis (standard meteorological heading).

! Modifications:
! 4/23/21 (J. Dostalek) Started this Fortran 90 version
!-------------------------------------------------------------------------------

  implicit none

  real, intent(in) :: x
  real, intent(in) :: y
  real, intent(out) :: r
  real, intent(out) :: theta

  real, parameter :: rtd = 57.296

  r = sqrt(x*x + y*y)

  if (r <= 0.0) then
    theta = 0.0
    return
  end if

  theta = rtd*acos(x/r)
  if (y < 0.0) theta = 360.0 - theta

! Convert theta to heading
  theta = 90.0 - theta
  if (theta < 0.0) theta = theta + 360.0

  return
end subroutine ctorh

subroutine rhtoc(r, thetah, x, y)
! This routine converts from radial coordinates
! to Cartesian coordinates, where theta is in
! degrees measured clockwise from
! the +y-axis (standard meteorological heading).

! Modifications:
! 4/23/21 (J. Dostalek) Started this Fortran 90 version
!-------------------------------------------------------------------------------

  implicit none

  real, intent(in) :: r
  real, intent(in) :: thetah
  real, intent(out) :: x
  real, intent(out) :: y

  real, parameter :: rtd = 57.296

  real :: theta

! Convert theta from heading to standard definition
  theta = 90.0 - thetah
  if (theta < 0.0) theta = theta + 360.0

  x = r*cos(theta/rtd)
  y = r*sin(theta/rtd)

  return
end subroutine rhtoc