Comdiv Subroutine

  subroutine Comdiv(ar, ai, br, bi, cr, ci, rc )
    !+ Performs complex division `c` = `a` / `b`
    real(double) ::  ar
      !+ Real part of numerator
    real(double) ::  ai
      !+ Imaginary part of numerator
    real(double) ::  br
      !+ Real part of denominator
    real(double) ::  bi
      !+ Imaginary part of denominator
    real(double) ::  cr
      !+ Real part of quotient
    real(double) ::  ci
      !+ Imaginary part of quotient
    integer      :: rc
      !+ return code

    real(double) :: tmp
    if (br == ZERO.AND.bi == ZERO) then
       rc = 1
       return
    endif
    if (dabs(br) > dabs(bi)) then
       tmp = bi / br
       br  = tmp * bi + br
       cr  = (ar + tmp * ai) / br
       ci  = (ai - tmp * ar) / br
    else
       tmp = br / bi
       bi  = tmp * br + bi
       cr  = (tmp * ar + ai) / bi
       ci  = (tmp * ai - ar) / bi
    endif
    rc = 0
  end subroutine Comdiv !Comdiv