MODULE module_bioemi_beis314
!
! BEIS3.14 Modified to include sesquiterpenes and MBO
! BEIS3.13 Emissions Module for WRF-Chem
! BEIS3.11 Written by Greg Frost 6/2004
!
! Modified to BEIS3.14 by Stu McKeen 11/9/09 - addition of sesquiterpene total emissions
! Modified to BEIS3.13 by Stu McKeen 9/27/07
! Note: species to vegetation assignment table (beis3_efac_v0.99_053105.txt) must be
! applied to front end for BEIS3.13 reference emission files - done outside of WRF-Chem
! BEIS3.13 also recommends 10meter or 2m air temp. for isoprene emissions, (no guidance)
! See: Changes to the Biogenic Emissions Inventory System Version 3 (BEIS3)
! Donna Schwede, George Pouliot, and Tom Pierce, presented at 2005 CMAS conference:
! available on the web at http://www/cmascenter.org/conference/2005/archive.cfm
!
! Using off-line gridded standard biogenic emissions
! for each model compound with such emissions,
! model shortwave solar flux (isoprene only),
! & air temperature, pressure, and density in lowest model level,
! calculates actual biogenic emissions of each compound.
! Based on hrbeis311.f from BEIS3.11 for SMOKE, with major
! surgery performed on original routines for use with WRF-Chem.
! This version assumes chemical mechanism is RACM.
! The following 16 RACM compounds have biogenic emissions:
! iso, no, oli, api, lim, xyl, hc3, ete, olt, ket, ald, hcho, eth, ora2, co, nr
!23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
CONTAINS
SUBROUTINE bio_emissions_beis314(id,config_flags,ktau,curr_secs, &
dtstep,julday,gmt,xlat,xlong,t_phy,p_phy,gsw, &
sebio_iso,sebio_oli,sebio_api,sebio_lim,sebio_xyl, &
sebio_hc3,sebio_ete,sebio_olt,sebio_ket,sebio_ald, &
sebio_hcho,sebio_eth,sebio_ora2,sebio_co,sebio_nr, &
sebio_sesq,sebio_mbo, &
noag_grow,noag_nongrow,nononag,slai, &
ebio_iso,ebio_oli,ebio_api,ebio_lim,ebio_xyl, &
ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_nr,ebio_no, &
ebio_sesq,ebio_mbo, & ! sesq and mbo are added
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
USE module_configure
USE module_state_description
IMPLICIT NONE
! .. Parameters ..
TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
! .. Indices ..
INTEGER, INTENT(IN ) :: id, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
! .. Passed variables ..
INTEGER, INTENT (IN) :: ktau, & ! time step number
julday ! current simulation julian day
REAL(KIND=8), INTENT(IN) :: curr_secs ! seconds into the simulation
REAL, INTENT (IN) :: gmt,dtstep
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
INTENT(IN ) :: &
t_phy, & !air T (K)
p_phy !P (Pa)
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(IN ) :: &
xlat, & !latitude (deg)
xlong, & !longitude (deg)
gsw !downward shortwave surface flux (W/m^2)
! Normalized biogenic emissions for standard conditions (moles compound/km^2/hr)
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(IN ) :: &
sebio_iso,sebio_oli,sebio_api,sebio_lim,sebio_xyl, &
sebio_hc3,sebio_ete,sebio_olt,sebio_ket,sebio_ald, &
sebio_hcho,sebio_eth,sebio_ora2,sebio_co,sebio_nr, &
sebio_sesq,sebio_mbo, &
noag_grow,noag_nongrow,nononag
! Leaf area index for isoprene
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(IN ) :: slai
! Actual biogenic emissions (moles compound/km^2/hr)
REAL, DIMENSION( ims:ime , jms:jme ), &
INTENT(INOUT ) :: &
ebio_iso,ebio_oli,ebio_api,ebio_lim,ebio_xyl, &
ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_nr,ebio_no, &
ebio_sesq,ebio_mbo
! .. Local Scalars ..
INTEGER :: i,j
! Variables for 1 element of I/O arrays
! met and phys input variables
REAL :: tair ! surface air temperature (K)
REAL :: tsolar ! downward shortwave surface flux (W/m^2)
REAL :: pres ! surface pressure (mb)
REAL :: ylat ! latitude (deg)
REAL :: ylong ! longitude (deg)
! normalized emissions (moles compound/km^2/hr)
REAL :: se_iso,se_oli,se_api,se_lim,se_xyl, &
se_hc3,se_ete,se_olt,se_ket,se_ald, &
se_hcho,se_eth,se_ora2,se_co,se_nr, &
se_mbo,se_sesq, &
growagno,ngrowagno,nonagno
! leaf area index for isoprene
REAL :: tlai
! actual emissions for NO
REAL :: e_no
! Other parameters needed in calculations
! Guenther's parameterizations: Guenther et al. JGR 98, 12609-12617, 1993
REAL :: ct, dt ! Guenther's temperature correction for isoprene
REAL :: cfno ! NO correction factor
REAL :: cfovoc ! non-isoprene correction factor
REAL :: par ! PAR = photosynthetically active radiation (micromole/m^2/s)
REAL :: csubl ! C sub l from Guenther
REAL :: zen ! zenith angle (radians)
REAL :: coszen ! cosine(zenith angle)
REAL :: pardb ! PAR direct beam
REAL :: pardif ! PAR diffuse
REAL :: gmtp ! current simulation time
! Error message variables
INTEGER , PARAMETER :: ldev = 6 ! unit number for log file
CHARACTER*256 :: mesg
! Functions called directly or indirectly
! clnew calculates csubl based on zenith angle, par, and lai
! cguen Guenther's equation for computing light correction
! fertilizer_adj computes fertlizer adjustment factor
! veg_adj computes vegatation adjustment factor
! growseason computes day of growing season
! Subroutines called directly or indirectly
! calc_zenithb calculates zenith angle from latitude, longitude, julian day, and GMT
! NOTE: longitude input for this routine is nonstandard: >0 for W, <0 for E!!
! getpar computes PAR (direct beam and diffuse) in umol/m2-sec from downward shortwave flux
! hrno algorithm to estimate NO emissions; does not include precipitation adjustment
!***************************************
! begin body of subroutine bio_emissions_beis314
! hour into integration
!The old gmtp method will break for runs longer than about 12 days with r4
! gmtp=float(ktau)*dtstep/3600.
gmtp = curr_secs/3600._8
!
gmtp=mod(gmt+gmtp,24.)
write(mesg,*) 'calculate beis314 emissions at gmtp = ',gmtp
call wrf_debug(15,mesg)
DO 100 j=jts,jte
DO 100 i=its,ite
tair = t_phy(i,kts,j)
pres = .01*p_phy(i,kts,j)
ylat = xlat(i,j)
ylong = xlong(i,j)
tsolar = gsw(i,j)
tlai = slai(i,j)
se_iso = sebio_iso(i,j)
se_oli = sebio_oli(i,j)
se_api = sebio_api(i,j)
se_lim = sebio_lim(i,j)
se_xyl = sebio_xyl(i,j)
se_hc3 = sebio_hc3(i,j)
se_ete = sebio_ete(i,j)
se_olt = sebio_olt(i,j)
se_ket = sebio_ket(i,j)
se_ald = sebio_ald(i,j)
se_hcho = sebio_hcho(i,j)
se_eth = sebio_eth(i,j)
se_ora2 = sebio_ora2(i,j)
se_co = sebio_co(i,j)
se_nr = sebio_nr(i,j)
! SESQ and MBO are added
se_sesq = sebio_sesq(i,j)
se_mbo = sebio_mbo(i,j)
growagno = noag_grow(i,j)
ngrowagno = noag_nongrow(i,j)
nonagno = nononag(i,j)
!....Perform checks on max and min bounds for temperature
IF (tair .LT. 200.0) THEN
! WRITE( mesg, 94010 )
! & 'tair=', tair,
! & 'too low at i,j= ',i,',',j
WRITE( ldev, * ) mesg
END IF
IF (tair .GT. 315.0 ) THEN
! WRITE( mesg, 94020 )
! & 'tair=', tair,
! & 'too high at i,j= ',i,',',j,
! & '...resetting to 315K'
tair = 315.0
! WRITE( ldev, * ) mesg
ENDIF
!... Isoprene emissions
!...... Calculate temperature correction term
dt = 28668.514 / tair
ct = EXP( 37.711 - 0.398570815 * dt ) / &
(1.0 + EXP( 91.301 - dt ) )
!...... Calculate zenith angle in radians
! NOTE: nonstandard longitude input here: >0 for W, <0 for E!!
CALL calc_zenithb(ylat,-ylong,julday,gmtp,zen)
coszen = COS(zen)
!...... Convert tsolar to PAR and find direct and diffuse fractions
CALL getpar( tsolar, pres, zen, pardb, pardif )
par = pardb + pardif
!...... Check max/min bounds of PAR and calculate
!...... biogenic isoprene
IF ( par .LT. 0.00 .OR. par .GT. 2600.0 ) THEN
! WRITE( mesg, 94010 )
! & 'PAR=', par,
! & 'out of range at i,j= ',i,',',j
! WRITE( ldev, * ) mesg
ENDIF
!...... Check max bound of LAI
IF ( tlai .GT. 10.0 ) THEN
! WRITE( mesg, 94010 )
! & 'LAI=', tlai,
! & 'out of range at i,j= ',i,',',j
! WRITE( ldev, * ) mesg
ENDIF
!...... Initialize csubl
csubl = 0.0
!...... If PAR < 0.01 or zenith angle > 89 deg, set isoprene emissions to 0.
IF ( par .LE. 0.01 .OR. coszen .LE. 0.02079483 ) THEN
ebio_iso(i,j) = 0.0
ELSE
!...... Calculate csubl including shading if LAI > 0.1
IF ( tlai .GT. 0.1 ) THEN
csubl = clnew( zen, pardb, pardif, tlai )
!...... Otherwise calculate csubl without considering LAI
ELSE ! keep this or not?
csubl = cguen( par )
ENDIF
ebio_iso(i,j) = se_iso * ct * csubl
ENDIF
!... Other biogenic emissions except NO:
!...... RACM: oli, api, lim, hc3, ete, olt, ket, ald, hcho, eth, ora2, co
cfovoc = EXP( 0.09 * ( tair - 303.0 ) )
ebio_oli(i,j) = se_oli * cfovoc
ebio_api(i,j) = se_api * cfovoc
ebio_lim(i,j) = se_lim * cfovoc
ebio_xyl(i,j) = se_xyl * cfovoc
ebio_hc3(i,j) = se_hc3 * cfovoc
ebio_ete(i,j) = se_ete * cfovoc
ebio_olt(i,j) = se_olt * cfovoc
ebio_ket(i,j) = se_ket * cfovoc
ebio_ald(i,j) = se_ald * cfovoc
ebio_hcho(i,j) = se_hcho * cfovoc
ebio_eth(i,j) = se_eth * cfovoc
ebio_ora2(i,j) = se_ora2 * cfovoc
ebio_co(i,j) = se_co * cfovoc
ebio_nr(i,j) = se_nr * cfovoc
! SESQ and MBO are added
ebio_sesq(i,j) = se_sesq * cfovoc
ebio_mbo(i,j) = se_mbo * cfovoc
!... NO emissions
CALL hrno( julday, growagno, ngrowagno, nonagno, tair, e_no)
ebio_no(i,j) = e_no
100 CONTINUE
RETURN
!****************** FORMAT STATEMENTS ******************************
!........... Informational (LOG) message formats... 92xxx
!........... Internal buffering formats............ 94xxx
94010 FORMAT( A, F10.2, 1X, A, I4, A1, I4)
94020 FORMAT( A, F10.2, 1X, A, I4, A1, I4, 1X, A)
!***************** CONTAINS ********************************************
CONTAINS
REAL FUNCTION clnew( zen, pardb, pardif, tlai )
!........ Function to calculate csubl based on zenith angle, PAR, and LAI
!******** Reference:CN98
! Campbell, G.S. and J.M. Norman. 1998. An Introduction to Environmental Biophysics,
! Springer-Verlag, New York.
IMPLICIT NONE
REAL, INTENT (IN) :: pardb ! direct beam PAR( umol/m2-s)
REAL, INTENT (IN) :: pardif ! diffuse PAR ( umol/m2-s)
REAL, INTENT (IN) :: zen ! solar zenith angle (radians)
REAL, INTENT (IN) :: tlai ! leaf area index for grid cell
!............. Local variables
REAL kbe ! extinction coefficient for direct beam
REAL ALPHA ! leave absorptivity
REAL KD ! extinction coefficient for diffuse radiation
REAL SQALPHA ! square root of alpha
REAL canparscat ! exponentially wtd scattered PAR (umol/m2-s)
REAL canpardif ! exponentially wtd diffuse PAR (umol/m2-s)
REAL parshade ! PAR on shaded leaves (umol/m2-s)
REAL parsun ! PAR on sunlit leaves (umol/m2-s)
REAL laisun ! LAI that is sunlit
REAL fracsun ! fraction of leaves that are sunlit
REAL fracshade ! fraction of leaves that are shaded
ALPHA = 0.8
SQALPHA = SQRT(0.8)
KD = 0.68
!........... CN98 - eqn 15.4, assume x=1
kbe = 0.5 * SQRT(1. + TAN( zen ) * TAN( zen ))
!.......... CN98 - p. 261 (this is usually small)
canparscat = 0.5 * pardb * (EXP(-1.*sqalpha*kbe*tlai) - &
EXP(-1.* kbe * tlai))
!.......... CN98 - p. 261 (assume exponentially wtd avg)
canpardif = pardif * (1.-EXP(-1.*sqalpha*kd*tlai)) / &
(sqalpha*kd*tlai)
!......... CN98 - p. 261 (for next 3 eqns)
parshade = canpardif + canparscat
parsun = kbe * pardb + parshade
laisun = (1. - EXP( -1. * kbe * tlai))/kbe
fracsun = laisun/tlai
fracshade = 1. - fracsun
!.......... cguen is guenther's eqn for computing light correction as a
!.......... function of PAR...fracSun should probably be higher since
!.......... sunlit leaves tend to be thicker than shaded leaves. But
!.......... since we need to make crude asmptns regarding leave
!.......... orientation (x=1), will not attempt to fix at the moment.
clnew =fracsun * cguen(parsun) + fracshade * cguen(parshade)
RETURN
END FUNCTION clnew
REAL FUNCTION cguen( partmp )
!.......... Guenther's equation for computing light correction
! Reference: Guenther, A., B. Baugh, G. Brasseur, J. Greenberg, P. Harley, L. Klinger,
! D. Serca, and L. Vierling, 1999: Isoprene emission estimates and uncertainties
! for the Central African EXPRESSO Study domain. J. Geophys. Res., 104, 30625-30639.
IMPLICIT NONE
REAL, INTENT (IN) :: partmp
REAL, PARAMETER :: alpha = 0.001
REAL, PARAMETER :: cl = 1.42
IF ( partmp .LE. 0.01) THEN
cguen = 0.0
ELSE
cguen = (alpha *cl * partmp) / &
SQRT(1. + alpha * alpha * partmp * partmp)
ENDIF
RETURN
END FUNCTION cguen
END SUBROUTINE bio_emissions_beis314
!=================================================================
SUBROUTINE calc_zenithb(lat,long,ijd,gmt,zenith)
! Based on calc_zenith from WRF-Chem module_phot_mad.F
! this subroutine calculates solar zenith angle for a
! time and location. must specify:
! input:
! lat - latitude in decimal degrees
! long - longitude in decimal degrees
! NOTE: Nonstandard convention for long: >0 for W, <0 for E!!
! gmt - greenwich mean time - decimal military eg.
! 22.75 = 45 min after ten pm gmt
! output
! zenith - in radians (GJF, 6/2004)
! remove azimuth angle calculation since not needed (GJF, 6/2004)
! .. Scalar Arguments ..
CHARACTER*256 :: mesg
REAL :: gmt, lat, long, zenith
INTEGER :: ijd
! .. Local Scalars ..
REAL :: csz, cw, d, decl, dr, ec, epsi, eqt, eyt, feqt, feqt1, &
feqt2, feqt3, feqt4, feqt5, feqt6, feqt7, lbgmt, lzgmt, ml, pepsi, &
pi, ra, rdecl, reqt, rlt, rml, rphi, rra, ssw, sw, tab, w, wr, &
yt, zpt, zr
INTEGER :: jd
! .. Intrinsic Functions ..
INTRINSIC acos, atan, cos, min, sin, tan
! convert to radians
pi = 3.1415926535590
dr = pi/180.
rlt = lat*dr
rphi = long*dr
! ???? + (yr - yref)
jd = ijd
d = jd + gmt/24.0
! calc geom mean longitude
ml = 279.2801988 + .9856473354*d + 2.267E-13*d*d
rml = ml*dr
! calc equation of time in sec
! w = mean long of perigee
! e = eccentricity
! epsi = mean obliquity of ecliptic
w = 282.4932328 + 4.70684E-5*d + 3.39E-13*d*d
wr = w*dr
ec = 1.6720041E-2 - 1.1444E-9*d - 9.4E-17*d*d
epsi = 23.44266511 - 3.5626E-7*d - 1.23E-15*d*d
pepsi = epsi*dr
yt = (tan(pepsi/2.0))**2
cw = cos(wr)
sw = sin(wr)
ssw = sin(2.0*wr)
eyt = 2.*ec*yt
feqt1 = sin(rml)*(-eyt*cw-2.*ec*cw)
feqt2 = cos(rml)*(2.*ec*sw-eyt*sw)
feqt3 = sin(2.*rml)*(yt-(5.*ec**2/4.)*(cw**2-sw**2))
feqt4 = cos(2.*rml)*(5.*ec**2*ssw/4.)
feqt5 = sin(3.*rml)*(eyt*cw)
feqt6 = cos(3.*rml)*(-eyt*sw)
feqt7 = -sin(4.*rml)*(.5*yt**2)
feqt = feqt1 + feqt2 + feqt3 + feqt4 + feqt5 + feqt6 + feqt7
eqt = feqt*13751.0
! convert eq of time from sec to deg
reqt = eqt/240.
! calc right ascension in rads
ra = ml - reqt
rra = ra*dr
! calc declination in rads, deg
tab = 0.43360*sin(rra)
rdecl = atan(tab)
decl = rdecl/dr
! calc local hour angle
lbgmt = 12.0 - eqt/3600. + long*24./360.
lzgmt = 15.0*(gmt-lbgmt)
zpt = lzgmt*dr
csz = sin(rlt)*sin(rdecl) + cos(rlt)*cos(rdecl)*cos(zpt)
if(csz.gt.1) then
write(mesg,*) 'calczen,csz ',csz
call wrf_debug(15,mesg)
endif
csz = min(1.,csz)
zr = acos(csz)
! zenith = zr/dr
! keep zenith angle in radians for later use (GJF 6/2004)
zenith = zr
RETURN
END SUBROUTINE calc_zenithb
!=================================================================
SUBROUTINE getpar( tsolar, pres, zen, pardb, pardif )
!***********************************************************************
! subroutine body starts at line
!
! DESCRIPTION:
!
! Based on code from Bart Brashers (10/2000), which was based on
! code from Weiss and Norman (1985).
!
!
! PRECONDITIONS REQUIRED:
! Solar radiation (W/m2) and pressure (mb)
!
! SUBROUTINES AND FUNCTIONS CALLED:
!
! REVISION HISTORY:
! 3/01 : Prototype by JMV
!
!***********************************************************************
!
! Project Title: Sparse Matrix Operator Kernel Emissions (SMOKE) Modeling
! System
! File: @(#)Id: getpar.f,v 1.1.1.1 2001/03/27 19:08:49 smith_w Exp
!
! COPYRIGHT (C) 2001, MCNC--North Carolina Supercomputing Center
! All Rights Reserved
!
! See file COPYRIGHT for conditions of use.
!
! MCNC-Environmental Programs Group
! P.O. Box 12889
! Research Triangle Park, NC 27709-2889
!
! env_progs@mcnc.org
!
! Pathname: Source: /env/proj/archive/cvs/jmv/beis3v0.9/getpar.f,v
! Last updated: Date: 2001/03/27 19:08:49
!
!***********************************************************************
IMPLICIT NONE
!........ Inputs
REAL , INTENT (IN) :: tsolar ! modeled or observed total radiation (W/m2)
REAL , INTENT (IN) :: pres ! atmospheric pressure (mb)
REAL , INTENT (IN) :: zen ! solar zenith angle (radians)
!........ Outputs
REAL, INTENT (OUT) :: pardb ! direct beam PAR( umol/m2-s)
REAL, INTENT (OUT) :: pardif ! diffuse PAR ( umol/m2-s)
!........... PARAMETERS and their descriptions:
REAL, PARAMETER :: watt2umol = 4.6 ! convert W/m^2 to umol/m^2-s (4.6)
!
REAL ratio ! transmission fraction for total radiation
REAL ot ! optical thickness
REAL rdvis ! possible direct visible beam (W/m^2)
REAL rfvis ! possible visible diffuse (W/m^2)
REAL wa ! water absorption in near-IR (W/m^2)
REAL rdir ! direct beam in near-IR (W/m^2)
REAL rfir ! diffuse near-IR (W/m^2)
REAL rvt ! total possible visible radiation (W/m^2)
REAL rirt ! total possible near-IR radiation (W/m^2)
REAL fvis ! fraction of visible to total
REAL fvb ! fraction of visible that is direct beam
REAL fvd ! fraction of visible that is diffuse
!***************************************
! begin body of subroutine
!............ Assume that PAR = 0 if zenith angle is greater than 87 degrees
!............ or if solar radiation is zero
IF (zen .GE. 1.51844 .OR. tsolar .LE. 0.) THEN
pardb = 0.
pardif = 0.
RETURN
ENDIF
!............ Compute clear sky (aka potential) radiation terms
ot = pres / 1013.25 / COS(zen) !Atmospheric Optical thickness
rdvis = 600. * EXP(-.185 * ot) * COS(zen) !Direct visible beam, eqn (1)
rfvis = 0.42 * (600 - rdvis) * COS(zen) !Visible Diffuse, eqn (3)
wa = 1320 * .077 * (2. * ot)**0.3 !water absorption in near-IR, eqn (6)
rdir = (720. * EXP(-0.06 * ot)-wa) * COS(zen) !Direct beam near-IR, eqn (4)
rfir = 0.65 * (720. - wa - rdir) * COS(zen) !Diffuse near-IR, eqn (5)
rvt = rdvis + rfvis !Total visible radiation, eqn (9)
rirt = rdir + rfir !Total near-IR radiation, eqn (10)
fvis = rvt/(rirt + rvt) !Fraction of visible to total radiation, eqn 7
ratio = tsolar /(rirt + rvt) !Ratio of "actual" to clear sky solar radiation
!............ Compute fraction of visible that is direct beam
IF (ratio .GE. 0.89) THEN
fvb = rdvis/rvt * 0.941124
ELSE IF (ratio .LE. 0.21) THEN
fvb = rdvis/rvt * 9.55E-3
ELSE
fvb = rdvis/rvt * (1.-((0.9 - ratio)/0.7)**0.666667)
ENDIF
fvd = 1. - fvb
!............ Compute PAR (direct beam and diffuse) in umol/m2-sec
pardb = tsolar * fvis * fvb * watt2umol
pardif = tsolar * fvis * fvd * watt2umol
RETURN
!****************** FORMAT STATEMENTS ******************************
!........... Informational (LOG) message formats... 92xxx
!........... Internal buffering formats............ 94xxx
END SUBROUTINE getpar
SUBROUTINE hrno( julday, growagno, ngrowagno, nonagno, tairin, e_no)
!***********************************************************************
! subroutine body starts at line 150
!
! DESCRIPTION:
!
! Uses new NO algorithm NO = Normalized*Tadj*Fadj*Cadj
! to estimate NO emissions
! Information needed to estimate NO emissions
! Julian Day (integer) julday
! Surface Temperature (MCIP field) tair (K)
! Note: Precipitation adjustment not used in the WRF-Chem implementation of BEIS3.11
! because of differences in soil categories between BEIS and WRF-Chem
!
! The calculation are based on the following paper by J.J. Yienger and H. Levy II
! J.J. Yienger and H. Levy II, Journal of Geophysical Research, vol 100,11447-11464,1995
!
! Also see the following paper for more information:
! Proceedings of the Air and Waste Management Association/U.S. Environmental Protection
! Agency EMission Inventory Conference, Raleigh October 26-28, 1999 Raleigh NC
! by Tom Pierce and Lucille Bender
!
! REFERENCES
!
! JACQUEMIN B. AND NOILHAN J. (1990), BOUND.-LAYER METEOROL., 52, 93-134.
! J.J. Yienger and H. Levy II, Journal of Geophysical Research, vol 100,11447-11464,1995
! T. Pierce and L. Bender, Examining the Temporal Variability of Ammonia and Nitric Oxide Emissions from Agricultural Processes
! Proceedings of the Air and Waste Management Association/U.S. Environmental Protection
! Agency EMission Inventory Conference, Raleigh October 26-28, 1999 Raleigh NC
!
! PRECONDITIONS REQUIRED:
! Normalized NO emissions, Surface Temperature
!
! SUBROUTINES AND FUNCTIONS CALLED (directly or indirectly):
! fertilizer_adj computes fertlizer adjustment factor
! veg_adj computes vegatation adjustment factor
! growseason computes day of growing season
!
!
! REVISION HISTORY:
! 10/01 : Prototype by GAP
!
!***********************************************************************
!
! Project Title: BEIS3 Enhancements for NO emission calculation
! File: hrno.f
!
!
!***********************************************************************
IMPLICIT NONE
!........... ARGUMENTS and their descriptions:
INTEGER, INTENT (IN) :: julday ! current julian day
REAL, INTENT (IN) :: tairin ! air temperature (K)
REAL, INTENT (IN) :: growagno ! norm NO emissions, agricultural, growing
REAL, INTENT (IN) :: ngrowagno ! norm NO emissions, agricultural, not growing
REAL, INTENT (IN) :: nonagno ! norm NO emissions, non-agricultural
REAL, INTENT (OUT) :: e_no ! output NO emissions
!........... SCRATCH LOCAL VARIABLES and their descriptions:
REAL cfno ! NO correction factor
REAL cfnograss ! NO correction factor for grasslands
REAL tsoi ! soil temperature
REAL tair ! air temperature
REAL :: cfnowet, cfnodry
INTEGER gday
!***********************************************************************
tair = tairin
!............. calculate NO emissions by going thru temperature cases
gday = growseason(julday)
! Control of growing season should be done in the input files for BEIS3.13
! gday = 91
IF (gday .eq. 0) THEN !not growing season
IF ( tair .GT. 303.00 ) tair = 303.00
IF ( tair .GT. 268.8690 ) THEN
cfno = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
ELSE
cfno = 0.0
ENDIF
e_no = &
ngrowagno * cfno & !agriculture
+ nonagno * cfno ! non-agriculture
ELSE
tsoi = 0.72*tair+82.28
IF (tsoi .LE. 273.16) tsoi = 273.16
IF (tsoi .GE. 303.16) tsoi = 303.16
cfnodry = (1./3.)*(1./30.)*(tsoi-273.16) ! see YL 1995 Equa 9a p. 11452
IF (tsoi .LE. 283.16) THEN ! linear cold case
cfnowet = (tsoi-273.16)*EXP(-0.103*30.0)*0.28 ! see YL 1995 Equ 7b
ELSE ! exponential case
cfnowet = EXP(0.103*(tsoi-273.16)) &
*EXP(-0.103*30.0)
ENDIF
cfno = 0.5*cfnowet + 0.5*cfnodry
IF ( tair .GT. 303.00 ) tair = 303.00
IF ( tair .GT. 268.8690 ) THEN
cfnograss = EXP( 0.04686 * tair - 14.30579 ) ! grass (from BEIS2)
ELSE
cfnograss = 0.0
ENDIF
e_no = growagno * cfno *fertilizer_adj(julday)*veg_adj(julday) &
+ nonagno * cfnograss
ENDIF
RETURN
!***************** CONTAINS ********************************************
CONTAINS
REAL FUNCTION fertilizer_adj(julday)
!*****************************************************************
!
! SUMMARY:
! computes fertilizer adjustment factor from Julian day
!
! FUNCTION CALLS:
! growseason computes day of growing season
!
! NOTE: julday = Julian day format
!
!*****************************************************************
implicit none
integer julday
!
!******** local scratch variables
!
integer gday
!
!******** function calls
!
gday = growseason(julday)
! Control of growing season should be done in the input files for BEIS3.13
! gday = 91
IF (gday .EQ. 0) THEN
fertilizer_adj = 0.0
ELSEIF ((gday .GE. 1) .AND. (gday .LT. 30)) THEN ! first month of growing season
fertilizer_adj = 1.0
ELSEIF (gday .GE. 30) THEN
fertilizer_adj = 1.0+30.0/184.0-float(gday)/184.0
ELSE
write (*,*) 'ERROR: invalid Julian day'
write (*,*) 'julday = ', julday
write (*,*) 'growing season day = ',gday
CALL wrf_error_fatal ( 'INVALID GROWING SEASON DAY')
ENDIF
RETURN
END FUNCTION fertilizer_adj
REAL FUNCTION veg_adj(julday)
!*****************************************************************
!
! SUMMARY:
! computes vegetation adjustment factor from Julian day
!
! FUNCTION CALLS:
! growseason computes day of growing season
!
! NOTE: julday = Julian day format
!
!*****************************************************************
implicit none
integer julday
!
!******** locals
!
integer gday
!
!******* function calls
!
gday = growseason(julday)
! Control of growing season should be done in the input files for BEIS3.13
!gday = 91
IF (gday .LE. 30) THEN
veg_adj = 1.0
ELSEIF ((gday .GT. 30) .AND. (gday .LT. 60)) THEN
veg_adj = 1.5-(float(gday)/60.0)
ELSEIF (gday .GE. 60) THEN
veg_adj = 0.5
ELSE
write (*,*) 'ERROR: invalid Julian day'
write (*,*) 'julday = ', julday
write (*,*) 'growing season day = ',gday
CALL wrf_error_fatal ( 'veg_adj: INVALID GROWING SEASON DAY' )
ENDIF
RETURN
END FUNCTION veg_adj
END SUBROUTINE hrno
INTEGER FUNCTION growseason(julday)
!*****************************************************************
!
! SUMMARY:
! computes day of growing season from Julian day
!
! NOTE: julday = Julian day format
!
!*****************************************************************
implicit none
integer julday
!*******
!
!
! given Julian day, compute day of growing season
!
!
!
!******** locals
integer gsjulian_start
integer gsjulian_end
data gsjulian_start /91/ !=April 1 in non-leap-year
data gsjulian_end /304/ !=Oct 31 in non-leap-year
IF ((julday .GE. gsjulian_start) &
.AND. (julday .LE. gsjulian_end)) THEN ! growing season
growseason = julday-gsjulian_start+1
ELSEIF ((julday .GE. 1) & ! before or after growing season
.AND. (julday .LE. 366)) THEN
growseason = 0
ELSE
write (*,*) 'ERROR: Invalid julday '
write (*,*) 'julday = ',julday
CALL wrf_error_fatal ( 'growseason: INVALID JULIAN DAY')
ENDIF
RETURN
END FUNCTION growseason
END MODULE module_bioemi_beis314