1 2 3 SUBROUTINE ACALC 4 !*********************************************************************** 5 ! ACALC Module of the AMS/EPA Regulatory Model - AERMOD 6 ! 7 ! PURPOSE: Calculates concentration or deposition values 8 ! for AREA sources utilizing an integrated line source. 9 ! 10 ! PROGRAMMER: Jeff Wang, Roger Brode 11 ! 12 ! DATE: March 2, 1992 13 ! 14 ! MODIFIED: Modified to include initialization of __VAL arrays 15 ! at end of receptor loop. 16 ! R. W. Brode, MACTEC (f/k/a PES), Inc., 10/26/04 17 ! 18 ! MODIFIED: To include the PVMRM and OLM options for 19 ! modeling conversion of NOx to NO2. 20 ! Added debug statement based on ENSR code. 21 ! R. W. Brode, MACTEC (f/k/a PES), Inc., 07/27/04 22 ! 23 ! MODIFIED: To include tilted plume for point source 24 ! approximation of particle emissions. 25 ! R. W. Brode, MACTEC (f/k/a PES), Inc., 07/23/04 26 ! 27 ! MODIFIED: To allow TOXICS option for AREAPOLY sources. 28 ! R. W. Brode, MACTEC (f/k/a PES), Inc., 05/12/04 29 ! 30 ! MODIFIED: To assign value to XDIST for use in dry depletion. 31 ! R. W. Brode, MACTEC (f/k/a PES), Inc., 03/19/04 32 ! 33 ! MODIFIED: To avoid potential math errors for AREAPOLY sources 34 ! R.W. Brode, PES, Inc. - 02/25/02 35 ! 36 ! MODIFIED: To incorporate numerical integration algorithm 37 ! for AREA source - 7/7/93 38 ! 39 ! INPUTS: Source Parameters for Specific Source 40 ! Arrays of Receptor Locations 41 ! Meteorological Variables for One Hour 42 ! 43 ! OUTPUTS: Array of 1-hr CONC or DEPOS Values for Each Source/Receptor 44 ! 45 ! CALLED FROM: CALC 46 !*********************************************************************** 47 48 ! Variable Declarations 49 USE MAIN1 50 IMPLICIT NONE 51 CHARACTER MODNAM*12 52 INTEGER :: I , J 53 REAL :: XDEP , WIDTH , LENGTH , XMAXR , QTKSAV , XPOINT , ADJ 54 REAL :: AEROUT(NUMTYP) , FYOUT 55 56 SAVE 57 REAL XSPA(NVMAX) , YSPA(NVMAX) 58 59 ! Variable Initializations 60 MODNAM = 'ACALC' ! 9132 61 WAKE = .FALSE. 62 63 DO ITYP = 1 , NUMTYP 64 HRVAL(ITYP) = 0.0 ! 9132 65 HRVALD(ITYP) = 0.0 66 AERVAL(ITYP) = 0.0 67 AERVALD(ITYP) = 0.0 68 AEROUT(ITYP) = 0.0 69 ENDDO 70 71 ! Set the Source Variables for This Source --- CALL SETSRC 72 CALL SETSRC ! 9132 73 74 ! Apply Variable Emission Rate and Unit Factors --- CALL EMFACT 75 CALL EMFACT(QS) 76 77 ! Initialize 'ARC' Arrays for EVALFILE Output --- CALL EVLINI 78 IF ( EVAL(ISRC) ) CALL EVLINI 79 80 IF ( QTK.NE.0.0 ) THEN 81 ! Set Mixing Height and Profiles for Urban Option if Needed 82 IF ( STABLE .AND. URBAN ) THEN ! 9132 83 IF ( URBSRC(ISRC).EQ.'Y' ) THEN ! 0 84 URBSTAB = .TRUE. ! 0 85 ZI = AMAX1(ZIURB,ZIMECH) 86 GRIDSV = GRDSVU 87 GRIDSW = GRDSWU 88 GRIDTG = GRDTGU 89 GRIDPT = GRDPTU 90 OBULEN = ABS(URBOBULEN) 91 USTAR = URBUSTR 92 ELSEIF ( URBSRC(ISRC).EQ.'N' ) THEN 93 URBSTAB = .FALSE. ! 0 94 ZI = ZIRUR 95 GRIDSV = GRDSVR 96 GRIDSW = GRDSWR 97 GRIDTG = GRDTGR 98 GRIDPT = GRDPTR 99 OBULEN = RUROBULEN 100 USTAR = RURUSTR 101 ENDIF 102 ELSE 103 URBSTAB = .FALSE. ! 9132 104 ENDIF 105 106 ! Initialize meteorological variables --- CALL METINI 107 CALL METINI ! 9132 108 ! Initialize miscellaneous variables 109 FB = 0.0 110 FM = 0.0 111 PPF = 0.0 112 HSP = HS 113 DHP = 0.0 114 DHP1 = 0.0 115 DHP2 = 0.0 116 DHP3 = 0.0 117 DHCRIT = 0.0 118 XFINAL = 0.0 119 XMIXED = ZI*UAVG/SWAVG 120 IF ( XMIXED.LT.XFINAL ) XMIXED = XFINAL 121 ZMIDMX = 0.5*ZI 122 123 !DEP Initialize PDF parameters for use in calculating ZSUBP 124 IF ( UNSTAB .AND. (HS.LT.ZI) ) CALL PDF 125 ! Set Deposition Variables for this Source 126 ! Calculate Deposition Velocities for this Source --- CALL VDP 127 IF ( LDPART .OR. LDGAS ) CALL VDP 128 IF ( LWPART .OR. LWGAS ) THEN 129 !PES Set value of ZSUBP = MAX( ZI, TOP OF PLUME ), where 130 !PES TOP OF PLUME is defined as plume height (HE) plus 2.15*SZ, 131 !PES evaluated at a distance of 20 kilometers downwind. 132 !PES Apply minimum value of 500m and maximum value of 10,000m. 133 IF ( STABLE .OR. (UNSTAB .AND. HS.GE.ZI) ) THEN ! 0 134 CALL SIGZ(20000.) ! 0 135 ZSUBP = MAX(500.,ZI,HS+SZCOEF*SZAS) 136 ELSEIF ( UNSTAB ) THEN 137 CALL SIGZ(20000.) ! 0 138 ZSUBP = MAX(500.,ZI,HS+SZCOEF*(SZAD1+SZAD2)/2.) 139 ENDIF 140 ZSUBP = MIN(10000.,ZSUBP) ! 0 141 ! Calculate Scavenging Ratios for this Source --- CALL SCAVRAT 142 CALL SCAVRAT 143 ENDIF 144 145 ! Begin Receptor LOOP 146 RECEPTOR_LOOP:DO IREC = 1 , NUMREC ! 9132 147 ! Calculate Down and Crosswind Distances --- CALL ARDIST 148 IF ( EVONLY ) THEN !1315008 149 CALL ARDIST(IEVENT,XDEP,WIDTH,LENGTH,XMAXR) ! 0 150 ELSE 151 CALL ARDIST(IREC,XDEP,WIDTH,LENGTH,XMAXR) !1315008 152 ENDIF 153 154 ! Check to see if receptor is upwind of area source 155 IF ( XMAXR.LT.1.0 ) GOTO 50 !1315008 156 157 ! Check to see if receptor is more than 80km from source. 158 IF ( TOXICS .AND. DISTR.GT.80000. ) GOTO 50 ! 740558 159 160 ! Set initial effective parameters 161 UEFF = US ! 740558 162 SVEFF = SVS 163 SWEFF = SWS 164 TGEFF = TGS 165 IF ( UNSTAB .AND. (HS.LT.ZI) ) THEN 166 UEFFD = US ! 121272 167 SVEFFD = SVS 168 SWEFFD = SWS 169 UEFFN = US 170 SVEFFN = SVS 171 SWEFFN = SWS 172 ENDIF 173 174 ! Check to see if receptor is beyond edge of plume laterally. 175 IF ( (ABS(Y)-0.5*WIDTH).GT.0. ) THEN ! 740558 176 CALL ADISY(XMAXR) ! 712088 177 IF ( (ABS(Y)-0.5*WIDTH).GE.4.*SY ) GOTO 50 178 ENDIF 179 180 IF ( DEBUG ) THEN ! 271432 181 WRITE (DBGUNT,6015) UEFF , SVEFF , SWEFF ! 0 182 ! ENSR ENHANCEMENT OF WRITE STATEMENT 183 6015 FORMAT (//,'AERMOD AREA SOURCE COMPONENT',/,5X, & 184 & 'Initial effective parameters for ', & 185 & 'stable or direct convective ','plume:',//,5x, & 186 & 'Ueff = ',F7.2,' m/s; ','SVeff = ',F7.2, & 187 & ' m/s; SWeff = ',F7.2,' m/s.',/) 188 ENDIF 189 190 ! Determine the CRITical Dividing Streamline --- CALL CRITDS 191 CALL CRITDS(HE) ! 271432 192 193 ! Set distance factor for point source approx. for TOXICS based 194 ! on "equivalent" PG stability class (KST) 195 IF ( URBSTAB ) THEN 196 VP_FACT = VIRTPNT_URB(KST) ! 0 197 ELSE 198 VP_FACT = VIRTPNT_RUR(KST) ! 271432 199 ENDIF 200 201 ! Calculate distance for switch to point source approx. for TOXICS 202 XPOINT = 1.5*LENGTH + VP_FACT*WIDTH ! 271432 203 204 ! --- Assign XDIST for use in dry depletion (FUNCTION F2INT) 205 XDIST = X 206 207 IF ( .NOT.TOXICS .OR. (TOXICS .AND. X.LT.XPOINT) ) THEN 208 IF ( ARDPLETE ) THEN ! 271432 209 210 IF ( LDGAS .OR. LWGAS ) THEN ! 0 211 CALL PDEPG(XDEP) ! 0 212 ELSE 213 DQCORG = 1.0 ! 0 214 WQCORG = 1.0 215 ENDIF 216 IF ( LDPART .OR. LWPART ) THEN ! 0 217 CALL PDEP(XDEP) ! 0 218 ELSEIF ( NPD.GT.0 ) THEN 219 DO J = 1 , NPD ! 0 220 DQCOR(J) = 1.0 ! 0 221 WQCOR(J) = 1.0 222 ENDDO 223 ENDIF 224 225 ENDIF 226 227 DO ITYP = 1 , NUMTYP ! 271432 228 ! Calculate Area Source Integral --- CALL AREAIN 229 CALL AREAIN ! 271432 230 ENDDO 231 ELSE 232 ! Use point source approximation 233 ! Save emissions per unit area and calculate total emissions 234 QTKSAV = QTK ! 0 235 IF ( SRCTYP(ISRC).EQ.'AREA' .OR. SRCTYP(ISRC) & 236 & .EQ.'AREAPOLY' ) THEN 237 ! Note that XINIT and YINIT are equivalent values for AREAPOLY 238 QTK = QTK*XINIT*YINIT ! 0 239 ELSEIF ( SRCTYP(ISRC).EQ.'AREACIRC' ) THEN 240 QTK = QTK*PI*RADIUS(ISRC)*RADIUS(ISRC) ! 0 241 ENDIF 242 SYINIT = 0.0 ! 0 243 244 ! Determine Deposition Correction Factors 245 IF ( LDGAS .OR. LWGAS ) THEN 246 CALL PDEPG(X) ! 0 247 ELSE 248 DQCORG = 1.0 ! 0 249 WQCORG = 1.0 250 ENDIF 251 IF ( LDPART .OR. LWPART ) THEN ! 0 252 CALL PDEP(X) ! 0 253 ELSEIF ( NPD.GT.0 ) THEN 254 DO J = 1 , NPD ! 0 255 DQCOR(J) = 1.0 ! 0 256 WQCOR(J) = 1.0 257 ENDDO 258 ENDIF 259 260 ! Define plume centroid height (CENTER) for use in 261 ! inhomogeniety calculations 262 CALL CENTROID(X) ! 0 263 264 ! If the atmosphere is unstable and the stack 265 ! top is below the mixing height, calculate 266 ! the CBL PDF coefficients --- CALL PDF 267 IF ( UNSTAB .AND. (HS.LT.ZI) ) CALL PDF 268 269 ! Determine Effective Plume Height --- CALL HEFF 270 CALL HEFF(X) 271 272 ! Compute effective parameters using an 273 ! iterative average through plume rise layer 274 CALL IBLVAL(X) 275 276 ! Call PDF & HEFF again for final CBL plume heights 277 IF ( UNSTAB .AND. (HS.LT.ZI) ) THEN 278 CALL PDF ! 0 279 CALL HEFF(X) 280 ENDIF 281 282 ! Determine Dispersion Parameters --- CALL VDIS 283 CALL VDIS(X) ! 0 284 285 ! Calculate the 'y-term' contribution to 286 ! dispersion, FSUBY, for coherent plume --- CALL FYPLM 287 CALL FYPLM(SY,FYOUT) 288 FSUBY = FYOUT 289 FSUBYD = FSUBY 290 FSUBYN = FSUBYD 291 292 ! Set lateral term = 0.0 for penetrated source 293 FSUBY3 = 0.0 294 295 ! Check for zero "y-terms"; if zero then skip calculations 296 ! and go to next receptor. 297 IF ( FSUBY.EQ.0.0 .AND. FSUBY3.EQ.0.0 ) THEN 298 DO ITYP = 1 , NUMTYP ! 0 299 HRVAL(ITYP) = 0.0 ! 0 300 IF ( WETSCIM ) HRVALD(ITYP) = 0.0 301 ENDDO 302 303 ELSE 304 305 IF ( NPD.EQ.0 ) THEN ! 0 306 ! Perform calculations for gases 307 ! Assign plume tilt, HV = 0.0 308 309 ADJ = DQCORG*WQCORG ! 0 310 311 IF ( STABLE .OR. (UNSTAB .AND. (HS.GE.ZI)) ) THEN 312 ! Calculate height of the "effective reflecting surface" 313 CALL REFL_HT(HE,X,0.0,VSIGZ,HSBL) ! 0 314 ELSEIF ( UNSTAB ) THEN 315 HSBL = 0.0 ! 0 316 ENDIF 317 318 ! Determine the CRITical Dividing Streamline--- CALL CRITDS 319 CALL CRITDS(HE) ! 0 320 321 ! Calculate the fraction of plume below 322 ! HCRIT, PHEE --- CALL PFRACT 323 CALL PFRACT(HE) 324 325 ! Calculate FOPT = f(PHEE) --- CALL FTERM 326 CALL FTERM 327 328 ! Calculate Conc. for Virtual Point Source --- CALL AER_PCHI 329 CALL AER_PCHI(X,ADJ,VDEPG,0,AEROUT) 330 HRVAL = AEROUT 331 332 ELSE 333 ! Perform calculations for particles, loop through particle sizes 334 335 ! Begin loop over particle sizes 336 DO J = 1 , NPD ! 0 337 338 ! Calculate Plume Tilt Due to Settling, HV 339 HV = (X/US)*VGRAV(J) ! 0 340 341 ! Adjust Jth contribution by mass fraction and source 342 ! depletion 343 ADJ = PHI(J)*DQCOR(J)*WQCOR(J) 344 345 IF ( STABLE .OR. (UNSTAB .AND. (HS.GE.ZI)) ) & 346 & THEN 347 ! Calculate height of the "effective reflecting surface" 348 ! Calculate Settled Plume Height(s), HESETL 349 HESETL = MAX(0.0,HE-HV) ! 0 350 CALL REFL_HT(HESETL,X,0.0,VSIGZ,HSBL) 351 ELSEIF ( UNSTAB ) THEN 352 ! Calculate Settled Plume Height(s), HESETL 353 HESETL = MAX(0.0,0.5*(HED1+HED2)-HV) ! 0 354 HSBL = 0.0 355 ENDIF 356 357 ! Determine the CRITical Dividing Streamline--- CALL CRITDS 358 CALL CRITDS(HESETL) ! 0 359 360 ! Calculate the fraction of plume below 361 ! HCRIT, PHEE --- CALL PFRACT 362 CALL PFRACT(HESETL) 363 364 ! Calculate FOPT = f(PHEE) --- CALL FTERM 365 CALL FTERM 366 367 ! Calculate Conc. for Virtual Point Source --- CALL AER_PCHI 368 CALL AER_PCHI(X,ADJ,VDEP(J),J,AEROUT) 369 HRVAL = HRVAL + AEROUT 370 371 ENDDO 372 ! End loop over particle sizes 373 374 ENDIF 375 ENDIF 376 QTK = QTKSAV ! 0 377 ENDIF 378 379 IF ( PVMRM .AND. .NOT.O3MISS ) THEN ! 271432 380 ! --- Store data by source and receptor for PVMRM option 381 DO ITYP = 1 , NUMTYP ! 0 382 CHI(IREC,ISRC,ITYP) = HRVAL(ITYP) ! 0 383 ENDDO 384 IF ( STABLE .OR. (UNSTAB .AND. (HS.GE.ZI)) ) THEN ! 0 385 HECNTR(IREC,ISRC) = HE ! 0 386 UEFFS(IREC,ISRC) = UEFF 387 ELSE 388 HECNTR(IREC,ISRC) = CENTER ! 0 389 UEFFS(IREC,ISRC) = UEFFD 390 ENDIF 391 IF ( PPF.GT.0.0 ) THEN ! 0 392 HECNTR3(IREC,ISRC) = HE3 ! 0 393 PPFACT(ISRC) = PPF 394 UEFF3S(IREC,ISRC) = UEFF3 395 ELSE 396 PPFACT(ISRC) = 0.0 ! 0 397 ENDIF 398 FOPTS(IREC,ISRC) = FOPT ! 0 399 ! Cycle to next receptor & skip call to SUMVAL (will be done later) 400 GOTO 50 401 ELSEIF ( OLM .AND. .NOT.O3MISS ) THEN 402 ! --- Store data by source and receptor for OLM option 403 DO ITYP = 1 , NUMTYP ! 0 404 CHI(IREC,ISRC,ITYP) = HRVAL(ITYP) ! 0 405 ENDDO 406 ! Cycle to next receptor & skip call to SUMVAL (will be done later) 407 GOTO 50 ! 0 408 ENDIF 409 410 ! Sum HRVAL to AVEVAL and ANNVAL Arrays --- CALL SUMVAL 411 IF ( EVONLY ) THEN ! 271432 412 CALL EV_SUMVAL ! 0 413 ELSE 414 CALL SUMVAL ! 271432 415 ENDIF 416 417 ! Initialize __VAL arrays 418 DO ITYP = 1 , NUMTYP ! 271432 419 HRVAL(ITYP) = 0.0 ! 271432 420 HRVALD(ITYP) = 0.0 421 AERVAL(ITYP) = 0.0 422 AERVALD(ITYP) = 0.0 423 AEROUT(ITYP) = 0.0 424 ENDDO 425 426 50 ENDDO RECEPTOR_LOOP 427 ! End Receptor LOOP 428 429 ! Output 'ARC' Values for EVALFILE --- CALL EVALFL 430 IF ( EVAL(ISRC) ) CALL EVALFL ! 9132 431 432 ENDIF 433 434 CONTINUE ! 9132 435 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