chem/module_bioemi_beis314.F

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