COMPLEX (8) function cgamma(mo,z)
!SUBROUTINE cgamma(mo, z, w)
!-----------------------------------------------------------------------

!        EVALUATION OF THE COMPLEX GAMMA AND LOGGAMMA FUNCTIONS

!                        ---------------

!     MO IS AN INTEGER, Z A COMPLEX ARGUMENT, AND W A COMPLEX VARIABLE.

!                 W = GAMMA(Z)       IF MO = 0
!                 W = LN(GAMMA(Z))   OTHERWISE

!-----------------------------------------------------------------------
!     WRITTEN BY ALFRED H. MORRIS, JR.
!        NAVAL SURFACE WARFARE CENTER
!        DAHLGREN, VIRGINIA

!     This version, in a subset of Fortran 90, prepared by
!     Alan.Miller @ vic.cmis.csiro.au
!     http://www.ozemail.com.au/~milleraj

!     This version is accurate to within 5 in the 14th significant
!     decimal digit.
!-----------------------------------------------------------------------

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

INTEGER, INTENT(IN)       :: mo
COMPLEX (dp), INTENT(IN)  :: z
COMPLEX (dp) :: w !, INTENT(OUT) :: w

! Local variables
COMPLEX (dp) :: eta, eta2, sum
REAL (dp), PARAMETER :: c0(12) = (/ .833333333333333E-01_dp,  &
          -.277777777777778E-02_dp, .793650793650794E-03_dp,  &
          -.595238095238095E-03_dp, .841750841750842E-03_dp,  &
          -.191752691752692E-02_dp, .641025641025641E-02_dp,  &
          -.295506535947712E-01_dp, .179644372368831_dp,      &
          -1.39243221690590_dp,     13.4028640441684_dp,      &
          -156.848284626002_dp /), pi = 3.14159265358979_dp,  &
    pi2  = 6.28318530717959_dp, alpi = 1.14472988584940_dp,  &
    hl2p = .918938533204673_dp, half = 0.5_dp
REAL (dp)  :: a, a1, a2, c, cn, cut, d, eps, et, e2t, h1, h2, s, sn, &
              s1, s2, t, t1, t2, u, u1, u2, v1, v2, w1, w2, x, y, y2
INTEGER    :: j, k, l, m, max, n, nm1
!---------------------------
!     ALPI = LOG(PI)
!     HL2P = 0.5 * LOG(2*PI)
!---------------------------

!     ****** MAX AND EPS ARE MACHINE DEPENDENT CONSTANTS.
!            MAX IS THE LARGEST POSITIVE INTEGER THAT MAY
!            BE USED, AND EPS IS THE SMALLEST REAL NUMBER
!            SUCH THAT 1.0 + EPS > 1.0.

!                      MAX = IPMPAR(3)
max = HUGE(3)
eps = EPSILON(1.0_dp)

!---------------------------
x = REAL(z, KIND=dp)
y = AIMAG(z)
IF (x < 0.0_dp) THEN
!-----------------------------------------------------------------------
!            CASE WHEN THE REAL PART OF Z IS NEGATIVE
!-----------------------------------------------------------------------
  y = ABS(y)
  t = -pi * y
  et = EXP(t)
  e2t = et * et

!     SET  A1 = (1 + E2T)/2  AND  A2 = (1 - E2T)/2

  a1 = half * (1.0_dp + e2t)
  t2 = t + t
  IF (t2 >= -0.15_dp) THEN
    a2 = -half * rexp(t2)
  ELSE
    a2 = half * (half + (half - e2t))
  END IF

!     COMPUTE SIN(PI*X) AND COS(PI*X)

  IF (ABS(x) >= MIN(REAL(MAX), 1.0_dp/eps)) GO TO 70
  k = ABS(x)
  u = x + k
  k = MOD(k,2)
  IF (u <= -half) THEN
    u = half + (half + u)
    k = k + 1
  END IF
  u = pi * u
  sn = SIN(u)
  cn = COS(u)
  IF (k == 1) THEN
    sn = -sn
    cn = -cn
  END IF

!     SET  H1 + H2*I  TO  PI/SIN(PI*Z)  OR  LOG(PI/SIN(PI*Z))

  a1 = sn * a1
  a2 = cn * a2
  a = a1 * a1 + a2 * a2
  IF (a == 0.0_dp) GO TO 70
  IF (mo == 0) THEN

    h1 = a1 / a
    h2 = -a2 / a
    c = pi * et
    h1 = c * h1
    h2 = c * h2
  ELSE

    h1 = (alpi+t) - half * LOG(a)
    h2 = -ATAN2(a2,a1)
  END IF
  IF (AIMAG(z) >= 0.0_dp) THEN
    x = 1.0_dp - x
    y = -y
  ELSE
    h2 = -h2
    x = 1.0_dp - x
  END IF
END IF
!-----------------------------------------------------------------------
!           CASE WHEN THE REAL PART OF Z IS NONNEGATIVE
!-----------------------------------------------------------------------
w1 = 0.0_dp
w2 = 0.0_dp
n = 0
t = x
y2 = y * y
a = t * t + y2
cut = 36.0_dp
IF (eps > 1.e-8_dp) cut = 16.0_dp
IF (a < cut) THEN
  IF (a == 0.0_dp) GO TO 70
 10    n = n + 1
  t = t + 1.0_dp
  a = t * t + y2
  IF (a < cut) GO TO 10

