SUBROUTINE WOFZ (XI, YI, U, V, FLAG)
IMPLICIT NONE
LOGICAL A, B, FLAG
REAL(8)::FACTOR,RMAXREAL,RMAXEXP,RMAXGONI
integer :: N,J,I,kapn,nu,NP1
real(8) :: XI,yi,u,v,xabs,yabs,x,y,QRHO,XABSQ,XQUAD,YQUAD,xsum,ysum,XAUX
real(8) :: u1,v1,daux,u2,v2,h,h2,qlambda,rx,ry,sx,sy,tx,ty,C,w1
PARAMETER (FACTOR = 1.12837916709551257388D0,&
RMAXREAL = 1.340780792994259D+154,&
RMAXEXP = 709.0895657128241D0,&
RMAXGONI = 0.6746518850690209D10)
FLAG = .FALSE.
XABS = DABS(XI)
YABS = DABS(YI)
X = XABS/6.3
Y = YABS/4.4
IF ((XABS.GT.RMAXREAL).OR.(YABS.GT.RMAXREAL)) GOTO 100
QRHO = X**2 + Y**2
XABSQ = XABS**2
XQUAD = XABSQ - YABS**2
YQUAD = 2*XABS*YABS
A = QRHO.LT.0.085264D0
IF (A) THEN
QRHO = (1-0.85*Y)*DSQRT(QRHO)
N = IDNINT(6 + 72*QRHO)
J = 2*N+1
XSUM = 1.0/J
YSUM = 0.0D0
DO I=N, 1, -1
J = J - 2
XAUX = (XSUM*XQUAD - YSUM*YQUAD)/I
YSUM = (XSUM*YQUAD + YSUM*XQUAD)/I
XSUM = XAUX + 1.0/J
enddo
U1 = -FACTOR*(XSUM*YABS + YSUM*XABS) + 1.0
V1 = FACTOR*(XSUM*XABS - YSUM*YABS)
DAUX = DEXP(-XQUAD)
U2 = DAUX*DCOS(YQUAD)
V2 = -DAUX*DSIN(YQUAD)
U = U1*U2 - V1*V2
V = U1*V2 + V1*U2
ELSE
IF (QRHO.GT.1.0) THEN
H = 0.0D0
KAPN = 0
QRHO = DSQRT(QRHO)
NU = IDINT(3 + (1442/(26*QRHO+77)))
ELSE
QRHO = (1-Y)*DSQRT(1-QRHO)
H = 1.88*QRHO
H2 = 2*H
KAPN = IDNINT(7 + 34*QRHO)
NU = IDNINT(16 + 26*QRHO)
ENDIF
B = (H.GT.0.0)
IF (B) QLAMBDA = H2**KAPN
RX = 0.0
RY = 0.0
SX = 0.0
SY = 0.0
DO N=NU, 0, -1
NP1 = N + 1
TX = YABS + H + NP1*RX
TY = XABS - NP1*RY
C = 0.5/(TX**2 + TY**2)
RX = C*TX
RY = C*TY
IF ((B).AND.(N.LE.KAPN)) THEN
TX = QLAMBDA + SX
SX = RX*TX - RY*SY
SY = RY*TX + RX*SY
QLAMBDA = QLAMBDA/H2
ENDIF
enddo
IF (H.EQ.0.0) THEN
U = FACTOR*RX
V = FACTOR*RY
ELSE
U = FACTOR*SX
V = FACTOR*SY
END IF
IF (YABS.EQ.0.0) U = DEXP(-XABS**2)
END IF
IF (YI.LT.0.0) THEN
IF (A) THEN
U2 = 2*U2
V2 = 2*V2
ELSE
XQUAD = -XQUAD
IF ((YQUAD.GT.RMAXGONI).OR.(XQUAD.GT.RMAXEXP)) GOTO 100
W1 = 2*DEXP(XQUAD)
U2 = W1*DCOS(YQUAD)
V2 = -W1*DSIN(YQUAD)
END IF
U = U2 - U
V = V2 - V
IF (XI.GT.0.0) V = -V
ELSE
IF (XI.LT.0.0) V = -V
END IF
RETURN
100 FLAG = .TRUE.
RETURN
END SUBROUTINE WOFZ