1 2 SUBROUTINE PVMRM_CALC 3 !*********************************************************************** 4 ! PVMRM_CALC Module of the AMS/EPA Regulatory Model - AERMOD 5 ! 6 ! PURPOSE: Processes Hourly Results for PVMRM Option 7 ! 8 ! PROGRAMMER: Roger W. Brode, MACTEC, Inc. (f/k/a PES, Inc.) 9 ! 10 ! DATE: May 12, 2004 11 ! 12 ! INPUTS: 13 ! 14 ! 15 ! OUTPUTS: 16 ! 17 ! 18 ! CALLED FROM: HRLOOP 19 !*********************************************************************** 20 21 ! Variable Declarations 22 USE MAIN1 23 IMPLICIT NONE 24 CHARACTER MODNAM*12 25 CHARACTER LEADIN*25 26 REAL , PARAMETER :: EPS = 0.01 27 REAL :: MAXCHI(NUMREC,NUMTYP) , MAXCONC , CWDIST , CWMAX , CWMIN 28 REAL :: DWDIST , DWMAX , DWMIN , CWDELT , DWDELT 29 REAL :: HMNH , HMXH , HMNT , HMXT , HMNH3 , HMXH3 , HMNT3 , & 30 & HMXT3 , BVERT , BVERT3 , BLONG , BHORIZ , VOLDOM , & 31 & VOLSUM , QSUM , XDOM , YDOM , DISTDOM , UDOM , VOLDOM3 , & 32 & VOLSUM3 , FDOM 33 REAL :: DOMO3MOLES , SUMO3MOLES , DOMNOXMOLES , SUMNOXMOLES , & 34 & DOMCONVERTED , SUMCONVERTED , PERCENTNO2 35 REAL :: AVE_NO2RATIO , SUM_NO2RAT , QAREA 36 REAL :: ZCORR , PCORR , TCORR , PTZ , TAP 37 REAL :: XDEP , WIDTH , LENGTH , XMAXR 38 INTEGER :: DOMIDX(NUMREC,NUMTYP) , IDOM , NDXBLZ , NUMCONT 39 40 SAVE 41 42 ! Variable Initializations 43 MODNAM = 'PVMRM_CALC' ! 0 44 45 ! Determine maximum values and corresponding index by receptor 46 MAXCHI = MAXVAL(CHI,DIM=2) 47 DOMIDX = MAXLOC(CHI,DIM=2) 48 49 ! Begin Receptor LOOP 50 RECEPTOR_LOOP:DO IREC = 1 , NUMREC 51 52 ! Get max concentration and index for dominant source for this receptor 53 MAXCONC = MAXCHI(IREC,1) ! 0 54 IF ( MAXCONC.EQ.0.0 ) GOTO 100 55 56 ! Set index for dominant source, IDOM 57 IDOM = DOMIDX(IREC,1) ! 0 58 UDOM = UEFFS(IREC,IDOM) 59 FDOM = FOPTS(IREC,IDOM) 60 61 ! Determine extent of major contributing sources 62 HMNT = 1.0E10 63 HMXT = 0.0 64 HMNH = 1.0E10 65 HMXH = 0.0 66 HMNT3 = 1.0E10 67 HMXT3 = 0.0 68 HMNH3 = 1.0E10 69 HMXH3 = 0.0 70 CWMIN = 1.0E20 71 CWMAX = -1.0E20 72 DWMIN = 1.0E20 73 DWMAX = -1.0E20 74 75 NUMCONT = 0 76 ! Loop through sources to obtain major contributing sources 77 DO ISRC = 1 , NUMSRC 78 IF ( CHI(IREC,ISRC,1).GE.0.5*MAXCONC ) THEN ! 0 79 NUMCONT = NUMCONT + 1 ! 0 80 IF ( SRCTYP(ISRC)(1:4).EQ.'AREA' ) THEN 81 CALL SETSRC ! 0 82 CALL ARDIST(IREC,XDEP,WIDTH,LENGTH,XMAXR) 83 ! Store max and min values of crosswind and downwind distances 84 CWMAX = MAX(CWMAX,Y+0.5*WIDTH) 85 CWMIN = MIN(CWMIN,Y-0.5*WIDTH) 86 DWMAX = MAX(DWMAX,X+0.5*LENGTH) 87 DWMIN = MIN(DWMIN,X-0.5*LENGTH) 88 ELSE 89 ! Calculate crosswind distance from source to receptor 90 CWDIST = (AXR(IREC)-AXS(ISRC)) & 91 & *WDCOS - (AYR(IREC)-AYS(ISRC))*WDSIN 92 ! Calculate downwind distance from source to receptor 93 DWDIST = -((AXR(IREC)-AXS(ISRC))*WDSIN+(AYR(IREC)-AYS(& 94 & ISRC))*WDCOS) 95 ! Store max and min values of crosswind and downwind distances 96 CWMAX = MAX(CWMAX,CWDIST) 97 CWMIN = MIN(CWMIN,CWDIST) 98 DWMAX = MAX(DWMAX,DWDIST) 99 DWMIN = MIN(DWMIN,DWDIST) 100 ENDIF 101 ! Assign receptor height above stack base for dominant source 102 ZRT = AZELEV(IREC) - AZS(ISRC) + AZFLAG(IREC) ! 0 103 ! Check plume height ranges for horizontal and terrain responding 104 IF ( HECNTR(IREC,ISRC).LT.HMNT ) HMNT = HECNTR(IREC,ISRC) 105 IF ( HECNTR(IREC,ISRC).GT.HMXT ) HMXT = HECNTR(IREC,ISRC) 106 IF ( HECNTR(IREC,ISRC)-ZRT.LT.HMNH ) & 107 & HMNH = HECNTR(IREC,ISRC) - ZRT 108 IF ( HECNTR(IREC,ISRC)-ZRT.GT.HMXH ) & 109 & HMXH = HECNTR(IREC,ISRC) - ZRT 110 IF ( PPFACT(ISRC).GT.0.0 ) THEN 111 IF ( HECNTR3(IREC,ISRC).LT.HMNT3 ) & 112 & HMNT3 = HECNTR3(IREC,ISRC) 113 IF ( HECNTR3(IREC,ISRC).GT.HMXT3 ) & 114 & HMXT3 = HECNTR3(IREC,ISRC) 115 IF ( HECNTR3(IREC,ISRC)-ZRT.LT.HMNH3 ) & 116 & HMNH3 = HECNTR3(IREC,ISRC) - ZRT 117 IF ( HECNTR3(IREC,ISRC)-ZRT.GT.HMXH3 ) & 118 & HMXH3 = HECNTR3(IREC,ISRC) - ZRT 119 ENDIF 120 ENDIF 121 ENDDO 122 123 ! Set vertical dimensions of "box" for major cont. sources. 124 ! Use terrain weighting factor for dominant source (FDOM) to 125 ! combine heights for horizontal and terrain responding plumes. 126 BVERT = FDOM*(HMXH-HMNH) + (1.-FDOM)*(HMXT-HMNT) ! 0 127 IF ( UNSTAB .AND. PPFACT(IDOM).GT.0.0 ) THEN 128 BVERT3 = FDOM*(HMXH3-HMNH3) + (1.-FDOM)*(HMXT3-HMNT3) ! 0 129 ELSE 130 BVERT3 = 0.0 ! 0 131 ENDIF 132 133 ! Set horizontal dimensions of "box" for major cont. sources. 134 IF ( CWMAX.GT.CWMIN ) THEN ! 0 135 BHORIZ = CWMAX - CWMIN ! 0 136 ELSE 137 BHORIZ = 0.0 ! 0 138 ENDIF 139 140 ! Calculate Distance from Dominant Source 141 IF ( SRCTYP(IDOM)(1:4).EQ.'AREA' ) THEN ! 0 142 ISRC = IDOM ! 0 143 CALL SETSRC 144 CALL ARDIST(IREC,XDEP,WIDTH,LENGTH,XMAXR) 145 ! Assign downwind distance based on most upwind vertex of area source 146 XDOM = X 147 YDOM = Y 148 DISTDOM = XMAXR 149 ! Calculate volume of dominant plume, passing WIDTH for BHORIZ 150 CALL PLUME_VOL(DISTDOM,IDOM,0.0,0.0,WIDTH,VOLDOM) 151 ELSE 152 XDOM = -((AXR(IREC)-AXS(IDOM))*WDSIN+(AYR(IREC)-AYS(IDOM)) & 153 & *WDCOS) 154 YDOM = (AXR(IREC)-AXS(IDOM))*WDCOS - (AYR(IREC)-AYS(IDOM)) & 155 & *WDSIN 156 DISTDOM = SQRT(XDOM*XDOM+YDOM*YDOM) 157 ! Calculate volume of dominant plume 158 CALL PLUME_VOL(DISTDOM,IDOM,0.0,0.0,0.0,VOLDOM) 159 ENDIF 160 161 IF ( BVERT.GT.0.0 .OR. BHORIZ.GT.0.0 ) THEN ! 0 162 ! Calculate volume of dominant plus major contributing sources 163 CALL PLUME_VOL(DISTDOM,IDOM,BVERT,BVERT3,BHORIZ,VOLSUM) ! 0 164 ELSE 165 VOLSUM = VOLDOM ! 0 166 ENDIF 167 168 ! Correct plume volume to conditions of standard temperature 169 ! and pressure (Hanrahan, 1999) 170 ! 0.028966 is the average kg/mole of air 171 172 ! First obtain a height to use in the correction 173 ZCORR = AZS(IDOM) + HECNTR(IREC,IDOM) ! 0 174 175 ! Calculate ambient temperature at plume height from pot. temp. profile 176 CALL LOCATE(GRIDHT,1,MXGLVL,ZCORR,NDXBLZ) 177 IF ( NDXBLZ.GE.1 ) THEN 178 !---- Potential temperature 179 CALL GINTRP(GRIDHT(NDXBLZ),GRIDPT(NDXBLZ),GRIDHT(NDXBLZ+1), & 180 & GRIDPT(NDXBLZ+1),ZCORR,PTZ) 181 ELSE 182 ! Use GRID value for lowest level 183 PTZ = GRIDPT(1) ! 0 184 ENDIF 185 TAP = PTZ - GOVRCP*ZCORR ! 0 186 187 PCORR = EXP(-G*0.028966*ZCORR/(RGAS*TAP)) 188 TCORR = 273.15/TAP 189 190 VOLDOM = VOLDOM*PCORR*TCORR 191 VOLSUM = VOLSUM*PCORR*TCORR 192 193 ! get moles of ozone in dominant plume 194 ! O3 expressed in ug/m3 195 DOMO3MOLES = VOLDOM*O3CONC/(10**6*48.) ! O3 is ug/m3 196 197 ! get moles of ozone in the combined plume 198 SUMO3MOLES = VOLSUM*O3CONC/(10**6*48.) ! O3 is ug/m3 199 200 ! get moles of NOx in the dominant plume 201 ! Use molecular weight of NO2 (46) since emission calcs are based on NO2 202 IF ( SRCTYP(IDOM).EQ.'AREA' .OR. SRCTYP(IDOM).EQ.'AREAPOLY' ) & 203 & THEN 204 ISRC = IDOM ! 0 205 CALL SETSRC 206 DOMNOXMOLES = AQS(IDOM)*XINIT*YINIT*DISTDOM/(UDOM*46.) 207 ELSEIF ( SRCTYP(IDOM).EQ.'AREACIRC' ) THEN 208 ISRC = IDOM ! 0 209 CALL SETSRC 210 DOMNOXMOLES = AQS(IDOM)*PI*RADIUS(IDOM)*RADIUS(IDOM) & 211 & *DISTDOM/(UDOM*46.) 212 ELSE 213 DOMNOXMOLES = AQS(IDOM)*DISTDOM/(UDOM*46.) ! 0 214 ENDIF 215 216 ! get moles of NOx in the merged plume 217 ! Sum emissions for sources within projected width of major 218 ! contributing sources 219 ! Assign small delta values for extent of box containing major 220 ! contributing sources to ensure inclusion of sources on the edge 221 CWDELT = MAX(0.1,EPS*(CWMAX-CWMIN)) ! 0 222 DWDELT = MAX(0.1,EPS*(DWMAX-DWMIN)) 223 QSUM = 0.0 224 SUM_NO2RAT = 0.0 225 DO ISRC = 1 , NUMSRC 226 IF ( SRCTYP(ISRC)(1:4).EQ.'AREA' ) THEN ! 0 227 CALL SETSRC ! 0 228 CALL ARDIST(IREC,XDEP,WIDTH,LENGTH,XMAXR) 229 ! Assign CWDIST and DWDIST to Y and X based on distance from 230 ! center of area source to receptor. If center of area source 231 ! is within "box" of major contributing sources, then include 232 ! it's emissions in calculation of NOx moles. 233 CWDIST = Y 234 DWDIST = X 235 ELSE 236 ! Calculate crosswind distance from source to receptor 237 CWDIST = (AXR(IREC)-AXS(ISRC)) & 238 & *WDCOS - (AYR(IREC)-AYS(ISRC))*WDSIN 239 ! Calculate downwind distance from source to receptor 240 DWDIST = - & 241 & ((AXR(IREC)-AXS(ISRC))*WDSIN+(AYR(IREC)-AYS(ISRC& 242 & ))*WDCOS) 243 ENDIF 244 ! Check for crosswind disance between MIN and MAX of projected width 245 IF ( CWDIST.GE.CWMIN-CWDELT .AND. CWDIST.LE.CWMAX+CWDELT ) & 246 & THEN 247 IF ( DWDIST.GE.DWMIN-DWDELT .AND. DWDIST.LE.DWMAX+ & 248 & DWDELT ) THEN 249 IF ( SRCTYP(ISRC).EQ.'AREA' .OR. SRCTYP(ISRC) & 250 & .EQ.'AREAPOLY' ) THEN 251 QAREA = AQS(ISRC)*XINIT*YINIT ! 0 252 QSUM = QSUM + QAREA 253 SUM_NO2RAT = SUM_NO2RAT + ANO2_RATIO(ISRC)*QAREA 254 ELSEIF ( SRCTYP(ISRC).EQ.'AREACIRC' ) THEN 255 QAREA = AQS(ISRC)*PI*RADIUS(ISRC)*RADIUS(ISRC) ! 0 256 QSUM = QSUM + QAREA 257 SUM_NO2RAT = SUM_NO2RAT + ANO2_RATIO(ISRC)*QAREA 258 ELSE 259 QSUM = QSUM + AQS(ISRC) ! 0 260 SUM_NO2RAT = SUM_NO2RAT + ANO2_RATIO(ISRC) & 261 & *AQS(ISRC) 262 ENDIF 263 ENDIF 264 ENDIF 265 ENDDO 266 267 !TMP Check for QSUM = 0.0 error (i.e., dominant source not inc. in sum) 268 !TMP This condition should never be encountered 269 IF ( QSUM.EQ.0.0 ) THEN ! 0 270 WRITE (IOUNIT,*) & 271 & 'Error in PVMRM_CALC; QSUM = 0.0. Aborting!!!' 272 WRITE (*,*) 'Error in PVMRM_CALC; QSUM = 0.0. Aborting!!!' 273 STOP 274 ENDIF 275 276 ! Calculate NOx moles for combined plume, and average in-stack ratio 277 ! Use molecular weight of NO2 (46) since emission calcs are based on NO2 278 SUMNOXMOLES = QSUM*DISTDOM/(UDOM*46.) ! 0 279 AVE_NO2RATIO = SUM_NO2RAT/QSUM 280 281 ! Calculate NOx conversion ratios for dominant plume and combined plume 282 DOMCONVERTED = DOMO3MOLES/(DOMNOXMOLES) 283 SUMCONVERTED = SUMO3MOLES/(SUMNOXMOLES) 284 285 ! 286 ! Find which is more important -- the dominant plume or the combined plume 287 ! Minimum used as the most controlling plumes will have the least conversion 288 IF ( DOMCONVERTED.LT.SUMCONVERTED ) THEN 289 PERCENTNO2 = DOMCONVERTED + ANO2_RATIO(IDOM) ! 0 290 ELSE 291 PERCENTNO2 = SUMCONVERTED + AVE_NO2RATIO ! 0 292 ENDIF 293 294 ! Limit to equilibrium concentration of NO2 (default set at 90 percent) 295 IF ( PERCENTNO2.GT.NO2EQUIL ) PERCENTNO2 = NO2EQUIL ! 0 296 297 ! Write data to PVMRM.TXT debugging file 298 IF ( DEBUG .AND. DOMCONVERTED.LT.SUMCONVERTED ) THEN 299 WRITE (50,9987) KURDAT , IREC , SRCID(IDOM) , DISTDOM , & 300 & MAXCONC , NUMCONT , O3CONC , DOMO3MOLES , & 301 & DOMNOXMOLES , BHORIZ , BVERT , VOLDOM , & 302 & PERCENTNO2 303 9987 FORMAT (1x,'DOM: ',i8,i5,1x,a8,2F11.4,1x,i4,1x,f10.3, & 304 & 2(1x,f12.3),2(1x,f10.3),1x,e12.4,1x,f8.3) 305 ELSEIF ( DEBUG ) THEN 306 WRITE (50,9988) KURDAT , IREC , SRCID(IDOM) , DISTDOM , & 307 & MAXCONC , NUMCONT , O3CONC , SUMO3MOLES , & 308 & SUMNOXMOLES , BHORIZ , BVERT , VOLSUM , & 309 & PERCENTNO2 310 9988 FORMAT (1x,'SUM: ',i8,i5,1x,a8,2F11.4,1x,i4,1x,f10.3, & 311 & 2(1x,f12.3),2(1x,f10.3),1x,e12.4,1x,f8.3) 312 ENDIF 313 314 ! Update HRVAL, AVEVAL and ANNVAL Arrays 315 DO ISRC = 1 , NUMSRC ! 0 316 DO ITYP = 1 , NUMTYP ! 0 317 HRVAL(ITYP) = CHI(IREC,ISRC,ITYP)*PERCENTNO2 ! 0 318 ! Sum HRVAL to AVEVAL and ANNVAL Arrays --- CALL SUMVAL 319 ENDDO 320 IF ( EVONLY ) THEN ! 0 321 CALL EV_SUMVAL ! 0 322 ELSE 323 CALL SUMVAL ! 0 324 ENDIF 325 ! Check ARC centerline values for EVALFILE 326 ! output --- CALL EVALCK 327 IF ( EVAL(ISRC) ) CALL EVALCK ! 0 328 ENDDO 329 330 ! Initialize __VAL arrays 331 DO ITYP = 1 , NUMTYP ! 0 332 HRVAL(ITYP) = 0.0 ! 0 333 HRVALD(ITYP) = 0.0 334 AERVAL(ITYP) = 0.0 335 AERVALD(ITYP) = 0.0 336 ENDDO 337 338 100 ENDDO RECEPTOR_LOOP 339 ! End Receptor LOOP 340 341 CONTINUE ! 0 342 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