Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
29276 Discussions

error #6637: This actual argument must be the name of an external user function or the name of an intrinsic function. [FUN1]

Leigh_Wardle
Beginner
‎09-21-2009 06:57 PM
2,246 Views
Solved Jump to solution 
 Hi,

I get the following compilation error on Lines 962, 997 and 1009 of SUBROUTINE HYBFHT (copy below):

error #6637: This actual argument must be the name of an external user function or the name of an intrinsic function.   [FUN1]

Each of these lines is of the form:

          CALL BESQUX(A,B,TERMR,TERMI,NG,NW,ORDER,R,FUN1,REL1,
     1  FWORK,MAXREL)

But as far as I can see I have declared FUN1 both in the calling routine (HYBFHT) and the routine being called (BESQUX).
The declaration in both cases is:

EXTERNAL COMPLEX FUN1

Any suggestions would be most appreciated.

Regards,
Leigh

SUBROUTINE HYBFHT(B,NB,NREL,TOL,NTOL,RERR,AERR,NORD,FUN1,REL1,
     *  ZWORK,FWORK,BWORK,NPCS, IPATH,MAXDIM,MAXREL, ZANS,NFUN1,IERR)
C=======================================================================
C     BY W.L.ANDERSON, U.S.GEOLOGICAL SURVEY, DENVER, COLORADO, USA.
C=======================================================================
C
C     HYBFHT - A HYBRID FAST HANKEL TRANSFORM SUBPROGRAM FOR RAPID AND
C     ACCURATE EVALUATION OF HANKEL TRANSFORMS OF ORDERS 0 AND/OR 1.
C     SUB-METHODS WITHIN HYBFHT ARE RELATED OR LAGGED CONVOLUTIONS USING
C     ANDERSON'S (1982) DIGITAL FILTERING SUBPROGRAM (HANKL), AND/OR
C     CHAVE'S (1983) DIRECT QUADRATURE SUBPROGRAM (BESAUX), WHERE THE
C     COMPLEX KERNEL FUNCTIONS CAN BE HIGHLY OSCILLATORY AS WELL AS
C     DECREASING, INCREASING, OR OF BOUNDED VARIATION.
C
C     SUBPROGRAM HYBFHT HAS SEVERAL OPTIONS THAT THE USER MAY SELECT
C     TO CONTROL HOW THE ALGORITHM PERFORMS IN SPEED AND/OR ACCURACY
C     (SEE DEFINITION OF PARAMETER IPATH).  THE IDEA WAS TO USE TWO
C     DIFFERENT ALGORITHMS IN SUCH A WAY AS TO TAKE ADVANTAGE OF THE
C     BEST FEATURES OF EACH ALGORITHM FOR TOTAL PERFORMANCE, AND TO
C     SOLVE A WIDER CLASS OF PROBLEMS INVOLVING HANKEL TRANSFORMS.
C
C     NOTE THAT THE ORIGINAL SUB-METHOD SUBPROGRAMS HAVE BEEN RENAMED
C     AND SOMEWHAT MODIFIED TO COEXIST UNDER CONTROL OF HYBFHT IN
C     SINGLE-PRECISION COMPLEX ARITHMETIC.  (A DOUBLE-PRECISION VERSION
C     USING COMPLEX*16, CALLED HY2FHT, IS ALSO AVAILABLE.)
C
C=======================================================================
C
C  REFERENCES:
C
C  ANDERSON, W.L., 1982, FAST HANKEL TRANSFORMS USING RELATED AND LAGGED
C    CONVOLUTIONS: ACM TRANS. ON MATH. SOFTWARE, VOL.8, NO.4, P.344-368.
C
C  --------, 1984, DISCUSSION ON: "NUMERICAL INTEGRATION OF RELATED
C    HANKEL TRANSFORMS BY QUADRATURE AND CONTINUED FRACTION EXPANSION"
C    (BY A.D. CHAVE), WITH REPLY BY AUTHOR: GEOPHYSICS, VOL.49, NO.10,
C    P.1811-1813.
C
C  CHAVE, A.D., 1983, NUMERICAL INTEGRATION OF RELATED HANKEL TRANSFORMS
C    BY QUADRATURE AND CONTINUED FRACTION EXPANSION: GEOPHYSICS, VOL.48,
C    NO.12, P.1671-1686.
C
C=======================================================================
C
C  INPUT PARAMETERS:
C
C      B() - (REAL) ARRAY OF HANKEL TRANSFORM ARGUMENTS B(I), I=1,..,NB.
C            FOR ONLY A SINGLE ARGUMENT B(1), USE NB=1; NB>1 IMPLIES
C            SUBROUTINES HANKL AND/OR BESAUX TO BE INVOKED AS REQUIRED
C            USING THE SPECIFIC B-ARGUMENTS AS LISTED IN THE INPUT ARRAY
C            B(I), WHERE B(I)>B(I-1)>0.0 FOR ALL I=2,...,NB. NOTE THAT
C            NB>1 OPTION IS FURTHER CONTROLLED BY PARAMETER IPATH BELOW.
C            B() IS ASSUMED DIMENSIONED B(MAXDIM), WHERE MAXDIM IS
C            DEFINED IN HYBFHT AS AN ADJUSTABLE DIMENSION.  WHEN HANKL
C            IS USED (IPATH=+-1), ARRAY B() IS USED TO INTERPOLATE
C            A CUBIC SPLINE AFTER LAGGED CONVOLUTIONS ARE PERFORMED;
C            WHEN BESAUX IS USED (IPATH=+-2), EACH B() IS USED DIRECTLY
C            TO OBTAIN RESULTS VIA A SEPARATE CALL TO BESAUX.
C       NB - (INTEGER) NUMBER OF B'S GIVEN IN ARRAY B(I), I=1,NB>=1.
C            NOTE THAT NB IS THE EXACT NUMBER OF INCREASING ARGUMENTS,
C            AND IS NOT THE NUMBER OF LAGGED CONVOLUTIONS (MB) AS
C            DEFINED IN SUBROUTINE HANKL.  IN GENERAL, 1<=NB<=NREL<=MAXREL IS REQUIRED (SEE MAXREL BELOW).
C      TOL - (REAL) TOLERANCE FACTOR PASSED TO SUBROUTINE HANKL, AS
C            REQUIRED, DEPENDING ON PARAMETER IPATH SELECTED.  NORMALLY,
C            TOL=.1E-8 IS THE SUGGESTED TOLERANCE THAT SHOULD BE USED TO
C            YIELD THE APPROXIMATE RELATIVE ERROR (.1E-6) FROM HANKL
C            USING 32-BIT FLOATING-POINT WORDS IN SINGLE-PRECISION.
C            TOL=0.0 CAN BE SET FOR THE MAXIMUM ACCURACY POSSIBLE.
C            (SEE MORE INFORMATION ON "TOL" IN SUBROUTINE HANKL.)
C     NTOL - (INTEGER) NUMBER OF CONSECUTIVE TIMES THE TOLERANCE IS MET
C            IN SUBROUTINE HANKL, AS REQUIRED, DEPENDING ON PARAMETER
C            IPATH SELECTED.  NORMALLY, FOR SMOOTH NON-OSCILLATORY
C            KERNELS, NTOL=1 IS SUFFICIENT, AND NTOL=2 OR 3 POSSIBLY FOR
C            OSCILLATING KERNELS TO AVOID A PREMATURE CUT-OFF OF THE
C            ADAPTIVE CONVOLUTION ALGORITHM IN HANKL.
C            (SEE MORE INFORMATION ON "NTOL" IN SUBROUTINE HANKL.)
C     RERR - (REAL) REQUESTED RELATIVE ERROR PASSED TO SUBROUTINE BESAUX
C            AS REQUIRED, DEPENDING ON PARAMETER IPATH SELECTED.  RERR
C            SHOULD NOT BE CHOSEN MUCH SMALLER THAN ABOUT .1E-6 BECAUSE
C            THIS IS A 32-BIT FLOATING-POINT COMPLEX VERSION OF BESAUX
C            TO COEXIST WITH HANKL IN SINGLE-PRECISION COMPLEX.
C            (SEE MORE INFORMATION ON "RERR" IN SUBROUTINE BESAUX.)
C     AERR - (REAL) REQUESTED ABSOLUTE ERROR PASSED TO SUBROUTINE BESAUX
C            AS REQUIRED, DEPENDING ON PARAMETER IPATH SELECTED.  AERR
C            SHOULD NOT BE CHOSEN TOO CLOSE TO ZERO, BUT THIS DEPENDS
C            ON THE PARTICULAR "FUN1" AND/OR "REL1" BEING USED.  A GOOD
C            RULE OF THUMB IS TO GENERALLY CHOOSE AERR<=1.E-2*RERR.
C            (SEE MORE INFORMATION ON "AERR" IN SUBROUTINE BESAUX, AND
C            ALSO IN CHAVE, 1983, P.1673.)
C   NORD() - (INTEGER) ARRAY OF HANKEL TRANSFORM ORDERS NORD(I), I=1,
C            NREL>=1.  NOTE THAT NORD(I)=0 OR 1 IS ALWAYS REQUIRED FOR
C            I=1,NREL. ONLY ORDERS 0 OR 1 MAY BE USED IN HYBFHT, IN ANY
C            COMBINATION; HOWEVER, NOT ALL NREL TRANSFORMS CAN BE FOUND
C            AS RELATED TRANSFORMS IF ORDERS ARE UNEQUAL AND WHEN IPATH
C            IS SELECTED SUCH THAT SUBROUTINE BESAUX IS CALLED.  THIS
C            IS NOT THE CASE IF SUBROUTINE HANKL IS CALLED, SINCE THE
C            DIGITAL FILTERS CAN OBTAIN RELATED TRANSFORMS OF EITHER
C            ORDERS 0 OR 1 IN ANY NORD() SEQUENCE DESIRED.
C            NORD() IS ASSUMED DIMENSIONED NORD(MAXREL), WHERE MAXREL
C            IS DEFINED IN HYBFHT AS AN ADJUSTABLE DIMENSION.
C     FUN1 - NAME OF AN EXTERNAL DECLARED COMPLEX FUNCTION WRITTEN BY
C            THE USER TO EVALUATE FUN1(X) FOR X>0.0.  AN EXTERNAL
C            STATEMENT MUST APPEAR IN THE USER'S CALLING PROGRAM AS:
C                  EXTERNAL FUN1
C            AND THE COMPLEX FUNCTION MUST HAVE THE GENERAL FORM:
C            (REPLACE "<...>" BELOW AS REQUIRED.)
C
C                  COMPLEX FUNCTION FUN1(X)
C                  COMMON/PASS/<...INCLUDE PARAMETERS AS NEEDED...>
C                  <...>
C                  FUN1=CMPLX(<...>,<...>)
C                  RETURN
C                  END
C
C            NOTES:  IF PARAMETERS ARE NEEDED IN FUN1 (FOR EXAMPLE A,C),
C            THEN USE "COMMON/PASS/A,C" IN CALLING PROGRAM AND FUN1.
C            FUN1 MAY BE OSCILLATING, INCREASING, DECREASING, OR OF
C            BOUNDED VARIATION AS THE ARGUMENT X>0 INCREASES. (SEE MORE
C            INFORMATION ON "FUN1" IN SUBROUTINES HANKL AND BESAUX.)
C     REL1 - NAME OF AN EXTERNAL DECLARED SUBROUTINE WRITTEN BY THE USER
C            TO EVALUATE ANY (AND ALL) NREL RELATED KERNELS FJ(NREL)
C            AS A FUNCTION OF F1=FUN1(X), ASSUMED PREVIOUSLY EVALUATED.
C            AN EXTERNAL STATEMENT MUST APPEAR IN THE USER'S PROGRAM AS:
C                  EXTERNAL REL1
C            AND THE SUBROUTINE MUST HAVE THE GENERAL FORM:
C            (REPLACE "<...>" BELOW AS REQUIRED.)
C
C                  SUBROUTINE REL1(X,F1,FJ,J)
C                  COMPLEX F1,FJ(1)
C                  COMMON/PASS/<...INCLUDE PARAMETERS AS NEEDED...>
C                  IF(J.LT.0) GO TO (1,2,<...UP TO NREL...>), IABS(J)
C            1     FJ(1)=F1
C                    IF(J.LT.0) RETURN
C            2     FJ(2)=  <...EXPRESSION IN TERMS OF X,F1, ETC...>
C                    IF(J.LT.0) RETURN
C                  <...>
C                  RETURN
C                  END
C
C        DEFINITION OF REL1(X,F1,FJ,J) PARAMETERS:
C        ----------------------------------------
C         X = REAL INPUT VARIABLE (CONTROLLED BY HYBFHT).
C        F1 = COMPLEX INPUT FUNCTION FUN1(X) VALUE GIVEN BY HYBFHT.
C       FJ()= COMPLEX OUTPUT ARRAY OF RELATED KERNEL FUNCTIONS, WHICH
C       ****  THE USER MUST CODE;  NOTE THAT EACH FJ(J) MAY OR MAY NOT
C             BE A SIMPLE FUNCTION OF X, F1 AND/OR OTHER PARAMETERS AS
C             PASSED IN COMMON/PASS/ FROM THE USER'S CALLING PROGRAM.
C         J = INTEGER INPUT VARIABLE (CONTROLLED BY SUBROUTINE HYBFHT).
C             OPTION: WHEN J<0, THEN ONLY THE IABS(J) RELATED KERNEL
C             IS RETURNED IN FJ(IABS(J)); WHEN J>0, ALL RELATED KERNELS
C             FJ(I), I=1,J ARE RETURNED, WHERE FJ(1)=F1 ALWAYS REQUIRED.
C             (SEE REL1 EXAMPLE INCLUDED WITH THE TEST DRIVER OF HYBFHT)
C  ZWORK() - (COMPLEX) ARRAY OF DIMENSION ZWORK(283,MAXREL), WHERE
C            MAXREL IS DEFINED IN HYBFHT AS AN ADJUSTABLE DIMENSION.
C            ZWORK() IS ONLY USED BY SUBROUTINE HANKL, AS REQUIRED,
C            DEPENDING ON PARAMETER IPATH SELECTION.
C  FWORK() - (COMPLEX) ARRAY OF DIMENSION FWORK(MAXREL), WHERE
C            MAXREL IS DEFINED IN HYBFHT AS AN ADJUSTABLE DIMENSION.
C            FWORK() IS USED BY EITHER SUBROUTINES HANKL OR HYBFHT,
C            AS REQUIRED, DEPENDING ON PARAMETER IPATH SELECTION.
C  BWORK() - (REAL) ARRAY OF DIMENSION BWORK(MAXDIM,9), WHERE MAXDIM
C            IS DEFINED IN HYBFHT AS AN ADJUSTABLE DIMENSION.  BWORK()
C            IS ONLY USED WHENEVER SUBROUTINE HANKL IS CALLED, AS
C            REQUIRED, DEPENDING ON PARAMETER IPATH SELECTION.
C     NPCS - (INTEGER) NUMBER OF PIECES FOR PARTIAL INTEGRATION AS
C            USED BY SUBROUTINE BESAUX, AS REQUIRED, DEPENDING ON
C            PARAMETER IPATH SELECTION.  NORMALLY SET NPCS=1, BUT NPCS
C            CAN BE SET GREATER THAN 1 FOR DIFFICULT KERNELS.
C            (SEE MORE INFORMATION ON "NPCS" IN SUBROUTINE BESAUX.)
C    IPATH - (INTEGER) OPTION CONTROL PATH IN SUBROUTINE HYBFHT DEFINED
C            AS FOLLOWS:      ---------------------------------
C        = 1 TO FORCE USE OF ONLY SUBROUTINE HANKL, AND AS CONTROLLED
C            BY INPUT PARAMETERS ABOVE RELATED TO SUBROUTINE HANKL.
C        = 2 TO FORCE USE OF ONLY SUBROUTINE BESAUX, AND AS CONTROLLED
C            BY INPUT PARAMETERS ABOVE RELATED TO SUBROUTINE BESAUX.
C       = -1 TO BEGIN WITH SUBROUTINE HANKL, BUT CHECK DURING ADAPTIVE
C            CONVOLUTION IF A HIGHLY OSCILLATORY KERNEL IS BEING USED
C            SUCH THAT THE OSCILLATIONS EXCEED THE FILTER'S NYQUIST
C            FREQUENCY AS TRANSFORMED TO LOG-SPACE, I.E., 1/(2*0.2).
C            IF THIS OCCURS, THEN HYBFHT WILL AUTOMATICALLY ABORT HANKL
C            AND SWITCH TO SUBROUTINE BESAUX.  THE NUMBER OF FUNCTION
C            FUN1 EVALUATIONS WILL BE NOTED IN OUTPUT ARRAY NFUN1(1) FOR
C            HANKL (UNTIL ABORTED) AND NFUN1(2) AT EXIT.  IF NFUN1(2)=0
C            AT EXIT, THEN NO ABORT OCCURRED, AND SUBROUTINE HANKL WAS
C            ONLY USED.
C       = -2 PERFORM A RELATIVE ERROR CHECK USING BOTH SUBROUTINES HANKL
C            AND BESAUX AT THE END POINTS B(1) AND B(NB) ONLY WHEN
C            NB>2 IS SELECTED.  THE ASSUMPTION IS THAT BESAUX IS MORE
C            ACCURATE, IN GENERAL, AND HENCE CAN BE USED TO DETERMINE
C            IF HANKL IS SUFFICIENTLY ACCURATE FOR LAGGED CONVOLUTION
C            TO BE PERFORMED FOR THE ENTIRE RANGE OF B(I),I=1,NB>2.
C            IF THE RELATIVE ERRORS AT B(1) OR B(NB) ARE LESS THAN OR
C            EQUAL TO PARAMETER RERR*10, THEN SUBROUTINE HANKL IS USED;
C            ELSE SUBROUTINE BESAUX IS USED.  THE NUMBER OF FUNCTION
C            FUN1 EVALUATIONS WILL BE NOTED IN OUTPUT ARRAY NFUN1(1)
C            AND NFUN1(2), WHERE SOME EXTRA EVALUATIONS ARE SOMETIMES
C            NECESSARY TO PERFORM THE RELATIVE ERROR TESTS, AND WHEN
C            HANKL IS FOLLOWED TO COMPLETION.  HOWEVER, IF NB>>2,
C            THEN MANY FUN1 EVALUATIONS ARE SAVED WITH HANKL THAN WITH
C            BESAUX DUE TO LAGGED CONVOLUTIONS.  THUS IPATH=-2 IS SAFE
C            AND USUALLY NOT TIME-CONSUMING IN MOST CASES WHEN NB>>2.
C            (IF IPATH=-2 AND NB<=2, THEN ONLY BESAUX IS USED.)
C   MAXDIM - (INTEGER) MAXIMUM DIMENSION USED IN PROGRAM THAT CALLS
C            SUBROUTINE HYBFHT FOR ARRAYS: B(MAXDIM),BWORK(MAXDIM,9),
C            ZANS(MAXDIM,MAXREL).  NB<=MAXDIM IS REQUIRED IF IPATH=+-2
C            AND MB=AINT(5.*ALOG(B(NB)/B(1)))+5<=MAXDIM IF IPATH=+-1.
C            IT IS RECOMMENDED THE LATTER RELATION ALWAYS BE USED SUCH
C            THAT MAXDIM>>NB. (THIS WILL ALLOW FOR ANY IPATH OPTION.)
C   MAXREL - (INTEGER) MAXIMUM DIMENSION USED IN PROGRAM THAT CALLS
C            SUBROUTINE HYBFHT FOR ARRAYS: ZWORK(283,MAXREL),
C            ZANS(MAXDIM,MAXREL), AND FWORK(MAXREL).  NOTE THAT
C            NREL<=MAXREL IS REQUIRED.
C
C-----------------------------------------------------------------------
C
C  OUTPUT PARAMETERS:
C
C   ZANS() - (COMPLEX) ARRAY OF HANKEL TRANSFORM RESULTS ZANS(I,J),
C            I=1,NB AND J=1,NREL.  THE COMPLEX ARRAY IS ASSUMED TO BE
C            DIMENSIONED IN THE CALLING PROGRAM AS: ZANS(MAXDIM,MAXREL).
C            NOTE THAT 1<=NB<=MAXDIM AND 1<=NREL<=MAXREL ARE REQUIRED.
C  NFUN1() - (INTEGER) ARRAY OF NUMBER OF FUNCTION FUN1 EVALUATIONS
C            TAKEN BY EITHER SUBROUTINE HANKL (STORED IN NFUN1(1))
C            AND/OR SUBROUTINE BESAUX (STORED IN NFUN1(2)).  USUALLY,
C            NFUN1(1) IS NOT MUCH LARGER THAN THE NUMBER OF WEIGHTS
C            NEEDED FOR A SINGLE DIRECT CONVOLUTION (FOR ANY NB, NREL).
C            WITH REGARD TO IPATH SELECTED, NFUN1(2) IS GENERALLY LARGER
C            THAN NFUN1(1), BUT THIS DEPENDS ON PARAMETERS TOL, RERR,
C            AERR SELECTED.  NOTE THAT NFUN1() REFLECTS THE OUTCOME OF
C            PARAMETER IPATH, AND AS A RESULT OF THE FUN1 AND REL1
C            SUBPROGRAMS USED IN A CALL TO SUBROUTINE HYBFHT.
C     IERR - (INTEGER) ERROR RETURN CODE DEFINED AS FOLLOWS:
C        = 0 NO ERROR DETECTED IN INPUT PARAMETERS.
C        = 1 ERROR RETURN: NB<1 OR NB>MAXDIM NOT ALLOWED.
C        = 2 ERROR RETURN: NREL<1 OR NREL>MAXREL NOT ALLOWED.
C        = 3 ERROR RETURN: B(I-1)<=0 OR B(I)<=B(I-1), I=2,NB.
C            (ARRAY B(I)=0.0 OR NOT INCREASING FOR I=1,NB.)
C        = 4 ERROR IPATH=+-1: MB=AINT(5.*ALOG(B(NB)/B(1)))+5>MAXDIM.
C            (MB TOO LARGE, OR MAXDIM TOO SMALL HAS BEEN SELECTED.)
C        = 5 ERROR IPATH=+-1: SPLIN1 ERROR RETURN AFTER CALL SPLIN1.
C            (THIS SHOULD NEVER OCCUR IN PRACTICE DUE TO PRIOR CHECKS.)
C        = 6 ERROR RETURN: IPATH<-2, IPATH=0, OR IPATH>2 NOT ALLOWED.
C       = 11 ERROR IPATH=+-1 OR -2: ERROR RETURN AFTER CALL HANKL.
C            (SEE SUBROUTINE HANKL IERR=1 CAUSES AND CORRECTIVE ACTION.)
C       = 21 ERROR IPATH=2: ERROR RETURN AFTER CALL BESAUX.
C            (SEE SUBROUTINE BESAUX IERR=1 CAUSES; POSSIBLE CORRECTION
C             IS TO TRY USING NPCS>1 AND INCREASE RERR, AERR.)
C       = 31 ERROR IPATH=-1 OR -2: ERROR RETURN AFTER CALL BESAUX.
C            (SEE SUBROUTINE BESAUX IERR=1 CAUSES; POSSIBLE CORRECTION
C             IS TO TRY USING NPCS>1 AND INCREASE RERR, AERR.)
C
C=======================================================================
C
C  HYBFHT USAGE:  (REPLACE "<...>" BELOW AS REQUIRED.)
C
C    THE USER'S PROGRAM SHOULD CONTAIN THE FOLLOWING STATEMENTS:
C    ----------------------------------------------------------
C
C     PARAMETER (MAXDIM= <...>,MAXREL= <...>)
C     INTEGER NORD(MAXREL),NFUN1(2)
C     REAL B(MAXDIM),BWORK(MAXDIM,9)
C     COMPLEX ZWORK(283,MAXREL),ZANS(MAXDIM,MAXREL),FWORK(MAXREL)
C     EXTERNAL FUN1,REL1 <,...POSSIBLY INCLUDE FUN2,REL2, ETC...>
C     COMMON/PASS/<...USE TO PASS OTHER PARAMETERS TO FUN1, ETC...>
C     <...>
C     CALL HYBFHT(<...SEE CALLING SEQUENCE ABOVE...>)
C     IF(IERR.NE.0) STOP <...>
C     <...>
C     END
C     COMPLEX FUNCTION FUN1(X)
C     <...SEE FORM OF FUN1 ABOVE...>
C     END
C     SUBROUTINE REL1(X,F1,FJ,J)
C     <...SEE FORM OF REL1 ABOVE...>
C     END
C
C=======================================================================
C
      INTEGER NORD(MAXREL),NFUN1(2)
      REAL B(MAXDIM),BWORK(MAXDIM,9),DER(2)
      COMPLEX ZWORK(283,MAXREL),ZANS(MAXDIM,MAXREL),FWORK(MAXREL)
      DOUBLE PRECISION BESR,BESI
      EXTERNAL COMPLEX FUN1 !20090922
      EXTERNAL REL1
      COMMON/CONTRL/JPATH,J
      COMMON/TEST/NG,NF2,NI
      
      DATA DER/2*0.0/