!     LET S1 + S2*I BE THE PRODUCT OF THE TERMS (Z+J)/(Z+N)

  u1 = (x*t+y2) / a
  u2 = y / a
  s1 = u1
  s2 = n * u2
  IF (n >= 2) THEN
    u = t / a
    nm1 = n - 1
    DO j = 1, nm1
      v1 = u1 + j * u
      v2 = (n-j) * u2
      c = s1 * v1 - s2 * v2
      d = s1 * v2 + s2 * v1
      s1 = c
      s2 = d
    END DO
  END IF

!     SET  W1 + W2*I = LOG(S1 + S2*I)  WHEN MO IS NONZERO

  s = s1 * s1 + s2 * s2
  IF (mo /= 0) THEN
    w1 = half * LOG(s)
    w2 = ATAN2(s2,s1)
  END IF
END IF

!     SET  V1 + V2*I = (Z - 0.5) * LOG(Z + N) - Z

t1 = half * LOG(a) - 1.0_dp
t2 = ATAN2(y,t)
u = x - half
v1 = (u*t1-half) - y * t2
v2 = u * t2 + y * t1

!     LET A1 + A2*I BE THE ASYMPTOTIC SUM

eta = CMPLX(t/a, -y/a, KIND=dp)
eta2 = eta * eta
m = 12
IF (a >= 289.0_dp) m = 6
IF (eps > 1.e-8) m = m / 2
sum = CMPLX(c0(m), 0.0_dp, KIND=dp)
l = m
DO j = 2, m
  l = l - 1
  sum = CMPLX(c0(l), 0.0_dp, KIND=dp) + sum * eta2
END DO
sum = sum * eta
a1 = REAL(sum, KIND=dp)
a2 = AIMAG(sum)
!-----------------------------------------------------------------------
!                 GATHERING TOGETHER THE RESULTS
!-----------------------------------------------------------------------
w1 = (((a1 + hl2p) - w1) + v1) - n
w2 = (a2 - w2) + v2
IF (REAL(z, KIND=dp) < 0.0_dp) GO TO 50
IF (mo == 0) THEN

!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO = 0

  a = EXP(w1)
  w1 = a * COS(w2)
  w2 = a * SIN(w2)
  IF (n == 0) GO TO 60
  c = (s1*w1 + s2*w2) / s
  d = (s1*w2 - s2*w1) / s
  w1 = c
  w2 = d
  GO TO 60
END IF

!     CASE WHEN THE REAL PART OF Z IS NONNEGATIVE AND MO IS NONZERO.
!     THE ANGLE W2 IS REDUCED TO THE INTERVAL -PI < W2 <= PI.

 40    IF (w2 <= pi) THEN
  k = half - w2 / pi2
  w2 = w2 + pi2 * k
  GO TO 60
END IF
k = w2 / pi2 - half
w2 = w2 - pi2 * REAL(k+1)
IF (w2 <= -pi) w2 = pi
GO TO 60

!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO IS NONZERO

 50    IF (mo /= 0) THEN
  w1 = h1 - w1
  w2 = h2 - w2
  GO TO 40
END IF

!     CASE WHEN THE REAL PART OF Z IS NEGATIVE AND MO = 0

a = EXP(-w1)
t1 = a * COS(-w2)
t2 = a * SIN(-w2)
w1 = h1 * t1 - h2 * t2
w2 = h1 * t2 + h2 * t1
IF (n /= 0) THEN
  c = w1 * s1 - w2 * s2
  d = w1 * s2 + w2 * s1
  w1 = c
  w2 = d
END IF

!     TERMINATION

 60    w = CMPLX(w1, w2, KIND=dp)
cgamma = w
RETURN
!-----------------------------------------------------------------------
!             THE REQUESTED VALUE CANNOT BE COMPUTED
!-----------------------------------------------------------------------
 70    w = (0.0_dp, 0.0_dp)
cgamma = w
RETURN

CONTAINS


FUNCTION rexp(x) RESULT(fn_val)
!-----------------------------------------------------------------------
!            EVALUATION OF THE FUNCTION EXP(X) - 1
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: x
REAL (dp)             :: fn_val

! Local variables
REAL (dp), PARAMETER           :: p1 =  .914041914819518E-09_dp,  &
    p2 = .238082361044469E-01_dp, q1 = -.499999999085958_dp,      &
    q2 = .107141568980644_dp,     q3 = -.119041179760821E-01_dp,  &
    q4 = .595130811860248E-03_dp
REAL (dp) :: e
!-----------------------
IF (ABS(x) <= 0.15_dp) THEN
  fn_val = x * (((p2*x + p1)*x + 1.0_dp) /  &
           ((((q4*x + q3)*x + q2)*x + q1)*x + 1.0_dp))
  RETURN
END IF

IF (x >= 0.0_dp) THEN
  e = EXP(x)
  fn_val = e * (half + (half - 1.0_dp/e))
  RETURN
END IF
IF (x >= -37.0_dp) THEN
  fn_val = (EXP(x) - half) - half
  RETURN
END IF
fn_val = -1.0_dp
RETURN
END FUNCTION rexp

!END SUBROUTINE cgamma
end function cgamma
