module_bioemi_beis311.F
References to this file elsewhere.
1 MODULE module_bioemi_beis311
2
3 ! BEIS3.11 Emissions Module for WRF-Chem
4 ! Written by Greg Frost 6/2004
5 ! Using off-line gridded standard biogenic emissions
6 ! for each model compound with such emissions,
7 ! model shortwave solar flux (isoprene only),
8 ! & air temperature, pressure, and density in lowest model level,
9 ! calculates actual biogenic emissions of each compound.
10 ! Based on hrbeis311.f from BEIS3.11 for SMOKE, with major
11 ! surgery performed on original routines for use with WRF-Chem.
12 ! This version assumes chemical mechanism is RACM.
13 ! The following 16 RACM compounds have biogenic emissions:
14 ! iso, no, oli, api, lim, xyl, hc3, ete, olt, ket, ald, hcho, eth, ora2, co, nr
15 !23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
16
17 CONTAINS
18 SUBROUTINE bio_emissions_beis311(id,config_flags,ktau,dtstep, &
19 julday,gmt,xlat,xlong,t_phy,p_phy,gsw, &
20 sebio_iso,sebio_oli,sebio_api,sebio_lim,sebio_xyl, &
21 sebio_hc3,sebio_ete,sebio_olt,sebio_ket,sebio_ald, &
22 sebio_hcho,sebio_eth,sebio_ora2,sebio_co,sebio_nr, &
23 noag_grow,noag_nongrow,nononag,slai, &
24 ebio_iso,ebio_oli,ebio_api,ebio_lim,ebio_xyl, &
25 ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
26 ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_nr,ebio_no, &
27 ids,ide, jds,jde, kds,kde, &
28 ims,ime, jms,jme, kms,kme, &
29 its,ite, jts,jte, kts,kte )
30
31 USE module_configure
32 USE module_state_description
33
34 IMPLICIT NONE
35
36 ! .. Parameters ..
37 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
38
39 ! .. Indices ..
40 INTEGER, INTENT(IN ) :: id,ktau, &
41 ids,ide, jds,jde, kds,kde, &
42 ims,ime, jms,jme, kms,kme, &
43 its,ite, jts,jte, kts,kte
44 ! .. Passed variables ..
45 INTEGER, INTENT (IN) :: julday ! current simulation julian day
46
47 REAL, INTENT (IN) :: gmt,dtstep
48
49 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
50 INTENT(IN ) :: &
51 t_phy, & !air T (K)
52 p_phy !P (Pa)
53
54 REAL, DIMENSION( ims:ime , jms:jme ), &
55 INTENT(IN ) :: &
56 xlat, & !latitude (deg)
57 xlong, & !longitude (deg)
58 gsw !downward shortwave surface flux (W/m^2)
59
60 ! Normalized biogenic emissions for standard conditions (moles compound/km^2/hr)
61 REAL, DIMENSION( ims:ime , jms:jme ), &
62 INTENT(IN ) :: &
63 sebio_iso,sebio_oli,sebio_api,sebio_lim,sebio_xyl, &
64 sebio_hc3,sebio_ete,sebio_olt,sebio_ket,sebio_ald, &
65 sebio_hcho,sebio_eth,sebio_ora2,sebio_co,sebio_nr, &
66 noag_grow,noag_nongrow,nononag
67
68 ! Leaf area index for isoprene
69 REAL, DIMENSION( ims:ime , jms:jme ), &
70 INTENT(IN ) :: slai
71
72 ! Actual biogenic emissions (moles compound/km^2/hr)
73 REAL, DIMENSION( ims:ime , jms:jme ), &
74 INTENT(INOUT ) :: &
75 ebio_iso,ebio_oli,ebio_api,ebio_lim,ebio_xyl, &
76 ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
77 ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_nr,ebio_no
78
79
80 ! .. Local Scalars ..
81
82 INTEGER :: i,j
83
84 ! Variables for 1 element of I/O arrays
85 ! met and phys input variables
86 REAL :: tair ! surface air temperature (K)
87 REAL :: tsolar ! downward shortwave surface flux (W/m^2)
88 REAL :: pres ! surface pressure (mb)
89 REAL :: ylat ! latitude (deg)
90 REAL :: ylong ! longitude (deg)
91 ! normalized emissions (moles compound/km^2/hr)
92 REAL :: se_iso,se_oli,se_api,se_lim,se_xyl, &
93 se_hc3,se_ete,se_olt,se_ket,se_ald, &
94 se_hcho,se_eth,se_ora2,se_co,se_nr, &
95 growagno,ngrowagno,nonagno
96 ! leaf area index for isoprene
97 REAL :: tlai
98 ! actual emissions for NO
99 REAL :: e_no
100
101 ! Other parameters needed in calculations
102 ! Guenther's parameterizations: Guenther et al. JGR 98, 12609-12617, 1993
103 REAL :: ct, dt ! Guenther's temperature correction for isoprene
104 REAL :: cfno ! NO correction factor
105 REAL :: cfovoc ! non-isoprene correction factor
106 REAL :: par ! PAR = photosynthetically active radiation (micromole/m^2/s)
107 REAL :: csubl ! C sub l from Guenther
108 REAL :: zen ! zenith angle (radians)
109 REAL :: coszen ! cosine(zenith angle)
110 REAL :: pardb ! PAR direct beam
111 REAL :: pardif ! PAR diffuse
112 REAL :: gmtp ! current simulation time
113
114 ! Error message variables
115 INTEGER , PARAMETER :: ldev = 6 ! unit number for log file
116 CHARACTER*256 :: mesg
117
118 ! Functions called directly or indirectly
119 ! clnew calculates csubl based on zenith angle, par, and lai
120 ! cguen Guenther's equation for computing light correction
121 ! fertilizer_adj computes fertlizer adjustment factor
122 ! veg_adj computes vegatation adjustment factor
123 ! growseason computes day of growing season
124
125 ! Subroutines called directly or indirectly
126 ! calc_zenithb calculates zenith angle from latitude, longitude, julian day, and GMT
127 ! NOTE: longitude input for this routine is nonstandard: >0 for W, <0 for E!!
128 ! getpar computes PAR (direct beam and diffuse) in umol/m2-sec from downward shortwave flux
129 ! hrno algorithm to estimate NO emissions; does not include precipitation adjustment
130
131 !***************************************
132 ! begin body of subroutine bio_emissions_beis311
133
134 ! hour into integration
135 gmtp=float(ktau)*dtstep/3600.
136 !
137 gmtp=mod(gmt+gmtp,24.)
138 write(mesg,*) 'calculate beis311 emissions at gmtp = ',gmtp
139 call wrf_debug(15,mesg)
140 DO 100 j=jts,jte
141 DO 100 i=its,ite
142
143 tair = t_phy(i,kts,j)
144 pres = .01*p_phy(i,kts,j)
145 ylat = xlat(i,j)
146 ylong = xlong(i,j)
147 tsolar = gsw(i,j)
148 tlai = slai(i,j)
149 se_iso = sebio_iso(i,j)
150 se_oli = sebio_oli(i,j)
151 se_api = sebio_api(i,j)
152 se_lim = sebio_lim(i,j)
153 se_xyl = sebio_xyl(i,j)
154 se_hc3 = sebio_hc3(i,j)
155 se_ete = sebio_ete(i,j)
156 se_olt = sebio_olt(i,j)
157 se_ket = sebio_ket(i,j)
158 se_ald = sebio_ald(i,j)
159 se_hcho = sebio_hcho(i,j)
160 se_eth = sebio_eth(i,j)
161 se_ora2 = sebio_ora2(i,j)
162 se_co = sebio_co(i,j)
163 se_nr = sebio_nr(i,j)
164 growagno = noag_grow(i,j)
165 ngrowagno = noag_nongrow(i,j)
166 nonagno = nononag(i,j)
167
168 !....Perform checks on max and min bounds for temperature
169
170 IF (tair .LT. 200.0) THEN
171 ! WRITE( mesg, 94010 )
172 ! & 'tair=', tair,
173 ! & 'too low at i,j= ',i,',',j
174 WRITE( ldev, * ) mesg
175 END IF
176
177 IF (tair .GT. 315.0 ) THEN
178 ! WRITE( mesg, 94020 )
179 ! & 'tair=', tair,
180 ! & 'too high at i,j= ',i,',',j,
181 ! & '...resetting to 315K'
182 tair = 315.0
183 ! WRITE( ldev, * ) mesg
184 ENDIF
185
186 !... Isoprene emissions
187 !...... Calculate temperature correction term
188 dt = 28668.514 / tair
189 ct = EXP( 37.711 - 0.398570815 * dt ) / &
190 (1.0 + EXP( 91.301 - dt ) )
191
192 !...... Calculate zenith angle in radians
193 ! NOTE: nonstandard longitude input here: >0 for W, <0 for E!!
194 CALL calc_zenithb(ylat,-ylong,julday,gmtp,zen)
195 coszen = COS(zen)
196
197 !...... Convert tsolar to PAR and find direct and diffuse fractions
198 CALL getpar( tsolar, pres, zen, pardb, pardif )
199 par = pardb + pardif
200
201 !...... Check max/min bounds of PAR and calculate
202 !...... biogenic isoprene
203 IF ( par .LT. 0.00 .OR. par .GT. 2600.0 ) THEN
204 ! WRITE( mesg, 94010 )
205 ! & 'PAR=', par,
206 ! & 'out of range at i,j= ',i,',',j
207 ! WRITE( ldev, * ) mesg
208 ENDIF
209
210 !...... Check max bound of LAI
211 IF ( tlai .GT. 10.0 ) THEN
212 ! WRITE( mesg, 94010 )
213 ! & 'LAI=', tlai,
214 ! & 'out of range at i,j= ',i,',',j
215 ! WRITE( ldev, * ) mesg
216 ENDIF
217
218 !...... Initialize csubl
219 csubl = 0.0
220
221 !...... If PAR < 0.01 or zenith angle > 89 deg, set isoprene emissions to 0.
222 IF ( par .LE. 0.01 .OR. coszen .LE. 0.02079483 ) THEN
223 ebio_iso(i,j) = 0.0
224
225 ELSE
226
227 !...... Calculate csubl including shading if LAI > 0.1
228 IF ( tlai .GT. 0.1 ) THEN
229 csubl = clnew( zen, pardb, pardif, tlai )
230
231 !...... Otherwise calculate csubl without considering LAI
232 ELSE ! keep this or not?
233 csubl = cguen( par )
234
235 ENDIF
236
237 ebio_iso(i,j) = se_iso * ct * csubl
238
239 ENDIF
240
241
242 !... Other biogenic emissions except NO:
243 !...... RACM: oli, api, lim, hc3, ete, olt, ket, ald, hcho, eth, ora2, co
244
245 cfovoc = EXP( 0.09 * ( tair - 303.0 ) )
246
247 ebio_oli(i,j) = se_oli * cfovoc
248 ebio_api(i,j) = se_api * cfovoc
249 ebio_lim(i,j) = se_lim * cfovoc
250 ebio_xyl(i,j) = se_xyl * cfovoc
251 ebio_hc3(i,j) = se_hc3 * cfovoc
252 ebio_ete(i,j) = se_ete * cfovoc
253 ebio_olt(i,j) = se_olt * cfovoc
254 ebio_ket(i,j) = se_ket * cfovoc
255 ebio_ald(i,j) = se_ald * cfovoc
256 ebio_hcho(i,j) = se_hcho * cfovoc
257 ebio_eth(i,j) = se_eth * cfovoc
258 ebio_ora2(i,j) = se_ora2 * cfovoc
259 ebio_co(i,j) = se_co * cfovoc
260 ebio_nr(i,j) = se_nr * cfovoc
261
262 !... NO emissions
263
264 CALL hrno( julday, growagno, ngrowagno, nonagno, tair, e_no)
265
266 ebio_no(i,j) = e_no
267
268 100 CONTINUE
269
270 RETURN
271
272
273 !****************** FORMAT STATEMENTS ******************************
274
275 !........... Informational (LOG) message formats... 92xxx
276
277
278 !........... Internal buffering formats............ 94xxx
279
280
281 94010 FORMAT( A, F10.2, 1X, A, I4, A1, I4)
282 94020 FORMAT( A, F10.2, 1X, A, I4, A1, I4, 1X, A)
283
284
285 !***************** CONTAINS ********************************************
286 CONTAINS
287
288 REAL FUNCTION clnew( zen, pardb, pardif, tlai )
289
290 !........ Function to calculate csubl based on zenith angle, PAR, and LAI
291
292 IMPLICIT NONE
293
294 REAL, INTENT (IN) :: pardb ! direct beam PAR( umol/m2-s)
295 REAL, INTENT (IN) :: pardif ! diffuse PAR ( umol/m2-s)
296 REAL, INTENT (IN) :: zen ! solar zenith angle (radians)
297 REAL, INTENT (IN) :: tlai ! leaf area index for grid cell
298 REAL kbe ! extinction coefficient for direct beam
299 REAL canparscat ! exponentially wtd scattered PAR (umol/m2-s)
300 REAL canpardif ! exponentially wtd diffuse PAR (umol/m2-s)
301 REAL parshade ! PAR on shaded leaves (umol/m2-s)
302 REAL parsun ! PAR on sunlit leaves (umol/m2-s)
303 REAL laisun ! LAI that is sunlit
304 REAL fracsun ! fraction of leaves that are sunlit
305 REAL fracshade ! fraction of leaves that are shaded
306
307 !........... CN98 - eqn 15.4, assume x=1
308
309 kbe = 0.5 * SQRT(1. + TAN( zen ) * TAN( zen ))
310
311 !.......... CN98 - p. 261 (this is usually small)
312
313 canparscat = 0.5 * pardb * (EXP(-0.894 * kbe * tlai) - &
314 EXP(-1.* kbe * tlai))
315
316 !.......... CN98 - p. 261 (assume exponentially wtd avg)
317
318 canpardif = pardif * (1. - EXP(-0.61 * tlai))/(0.61 * tlai)
319
320 !......... CN98 - p. 261 (for next 3 eqns)
321
322 parshade = canpardif + canparscat
323 parsun = kbe * (pardb + pardif) + parshade
324 laisun = (1. - EXP( -kbe * tlai))/kbe
325 fracsun = laisun/tlai
326 fracshade = 1. - fracsun
327
328 !.......... cguen is guenther's eqn for computing light correction as a
329 !.......... function of PAR...fracSun should probably be higher since
330 !.......... sunlit leaves tend to be thicker than shaded leaves. But
331 !.......... since we need to make crude asmptns regarding leave
332 !.......... orientation (x=1), will not attempt to fix at the moment.
333
334 clnew =fracsun * cguen(parsun) + fracshade * cguen(parshade)
335
336 RETURN
337 END FUNCTION clnew
338
339 REAL FUNCTION cguen( partmp )
340
341 !.......... Guenther's equation for computing light correction
342
343 IMPLICIT NONE
344 REAL, INTENT (IN) :: partmp
345 REAL, PARAMETER :: alpha2 = 0.00000729
346
347 IF ( partmp .LE. 0.01) THEN
348 cguen = 0.0
349 ELSE
350 cguen = (0.0028782 * partmp) / &
351 SQRT(1. + alpha2 * partmp * partmp)
352 ENDIF
353
354 RETURN
355 END FUNCTION cguen
356
357 END SUBROUTINE bio_emissions_beis311
358
359 !=================================================================
360
361 SUBROUTINE calc_zenithb(lat,long,ijd,gmt,zenith)
362 ! Based on calc_zenith from WRF-Chem module_phot_mad.F
363 ! this subroutine calculates solar zenith angle for a
364 ! time and location. must specify:
365 ! input:
366 ! lat - latitude in decimal degrees
367 ! long - longitude in decimal degrees
368 ! NOTE: Nonstandard convention for long: >0 for W, <0 for E!!
369 ! gmt - greenwich mean time - decimal military eg.
370 ! 22.75 = 45 min after ten pm gmt
371 ! output
372 ! zenith - in radians (GJF, 6/2004)
373 ! remove azimuth angle calculation since not needed (GJF, 6/2004)
374 ! .. Scalar Arguments ..
375 REAL :: gmt, lat, long, zenith
376 INTEGER :: ijd
377 ! .. Local Scalars ..
378 REAL :: csz, cw, d, decl, dr, ec, epsi, eqt, eyt, feqt, feqt1, &
379 feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
380 pi, ra, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, &
381 yt, zpt, zr
382 INTEGER :: jd
383 CHARACTER*256 :: mesg
384 ! .. Intrinsic Functions ..
385 INTRINSIC acos, atan, cos, min, sin, tan
386 ! convert to radians
387 pi = 3.1415926535590
388 dr = pi/180.
389 rlt = lat*dr
390 rphi = long*dr
391
392 ! ???? + (yr - yref)
393
394 jd = ijd
395
396 d = jd + gmt/24.0
397 ! calc geom mean longitude
398 ml = 279.2801988 + .9856473354*d + 2.267E-13*d*d
399 rml = ml*dr
400
401 ! calc equation of time in sec
402 ! w = mean long of perigee
403 ! e = eccentricity
404 ! epsi = mean obliquity of ecliptic
405 w = 282.4932328 + 4.70684E-5*d + 3.39E-13*d*d
406 wr = w*dr
407 ec = 1.6720041E-2 - 1.1444E-9*d - 9.4E-17*d*d
408 epsi = 23.44266511 - 3.5626E-7*d - 1.23E-15*d*d
409 pepsi = epsi*dr
410 yt = (tan(pepsi/2.0))**2
411 cw = cos(wr)
412 sw = sin(wr)
413 ssw = sin(2.0*wr)
414 eyt = 2.*ec*yt
415 feqt1 = sin(rml)*(-eyt*cw-2.*ec*cw)
416 feqt2 = cos(rml)*(2.*ec*sw-eyt*sw)
417 feqt3 = sin(2.*rml)*(yt-(5.*ec**2/4.)*(cw**2-sw**2))
418 feqt4 = cos(2.*rml)*(5.*ec**2*ssw/4.)
419 feqt5 = sin(3.*rml)*(eyt*cw)
420 feqt6 = cos(3.*rml)*(-eyt*sw)
421 feqt7 = -sin(4.*rml)*(.5*yt**2)
422 feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
423 eqt = feqt*13751.0
424
425 ! convert eq of time from sec to deg
426 reqt = eqt/240.
427 ! calc right ascension in rads
428 ra = ml - reqt
429 rra = ra*dr
430 ! calc declination in rads, deg
431 tab = 0.43360*sin(rra)
432 rdecl = atan(tab)
433 decl = rdecl/dr
434 ! calc local hour angle
435 lbgmt = 12.0 - eqt/3600. + long*24./360.
436 lzgmt = 15.0*(gmt-lbgmt)
437 zpt = lzgmt*dr
438 csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
439 if(csz.gt.1) then
440 write(mesg,*) 'calczen,csz ',csz
441 call wrf_debug(15,mesg)
442 endif
443 csz = min(1.,csz)
444 zr = acos(csz)
445 ! zenith = zr/dr
446 ! keep zenith angle in radians for later use (GJF 6/2004)
447 zenith = zr
448
449 RETURN
450
451 END SUBROUTINE calc_zenithb
452
453 !=================================================================
454
455
456 SUBROUTINE getpar( tsolar, pres, zen, pardb, pardif )
457
458 !***********************************************************************
459 ! subroutine body starts at line
460 !
461 ! DESCRIPTION:
462 !
463 ! Based on code from Bart Brashers (10/2000), which was based on
464 ! code from Weiss and Norman (1985).
465 !
466 !
467 ! PRECONDITIONS REQUIRED:
468 ! Solar radiation (W/m2) and pressure (mb)
469 !
470 ! SUBROUTINES AND FUNCTIONS CALLED:
471 !
472 ! REVISION HISTORY:
473 ! 3/01 : Prototype by JMV
474 !
475 !***********************************************************************
476 !
477 ! Project Title: Sparse Matrix Operator Kernel Emissions (SMOKE) Modeling
478 ! System
479 ! File: @(#)Id: getpar.f,v 1.1.1.1 2001/03/27 19:08:49 smith_w Exp
480 !
481 ! COPYRIGHT (C) 2001, MCNC--North Carolina Supercomputing Center
482 ! All Rights Reserved
483 !
484 ! See file COPYRIGHT for conditions of use.
485 !
486 ! MCNC-Environmental Programs Group
487 ! P.O. Box 12889
488 ! Research Triangle Park, NC 27709-2889
489 !
490 ! env_progs@mcnc.org
491 !
492 ! Pathname: Source: /env/proj/archive/cvs/jmv/beis3v0.9/getpar.f,v
493 ! Last updated: Date: 2001/03/27 19:08:49
494 !
495 !***********************************************************************
496
497 IMPLICIT NONE
498
499 !........ Inputs
500
501 REAL , INTENT (IN) :: tsolar ! modeled or observed total radiation (W/m2)
502 REAL , INTENT (IN) :: pres ! atmospheric pressure (mb)
503 REAL , INTENT (IN) :: zen ! solar zenith angle (radians)
504
505 !........ Outputs
506
507 REAL, INTENT (OUT) :: pardb ! direct beam PAR( umol/m2-s)
508 REAL, INTENT (OUT) :: pardif ! diffuse PAR ( umol/m2-s)
509
510 !........... PARAMETERS and their descriptions:
511
512 REAL, PARAMETER :: watt2umol = 4.6 ! convert W/m^2 to umol/m^2-s (4.6)
513
514 !
515 REAL ratio ! transmission fraction for total radiation
516 REAL ot ! optical thickness
517 REAL rdvis ! possible direct visible beam (W/m^2)
518 REAL rfvis ! possible visible diffuse (W/m^2)
519 REAL wa ! water absorption in near-IR (W/m^2)
520 REAL rdir ! direct beam in near-IR (W/m^2)
521 REAL rfir ! diffuse near-IR (W/m^2)
522 REAL rvt ! total possible visible radiation (W/m^2)
523 REAL rirt ! total possible near-IR radiation (W/m^2)
524 REAL fvis ! fraction of visible to total
525 REAL fvb ! fraction of visible that is direct beam
526 REAL fvd ! fraction of visible that is diffuse
527
528 !***************************************
529 ! begin body of subroutine
530
531 !............ Assume that PAR = 0 if zenith angle is greater than 87 degrees
532 !............ or if solar radiation is zero
533
534 IF (zen .GE. 1.51844 .OR. tsolar .LE. 0.) THEN
535 pardb = 0.
536 pardif = 0.
537 RETURN
538 ENDIF
539
540 !............ Compute clear sky (aka potential) radiation terms
541
542 ot = pres / 1013.25 / COS(zen) !Atmospheric Optical thickness
543 rdvis = 600. * EXP(-.185 * ot) * COS(zen) !Direct visible beam, eqn (1)
544 rfvis = 0.42 * (600 - rdvis) * COS(zen) !Visible Diffuse, eqn (3)
545 wa = 1320 * .077 * (2. * ot)**0.3 !water absorption in near-IR, eqn (6)
546 rdir = (720. * EXP(-0.06 * ot)-wa) * COS(zen) !Direct beam near-IR, eqn (4)
547 rfir = 0.65 * (720. - wa - rdir) * COS(zen) !Diffuse near-IR, eqn (5)
548
549 rvt = rdvis + rfvis !Total visible radiation, eqn (9)
550 rirt = rdir + rfir !Total near-IR radiation, eqn (10)
551 fvis = rvt/(rirt + rvt) !Fraction of visible to total radiation, eqn 7
552 ratio = tsolar /(rirt + rvt) !Ratio of "actual" to clear sky solar radiation
553
554 !............ Compute fraction of visible that is direct beam
555
556 IF (ratio .GE. 0.89) THEN
557 fvb = rdvis/rvt * 0.941124
558 ELSE IF (ratio .LE. 0.21) THEN
559 fvb = rdvis/rvt * 9.55E-3
560 ELSE
561 fvb = rdvis/rvt * (1.-((0.9 - ratio)/0.7)**0.666667)
562 ENDIF
563 fvd = 1. - fvb
564
565 !............ Compute PAR (direct beam and diffuse) in umol/m2-sec
566
567 pardb = tsolar * fvis * fvb * watt2umol
568 pardif = tsolar * fvis * fvd * watt2umol
569
570
571 RETURN
572
573 !****************** FORMAT STATEMENTS ******************************
574
575 !........... Informational (LOG) message formats... 92xxx
576
577
578 !........... Internal buffering formats............ 94xxx
579
580 END SUBROUTINE getpar
581
582 SUBROUTINE hrno( julday, growagno, ngrowagno, nonagno, tairin, e_no)
583
584 !***********************************************************************
585 ! subroutine body starts at line 150
586 !
587 ! DESCRIPTION:
588 !
589 ! Uses new NO algorithm NO = Normalized*Tadj*Fadj*Cadj
590 ! to estimate NO emissions
591 ! Information needed to estimate NO emissions
592 ! Julian Day (integer) julday
593 ! Surface Temperature (MCIP field) tair (K)
594 ! Note: Precipitation adjustment not used in the WRF-Chem implementation of BEIS3.11
595 ! because of differences in soil categories between BEIS and WRF-Chem
596 !
597 ! The calculation are based on the following paper by J.J. Yienger and H. Levy II
598 ! J.J. Yienger and H. Levy II, Journal of Geophysical Research, vol 100,11447-11464,1995
599 !
600 ! Also see the following paper for more information:
601 ! Proceedings of the Air and Waste Management Association/U.S. Environmental Protection
602 ! Agency EMission Inventory Conference, Raleigh October 26-28, 1999 Raleigh NC
603 ! by Tom Pierce and Lucille Bender
604 !
605 ! REFERENCES
606 !
607 ! JACQUEMIN B. AND NOILHAN J. (1990), BOUND.-LAYER METEOROL., 52, 93-134.
608 ! J.J. Yienger and H. Levy II, Journal of Geophysical Research, vol 100,11447-11464,1995
609 ! T. Pierce and L. Bender, Examining the Temporal Variability of Ammonia and Nitric Oxide Emissions from Agricultural Processes
610 ! Proceedings of the Air and Waste Management Association/U.S. Environmental Protection
611 ! Agency EMission Inventory Conference, Raleigh October 26-28, 1999 Raleigh NC
612 !
613 ! PRECONDITIONS REQUIRED:
614 ! Normalized NO emissions, Surface Temperature
615 !
616 ! SUBROUTINES AND FUNCTIONS CALLED (directly or indirectly):
617 ! fertilizer_adj computes fertlizer adjustment factor
618 ! veg_adj computes vegatation adjustment factor
619 ! growseason computes day of growing season
620 !
621 !
622 ! REVISION HISTORY:
623 ! 10/01 : Prototype by GAP
624 !
625 !***********************************************************************
626 !
627 ! Project Title: BEIS3 Enhancements for NO emission calculation
628 ! File: hrno.f
629 !
630 !
631 !***********************************************************************
632
633 IMPLICIT NONE
634
635 !........... ARGUMENTS and their descriptions:
636
637 INTEGER, INTENT (IN) :: julday ! current julian day
638
639
640 REAL, INTENT (IN) :: tairin ! air temperature (K)
641 REAL, INTENT (IN) :: growagno ! norm NO emissions, agricultural, growing
642 REAL, INTENT (IN) :: ngrowagno ! norm NO emissions, agricultural, not growing
643 REAL, INTENT (IN) :: nonagno ! norm NO emissions, non-agricultural
644
645 REAL, INTENT (OUT) :: e_no ! output NO emissions
646
647 !........... SCRATCH LOCAL VARIABLES and their descriptions:
648
649 REAL cfno ! NO correction factor
650 REAL cfnograss ! NO correction factor for grasslands
651 REAL tsoi ! soil temperature
652 REAL tair ! air temperature
653
654 REAL :: cfnowet, cfnodry
655
656 INTEGER growseason, gday
657 EXTERNAL growseason
658 !***********************************************************************
659
660 tair = tairin
661
662 !............. calculate NO emissions by going thru temperature cases
663
664 ! gday = growseason(julday)
665 gday = 91
666 IF (gday .eq. 0) THEN !not growing season
667 IF ( tair .GT. 303.00 ) tair = 303.00
668
669 IF ( tair .GT. 268.8690 ) THEN
670 cfno = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
671 ELSE
672 cfno = 0.0
673 ENDIF
674
675 e_no = &
676 ngrowagno * cfno & !agriculture
677 + nonagno * cfno ! non-agriculture
678
679 ELSE
680
681 tsoi = 0.72*tair+82.28
682 IF (tsoi .LE. 273.16) tsoi = 273.16
683 IF (tsoi .GE. 303.16) tsoi = 303.16
684
685 cfnodry = (1./3.)*(1./30.)*(tsoi-273.16) ! see YL 1995 Equa 9a p. 11452
686 IF (tsoi .LE. 283.16) THEN ! linear cold case
687 cfnowet = (tsoi-273.16)*EXP(-0.103*30.0)*0.28 ! see YL 1995 Equ 7b
688 ELSE ! exponential case
689 cfnowet = EXP(0.103*(tsoi-273.16)) &
690 *EXP(-0.103*30.0)
691 ENDIF
692 cfno = 0.5*cfnowet + 0.5*cfnodry
693
694 IF ( tair .GT. 303.00 ) tair = 303.00
695
696 IF ( tair .GT. 268.8690 ) THEN
697 cfnograss = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
698 ELSE
699 cfnograss = 0.0
700 ENDIF
701
702 e_no = growagno * cfno *fertilizer_adj(julday)*veg_adj(julday) &
703 + nonagno * cfnograss
704
705 ENDIF
706
707 RETURN
708
709 !***************** CONTAINS ********************************************
710 CONTAINS
711
712 REAL FUNCTION fertilizer_adj(julday)
713 !*****************************************************************
714 !
715 ! SUMMARY:
716 ! computes fertilizer adjustment factor from Julian day
717 !
718 ! FUNCTION CALLS:
719 ! growseason computes day of growing season
720 !
721 ! NOTE: julday = Julian day format
722 !
723 !*****************************************************************
724 implicit none
725 integer julday
726 !
727 !******** local scratch variables
728 !
729 integer gday
730 !
731 !******** function calls
732 !
733 INTEGER growseason
734 EXTERNAL growseason
735
736 ! gday = growseason(julday)
737 gday = 91
738
739 IF (gday .EQ. 0) THEN
740 fertilizer_adj = 0.0
741 ELSEIF ((gday .GE. 1) .AND. (gday .LT. 30)) THEN ! first month of growing season
742 fertilizer_adj = 1.0
743 ELSEIF (gday .GE. 30) THEN
744 fertilizer_adj = 1.0+30.0/184.0-float(gday)/184.0
745 ELSE
746 write (*,*) 'ERROR: invalid Julian day'
747 write (*,*) 'julday = ', julday
748 write (*,*) 'growing season day = ',gday
749 CALL wrf_error_fatal ( 'INVALID GROWING SEASON DAY')
750 ENDIF
751
752 RETURN
753
754 END FUNCTION fertilizer_adj
755
756
757 REAL FUNCTION veg_adj(julday)
758 !*****************************************************************
759 !
760 ! SUMMARY:
761 ! computes vegetation adjustment factor from Julian day
762 !
763 ! FUNCTION CALLS:
764 ! growseason computes day of growing season
765 !
766 ! NOTE: julday = Julian day format
767 !
768 !*****************************************************************
769 implicit none
770
771 integer julday
772
773
774 !
775 !******** locals
776 !
777 integer gday
778
779 !
780 !******* function calls
781 !
782 INTEGER growseason
783 EXTERNAL growseason
784
785 !gday = growseason(julday)
786 gday = 91
787
788 IF (gday .LE. 30) THEN
789 veg_adj = 1.0
790 ELSEIF ((gday .GT. 30) .AND. (gday .LT. 60)) THEN
791 veg_adj = 1.5-(float(gday)/60.0)
792 ELSEIF (gday .GE. 60) THEN
793 veg_adj = 0.5
794 ELSE
795 write (*,*) 'ERROR: invalid Julian day'
796 write (*,*) 'julday = ', julday
797 write (*,*) 'growing season day = ',gday
798 CALL wrf_error_fatal ( 'veg_adj: INVALID GROWING SEASON DAY' )
799 ENDIF
800
801
802 RETURN
803
804
805 END FUNCTION veg_adj
806
807 END SUBROUTINE hrno
808
809 INTEGER FUNCTION growseason(julday)
810 !*****************************************************************
811 !
812 ! SUMMARY:
813 ! computes day of growing season from Julian day
814 !
815 ! NOTE: julday = Julian day format
816 !
817 !*****************************************************************
818 implicit none
819 integer julday
820
821 !*******
822 !
823 !
824 ! given Julian day, compute day of growing season
825 !
826 !
827 !
828 !******** locals
829
830 integer gsjulian_start
831 integer gsjulian_end
832
833 data gsjulian_start /91/ !=April 1 in non-leap-year
834 data gsjulian_end /304/ !=Oct 31 in non-leap-year
835
836 IF ((julday .GE. gsjulian_start) &
837 .AND. (julday .LE. gsjulian_end)) THEN ! growing season
838
839 growseason = julday-gsjulian_start+1
840
841
842 ELSEIF ((julday .GE. 1) & ! before or after growing season
843 .AND. (julday .LE. 366)) THEN
844
845 growseason = 0
846
847 ELSE
848 write (*,*) 'ERROR: Invalid julday '
849 write (*,*) 'julday = ',julday
850 CALL wrf_error_fatal ( 'growseason: INVALID JULIAN DAY')
851 ENDIF
852
853
854 RETURN
855 END FUNCTION growseason
856
857
858 END MODULE module_bioemi_beis311