C
C  PRESET AND TEST SOME PARAMETERS (OTHERS TESTED IN SUBR. CALLS)
      IERR=0
      NFUN1(1)=0
      NFUN1(2)=0
      IF(NB.LT.1.OR.NB.GT.MAXDIM) THEN
        IERR=1
        RETURN
      ELSE IF(NREL.LT.1.OR.NREL.GT.MAXREL) THEN
        IERR=2
        RETURN
      ELSE IF(NB.GT.1) THEN
        DO 5 I=2,NB
          IF(B(I).LE.B(I-1).OR.B(I-1).LE.0.0) THEN
            IERR=3
            RETURN
          ENDIF
5       CONTINUE
      ENDIF
C  PROCESS VIA IPATH SELECTION (SEE IPATH DEFINITIONS)
      JPATH=IPATH
10    KPATH=IABS(JPATH)
C--BEGINNING OF "IF(KPATH.EQ.1) THEN" LOOP:
      IF(KPATH.EQ.1) THEN
        IF(NB.GT.1) THEN
          BMAX=B(NB)
          MB=AINT(5.*ALOG(BMAX/B(1)))+5
          IF(MB.GT.MAXDIM) THEN
            IERR=4
            RETURN
          ENDIF
        ELSE
          BMAX=B(1)
          MB=1
        ENDIF
