1 MODULE module_bioemi_beis311
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
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 )
32 USE module_state_description
37 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
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
47 REAL, INTENT (IN) :: gmt,dtstep
49 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
54 REAL, DIMENSION( ims:ime , jms:jme ), &
56 xlat, & !latitude (deg)
57 xlong, & !longitude (deg)
58 gsw !downward shortwave surface flux (W/m^2)
60 ! Normalized biogenic emissions for standard conditions (moles compound/km^2/hr)
61 REAL, DIMENSION( ims:ime , jms:jme ), &
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
68 ! Leaf area index for isoprene
69 REAL, DIMENSION( ims:ime , jms:jme ), &
72 ! Actual biogenic emissions (moles compound/km^2/hr)
73 REAL, DIMENSION( ims:ime , jms:jme ), &
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
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
98 ! actual emissions for NO
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
114 ! Error message variables
115 INTEGER , PARAMETER :: ldev = 6 ! unit number for log file
116 CHARACTER*256 :: mesg
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
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
131 !***************************************
132 ! begin body of subroutine bio_emissions_beis311
134 ! hour into integration
135 gmtp=float(ktau)*dtstep/3600.
137 gmtp=mod(gmt+gmtp,24.)
138 write(mesg,*) 'calculate beis311 emissions at gmtp = ',gmtp
139 call wrf_debug(15,mesg)
143 tair = t_phy(i,kts,j)
144 pres = .01*p_phy(i,kts,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)
168 !....Perform checks on max and min bounds for temperature
170 IF (tair .LT. 200.0) THEN
171 ! WRITE( mesg, 94010 )
173 ! & 'too low at i,j= ',i,',',j
174 WRITE( ldev, * ) mesg
177 IF (tair .GT. 315.0 ) THEN
178 ! WRITE( mesg, 94020 )
180 ! & 'too high at i,j= ',i,',',j,
181 ! & '...resetting to 315K'
183 ! WRITE( ldev, * ) mesg
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 ) )
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)
197 !...... Convert tsolar to PAR and find direct and diffuse fractions
198 CALL getpar( tsolar, pres, zen, pardb, pardif )
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 )
206 ! & 'out of range at i,j= ',i,',',j
207 ! WRITE( ldev, * ) mesg
210 !...... Check max bound of LAI
211 IF ( tlai .GT. 10.0 ) THEN
212 ! WRITE( mesg, 94010 )
214 ! & 'out of range at i,j= ',i,',',j
215 ! WRITE( ldev, * ) mesg
218 !...... Initialize csubl
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
227 !...... Calculate csubl including shading if LAI > 0.1
228 IF ( tlai .GT. 0.1 ) THEN
229 csubl = clnew( zen, pardb, pardif, tlai )
231 !...... Otherwise calculate csubl without considering LAI
232 ELSE ! keep this or not?
237 ebio_iso(i,j) = se_iso * ct * csubl
242 !... Other biogenic emissions except NO:
243 !...... RACM: oli, api, lim, hc3, ete, olt, ket, ald, hcho, eth, ora2, co
245 cfovoc = EXP( 0.09 * ( tair - 303.0 ) )
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
264 CALL hrno( julday, growagno, ngrowagno, nonagno, tair, e_no)
273 !****************** FORMAT STATEMENTS ******************************
275 !........... Informational (LOG) message formats... 92xxx
278 !........... Internal buffering formats............ 94xxx
281 94010 FORMAT( A, F10.2, 1X, A, I4, A1, I4)
282 94020 FORMAT( A, F10.2, 1X, A, I4, A1, I4, 1X, A)
285 !***************** CONTAINS ********************************************
288 REAL FUNCTION clnew( zen, pardb, pardif, tlai )
290 !........ Function to calculate csubl based on zenith angle, PAR, and LAI
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
307 !........... CN98 - eqn 15.4, assume x=1
309 kbe = 0.5 * SQRT(1. + TAN( zen ) * TAN( zen ))
311 !.......... CN98 - p. 261 (this is usually small)
313 canparscat = 0.5 * pardb * (EXP(-0.894 * kbe * tlai) - &
314 EXP(-1.* kbe * tlai))
316 !.......... CN98 - p. 261 (assume exponentially wtd avg)
318 canpardif = pardif * (1. - EXP(-0.61 * tlai))/(0.61 * tlai)
320 !......... CN98 - p. 261 (for next 3 eqns)
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
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.
334 clnew =fracsun * cguen(parsun) + fracshade * cguen(parshade)
339 REAL FUNCTION cguen( partmp )
341 !.......... Guenther's equation for computing light correction
344 REAL, INTENT (IN) :: partmp
345 REAL, PARAMETER :: alpha2 = 0.00000729
347 IF ( partmp .LE. 0.01) THEN
350 cguen = (0.0028782 * partmp) / &
351 SQRT(1. + alpha2 * partmp * partmp)
357 END SUBROUTINE bio_emissions_beis311
359 !=================================================================
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:
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
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
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, &
383 CHARACTER*256 :: mesg
384 ! .. Intrinsic Functions ..
385 INTRINSIC acos, atan, cos, min, sin, tan
397 ! calc geom mean longitude
398 ml = 279.2801988 + .9856473354*d + 2.267E-13*d*d
401 ! calc equation of time in sec
402 ! w = mean long of perigee
404 ! epsi = mean obliquity of ecliptic
405 w = 282.4932328 + 4.70684E-5*d + 3.39E-13*d*d
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
410 yt = (tan(pepsi/2.0))**2
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
425 ! convert eq of time from sec to deg
427 ! calc right ascension in rads
430 ! calc declination in rads, deg
431 tab = 0.43360*sin(rra)
434 ! calc local hour angle
435 lbgmt = 12.0 - eqt/3600. + long*24./360.
436 lzgmt = 15.0*(gmt-lbgmt)
438 csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
440 write(mesg,*) 'calczen,csz ',csz
441 call wrf_debug(15,mesg)
446 ! keep zenith angle in radians for later use (GJF 6/2004)
451 END SUBROUTINE calc_zenithb
453 !=================================================================
456 SUBROUTINE getpar( tsolar, pres, zen, pardb, pardif )
458 !***********************************************************************
459 ! subroutine body starts at line
463 ! Based on code from Bart Brashers (10/2000), which was based on
464 ! code from Weiss and Norman (1985).
467 ! PRECONDITIONS REQUIRED:
468 ! Solar radiation (W/m2) and pressure (mb)
470 ! SUBROUTINES AND FUNCTIONS CALLED:
473 ! 3/01 : Prototype by JMV
475 !***********************************************************************
477 ! Project Title: Sparse Matrix Operator Kernel Emissions (SMOKE) Modeling
479 ! File: @(#)Id: getpar.f,v 1.1.1.1 2001/03/27 19:08:49 smith_w Exp
481 ! COPYRIGHT (C) 2001, MCNC--North Carolina Supercomputing Center
482 ! All Rights Reserved
484 ! See file COPYRIGHT for conditions of use.
486 ! MCNC-Environmental Programs Group
488 ! Research Triangle Park, NC 27709-2889
492 ! Pathname: Source: /env/proj/archive/cvs/jmv/beis3v0.9/getpar.f,v
493 ! Last updated: Date: 2001/03/27 19:08:49
495 !***********************************************************************
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)
507 REAL, INTENT (OUT) :: pardb ! direct beam PAR( umol/m2-s)
508 REAL, INTENT (OUT) :: pardif ! diffuse PAR ( umol/m2-s)
510 !........... PARAMETERS and their descriptions:
512 REAL, PARAMETER :: watt2umol = 4.6 ! convert W/m^2 to umol/m^2-s (4.6)
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
528 !***************************************
529 ! begin body of subroutine
531 !............ Assume that PAR = 0 if zenith angle is greater than 87 degrees
532 !............ or if solar radiation is zero
534 IF (zen .GE. 1.51844 .OR. tsolar .LE. 0.) THEN
540 !............ Compute clear sky (aka potential) radiation terms
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)
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
554 !............ Compute fraction of visible that is direct beam
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
561 fvb = rdvis/rvt * (1.-((0.9 - ratio)/0.7)**0.666667)
565 !............ Compute PAR (direct beam and diffuse) in umol/m2-sec
567 pardb = tsolar * fvis * fvb * watt2umol
568 pardif = tsolar * fvis * fvd * watt2umol
573 !****************** FORMAT STATEMENTS ******************************
575 !........... Informational (LOG) message formats... 92xxx
578 !........... Internal buffering formats............ 94xxx
580 END SUBROUTINE getpar
582 SUBROUTINE hrno( julday, growagno, ngrowagno, nonagno, tairin, e_no)
584 !***********************************************************************
585 ! subroutine body starts at line 150
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
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
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
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
613 ! PRECONDITIONS REQUIRED:
614 ! Normalized NO emissions, Surface Temperature
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
623 ! 10/01 : Prototype by GAP
625 !***********************************************************************
627 ! Project Title: BEIS3 Enhancements for NO emission calculation
631 !***********************************************************************
635 !........... ARGUMENTS and their descriptions:
637 INTEGER, INTENT (IN) :: julday ! current julian day
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
645 REAL, INTENT (OUT) :: e_no ! output NO emissions
647 !........... SCRATCH LOCAL VARIABLES and their descriptions:
649 REAL cfno ! NO correction factor
650 REAL cfnograss ! NO correction factor for grasslands
651 REAL tsoi ! soil temperature
652 REAL tair ! air temperature
654 REAL :: cfnowet, cfnodry
656 INTEGER growseason, gday
658 !***********************************************************************
662 !............. calculate NO emissions by going thru temperature cases
664 ! gday = growseason(julday)
666 IF (gday .eq. 0) THEN !not growing season
667 IF ( tair .GT. 303.00 ) tair = 303.00
669 IF ( tair .GT. 268.8690 ) THEN
670 cfno = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
676 ngrowagno * cfno & !agriculture
677 + nonagno * cfno ! non-agriculture
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
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)) &
692 cfno = 0.5*cfnowet + 0.5*cfnodry
694 IF ( tair .GT. 303.00 ) tair = 303.00
696 IF ( tair .GT. 268.8690 ) THEN
697 cfnograss = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
702 e_no = growagno * cfno *fertilizer_adj(julday)*veg_adj(julday) &
703 + nonagno * cfnograss
709 !***************** CONTAINS ********************************************
712 REAL FUNCTION fertilizer_adj(julday)
713 !*****************************************************************
716 ! computes fertilizer adjustment factor from Julian day
719 ! growseason computes day of growing season
721 ! NOTE: julday = Julian day format
723 !*****************************************************************
727 !******** local scratch variables
731 !******** function calls
736 ! gday = growseason(julday)
739 IF (gday .EQ. 0) THEN
741 ELSEIF ((gday .GE. 1) .AND. (gday .LT. 30)) THEN ! first month of growing season
743 ELSEIF (gday .GE. 30) THEN
744 fertilizer_adj = 1.0+30.0/184.0-float(gday)/184.0
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')
754 END FUNCTION fertilizer_adj
757 REAL FUNCTION veg_adj(julday)
758 !*****************************************************************
761 ! computes vegetation adjustment factor from Julian day
764 ! growseason computes day of growing season
766 ! NOTE: julday = Julian day format
768 !*****************************************************************
780 !******* function calls
785 !gday = growseason(julday)
788 IF (gday .LE. 30) THEN
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
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' )
809 INTEGER FUNCTION growseason(julday)
810 !*****************************************************************
813 ! computes day of growing season from Julian day
815 ! NOTE: julday = Julian day format
817 !*****************************************************************
824 ! given Julian day, compute day of growing season
830 integer gsjulian_start
833 data gsjulian_start /91/ !=April 1 in non-leap-year
834 data gsjulian_end /304/ !=Oct 31 in non-leap-year
836 IF ((julday .GE. gsjulian_start) &
837 .AND. (julday .LE. gsjulian_end)) THEN ! growing season
839 growseason = julday-gsjulian_start+1
842 ELSEIF ((julday .GE. 1) & ! before or after growing season
843 .AND. (julday .LE. 366)) THEN
848 write (*,*) 'ERROR: Invalid julday '
849 write (*,*) 'julday = ',julday
850 CALL wrf_error_fatal ( 'growseason: INVALID JULIAN DAY')
855 END FUNCTION growseason
858 END MODULE module_bioemi_beis311