      SUBROUTINE INITAL
C...
C...  TWO DIMENSIONAL INTEGRATION
C...
C...  THE FOLLOWING CODING ILLUSTRATES THE EVALUATION OF AN INTEGRAL
C...  OF THE FORM
C...
C...          YU    XU
C...     I = INT ( INT F(X,Y) DX) DY
C...          YL    XL
C...
C...  WHERE
C...
C...     I       INTEGRAL TO BE EVALUATED
C...
C...     F(X,Y)  GIVEN INTEGRAND
C...
C...     YU,YL   UPPER AND LOWER LIMITS FOR Y
C...
C...     XU,XL   UPPER AND LOWER LIMITS FOR X
C...
C...  XU AND XL ARE GENERALLY OF THE FORM XU = F(Y), XL = G(Y).
C...
C...  IN THE FOLLOWING CODING, THE INTEGRATION IN Y IS DONE BY THE USUAL
C...  ODE INTEGRATION THROUGH COMMON/T/, /Y/ AND /F/.  THE INTEGRATION IN
C...  X IS DONE BY SUBROUTINE INTX CALLED IN DERV.  INTX IMPLEMENTS THE
C...  SECOND ORDER ADAMS BASHFORTH METHOD WITH STARTING BY THE MODIFIED
C...  EULER METHOD.
C...
C...  F(X,Y) IS DEFINED IN A USER-SUPPLIED FUNCTION CALL BY INTX. TO
C...  TEST THE METHOD, A SAMPLE PROBLEM IS CODED FOR F(X,Y) = 1 OVER
C...  A QUARTER OF THE CIRCLE WITH UNIT RADIUS, FOR WHICH THE CORRECT
C...  ANSWER (THE AREA OF A QUARTER CIRCLE) IS PI/4.  THIS RESULT CAN
C...  BE DEMONSTRATED IN THE FOLLOWING WAY, WHICH GIVES SOME INSIGHT
C...  INTO THE OPERATION OF THE FOLLOWING CODING (XU = (1 - Y**2)**0.5
C...  FOR A CIRCLE WITH THE EQUATION X**2 + Y**2 = 1)
C...
C...       1     XU             1                 1
C...  I = INT ( INT 1 DX) DY = INT (XU - 0) DY = INT (1 - Y**2)**0.5 DY
C...       0     0              0                 0
C...
C...                                          1
C...    = (1/2)(Y*(1 - Y**2)**0.5 + ARCSIN(Y))  = PI/4
C...                                          0
      COMMON/T/Y
     +      /Y/I
     +      /F/DIDT
      REAL I
C...
C...  INITIAL CONDITION
      I=0.
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/Y
     +      /Y/I
     +      /F/DIDT
      REAL I
C...
C...  EXTERNAL THE FUNCTION CALLED BY INTX TO EVALUATE THE INTEGRAND
      EXTERNAL FXY
C...
C...  EVALUATE THE INTEGRAND
C...     LIMITS FOR X (BECAUSE OF ROUNDING, Y CAN BECOME SLIGHTLY GREATER
C...     THAN ONE WHICH CAUSES PROBLEMS WITH THE SQUARE ROOT UNLESS THE
C...     ABS IS USED)
            XL=0.
            XU=SQRT(ABS(1.0-Y**2))
C...
C...     NUMBER OF INTEGRATION INTERVALS IN X
            NX=100
C...
C...     INTEGRATION IN X
            CALL INTX(FXY,XL,XU,NX,Y,SUMX)
C...
C...     INTEGRAL IN X BECOMES THE INTEGRAND FOR THE Y INTEGRATION
            DIDT=SUMX
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/Y
     +      /Y/I
     +      /F/DIDT
      REAL I
C...
C...  WRITE THE SOLUTION (AFTER MULTIPLICATION BY 4 SO THE ANSWER IS
C...  PI, THE AREA OF A CIRCLE WITH UNIT RADIUS)
      WRITE(NO,1)Y,4.0*I