C  FHT FILTER METHOD (JPATH=-1 OR +1)
C  **********************************
        CALL HANKL(BMAX,MB,NREL,TOL,NTOL,NORD,FUN1,REL1,ZWORK,FWORK,
     1 MAXDIM,MAXREL, ZANS,BWORK(1,1),NF1,IER)
        NFUN1(1)=NFUN1(1)+NF1
        IF(IER.NE.0) THEN
          IERR=IER+10
          RETURN
        ENDIF
C  JPATH=-2 CAN BE SET AUTOMATICALLY BY HANKL IF IPATH=-1 WAS SET AND
C  POOR ACCURACY OBTAINED BY HANKL DUE TO HIGHLY OSCILLATING KERNELS.
        IF(JPATH.EQ.-2) THEN
          I1=1
          INB=NB
          GO TO 110
        ENDIF
        IF(NB.EQ.1) RETURN
C  NB>1, GET SPLINE INTERPOLATED ZANS(I,NREL) AT B(I),I=1,NB.
        DO 100 J=1,NREL
          K=8
          DO 20 I=1,MB
            BWORK(I,8)=0.0
            BWORK(I,9)=0.0
            IF(J.EQ.1) BWORK(I,1)=ALOG(BWORK(I,1))
            BWORK(I,2)=REAL(ZANS(I,J))
            IF(AIMAG(ZANS(I,J)).NE.0.0) K=-8
20        CONTINUE
30        KK=IABS(K)
          CALL SPLIN1(MB,0.0,BWORK(1,1),BWORK(1,2),BWORK(1,3),
     1     BWORK(1,4),BWORK(1,5),0,DER,BWORK(1,6),BWORK(1,7))
          IF(MB.LT.0) THEN
            IERR=5
            RETURN
          ENDIF
          DO 40 I=1,NB
            BLOG=ALOG(B(I))
            IF(BLOG.LT.BWORK(1,1)) THEN
              BWORK(I,KK)=BWORK(I,2)
            ELSE IF(BLOG.GT.BWORK(MB,1)) THEN
              BWORK(I,KK)=BWORK(MB,2)
            ELSE
              CALL SPOINT(MB,BWORK(1,1),BWORK(1,2),BWORK(1,3),
     1         BWORK(1,4),BWORK(1,5),BLOG,BWORK(I,KK))
            ENDIF
40        CONTINUE
          IF(K.EQ.-8) THEN
            K=9
            DO 50 I=1,MB
              BWORK(I,2)=AIMAG(ZANS(I,J))
50          CONTINUE
            GO TO 30
          ELSE
            DO 60 I=1,NB
              ZANS(I,J)=CMPLX(BWORK(I,8),BWORK(I,9))
60          CONTINUE
          ENDIF
100     CONTINUE
C  DONE WITH SPLINE INTERPOLATION WHEN NB>1 AND NREL>=1.
        RETURN
C--REST OF "IF(KPATH.EQ.1) THEN" LOOP ABOVE:
      ELSE IF(KPATH.EQ.2) THEN
        I1=1
        INB=NB
        IF(IPATH.EQ.2.OR.NB.LT.3) GO TO 110
C  SPECIAL PREPROCESSING ONLY WHEN IPATH=-2 AND NB>2
        RERR10=10.*RERR
        DO 105 I=NB,1,1-NB
          CALL HANKL(B(I),1,NREL,TOL,NTOL,NORD,FUN1,REL1,ZWORK,FWORK,
     1   MAXDIM,MAXREL, ZANS,BWORK(1,1),NF1,IER)
          NFUN1(1)=NFUN1(1)+NF1
          IF(IER.NE.0) THEN
            IERR=IER+10
            RETURN
          ENDIF
          DO 105 J=1,NREL
            IF(J.GT.1.AND.NORD(J).EQ.NORD(1)) THEN
              NEW=2
              CALL RELATE(FWORK,MAXREL,REL1,J,0)
            ELSE
              IF(NREL.GT.1) THEN
                NEW=1
              ELSE
                NEW=0
              ENDIF
            ENDIF
            TESTR=REAL(ZANS(1,J))
            TESTI=AIMAG(ZANS(1,J))
C  GET ASSUMED TRUE VALUES FROM BESAUX AT THIS (I,J),I=NB OR 1, J=1,NREL
C  AND COMPARE WITH HANKL TEST VALUES AT THE END POINTS B(NB), B(1).
C  SWITCH TO IPATH=1 OR 2 DEPENDING ON RELATIVE ERROR TESTS.
            CALL BESAUX(BESR,BESI,DFLOAT(NORD(J)),1,7,DBLE(B(I)),FUN1,
     1     REL1, FWORK,MAXREL,DBLE(RERR),DBLE(AERR),NPCS,NEW,IER)
            NFUN1(2)=NFUN1(2)+NF2
            IF(IER.NE.0) THEN
              IERR=IER+30
C             RETURNS WITH IERR=31
            ENDIF
            ZANS(I,J)=CMPLX(SNGL(BESR),SNGL(BESI))
            IF(J.EQ.1.AND.NREL.GT.1) CALL RELATE(FWORK,MAXREL,REL1,J,1)
            TRUE=SNGL(BESR)
            IF(TRUE.NE.0.0) THEN
              RELERR=ABS((TESTR-TRUE)/TRUE)
              IF(RELERR.GT.RERR10) THEN
                IF(I.EQ.NB) THEN
                  INB=NB-1
                ELSE IF(I.EQ.1) THEN
                  I1=2
                  INB=NB-1
                ENDIF
                JPATH=2
                IF(J.EQ.NREL) GO TO 110
              ENDIF
            ENDIF
            TRUE=SNGL(BESI)
            IF(TRUE.NE.0.0) THEN
              RELERR=ABS((TESTI-TRUE)/TRUE)
              IF(RELERR.GT.RERR10) THEN
                IF(I.EQ.NB) THEN
                  INB=NB-1
                ELSE IF(I.EQ.1) THEN
                  I1=2
                  INB=NB-1
                ENDIF
                JPATH=2
                IF(J.EQ.NREL) GO TO 110
              ENDIF
            ENDIF
105     CONTINUE
        JPATH=1
        GO TO 10
C--REST OF "IF(KPATH.EQ.1) THEN" LOOP ABOVE:
      ELSE
        IERR=6
        RETURN
      ENDIF
C  QUADRATURE METHOD (JPATH=-2 OR +2)
C  **********************************
110   DO 200 I=I1,INB
      DO 200 J=1,NREL
        IF(J.GT.1.AND.NORD(J).EQ.NORD(1)) THEN
          NEW=2
          CALL RELATE(FWORK,MAXREL,REL1,J,0)
        ELSE
          IF(NREL.GT.1) THEN
            NEW=1
          ELSE
            NEW=0
          ENDIF
        ENDIF
        CALL BESAUX(BESR,BESI,DFLOAT(NORD(J)),1,7,DBLE(B(I)),FUN1,
     1 REL1, FWORK,MAXREL,DBLE(RERR),DBLE(AERR),NPCS,NEW,IER)
        NFUN1(2)=NFUN1(2)+NF2
        IF(IER.NE.0) THEN
          IF(JPATH.EQ.-2) THEN
            IERR=IER+30
          ELSE
            IERR=IER+20
          ENDIF
C         RETURNS WITH IERR>=21
        ENDIF
        ZANS(I,J)=CMPLX(SNGL(BESR),SNGL(BESI))
        IF(J.EQ.1.AND.NREL.GT.1) CALL RELATE(FWORK,MAXREL,REL1,J,1)
200   CONTINUE
      RETURN
      END
      SUBROUTINE RELATE(FWORK,MAXREL,REL1,JREL,ISAVE)
C----------------------------------------------------------------------
C  GET RELATED KERNELS AS REQUIRED FOR SUBROUTINE BESAUX VIA CALL REL1.
C  NOTE ALL I/O PARAMETERS CONTROLLED BY HYBFHT.  THE USER ONLY NEEDS
C  TO CODE EXTERNAL SUBROUTINE REL1 TO COMPUTE ANY (OR ALL) RELATED
C  KERNELS. (SEE COMMENTS ABOUT SUBROUTINE REL1 IN HYBFHT.)
C----------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,K,O-Z)
C  NTERM IS MAXIMUM NUMBER OF BESSEL FUNCTION LOOPS STORED IF NEW.NE.0
      PARAMETER (NTERM=50)
      SAVE KSAVE,KASAVE,KBSAVE,NKSAVE,NPSAVE,NPSAV
      REAL X,FRE,FIM
      DIMENSION KARG(255,NTERM),KERN(510,NTERM),NK(NTERM),
     1 KBES(510,NTERM),KSAVE(510,NTERM),KASAVE(255,NTERM),
     2 KBSAVE(510,NTERM),NKSAVE(NTERM)
      COMPLEX FWORK(MAXREL)
      EXTERNAL REL1
      COMMON/BESINT/NK,NP,NPS,KARG,KERN,KBES
      IF(ISAVE.EQ.1) THEN
        NPSAV=NP
        NPSAVE=NPS
        DO 1 I=1,MIN(NPS,100)
          NKSAVE(I)=NK(I)
1       CONTINUE
      ELSE
        NP=NPSAV
        NPS=NPSAVE
        DO 2 I=1,MIN(NPS,100)
          NK(I)=NKSAVE(I)
2       CONTINUE
      ENDIF
C
      DO 10 I=1,MIN(NPS,100)
        LK=1
        DO 10 J=1,2**(NK(I)+1)-1
          IF(ISAVE.EQ.1) THEN
            TBES=KBES(J,I)
            IF(TBES.NE.0.0) THEN
              RBES=1.0/TBES
            ELSE
              RBES=0.0
            ENDIF
            KSAVE(LK,I)=KERN(LK,I)*RBES
            KSAVE(LK+1,I)=KERN(LK+1,I)*RBES
            KASAVE(J,I)=KARG(J,I)
            KBSAVE(J,I)=KBES(J,I)
            LK=LK+2
            GO TO 10
          ENDIF
C--GET SAVED ARGUMENT, BESSEL, AND KERNEL FUNCTION VALUE FOR F1.
          X=KASAVE(J,I)
          TBES=KBSAVE(J,I)
          FRE=KSAVE(LK,I)
          FIM=KSAVE(LK+1,I)
