1 2 !----------------------------------------------------------------------- 3 SUBROUTINE QATR3(XL,XU,EPS,EPST,NDIM,IMIN,Y,YD,IER,I,AUX) 4 !----------------------------------------------------------------------- 5 ! 6 ! --- ISCST3 QATR3 7 ! 8 ! PURPOSE: Integration routine adapted from the IBM SSP program 9 ! DQATR. Modified for single precision. This is a COPY 10 ! of QATR for use in area source integrations. 11 ! 12 ! ARGUMENTS: 13 ! PASSED: xl,xu lower and upper limits of integration [r] 14 ! eps fractional error used to define convergence [r] 15 ! epst lower theshold check for value of integral [r] 16 ! ndim dimension of array aux (max no. of steps) [i] 17 ! imin minimum number of "steps" for integral [i] 18 ! fct external function (integrand) [r] 19 ! aux working array, passed to allow variable dim. [r] 20 ! RETURNED: y value of integral [r] 21 ! yd value of integral without wet depletion [r] 22 ! ier status flag at termination [i] 23 ! i number of subdivision steps [i] 24 ! 25 ! CALLING ROUTINES: PSIDE, PSIDE2 26 ! 27 ! EXTERNAL ROUTINES: none 28 !----------------------------------------------------------------------- 29 30 ! NOTES: status flags denote the following -- 31 ! ier=0 value of integral converged to within eps 32 ! ier=1 value of integral is diverging (not used) 33 ! ier=2 value of integral did not converge to within 34 ! eps before ndim limit was reached 35 36 IMPLICIT NONE 37 38 INTEGER :: NDIM , IMIN , IER 39 REAL :: Y , YD , EPS , EPST , XU , XL , H , HH , DELT2 , P , & 40 & DELT1 , HD , X , SM , SMD , Q , AUX(NDIM) , AUXD(NDIM) , & 41 & P1 , P2 , P1D , P2D 42 INTEGER :: I , II , JI , J , JJ 43 44 DO I = 1 , NDIM !1037136 45 AUX(I) = 0.0 ! 12445K 46 AUXD(I) = 0.0 47 ENDDO 48 49 !--- Preparations for Romberg loop 50 CALL PLUMEF(XL,P1,P1D) !1037136 51 CALL PLUMEF(XU,P2,P2D) 52 AUX(1) = 0.5*(P1+P2) 53 AUXD(1) = 0.5*(P1D+P2D) 54 H = XU - XL 55 56 IF ( H.EQ.0.0 .OR. AUX(1).EQ.0.0 ) THEN 57 IER = 0 ! 3544 58 Y = 0.0 59 YD = 0.0 60 RETURN 61 ENDIF 62 63 HH = H !1033592 64 DELT2 = 0. 65 P = 1. 66 JJ = 1 67 68 DO I = 2 , NDIM 69 Y = AUX(1) !4135462 70 DELT1 = DELT2 71 HD = HH 72 HH = 0.5*HH 73 P = 0.5*P 74 X = XL + HH 75 SM = 0. 76 SMD = 0. 77 78 DO J = 1 , JJ 79 CALL PLUMEF(X,P1,P1D) ! 15521K 80 SM = SM + P1 81 SMD = SMD + P1D 82 X = X + HD 83 ENDDO 84 85 !---- A new approximation to the integral is computed by means 86 ! of the trapezoidal rule 87 AUX(I) = 0.5*AUX(I-1) + P*SM !4135462 88 AUXD(I) = 0.5*AUXD(I-1) + P*SMD 89 90 !---- Start of Rombergs extrapolation method 91 92 Q = 1. 93 JI = I - 1 94 DO J = 1 , JI 95 II = I - J ! 10341K 96 Q = Q + Q 97 Q = Q + Q 98 AUX(II) = AUX(II+1) + (AUX(II+1)-AUX(II))/(Q-1.) 99 AUXD(II) = AUXD(II+1) + (AUXD(II+1)-AUXD(II))/(Q-1.) 100 ENDDO 101 102 !---- End of Romberg step 103 104 ! Compute absolute error, delt2 105 DELT2 = ABS(Y-AUX(1)) !4135462 106 107 IF ( I.GE.IMIN ) THEN 108 ! Check for covergence of algorithm 109 IF ( ABS(AUX(1)).LT.EPST ) THEN !1034686 110 ! Lower threshold convergence test 111 IER = 0 ! 89414 112 Y = H*AUX(1) 113 YD = H*AUXD(1) 114 RETURN 115 ELSEIF ( DELT2.LE.EPS*ABS(AUX(1)) ) THEN 116 ! Relative error convergence test 117 IER = 0 ! 942796 118 Y = H*AUX(1) 119 YD = H*AUXD(1) 120 RETURN 121 ELSEIF ( ABS(HH).LT.1.0 ) THEN 122 ! Minimum "delta-x" convergence test; < 1.0m 123 IER = 0 ! 1382 124 Y = H*AUX(1) 125 YD = H*AUXD(1) 126 RETURN 127 ENDIF 128 ENDIF 129 130 JJ = JJ + JJ !3101870 131 ENDDO 132 133 ! Convergence not reached within maximum number of steps 134 IER = 2 ! 0 135 Y = H*AUX(1) 136 YD = H*AUXD(1) 137 138 CONTINUE 139 END
HyperKWIC - Version 1.00DD executed at 20:00 on 1 Mar 2018 | Personal or Academic or Evaluation User | Free for Non-Commercial, Non-Government Use