1 2 SUBROUTINE PLUME_VOL(XARG,ISDX,BVARG,BVARG3,BHARG,VOLOUT) 3 !*********************************************************************** 4 ! PLUME_VOL Module of the AMS/EPA Regulatory Model - AERMOD 5 ! 6 ! PURPOSE: Calculates plume volume for PVMRM option 7 ! 8 ! PROGRAMMER: Roger W. Brode, PES, Inc. 9 ! 10 ! DATE: May 6, 2002 11 ! 12 ! INPUTS: 13 ! 14 ! 15 ! OUTPUTS: 16 ! 17 ! 18 ! CALLED FROM: PVMRM_CALC 19 !*********************************************************************** 20 21 ! Variable Declarations 22 USE MAIN1 23 IMPLICIT NONE 24 CHARACTER MODNAM*12 25 INTEGER , PARAMETER :: NTAB = 79 26 REAL , PARAMETER :: NSUBZ = 4.00 , SRMIN = 5.0 27 REAL :: XARG , BVARG , BVARG3 , BHARG , VERT , VOLOUT , SUM , SUM3 28 REAL :: XTABLE(NTAB) , SR(NTAB) , SR3(NTAB) , SIGR , SIGR3 29 REAL :: XTAB , DELTAX , BVTMP , SIGRZTMP 30 INTEGER :: I , ISDX 31 INTEGER :: KITER , NDXHE , NDXZMID , NDXZPL , NDXHE3 32 INTEGER :: IPOSITN , NDXBH , N1 , N2 , N , IS 33 REAL :: HSPRIM , ZPLM , DHFOLD , SVPM , UPM , TGPM , PTPM , PTP , & 34 & ZMID 35 36 SAVE 37 38 ! Variable Initializations 39 DATA XTABLE/10. , 20. , 30. , 40. , 50. , 60. , 70. , 80. , 90. , & 40 & 100. , 110. , 120. , 130. , 140. , 150. , 160. , 170. , & 41 & 180. , 190. , 200. , 210. , 220. , 230. , 240. , 250. , & 42 & 260. , 300. , 400. , 500. , 600. , 700. , 800. , 900. , & 43 & 1000. , 1100. , 1200. , 1300. , 1400. , 1500. , 1600. , & 44 & 1700. , 1800. , 1900. , 2000. , 2100. , 2200. , 2300. , & 45 & 2400. , 2500. , 2600. , 2700. , 2800. , 2900. , 3000. , & 46 & 3100. , 3300. , 3500. , 3700. , 4000. , 5000. , 5500. , & 47 & 7000. , 8000. , 9000. , 10000. , 12000. , 14000. , 16000. , & 48 & 18000. , 20000. , 23000. , 25000. , 27500. , 30000. , & 49 & 35000. , 40000. , 45000. , 50000. , 500000./ 50 DATA SR/NTAB*0.0/ , SR3/NTAB*0.0/ 51 52 MODNAM = 'PLUME_VOL' ! 0 53 54 ! Assign source index to global variable 55 ISRC = ISDX 56 SUM = 0.0 57 SUM3 = 0.0 58 59 ! Set Mixing Height and Profiles for Urban Option if Needed 60 IF ( STABLE .AND. URBAN ) THEN 61 IF ( URBSRC(ISRC).EQ.'Y' ) THEN ! 0 62 URBSTAB = .TRUE. ! 0 63 ZI = AMAX1(ZIURB,ZIMECH) 64 GRIDSV = GRDSVU 65 GRIDSW = GRDSWU 66 GRIDTG = GRDTGU 67 GRIDPT = GRDPTU 68 ELSEIF ( URBSRC(ISRC).EQ.'N' ) THEN 69 URBSTAB = .FALSE. ! 0 70 ZI = ZIRUR 71 GRIDSV = GRDSVR 72 GRIDSW = GRDSWR 73 GRIDTG = GRDTGR 74 GRIDPT = GRDPTR 75 ENDIF 76 ELSE 77 URBSTAB = .FALSE. ! 0 78 ENDIF 79 80 ! Set the Source Variables for This Source --- CALL SETSRC 81 CALL SETSRC ! 0 82 83 ! Set Deposition Variables for this Source 84 ! Calculate Deposition Velocities for this Source --- CALL VDP 85 IF ( LDPART .OR. LWPART .OR. LDGAS ) CALL VDP 86 87 ! Calculate Scavenging Ratios for this Source --- CALL SCAVRAT 88 IF ( LWPART .OR. LWGAS ) CALL SCAVRAT 89 90 ! Calculate the initial meteorological variables --- CALL METINI 91 CALL METINI 92 93 IF ( SRCTYP(ISRC).EQ.'VOLUME' .OR. SRCTYP(ISRC)(1:4).EQ.'AREA' ) & 94 & THEN 95 FB = 0.0 ! 0 96 FM = 0.0 97 PPF = 0.0 98 HSP = HS 99 DHP = 0.0 100 DHP1 = 0.0 101 DHP2 = 0.0 102 DHP3 = 0.0 103 DHCRIT = 0.0 104 XFINAL = 0.0 105 XMIXED = ZI*UAVG/SWAVG 106 IF ( XMIXED.LT.XFINAL ) XMIXED = XFINAL 107 ZMIDMX = 0.5*ZI 108 109 ELSEIF ( SRCTYP(ISRC).EQ.'POINT' ) THEN 110 ! Calculate Buoyancy and Momentum Fluxes --- CALL FLUXES 111 CALL FLUXES ! 0 112 113 ! Set Wake and Building Type Switches --- CALL WAKFLG 114 ! --- NOTE: WAKFLG sets building dimensions based on wind 115 ! direction at stack top. 116 WAKE = .FALSE. 117 118 ! Define temporary values of CENTER and SURFAC based on HS 119 CENTER = HS 120 IF ( CENTER.LT.0.1*ZI ) THEN 121 SURFAC = .TRUE. ! 0 122 ELSE 123 SURFAC = .FALSE. ! 0 124 ENDIF 125 126 ! Check for stack-tip downwash option and adjust if necessary 127 IF ( NOSTD ) THEN ! 0 128 ! No stack-tip downwash, no adjustments necessary 129 HSP = HS ! 0 130 ELSE 131 ! Make adjustments for stack-tip downwash 132 HSP = HSPRIM(US,VS,HS,DS) ! 0 133 ENDIF 134 135 ! Calculate Distance to Final Rise --- CALL DISTF 136 CALL DISTF ! 0 137 138 IF ( DEBUG ) THEN 139 WRITE (DBGUNT,6000) DHFAER , UP , TGS ! 0 140 6000 FORMAT (/,5X,'INITIAL PLUME RISE ESTIMATE: DELH = ',F6.1, & 141 & ' M; Uplume = ',F5.2,' M/S; DTHDZ = ',F7.4, & 142 & ' DEG K/M') 143 ENDIF 144 145 IF ( STABLE .OR. (UNSTAB .AND. (HS.GE.ZI)) ) THEN ! 0 146 ! Use iterative approach to stable plume rise calculations 147 KITER = 0 ! 0 148 50 ZPLM = HSP + 0.5*DHFAER ! 0 149 DHFOLD = DHFAER 150 151 !---- Locate index below ZPLM 152 153 CALL LOCATE(GRIDHT,1,MXGLVL,ZPLM,NDXZPL) 154 155 !---- Get Wind speed at ZPLM; replace UP. Also, replace TGP, 156 ! vertical potential temperature gradient, if stable. 157 158 CALL GINTRP(GRIDHT(NDXZPL),GRIDSV(NDXZPL),GRIDHT(NDXZPL+1), & 159 & GRIDSV(NDXZPL+1),ZPLM,SVPM) 160 CALL GINTRP(GRIDHT(NDXZPL),GRIDWS(NDXZPL),GRIDHT(NDXZPL+1), & 161 & GRIDWS(NDXZPL+1),ZPLM,UPM) 162 SVPM = AMAX1(SVPM,SVMIN,0.05*UPM) 163 UPM = SQRT(UPM*UPM+2.*SVPM*SVPM) 164 !RWB Use average of stack top and midpoint wind speeds. 165 UP = 0.5*(US+UPM) 166 167 CALL GINTRP(GRIDHT(NDXZPL),GRIDTG(NDXZPL),GRIDHT(NDXZPL+1), & 168 & GRIDTG(NDXZPL+1),ZPLM,TGPM) 169 CALL GINTRP(GRIDHT(NDXZPL),GRIDPT(NDXZPL),GRIDHT(NDXZPL+1), & 170 & GRIDPT(NDXZPL+1),ZPLM,PTPM) 171 !RWB Use average of stack top and midpoint temperature gradients. 172 TGP = 0.5*(TGS+TGPM) 173 PTP = 0.5*(PTS+PTPM) 174 BVF = SQRT(G*TGP/PTP) 175 IF ( BVF.LT.1.0E-10 ) BVF = 1.0E-10 176 BVPRIM = 0.7*BVF 177 178 CALL DISTF 179 180 KITER = KITER + 1 181 182 !RJP Add temporary debugging statements 183 184 IF ( DEBUG ) THEN 185 WRITE (DBGUNT,6001) KITER , DHFOLD , DHFAER , ZPLM , UP ,& 186 & TGP 187 6001 FORMAT (/,5X,'OPTH2 ITER. #',I1,': OLD DELH = ',F6.1, & 188 & ' M; NEW DELH = ',F6.1,' M; MET LEVEL = ',F6.1, & 189 & ' M; NEW Upl = ',F5.2,' M/S; NEW DTHDZ = ',F7.4, & 190 & ' K/M') 191 ENDIF 192 IF ( ABS((DHFOLD-DHFAER)/DHFAER).LT.0.01 ) GOTO 60 ! 0 193 IF ( KITER.GE.5 ) THEN ! 0 194 DHFAER = 0.5*(DHFAER+DHFOLD) ! 0 195 IF ( DEBUG ) WRITE (DBGUNT,6002) DHFAER 196 6002 FORMAT (/,5X,'OPTH2 ITERATION FAILED TO CONVERGE; PLUME',& 197 & ' RISE SET AT ',F6.1,' METERS.',/) 198 GOTO 60 199 ELSE 200 GOTO 50 ! 0 201 ENDIF 202 203 60 CONTINUE ! 0 204 205 !RWB After completing iteration, reset UP and TGP to stack top 206 !RWB values for subsequent distance-dependent plume rise calcs. 207 UP = US 208 TGP = TGS 209 PTP = PTS 210 BVF = SQRT(G*TGP/PTP) 211 IF ( BVF.LT.1.0E-10 ) BVF = 1.0E-10 212 BVPRIM = 0.7*BVF 213 ENDIF 214 215 ! Initialize FSTREC Logical Switch for First Receptor of Loop; 216 FSTREC = .TRUE. ! 0 217 PRM_FSTREC = .TRUE. 218 219 ZMIDMX = 0.5*ZI 220 221 !RJP 222 !RJP Calculate distance to uniformly mixed plume within the 223 !RJP boundary layer (XMIXED) after Turner's Workbook (1970), page 7: 224 !RJP distance is approximately (Zi * UAVG)/SWAVG, where UAVG 225 !RJP and SWAVG are wind speed and sigma-w averaged over the depth 226 !RJP between the ground and Zi (or the plume height, if higher in 227 !RJP stable conditions); this height is denoted as 2 * ZMIDMX. 228 !RJP 229 !RJP First, get refined estimate of final rise and distance to final 230 !RJP rise if downwash conditions prevail. 231 !RJP 232 XFINAL = XMAX 233 DHCRIT = DHFAER 234 XMIXED = ZI*UAVG/SWAVG 235 IF ( UNSTAB .AND. HS.LT.ZI ) THEN 236 ! Check for XMIXED smaller than 1.25*XFINAL 237 IF ( XMIXED.LT.1.25*XFINAL ) THEN ! 0 238 XFINAL = 0.8*XMIXED ! 0 239 CALL CBLPRD(XFINAL) 240 DHCRIT = DHP1 241 ENDIF 242 ENDIF 243 244 ENDIF 245 246 ! First build table of relative dispersion coefficients for 247 ! dominant source (source index = ISDX) 248 DO I = 1 , NTAB ! 0 249 250 XTAB = XTABLE(I) ! 0 251 IF ( I.GT.1 ) THEN 252 DELTAX = XTABLE(I) - XTABLE(I-1) ! 0 253 ELSE 254 DELTAX = XTABLE(1) ! 0 255 ENDIF 256 257 IF ( XARG.GT.XTAB ) THEN ! 0 258 259 ! Define plume centroid height (CENTER) for use in 260 ! inhomogeneity calculations 261 CALL CENTROID(XTAB) ! 0 262 263 ! Calculate the plume rise --- CALL DELTAH 264 IF ( SRCTYP(ISRC).EQ.'POINT' ) CALL DELTAH(XTAB) 265 266 ! If the atmosphere is unstable and the stack 267 ! top is below the mixing height, calculate 268 ! the CBL PDF coefficients --- CALL PDF 269 IF ( UNSTAB .AND. (HS.LT.ZI) ) CALL PDF 270 271 ! Determine Effective Plume Height --- CALL HEFF 272 CALL HEFF(XTAB) 273 274 ! Compute effective parameters using an 275 ! average through plume rise layer 276 CALL IBLVAL(XTAB) 277 278 ! Call PDF & HEFF again for final CBL plume heights 279 IF ( UNSTAB .AND. (HS.LT.ZI) ) THEN 280 CALL PDF ! 0 281 CALL HEFF(XTAB) 282 ENDIF 283 284 IF ( SRCTYP(ISRC).EQ.'POINT' ) THEN ! 0 285 ! Call BID to get buoyancy-induced dispersion terms 286 CALL BID ! 0 287 ELSE 288 SZB = 0.0 ! 0 289 SZBD = 0.0 290 SZB3 = 0.0 291 ENDIF 292 293 ! Determine Relative Dispersion Parameters --- CALL RELDISP 294 CALL RELDISP(XTAB,SR(I),SR3(I)) ! 0 295 296 IF ( I.GT.1 ) THEN 297 SIGR = MAX(SRMIN,0.5*(SR(I)+SR(I-1))) ! 0 298 ELSE 299 SIGR = MAX(SRMIN,0.5*SR(I)) ! 0 300 ENDIF 301 302 ! Calculate vertical dimension (VERT) taking into account ZI limit 303 VERT = BVARG + 2.*NSUBZ*SIGR ! 0 304 IF ( UNSTAB .AND. VERT.GT.ZI ) THEN 305 VERT = ZI ! 0 306 BVTMP = MAX(0.0,ZI-2.*NSUBZ*SIGR) 307 IF ( 2.*NSUBZ*SIGR.GT.ZI ) THEN 308 SIGRZTMP = ZI/(2.*NSUBZ) ! 0 309 ELSE 310 SIGRZTMP = SIGR ! 0 311 ENDIF 312 ELSE 313 BVTMP = BVARG ! 0 314 SIGRZTMP = SIGR 315 ENDIF 316 317 318 ! Plume volume calculation based on rectangle with rounded corners 319 ! First component is for major contributing plume only 320 SUM = SUM + (PI*(NSUBZ*SIGR)*(NSUBZ*SIGRZTMP)+VERT*BHARG+ & 321 & 2.*NSUBZ*SIGR*BVTMP)*DELTAX 322 323 IF ( UNSTAB .AND. PPFACT(ISRC).GT.0.0 ) THEN 324 325 IF ( I.GT.1 ) THEN ! 0 326 SIGR3 = MAX(SRMIN,0.5*(SR3(I)+SR3(I-1))) ! 0 327 ELSE 328 SIGR3 = MAX(SRMIN,0.5*SR3(I)) ! 0 329 ENDIF 330 331 VERT = BVARG3 + 2.*NSUBZ*SIGR3 ! 0 332 333 SUM3 = SUM3 + (PI*(NSUBZ*SIGR3)*(NSUBZ*SIGR3)+VERT*BHARG+& 334 & 2.*NSUBZ*SIGR3*BVARG3)*DELTAX 335 336 ENDIF 337 338 ELSE 339 340 ! Define plume centroid height (CENTER) for use in 341 ! inhomogeneity calculations 342 CALL CENTROID(XARG) ! 0 343 344 ! Calculate the plume rise --- CALL DELTAH 345 IF ( SRCTYP(ISRC).EQ.'POINT' ) CALL DELTAH(XARG) 346 347 ! If the atmosphere is unstable and the stack 348 ! top is below the mixing height, calculate 349 ! the CBL PDF coefficients --- CALL PDF 350 IF ( UNSTAB .AND. (HS.LT.ZI) ) CALL PDF 351 352 ! Determine Effective Plume Height --- CALL HEFF 353 CALL HEFF(XARG) 354 355 ! Compute effective parameters using an 356 ! average through plume rise layer 357 CALL IBLVAL(XARG) 358 359 ! Call PDF & HEFF again for final CBL plume heights 360 IF ( UNSTAB .AND. (HS.LT.ZI) ) THEN 361 CALL PDF ! 0 362 CALL HEFF(XARG) 363 ENDIF 364 365 IF ( SRCTYP(ISRC).EQ.'POINT' ) THEN ! 0 366 ! Call BID to get buoyancy-induced dispersion terms 367 CALL BID ! 0 368 ELSE 369 SZB = 0.0 ! 0 370 SZBD = 0.0 371 SZB3 = 0.0 372 ENDIF 373 374 ! Determine Relative Dispersion Parameters --- CALL RELDISP 375 CALL RELDISP(XARG,SR(I),SR3(I)) ! 0 376 377 DELTAX = XARG - XTABLE(I-1) 378 379 IF ( I.GT.1 ) THEN 380 SIGR = MAX(SRMIN,0.5*(SR(I)+SR(I-1))) ! 0 381 ELSE 382 SIGR = MAX(SRMIN,0.5*SR(I)) ! 0 383 ENDIF 384 385 ! Calculate vertical dimension (VERT) taking into account ZI limit 386 VERT = BVARG + 2.*NSUBZ*SIGR ! 0 387 IF ( UNSTAB .AND. VERT.GT.ZI ) THEN 388 VERT = ZI ! 0 389 BVTMP = MAX(0.0,ZI-2.*NSUBZ*SIGR) 390 IF ( 2.*NSUBZ*SIGR.GT.ZI ) THEN 391 SIGRZTMP = ZI/(2.*NSUBZ) ! 0 392 ELSE 393 SIGRZTMP = SIGR ! 0 394 ENDIF 395 ELSE 396 BVTMP = BVARG ! 0 397 SIGRZTMP = SIGR 398 ENDIF 399 400 ! Plume volume calculation based on rectangle with rounded corners 401 ! First component is for major contributing plume only 402 SUM = SUM + (PI*(NSUBZ*SIGR)*(NSUBZ*SIGRZTMP)+VERT*BHARG+ & 403 & 2.*NSUBZ*SIGR*BVTMP)*DELTAX 404 405 IF ( UNSTAB .AND. PPFACT(ISRC).GT.0.0 ) THEN 406 407 IF ( I.GT.1 ) THEN ! 0 408 SIGR3 = MAX(SRMIN,0.5*(SR3(I)+SR3(I-1))) ! 0 409 ELSE 410 SIGR3 = MAX(SRMIN,0.5*SR3(I)) ! 0 411 ENDIF 412 413 VERT = BVARG3 + 2.*NSUBZ*SIGR3 ! 0 414 415 SUM3 = SUM3 + (PI*(NSUBZ*SIGR3)*(NSUBZ*SIGR3)+VERT*BHARG+& 416 & 2.*NSUBZ*SIGR3*BVARG3)*DELTAX 417 418 ENDIF 419 420 GOTO 100 ! 0 421 422 ENDIF 423 424 ENDDO 425 426 ! Combine volume for "direct" plume volume and penetrated plume volume 427 100 VOLOUT = SUM*(1.-PPFACT(ISRC)) + SUM3*PPFACT(ISRC) ! 0 428 429 CONTINUE 430 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