1 2 SUBROUTINE PSIDE_TOX(U1,U2,DVAL,DVALD) 3 !*********************************************************************** 4 ! PSIDE_TOX Module of the AMS/EPA Regulatory Model - AERMOD 5 ! 6 ! Special version of PSIDE optimized for TOXICS applications. 7 ! Utilizes Romberg Integration (QATR3) or Gaussian Quadrature (QG2) 8 ! depending on the source receptor geometry. 9 ! 10 ! PURPOSE: INTEGRATES SIDE K of POLYGON 11 ! int f(u)*CNF(v(u)/sig(u))=f(u)*vn(u) from u1 to u2 12 ! CNF = cumulative normal distribution 13 ! Computes W(1),W(2)--normalized plume width at u1 u2 14 ! Checks for w(i) outside of -3.9,3.9 with i+, i- 15 ! L=-3.9 U=3.9 = bounds for testing 16 ! Integrates according to case encountered: 17 ! situation CASE iplus iminus integral limits 18 ! L<w1,w2<U 1 0 0 u1,u2 19 ! w1,w2<L 2 0 1+2 don't compute 20 ! w1,w2>U 3 1+2 0 u1,u2 21 ! w1<L<w2<U 4 0 1 u-,u2 22 ! w2<L<w1<U 5 0 2 u1,u- 23 ! L<w1<U<w2 6 2 0 u1,u+ u+,u2 24 ! L<w2<U<w1 7 1 0 u1,u+ u+,u2 25 ! w1<L<U<w2 8 2 1 u-,u+ u+,u2 26 ! w2<L<U<w1 9 1 2 u1,u+ u+,u- 27 ! 28 ! u+ = value of u such that w(u)=U=3.9 29 ! u- = " w(u)=L=-3.9 30 ! u+,u- computed with Brent's Algorithm 31 ! 32 ! IF uplus >0, part of side is outside plume 33 ! but is integrated anyway, unless there is 34 ! a corresponding part on another side that will 35 ! cause cancellation. This is determined in 36 ! PSIDE2; 37 ! 38 ! PROGRAMMER: Jeff Wang, Roger Brode 39 ! Adapted From Codes By Richard Strelitz, CSC 40 ! 41 ! DATE: July 7, 1993 42 ! 43 ! MODIFIED by Roger Brode, PES, Inc. to use QATR routine for 44 ! integration, and to pass VN(1) and VN(2) to 45 ! ZBRENT. Also changed to use -3.9 to +3.9 for 46 ! width of plume for internal consistency with 47 ! SUBROUTINE PWIDTH. - 12/14/98 48 ! 49 ! MODIFIED by Roger Brode, PES, Inc. to correct lower integration 50 ! limit for Case 4, and to remove extraneous calls 51 ! to XWIDTH and PWIDTH after calls to ZBRENT. - 7/29/94 52 ! 53 ! INPUTS: End Points of The Segments 54 ! 55 ! OUTPUTS: Integral Value (if any) for Segment 56 ! 57 ! CALLED FROM: AREAIN 58 !*********************************************************************** 59 60 ! Variable Declarations 61 USE MAIN1 62 IMPLICIT NONE 63 CHARACTER MODNAM*12 64 65 SAVE 66 67 !---- Set convergence criteria for calls to QATR3: 68 ! NDIM = maximum number of integration steps 69 ! IMIN = minimum number of integration steps 70 ! EPSR = relative error tolerance for integral 71 ! EPST = minimum value threshold for integral 72 !---- 73 INTEGER , PARAMETER :: NDIM = 10 , IMIN = 4 74 REAL , PARAMETER :: EPSR = 2.0E-2 , EPST = 1.0E-5 75 76 !---- Set distance factor for switching to Gaussian Quadrature, QG_FACT 77 ! If Xmax - Xmin is .LT. QG_FACT*Xmin, then use QG2, where 78 ! Xmax and Xmin are the distances to the endpoints of the side. 79 REAL , PARAMETER :: QG_FACT = 5.0 80 81 INTEGER :: I , KS , IMINUS , IPLUS , NOUT , ICON 82 REAL :: DVAL , DVALD , U1 , U2 , UMINUS , UPLUS , AUX(NDIM) , & 83 & U(2) , V1(2) , VN(2) , W(2) 84 85 ! Variable Initializations 86 MODNAM = 'PSIDE_TOX' ! 0 87 88 ! NSEG = number of segments; set to 0 in AREAIN 89 ! for each source/rcvr/time step 90 DVAL = 0.0 91 DVALD = 0.0 92 93 DO I = 1 , 2 94 KS = IVERT + I - 1 ! 0 95 U(I) = UVERT(KS) 96 V1(I) = VVERT(KS) 97 VN(I) = VNVERT(KS) 98 W(I) = WVERT(KS) 99 ENDDO 100 101 IMINUS = 0 ! 0 102 IPLUS = 0 103 UMINUS = -1.0 104 UPLUS = -1.0 105 DO I = 1 , 2 106 IF ( VN(I).LT.-3.9 ) IMINUS = I + IMINUS ! 0 107 IF ( VN(I).GT.3.9 ) IPLUS = I + IPLUS 108 ENDDO 109 110 IF ( IPLUS.EQ.1 .OR. IPLUS.EQ.2 ) & 111 & CALL ZBRENT(1,U(1),U(2),VN(1),VN(2),1.0,UPLUS) 112 IF ( IMINUS.EQ.1 .OR. IMINUS.EQ.2 ) & 113 & CALL ZBRENT(-1,U(1),U(2),VN(1),VN(2),1.0,UMINUS) 114 115 !---- CASE DEPENDs on iplus, iminus 116 IF ( IPLUS.EQ.0 ) THEN 117 IF ( IMINUS.EQ.0 ) THEN ! 0 118 ! iplus iminus case 119 ! 0 0 1 120 IF ( ABS(U2-U1).LT.QG_FACT*MIN(U1,U2) ) THEN ! 0 121 CALL QG2(U1,U2,DVAL,DVALD) ! 0 122 ELSE 123 CALL QATR3(U1,U2,EPSR,EPST,NDIM,IMIN,DVAL,DVALD,ICON, & 124 & NOUT,AUX) 125 ENDIF 126 ELSEIF ( IMINUS.EQ.3 ) THEN 127 ! 0 3 2 128 DVAL = 0 ! 0 129 ELSEIF ( IMINUS.EQ.1 ) THEN 130 ! 0 1 4 131 IF ( ABS(U2-UMINUS).LT.QG_FACT*MIN(UMINUS,U2) ) THEN ! 0 132 CALL QG2(UMINUS,U2,DVAL,DVALD) ! 0 133 ELSE 134 CALL QATR3(UMINUS,U2,EPSR,EPST,NDIM,IMIN,DVAL,DVALD,ICON,& 135 & NOUT,AUX) 136 ENDIF 137 ELSE 138 ! 0 2 5 139 IF ( ABS(UMINUS-U1).LT.QG_FACT*MIN(U1,UMINUS) ) THEN ! 0 140 CALL QG2(U1,UMINUS,DVAL,DVALD) ! 0 141 ELSE 142 CALL QATR3(U1,UMINUS,EPSR,EPST,NDIM,IMIN,DVAL,DVALD,ICON,& 143 & NOUT,AUX) 144 ENDIF 145 ENDIF 146 147 ELSEIF ( IPLUS.EQ.1 ) THEN 148 NSEGS = NSEGS + 1 ! 0 149 UASEGS(NSEGS) = U1 150 UBSEGS(NSEGS) = UPLUS 151 IF ( IMINUS.EQ.0 ) THEN 152 ! 1 0 7 153 IF ( ABS(U2-UPLUS).LT.QG_FACT*MIN(UPLUS,U2) ) THEN ! 0 154 CALL QG2(UPLUS,U2,DVAL,DVALD) ! 0 155 ELSE 156 CALL QATR3(UPLUS,U2,EPSR,EPST,NDIM,IMIN,DVAL,DVALD,ICON, & 157 & NOUT,AUX) 158 ENDIF 159 ELSE 160 ! 1 2 9 161 IF ( ABS(UMINUS-UPLUS).LT.QG_FACT*MIN(UPLUS,UMINUS) ) THEN ! 0 162 CALL QG2(UPLUS,UMINUS,DVAL,DVALD) ! 0 163 ELSE 164 CALL QATR3(UPLUS,UMINUS,EPSR,EPST,NDIM,IMIN,DVAL,DVALD, & 165 & ICON,NOUT,AUX) 166 ENDIF 167 ENDIF 168 169 ELSEIF ( IPLUS.EQ.2 ) THEN 170 NSEGS = NSEGS + 1 ! 0 171 UASEGS(NSEGS) = UPLUS 172 UBSEGS(NSEGS) = U2 173 IF ( IMINUS.EQ.0 ) THEN 174 ! 2 0 6 175 IF ( ABS(UPLUS-U1).LT.QG_FACT*MIN(U1,UPLUS) ) THEN ! 0 176 CALL QG2(U1,UPLUS,DVAL,DVALD) ! 0 177 ELSE 178 CALL QATR3(U1,UPLUS,EPSR,EPST,NDIM,IMIN,DVAL,DVALD,ICON, & 179 & NOUT,AUX) 180 ENDIF 181 ELSE 182 ! 2 1 8 183 IF ( ABS(UPLUS-UMINUS).LT.QG_FACT*MIN(UMINUS,UPLUS) ) THEN ! 0 184 CALL QG2(UMINUS,UPLUS,DVAL,DVALD) ! 0 185 ELSE 186 CALL QATR3(UMINUS,UPLUS,EPSR,EPST,NDIM,IMIN,DVAL,DVALD, & 187 & ICON,NOUT,AUX) 188 ENDIF 189 ENDIF 190 191 ELSE 192 ! 3 0 3 193 NSEGS = NSEGS + 1 ! 0 194 UASEGS(NSEGS) = U1 195 UBSEGS(NSEGS) = U2 196 ENDIF 197 198 CONTINUE ! 0 199 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