C--GET RELATED KERNEL FROM SAVED KERNEL AS REQUIRED VIA CALL REL1.
          CALL REL1(X,CMPLX(FRE,FIM),FWORK,-JREL)
          KERN(LK,I)=REAL(FWORK(JREL))*TBES
          KERN(LK+1,I)=AIMAG(FWORK(JREL))*TBES
          LK=LK+2
10    CONTINUE
      RETURN
      END
      SUBROUTINE SPLIN1(M,H,X,Y,A,B,C,IT,D,P,S)
C--ONE DIMENSIONAL CUBIC SPLINE COEFFICIENT DETERMINATION.
C
C        BY  W.L.ANDERSON, U.S. GEOLOGICAL SURVEY, DENVER, COLORADO
C
C  PARMS--- M= NUMBER OF DATA POINTS .GT. 2
C           H= EQUAL INTERVAL OPTION WHEN H.GT.0. (USE DUMMY X HERE),
C              UNEQUAL INTERVALS IF H=0. (X REQUIRED STORAGE)
C           X= INDEP.VAR WHEN H=0. (DIM .GE. M).
C           Y= DEPENDENT VARIABLE  (DIM .GE. M).
C           A,B,C=COEFF.ARRAYS (EACH DIM .GE. M)
C                 RESULTS ARE RETURNED IN 1ST(M-1) ELEMENTS OF A,B,&C.
C                 ALSO USED AS WORK ARRAYS DURING EXECUTION.
C           IT= TYPE OF BOUNDARY CONDITION SUPPLIED IN D ARRAY. USE
C              IT=1 IF 1ST DERIVATIVES GIVEN AT END POINTS, OR
C              IT=0 IF 2ND DERIVATIVES GIVEN AT END POINTS.
C           D= BOUNDARY ARRAY (DIM 2) AT POINT 1 AND M RESPECTIVELY.
C           P,S= WORK ARRAYS (EACH DIM=M).
C--ERROR RETURN WITH M=-(ABS(M)) IF ANY PARM OUT OF RANGE.
C  THE RESULTING CUBIC SPLINE IS OF THE FORM:
C     Y=Y(I)+A(I)*(X-X(I))+B(I)*(X-X(I))**2+C(I)*(X-X(I))**3
C       FOR I=1,2,...,M-1
C
C
      REAL*4  X(1),Y(1),A(1),B(1),C(1),D(2),P(1),S(1),MUL
      IF(IT.LT.0.OR.IT.GT.1.OR.H.LT.0..OR.M.LT.3) GO TO 999
      N=M-1
      IF(IT.EQ.0) GO TO 20
C--1ST DERIVATIVE BOUNDARIES GIVEN
      NE=N-1
      IF(H) 999,11,1
C--EQUAL SPACING H .GT. 0. AND IT=1
    1 HH=3.0/H
      DO 2 I=1,NE
      B(I)=4.0
      C(I)=1.0
      A(I)=1.0
    2 P(I)=HH*(Y(I+2)-Y(I))
      P(1)=P(1)-D(1)
      P(NE)=P(NE)-D(2)
C--SOLUTION OF TRIDIAGONAL MATRIX EQ. OF ORDER NE
    3 C(1)=C(1)/B(1)
      P(1)=P(1)/B(1)
      DO 4 I=2,NE
      MUL=1.0/(B(I)-A(I)*C(I-1))
      C(I)=MUL*C(I)
    4 P(I)=MUL*(P(I)-A(I)*P(I-1))
C--OBTAIN SPLINE COEFFICIENTS
      A(NE+IT)=P(NE)
      I=NE-1
    5 A(I+IT)=P(I)-C(I)*A(I+IT+1)
      I=I-1
       IF(I.GE.1) GO TO 5
      IF(IT.EQ.0) GO TO 6
      A(1)=D(1)
      A(M)=D(2)
    6 IF(H.EQ.0.) GO TO 14
      HH=1.0/H
      DO 7 I=1,N
      MUL=HH*(Y(I+1)-Y(I))
      B(I)=HH*(3.0*MUL-(A(I+1)+2.0*A(I)))
    7 C(I)=HH*HH*(-2.0*MUL+A(I+1)+A(I))
      RETURN
C--UNEQUAL SPACING H=0.. AND IT=1
   11 DO 12 I=1,N
   12 S(I+1)=X(I+1)-X(I)
      DO 13 I=1,NE
      B(I)=2.0*(S(I+1)+S(I+2))
      C(I)=S(I+1)
      A(I)=S(I+2)
   13 P(I)=3.0*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
     $ (S(I+1)*S(I+2))
      P(1)=P(1)-S(3)*D(1)
      P(NE)=P(NE)-S(N)*D(2)
      GO TO 3
   14 DO 15 I=1,N
      HH=1.0/S(I+1)
      MUL=(Y(I+1)-Y(I))*HH**2
      B(I)=3.0*MUL-(A(I+1)+2.0*A(I))*HH
   15 C(I)=-2.0*MUL*HH+(A(I+1)+A(I))*HH**2
      RETURN
C--2ND DERIVATIVE BOUNDARIES GIVEN
   20 NE=N+1
      IF(H) 999,31,21
C--EQUAL SPACING H .GT. 0 AND IT=0
   21 HH=3.0/H
      DO 22 I=2,N
      B(I)=4.0
      C(I)=1.0
      A(I)=1.0
   22 P(I)=HH*(Y(I+1)-Y(I-1))
      B(1)=2.0
      B(NE)=2.0
      C(1)=1.0
      C(NE)=1.0
      A(NE)=1.0
      P(1)=HH*(Y(2)-Y(1))-0.5*H*D(1)
      P(NE)=HH*(Y(M)-Y(N))+0.5*H*D(2)
      GO TO 3
C--UNEQUAL SPACING H=0 AND IT=0
   31 DO 32 I=1,N
   32 S(I+1)=X(I+1)-X(I)
      N1=N-1
      DO 33 I=1,N1
      B(I+1)=2.0*(S(I+1)+S(I+2))
      C(I+1)=S(I+1)
      A(I+1)=S(I+2)
   33 P(I+1)=3.0*(S(I+1)**2*(Y(I+2)-Y(I+1))+S(I+2)**2*(Y(I+1)-Y(I)))/
     *     (S(I+1)*S(I+2))
      B(1)=2.0
      B(NE)=2.0
      C(1)=1.0
      C(NE)=1.0
      A(NE)=1.0
      P(1)=3.0*(Y(2)-Y(1))/S(2)-0.5*S(2)*D(1)
      P(NE)=3.0*(Y(M)-Y(N))/S(M)+0.5*S(M)*D(2)
      GO TO 3
  999 M=-IABS(M)
      RETURN
      END
      SUBROUTINE SPOINT(M,X,Y,A,B,C,XX,YY)