1     FORMAT(//,5H Y = ,F3.1,3X,4HI = ,F8.6)
      RETURN
      END
      FUNCTION FXY(X,Y)
C...
C...  FUNCTION FXY(X,Y) DEFINES THE INTEGRAND AS A FUNCTION OF X AND
C...  Y (IN THE FOLLOWING EXAMPLE, THE INTEGRAND IS CONSTANT, BUT IN
C...  GENERAL, IT WILL BE A FUNCTION OF THE INPUTS X AND Y)
C...
      FXY=1.0
      RETURN
      END
      SUBROUTINE INTX(FXY,XL,XU,NX,Y,SUMX)
C...
C...  SUBROUTINE INTX DOES A ONE DIMENSIONAL INTEGRATION OVER NX IN-
C...  TERVALS IN X AT A GIVEN VALUE OF Y.  THE COMPUTED INTEGRAL IS
C...  RETURNED AS SUMX.
C...
C...  EXTERNAL THE ROUTINE CALLED BY INTX TO EVALUATE THE INTEGRAND
      EXTERNAL FXY
C...
C...  SET THE INTEGRATION INTERVAL
      DX=(XU-XL)/FLOAT(NX)
C...
C...  BEGIN THE INTEGRATION OVER NX STEPS
C...
C...  NSTART = 1 - SECOND ORDER, RUNGE KUTTA (MODIFIED EULER) TO START
C...  THE INTEGRATION
      NSTART=1
      SUMX=0.
      DO 1 I=1,NX
C...
C...  TAKE A SECOND ORDER (NSTART = 1) RUNGE KUTTA STEP AT THE BEGINNING
C...  SINCE THE ADAMS BASHFORTH METHOD IS NOT SELF-STARTING
      IF(NSTART.EQ.1)THEN
C...
C...     FIRST ORDER, RUNGE KUTTA (EULER) STEP FROM N-1 = 1-1 = 0 TO
C...     N = 1
         FNM2=FXY(XL,Y)
         SUMX=SUMX+FNM2*DX
         X=XL+DX
C...
C...     THE DERIVATIVE AT N = 1 BECOMES THE DERIVATIVE AT N-1 FOR
C...     THE FIRST STEP OF THE ADAMS BASHFORTH INTEGRATION
         FNM1=FXY(X,Y)
C...
C...     SECOND ORDER, RUNGE KUTTA STEP
         SUMX=SUMX+0.5*(FNM1-FNM2)*DX
      END IF
C...
C...  TAKE A SECOND ORDER, ADAMS BASHFORTH STEP FROM N-1 TO N
      IF(NSTART.EQ.2)THEN
         SUMX=SUMX+DX*(1.5*FNM1-0.5*FNM2)
         X=X+DX
C...
C...     DERIVATIVES AT N-1 AND N-2 FOR NEXT STEP
         FNM2=FNM1
         FNM1=FXY(X,Y)
      END IF
      NSTART=2
C...
C...  THE FOLLOWING WRITE WAS USED TO CHECK THE OPERATION OF SUBROUTINE
C...  INTX.  WHEN ACTIVATED, IT ALSO GIVES A COMPLETE PICTURE OF HOW THE
C...  TWO DIMENSIONAL INTEGRATION WORKS.  HOWEVER, A WORD OF CAUTION -
C...  THIS WRITE WILL BE EXECUTED MANY TIMES UNLESS THE NUMBER OF INCRE-
C...  MENTS IN X AND Y IS SMALL, E.G., 10 IN EACH DIRECTION RESULTS IN
C...  100 CALLS TO INTX
C...  WRITE(*,*)' I = ',I,' Y = ',Y,' X = ',X,' SUMX = ',SUMX
C...
C...  TAKE THE NEXT STEP IN X
1     CONTINUE
C...
C...  INTEGRATION IS COMPLETE AT Y
      RETURN
      END
