      SUBROUTINE GWLOOKU3
C
C     * JUN 15/2013 - M.LAZARE.   NEW VERSION FOR GCM17+:
C     *                           - NUMERICAL RECIPE FUNCTIONS
C     *                             WHICH THIS ROUTINE USES
C     *                             HAVE BEEN UPGRADED FOR
C     *                             ACCURACY, SO NEW NAMES
C     *                             CALLED {EF3,KF3,RF3,RD3}.
C     * SEP 11/2006 - F.MAJAESS.  PREVIOUS VERSION FOR GCM15F:
C     *                           - HARD-CODED TO REAL*8 AND
C     *                             CALLS NEW FUNCTION NAMES
C     *                             NOW IS "MODEL" INSTEAD OF
C     *                             "COMM", TREATED SIMILARILY.
C     * MAR 24/1999 - J SCINOCCA. ORIGINAL VERSION GWLOOKU.
C
C     * THIS SUBROUTINE FILLS TWO LOOK-UP TABLES DP AND TP
C     * USED BY GWDFL14. DP AND TP ARE DERIVED BY NUMERICALLY
C     * INTEGRATING INCOMPLETE ELLIPTIC INTEGRALS.  THE NUMERICAL
C     * INTEGRATION IS PERFORMED USING THE FUNCTIONS EF2,KF2,RF2 
C     * AND RD2 OBTAINED FROM NUMERICAL RECIPES.
C     * OTHER QUANTITIES OUTPUT THROUGH THE COMMON GROUP GWDTABL ARE
C     * PI - ONE HALF THE CIRCUMFERENCE OF THE UNIT CIRCLE 
C     * NPSI,NGAM - THE DIMENSIONS OF DP AND TP
C     * DPSI,DGAM - GRID INCREMENTS IN EACH DIRECTION
C     * PSIST,GAMST - STARTING VALUE OF EACH DIMENSION
C
      IMPLICIT REAL*8 (A-H,O-Z),
     +INTEGER (I-N)

      COMMON /GWDTABL/PI,NPSI,NGAM,DPSI,DGAM,PSIST,GAMST,
     1                DP(100,100),TP(100,100)
      REAL*8 PI,DPSI,DGAM,PSIST,GAMST,DP,TP

      REAL*8 D3P,D3M,T3P,T3M,BP,BM,CP,CM,BB,CC,BPI,CPI
      REAL*8 CSP,SNP,FCT1,FCT2,GAMS
      REAL*8 EF3,KF3,RF3,RD3

      EXTERNAL EF3,KF3,RF3,RD3
C------------------------------------------------------------------
      PI=2*ASIN(1.)
      NPSI=100
      NGAM=100
      GAMST=0.
      DGAM=0.999/FLOAT(NGAM-1)
      PSIST=-PI/2
      DPSI=PI/FLOAT(NPSI-1)

      DO I=1,NGAM
        GAM=REAL(I-1)*DGAM+GAMST
        GAMS=GAM**2
        DO J=1,NPSI
          PSI=REAL(J-1)*DPSI+PSIST
          CSP=COS(PSI)
          SNP=SIN(PSI)
          FCT1=SQRT(CSP**2+GAMS*SNP**2)
          FCT2=SQRT(SNP**2+GAMS*CSP**2)
          APSI=ABS(PSI)
          BA1=SQRT(1.-GAM**2)
          B1=(EF3(BA1,APSI)-KF3(BA1,APSI)*GAM**2)/(BA1**2)
          B2=(EF3(BA1,PI/2.-APSI)-KF3(BA1,PI/2.-APSI)*GAM**2)
     1       /(BA1**2) 
          BB=B1+B2
          BPI=(EF3(BA1,PI/2.)-KF3(BA1,PI/2.)*GAM**2)/(BA1**2)
          C1=(KF3(BA1,APSI)-EF3(BA1,APSI))/(BA1**2)
          C2=(KF3(BA1,PI/2.-APSI)-EF3(BA1,PI/2.-APSI))/(BA1**2)
          CC=C1+C2
          CPI=(KF3(BA1,PI/2.)-EF3(BA1,PI/2.))/(BA1**2)
          IF(PSI.GE.0.) THEN
            BP=BB
            CP=CC
            BM=2.*BPI-BB
            CM=2.*CPI-CC
          ELSE
            BP=2.*BPI-BB
            CP=2.*CPI-CC
            BM=BB
            CM=CC
          ENDIF
          DP(I,J)=CSP**2*BP+GAMS*SNP**2*CP
     1           -2.*GAM*CSP*SNP*(FCT2-FCT1)/(GAMS-1.)
          TP(I,J)=CSP*SNP*BP-GAMS*CSP*SNP*CP
     1           +GAM*(CSP**2-SNP**2)*(FCT2-FCT1)/(GAMS-1.)
        ENDDO
      ENDDO

      RETURN
      END 
