module_bioemi_megan2.F
References to this file elsewhere.
1 MODULE module_bioemi_megan2
2
3 ! MEGAN v2.04 Emissions Module for WRF-Chem
4 !
5 ! Reference:
6 !
7 ! Estimates of global terrestial isoprene emissions using MEGAN
8 ! (Model of Emissions of Gases and Aerosols from Nature )
9 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, P.I. Palmer, and C. Geron
10 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006.
11 !
12 ! MEGAN v2.0 Documentation
13 !
14 !
15 ! August 2007
16 !
17 ! Serena H. Chung Washington State University
18 ! Tan Sakulyanontvittaya University of Colorado
19 ! Christine Wiedinmyer National Center for Atmospheric Research
20 !
21 !
22 ! 11/08/2007 SHC Took out some "if (ktau ==1) then ... end if " statements
23 !
24
25 CONTAINS
26
27 SUBROUTINE bio_emissions_megan2(id,config_flags,ktau,dtstep, &
28 julday,gmt,xlat,xlong,p_phy,rho_phy,dz8w, &
29 chem, ne_area, &
30 current_month, &
31 T2,swdown, &
32 nmegan, EFmegan, msebio_isop, &
33 mlai, &
34 pftp_bt, pftp_nt, pftp_sb, pftp_hb, &
35 mtsa, &
36 mswdown, &
37 mebio_isop, mebio_apin, mebio_bpin, mebio_bcar, &
38 mebio_acet, mebio_mbo, mebio_no, &
39 ebio_iso,ebio_oli,ebio_api,ebio_lim, &
40 ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
41 ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_no, &
42 e_bio, &
43 ids,ide, jds,jde, kds,kde, &
44 ims,ime, jms,jme, kms,kme, &
45 its,ite, jts,jte, kts,kte )
46
47 USE module_configure
48 USE module_state_description
49 USE module_data_megan2
50 USE module_data_mgn2mech
51 ! USE module_bioemi_beis311, ONLY : getpar, calc_zenithb
52
53 IMPLICIT NONE
54
55 ! Subroutine arguments
56
57 ! ...simulation parameters
58 TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
59
60 ! ...domain id, current time step counter, xyz indices ..
61 INTEGER, INTENT(IN ) :: id,ktau, &
62 ids,ide, jds,jde, kds,kde, &
63 ims,ime, jms,jme, kms,kme, &
64 its,ite, jts,jte, kts,kte
65
66 ! ...current julian day
67 INTEGER, INTENT (IN) :: julday
68 !...GTM hour of start of simulation, time step in seconds
69 REAL, INTENT(IN) :: gmt,dtstep
70
71 ! ...3rd dimension size of array e_bio
72 INTEGER, INTENT (IN) :: ne_area
73
74 !...pressure (Pa)
75 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), &
76 INTENT(IN) :: p_phy
77
78 !...latitude and longitude (degrees)
79 REAL, DIMENSION( ims:ime , jms:jme ), &
80 INTENT(IN) :: xlat, xlong
81
82 !... air density (kg air/m3)
83 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
84 INTENT(IN) :: rho_phy
85
86 !...full layer height (m)
87 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
88 INTENT(IN) :: dz8w
89
90 !...2-meter temperature (K)
91 REAL, DIMENSION( ims:ime , jms:jme ), &
92 INTENT(IN) :: T2
93
94 !...downward shortwave surface flux (W/m2)
95 REAL, DIMENSION( ims:ime , jms:jme ), &
96 INTENT(IN) :: swdown
97
98 !...Number of MEGAN v2.04 species as specified by the namelist
99 !...variable nmegan; nmegan should equal n_spca_spc (this will
100 !...be checked later.) Currently nmegan=n_spca_spc=138.
101 INTEGER, INTENT(IN) :: nmegan
102
103 !...Emissions factors for nmegan=n_spca_spc=138 MEGAN v2.04 species
104 REAL, DIMENSION (ims:ime, jms:jme , nmegan) , &
105 INTENT(INOUT) :: EFmegan
106
107 !...Emission factor for isoprene (read in from file
108 !...(wrfbiochemi_d<domain>)
109 !...(moles compound/km^2/hr)
110 REAL, DIMENSION( ims:ime , jms:jme ), &
111 INTENT(IN ) :: msebio_isop
112
113 !...Plant functional group percentage (read in from file
114 !...(wrfbiochemi_d<domain>)
115 REAL, DIMENSION ( ims:ime , jms:jme ), &
116 INTENT(IN) :: &
117 pftp_bt, pftp_nt, pftp_sb, pftp_hb
118
119 !..."Climatological" Leaf area index (read in from file
120 !...(wrfbiochemi_d<domain>)
121 REAL, DIMENSION( ims:ime , jms:jme , 12 ), &
122 INTENT(IN) :: mlai
123
124 !..."Climatological" surface air temperature (K) (read in from file
125 !...(wrfbiochemi_d<domain>)
126 REAL, DIMENSION( ims:ime , jms:jme , 12 ), &
127 INTENT(IN) :: mtsa
128
129 !..."Climatological" downward radiation (W/m2) (read in from file
130 !...(wrfbiochemi_d<domain>)
131 REAL, DIMENSION ( ims:ime , jms:jme , 12 ), &
132 INTENT(IN) :: mswdown
133
134 !...Actual emissions for a few selected species as diagnostics, using
135 !...MEGAN v2.0 classes of species classification
136 !...(mol km-2 hr-1)
137 REAL, DIMENSION( ims:ime , jms:jme ), &
138 INTENT(INOUT) :: &
139 mebio_isop, mebio_apin, mebio_bpin, mebio_bcar, &
140 mebio_acet, mebio_mbo, mebio_no
141
142 !...Actual biogenic emissions, converted to mechanisms species.
143 !...(ppm m/min)
144 REAL, DIMENSION( ims:ime, jms:jme, ne_area ), &
145 INTENT(INOUT ) :: e_bio
146
147 !...Actual biogenic emissions, converted to mechanisms species.
148 !...These variables were originally for BEIS3.11 biogenic emissions
149 !...modules.
150 !...(moles compound/km^2/hr)
151 REAL, DIMENSION( ims:ime , jms:jme ), &
152 INTENT(INOUT ) :: &
153 ebio_iso,ebio_oli,ebio_api,ebio_lim, &
154 ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald, &
155 ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_no
156
157 !...Array of chemical concentrations
158 !... in - concentrations before biogenic emissions
159 !... out - concentrations after biogeniec emissions
160 !... gas-phace concentrations are in ppm
161 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
162 INTENT(INOUT) :: chem
163
164 !...Current month
165 INTEGER, INTENT(IN) :: current_month
166
167 ! Local parameters
168
169 !...Below which set emissions rate to zero (mol km-2 hr-1)
170 REAL, PARAMETER :: min_emis = 0.001
171
172
173 !...number of days in each month
174 INTEGER, PARAMETER :: DaysInMonth(12) = (/ &
175 31,28,31,30,31,30,31,31,30,31,30,31 /)
176 !...conversion between radians and degrees
177 REAL, PARAMETER :: PI = 3.14159
178 REAL, PARAMETER :: D2RAD = PI/180.0
179
180
181 ! Local Scalars
182
183 CHARACTER(len=256) :: mesg
184 INTEGER :: i,j,k,i_class, i_spc, icount, p_in_chem
185 INTEGER :: previous_month
186
187 !...minutes since start of run to the middle of the
188 !...current times step (seconds included as decimals)
189 REAL :: xtime
190
191 !...the GMT hour of the middle of the current time step
192 !...(can be greater than 24)
193 INTEGER :: ixhour
194 REAL :: xhour
195
196 !...minutes past the previous hour mark, at the
197 !...middle of the current time step
198 REAL :: xmin
199
200 !...the GMT hour of the middle of the current time step
201 !...(between 0 and 24)
202 REAL :: gmtp
203
204 !...GMT hour plus minutes (in fractaionl hour) of the middle
205 !...of the current time step
206 REAL :: tmidh
207
208 !...Current and previous month leaf area index
209 REAL :: LAIc, LAIp
210
211 !...temperature((K) and pressure (mb)
212 REAL :: tsa, tsa24, pres
213
214 !...latitude and longitude (degrees)
215 REAL :: lat, lon
216
217 !...downward solar radiation, current and some 24-hour mean (W/m2)
218 REAL :: swd, swd24
219 !...photosynthetic photon flux density (i.e. PPDF or PAR)
220 !...(micromole m-2 s-1)
221 REAL :: par, par24, pardb, pardif
222
223 !...solar zenith angle (radians), cosine of zenith angle
224 REAL :: zen , coszen
225
226 !...days in the previous month
227 REAL :: tstlen
228
229 !...emissions factor (microgram m-2 hr-1)
230 REAL :: epsilon
231
232 !...MEGAN v2.04 emissions adjustment factors for leaf area, temperature,
233 !...light, leaf age, and soil moisture
234 !...(dimensionless)
235 REAL :: gam_LHT, gam_TMP, gam_PHO,gam_AGE, gam_SMT
236
237 !...normalized ratio accounting for production and loss within
238 !...plant canopies (dimensionless)
239 REAL :: rho
240
241 !...Some light-dependent factor (dimensionless)
242 REAL :: ldf
243
244 !...conversion factor to convert emissions rates in
245 !...mol km-2 hr-1 to concentrations in ppm
246 REAL :: emis2ppm
247
248 !...conversion factor from mol km-2 hr-1 to ppm m min-1
249 REAL :: convert2
250
251 !...emission rate converted to mechanism species in mol km-2 hr-1
252 REAL :: gas_emis
253
254 ! Local Arrays
255
256 !...emissions adjustment factors for n_mgn_spc=20 classes of
257 !...MEGAN v2.04 specie.
258 !...adjust_factor = [GAMMA]*[rho] (see comments later)
259 !...(dimensionless)
260 REAL, DIMENSION(n_mgn_spc) :: adjust_factor
261
262 !...plant functional type fractions
263 REAL :: pft_frac(n_pft)
264
265 !...actually emissions rates of n_spca_spc=138 MEGAN v2.04 species
266 !...(mol km-2 hr-2)
267 REAL, DIMENSION ( n_spca_spc ) :: E_megan2
268
269 ! End header ------------------------------------------------------
270
271
272 ! MEGAN v2.04 has nmegan=n_spca_spc=138 species, which are lumped
273 ! into n_mgn_spc=20 classes. The number, names and indices of
274 ! these classes and species are defined in module_data_megan2.F.
275 ! They need to follow a few rules
276
277 IF ( ktau .EQ. 1 ) THEN
278
279 ! The size of variable EFmegan(:,:,nmegan) is allocated based on
280 ! the value of namelist variable nmegan. nmegan should be equal
281 ! to n_spca_spc (though can be greater than to n_spca_spc).
282 IF ( nmegan .NE. n_spca_spc ) THEN
283 WRITE(mesg,*)'namelist variable nmegan does not match n_spca_spc'
284 CALL wrf_error_fatal(mesg)
285 END IF
286
287 ! For programming, the ordering of the species or classes of
288 ! species should not matter, except that isoprene should always
289 ! be first; therefore, imgn_isop=1 and is_isoprene=1 always.
290 IF ( (imgn_isop .NE. 1) .OR. (is_isoprene .NE. 1) ) THEN
291 WRITE(mesg,*)'imgn_isop and is_isoprene in bio_emissions_megan should be 1'
292 CALL wrf_error_fatal(mesg)
293 END IF
294
295 END IF
296
297
298 ! Initialize diagnostic output
299 ebio_iso ( its:ite , jts:jte ) = 0.0
300 ebio_oli ( its:ite , jts:jte ) = 0.0
301 ebio_api ( its:ite , jts:jte ) = 0.0
302 ebio_lim ( its:ite , jts:jte ) = 0.0
303 ebio_hc3 ( its:ite , jts:jte ) = 0.0
304 ebio_ete ( its:ite , jts:jte ) = 0.0
305 ebio_olt ( its:ite , jts:jte ) = 0.0
306 ebio_ket ( its:ite , jts:jte ) = 0.0
307 ebio_ald ( its:ite , jts:jte ) = 0.0
308 ebio_hcho ( its:ite , jts:jte ) = 0.0
309 ebio_eth ( its:ite , jts:jte ) = 0.0
310 ebio_ora2 ( its:ite , jts:jte ) = 0.0
311 ebio_co ( its:ite , jts:jte ) = 0.0
312 ebio_no ( its:ite , jts:jte ) = 0.0
313 e_bio ( its:ite , jts:jte , 1:ne_area) = 0.0
314
315 !...the following is redundant if there is no
316 !...bug in the subroutine
317 mebio_isop ( its:ite , jts:jte ) = 0.0
318 mebio_apin ( its:ite , jts:jte ) = 0.0
319 mebio_bpin ( its:ite , jts:jte ) = 0.0
320 mebio_bcar ( its:ite , jts:jte ) = 0.0
321 mebio_acet ( its:ite , jts:jte ) = 0.0
322 mebio_mbo ( its:ite , jts:jte ) = 0.0
323 mebio_no ( its:ite , jts:jte ) = 0.0
324
325
326 ! Extract climatological values for relevant months.
327 !
328 ! In MEGAN v2.04, emissions rates dependent on ambient conditions
329 ! of the past 24 hours to the past month or so. The implementation
330 ! of MEGAN v2.04 here uses monthly-mean values of the previous
331 ! month for any past history required by the model. The monthly-
332 ! -mean values should be provided as input in the
333 ! wrfbiochemi_d<domain> file. Fully implementation (not done here)
334 ! require online calculations of 24-hour and 240-hour mean of
335 ! surface air temperature and downward PAR
336 !
337 ! MEGAN v2.04 also requires time-dependent leaf area index to
338 ! estimate leaf age. Here, leaf area indices of the current
339 ! and the previous months are used. The data should be
340 ! provided in wrfbiochemi_d<domain> file.
341
342 IF (current_month > 1) THEN
343 previous_month = current_month -1
344 ELSE
345 previous_month = 12
346 END IF
347
348
349 ! Following module_phot_fastj.F, determine current
350 ! time of day in GMT at the middle of the current
351 ! time step, tmidh.
352 ! ktau - time step counter
353 ! dstep - time step in seconds
354 ! gmt - starting hour (in GMT) of the simulation
355
356 !...minutes since start of run to the middle of the
357 !...current times step (seconds included as decimals)
358 xtime=(ktau-1)*dtstep/60. + dtstep/120.
359 !...the GMT hour of the middle of the current time step
360 !...(can be greater than 24)
361 ixhour = ifix(gmt + 0.01) + ifix(xtime/60.)
362 xhour=float(ixhour)
363 !...minutes past the previous hour mark, at the
364 !...middle of the current time step
365 xmin = 60.*gmt + (xtime-xhour*60)
366 !...the GMT hour of the middle of the current time step
367 !...(between 0 and 24)
368 gmtp=MOD(xhour,24.)
369 !...GMT hour plus minutes (in fractaionl hour) of the middle
370 !...of the current time step
371 tmidh= gmtp + xmin/60.
372
373 WRITE(mesg,*) 'calculate MEGAN emissions at ktau, gmtp, tmidh = ',ktau, gmtp, tmidh
374 CALL wrf_message(mesg)
375
376
377 ! Get the mechanism converstion table
378 ! ( Even though the mechanism converstion table is time-independent,
379 ! do this for all time steps to be sure there will be no issue with
380 ! restart runs. This should be edited eventually to reduce
381 ! redundant calculations.)
382 ! SHC (11/08/2007)
383 GAS_MECH_SELECT1: SELECT CASE (config_flags%chem_opt)
384
385 CASE (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_KPP)
386 ! get p_of_radm2cbmz(:), p_of_radm2(:), and radm2_per_megan(:)
387 CALL get_megan2radm2_table
388
389 CASE (RACM, RACMSORG,RACM_KPP,RACMSORG_KPP, RACM_MIM_KPP )
390
391 ! get p_of_megan2racm(:), p_of_racm(:), and racm_per_megan(:)
392 CALL get_megan2racm_table
393
394 CASE (CBMZ, CBMZ_BB, CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ)
395
396 ! get p_of_megan2cbmz(:), p_of_cbmz(:), and cbmz_per_megan(:)
397 CALL get_megan2cbmz_table
398
399 CASE DEFAULT
400
401 CALL wrf_error_fatal('Species conversion table for MEGAN v2.04 not available. ')
402
403 END SELECT GAS_MECH_SELECT1
404
405 ! Calcuate biogenic emissions grid by grid
406
407 j_loop: DO j = jts, jte
408 i_loop: DO i = its, ite
409
410
411 ! Put variables of ambient conditions into scalar variables
412
413 tsa = T2(i,j) ! air temperature at 2-meter (K)
414 pres = 0.01*p_phy(i,kts,j) ! surface pressure (mb)
415 lat = xlat(i,j) ! latitude (degree)
416 lon = xlong(i,j) ! longitude (degress)
417 swd = swdown(i,j) ! downward solar radiation (W/m2)
418 LAIc = mlai(i,j,current_month) ! current month leaf area index
419 LAIp = mlai(i,j,previous_month) ! previous month leaf area index
420
421 !...Emissions are dependent on the ambient conditions in the last
422 !...24 to 240 hours; here, use input data for monthly mean of the
423 !...the previous month
424 tsa24 = mtsa (i,j,previous_month) ! [=]K
425 swd24 = mswdown (i,j,previous_month) ! [=] W/m2
426
427 !...Perform checks on max and min bounds for temperature
428 IF (tsa .LT. 200.0) THEN
429 WRITE (mesg,'("temperature too low at i=",i3," ,j=",i3 )')i,j
430 CALL wrf_message(mesg)
431 END IF
432 IF (tsa .GT. 315.0 ) THEN
433 WRITE (mesg,'("temperature too high at i=",i3," ,j=",i3," ;resetting to 315K" )')i,j
434 CALL wrf_message(mesg)
435 tsa = 315.0
436 END IF
437
438 ! !...Calculate zenith angle (in radians)
439 ! !...(NOTE: nonstandard longitude input here: >0 for W, <0 for E!!)
440 ! !...(subroutine calc_zenithb is in module_bioemis_beis311.F)
441 ! CALL calc_zenithb(lat,-lon,julday,tmidh,zen)
442 ! coszen = COS(zen)
443
444 !....Convert downward solar radiation to photosynthetically
445 !....active radiation
446
447 ! !......(subroutine getpar is in module_bioemis_beis311.F)
448 ! CALL getpar( swd, pres, zen, pardb, pardif )
449 ! par = pardb + pardif ! micro-mole/m2/s
450
451 !......assume 4.766 (umol m-2 s-1) per (W m-2)
452 !......assume 1/2 of swd is in 400-700 nm band
453 par = 4.766 * 0.5 * swd
454
455 !......Check max/min bounds of PAR
456 IF ( par .LT. 0.00 .OR. par .GT. 2600.0 ) THEN
457 WRITE (mesg,'("par out of range at i=",i3," ,j=",i3," par =",f8.3 )')i,j,par
458 CALL wrf_message(mesg)
459 END IF
460
461 !......For the 24-avg PAR,
462 !......assume 4.766 (umol m-2 s-1) per (W m-2)
463 !......assume 1/2 of swd is in 400-700 band
464 par24 = swd24 * 4.766 * 0.5
465
466 ! ------------------------------------------------------------
467 !
468 ! MEGAN v2.04 Box Model
469 !
470 ! Reference:
471 !
472 ! Estimates of global terrestial isoprene emissions using MEGAN
473 ! (Model of Emissions of Gases and Aerosols from Nature )
474 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer,
475 ! P.I. Palmer, and C. Geron
476 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006
477 !
478 ! MEGAN v2.0 Documentation
479 !
480 !
481 ! The following code is based on Tan's megan.F dated 11/21/2006
482 !
483 ! Scientific algorithm
484 !
485 ! Emission = [epsilon][gamma][rho]
486 !
487 ! where [epsilon] = emission factor (usually um m-2 hr or mole km-2 hr-1)
488 ! [gamma] = emission activity factor (dimensionless)
489 ! [rho] = production and loss within plant canopies
490 ! (dimensionless)
491 !
492 ! [gamma] = [gamma_CE][gamma_age][gamma_SM]
493 !
494 ! where [gamma_CE] = canopy correction factor
495 ! [gamma_age] = leaf age correction factor
496 ! [gamma_SM] = soil moisture correction factor
497 !
498 ! [gamma_CE] = [gamma_LAI][gamma_T]((1-LDF) + LDF*[gamma_P] )
499 !
500 ! where [gamma_LAI] = leaf area index factor
501 ! [gamma_P] = PPFD emission activity factor
502 ! [gamma_T] = temperature response factor
503 ! LDF =
504 !
505 ! Emission = [epsilon][gamma_LAI][gamma_T][gamma_age]
506 ! x ((1-LDF) + LDF*[gamma_P] )[rho]
507 !
508 ! or
509 !
510 ! Emission = adjust_factor [epsilon]
511 !
512 ! where
513 !
514 ! adjust_fact = [gamma_LAI][gamma_T][gamma_age]((1-LDF) + LDF*[gamma_P] )[rho]
515
516
517 !
518 ! Calculate the dimensionless emission adjustment factor.
519 ! MEGAN v2.04 has n_spca_spc = 138 species. These species are
520 ! lumped into n_mgn_spc=20 classes. The emission adjustment
521 ! factors are different for the 20 classes of species.
522 !
523 ! NOTE: This version of the code contains the corrected equation for
524 ! gamma P (based on a revised version of equation 11b from Guenther et al., 2006)
525 ! CW (08/16/2007)
526 !
527 ! NOTE: This version of the code contains the updated emission factors (static)
528 ! and beta values based on Alex's V2.1 notes sent out on August 13, 2007
529 ! CW (08/16/2007)
530 !
531 ! NOTE: This version of the code applies the same gamma T equation to the
532 ! emissions of all compounds other than isoprene. This occurs regardless
533 ! of whether the emissions are light dependent or not. This is NOT the same
534 ! as what Alex has in his code. In his code, the light-dependent emissions
535 ! are also given the isoprene gamma T. Because all emissions (other than isoprene
536 ! are assigned the same gamma T, this could lead to overestimates of emissions
537 ! at high temperatures (>40C). Light-dependent emisisons (e.g., SQTs) should
538 ! fall off at high temperatures. (CW, 08/16/2007)
539
540 !...Calculate adjustment factor components that are species-independent
541
542 !......Get the leaf area index factor gam_LHT
543 CALL GAMMA_LAI( LAIc, gam_LHT)
544
545 !......Get the light emission activity factor gam_PHO
546 CALL GAMMA_P( julday, tmidh, lat, lon, par, par24, gam_PHO )
547
548 !......Get the soil moisture factor gam_SMT
549 !......for now, set = 1.0
550 gam_SMT = 1.0
551
552 !...Calculate the overall emissions adjustment factors, for
553 !...each of the n_mgn_spc=20 classes of compounds
554
555 DO i_class = 1, n_mgn_spc
556
557 ! Get the temperature response factor gam_TMP
558 ! One algorithm for isoprene, and one for non-isoprene
559
560 IF ( i_class == imgn_isop ) THEN
561 CALL GAMMA_TISOP( tsa, tsa24, gam_TMP )
562 ELSE
563 CALL GAMMA_TNISP( i_class , tsa, gam_TMP )
564 END IF
565
566 ! Get the leaf age correction factor gam_AGE
567
568 !...Time step (days) between LAIc and LAIp:
569 !...Since monthly mean LAI is used,
570 !...use # of days in the previous month
571 tstlen = REAL(DaysInMonth(previous_month),KIND(1.0))
572
573 CALL GAMMA_A( i_class , LAIp, LAIc, TSTLEN, tsa24, gam_AGE )
574
575 ! rho - normalized ratio accounting for production and
576 ! oss within plant canopies; rho_fct is defined in
577 ! module_data_megan2.F; currently rho_fct = 1.0 for all
578 ! species (dimensionless)
579 rho = rho_fct(i_class)
580
581 ! Fraction of emission to apply light-dependence factor
582 ! ldf_fct is defined in module_data_megan2.F
583 ! (dimensionless)
584 ldf = ldf_fct(i_class)
585
586 ! The overall emissions adjustment factor
587 ! (dimensionless)
588
589 adjust_factor(i_class) = gam_TMP * gam_AGE * gam_LHT * gam_SMT * rho * &
590 ( (1.0-LDF) + gam_PHO*LDF )
591
592 END DO !i_class = 1, n_mgn_spc (loop over classes of MEGAN species )
593
594
595 ! For isoprene, the emission factor is already read in from
596 ! wrfbiochemi_d<domain> file; therefore, actual emissions rate
597 ! can be calculated here already.
598 ! (mol km-2 hr-1)
599 E_megan2(is_isoprene) = adjust_factor(imgn_isop)*msebio_isop(i,j)
600 IF ( E_megan2(is_isoprene) .LT. min_emis ) E_megan2(is_isoprene)=0.
601
602
603 ! Calculate emissions for all n_spca_spc=nmegan=138 MEGAN v2.04
604 ! species, except for isoprene. For non-isoprene emissions,
605 ! the emission factor [epsilon] has to be calculated
606 ! for the first time step.
607
608 !...Loop over species, because emission factor [epsilon] is
609 !...different for each species
610 !...( i_spc = 1 is skipped in the do loop below to skip
611 !...isoprene; this works because is_isoprene = 1 )
612 DO i_spc = 2, n_spca_spc
613
614 ! The lumped class in which the current species is a member
615 i_class = mg20_map (i_spc)
616
617 ! Calculate emission factor (microgram m-2 hr-1) for
618 ! species i_spc
619 ! ( Even though EFmegan is time invariant, for now calculate
620 ! EFmegan for every time step to be sure there will be
621 ! no issue with restart runs.
622 ! SHC (11/08/2007)
623 ! IF ( ktau .EQ. 1 ) THEN
624
625 ! Grab plant functional type fractions for current grid
626 ! cell (pftp_* is the plant functional type % and was
627 ! read in from wrfbiochemi_d<domain> file.)
628 pft_frac(k_bt) = 0.01*pftp_bt(i,j)
629 pft_frac(k_nt) = 0.01*pftp_nt(i,j)
630 pft_frac(k_sb) = 0.01*pftp_sb(i,j)
631 pft_frac(k_hb) = 0.01*pftp_hb(i,j)
632
633 ! Sum up emissions factor over plant functional types
634 epsilon = 0.0
635 DO k = 1, n_pft !loop over plant functional types
636 epsilon = epsilon + &
637 pft_frac(k)*EF(i_class,k)*EF_frac(i_spc,k)
638 END DO
639
640 ! Store emission factor to variable EFmegan (which is
641 ! declared in Registry/registry.chem)
642 ! (migrogram m-2 hr-1)
643 EFmegan(i,j,i_spc) = epsilon
644
645 ! END IF ! ( ktau .EQ. 1 )
646
647 ! Calculate actual emission rate for species i_spc;
648 ! also, convert units from (microgram m-2 hr-1) to
649 ! (mol km-2 hr-1)
650 E_megan2(i_spc) = EFmegan(i,j,i_spc)* &
651 adjust_factor(i_class)/spca_mwt(i_spc)
652 IF ( E_megan2(i_spc) .LT. min_emis ) E_megan2(i_spc)=0.
653
654 END DO !i_spc = 2, n_spca_spc, loop over all non-isoprene species
655
656
657 ! Output emissions for some species as diagnostics
658 ! (mol km-2 hr-1)
659 mebio_isop (i,j) = E_megan2 ( is_isoprene )
660 mebio_apin (i,j) = E_megan2 ( is_pinene_a )
661 mebio_bpin (i,j) = E_megan2 ( is_pinene_b )
662 mebio_bcar (i,j) = E_megan2 ( is_caryophyllene_b )
663 mebio_acet (i,j) = E_megan2 ( is_acetone )
664 mebio_mbo (i,j) = E_megan2 ( is_MBO_2m3e2ol )
665 mebio_no (i,j) = E_megan2 ( is_nitric_OXD )
666
667
668 ! Speciate the n_spca_spc=nmegan=138 species into
669 ! the gas-phase mechanism species
670
671
672 !...conversion factor to convert emissions rates in
673 !...mol km-2 hr-1 to concentrations in ppm
674 ! 0.02897 kg/mol is molecular of air
675 ! rho_phy is air density in kg air/m3
676 ! dz8w is layer height in meters
677 ! dtstep is time step in seconds
678 emis2ppm = 0.02897*dtstep/(rho_phy(i,kts,j)*dz8w(i,kts,j)*3600.)
679
680 !...conversion factor from mol km-2 hr-1 to ppm m min-1
681 !...(e_bio is in units of ppm m min-1)
682 convert2 = 0.02897/(rho_phy(i,kts,j)*60.)
683
684
685 !...
686 GAS_MECH_SELECT: SELECT CASE (config_flags%chem_opt)
687
688 CASE (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_KPP)
689
690 DO icount = 1, n_megan2radm2
691
692 IF ( p_of_radm2(icount) .NE. non_react ) THEN
693
694 ! Get index to chem array for the corresponding RADM2
695 ! species.
696 p_in_chem = p_of_radm2(icount)
697
698 ! Check if the species is actually in the mechanism
699 IF ( p_in_chem > param_first_scalar ) THEN
700
701 ! Emission rate for mechanism species in mol km-2 hr-1
702 gas_emis = radm2_per_megan(icount) * E_megan2(p_of_megan2radm2(icount))
703
704 ! Increase gas-phase concentrations (in ppmv) due to
705 ! biogenic emissions
706 chem(i,kts,j,p_in_chem) = chem(i,kts,j,p_in_chem) + gas_emis*emis2ppm
707
708 ! Add emissions to diagnostic output variables.
709 ! ebio_xxx (mol km-2 hr-1) were originally used by the
710 ! BEIS3.11 biogenic emissions module.
711 ! I have also borrowed variable e_bio (ppm m min-1).
712 IF ( p_in_chem .EQ. p_iso ) THEN
713 ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
714 e_bio(i,j,p_iso-1) = e_bio(i,j,p_iso-1) + gas_emis*convert2
715 ELSE IF ( p_in_chem .EQ. p_oli) THEN
716 ebio_oli(i,j) = ebio_oli(i,j) + gas_emis
717 e_bio(i,j,p_oli-1) = e_bio(i,j,p_oli-1) + gas_emis*convert2
718 ELSE IF ( p_in_chem .EQ. p_hc3) THEN
719 ebio_hc3(i,j) = ebio_hc3(i,j) + gas_emis
720 e_bio(i,j,p_hc3-1) = e_bio(i,j,p_hc3-1) + gas_emis*convert2
721 ELSE IF ( p_in_chem .EQ. p_olt) THEN
722 ebio_olt(i,j) = ebio_olt(i,j) + gas_emis
723 e_bio(i,j,p_olt-1) = e_bio(i,j,p_olt-1) + gas_emis*convert2
724 ELSE IF ( p_in_chem .EQ. p_ket) THEN
725 ebio_ket(i,j) = ebio_ket(i,j) + gas_emis
726 e_bio(i,j,p_ket-1) = e_bio(i,j,p_ket-1) + gas_emis*convert2
727 ELSE IF ( p_in_chem .EQ. p_ald) THEN
728 ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
729 e_bio(i,j,p_ald-1) = e_bio(i,j,p_ald-1) + gas_emis*convert2
730 ELSE IF ( p_in_chem .EQ. p_hcho) THEN
731 ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
732 e_bio(i,j,p_hcho-1) = e_bio(i,j,p_hcho-1) + gas_emis*convert2
733 ELSE IF ( p_in_chem .EQ. p_eth) THEN
734 ebio_eth(i,j) = ebio_eth(i,j) + gas_emis
735 e_bio(i,j,p_eth-1) = e_bio(i,j,p_eth-1) + gas_emis*convert2
736 ELSE IF ( p_in_chem .EQ. p_ora2) THEN
737 ebio_ora2(i,j) = ebio_ora2(i,j) + gas_emis
738 e_bio(i,j,p_ora2-1) = e_bio(i,j,p_ora2-1) + gas_emis*convert2
739 ELSE IF ( p_in_chem .EQ. p_co) THEN
740 ebio_co(i,j) = ebio_co(i,j) + gas_emis
741 e_bio(i,j,p_co-1) = e_bio(i,j,p_co-1) + gas_emis*convert2
742 ELSE IF ( p_in_chem .EQ. p_no) THEN
743 ebio_no(i,j) = ebio_no(i,j) + gas_emis
744 e_bio(i,j,p_no-1) = e_bio(i,j,p_no-1) + gas_emis*convert2
745 ELSE IF ( p_in_chem .EQ. p_ol2) THEN
746 e_bio(i,j,p_ol2-1) = e_bio(i,j,p_ol2-1) + gas_emis*convert2
747 ELSE IF ( p_in_chem .EQ. p_hc5) THEN
748 e_bio(i,j,p_hc5-1) = e_bio(i,j,p_hc5-1) + gas_emis*convert2
749 ELSE IF ( p_in_chem .EQ. p_hc8) THEN
750 e_bio(i,j,p_hc8-1) = e_bio(i,j,p_hc8-1) + gas_emis*convert2
751 ELSE IF ( p_in_chem .EQ. p_ora1) THEN
752 e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
753 END IF
754
755 END IF !( p_in_chem > param_first_scalar )
756
757 END IF !( p_of_ramd2(icount) .NE. non_react )
758
759 END DO
760
761 CASE (RACM, RACMSORG,RACM_KPP,RACMSORG_KPP, RACM_MIM_KPP )
762
763 DO icount = 1, n_megan2racm
764
765 IF ( p_of_racm(icount) .NE. non_react ) THEN
766
767 ! Get index to chem array for the corresponding RACM
768 ! species.
769 p_in_chem = p_of_racm(icount)
770
771 ! Check if the species is actually in the mechanism
772 IF( p_in_chem > param_first_scalar ) THEN
773
774 ! Emission rate of mechanism species in mol km-2 hr-1
775 gas_emis = racm_per_megan(icount) * E_megan2(p_of_megan2racm(icount))
776
777 ! Increase gas-phase concentrations (in ppmv) due to
778 ! biogenic emissions
779 chem(i,kts,j,p_in_chem) = chem(i,kts,j,p_in_chem) + gas_emis * emis2ppm
780
781 ! Add emissions to diagnostic output variables.
782 ! ebio_xxx (mol km-2 hr-1) were originally used by the
783 ! BEIS3.11 biogenic emissions module.
784 ! I have also borrowed variable e_bio (ppm m min-1).
785 IF ( p_in_chem .EQ. p_iso ) THEN
786 ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
787 e_bio(i,j,p_iso-1) = e_bio(i,j,p_iso-1) + gas_emis*convert2
788 ELSE IF ( p_in_chem .EQ. p_oli) THEN
789 ebio_oli(i,j) = ebio_oli(i,j) + gas_emis
790 e_bio(i,j,p_oli-1) = e_bio(i,j,p_oli-1) + gas_emis*convert2
791 ELSE IF ( p_in_chem .EQ. p_api) THEN
792 ebio_api(i,j) = ebio_api(i,j) + gas_emis
793 e_bio(i,j,p_api-1) = e_bio(i,j,p_api-1) + gas_emis*convert2
794 ELSE IF ( p_in_chem .EQ. p_lim) THEN
795 ebio_lim(i,j) = ebio_lim(i,j) + gas_emis
796 e_bio(i,j,p_lim-1) = e_bio(i,j,p_lim-1) + gas_emis*convert2
797 ELSE IF ( p_in_chem .EQ. p_hc3) THEN
798 ebio_hc3(i,j) = ebio_hc3(i,j) + gas_emis
799 e_bio(i,j,p_hc3-1) = e_bio(i,j,p_hc3-1) + gas_emis*convert2
800 ELSE IF ( p_in_chem .EQ. p_ete) THEN
801 ebio_ete(i,j) = ebio_ete(i,j) + gas_emis
802 e_bio(i,j,p_ete-1) = e_bio(i,j,p_ete-1) + gas_emis*convert2
803 ELSE IF ( p_in_chem .EQ. p_olt) THEN
804 ebio_olt(i,j) = ebio_olt(i,j) + gas_emis
805 e_bio(i,j,p_olt-1) = e_bio(i,j,p_olt-1) + gas_emis*convert2
806 ELSE IF ( p_in_chem .EQ. p_ket) THEN
807 ebio_ket(i,j) = ebio_ket(i,j) + gas_emis
808 e_bio(i,j,p_ket-1) = e_bio(i,j,p_ket-1) + gas_emis*convert2
809 ELSE IF ( p_in_chem .EQ. p_ald) THEN
810 ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
811 e_bio(i,j,p_ald-1) = e_bio(i,j,p_ald-1) + gas_emis*convert2
812 ELSE IF ( p_in_chem .EQ. p_hcho) THEN
813 ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
814 e_bio(i,j,p_hcho-1) = e_bio(i,j,p_hcho-1) + gas_emis*convert2
815 ELSE IF ( p_in_chem .EQ. p_eth) THEN
816 ebio_eth(i,j) = ebio_eth(i,j) + gas_emis
817 e_bio(i,j,p_eth-1) = e_bio(i,j,p_eth-1) + gas_emis*convert2
818 ELSE IF ( p_in_chem .EQ. p_ora2) THEN
819 ebio_ora2(i,j) = ebio_ora2(i,j) + gas_emis
820 e_bio(i,j,p_ora2-1) = e_bio(i,j,p_ora2-1) + gas_emis*convert2
821 ELSE IF ( p_in_chem .EQ. p_co) THEN
822 ebio_co(i,j) = ebio_co(i,j) + gas_emis
823 e_bio(i,j,p_co-1) = e_bio(i,j,p_co-1) + gas_emis*convert2
824 ELSE IF ( p_in_chem .EQ. p_no) THEN
825 ebio_no(i,j) = ebio_no(i,j) + gas_emis
826 e_bio(i,j,p_no-1) = e_bio(i,j,p_no-1) + gas_emis*convert2
827 ELSE IF ( p_in_chem .EQ. p_hc5) THEN
828 e_bio(i,j,p_hc5-1) = e_bio(i,j,p_hc5-1) + gas_emis*convert2
829 ELSE IF ( p_in_chem .EQ. p_hc8) THEN
830 e_bio(i,j,p_hc8-1) = e_bio(i,j,p_hc8-1) + gas_emis*convert2
831 ELSE IF ( p_in_chem .EQ. p_ora1) THEN
832 e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
833 END IF
834
835 END IF !( p_in_chem > param_first_scalar )
836
837
838 END IF !( p_of_racm(icount) .NE. non_react )
839
840 END DO
841
842 CASE (CBMZ, CBMZ_BB, CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ)
843
844 DO icount = 1, n_megan2cbmz
845
846 IF ( p_of_cbmz (icount) .NE. non_react ) THEN
847
848 ! Get index to chem array for the corresponding CBMZ
849 ! species.
850 p_in_chem = p_of_cbmz(icount)
851
852 ! Check if the species is actually in the mechanism
853 ! (e.g., NH3 is in the mechanism only if aerosols
854 ! are simulated)
855 IF( p_in_chem > param_first_scalar ) THEN
856
857 ! Emission rate of mechanism species in mol km-2 hr-1
858 gas_emis = cbmz_per_megan(icount) * E_megan2(p_of_megan2cbmz(icount))
859
860 ! Increase gas-phase concentrations (in ppmv) due to
861 ! biogenic emissions
862 chem(i,kts,j,p_in_chem) = chem(i,kts,j,p_in_chem) + gas_emis * emis2ppm
863
864 ! Add emissions to diagnostic output variables.
865 ! ebio_xxx (mol km-2 hr-1) were originally used by the
866 ! BEIS3.11 biogenic emissions module.
867 ! I have also borrowed variable e_bio (ppm m min-1).
868 IF ( p_in_chem .EQ. p_iso ) THEN
869 ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
870 e_bio(i,j,p_iso-1) = e_bio(i,j,p_iso-1) + gas_emis*convert2
871 ELSE IF ( p_in_chem .EQ. p_oli) THEN
872 ebio_oli(i,j) = ebio_oli(i,j) + gas_emis
873 e_bio(i,j,p_oli-1) = e_bio(i,j,p_oli-1) + gas_emis*convert2
874 ELSE IF ( p_in_chem .EQ. p_olt) THEN
875 ebio_olt(i,j) = ebio_olt(i,j) + gas_emis
876 e_bio(i,j,p_olt-1) = e_bio(i,j,p_olt-1) + gas_emis*convert2
877 ELSE IF ( p_in_chem .EQ. p_ket) THEN
878 ebio_ket(i,j) = ebio_ket(i,j) + gas_emis
879 e_bio(i,j,p_ket-1) = e_bio(i,j,p_ket-1) + gas_emis*convert2
880 ELSE IF ( p_in_chem .EQ. p_ald) THEN
881 ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
882 e_bio(i,j,p_ald-1) = e_bio(i,j,p_ald-1) + gas_emis*convert2
883 ELSE IF ( p_in_chem .EQ. p_hcho) THEN
884 ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
885 e_bio(i,j,p_hcho-1) = e_bio(i,j,p_hcho-1) + gas_emis*convert2
886 ELSE IF ( p_in_chem .EQ. p_eth) THEN
887 ebio_eth(i,j) = ebio_eth(i,j) + gas_emis
888 e_bio(i,j,p_eth-1) = e_bio(i,j,p_eth-1) + gas_emis*convert2
889 ELSE IF ( p_in_chem .EQ. p_ora2) THEN
890 ebio_ora2(i,j) = ebio_ora2(i,j) + gas_emis
891 e_bio(i,j,p_ora2-1) = e_bio(i,j,p_ora2-1) + gas_emis*convert2
892 ELSE IF ( p_in_chem .EQ. p_co) THEN
893 ebio_co(i,j) = ebio_co(i,j) + gas_emis
894 e_bio(i,j,p_co-1) = e_bio(i,j,p_co-1) + gas_emis*convert2
895 ELSE IF ( p_in_chem .EQ. p_no) THEN
896 ebio_no(i,j) = ebio_no(i,j) + gas_emis
897 e_bio(i,j,p_no-1) = e_bio(i,j,p_no-1) + gas_emis*convert2
898 ELSE IF ( p_in_chem .EQ. p_ol2) THEN
899 e_bio(i,j,p_ol2-1) = e_bio(i,j,p_ol2-1) + gas_emis*convert2
900 ELSE IF ( p_in_chem .EQ. p_ora1) THEN
901 e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
902 END IF
903
904
905 END IF !( p_in_chem > param_first_scalar )
906
907
908 END IF ! ( p_of_cbmz (icount) .NE. non_react )
909
910 END DO
911
912 CASE DEFAULT
913
914 CALL wrf_error_fatal('Species conversion table for MEGAN v2.04 not available. ')
915
916 END SELECT GAS_MECH_SELECT
917
918
919
920 END DO i_loop ! i = its, ite
921 END DO j_loop ! j = jts, jte
922
923
924 CONTAINS
925
926 ! -----------------------------------------------------------------
927 ! SUBROUTINE GAMMA_TISOP returns the GAMMA_T value for isoprene
928 ! Orginally from Tan's gamma_etc.F
929 ! -----------------------------------------------------------------
930
931 SUBROUTINE GAMMA_TISOP( TEMP, D_TEMP, gam_T )
932 !
933 ! Description :
934 !
935 ! MEGAN biogenic emissions adjustment factor for temperature
936 ! for isoprene
937 !
938 ! Reference:
939 !
940 ! Estimates of global terrestial isoprene emissions using MEGAN
941 ! (Model of Emissions of Gases and Aerosols from Nature )
942 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer,
943 ! P.I. Palmer, and C. Geron
944 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006 !
945 !
946
947 IMPLICIT NONE
948
949 ! hourly surface air temperature (K)
950 ! (here use instantaneous temperature
951 REAL, INTENT(IN) :: TEMP
952 ! daily-mean surface airtemperature (K)
953 ! (here use the previous month's monthly mean)
954 REAL, INTENT(IN) :: D_TEMP
955 !temperature adjustment factor
956 REAL, INTENT(OUT) :: gam_T
957
958 ! Local parameters
959 REAL :: Eopt, Topt, X
960 REAL :: AAA, BBB
961 REAL, PARAMETER :: CT1 = 80.0
962 REAL, PARAMETER :: CT2 = 200.0
963
964 ! End header ----------------------------------------------------
965
966 ! Below Eqn (14) of Guenther et al. [2006]
967 ! (assuming T_daily = D_TEMP)
968 Eopt = 1.75 * EXP(0.08*(D_TEMP-297.0))
969
970 ! Eqn (8) of Guenther et al. [2006]
971 ! (assuming T_daily = D_TEMP)
972 Topt = 313.0 + ( 0.6*(D_TEMP-297.0) )
973
974 ! Eqn (5) of Guenther et al. [2006]
975 X = ( (1.0/Topt)-(1.0/TEMP) ) / 0.00831
976 AAA = Eopt*CT2*EXP(CT1*X)
977 BBB = ( CT2-CT1*( 1.0-EXP(CT2*X) ) )
978 gam_T = AAA/BBB
979
980 END SUBROUTINE GAMMA_TISOP
981
982 ! -----------------------------------------------------------------
983 ! SUBROUITINE GAMMA_TNISP returns the GAMMA_T value for
984 ! non-isoprene species
985 ! Originally from Tan's gamma_etc.F
986 !------------------------------------------------------------------
987
988 SUBROUTINE GAMMA_TNISP( SPCNUM, TEMP, gam_T )
989 !
990 ! Description :
991 !
992 ! MEGAN biogenic emissions adjustment factor for temperature
993 ! for non-isoprene species.
994 !
995 ! Reference:
996 !
997 ! MEGAN v2.0 Documentation
998 !
999 ! Method:
1000 !
1001 ! GAMMA_T = exp[BETA*(T-Ts)]
1002 ! where BETA = temperature dependent parameter
1003 ! Ts = standard temperature (normally 303K, 30C)
1004 !
1005
1006 IMPLICIT NONE
1007
1008 INTEGER, INTENT(IN) :: SPCNUM ! Species number
1009 REAL, INTENT(IN) :: TEMP
1010 REAL, INTENT(OUT) :: gam_T
1011 REAL, PARAMETER :: Ts = 303.0
1012
1013 ! End header ----------------------------------------------------
1014
1015 ! TDF_PRM is defined in module_data_megan2.F
1016 gam_T = EXP( TDF_PRM(SPCNUM)*(TEMP-Ts) )
1017
1018 END SUBROUTINE GAMMA_TNISP
1019
1020
1021 ! --------------------------------------------------------------------
1022 ! SUBROUTINE GAMMA_LAI
1023 ! Originally from Tan's gamma_etc.F
1024 ! --------------------------------------------------------------------
1025
1026 SUBROUTINE GAMMA_LAI(LAI, gam_L )
1027 ! Description :
1028 !
1029 ! MEGAN biogenic emissions adjustment factor for leaf area
1030 ! index
1031 !
1032 ! Reference:
1033 !
1034 ! Estimates of global terrestial isoprene emissions using MEGAN
1035 ! (Model of Emissions of Gases and Aerosols from Nature )
1036 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer,
1037 ! P.I. Palmer, and C. Geron
1038 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006 !
1039 !
1040 ! Method:
1041 ! 0.49[LAI]
1042 ! GAMMA_LAI = ---------------- (dimensionless)
1043 ! (1+0.2LAI^2)^0.5
1044 !
1045
1046 IMPLICIT NONE
1047 REAL, INTENT(IN) :: LAI
1048 REAL, INTENT(OUT) :: gam_L
1049
1050 ! End header ----------------------------------------------------
1051
1052
1053 ! Eqn (15) of Guenther et al. [2006]
1054 gam_L = (0.49*LAI) / ( SQRT(1.0+0.2*(LAI**2)) )
1055
1056 RETURN
1057 END SUBROUTINE GAMMA_LAI
1058
1059 !-------------------------------------------------------------------
1060 ! SUBROUTINE GAMMA_P
1061 ! Originally from Tan's gamma_etc.F
1062 !-------------------------------------------------------------------
1063 SUBROUTINE GAMMA_P( &
1064 DOY_in, tmidh, LAT, LONG, &
1065 PPFD, D_PPFD, gam_P &
1066 )
1067 !
1068 ! Description :
1069 !
1070 ! MEGAN biogenic emissions adjustment factor for
1071 ! photosynthetic photon flux density (PPFD or PAR)
1072 !
1073 ! Reference:
1074 !
1075 ! Estimates of global terrestial isoprene emissions using MEGAN
1076 ! (Model of Emissions of Gases and Aerosols from Nature )
1077 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer,
1078 ! P.I. Palmer, and C. Geron
1079 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006
1080 !
1081 ! Method:
1082 !
1083 ! GAMMA_P = 0.0 sin(a)<=0
1084 !
1085 ! GAMMA_P = sin(a)[2.46*0.9*PHI^3*(1+0.0005(Pdaily-400))]
1086 ! 0<a<180
1087 ! where PHI = above canopy PPFD transmission (non-dimension)
1088 ! Pdaily = daily average above canopy PPFD (umol/m2s)
1089 ! a = solar angle (degree)
1090 !
1091 ! Note: AAA = 2.46*BBB*PHI-0.9*PHI^2
1092 ! BBB = (1+0.0005(Pdaily-400))
1093 ! GAMMA_P = sin(a)*AAA
1094 !
1095 ! Pac
1096 ! PHI = -----------
1097 ! sin(a)*Ptoa
1098 !
1099 ! where Pac = above canopy PPFD (umol/m2s)
1100 ! Ptoa = PPFD at the top of atmosphere (umol/m2s)
1101 !
1102 ! Pac = SRAD * 4.766 mmmol/m2-s * 0.5
1103 !
1104 ! Ptoa = 3000 + 99*cos[2*3.14-( DOY-10)/365 )]
1105 ! where DOY = day of year (julian day)
1106 !
1107 ! NOTE: This code has been corrected. The gamma P equation as defined in the
1108 ! original Guenther et al., 2006 (equation 11b) is incorrect. This has the
1109 ! corrected algorithm. (CW, 08/16/2007)
1110 !-----------------------------------------------------------------
1111
1112 IMPLICIT NONE
1113
1114 INTEGER, INTENT(IN) :: DOY_in ! julian day at GMT
1115
1116 ! GMT hour plus minutes (in fractaionl hour) of the middle
1117 ! of the current time step
1118 REAL, INTENT(IN) :: tmidh
1119 REAL, INTENT(IN) :: LAT ! Latitude [=] degrees
1120 REAL, INTENT(IN) :: LONG ! Longitude [=] degrees
1121 REAL, INTENT(IN) :: PPFD ! Photosynthetic Photon Flus Density
1122 REAL, INTENT(IN) :: D_PPFD ! Daily PPFD
1123 REAL, INTENT(OUT) :: gam_P ! GAMMA_P
1124
1125
1126 !... Local scalars
1127 INTEGER :: DOY ! local julian day
1128 REAL :: HOUR ! solar hour
1129 REAL :: AAA, BBB
1130 REAL :: SIN_solarangle ! sin(solar angle)
1131 REAL :: Ptoa, Pac, PHI
1132
1133 ! End header ----------------------------------------------------
1134
1135 ! Convert time of the middle of the current time step
1136 ! from GMT to solar hour (include minutes in decimals)
1137 DOY = DOY_in
1138 HOUR = tmidh + long/15.
1139 IF ( HOUR .LT. 0.0 ) THEN
1140 HOUR = HOUR + 24.0
1141 DOY = DOY - 1
1142 ENDIF
1143
1144 ! Above canopy photosynthetic photo flux density (PPFD)
1145 ! ( micromole/m2/s )
1146 Pac = PPFD
1147
1148 ! Get sin of solar elevation angle
1149 CALL SOLARANGLE( DOY, HOUR, LAT, SIN_solarangle )
1150
1151 ! Calculate gamma_p in Eqn (10) of Guenther et al. [2006]
1152 IF ( SIN_solarangle .LE. 0.0 ) THEN
1153 ! Eqn (11a) of Guenther et al. [2006]
1154 gam_P = 0.0
1155 ELSE
1156 ! PPFD at top of the atmosphere
1157 ! Eqn (13) of Guenther et al. [2006]
1158 ! ( micromole/m2/s )
1159 Ptoa = 3000.0 + 99.0 * COS( 2.*3.14*(DOY-10.)/365. )
1160 ! Above canopy PPFD transmission
1161 ! Eqn (12) of Guenther et al. [2006]
1162 ! (nondimensional)
1163 PHI = Pac/(SIN_solarangle * Ptoa)
1164 ! Eqn (11b) of Guenther et al. [2006]
1165 ! (Note: typo in the paper; correction made 08/06/2007)
1166 BBB = 1. + 0.0005*( D_PPFD-400. )
1167 AAA = 2.46 * BBB * PHI - 0.9 * (PHI**2)
1168 gam_P = SIN_solarangle * AAA
1169
1170 ENDIF
1171 ! Screening the unforced errors
1172 ! IF solar elevation angle is less than 1 THEN
1173 ! gamma_p can not be greater than 0.1.
1174 IF (SIN_solarangle .LE. 0.0175 .AND. gam_P .GT. 0.1) THEN
1175 gam_P = 0.1
1176 ENDIF
1177
1178
1179 END SUBROUTINE GAMMA_P
1180
1181 ! ----------------------------------------------------------------
1182 ! SUBROUTINE GAMMA_A returns GAMMA_A
1183 ! Originally from Tan's gamma_etc.F
1184 !------------------------------------------------------------------
1185 SUBROUTINE GAMMA_A( i_spc, LAIp, LAIc, TSTLEN, D_TEMP, gam_A )
1186 ! Description :
1187 !
1188 ! MEGAN biogenic emissions adjustment factor for leaf age
1189 !
1190 ! Reference:
1191 !
1192 ! Estimates of global terrestial isoprene emissions using MEGAN
1193 ! (Model of Emissions of Gases and Aerosols from Nature )
1194 ! A. Guenther, T. Karl, P. Harley, C. Wiedinmyer,
1195 ! P.I. Palmer, and C. Geron
1196 ! Atmospheric Chemistry and Physics, 6, 3181-3210, 2006
1197 !
1198 ! MEGAN v2.0 Documentation
1199 !
1200 !
1201 ! Method:
1202 !
1203 ! GAMMA_age = Fnew*Anew + Fgro*Agro + Fmat*Amat + Fold*Aold
1204 ! where Fnew = new foliage fraction
1205 ! Fgro = growing foliage fraction
1206 ! Fmat = mature foliage fraction
1207 ! Fold = old foliage fraction
1208 ! Anew = relative emission activity for new foliage
1209 ! Agro = relative emission activity for growing foliage
1210 ! Amat = relative emission activity for mature foliage
1211 ! Aold = relative emission activity for old foliage
1212 !
1213 !
1214 ! For foliage fraction
1215 ! Case 1) LAIc = LAIp
1216 ! Fnew = 0.0 , Fgro = 0.1 , Fmat = 0.8 , Fold = 0.1
1217 !
1218 ! Case 2) LAIp > LAIc
1219 ! Fnew = 0.0 , Fgro = 0.0
1220 ! Fmat = 1-Fold
1221 ! Fold = (LAIp-LAIc)/LAIp
1222 !
1223 ! Case 3) LAIp < LAIc
1224 ! Fnew = 1-(LAIp/LAIc) t <= ti
1225 ! = (ti/t) * ( 1-(LAIp/LAIc) ) t > ti
1226 !
1227 ! Fmat = LAIp/LAIc t <= tm
1228 ! = (LAIp/LAIc) +
1229 ! ( (t-tm)/t ) * ( 1-(LAIp/LAIc) ) t > tm
1230 !
1231 ! Fgro = 1 - Fnew - Fmat
1232 ! Fold = 0.0
1233 !
1234 ! where
1235 ! ti = 5 + (0.7*(300-Tt)) Tt <= 303
1236 ! = 2.9 Tt > 303
1237 ! tm = 2.3*ti
1238 !
1239 ! t = length of the time step (days)
1240 ! ti = number of days between budbreak and the induction of
1241 ! emission
1242 ! tm = number of days between budbreak and the initiation of
1243 ! peak emissions rates
1244 ! Tt = average temperature (K) near top of the canopy during
1245 ! current time period (daily ave temp for this case)
1246 !
1247 !
1248 ! For relative emission activity
1249 ! Case 1) Constant
1250 ! Anew = 1.0 , Agro = 1.0 , Amat = 1.0 , Aold = 1.0
1251 !
1252 ! Case 2) Monoterpenes
1253 ! Anew = 2.0 , Agro = 1.8 , Amat = 0.95 , Aold = 1.0
1254 !
1255 ! Case 3) Sesquiterpenes
1256 ! Anew = 0.4 , Agro = 0.6 , Amat = 1.075, Aold = 1.0
1257 !
1258 ! Case 4) Methanol
1259 ! Anew = 3.0 , Agro = 2.6 , Amat = 0.85 , Aold = 1.0
1260 !
1261 ! Case 5) Isoprene
1262 ! Anew = 0.05 , Agro = 0.6 , Amat = 1.125, Aold = 1.0
1263
1264
1265 IMPLICIT NONE
1266
1267 ! SUBROUTINE arguments
1268
1269 !..."Pointer" for class of species
1270 INTEGER, INTENT(IN) :: i_spc
1271 !...average temperature of the previous 24-hours
1272 REAL, INTENT(IN) :: D_TEMP
1273 !...leaf area index of the current and previous
1274 !...month
1275 REAL, INTENT(IN) :: LAIp, LAIc
1276 !...time step between LAIc and LAIp (days)
1277 REAL, INTENT(IN) :: TSTLEN
1278 !...emissions adjustment factor accounting for leaf age
1279 REAL, INTENT(OUT) :: gam_A
1280
1281 ! Local scalars
1282
1283 !...leaf age fractions
1284 REAL :: Fnew, Fgro, Fmat, Fold
1285 !...relative emission activity index
1286 INTEGER :: AINDX
1287 !...time step between LAIC and LAIp (days)
1288 INTEGER :: t
1289 !...number of days between budbreak and the induction emission
1290 REAL ti
1291 !...number of days between budbreak and the initiation of peak
1292 !...emissions rates
1293 REAL tm
1294 !
1295
1296 REAL Tt ! average temperature (K)
1297 ! daily ave temp
1298
1299 ! End header ----------------------------------------------------
1300
1301 ! Choose relative emission activity class
1302 ! See Table 2 of MEGAN v2.0 Documentation
1303 !
1304
1305 IF ( (i_spc==imgn_acto) .OR. (i_spc==imgn_acta) .OR. (i_spc==imgn_form) &
1306 .OR. (i_spc==imgn_ch4) .OR. (i_spc==imgn_no) .OR. (i_spc==imgn_co) &
1307 ) THEN
1308 AINDX = 1
1309
1310 ELSE IF ( (i_spc==imgn_myrc) .OR. (i_spc==imgn_sabi) .OR. (i_spc==imgn_limo) &
1311 .OR. (i_spc==imgn_3car) .OR. (i_spc==imgn_ocim) .OR. (i_spc==imgn_bpin) &
1312 .OR. (i_spc==imgn_apin) .OR. ( i_spc==imgn_omtp) &
1313 ) THEN
1314 AINDX = 2
1315
1316 ELSE IF ( (i_spc==imgn_afarn) .OR. (i_spc==imgn_bcar) .OR. (i_spc==imgn_osqt) &
1317 ) THEN
1318 AINDX = 3
1319
1320 ELSE IF (i_spc==imgn_meoh) THEN
1321 aindx = 4
1322
1323 ELSE IF ( (i_spc==imgn_isop) .OR. (i_spc==imgn_mbo) ) THEN
1324 aindx = 5
1325 ELSE
1326 WRITE(mesg,fmt = '("Invalid i_spc in SUBROUTINE GAMMA_A; i_spc = ", I3)') i_spc
1327 CALL wrf_error_fatal(mesg)
1328 END IF
1329
1330
1331
1332 ! Time step between LAIp and LAIc (days)
1333 t = TSTLEN
1334 ! Tt is the average ambient air temperature (K) of the preceeding time
1335 ! interval. Here, use the monthly-mean surface air temperature
1336 Tt = D_TEMP
1337
1338 ! Calculate foliage fraction
1339 ! (section 3.2.2 of Guenther et al. [2006])
1340 IF (LAIp .EQ. LAIc) THEN
1341 Fnew = 0.0
1342 Fgro = 0.1
1343 Fmat = 0.8
1344 Fold = 0.1
1345 ELSEIF (LAIp .GT. LAIc) THEN
1346 Fnew = 0.0
1347 Fgro = 0.0
1348 Fold = ( LAIp-LAIc ) / LAIp
1349 Fmat = 1.0-Fold
1350 ELSE ! LAIp < LAIc
1351 ! Calculate ti, which is the number of days between budbreak and
1352 ! the induction of isoprene emission.
1353 IF (Tt .LE. 303.0) THEN
1354 ! Eqn (18a) of Guenther et al. [2006]
1355 ti = 5.0 + 0.7*(300.0-Tt)
1356 ELSE
1357 ! Eqn (18b) of Guenther et al. [2006]
1358 ti = 2.9
1359 ENDIF
1360 ! tm is the number of days between budbreak and the initiation
1361 ! of peak isoprene emissions rates.
1362 ! Eqn (19) of Guenther et al. [2006]
1363 tm = 2.3*ti
1364
1365 ! Calculate Fnew and Fmat, then Fgro and Fold
1366 ! Fnew
1367 IF (t .LE. ti) THEN
1368 ! Eqn (17a) of Guenther et al. [2006]
1369 Fnew = 1.0 - (LAIp/LAIc)
1370 ELSE
1371 ! Eqn (17b) of Guenther et al. [2006]
1372 Fnew = (ti/t) * ( 1-(LAIp/LAIc) )
1373 ENDIF
1374
1375 ! Fmat
1376 IF (t .LE. tm) THEN
1377 ! Eqn (17c) of Guenther et al. [2006]
1378 Fmat = LAIp/LAIc
1379 ELSE
1380 ! Eqn (17d) of Guenther et al. [2006]
1381 Fmat = (LAIp/LAIc) + ( (t-tm)/t ) * ( 1-(LAIp/LAIc) )
1382 ENDIF
1383
1384 Fgro = 1.0 - Fnew - Fmat
1385 Fold = 0.0
1386
1387 ENDIF
1388
1389 !Calculate GAMMA_A
1390 ! Anew, Agro, Amat, Aold are defined in module_data_megan2.F
1391 gam_A = Fnew*Anew(AINDX) + Fgro*Agro(AINDX) &
1392 + Fmat*Amat(AINDX) + Fold*Aold(AINDX)
1393
1394
1395 END SUBROUTINE GAMMA_A
1396
1397 ! ----------------------------------------------------------------
1398 ! SUBROUTINE SOLARANGLE calculates the solar angle
1399 ! Originally from Tan's solarangle.F
1400 !------------------------------------------------------------------
1401 SUBROUTINE SOLARANGLE( DAY, SHOUR, LAT, SIN_solarangle )
1402 !
1403 !
1404 ! Input:
1405 ! 1) Day of year
1406 ! 2) Latitude
1407 ! 3) Hour
1408 !
1409 ! Output: sin of solar angle
1410 !
1411
1412 IMPLICIT NONE
1413
1414 ! Arguments
1415 INTEGER, INTENT(IN) :: DAY ! DOY or julian day
1416 REAL, INTENT(IN) :: SHOUR ! Solar hour
1417 REAL, INTENT(IN) :: LAT ! Latitude
1418 REAL, INTENT(OUT) :: SIN_solarangle
1419
1420 ! Local scalars
1421 REAL :: sindelta, cosdelta, A, B
1422
1423 ! End header -----------------------------------------------------
1424
1425 sindelta = -SIN(0.40907) * COS( 6.28*(REAL(DAY,KIND(0.))+10.)/365. )
1426 cosdelta = SQRT(1.-sindelta**2.)
1427
1428 A = SIN( LAT*D2RAD ) * sindelta
1429 B = COS( LAT*D2RAD ) * cosdelta
1430
1431 SIN_solarangle = A + B * COS(2.*PI*(SHOUR-12.)/24.)
1432
1433
1434 END SUBROUTINE SOLARANGLE
1435
1436
1437
1438
1439 END SUBROUTINE bio_emissions_megan2
1440
1441 END MODULE module_bioemi_megan2