C--GIVEN CUBIC SPLINE COEFF'S A,B,C,AND M OBS.DATA ARRAYS X,Y
C  SPOINT EVALUATES THE PIECEWISE CUBIC SPLINE ORDINATE YY AT THE
C  ABSCISSA XX, WHERE XX IS IN THE CLOSED INTERVAL (X(1),X(M)).
C  NOTE: IF COMPUTING OVER EQUAL INTERVALS, USE THE SUBR 'CUBIC'
C  WHICH REQUIRES ONLY ONE CALL.
C
      DIMENSION X(1),Y(1),A(1),B(1),C(1)
      IF(XX.LT.X(1).OR.XX.GT.X(M)) GO TO 9
      M1=M-1
      DO 1 I=1,M1
      J=I
      IF(XX.LE.X(I+1)) GO TO 2
    1 CONTINUE
    9 WRITE(6,60) XX,X(1),X(M)
   60 FORMAT('0ERROR IN SPOINT CALL--XX=',E16.8,' NOT IN CLOSED INTERVAL
     * (',E16.8,',',E16.8,')')
      RETURN
    2 Z=XX-X(J)
      YY=Y(J)+((C(J)*Z+B(J))*Z+A(J))*Z
      RETURN
      END
      SUBROUTINE BESAUX(BESR,BESI,ORDER,NL,NU,R,FUN1,REL1,FWORK,MAXREL,
     $  RERR,AERR, NPCS,NEW,IERR)
C
C  PERFORMS AUTOMATIC CALCULATION OF BESSEL TRANSFORM TO SPECIFIED
C  RELATIVE AND ABSOLUTE ERROR
C
C  ARGUMENT LIST:
C
C  BESR,BESI-REAL AND IMAGINARY PARTS RETURNED BY BESAUX
C  ORDER-ORDER OF THE BESSEL FUNCTION
C  NL-LOWER LIMIT FOR GAUSS ORDER TO START COMPUTATION
C  NU-UPPER LIMIT FOR GAUSS ORDER
C     NU,NL=1,...7 SELECTS 3,7,15,31,63,127,AND 255 POINT GAUSS
C     QUADRATURE BETWEEN THE ZERO CROSSINGS OF THE BESSEL FUNCTION
C  R-ARGUMENT OF THE BESSEL FUNCTION
C///////////////////////////////////////////////////////////////////////
C     THE FOLLOWING PARAMETERS ARE INCLUDED SO THAT BESAUX CAN BE RUN
C     UNDER THE GENERAL CONTROL OF SUBROUTINE HYBFHT:
C///////////////////////////////////////////////////////////////////////
C     FUN1 - NAME OF AN EXTERNAL DECLARED COMPLEX FUNCTION WRITTEN BY
C            THE USER TO EVALUATE FUN1(G) FOR G>0.0.  AN EXTERNAL
C            STATEMENT MUST APPEAR IN THE USER'S CALLING PROGRAM AS:
C                  EXTERNAL FUN1
C            AND THE COMPLEX FUNCTION MUST HAVE THE GENERAL FORM:
C            (REPLACE "<...>" BELOW AS REQUIRED.)
C
C                  COMPLEX FUNCTION FUN1(G)
C                  COMMON/PASS/<...INCLUDE PARAMETERS AS NEEDED...>
C                  <...>
C                  FUN1=CMPLX(<...>,<...>)
C                  RETURN
C                  END
C
C        NOTES: SEE COMMENTS IN SUBROUTINES HYBFHT AND HANKL ON FUN1.
C        THE EXTERNAL FUNCT CALL WAS CHANGED FROM THE FORM:
C        CALL FUN1(X,YR,YI)  --  ORIGINAL METHOD IN BESAUT -- TO
C        Z=FUN1(X),  WHERE COMPLEX FUN1 USED IN BESAUX, BESTRX, AND
C        BESQUX.  THIS METHOD WAS USED HERE TO BE COMPATIBLE WITH
C        EXTERNAL COMPLEX FUN1 USED IN CALL HANKL AND HYBFHT.
C     REL1 - NAME OF AN EXTERNAL DECLARED SUBROUTINE WRITTEN BY THE USER
C            TO EVALUATE ANY (AND ALL) NREL RELATED KERNELS FJ(NREL)
C            AS A FUNCTION OF F1=FUN1(G), ASSUMED PREVIOUSLY EVALUATED.
C            AN EXTERNAL STATEMENT MUST APPEAR IN THE USER'S PROGRAM AS:
C                  EXTERNAL REL1
C            AND THE SUBROUTINE MUST HAVE THE GENERAL FORM:
C            (REPLACE "<...>" BELOW AS REQUIRED.)
C
C                  SUBROUTINE REL1(G,F1,FJ,J)
C                  COMPLEX F1,FJ(1)
C                  COMMON/PASS/<...INCLUDE PARAMETERS AS NEEDED...>
C                  IF(J.LT.0) GO TO (1,2,<...UP TO NREL...>), IABS(J)
C            1     FJ(1)=F1
C                    IF(J.LT.0) RETURN
C            2     FJ(2)=  <...EXPRESSION IN TERMS OF G,F1, ETC...>
C                    IF(J.LT.0) RETURN
C                  <...>
C                  RETURN
C                  END
C           NOTES:
C             SEE COMMENTS IN SUBROUTINES HYBFHT AND HANKL ON REL1.
C             REL1 PARAMETERS G,F1,J ARE CONTROLLED BY SUBROUTINE BESAUX
C             WHERE THE USER MUST ONLY CODE THE FJ(J),J=1,NREL RELATED
C             KERNELS AS A FUNCTION OF G,F1, AND ANY OTHER PARAMETERS
C             THAT CAN BE PASSED IN COMMON/PASS/, AS REQUIRED.
C     FWORK - COMPLEX WORK ARRAY FWORK(MAXREL), WHICH IS USED IN CALL
C             TO REL1 AND PASSED TO BESAUX BY USER'S CALLING PROGRAM.
C    MAXREL - MAXIMUM DIMENSION OF FWORK(MAXREL) IN CALLING PROGRAM.
C///////////////////////////////////////////////////////////////////////
C  RERR,AERR-RELATIVE AND ABSOLUTE ERROR FOR TERMINATION
C       BESAUX TERMINATES WHEN INCREASING THE GAUSS ORDER DOES NOT
C       CHANGE THE RESULT BY MORE THAN RERR OR WHEN THE ABSOLUTE ERROR
C       IS LESS THAN AERR OR WHEN A GAUSS ORDER OF NU IS REACHED.
C  NPCS-NUMBER OF PIECES INTO WHICH EACH PARTIAL INTEGRAND IS DIVIDED,
C       ORDINARILY SET TO ONE. FOR VERY SMALL VALUES OF R WHERE
C       THE KERNEL FUNCTION IS APPRECIABLE ONLY OVER THE FIRST FEW
C       LOOPS OF THE BESSEL FUNCTION, NPCS MAY BE INCREASED TO ACHIEVE
C       REASONABLE ACCURACY.
C  NEW  IF NEW=1, THE INTEGRANDS ARE COMPUTED AND SAVED AT EACH GAUSS
C       ORDER. IF NEW=2, PREVIOUSLY COMPUTED INTEGRANDS ARE USED. NOTE
C       THAT ORDER,R, AND NPCS MUST NOT BE CHANGED WHEN SETTING NEW=2.
C  IERR-ERROR PARAMETER
C       IERR=0--NORMAL RETURN
C       IERR=1--RESULT NOT ACCURATE TO RERR DUE TO TOO LOW A GAUSS
C       ORDER OR CONVERGENCE NOT ACHIEVED IN BESTRX
C
C  SUBROUTINES REQUIRED:
C    BESTRX,BESQUX,JBESS,PADECF,CF,ZEROJ,DOT
C
C  A.CHAVE IGPP/UCSD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XSUM(1)
      PARAMETER (NSUM=0)
C  COMMON /TEST/ CONTAINS RUN STATISTICS
C  NG=HIGHEST GAUSS ORDER USED
C  NF=TOTAL FUNCTION EVALUATIONS  FOR ALL CALLS TO BESTRX
C  NI=TOTAL NUMBER OF PARTIAL INTEGRANDS ON LAST CALL
C  TO BESTRX ONLY
      COMMON /TEST/NG,NF,NI
C//
      EXTERNAL COMPLEX FUN1 !20090922
      EXTERNAL REL1
      COMPLEX FWORK(MAXREL)
C//
      NF=0
      IF(NL.GT.NU)THEN
        BESR=0.
        BESI=0.
        IERR=1
        RETURN
      ENDIF
      NW=MAX(NEW,1)
      CALL BESTRX(BESR,BESI,ORDER,NL,R,FUN1,REL1,FWORK,MAXREL,
     $ .1*RERR,.1*AERR, NPCS,XSUM,NSUM,NW,IERR)
      IF((IERR.NE.0).AND.(NL.EQ.7))THEN
        NG=NL
        RETURN
      ELSE
        OLDR=BESR
        OLDI=BESI
        DO 10 N=NL+1,NU
          CALL BESTRX(BESR,BESI,ORDER,N,R,FUN1,REL1,FWORK,MAXREL,
     $     .1*RERR,.1*AERR, NPCS,XSUM,NSUM,2,IERR)
          IF((IERR.NE.0).AND.(N.EQ.7))THEN
            BESR=OLDR
            BESI=OLDI
            NG=N
            RETURN
          ELSEIF((ABS(BESR-OLDR).LE.RERR*ABS(BESR)+AERR).AND.(ABS(BESI-
     $           OLDI).LE.RERR*ABS(BESI)+AERR))THEN
            NG=N
            RETURN
          ELSE
            OLDR=BESR
            OLDI=BESI
          ENDIF
   10   CONTINUE
      ENDIF
      NG=7
      IERR=1
      RETURN
      END
      SUBROUTINE BESTRX(BESR,BESI,ORDER,NG,R,FUN1,REL1,FWORK,MAXREL,
     $  RERR,AERR,NPCS, XSUM,NSUM,NEW,IERR)
C
C  COMPUTES BESSEL TRANSFORM OF SPECIFIED ORDER DEFINED AS
C  INTEGRAL(FUNCT(X)*J-SUB-ORDER(X*R)*DX) FROM X=0 TO INFINITY
C  COMPUTATION IS ACHIEVED BY INTEGRATION BETWEEN THE ASYMPTOTIC
C  ZERO CROSSINGS OF THE BESSEL FUNCTION USING GAUSS QUADRATURE.
C  THE RESULTING SERIES OF PARTIAL INTEGRANDS IS SUMMED BY CALCULATING
C  THE PADE APPROXIMANTS TO SPEED UP CONVERGENCE.
C
C  ARGUMENT LIST:
C
C  BESR,BESI  REAL AND IMAGINARY PARTS RETURNED BY BESTRX
C  ORDER  ORDER OF THE BESSEL FUNCTION
C  NG  NUMBER OF GAUSS POINTS TO USE IN THE QUADRATURE ROUTINE.
C          NG=1 THROUGH 7 SELECTS 3,7,15,31,63,126,AND 255 TERMS.
C  R     ARGUMENT OF THE BESSEL FUNCTION
C///////////////////////////////////////////////////////////////////////
C   PARAMETERS "FUN1,REL1,FWORK,MAXREL" ARE INCLUDED SO THAT BESAUX CAN
C   BE RUN UNDER THE GENERAL CONTROL OF SUBROUTINE HYBFHT.
C   (SEE SUBROUTINES BESAUX AND HYBFHT FOR DETAILS ON THESE PARAMETERS.)
C///////////////////////////////////////////////////////////////////////
C  RERR,AERR  SPECIFIED RELATIVE AND ABSOLUTE ERROR FOR THE
C          CALCULATION.  THE INTEGRATION
C          TERMINATES WHEN AN ADDITIONAL TERM DOES NOT CHANGE THE
C          RESULT BY MORE THAN RERR*RESULT+AERR
C  NPCS  NUMBER OF PIECES INTO WHICH EACH PARTIAL INTEGRAND IS DIVIDED,
C        ORDINARILY SET TO ONE. FOR VERY SMALL VALUES OF RANGE WHERE
C        THE KERNEL FUNCTION IS APPRECIABLE ONLY OVER THE FIRST FEW
C        LOOPS OF THE BESSEL FUNCTION, NPCS MAY BE INCREASED TO ACHIEVE
C        REASONABLE ACCURACY. NOTE THAT NPCS AFFECTS ONLY THE PADE
C        SUM PORTION OF THE INTEGRATION, OVER X(NSUM) TO INFINITY.
C  XSUM  VECTOR OF VALUES OF THE KERNEL ARGUMENT OF FUNCT FOR WHICH
C        EXPLICIT CALCULATION OF THE INTEGRAL IS DESIRED, SO THAT THE
C        INTEGRAL OVER 0 TO XSUM(NSUM) IS ADDED TO THE INTEGRAL OVER
C        XSUM(NSUM) TO INFINITY WITH THE PADE METHOD INVOKED ONLY FOR
C        THE LATTER. THIS ALLOWS THE PADE SUMMATION METHOD TO BE
C        OVERRIDDEN AND SOME TYPES OF SINGULARITIES TO BE HANDLED.
C  NSUM  NUMBER OF VALUES IN XSUM, MAY BE ZERO.
C  NEW   DETERMINES METHOD OF KERNEL CALCULATION
C        NEW=0 MEANS  CALCULATE BUT DO NOT SAVE INTEGRANDS
C        NEW=1 MEANS CALCULATE KERNEL BY CALLING FUNCT-SAVE KERNEL
C              TIMES BESSEL FUNCTION
C        NEW=2 MEANS USE SAVED KERNELS TIMES BESSEL FUNCTIONS IN
C              COMMON /BESINT/. NOTE THAT ORDER,R,NPCS,XSUM, AND
C              NSUM MAY NOT BE CHANGED WHEN SETTING NEW=2.
C  IERR  ERROR PARAMETER
C        0 NORMAL RETURN-INTEGRAL CONVERGED
C        1 MEANS NO CONVERGENCE AFTER NSTOP TERMS IN THE PADE SUM
C
C  SUBROUTINES REQUIRED:
C  BESQUX,PADECF,CF,ZEROJ,DOT,JBESS
C
C  A.CHAVE IGPP/UCSD
C
C  NTERM IS MAXIMUM NUMBER OF BESSEL FUNCTION LOOPS STORED IF NEW.NE.0
C  NSTOP IS MAXIMUM NUMBER OF PADE TERMS
C
      PARAMETER (NTERM=50, NSTOP=100)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION KARG,KERN,LASTR,LASTI,KBES
      DIMENSION KARG(255,NTERM),KERN(510,NTERM),SR(NSTOP),SI(NSTOP),
     $          NK(NTERM),XSUM(1),KBES(510,NTERM)
C  COMMON /TEST/ CONTAINS RUN STATISTICS
C  NGAUSS=HIGHEST GAUSS ORDER USED
C  NF=TOTAL FUNCTION EVALUATIONS FOR ALL CALLS TO BESTRX
C  NI=TOTAL NUMBER OF PARTIAL INTEGRANDS USED IN PADECF ON LAST CALL
C  TO BESTRX
      COMMON /TEST/NGAUSS,NF,NI
C  IF NEW.GT.0, COMMON /BESINT/ CONTAINS THE KERNEL FUNCTION ARGUMENTS
C  AND VALUES KARG(I,N),KERN(I,N) FOR N=1 TO NPS AND I=1 TO
C  2**(NK(N)+1)-1.
C  KERN(I,N), KERN(I+1,N) CONTAIN THE REAL AND IMAGINARY PARTS OF THE
C  KERNEL FUNCTION FOR THE CORRESPONDING ARGUMENT VALUE KARG(I,N).
C  EACH VALUE OF N CORRESPONDS TO ONE ENTRY IN XSUM OR TO ONE INTERVAL
C  BETWEEN ZERO CROSSINGS OF THE BESSEL FUNCTION DIVIDED BY NPCS.
C  NOTE THAT UP TO NTERM VALUES ARE STORED THE INTEGRATION PROCEEDS
C  FOR N.GT.NTERM AS IF NEW=0. NP IS AN INDEX PASSED INTERNALLY BY
C  BESTRX. (KBES ADDED FOR HYBFHT CONTROL.)
      COMMON /BESINT/NK,NP,NPS,KARG,KERN,KBES
C//
      EXTERNAL COMPLEX FUN1 !20090922
      EXTERNAL REL1
      COMPLEX FWORK(MAXREL)
C//
      IF(NEW.EQ.2)THEN
        NPO=NPS
      ELSE
        DO 5 I=1,NTERM
          NK(I)=0
    5   CONTINUE
        NPS=0
        NPO=NTERM
      ENDIF
C  CHECK FOR TRIVIAL CASE
      IF((ORDER.NE.0.0).AND.(R.EQ.0.))THEN
        BESR=0.
        BESI=0.
        IERR=0
        RETURN
      ENDIF
      NI=0
      NW=NEW
      NP=1
      NPB=1
      L=1
      B=0.0
      SUMR=0.0
      SUMI=0.0
      XSUMR=0.0
      XSUMI=0.0
      IF(NSUM.GT.0)THEN
C  COMPUTE BESSEL TRANSFORM EXPLICITLY ON (0,XSUM(NSUM))
        LASTR=0.0
        LASTI=0.0
        DO 10 N=1,NSUM
          IF((NW.EQ.2).AND.(NP.GT.NPO))NW=1
          IF(NP.GT.NTERM)NW=0
          A=B
          B=XSUM(N)
          CALL BESQUX(A,B,TERMR,TERMI,NG,NW,ORDER,R,FUN1,REL1,
     1  FWORK,MAXREL)
          XSUMR=XSUMR+TERMR
          XSUMI=XSUMI+TERMI
          IF((ABS(XSUMR-LASTR).LE.RERR*ABS(XSUMR)+AERR).AND.
     $      (ABS(XSUMI-LASTI).LE.RERR*ABS(XSUMI)+AERR))THEN
            BESR=XSUMR
            BESI=XSUMI
            IERR=0
            NPS=MAX(NP,NPS)
            RETURN
          ELSE
            NP=NP+1
            LASTR=XSUMR
            LASTI=XSUMI
          ENDIF
   10   CONTINUE
C  FIND FIRST ZERO CROSSING OF BESSEL FUNCTION BEYOND XSUM(NSUM)
   15   CONTINUE
        IF(ZEROJ(NPB,ORDER).GT.XSUM(NSUM)*R)GOTO 20
        NPB=NPB+1
        GOTO 15
      ENDIF
C  ENTRY POINT FOR PADE SUMMATION OF PARTIAL INTEGRANDS
   20 CONTINUE
      LASTR=0.0
      LASTI=0.0
      A=B
      B=ZEROJ(NPB,ORDER)/R
C  CALCULATE TERMS AND SUM WITH PADECF, QUITTING WHEN CONVERGENCE IS
C  OBTAINED
      DO 30 N=1,NSTOP
        IF(NPCS.EQ.1)THEN
          IF((NW.EQ.2).AND.(NP.GT.NPO))NW=1
          IF(NP.GT.NTERM)NW=0
          CALL BESQUX(A,B,TERMR,TERMI,NG,NW,ORDER,R,FUN1,REL1,
     1  FWORK,MAXREL)
          NP=NP+1
        ELSE
          TERMR=0.
          TERMI=0.
          XINC=(B-A)/NPCS
          AA=A
          BB=A+XINC
          DO 25 I=1,NPCS
            IF((NW.EQ.2).AND.(NP.GT.NPO))NW=1
            IF(NP.GT.NTERM)NW=0
            CALL BESQUX(AA,BB,TR,TI,NG,NW,ORDER,R,FUN1,REL1,
     1    FWORK,MAXREL)
            TERMR=TERMR+TR
            TERMI=TERMI+TI
            AA=BB
            BB=BB+XINC
            NP=NP+1
   25     CONTINUE
        ENDIF
        NI=NI+1
        SR(L)=TERMR
        SI(L)=TERMI
        CALL PADECF(SUMR,SUMI,SR,SI,L)
        IF((ABS(SUMR-LASTR).LE.RERR*ABS(SUMR)+AERR).AND.
     $     (ABS(SUMI-LASTI).LE.RERR*ABS(SUMI)+AERR))THEN
          BESR=XSUMR+SUMR
          BESI=XSUMI+SUMI
          IERR=0
          NPS=MAX(NP-1,NPS)
          RETURN
        ELSE
          LASTR=SUMR
          LASTI=SUMI
          NPB=NPB+1
          A=B
          B=ZEROJ(NPB,ORDER)/R
          L=L+1
        ENDIF
   30 CONTINUE
      BESR=XSUMR+SUMR
      BESI=XSUMI+SUMI
      IERR=1
      NPS=MAX(NP,NPS)
      RETURN
      END
      SUBROUTINE BESQUX(A,B,BESR,BESI,NG,NEW,ORDER,R,FUN1,REL1,
     1 FWORK,MAXREL)
C
C  THIS SUBROUTINE CALCULATES THE INTEGRAL OF F(X)*J-SUB-N(X*R)
C  OVER THE INTERVAL A TO B AT A SPECIFIED GAUSS ORDER
C  THE RESULT IS OBTAINED USING A SEQUENCE OF 1,3,7,15,31,63,
C  127, AND 255 POINT INTERLACING GAUSS FORMULAE SO THAT NO INTEGRAND
C  EVALUATIONS ARE WASTED. THE KERNEL FUNCTIONS MAY BE SAVED SO THAT
C  BESSEL TRANSFORMS OF SIMILAR KERNELS ARE COMPUTED WITHOUT NEW
C  EVALUATION OF THE KERNEL.
C  DETAILS ON THE FORMULAE ARE GIVEN IN 'THE OPTIMUM ADDITION OF POINTS
C  TO QUADRATURE FORMULAE' BY T.N.L. PATTERSON, MATHS.COMP.
C  22,847-856 (1968). GAUSS WEIGHTS TAKEN FROM
C  COMM. A.C.M. 16,694-699 (1973)
C
C ARGUMENT LIST:
C
C A      LOWER LIMIT OF INTEGRATION
C B      UPPER LIMIT OF INTEGRATION
C BESR,BESI RETURNED INTEGRAL VALUE REAL AND IMAGINARY PARTS
C NG NUMBER OF POINTS IN THE GAUSS FORMULA.  NG=1,...7
C        SELECTS 3,7,15,31,63,127,AND 255 POINT QUADRATURE.
C NEW SELECTS METHOD OF KERNEL EVALUATION
C     NEW=0 CALCULATES KERNELS BY CALLING F - NOTHING SAVED
C     NEW=1 CALCULATES KERNELS BY CALLING F AND SAVES KERNEL TIMES
C            BESSEL FUNCTION IN COMMON /BESINT/
C     NEW=2 USES SAVED KERNEL TIMES BESSEL FUNCTIONS IN
C           COMMON /BESINT/
C ORDER ORDER OF THE BESSEL FUNCTION
C R    ARGUMENT OF THE BESSEL FUNCTION
C///////////////////////////////////////////////////////////////////////
C   PARAMETERS "FUN1,REL1,FWORK,MAXREL" ARE INCLUDED SO THAT BESAUX CAN
C   BE RUN UNDER THE GENERAL CONTROL OF SUBROUTINE HYBFHT.
C   (SEE SUBROUTINES BESAUX AND HYBFHT FOR DETAILS ON THESE PARAMETERS.)
C///////////////////////////////////////////////////////////////////////
C//
C//  NOTE:
C//  SINGLE-PRESISION COMPLEX FOR BESAUX AT FUN1 LEVEL ONLY/
C//
C//  (SEE USE OF ZFUN1 AND FUN1 BELOW AS BOUNDED BY C// LINES.)
C///////////////////////////////////////////////////////////////////////
C
C  A.CHAVE IGPP/UCSD
C
C  MAXIMUM NUMBER OF BESSEL FUNCTION LOOPS THAT CAN BE SAVED
      PARAMETER (NTERM=50)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C//
      COMPLEX ZFUN1,FWORK(MAXREL) !20090922
      COMMON/CONTRL/JPATH,JREL
      EXTERNAL COMPLEX FUN1 !20090922
      EXTERNAL REL1
C//
      DOUBLE PRECISION KARG,KERN,JBESS,KBES
      DIMENSION WT(254),WA(127),NWA(7),NWT(7),KARG(255,NTERM),
     $          KERN(510,NTERM),FUNCT(254),FR1(64),FI1(64),FR2(64),
     $          FI2(64),BES1(64),BES2(64),NK(NTERM),KBES(510,NTERM)
C  COMMON /TEST/ CONTAINS RUN STATISTICS
C  NGAUSS=HIGHEST GAUSS ORDER USED
C  NF=TOTAL FUNCTION EVALUATIONS  FOR ALL CALLS TO BESTRX
C  NI=TOTAL NUMBER OF PARTIAL INTEGRANDS ON LAST CALL
C  TO BESTRX ONLY
      COMMON /TEST/NGAUSS,NF,NI
      COMMON /BESINT/NK,NP,NPS,KARG,KERN,KBES
C  STARTING ADDRESSES IN ARRAYS WT AND WA AT EACH ORDER
      DATA NWT/1,3,7,15,31,63,127/,NWA/1,2,4,8,16,32,64/
C  ARRAY WT CONTAINS GAUSS QUADRATURE WEIGHTS STORED AS NEEDED AT
C  EACH ORDER.  ARRAY WA CONTAINS THE CORRESPONDING ARGUMENT WEIGHTS.
      DATA (WT(I),I=1,20)/
     $0.55555555555555555556D+00,0.88888888888888888889D+00,
     $0.26848808986833344073D+00,0.10465622602646726519D+00,
     $0.40139741477596222291D+00,0.45091653865847414235D+00,
     $0.13441525524378422036D+00,0.51603282997079739697D-01,
     $0.20062852937698902103D+00,0.17001719629940260339D-01,
     $0.92927195315124537686D-01,0.17151190913639138079D+00,
     $0.21915685840158749640D+00,0.22551049979820668739D+00,
     $0.67207754295990703540D-01,0.25807598096176653565D-01,
     $0.10031427861179557877D+00,0.84345657393211062463D-02,
     $0.46462893261757986541D-01,0.85755920049990351154D-01/
      DATA (WT(I),I=21,40)/
     $0.10957842105592463824D+00,0.25447807915618744154D-02,
     $0.16446049854387810934D-01,0.35957103307129322097D-01,
     $0.56979509494123357412D-01,0.76879620499003531043D-01,
     $0.93627109981264473617D-01,0.10566989358023480974D+00,
     $0.11195687302095345688D+00,0.11275525672076869161D+00,
     $0.33603877148207730542D-01,0.12903800100351265626D-01,
     $0.50157139305899537414D-01,0.42176304415588548391D-02,
     $0.23231446639910269443D-01,0.42877960025007734493D-01,
     $0.54789210527962865032D-01,0.12651565562300680114D-02,
     $0.82230079572359296693D-02,0.17978551568128270333D-01/
      DATA (WT(I),I=41,60)/
     $0.28489754745833548613D-01,0.38439810249455532039D-01,
     $0.46813554990628012403D-01,0.52834946790116519862D-01,
     $0.55978436510476319408D-01,0.36322148184553065969D-03,
     $0.25790497946856882724D-02,0.61155068221172463397D-02,
     $0.10498246909621321898D-01,0.15406750466559497802D-01,
     $0.20594233915912711149D-01,0.25869679327214746911D-01,
     $0.31073551111687964880D-01,0.36064432780782572640D-01,
     $0.40715510116944318934D-01,0.44914531653632197414D-01,
     $0.48564330406673198716D-01,0.51583253952048458777D-01,
     $0.53905499335266063927D-01,0.55481404356559363988D-01/
      DATA (WT(I),I=61,80)/
     $0.56277699831254301273D-01,0.56377628360384717388D-01,
     $0.16801938574103865271D-01,0.64519000501757369228D-02,
     $0.25078569652949768707D-01,0.21088152457266328793D-02,
     $0.11615723319955134727D-01,0.21438980012503867246D-01,
     $0.27394605263981432516D-01,0.63260731936263354422D-03,
     $0.41115039786546930472D-02,0.89892757840641357233D-02,
     $0.14244877372916774306D-01,0.19219905124727766019D-01,
     $0.23406777495314006201D-01,0.26417473395058259931D-01,
     $0.27989218255238159704D-01,0.18073956444538835782D-03,
     $0.12895240826104173921D-02,0.30577534101755311361D-02/
      DATA (WT(I),I=81,100)/
     $0.52491234548088591251D-02,0.77033752332797418482D-02,
     $0.10297116957956355524D-01,0.12934839663607373455D-01,
     $0.15536775555843982440D-01,0.18032216390391286320D-01,
     $0.20357755058472159467D-01,0.22457265826816098707D-01,
     $0.24282165203336599358D-01,0.25791626976024229388D-01,
     $0.26952749667633031963D-01,0.27740702178279681994D-01,
     $0.28138849915627150636D-01,0.50536095207862517625D-04,
     $0.37774664632698466027D-03,0.93836984854238150079D-03,
     $0.16811428654214699063D-02,0.25687649437940203731D-02,
     $0.35728927835172996494D-02,0.46710503721143217474D-02/
      DATA (WT(I),I=101,120)/
     $0.58434498758356395076D-02,0.70724899954335554680D-02,
     $0.83428387539681577056D-02,0.96411777297025366953D-02,
     $0.10955733387837901648D-01,0.12275830560082770087D-01,
     $0.13591571009765546790D-01,0.14893641664815182035D-01,
     $0.16173218729577719942D-01,0.17421930159464173747D-01,
     $0.18631848256138790186D-01,0.19795495048097499488D-01,
     $0.20905851445812023852D-01,0.21956366305317824939D-01,
     $0.22940964229387748761D-01,0.23854052106038540080D-01,
     $0.24690524744487676909D-01,0.25445769965464765813D-01,
     $0.26115673376706097680D-01,0.26696622927450359906D-01/
      DATA (WT(I),I=121,140)/
     $0.27185513229624791819D-01,0.27579749566481873035D-01,
     $0.27877251476613701609D-01,0.28076455793817246607D-01,
     $0.28176319033016602131D-01,0.28188814180192358694D-01,
     $0.84009692870519326354D-02,0.32259500250878684614D-02,
     $0.12539284826474884353D-01,0.10544076228633167722D-02,
     $0.58078616599775673635D-02,0.10719490006251933623D-01,
     $0.13697302631990716258D-01,0.31630366082226447689D-03,
     $0.20557519893273465236D-02,0.44946378920320678616D-02,
     $0.71224386864583871532D-02,0.96099525623638830097D-02,
     $0.11703388747657003101D-01,0.13208736697529129966D-01/
      DATA (WT(I),I=141,160)/
     $0.13994609127619079852D-01,0.90372734658751149261D-04,
     $0.64476204130572477933D-03,0.15288767050877655684D-02,
     $0.26245617274044295626D-02,0.38516876166398709241D-02,
     $0.51485584789781777618D-02,0.64674198318036867274D-02,
     $0.77683877779219912200D-02,0.90161081951956431600D-02,
     $0.10178877529236079733D-01,0.11228632913408049354D-01,
     $0.12141082601668299679D-01,0.12895813488012114694D-01,
     $0.13476374833816515982D-01,0.13870351089139840997D-01,
     $0.14069424957813575318D-01,0.25157870384280661489D-04,
     $0.18887326450650491366D-03,0.46918492424785040975D-03/
      DATA (WT(I),I=161,180)/
     $0.84057143271072246365D-03,0.12843824718970101768D-02,
     $0.17864463917586498247D-02,0.23355251860571608737D-02,
     $0.29217249379178197538D-02,0.35362449977167777340D-02,
     $0.41714193769840788528D-02,0.48205888648512683476D-02,
     $0.54778666939189508240D-02,0.61379152800413850435D-02,
     $0.67957855048827733948D-02,0.74468208324075910174D-02,
     $0.80866093647888599710D-02,0.87109650797320868736D-02,
     $0.93159241280693950932D-02,0.98977475240487497440D-02,
     $0.10452925722906011926D-01,0.10978183152658912470D-01,
     $0.11470482114693874380D-01,0.11927026053019270040D-01/
      DATA (WT(I),I=181,200)/
     $0.12345262372243838455D-01,0.12722884982732382906D-01,
     $0.13057836688353048840D-01,0.13348311463725179953D-01,
     $0.13592756614812395910D-01,0.13789874783240936517D-01,
     $0.13938625738306850804D-01,0.14038227896908623303D-01,
     $0.14088159516508301065D-01,0.69379364324108267170D-05,
     $0.53275293669780613125D-04,0.13575491094922871973D-03,
     $0.24921240048299729402D-03,0.38974528447328229322D-03,
     $0.55429531493037471492D-03,0.74028280424450333046D-03,
     $0.94536151685852538246D-03,0.11674841174299594077D-02,
     $0.14049079956551446427D-02,0.16561127281544526052D-02/
      DATA (WT(I),I=201,220)/
     $0.19197129710138724125D-02,0.21944069253638388388D-02,
     $0.24789582266575679307D-02,0.27721957645934509940D-02,
     $0.30730184347025783234D-02,0.33803979910869203823D-02,
     $0.36933779170256508183D-02,0.40110687240750233989D-02,
     $0.43326409680929828545D-02,0.46573172997568547773D-02,
     $0.49843645647655386012D-02,0.53130866051870565663D-02,
     $0.56428181013844441585D-02,0.59729195655081658049D-02,
     $0.63027734490857587172D-02,0.66317812429018878941D-02,
     $0.69593614093904229394D-02,0.72849479805538070639D-02,
     $0.76079896657190565832D-02,0.79279493342948491103D-02/
      DATA (WT(I),I=221,240)/
     $0.82443037630328680306D-02,0.85565435613076896192D-02,
     $0.88641732094824942641D-02,0.91667111635607884067D-02,
     $0.94636899938300652943D-02,0.97546565363174114611D-02,
     $0.10039172044056840798D-01,0.10316812330947621682D-01,
     $0.10587167904885197931D-01,0.10849844089337314099D-01,
     $0.11104461134006926537D-01,0.11350654315980596602D-01,
     $0.11588074033043952568D-01,0.11816385890830235763D-01,
     $0.12035270785279562630D-01,0.12244424981611985899D-01,
     $0.12443560190714035263D-01,0.12632403643542078765D-01,
     $0.12810698163877361967D-01,0.12978202239537399286D-01/
      DATA (WT(I),I=241,254)/
     $0.13134690091960152836D-01,0.13279951743930530650D-01,
     $0.13413793085110098513D-01,0.13536035934956213614D-01,
     $0.13646518102571291428D-01,0.13745093443001896632D-01,
     $0.13831631909506428676D-01,0.13906019601325461264D-01,
     $0.13968158806516938516D-01,0.14017968039456608810D-01,
     $0.14055382072649964277D-01,0.14080351962553661325D-01,
     $0.14092845069160408355D-01,0.14094407090096179347D-01/
      DATA (WA(I),I=1,20)/
     $0.77459666924148337704D+00,0.96049126870802028342D+00,
     $0.43424374934680255800D+00,0.99383196321275502221D+00,
     $0.88845923287225699889D+00,0.62110294673722640294D+00,
     $0.22338668642896688163D+00,0.99909812496766759766D+00,
     $0.98153114955374010687D+00,0.92965485742974005667D+00,
     $0.83672593816886873550D+00,0.70249620649152707861D+00,
     $0.53131974364437562397D+00,0.33113539325797683309D+00,
     $0.11248894313318662575D+00,0.99987288812035761194D+00,
     $0.99720625937222195908D+00,0.98868475754742947994D+00,
     $0.97218287474858179658D+00,0.94634285837340290515D+00/
      DATA (WA(I),I=21,40)/
     $0.91037115695700429250D+00,0.86390793819369047715D+00,
     $0.80694053195021761186D+00,0.73975604435269475868D+00,
     $0.66290966002478059546D+00,0.57719571005204581484D+00,
     $0.48361802694584102756D+00,0.38335932419873034692D+00,
     $0.27774982202182431507D+00,0.16823525155220746498D+00,
     $0.56344313046592789972D-01,0.99998243035489159858D+00,
     $0.99959879967191068325D+00,0.99831663531840739253D+00,
     $0.99572410469840718851D+00,0.99149572117810613240D+00,
     $0.98537149959852037111D+00,0.97714151463970571416D+00,
     $0.96663785155841656709D+00,0.95373000642576113641D+00/
      DATA (WA(I),I=41,60)/
     $0.93832039777959288365D+00,0.92034002547001242073D+00,
     $0.89974489977694003664D+00,0.87651341448470526974D+00,
     $0.85064449476835027976D+00,0.82215625436498040737D+00,
     $0.79108493379984836143D+00,0.75748396638051363793D+00,
     $0.72142308537009891548D+00,0.68298743109107922809D+00,
     $0.64227664250975951377D+00,0.59940393024224289297D+00,
     $0.55449513263193254887D+00,0.50768775753371660215D+00,
     $0.45913001198983233287D+00,0.40897982122988867241D+00,
     $0.35740383783153215238D+00,0.30457644155671404334D+00,
     $0.25067873030348317661D+00,0.19589750271110015392D+00/
      DATA (WA(I),I=61,80)/
     $0.14042423315256017459D+00,0.84454040083710883710D-01,
     $0.28184648949745694339D-01,0.99999759637974846462D+00,
     $0.99994399620705437576D+00,0.99976049092443204733D+00,
     $0.99938033802502358193D+00,0.99874561446809511470D+00,
     $0.99780535449595727456D+00,0.99651414591489027385D+00,
     $0.99483150280062100052D+00,0.99272134428278861533D+00,
     $0.99015137040077015918D+00,0.98709252795403406719D+00,
     $0.98351865757863272876D+00,0.97940628167086268381D+00,
     $0.97473445975240266776D+00,0.96948465950245923177D+00,
     $0.96364062156981213252D+00,0.95718821610986096274D+00/
      DATA (WA(I),I=81,100)/
     $0.95011529752129487656D+00,0.94241156519108305981D+00,
     $0.93406843615772578800D+00,0.92507893290707565236D+00,
     $0.91543758715576504064D+00,0.90514035881326159519D+00,
     $0.89418456833555902286D+00,0.88256884024734190684D+00,
     $0.87029305554811390585D+00,0.85735831088623215653D+00,
     $0.84376688267270860104D+00,0.82952219463740140018D+00,
     $0.81462878765513741344D+00,0.79909229096084140180D+00,
     $0.78291939411828301639D+00,0.76611781930376009072D+00,
     $0.74869629361693660282D+00,0.73066452124218126133D+00,
     $0.71203315536225203459D+00,0.69281376977911470289D+00/
      DATA (WA(I),I=101,120)/
     $0.67301883023041847920D+00,0.65266166541001749610D+00,
     $0.63175643771119423041D+00,0.61031811371518640016D+00,
     $0.58836243444766254143D+00,0.56590588542365442262D+00,
     $0.54296566649831149049D+00,0.51955966153745702199D+00,
     $0.49570640791876146017D+00,0.47142506587165887693D+00,
     $0.44673538766202847374D+00,0.42165768662616330006D+00,
     $0.39621280605761593918D+00,0.37042208795007823014D+00,
     $0.34430734159943802278D+00,0.31789081206847668318D+00,
     $0.29119514851824668196D+00,0.26424337241092676194D+00,
     $0.23705884558982972721D+00,0.20966523824318119477D+00/
      DATA (WA(I),I=121,127)/
     $0.18208649675925219825D+00,0.15434681148137810869D+00,
     $0.12647058437230196685D+00,0.98482396598119202090D-01,
     $0.70406976042855179063D-01,0.42269164765363603212D-01,
     $0.14093886410782462614D-01/
C  CHECK FOR TRIVIAL CASE
      IF(A.GE.B)THEN
        BESR=0.
        BESI=0.
        RETURN
      ENDIF
      IF(NEW.EQ.2)GOTO 200
C  SCALE FACTORS
      SUM=(B+A)/2.
      DIFF=(B-A)/2.
C  ONE POINT GAUSS
C//
C//      CALL F(SUM,FZEROR,FZEROI)
      ZFUN1=FUN1(SNGL(SUM))
      IF(NEW.GE.1) THEN
        CALL REL1(SNGL(SUM),ZFUN1,FWORK,-JREL)
        ZFUN1=FWORK(JREL)
      ENDIF
      FZEROR=REAL(ZFUN1)
      FZEROI=AIMAG(ZFUN1)
C//
      NF=NF+1
      BESF=JBESS(SUM*R,ORDER)
      FZEROR=BESF*FZEROR
      FZEROI=BESF*FZEROI
      IF(NEW.EQ.1)THEN
        KARG(1,NP)=SUM
        KERN(1,NP)=FZEROR
        KERN(2,NP)=FZEROI
C//
      KBES(1,NP)=BESF
C//
        LA=2
        LK=3
      ENDIF
      N=1
      KN=1
    5 CONTINUE
C  STEP THROUGH GAUSS ORDERS
      DO 40 K=KN,NG
C  COMPUTE NEW FUNCTION VALUES
        NA=NWA(K)
        DO 10 J=1,NWA(K)
          X=WA(NA)*DIFF
          NA=NA+1
          SUMP=SUM+X
          SUMM=SUM-X
C//
C//          CALL F(SUMP,FR1(J),FI1(J))
        ZFUN1=FUN1(SNGL(SUMP))
      IF(NEW.GE.1) THEN
        CALL REL1(SNGL(SUMP),ZFUN1,FWORK,-JREL)
        ZFUN1=FWORK(JREL)
      ENDIF
        FR1(J)=REAL(ZFUN1)
        FI1(J)=AIMAG(ZFUN1)
C//          CALL F(SUMM,FR2(J),FI2(J))
        ZFUN1=FUN1(SNGL(SUMM))
      IF(NEW.GE.1) THEN
        CALL REL1(SNGL(SUMM),ZFUN1,FWORK,-JREL)
        ZFUN1=FWORK(JREL)
      ENDIF
        FR2(J)=REAL(ZFUN1)
        FI2(J)=AIMAG(ZFUN1)
C//
          NF=NF+2
          BES1(J)=JBESS(SUMP*R,ORDER)
          BES2(J)=JBESS(SUMM*R,ORDER)
          IF(NEW.GE.1)THEN
            KARG(LA,NP)=SUMP
            KARG(LA+1,NP)=SUMM
C//
          KBES(LA,NP)=BES1(J)
          KBES(LA+1,NP)=BES2(J)
C//
            LA=LA+2
          ENDIF
   10   CONTINUE
C  COMPUTE PRODUCTS OF KERNELS AND BESSEL FUNCTIONS
        DO 20 J=1,NWA(K)
          FR1(J)=BES1(J)*FR1(J)
          FI1(J)=BES1(J)*FI1(J)
          FR2(J)=BES2(J)*FR2(J)
          FI2(J)=BES2(J)*FI2(J)
          FUNCT(N)=FR1(J)+FR2(J)
          FUNCT(N+1)=FI1(J)+FI2(J)
          N=N+2
   20   CONTINUE
        IF(NEW.GE.1)THEN
          DO 30 J=1,NWA(K)
            KERN(LK,NP)=FR1(J)
            KERN(LK+1,NP)=FI1(J)
            KERN(LK+2,NP)=FR2(J)
            KERN(LK+3,NP)=FI2(J)
            LK=LK+4
   30     CONTINUE
        ENDIF
   40 CONTINUE
C  COMPUTE DOT PRODUCT OF WEIGHTS WITH INTEGRAND VALUES
      NW=NWT(NG)
C  DOT SHOULD BE REPLACED WITH SDOT FROM SCILIB ON CRAY-1 IF D.P. IS
C  NOT USED
      ACUMR=DOT(NW,WT(NW),1,FUNCT(1),2)
      ACUMI=DOT(NW,WT(NW),1,FUNCT(2),2)
      BESR=(ACUMR+WT(2*NW)*FZEROR)*DIFF
      BESI=(ACUMI+WT(2*NW)*FZEROI)*DIFF
      IF(NP.LE.NTERM)NK(NP)=NG
      RETURN
  200 CONTINUE
C  CONSTRUCT FUNCT FROM SAVED KERNELS
      DIFF=(B-A)/2.
      FZEROR=KERN(1,NP)
      FZEROI=KERN(2,NP)
      K=MIN(NK(NP),NG)
      NW=NWT(K)
      LK=3
      DO 210 N=1,2*NW,2
        FUNCT(N)=KERN(LK,NP)+KERN(LK+2,NP)
        FUNCT(N+1)=KERN(LK+1,NP)+KERN(LK+3,NP)
        LK=LK+4
  210 CONTINUE
      IF(NK(NP).GE.NG)THEN
C  NO ADDITIONAL ORDERS REQUIRED-COMPUTE DOT PRODUCT OF WEIGHTS WITH
C  INTEGRAND VALUES AND RETURN
        ACUMR=DOT(NW,WT(NW),1,FUNCT(1),2)
        ACUMI=DOT(NW,WT(NW),1,FUNCT(2),2)
        BESR=(ACUMR+WT(2*NW)*FZEROR)*DIFF
        BESI=(ACUMI+WT(2*NW)*FZEROI)*DIFF
        RETURN
      ELSE
C  COMPUTE ADDITIONAL ORDERS BEFORE TAKING DOT PRODUCT
        SUM=(B+A)/2.
        KN=K+1
        LA=LK/2+1
        GOTO 5
      ENDIF
      END
      SUBROUTINE PADECF(SUMR,SUMI,SR,SI,N)
C
C  COMPUTES SUM(S(I)),I=1,...N BY COMPUTATION OF PADE APPROXIMANT
C  USING CONTINUED FRACTION EXPANSION.  FUNCTION IS DESIGNED TO BE
C  CALLED SEQUENTIALLY AS N IS INCREMENTED FROM 1 TO ITS FINAL VALUE.
C  THE NTH CONTINUED FRACTION COEFFICIENT IS CALCULATED AND
C  STORED AND THE NTH CONVERGENT RETURNED.  IT IS UP TO THE USER TO
C  STOP THE CALCULATION WHEN THE DESIRED ACCURACY IS ACHIEVED.
C  ALGORITHM FROM HANGGI ET AL., Z.NATURFORSCH. 33A,402-417 (1977)
C  IN THEIR NOTATION, VECTORS CFCOR,CFCOI ARE LOWER CASE D, VECTORS DR,
C  DI ARE UPPER CASE D, VECTORS XR,XI ARE X, AND VECTORS SR,SI ARE S
C
C  A.CHAVE IGPP/UCSD
C
      PARAMETER (NC=100)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE XR,XI,DR,DI,CFCOR,CFCOI
      DIMENSION SR(1),SI(1),XR(NC),XI(NC),DR(0:NC),DI(0:NC),CFCOR(NC),
     $          CFCOI(NC)
      DATA DR(0)/-1./,DI(0)/0./
      IF(N.LT.3)THEN
C  INITIALIZE FOR RECURSIVE CALCULATIONS
        IF(N.EQ.1)THEN
          DO 10 I=1,NC
            XR(I)=0.
            XI(I)=0.
   10     CONTINUE
        ENDIF
        DR(N)=SR(N)
        DI(N)=SI(N)
        DENOM=DR(N-1)**2+DI(N-1)**2
        CFCOR(N)=-(DR(N-1)*DR(N)+DI(N-1)*DI(N))/DENOM
        CFCOI(N)=-(DR(N-1)*DI(N)-DR(N)*DI(N-1))/DENOM
        CALL CF(SUMR,SUMI,CFCOR,CFCOI,N)
        RETURN
      ELSE
        L=2*INT((N-1)/2)
C  UPDATE X VECTORS FOR RECURSIVE CALCULATION OF CF COEFFICIENTS
        DO 20 K=L,4,-2
          XR(K)=XR(K-1)+CFCOR(N-1)*XR(K-2)-CFCOI(N-1)*XI(K-2)
          XI(K)=XI(K-1)+CFCOR(N-1)*XI(K-2)+CFCOI(N-1)*XR(K-2)
   20   CONTINUE
        XR(2)=XR(1)+CFCOR(N-1)
        XI(2)=XI(1)+CFCOI(N-1)
C  INTERCHANGE ODD AND EVEN PARTS
        DO 30 K=1,MAX(1,L-1),2
          T1=XR(K)
          T2=XI(K)
          XR(K)=XR(K+1)
          XI(K)=XI(K+1)
          XR(K+1)=T1
          XI(K+1)=T2
   30   CONTINUE
C  COMPUTE FIRST COEFFICIENTS
        DR(N)=SR(N)
        DI(N)=SI(N)
        DO 40 K=1,MAX(1,L/2)
          DR(N)=DR(N)+SR(N-K)*XR(2*K-1)-SI(N-K)*XI(2*K-1)
          DI(N)=DI(N)+SI(N-K)*XR(2*K-1)+SR(N-K)*XI(2*K-1)
   40   CONTINUE
C  COMPUTE NEW CF COEFFICIENT
        DENOM=DR(N-1)**2+DI(N-1)**2
        CFCOR(N)=-(DR(N)*DR(N-1)+DI(N)*DI(N-1))/DENOM
        CFCOI(N)=-(DR(N-1)*DI(N)-DR(N)*DI(N-1))/DENOM
C  EVALUATE CONTINUED FRACTION
        CALL CF(SUMR,SUMI,CFCOR,CFCOI,N)
        RETURN
      ENDIF
      END
      SUBROUTINE CF(RESR,RESI,CFCOR,CFCOI,N)
C
C  EVALUATES A COMPLEX CONTINUED FRACTION BY RECURSIVE DIVISION
C  STARTING AT THE BOTTOM, AS USED BY PADECF
C  RESR,RESI ARE REAL AND IMAGINARY PARTS RETURNED
C  CFCOR,CFCOI ARE REAL AND IMAGINARY VECTORS OF CONTINUED FRACTION
C  COEFFICIENTS
C
      PARAMETER (NC=100,NCM=NC-1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CFCOR(1),CFCOI(1),ONE(NC)
      DATA ONE/0.,NCM*1./
      RESR=ONE(N)+C
0 Kudos
1 Solution
Steven_L_Intel1
Employee
‎09-21-2009 09:03 PM
2,246 Views
Solved Jump to solution
The statement:

EXTERNAL COMPLEX FUN1

doesn't do what you think it does. It declares a name COMPLEXFUN1 as EXTERNAL. You could write:

COMPLEX,EXTERNAL :: FUN1

or

COMPLEX FUN1
EXTERNAL FUN1

You only get the error when generated interface checking is enabled.

View solution in original post

0 Kudos
Copy link

Link Copied

4 Replies
Leigh_Wardle
Beginner
‎09-21-2009 07:06 PM
2,246 Views
Solved Jump to solution
Quoting - Leigh Wardle
[plain] Hi,

I get the following compilation error on Lines 962, 997 and 1009 of SUBROUTINE HYBFHT (copy below):

error #6637: This actual argument must be the name of an external user function or the name of an intrinsic function.   [FUN1]

Each of these lines is of the form:

          CALL BESQUX(A,B,TERMR,TERMI,NG,NW,ORDER,R,FUN1,REL1,
     1  FWORK,MAXREL)

But as far as I can see I have declared FUN1 both in the calling routine (HYBFHT) and the routine being called (BESQUX).
The declaration in both cases is:

EXTERNAL COMPLEX FUN1

Any suggestions would be most appreciated.

Regards,
Leigh


[/plain]

I forgot to mention:

  • This code compiled succesfully with v10.1 - but in all declarations FUN1 was declared with the attribute COMPLEX.
  • I am now compiling with the latest version of the compiler, v.11.1.046 Update 2.
0 Kudos
Copy link
Steven_L_Intel1
Employee
‎09-21-2009 09:03 PM
2,247 Views
Solved Jump to solution
The statement:

EXTERNAL COMPLEX FUN1

doesn't do what you think it does. It declares a name COMPLEXFUN1 as EXTERNAL. You could write:

COMPLEX,EXTERNAL :: FUN1

or

COMPLEX FUN1
EXTERNAL FUN1

You only get the error when generated interface checking is enabled.
0 Kudos
Copy link
Leigh_Wardle
Beginner
‎09-23-2009 12:16 AM
2,246 Views
Solved Jump to solution
The statement:

EXTERNAL COMPLEX FUN1

doesn't do what you think it does. It declares a name COMPLEXFUN1 as EXTERNAL. You could write:

COMPLEX,EXTERNAL :: FUN1

or

COMPLEX FUN1
EXTERNAL FUN1

You only get the error when generated interface checking is enabled.

Thanks, Steve, for setting me straight.

I am fascinated that the Fortran declaration EXTERNAL COMPLEX FUN1 declares a name COMPLEXFUN1 as EXTERNAL.

Regards,
Leigh
0 Kudos
Copy link
Steven_L_Intel1
Employee
‎09-23-2009 07:03 AM
2,246 Views
Solved Jump to solution
Quoting - Leigh Wardle
I am fascinated that the Fortran declaration EXTERNAL COMPLEX FUN1 declares a name COMPLEXFUN1 as EXTERNAL.

It wouldn't if you used free-form source, which I recommend to all. I recognize that if you are dealing with existing fixed-form source that it may not make sense to change it, but in free form source this would have been caught as an error by the compiler.

In fixed-form sources, blanks are ignored in nearly all contexts. A very famous error, said to have brought down a spacecraft, was coding a DO loop like this:

DO I=1.10

where a dot instead of a comma was used. In fixed-form source, this is the same as:

DOI = 1.10

IMPLICIT NONE would also help you here.
0 Kudos
Copy link
Reply

Community support is provided Monday to Friday. Other contact methods are available here.

Intel does not verify all solutions, including but not limited to any file transfers that may appear in this community. Accordingly, Intel disclaims all express and implied warranties, including without limitation, the implied warranties of merchantability, fitness for a particular purpose, and non-infringement, as well as any warranty arising from course of performance, course of dealing, or usage in trade.

For more complete information about compiler optimizations, see our Optimization Notice.