wrfv2_fire/chem/module_bioemi_beis311.F

   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