Adjusting include paths for removal of redundant code
[WRF.git] / chem / module_bioemi_megan2.F
bloba078447b7850d07dfb74c0d2a2a0eaee0275909a
1 MODULE module_bioemi_megan2
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   !
25 CONTAINS
27   SUBROUTINE bio_emissions_megan2(id,config_flags,ktau,dtstep,         &
28        curr_secs,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        ebio_c10h16,ebio_tol,ebio_bigalk, ebio_ch3oh,ebio_acet,         &
43        ebio_nh3,ebio_no2,ebio_c2h5oh,ebio_ch3cooh,ebio_mek,            &
44        ebio_bigene,ebio_c2h6,ebio_c2h4,ebio_c3h6,ebio_c3h8,ebio_so2,   &
45        ebio_dms,ebio_hcn,                                              &
46        ebio_c5h8,ebio_apinene,ebio_bpinene,ebio_toluene,               &
47        ebio_ch3cho,ebio_ch3co2h,ebio_tbut2ene,ebio_c2h5cho, &
48        ebio_nc4h10, &
49        ebio_sesq, ebio_mbo, ebio_bpi, ebio_myrc,                       &
50        ebio_alk3, ebio_alk4, ebio_alk5, ebio_ole1, ebio_ole2,          &    
51        ebio_aro1, ebio_aro2, ebio_ccho, ebio_meoh,                     &    
52        ebio_ethene, ebio_hcooh, ebio_terp, ebio_bald,                  &    
53        ebio_cco_oh, ebio_rco_oh,                                       &    
54        e_bio,                                                          &
55        ids,ide, jds,jde, kds,kde,                                      &
56        ims,ime, jms,jme, kms,kme,                                      &
57        its,ite, jts,jte, kts,kte                                       )
59     USE module_configure
60     USE module_state_description
61     USE module_data_megan2
62     USE module_data_mgn2mech
63 !    USE module_bioemi_beis313, ONLY : getpar, calc_zenithb
65     IMPLICIT NONE
67     ! Subroutine arguments
69     ! ...simulation parameters
70     TYPE(grid_config_rec_type),  INTENT(IN)    :: config_flags
72     ! ...domain id, current time step counter, xyz indices ..
73     INTEGER,   INTENT(IN   ) :: id,ktau,                               &
74          ids,ide, jds,jde, kds,kde,                                    &
75          ims,ime, jms,jme, kms,kme,                                    &
76          its,ite, jts,jte, kts,kte
78     ! ...current julian day
79     INTEGER, INTENT (IN) :: julday   
80     !...GTM hour of start of simulation, time step in seconds
81     REAL, INTENT(IN) :: gmt,dtstep
83     ! ...number of seconds into the simulation
84     REAL(KIND=8), INTENT(IN) :: curr_secs
86     ! ...3rd dimension size of array e_bio
87     INTEGER, INTENT (IN) :: ne_area
89     !...pressure (Pa)
90     REAL,  DIMENSION( ims:ime , kms:kme , jms:jme ),                   &
91          INTENT(IN) :: p_phy
93     !...latitude and longitude (degrees)
94     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
95          INTENT(IN) :: xlat, xlong
97     !... air density (kg air/m3)
98     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                      &
99          INTENT(IN) :: rho_phy
101     !...full layer height (m)
102    REAL,  DIMENSION( ims:ime , kms:kme , jms:jme )         ,           &
103           INTENT(IN) :: dz8w
105     !...2-meter temperature (K)
106     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
107          INTENT(IN) :: T2
109     !...downward shortwave surface flux (W/m2)
110     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
111          INTENT(IN) :: swdown                                    
113     !...Number of MEGAN v2.04 species as specified by the namelist 
114     !...variable nmegan; nmegan should equal n_spca_spc (this will
115     !...be checked later.)  Currently nmegan=n_spca_spc=138.
116     INTEGER, INTENT(IN) :: nmegan
117     
118     !...Emissions factors for nmegan=n_spca_spc=138 MEGAN v2.04 species
119     REAL, DIMENSION (ims:ime, jms:jme , nmegan) ,                      &
120          INTENT(INOUT) :: EFmegan
122     !...Emission factor for isoprene (read in from file
123     !...(wrfbiochemi_d<domain>)
124     !...(moles compound/km^2/hr)
125     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
126          INTENT(IN ) :: msebio_isop
128     !...Plant functional group percentage  (read in from file
129     !...(wrfbiochemi_d<domain>)
130     REAL, DIMENSION ( ims:ime , jms:jme ),                             &
131          INTENT(IN) ::                                                 &
132          pftp_bt, pftp_nt, pftp_sb, pftp_hb
134     !..."Climatological" Leaf area index  (read in from file
135     !...(wrfbiochemi_d<domain>)
136     REAL,  DIMENSION( ims:ime , jms:jme , 12 ),                        &
137          INTENT(IN) :: mlai
139     !..."Climatological" surface air temperature (K) (read in from file
140     !...(wrfbiochemi_d<domain>)
141     REAL,  DIMENSION( ims:ime , jms:jme , 12 ),                        &
142          INTENT(IN) :: mtsa
144     !..."Climatological" downward radiation (W/m2) (read in from file
145     !...(wrfbiochemi_d<domain>)
146     REAL, DIMENSION ( ims:ime , jms:jme , 12 ),                        &
147          INTENT(IN) :: mswdown
149     !...Actual emissions for a few selected species as diagnostics, using
150     !...MEGAN v2.0 classes of species classification
151     !...(mol km-2 hr-1)
152     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
153          INTENT(INOUT) ::                                              &
154          mebio_isop, mebio_apin, mebio_bpin, mebio_bcar,               &
155          mebio_acet, mebio_mbo, mebio_no
157     !...Actual biogenic emissions, converted to mechanisms species.
158     !...(ppm m/min)
159    REAL, DIMENSION( ims:ime, jms:jme, ne_area ),                       &
160          INTENT(INOUT ) :: e_bio
162     !...Actual biogenic emissions, converted to mechanisms species.
163     !...These variables were originally for BEIS3.11 biogenic emissions
164     !...modules.
165     !...(moles compound/km^2/hr)
166     REAL,  DIMENSION( ims:ime , jms:jme ),                             &
167          INTENT(INOUT  ) ::                                            &
168          ebio_iso,ebio_oli,ebio_api,ebio_lim,                          &
169          ebio_hc3,ebio_ete,ebio_olt,ebio_ket,ebio_ald,                 &
170          ebio_hcho,ebio_eth,ebio_ora2,ebio_co,ebio_no,                 &
171          ebio_c10h16,ebio_tol,ebio_bigalk, ebio_ch3oh,ebio_acet,       &
172          ebio_nh3,ebio_no2,ebio_c2h5oh,ebio_ch3cooh,ebio_mek,          &
173          ebio_bigene,ebio_c2h6,ebio_c2h4,ebio_c3h6,ebio_c3h8,ebio_so2, &
174          ebio_dms,ebio_hcn,                                            &
175          ebio_c5h8,ebio_apinene,ebio_bpinene,ebio_toluene,             &
176          ebio_ch3cho,ebio_ch3co2h,ebio_tbut2ene,ebio_c2h5cho,          &
177          ebio_nc4h10,                                                  &
178          ebio_sesq,ebio_mbo,ebio_bpi,ebio_myrc,                        &
179          ebio_alk3, ebio_alk4, ebio_alk5, ebio_ole1, ebio_ole2,        &    
180          ebio_aro1, ebio_aro2, ebio_ccho, ebio_meoh,                   &    
181          ebio_ethene, ebio_hcooh, ebio_terp, ebio_bald,                &    
182          ebio_cco_oh, ebio_rco_oh
184     !...Array of chemical concentrations
185     !...  in  - concentrations before biogenic emissions
186     !...  out - concentrations after biogeniec emissions
187     !... gas-phace concentrations are in ppm
188     REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),            &
189          INTENT(INOUT) :: chem
191     !...Current month
192     INTEGER, INTENT(IN) :: current_month
194     ! Local parameters
196     !...Below which set emissions rate to zero (mol km-2 hr-1)
197     REAL, PARAMETER :: min_emis = 0.001
198     
200     !...number of days in each month
201     INTEGER, PARAMETER :: DaysInMonth(12) = (/   &
202          31,28,31,30,31,30,31,31,30,31,30,31 /)
203     !...conversion between radians and degrees
204     REAL, PARAMETER :: PI = 3.14159 
205     REAL, PARAMETER :: D2RAD = PI/180.0 
208     ! Local Scalars
210     CHARACTER(len=256)   ::   mesg
211     INTEGER :: i,j,k,i_class, i_spc, icount, p_in_chem
212     INTEGER :: previous_month
214     !...minutes since start of run to the middle of the
215     !...current times step (seconds included as decimals)
216     REAL(KIND=8) :: xtime
218     !...the GMT hour of the middle of the current time step
219     !...(can be greater than 24)
220     INTEGER :: ixhour
221     REAL(KIND=8) :: xhour
223     !...minutes past the previous hour mark, at the
224     !...middle of the current time step
225     REAL :: xmin
227     !...the GMT hour of the middle of the current time step
228     !...(between 0 and 24)
229     REAL :: gmtp
231     !...GMT hour plus minutes (in fractaionl hour) of the middle
232     !...of the current time step
233     REAL :: tmidh
235     !...Current and previous month leaf area index
236     REAL :: LAIc, LAIp
238     !...temperature((K) and pressure (mb)
239     REAL :: tsa, tsa24, pres
241     !...latitude and longitude (degrees)
242     REAL :: lat, lon
244     !...downward solar radiation, current and some 24-hour mean (W/m2)
245     REAL :: swd, swd24
246     !...photosynthetic photon flux density (i.e. PPDF or PAR)  
247     !...(micromole m-2 s-1)
248     REAL :: par, par24, pardb, pardif
250     !...solar zenith angle (radians), cosine of zenith angle
251     REAL :: zen , coszen
253     !...days in the previous month
254     REAL :: tstlen
256     !...emissions factor (microgram m-2 hr-1)
257     REAL :: epsilon
259     !...MEGAN v2.04 emissions adjustment factors for leaf area, temperature,
260     !...light, leaf age, and soil moisture
261     !...(dimensionless)
262     REAL :: gam_LHT, gam_TMP, gam_PHO,gam_AGE, gam_SMT
264     !...normalized ratio accounting for production and loss within 
265     !...plant canopies (dimensionless)
266     REAL :: rho
268     !...Some light-dependent factor (dimensionless)
269     REAL :: ldf
271     !...conversion factor from mol km-2 hr-1 to ppm m min-1
272     REAL :: convert2
274     !...emission rate converted to mechanism species in mol km-2 hr-1
275     REAL :: gas_emis
277     ! Local Arrays
279     !...emissions adjustment factors for n_mgn_spc=20 classes of
280     !...MEGAN v2.04 specie.
281     !...adjust_factor = [GAMMA]*[rho] (see comments later)
282     !...(dimensionless)
283     REAL, DIMENSION(n_mgn_spc) :: adjust_factor
285     !...plant functional type fractions
286     REAL :: pft_frac(n_pft)
288     !...actually emissions rates of n_spca_spc=138 MEGAN v2.04 species
289     !...(mol km-2 hr-2)
290     REAL, DIMENSION ( n_spca_spc ) :: E_megan2
292     ! End header ------------------------------------------------------
295     ! MEGAN v2.04 has nmegan=n_spca_spc=138 species, which are lumped 
296     ! into n_mgn_spc=20 classes.  The number, names and indices of
297     ! these classes and species are defined in module_data_megan2.F.
298     ! They need to follow a few rules
299     
300     IF ( ktau .EQ. 1 ) THEN
302        ! The size of variable EFmegan(:,:,nmegan) is allocated based on
303        ! the value of namelist variable nmegan.  nmegan should be equal
304        ! to n_spca_spc (though can be greater than to n_spca_spc).
305        IF ( nmegan .NE. n_spca_spc ) THEN
306           WRITE(mesg,*)'namelist variable nmegan does not match n_spca_spc'
307           CALL wrf_error_fatal(mesg)          
308        END IF
310        ! For programming, the ordering of the  species or classes of
311        ! species should not matter, except that isoprene should always
312        ! be first; therefore, imgn_isop=1 and is_isoprene=1 always.
313        IF ( (imgn_isop .NE. 1) .OR. (is_isoprene .NE. 1) ) THEN
314           WRITE(mesg,*)'imgn_isop and is_isoprene in bio_emissions_megan should be 1'
315           CALL wrf_error_fatal(mesg)          
316        END IF
318     END IF
321     ! Initialize diagnostic output
322     ebio_iso  ( its:ite , jts:jte ) = 0.0
323     ebio_oli  ( its:ite , jts:jte ) = 0.0
324     ebio_api  ( its:ite , jts:jte ) = 0.0
325     ebio_lim  ( its:ite , jts:jte ) = 0.0
326     ebio_hc3  ( its:ite , jts:jte ) = 0.0
327     ebio_ete  ( its:ite , jts:jte ) = 0.0
328     ebio_olt  ( its:ite , jts:jte ) = 0.0
329     ebio_ket  ( its:ite , jts:jte ) = 0.0
330     ebio_ald  ( its:ite , jts:jte ) = 0.0
331     ebio_hcho ( its:ite , jts:jte ) = 0.0
332     ebio_eth  ( its:ite , jts:jte ) = 0.0
333     ebio_ora2 ( its:ite , jts:jte ) = 0.0
334     ebio_co   ( its:ite , jts:jte ) = 0.0
335     ebio_no   ( its:ite , jts:jte ) = 0.0
336     ebio_c10h16( its:ite , jts:jte ) = 0.0
337     ebio_tol  ( its:ite , jts:jte ) = 0.0
338     ebio_bigalk( its:ite , jts:jte ) = 0.0
339     ebio_ch3oh ( its:ite , jts:jte ) = 0.0
340     ebio_acet  ( its:ite , jts:jte ) = 0.0
341     ebio_nh3   ( its:ite , jts:jte ) = 0.0
342     ebio_no2   ( its:ite , jts:jte ) = 0.0
343     ebio_c2h5oh( its:ite , jts:jte ) = 0.0
344     ebio_ch3cooh( its:ite , jts:jte ) = 0.0
345     ebio_mek   ( its:ite , jts:jte ) = 0.0
346     ebio_bigene( its:ite , jts:jte ) = 0.0
347     ebio_c2h4  ( its:ite , jts:jte ) = 0.0
348     ebio_c2h6  ( its:ite , jts:jte ) = 0.0
349     ebio_c3h6  ( its:ite , jts:jte ) = 0.0
350     ebio_c3h8  ( its:ite , jts:jte ) = 0.0
351     ebio_so2   ( its:ite , jts:jte ) = 0.0
352     ebio_dms   ( its:ite , jts:jte ) = 0.0
353     ebio_terp  ( its:ite , jts:jte ) = 0.0
354     ebio_c5h8   ( its:ite , jts:jte ) = 0.0
355     ebio_apinene   ( its:ite , jts:jte ) = 0.0
356     ebio_bpinene   ( its:ite , jts:jte ) = 0.0
357     ebio_toluene   ( its:ite , jts:jte ) = 0.0
358     ebio_hcooh   ( its:ite , jts:jte ) = 0.0
359     ebio_ch3cho   ( its:ite , jts:jte ) = 0.0
360     ebio_c2h5oh   ( its:ite , jts:jte ) = 0.0
361     ebio_ch3co2h   ( its:ite , jts:jte ) = 0.0
362     ebio_tbut2ene   ( its:ite , jts:jte ) = 0.0
363     ebio_c2h5cho   ( its:ite , jts:jte ) = 0.0
364     ebio_nc4h10   ( its:ite , jts:jte ) = 0.0
365     ebio_sesq  ( its:ite , jts:jte ) = 0.0
366     ebio_mbo   ( its:ite , jts:jte ) = 0.0
367     ebio_bpi   ( its:ite , jts:jte ) = 0.0
368     ebio_myrc  ( its:ite , jts:jte ) = 0.0
369     e_bio     ( its:ite , jts:jte , 1:ne_area) = 0.0
370     
371     !...the following is redundant if there is no
372     !...bug in the subroutine
373     mebio_isop ( its:ite , jts:jte ) = 0.0
374     mebio_apin ( its:ite , jts:jte ) = 0.0
375     mebio_bpin ( its:ite , jts:jte ) = 0.0
376     mebio_bcar ( its:ite , jts:jte ) = 0.0
377     mebio_acet ( its:ite , jts:jte ) = 0.0 
378     mebio_mbo  ( its:ite , jts:jte ) = 0.0
379     mebio_no   ( its:ite , jts:jte ) = 0.0
382     ! Extract climatological values for relevant months.
383     !
384     !  In MEGAN v2.04, emissions rates dependent on ambient conditions
385     !  of the past 24 hours to the past month or so. The implementation
386     !  of MEGAN v2.04 here uses monthly-mean values of the previous
387     !  month for any past history required by the model.  The monthly-
388     !  -mean values should be provided as input in the
389     !  wrfbiochemi_d<domain> file.  Fully implementation (not done here)
390     !  require online calculations of 24-hour and 240-hour mean of
391     !  surface air temperature and downward PAR
392     !
393     !  MEGAN v2.04 also requires time-dependent leaf area index to
394     !  estimate leaf age.  Here, leaf area indices of the current
395     !  and the previous months are used.  The data should be
396     !  provided in wrfbiochemi_d<domain> file.
398     IF (current_month > 1) THEN
399        previous_month = current_month -1
400     ELSE
401        previous_month = 12
402     END IF
405     ! Following module_phot_fastj.F, determine current
406     ! time of day in GMT at the middle of the current 
407     ! time step, tmidh.
408     !     ktau  - time step counter
409     !     dstep - time step in seconds
410     !     gmt   - starting hour (in GMT) of the simulation
412     !...minutes since start of run to the middle of the
413     !...current times step (seconds included as decimals)
414     !(old way in r4 this will fail in about 2 yrs)...
415 !    xtime=(ktau-1)*dtstep/60. + dtstep/120.
416     xtime = curr_secs/60._8 + real(dtstep/120.,8)
417     !...the GMT hour of the middle of the current time step
418     !...(can be greater than 24)
419     ixhour = int(gmt + 0.01) + int(xtime/60._8)
420     xhour=real(ixhour,8)
421     !...minutes past the previous hour mark, at the
422     !...middle of the current time step
423     xmin = 60.*gmt + real(xtime-xhour*60._8,8)
424     !...the GMT hour of the middle of the current time step
425     !...(between 0 and 24)
426     gmtp=MOD(xhour,24._8)
427     !...GMT hour plus minutes (in fractaionl hour) of the middle
428     !...of the current time step
429     tmidh= gmtp + xmin/60.
431     WRITE(mesg,*) 'calculate MEGAN emissions at ktau, gmtp, tmidh = ',ktau, gmtp, tmidh
432     CALL wrf_message(mesg)
435     ! Get the mechanism converstion table
436     ! ( Even though the mechanism converstion table is time-independent,
437     ! do this for all time steps to be sure there will be no issue with 
438     ! restart runs.  This should be edited eventually to reduce 
439     ! redundant calculations.)
440     ! SHC  (11/08/2007)
441     GAS_MECH_SELECT1: SELECT CASE (config_flags%chem_opt)
442              
443     CASE (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_AQ, RADM2SORG_AQCHEM, RADM2SORG_KPP,GOCARTRADM2)
444        ! get p_of_radm2cbmz(:), p_of_radm2(:), and radm2_per_megan(:)
445        CALL get_megan2radm2_table
447     CASE (RACMSORG_AQ, RACMSORG_AQCHEM_KPP, RACM_ESRLSORG_AQCHEM_KPP, RACM_ESRLSORG_KPP, RACM_KPP, GOCARTRACM_KPP, RACMSORG_KPP, &
448           RACM_MIM_KPP, RACMPM_KPP)
449              
450        ! get p_of_megan2racm(:), p_of_racm(:), and racm_per_megan(:)
451        CALL get_megan2racm_table
452     CASE (RACM_SOA_VBS_KPP,RACM_SOA_VBS_AQCHEM_KPP,RACM_SOA_VBS_HET_KPP)
454         !get p_of_megan2racm(:), p_of_racm(:), and racm_per_megan(:)
455         CALL get_megan2racmSOA_table
457     CASE (CBMZ, CBMZ_BB, CBMZ_BB_KPP, CBMZ_MOSAIC_KPP, &
458           CBMZ_MOSAIC_4BIN, & 
459           CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ, &
460           CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, &
461           CBMZ_MOSAIC_DMS_4BIN_AQ, CBMZ_MOSAIC_DMS_8BIN_AQ, CBMZSORG, CBMZSORG_AQ, &
462           CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM3_AQ, CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM7_AQ)
463         
464        ! get p_of_megan2cbmz(:), p_of_cbmz(:), and cbmz_per_megan(:)
465        CALL get_megan2cbmz_table
467     CASE (CB05_SORG_AQ_KPP)
468        CALL get_megan2cb05_table
470     CASE ( CB05_SORG_VBS_AQ_KPP)
471        CALL get_megan2cb05vbs_table
473     CASE ( MOZART_KPP, MOZCART_KPP )
474        ! get p_of_megan2mozcart(:), p_of_mozcart(:), and mozcart_per_megan(:)
475        CALL get_megan2mozcart_table
476     CASE (  T1_MOZCART_KPP )
477        CALL get_megan2t1_mozc_table
478     CASE (  MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP )
479        CALL get_megan2mozm_table
481     CASE (SAPRC99_KPP,SAPRC99_MOSAIC_4BIN_VBS2_KPP, &
482          SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP,SAPRC99_MOSAIC_8BIN_VBS2_KPP)!BSINGH(12/03/2013): Added SAPRC 8 bin and non-aq on (04/07/2014) ! FIX FOR SAPRC07A
483        CALL get_megan2saprcnov_table
485     CASE ( CRIMECH_KPP, CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP )
486        ! get p_of_megan2crimech(:), p_of_crimech(:), and crimech_per_megan(:)
487        CALL get_megan2crimech_table
489     CASE DEFAULT
490        
491        CALL wrf_error_fatal('Species conversion table for MEGAN v2.04 not available. ')
493     END SELECT GAS_MECH_SELECT1
495     ! Calcuate biogenic emissions grid by grid
497     j_loop: DO j = jts, jte
498        i_loop: DO i = its, ite
501           ! Put variables of ambient conditions into scalar variables
503           tsa   = T2(i,j)                     ! air temperature at 2-meter (K)
504           pres  = 0.01*p_phy(i,kts,j)         ! surface pressure (mb)
505           lat   = xlat(i,j)                   ! latitude (degree) 
506           lon   = xlong(i,j)                  ! longitude (degress)
507           swd   = swdown(i,j)                 ! downward solar radiation (W/m2)
508           LAIc  = mlai(i,j,current_month)     ! current month leaf area index
509           LAIp  = mlai(i,j,previous_month)    ! previous month leaf area index
511           !...Emissions are dependent on the ambient conditions in the last
512           !...24 to 240 hours; here, use input data for monthly mean of the
513           !...the previous month
514           tsa24 = mtsa    (i,j,previous_month) ! [=]K
515           swd24 = mswdown (i,j,previous_month) ! [=] W/m2
517           !...Perform checks on max and min bounds for temperature
518           IF (tsa .LT. 200.0) THEN
519              WRITE (mesg,'("temperature too low at i=",i3," ,j=",i3 )')i,j
520              CALL wrf_message(mesg)
521           END IF
522           IF (tsa .GT. 315.0 ) THEN
523              WRITE (mesg,'("temperature too high at i=",i3," ,j=",i3," ;resetting to 315K" )')i,j
524              CALL wrf_message(mesg)
525              tsa = 315.0
526           END IF
528 !          !...Calculate zenith angle (in radians)
529 !          !...(NOTE: nonstandard longitude input here: >0 for W, <0 for E!!)
530 !          !...(subroutine calc_zenithb is in module_bioemis_beis313.F)
531 !          CALL calc_zenithb(lat,-lon,julday,tmidh,zen)
532 !          coszen = COS(zen)
534           !....Convert downward solar radiation to photosynthetically
535           !....active radiation
537 !          !......(subroutine getpar is in module_bioemis_beis313.F)
538 !          CALL getpar( swd, pres, zen, pardb, pardif )
539 !          par = pardb + pardif ! micro-mole/m2/s
541           !......assume 4.766 (umol m-2 s-1) per (W m-2)
542           !......assume 1/2 of swd is in 400-700 nm band
543           par = 4.766 * 0.5 * swd
545           !......Check max/min bounds of PAR
546           IF ( par .LT. 0.00 .OR. par .GT. 2600.0 ) THEN
547              WRITE (mesg,'("par out of range at i=",i3," ,j=",i3," par =",f8.3 )')i,j,par
548              CALL wrf_message(mesg)
549           END IF
551           !......For the 24-avg PAR, 
552           !......assume 4.766 (umol m-2 s-1) per (W m-2)
553           !......assume 1/2 of swd is in 400-700 band
554           par24 = swd24 * 4.766 * 0.5
556           ! ------------------------------------------------------------
557           !
558           !  MEGAN v2.04 Box Model
559           !
560           !  Reference: 
561           !
562           !    Estimates of global terrestial isoprene emissions using MEGAN
563           !    (Model of Emissions of Gases and Aerosols from Nature )
564           !    A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, 
565           !    P.I. Palmer, and C. Geron
566           !    Atmospheric Chemistry and Physics, 6, 3181-3210, 2006      
567           !
568           !    MEGAN v2.0 Documentation
569           !
570           !
571           !  The following code is based on Tan's megan.F dated 11/21/2006
572           !
573           !   Scientific algorithm
574           !
575           !           Emission = [epsilon][gamma][rho]
576           !
577           !         where [epsilon] = emission factor (usually um m-2 hr or mole km-2 hr-1)
578           !               [gamma]   = emission activity factor (dimensionless)
579           !               [rho]     = production and loss within plant canopies
580           !                         (dimensionless)
581           !
582           !            [gamma]  = [gamma_CE][gamma_age][gamma_SM]
583           !
584           !         where [gamma_CE]  = canopy correction factor
585           !               [gamma_age] = leaf age correction factor
586           !               [gamma_SM]  = soil moisture correction factor
587           !
588           !            [gamma_CE] = [gamma_LAI][gamma_T]((1-LDF) + LDF*[gamma_P] )
589           !                      
590           !         where [gamma_LAI] = leaf area index factor
591           !               [gamma_P]   = PPFD emission activity factor
592           !               [gamma_T]   = temperature response factor
593           !               LDF         = 
594           !
595           !           Emission = [epsilon][gamma_LAI][gamma_T][gamma_age]
596           !                      x ((1-LDF) + LDF*[gamma_P] )[rho]
597           !
598           !     or
599           ! 
600           !           Emission = adjust_factor [epsilon]
601           !
602           !         where
603           !
604           !           adjust_fact = [gamma_LAI][gamma_T][gamma_age]((1-LDF) + LDF*[gamma_P] )[rho]
607           !
608           ! Calculate the dimensionless emission adjustment factor. 
609           ! MEGAN v2.04 has n_spca_spc = 138 species.  These species are 
610           ! lumped into n_mgn_spc=20 classes.  The emission adjustment
611           ! factors are different for the 20 classes of species.
612           !
613           ! NOTE: This version of the code contains the corrected equation for 
614           ! gamma P (based on a revised version of equation 11b from Guenther et al., 2006)
615           ! CW (08/16/2007)
616           !
617           ! NOTE: This version of the code contains the updated emission factors (static)
618           ! and beta values based on Alex's V2.1 notes sent out on August 13, 2007
619           ! CW (08/16/2007)
620           ! 
621           ! NOTE: This version of the code applies the same gamma T equation to the 
622           ! emissions of all compounds other than isoprene. This occurs regardless 
623           ! of whether the emissions are light dependent or not. This is NOT the same 
624           ! as what Alex has in his code. In his code, the light-dependent emissions
625           ! are also given the isoprene gamma T. Because all emissions (other than isoprene
626           ! are assigned the same gamma T, this could lead to overestimates of emissions
627           ! at high temperatures (>40C). Light-dependent emisisons (e.g., SQTs) should
628           ! fall off at high temperatures. (CW, 08/16/2007)
630           !...Calculate adjustment factor components that are species-independent
632           !......Get the leaf area index factor gam_LHT
633           CALL GAMMA_LAI( LAIc, gam_LHT)
635           !......Get the light emission activity factor gam_PHO
636           CALL GAMMA_P( julday, tmidh, lat, lon, par, par24, gam_PHO )
638           !......Get the soil moisture factor gam_SMT
639           !......for now, set = 1.0
640           gam_SMT = 1.0
642           !...Calculate the overall emissions adjustment factors, for
643           !...each of the n_mgn_spc=20 classes of compounds
645           DO i_class = 1, n_mgn_spc
646              
647              ! Get the temperature response factor gam_TMP
648              !  One algorithm for isoprene, and one for non-isoprene
650              IF ( i_class == imgn_isop ) THEN
651                 CALL GAMMA_TISOP( tsa, tsa24, gam_TMP )
652              ELSE
653                 CALL GAMMA_TNISP( i_class , tsa, gam_TMP  )
654              END IF
656              ! Get the leaf age correction factor gam_AGE
658              !...Time step (days) between LAIc and LAIp:
659              !...Since monthly mean LAI is used,
660              !...use # of days in the previous month
661              tstlen = REAL(DaysInMonth(previous_month),KIND(1.0))
663              CALL GAMMA_A( i_class , LAIp, LAIc, TSTLEN, tsa24, gam_AGE )
665              ! rho - normalized ratio accounting for production and
666              ! oss within plant canopies; rho_fct is defined in 
667              ! module_data_megan2.F; currently rho_fct = 1.0 for all 
668              ! species (dimensionless)
669              rho = rho_fct(i_class)
671              ! Fraction of emission to apply light-dependence factor
672              ! ldf_fct is defined in module_data_megan2.F
673              ! (dimensionless)
674              ldf = ldf_fct(i_class)
676              ! The overall emissions adjustment factor
677              ! (dimensionless)
679              adjust_factor(i_class) = gam_TMP * gam_AGE * gam_LHT * gam_SMT * rho * &
680                   ( (1.0-LDF) + gam_PHO*LDF )
682           END DO !i_class = 1, n_mgn_spc (loop over classes of MEGAN species )
685           ! For isoprene, the emission factor is already read in from
686           ! wrfbiochemi_d<domain> file; therefore, actual emissions rate
687           ! can be calculated here already.
688           ! (mol km-2 hr-1)
689           E_megan2(is_isoprene) = adjust_factor(imgn_isop)*msebio_isop(i,j)
690           IF ( E_megan2(is_isoprene) .LT. min_emis ) E_megan2(is_isoprene)=0.
693           ! Calculate emissions for all n_spca_spc=nmegan=138 MEGAN v2.04
694           ! species, except for isoprene.  For non-isoprene emissions,
695           ! the emission factor [epsilon] has to be calculated
696           ! for the first time step.
698           !...Loop over species, because emission factor [epsilon] is 
699           !...different for each species
700           !...( i_spc = 1 is skipped in the do loop below to skip 
701           !...isoprene; this works because is_isoprene = 1 )
702           DO i_spc = 2, n_spca_spc 
704              ! The lumped class in which the current species is a member
705              i_class = mg20_map (i_spc)
707              ! Calculate emission factor (microgram m-2 hr-1) for 
708              ! species i_spc 
709              ! ( Even though EFmegan is time invariant, for now calculate
710              ! EFmegan for every time step to be sure there will be
711              ! no issue with restart runs.
712              ! SHC  (11/08/2007)
713 !             IF ( ktau .EQ. 1 ) THEN
715                 ! Grab plant functional type fractions for current grid
716                 ! cell (pftp_* is the plant functional type % and was
717                 ! read in from wrfbiochemi_d<domain> file.)
718                 pft_frac(k_bt) = 0.01*pftp_bt(i,j)
719                 pft_frac(k_nt) = 0.01*pftp_nt(i,j)
720                 pft_frac(k_sb) = 0.01*pftp_sb(i,j)
721                 pft_frac(k_hb) = 0.01*pftp_hb(i,j)
723                 ! Sum up emissions factor over plant functional types
724                 epsilon = 0.0
725                 DO k = 1, n_pft !loop over plant functional types
726                    epsilon = epsilon +                             &
727                         pft_frac(k)*EF(i_class,k)*EF_frac(i_spc,k)
728                 END DO
730                 ! Store emission factor to variable EFmegan (which is
731                 ! declared in Registry/registry.chem)
732                 ! (migrogram m-2 hr-1)
733                 EFmegan(i,j,i_spc) = epsilon
735 !             END IF ! ( ktau .EQ. 1 )
737              ! Calculate actual emission rate for species i_spc;
738              ! also, convert units from (microgram m-2 hr-1) to 
739              ! (mol km-2 hr-1)
740              E_megan2(i_spc) = EFmegan(i,j,i_spc)*        &
741                   adjust_factor(i_class)/spca_mwt(i_spc)
742              IF ( E_megan2(i_spc) .LT. min_emis ) E_megan2(i_spc)=0.
744           END DO !i_spc = 2, n_spca_spc, loop over all non-isoprene species
747           ! Output emissions for some species as diagnostics
748           ! (mol km-2 hr-1)
749 !          print*,'is_isoprene',is_isoprene
750 !          print*,'is_pinene_a',is_pinene_a
751 !          print*,'is_pinene_b',is_pinene_b
753 !          if (E_megan2 (is_isoprene).gt.0) print*,'E_megan2 (is_isoprene)',E_megan2 (is_isoprene)
754 !          if (E_megan2 (is_pinene_a).gt.0) print*,'E_megan2 (is_pinene_a)',E_megan2 (is_pinene_a)
756           mebio_isop  (i,j) = E_megan2 ( is_isoprene        )
757           mebio_apin  (i,j) = E_megan2 ( is_pinene_a        )
758           mebio_bpin  (i,j) = E_megan2 ( is_pinene_b        )
759           mebio_bcar  (i,j) = E_megan2 ( is_caryophyllene_b )
760           mebio_acet  (i,j) = E_megan2 ( is_acetone         )
761           mebio_mbo   (i,j) = E_megan2 ( is_MBO_2m3e2ol     )
762           mebio_no    (i,j) = E_megan2 ( is_nitric_OXD      )
765           ! Speciate the n_spca_spc=nmegan=138 species into
766           ! the gas-phase mechanism species
769           !...conversion factor from mol km-2 hr-1 to ppm m min-1
770           !...(e_bio is in units of ppm m min-1)
771           convert2 = 0.02897/(rho_phy(i,kts,j)*60.)
774           !...
775           GAS_MECH_SELECT: SELECT CASE (config_flags%chem_opt)
777           CASE ( MOZART_KPP, MOZCART_KPP )
779              DO icount = 1, n_megan2mozcart
780 !-----------------------------------------------------------------------
781 ! Get index to chem array for the corresponding MOZCART species.  
782 !-----------------------------------------------------------------------
783                 p_in_chem = p_of_mozcart(icount)
784 use_megan_emission : &
785                 IF ( p_in_chem /= non_react ) THEN
786 !-----------------------------------------------------------------------
787 ! Check if the species is actually in the mechanism
788 !-----------------------------------------------------------------------
789 is_mozcart_species : &
790                    IF ( p_in_chem >= param_first_scalar ) THEN
791 !-----------------------------------------------------------------------
792 ! Emission rate for mechanism species in mol km-2 hr-1
793 !-----------------------------------------------------------------------
794                       gas_emis = mozcart_per_megan(icount) * E_megan2(p_of_megan2mozcart(icount))
795 !-----------------------------------------------------------------------
796 ! Add emissions to diagnostic output variables.
797 ! ebio_xxx (mol km-2 hr-1) were originally used by the 
798 ! BEIS3.11 biogenic emissions module. 
799 ! I have also borrowed variable e_bio (ppm m min-1).
800 !-----------------------------------------------------------------------
801                       IF ( p_in_chem == p_isopr ) THEN
802                          ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
803                          e_bio(i,j,p_isopr-1)   = e_bio(i,j,p_isopr-1)  + gas_emis*convert2
804                       ELSE IF ( p_in_chem == p_no ) THEN
805                          ebio_no(i,j)  = ebio_no(i,j) + gas_emis
806                          e_bio(i,j,p_no-1)   = e_bio(i,j,p_no-1)  + gas_emis*convert2
807                       ELSE IF ( p_in_chem == p_no2 ) THEN
808                          ebio_no2(i,j)  = ebio_no2(i,j) + gas_emis
809                          e_bio(i,j,p_no2-1)   = e_bio(i,j,p_no2-1)  + gas_emis*convert2
810                       ELSE IF ( p_in_chem == p_co ) THEN
811                          ebio_co(i,j)  = ebio_co(i,j) + gas_emis
812                          e_bio(i,j,p_co-1)   = e_bio(i,j,p_co-1)  + gas_emis*convert2
813                       ELSE IF ( p_in_chem == p_hcho ) THEN
814                          ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
815                          e_bio(i,j,p_hcho-1)   = e_bio(i,j,p_hcho-1)  + gas_emis*convert2
816                       ELSE IF ( p_in_chem == p_ald ) THEN
817                          ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
818                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
819                       ELSE IF ( p_in_chem == p_acet ) THEN
820                          ebio_acet(i,j) = ebio_acet(i,j) + gas_emis
821                          e_bio(i,j,p_acet-1)   = e_bio(i,j,p_acet-1)  + gas_emis*convert2
822                       ELSE IF ( p_in_chem == p_tol ) THEN
823                          ebio_tol(i,j) = ebio_tol(i,j) + gas_emis
824                          e_bio(i,j,p_tol-1)   = e_bio(i,j,p_tol-1)  + gas_emis*convert2
825                       ELSE IF ( p_in_chem == p_c10h16 ) THEN
826                          ebio_c10h16(i,j) = ebio_c10h16(i,j) + gas_emis
827                          e_bio(i,j,p_c10h16-1)   = e_bio(i,j,p_c10h16-1)  + gas_emis*convert2
828                       ELSE IF ( p_in_chem == p_so2 ) THEN
829                          ebio_so2(i,j) = ebio_so2(i,j) + gas_emis
830                          e_bio(i,j,p_so2-1)   = e_bio(i,j,p_so2-1)  + gas_emis*convert2
831                       ELSE IF ( p_in_chem == p_dms ) THEN
832                          ebio_dms(i,j) = ebio_dms(i,j) + gas_emis
833                          e_bio(i,j,p_dms-1)   = e_bio(i,j,p_dms-1)  + gas_emis*convert2
834                       ELSE IF ( p_in_chem == p_bigalk ) THEN
835                          ebio_bigalk(i,j) = ebio_bigalk(i,j) + gas_emis
836                          e_bio(i,j,p_bigalk-1)   = e_bio(i,j,p_bigalk-1)  + gas_emis*convert2
837                       ELSE IF ( p_in_chem == p_bigene ) THEN
838                          ebio_bigene(i,j) = ebio_bigene(i,j) + gas_emis
839                          e_bio(i,j,p_bigene-1)   = e_bio(i,j,p_bigene-1)  + gas_emis*convert2
840                       ELSE IF ( p_in_chem == p_nh3 ) THEN
841                          ebio_nh3(i,j) = ebio_nh3(i,j) + gas_emis
842                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
843                       ELSE IF ( p_in_chem == p_ch3oh ) THEN
844                          ebio_ch3oh(i,j) = ebio_ch3oh(i,j) + gas_emis
845                          e_bio(i,j,p_ch3oh-1)   = e_bio(i,j,p_ch3oh-1)  + gas_emis*convert2
846                       ELSE IF ( p_in_chem == p_c2h5oh ) THEN
847                          ebio_c2h5oh(i,j) = ebio_c2h5oh(i,j) + gas_emis
848                          e_bio(i,j,p_c2h5oh-1)   = e_bio(i,j,p_c2h5oh-1)  + gas_emis*convert2
849                       ELSE IF ( p_in_chem == p_ch3cooh ) THEN
850                          ebio_ch3cooh(i,j) = ebio_ch3cooh(i,j) + gas_emis
851                          e_bio(i,j,p_ch3cooh-1)   = e_bio(i,j,p_ch3cooh-1)  + gas_emis*convert2
852                       ELSE IF ( p_in_chem == p_mek ) THEN
853                          ebio_mek(i,j) = ebio_mek(i,j) + gas_emis
854                          e_bio(i,j,p_mek-1)   = e_bio(i,j,p_mek-1)  + gas_emis*convert2
855                       ELSE IF ( p_in_chem == p_c2h4 ) THEN
856                          ebio_c2h4(i,j) = ebio_c2h4(i,j) + gas_emis
857                          e_bio(i,j,p_c2h4-1)   = e_bio(i,j,p_c2h4-1)  + gas_emis*convert2
858                       ELSE IF ( p_in_chem == p_c2h6 ) THEN
859                          ebio_c2h6(i,j) = ebio_c2h6(i,j) + gas_emis
860                          e_bio(i,j,p_c2h6-1)   = e_bio(i,j,p_c2h6-1)  + gas_emis*convert2
861                       ELSE IF ( p_in_chem == p_c3h6 ) THEN
862                          ebio_c3h6(i,j) = ebio_c3h6(i,j) + gas_emis
863                          e_bio(i,j,p_c3h6-1)   = e_bio(i,j,p_c3h6-1)  + gas_emis*convert2
864                       ELSE IF ( p_in_chem == p_c3h8 ) THEN
865                          ebio_c3h8(i,j) = ebio_c3h8(i,j) + gas_emis
866                          e_bio(i,j,p_c3h8-1)   = e_bio(i,j,p_c3h8-1)  + gas_emis*convert2
867                       END IF
868                    END IF is_mozcart_species
869                 END IF use_megan_emission
870              END DO
872           CASE ( T1_MOZCART_KPP )
874              DO icount = 1, n_megan2t1_mozc
875 !-----------------------------------------------------------------------
876 ! Get index to chem array for the corresponding T1_MOZCART species.  
877 !-----------------------------------------------------------------------
878                 p_in_chem = p_of_t1_mozc(icount)
879 use_megan_emis_a : &
880                 IF ( p_in_chem /= non_react ) THEN
881 !-----------------------------------------------------------------------
882 ! Check if the species is actually in the mechanism
883 !-----------------------------------------------------------------------
884 is_t1_mozc_species : &
885                    IF ( p_in_chem >= param_first_scalar ) THEN
886 !-----------------------------------------------------------------------
887 ! Emission rate for mechanism species in mol km-2 hr-1
888 !-----------------------------------------------------------------------
889                       gas_emis = t1_mozc_per_megan(icount) * E_megan2(p_of_megan2t1_mozc(icount))
890 !-----------------------------------------------------------------------
891 ! Add emissions to diagnostic output variables.
892 ! ebio_xxx (mol km-2 hr-1) were originally used by the 
893 ! BEIS3.11 biogenic emissions module. 
894 ! I have also borrowed variable e_bio (ppm m min-1).
895 !-----------------------------------------------------------------------
896                       IF ( p_in_chem == p_isopr ) THEN
897                          ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
898                          e_bio(i,j,p_isopr-1)   = e_bio(i,j,p_isopr-1)  + gas_emis*convert2
899                       ELSE IF ( p_in_chem == p_apin ) THEN
900                          ebio_api(i,j) = ebio_api(i,j) + gas_emis
901                          e_bio(i,j,p_apin-1)   = e_bio(i,j,p_apin-1)  + gas_emis*convert2
902                       ELSE IF ( p_in_chem == p_bpin ) THEN
903                          ebio_bpi(i,j) = ebio_bpi(i,j) + gas_emis
904                          e_bio(i,j,p_bpin-1)   = e_bio(i,j,p_bpin-1)  + gas_emis*convert2
905                       ELSE IF ( p_in_chem == p_limon ) THEN
906                          ebio_lim(i,j) = ebio_lim(i,j) + gas_emis
907                          e_bio(i,j,p_limon-1)   = e_bio(i,j,p_limon-1)  + gas_emis*convert2
908                       ELSE IF ( p_in_chem == p_myrc ) THEN
909                          ebio_myrc(i,j) = ebio_myrc(i,j) + gas_emis
910                          e_bio(i,j,p_myrc-1)   = e_bio(i,j,p_myrc-1)  + gas_emis*convert2
911                       ELSE IF ( p_in_chem == p_bcary ) THEN
912                          ebio_sesq(i,j) = ebio_sesq(i,j) + gas_emis
913                          e_bio(i,j,p_bcary-1)   = e_bio(i,j,p_bcary-1)  + gas_emis*convert2
914                       ELSE IF ( p_in_chem == p_mbo ) THEN
915                          ebio_mbo(i,j) = ebio_mbo(i,j) + gas_emis
916                          e_bio(i,j,p_mbo-1)   = e_bio(i,j,p_mbo-1)  + gas_emis*convert2
917                       ELSE IF ( p_in_chem == p_ch3oh ) THEN
918                          ebio_ch3oh(i,j) = ebio_ch3oh(i,j) + gas_emis
919                          e_bio(i,j,p_ch3oh-1)   = e_bio(i,j,p_ch3oh-1)  + gas_emis*convert2
920                       ELSE IF ( p_in_chem == p_c2h5oh ) THEN
921                          ebio_c2h5oh(i,j) = ebio_c2h5oh(i,j) + gas_emis
922                          e_bio(i,j,p_c2h5oh-1)   = e_bio(i,j,p_c2h5oh-1)  + gas_emis*convert2
923                       ELSE IF ( p_in_chem == p_hcho ) THEN
924                          ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
925                          e_bio(i,j,p_hcho-1)   = e_bio(i,j,p_hcho-1)  + gas_emis*convert2
926                       ELSE IF ( p_in_chem == p_ald ) THEN
927                          ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
928                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
929                       ELSE IF ( p_in_chem == p_ch3cooh ) THEN
930                          ebio_ch3cooh(i,j) = ebio_ch3cooh(i,j) + gas_emis
931                          e_bio(i,j,p_ch3cooh-1)   = e_bio(i,j,p_ch3cooh-1)  + gas_emis*convert2
932                       ELSE IF ( p_in_chem == p_hcooh ) THEN
933                          ebio_hcooh(i,j) = ebio_hcooh(i,j) + gas_emis
934                          e_bio(i,j,p_hcooh-1)   = e_bio(i,j,p_hcooh-1)  + gas_emis*convert2
935                       ELSE IF ( p_in_chem == p_hcn ) THEN
936                          ebio_hcn(i,j) = ebio_hcn(i,j) + gas_emis
937                          e_bio(i,j,p_hcn-1)   = e_bio(i,j,p_hcn-1)  + gas_emis*convert2
938                       ELSE IF ( p_in_chem == p_nh3 ) THEN
939                          ebio_nh3(i,j) = ebio_nh3(i,j) + gas_emis
940                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
941                       ELSE IF ( p_in_chem == p_co ) THEN
942                          ebio_co(i,j)  = ebio_co(i,j) + gas_emis
943                          e_bio(i,j,p_co-1)   = e_bio(i,j,p_co-1)  + gas_emis*convert2
944                       ELSE IF ( p_in_chem == p_c2h4 ) THEN
945                          ebio_c2h4(i,j) = ebio_c2h4(i,j) + gas_emis
946                          e_bio(i,j,p_c2h4-1)   = e_bio(i,j,p_c2h4-1)  + gas_emis*convert2
947                       ELSE IF ( p_in_chem == p_c2h6 ) THEN
948                          ebio_c2h6(i,j) = ebio_c2h6(i,j) + gas_emis
949                          e_bio(i,j,p_c2h6-1)   = e_bio(i,j,p_c2h6-1)  + gas_emis*convert2
950                       ELSE IF ( p_in_chem == p_c3h6 ) THEN
951                          ebio_c3h6(i,j) = ebio_c3h6(i,j) + gas_emis
952                          e_bio(i,j,p_c3h6-1)   = e_bio(i,j,p_c3h6-1)  + gas_emis*convert2
953                       ELSE IF ( p_in_chem == p_c3h8 ) THEN
954                          ebio_c3h8(i,j) = ebio_c3h8(i,j) + gas_emis
955                          e_bio(i,j,p_c3h8-1)   = e_bio(i,j,p_c3h8-1)  + gas_emis*convert2
956                       ELSE IF ( p_in_chem == p_bigalk ) THEN
957                          ebio_bigalk(i,j) = ebio_bigalk(i,j) + gas_emis
958                          e_bio(i,j,p_bigalk-1)   = e_bio(i,j,p_bigalk-1)  + gas_emis*convert2
959                       ELSE IF ( p_in_chem == p_bigene ) THEN
960                          ebio_bigene(i,j) = ebio_bigene(i,j) + gas_emis
961                          e_bio(i,j,p_bigene-1)   = e_bio(i,j,p_bigene-1)  + gas_emis*convert2
962                       ELSE IF ( p_in_chem == p_tol ) THEN
963                          ebio_tol(i,j) = ebio_tol(i,j) + gas_emis
964                          e_bio(i,j,p_tol-1)   = e_bio(i,j,p_tol-1)  + gas_emis*convert2
965                       END IF
966                    END IF is_t1_mozc_species
967                 END IF use_megan_emis_a
968              END DO
970           CASE ( MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP )
972              DO icount = 1, n_megan2mozm
973 !-----------------------------------------------------------------------
974 ! Get index to chem array for the corresponding MOZCART species.  
975 !-----------------------------------------------------------------------
976                 p_in_chem = p_of_mozm(icount)
977 use_megan_emis : &
978                 IF ( p_in_chem /= non_react ) THEN
979 !-----------------------------------------------------------------------
980 ! Check if the species is actually in the mechanism
981 !-----------------------------------------------------------------------
982 is_mozm_species : &
983                    IF ( p_in_chem >= param_first_scalar ) THEN
984 !-----------------------------------------------------------------------
985 ! Emission rate for mechanism species in mol km-2 hr-1
986 !-----------------------------------------------------------------------
987                       gas_emis = mozm_per_megan(icount) * E_megan2(p_of_megan2mozm(icount))
988 !-----------------------------------------------------------------------
989 ! Add emissions to diagnostic output variables.
990 ! ebio_xxx (mol km-2 hr-1) were originally used by the 
991 ! BEIS3.11 biogenic emissions module. 
992 ! I have also borrowed variable e_bio (ppm m min-1).
993 !-----------------------------------------------------------------------
994                       IF ( p_in_chem == p_isopr ) THEN
995                          ebio_iso(i,j) = ebio_iso(i,j) + gas_emis
996                          e_bio(i,j,p_isopr-1)   = e_bio(i,j,p_isopr-1)  + gas_emis*convert2
997                       ELSE IF ( p_in_chem == p_no ) THEN
998                          ebio_no(i,j)  = ebio_no(i,j) + gas_emis
999                          e_bio(i,j,p_no-1)   = e_bio(i,j,p_no-1)  + gas_emis*convert2
1000                       ELSE IF ( p_in_chem == p_no2 ) THEN
1001                          ebio_no2(i,j)  = ebio_no2(i,j) + gas_emis
1002                          e_bio(i,j,p_no2-1)   = e_bio(i,j,p_no2-1)  + gas_emis*convert2
1003                       ELSE IF ( p_in_chem == p_co ) THEN
1004                          ebio_co(i,j)  = ebio_co(i,j) + gas_emis
1005                          e_bio(i,j,p_co-1)   = e_bio(i,j,p_co-1)  + gas_emis*convert2
1006                       ELSE IF ( p_in_chem == p_hcho ) THEN
1007                          ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
1008                          e_bio(i,j,p_hcho-1)   = e_bio(i,j,p_hcho-1)  + gas_emis*convert2
1009                       ELSE IF ( p_in_chem == p_ald ) THEN
1010                          ebio_ald(i,j) = ebio_ald(i,j) + gas_emis
1011                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
1012                       ELSE IF ( p_in_chem == p_acet ) THEN
1013                          ebio_acet(i,j) = ebio_acet(i,j) + gas_emis
1014                          e_bio(i,j,p_acet-1)   = e_bio(i,j,p_acet-1)  + gas_emis*convert2
1015                       ELSE IF ( p_in_chem == p_tol ) THEN
1016                          ebio_tol(i,j) = ebio_tol(i,j) + gas_emis
1017                          e_bio(i,j,p_tol-1)   = e_bio(i,j,p_tol-1)  + gas_emis*convert2
1018                       ELSE IF ( p_in_chem == p_apin ) THEN
1019                          ebio_api(i,j) = ebio_api(i,j) + gas_emis
1020                          e_bio(i,j,p_apin-1)   = e_bio(i,j,p_apin-1)  + gas_emis*convert2
1021                       ELSE IF ( p_in_chem == p_bpin ) THEN
1022                          ebio_bpi(i,j) = ebio_bpi(i,j) + gas_emis
1023                          e_bio(i,j,p_bpin-1)   = e_bio(i,j,p_bpin-1)  + gas_emis*convert2
1024                       ELSE IF ( p_in_chem == p_limon ) THEN
1025                          ebio_lim(i,j) = ebio_lim(i,j) + gas_emis
1026                          e_bio(i,j,p_limon-1)   = e_bio(i,j,p_limon-1)  + gas_emis*convert2
1027                       ELSE IF ( p_in_chem == p_mbo ) THEN
1028                          ebio_mbo(i,j) = ebio_mbo(i,j) + gas_emis
1029                          e_bio(i,j,p_mbo-1)   = e_bio(i,j,p_mbo-1)  + gas_emis*convert2
1030                       ELSE IF ( p_in_chem == p_myrc ) THEN
1031                          ebio_myrc(i,j) = ebio_myrc(i,j) + gas_emis
1032                          e_bio(i,j,p_myrc-1)   = e_bio(i,j,p_myrc-1)  + gas_emis*convert2
1033                       ELSE IF ( p_in_chem == p_bcary ) THEN
1034                          ebio_sesq(i,j) = ebio_sesq(i,j) + gas_emis
1035                          e_bio(i,j,p_bcary-1)   = e_bio(i,j,p_bcary-1)  + gas_emis*convert2
1036                       ELSE IF ( p_in_chem == p_so2 ) THEN
1037                          ebio_so2(i,j) = ebio_so2(i,j) + gas_emis
1038                          e_bio(i,j,p_so2-1)   = e_bio(i,j,p_so2-1)  + gas_emis*convert2
1039                       ELSE IF ( p_in_chem == p_dms ) THEN
1040                          ebio_dms(i,j) = ebio_dms(i,j) + gas_emis
1041                          e_bio(i,j,p_dms-1)   = e_bio(i,j,p_dms-1)  + gas_emis*convert2
1042                       ELSE IF ( p_in_chem == p_bigalk ) THEN
1043                          ebio_bigalk(i,j) = ebio_bigalk(i,j) + gas_emis
1044                          e_bio(i,j,p_bigalk-1)   = e_bio(i,j,p_bigalk-1)  + gas_emis*convert2
1045                       ELSE IF ( p_in_chem == p_bigene ) THEN
1046                          ebio_bigene(i,j) = ebio_bigene(i,j) + gas_emis
1047                          e_bio(i,j,p_bigene-1)   = e_bio(i,j,p_bigene-1)  + gas_emis*convert2
1048                       ELSE IF ( p_in_chem == p_nh3 ) THEN
1049                          ebio_nh3(i,j) = ebio_nh3(i,j) + gas_emis
1050                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
1051                       ELSE IF ( p_in_chem == p_ch3oh ) THEN
1052                          ebio_ch3oh(i,j) = ebio_ch3oh(i,j) + gas_emis
1053                          e_bio(i,j,p_ch3oh-1)   = e_bio(i,j,p_ch3oh-1)  + gas_emis*convert2
1054                       ELSE IF ( p_in_chem == p_c2h5oh ) THEN
1055                          ebio_c2h5oh(i,j) = ebio_c2h5oh(i,j) + gas_emis
1056                          e_bio(i,j,p_c2h5oh-1)   = e_bio(i,j,p_c2h5oh-1)  + gas_emis*convert2
1057                       ELSE IF ( p_in_chem == p_ch3cooh ) THEN
1058                          ebio_ch3cooh(i,j) = ebio_ch3cooh(i,j) + gas_emis
1059                          e_bio(i,j,p_ch3cooh-1)   = e_bio(i,j,p_ch3cooh-1)  + gas_emis*convert2
1060                       ELSE IF ( p_in_chem == p_mek ) THEN
1061                          ebio_mek(i,j) = ebio_mek(i,j) + gas_emis
1062                          e_bio(i,j,p_mek-1)   = e_bio(i,j,p_mek-1)  + gas_emis*convert2
1063                       ELSE IF ( p_in_chem == p_c2h4 ) THEN
1064                          ebio_c2h4(i,j) = ebio_c2h4(i,j) + gas_emis
1065                          e_bio(i,j,p_c2h4-1)   = e_bio(i,j,p_c2h4-1)  + gas_emis*convert2
1066                       ELSE IF ( p_in_chem == p_c2h6 ) THEN
1067                          ebio_c2h6(i,j) = ebio_c2h6(i,j) + gas_emis
1068                          e_bio(i,j,p_c2h6-1)   = e_bio(i,j,p_c2h6-1)  + gas_emis*convert2
1069                       ELSE IF ( p_in_chem == p_c3h6 ) THEN
1070                          ebio_c3h6(i,j) = ebio_c3h6(i,j) + gas_emis
1071                          e_bio(i,j,p_c3h6-1)   = e_bio(i,j,p_c3h6-1)  + gas_emis*convert2
1072                       ELSE IF ( p_in_chem == p_c3h8 ) THEN
1073                          ebio_c3h8(i,j) = ebio_c3h8(i,j) + gas_emis
1074                          e_bio(i,j,p_c3h8-1)   = e_bio(i,j,p_c3h8-1)  + gas_emis*convert2
1075                       END IF
1076                    END IF is_mozm_species
1077                 END IF use_megan_emis
1078              END DO
1080           CASE (RADM2, RADM2_KPP, RADM2SORG, RADM2SORG_AQ, RADM2SORG_AQCHEM, RADM2SORG_KPP,GOCARTRADM2)
1082              DO icount = 1, n_megan2radm2
1084                 IF ( p_of_radm2(icount) .NE. non_react ) THEN
1085                 
1086                    ! Get index to chem array for the corresponding RADM2
1087                    ! species.  
1088                    p_in_chem = p_of_radm2(icount)
1090                    ! Check if the species is actually in the mechanism
1091                    IF ( p_in_chem >= param_first_scalar ) THEN
1092                       
1093                       ! Emission rate for mechanism species in mol km-2 hr-1
1094                       gas_emis = radm2_per_megan(icount) * E_megan2(p_of_megan2radm2(icount))
1096                       ! Add emissions to diagnostic output variables.
1097                       ! ebio_xxx (mol km-2 hr-1) were originally used by the 
1098                       ! BEIS3.11 biogenic emissions module. 
1099                       ! I have also borrowed variable e_bio (ppm m min-1).
1100                       IF ( p_in_chem .EQ. p_iso ) THEN
1101                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1102                          e_bio(i,j,p_iso-1)   = e_bio(i,j,p_iso-1)  + gas_emis*convert2
1103                       ELSE IF ( p_in_chem .EQ. p_oli) THEN
1104                          ebio_oli(i,j)        = ebio_oli(i,j)       + gas_emis
1105                          e_bio(i,j,p_oli-1)   = e_bio(i,j,p_oli-1)  + gas_emis*convert2
1106                       ELSE IF ( p_in_chem .EQ. p_hc3) THEN
1107                          ebio_hc3(i,j)        = ebio_hc3(i,j)       + gas_emis
1108                          e_bio(i,j,p_hc3-1)   = e_bio(i,j,p_hc3-1)  + gas_emis*convert2
1109                       ELSE IF ( p_in_chem .EQ. p_olt) THEN
1110                          ebio_olt(i,j)        = ebio_olt(i,j)       + gas_emis
1111                          e_bio(i,j,p_olt-1)   = e_bio(i,j,p_olt-1)  + gas_emis*convert2
1112                       ELSE IF ( p_in_chem .EQ. p_ket) THEN
1113                          ebio_ket(i,j)        = ebio_ket(i,j)       + gas_emis
1114                          e_bio(i,j,p_ket-1)   = e_bio(i,j,p_ket-1)  + gas_emis*convert2
1115                       ELSE IF ( p_in_chem .EQ. p_ald) THEN
1116                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1117                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
1118                       ELSE IF ( p_in_chem .EQ. p_hcho) THEN
1119                          ebio_hcho(i,j)       = ebio_hcho(i,j)      + gas_emis
1120                          e_bio(i,j,p_hcho-1)  = e_bio(i,j,p_hcho-1) + gas_emis*convert2
1121                       ELSE IF ( p_in_chem .EQ. p_eth) THEN
1122                          ebio_eth(i,j)        = ebio_eth(i,j)       + gas_emis
1123                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1124                       ELSE IF ( p_in_chem .EQ. p_ora2) THEN
1125                          ebio_ora2(i,j)       = ebio_ora2(i,j)      + gas_emis
1126                          e_bio(i,j,p_ora2-1)  = e_bio(i,j,p_ora2-1) + gas_emis*convert2
1127                       ELSE IF ( p_in_chem .EQ. p_co) THEN
1128                          ebio_co(i,j)         = ebio_co(i,j)        + gas_emis
1129                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)   + gas_emis*convert2
1130                       ELSE IF ( p_in_chem .EQ. p_no) THEN
1131                          ebio_no(i,j)         = ebio_no(i,j)        + gas_emis   
1132                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)   + gas_emis*convert2
1133                       ELSE IF ( p_in_chem .EQ. p_ol2) THEN
1134                           e_bio(i,j,p_ol2-1)  = e_bio(i,j,p_ol2-1)  + gas_emis*convert2
1135                       ELSE IF ( p_in_chem .EQ. p_hc5) THEN
1136                           e_bio(i,j,p_hc5-1)  = e_bio(i,j,p_hc5-1)  + gas_emis*convert2
1137                       ELSE IF ( p_in_chem .EQ. p_hc8) THEN
1138                           e_bio(i,j,p_hc8-1)  = e_bio(i,j,p_hc8-1)  + gas_emis*convert2
1139                       ELSE IF ( p_in_chem .EQ. p_ora1) THEN
1140                           e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
1141                       END IF
1143                    END IF !( p_in_chem >= param_first_scalar )
1145                 END IF !( p_of_ramd2(icount) .NE. non_react )
1146                 
1147              END DO
1149            CASE (RACMSORG_AQ, RACMSORG_AQCHEM_KPP, RACM_ESRLSORG_AQCHEM_KPP, RACM_ESRLSORG_KPP, RACM_KPP, GOCARTRACM_KPP, &
1150                  RACMSORG_KPP, RACM_MIM_KPP, RACMPM_KPP)
1152              DO icount = 1, n_megan2racm
1154                 IF ( p_of_racm(icount) .NE. non_react ) THEN
1156                    ! Get index to chem array for the corresponding RACM
1157                    ! species.  
1158                    p_in_chem = p_of_racm(icount)
1159                    
1160                    ! Check if the species is actually in the mechanism
1161                    IF( p_in_chem >= param_first_scalar ) THEN
1163                       ! Emission rate of mechanism species in mol km-2 hr-1
1164                       gas_emis =  racm_per_megan(icount) * E_megan2(p_of_megan2racm(icount))
1166                       ! Add emissions to diagnostic output variables.
1167                       ! ebio_xxx (mol km-2 hr-1) were originally used by the 
1168                       ! BEIS3.11 biogenic emissions module. 
1169                       ! I have also borrowed variable e_bio (ppm m min-1).
1170                       IF ( p_in_chem .EQ. p_iso ) THEN
1171                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1172                          e_bio(i,j,p_iso-1)   = e_bio(i,j,p_iso-1)  + gas_emis*convert2
1173                       ELSE IF ( p_in_chem .EQ. p_oli) THEN
1174                          ebio_oli(i,j)        = ebio_oli(i,j)       + gas_emis
1175                          e_bio(i,j,p_oli-1)   = e_bio(i,j,p_oli-1)  + gas_emis*convert2
1176                       ELSE IF ( p_in_chem .EQ. p_api) THEN
1177                          ebio_api(i,j)        = ebio_api(i,j)       + gas_emis
1178                          e_bio(i,j,p_api-1)   = e_bio(i,j,p_api-1)  + gas_emis*convert2
1179                       ELSE IF ( p_in_chem .EQ. p_lim) THEN
1180                          ebio_lim(i,j)        = ebio_lim(i,j)       + gas_emis
1181                          e_bio(i,j,p_lim-1)   = e_bio(i,j,p_lim-1)  + gas_emis*convert2
1182                       ELSE IF ( p_in_chem .EQ. p_hc3) THEN
1183                          ebio_hc3(i,j)        = ebio_hc3(i,j)       + gas_emis
1184                          e_bio(i,j,p_hc3-1)   = e_bio(i,j,p_hc3-1)  + gas_emis*convert2
1185                       ELSE IF ( p_in_chem .EQ. p_ete) THEN
1186                          ebio_ete(i,j)        = ebio_ete(i,j)       + gas_emis
1187                          e_bio(i,j,p_ete-1)   = e_bio(i,j,p_ete-1)  + gas_emis*convert2
1188                       ELSE IF ( p_in_chem .EQ. p_olt) THEN
1189                          ebio_olt(i,j)        = ebio_olt(i,j)       + gas_emis
1190                          e_bio(i,j,p_olt-1)   = e_bio(i,j,p_olt-1)  + gas_emis*convert2
1191                       ELSE IF ( p_in_chem .EQ. p_ket) THEN
1192                          ebio_ket(i,j)        = ebio_ket(i,j)       + gas_emis
1193                          e_bio(i,j,p_ket-1)   = e_bio(i,j,p_ket-1)  + gas_emis*convert2
1194                       ELSE IF ( p_in_chem .EQ. p_ald) THEN
1195                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1196                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
1197                       ELSE IF ( p_in_chem .EQ. p_hcho) THEN
1198                          ebio_hcho(i,j)       = ebio_hcho(i,j)      + gas_emis
1199                          e_bio(i,j,p_hcho-1)  = e_bio(i,j,p_hcho-1) + gas_emis*convert2
1200                       ELSE IF ( p_in_chem .EQ. p_eth) THEN
1201                          ebio_eth(i,j)        = ebio_eth(i,j)       + gas_emis
1202                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1203                       ELSE IF ( p_in_chem .EQ. p_ora2) THEN
1204                          ebio_ora2(i,j)       = ebio_ora2(i,j)      + gas_emis
1205                          e_bio(i,j,p_ora2-1)  = e_bio(i,j,p_ora2-1) + gas_emis*convert2
1206                       ELSE IF ( p_in_chem .EQ. p_co) THEN
1207                          ebio_co(i,j)         = ebio_co(i,j)        + gas_emis
1208                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)   + gas_emis*convert2
1209                       ELSE IF ( p_in_chem .EQ. p_no) THEN
1210                          ebio_no(i,j)         = ebio_no(i,j)        + gas_emis   
1211                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)   + gas_emis*convert2
1212                       ELSE IF ( p_in_chem .EQ. p_hc5) THEN
1213                           e_bio(i,j,p_hc5-1)  = e_bio(i,j,p_hc5-1)  + gas_emis*convert2
1214                       ELSE IF ( p_in_chem .EQ. p_hc8) THEN
1215                           e_bio(i,j,p_hc8-1)  = e_bio(i,j,p_hc8-1)  + gas_emis*convert2
1216                       ELSE IF ( p_in_chem .EQ. p_ora1) THEN
1217                           e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
1218                       END IF
1220                    END IF !( p_in_chem > param_first_scalar )
1221                    
1223                 END IF !( p_of_racm(icount) .NE. non_react )
1225              END DO
1227           CASE (CB05_SORG_AQ_KPP)
1229              DO icount = 1, n_megan2cb05
1230                 IF ( p_of_cb05 (icount) .NE. non_react ) THEN
1231                    ! Get index to chem array for the corresponding CB05
1232                    ! species.
1233                    p_in_chem = p_of_cb05(icount)
1235                    ! Check if the species is actually in the mechanism
1236                    ! (e.g., NH3 is in the mechanism only if aerosols
1237                    ! are simulated)
1238                    ! Check if the species is actually in the mechanism
1239                    IF ( p_in_chem >= param_first_scalar ) THEN
1241                       ! Emission rate for mechanism species in mol km-2 hr-1
1242                       gas_emis = cb05_per_megan(icount) * E_megan2(p_of_megan2cb05(icount))
1244                       ! Add emissions to diagnostic output variables.
1245                       ! ebio_xxx (mol km-2 hr-1) were originally used by the
1246                       ! BEIS3.11 biogenic emissions module.
1247                       ! I have also borrowed variable e_bio (ppm m min-1).
1248                       IF ( p_in_chem .EQ. p_isop ) THEN
1249                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1250                          e_bio(i,j,p_isop-1)   = e_bio(i,j,p_isop-1)  + gas_emis*convert2
1251                       END IF
1252                       IF ( p_in_chem .EQ. p_aacd ) THEN
1253                          e_bio(i,j,p_aacd-1)  = e_bio(i,j,p_aacd-1)  + gas_emis*convert2
1254                       END IF
1255                       IF ( p_in_chem .EQ. p_ald2 ) THEN
1256                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1257                          e_bio(i,j,p_ald2-1)  = e_bio(i,j,p_ald2-1)  + gas_emis*convert2
1258                       END IF
1259                       IF ( p_in_chem .EQ. p_aldx ) THEN
1260                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1261                          e_bio(i,j,p_aldx-1)  = e_bio(i,j,p_aldx-1)  + gas_emis*convert2
1262                       END IF
1263                       IF ( p_in_chem .EQ. p_apin ) THEN
1264                          ebio_api(i,j)        = ebio_api(i,j)       + gas_emis
1265                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1266                       END IF
1267                       IF ( p_in_chem .EQ. p_bpin ) THEN
1268                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1269                       END IF
1270                       IF ( p_in_chem .EQ. p_ch4 ) THEN
1271                          e_bio(i,j,p_ch4-1)   = e_bio(i,j,p_ch4-1)  + gas_emis*convert2
1272                       END IF
1273                       IF ( p_in_chem .EQ. p_co ) THEN
1274                          ebio_co(i,j)        = ebio_co(i,j)       + gas_emis
1275                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)  + gas_emis*convert2
1276                       END IF
1277                       IF ( p_in_chem .EQ. p_eth ) THEN
1278                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1279                       END IF
1280                       IF ( p_in_chem .EQ. p_etha ) THEN
1281                          e_bio(i,j,p_etha-1)  = e_bio(i,j,p_etha-1)  + gas_emis*convert2
1282                       END IF
1283                       IF ( p_in_chem .EQ. p_etoh ) THEN
1284                          e_bio(i,j,p_etoh-1)  = e_bio(i,j,p_etoh-1)  + gas_emis*convert2
1285                       END IF
1286                       IF ( p_in_chem .EQ. p_facd ) THEN
1287                          e_bio(i,j,p_facd-1)  = e_bio(i,j,p_facd-1)  + gas_emis*convert2
1288                       END IF
1289                       IF ( p_in_chem .EQ. p_form ) THEN
1290                          ebio_hcho(i,j)        = ebio_hcho(i,j)       + gas_emis
1291                          e_bio(i,j,p_form-1)  = e_bio(i,j,p_form-1)  + gas_emis*convert2
1292                       END IF
1293                       IF ( p_in_chem .EQ. p_hum ) THEN
1294                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1295                       END IF
1296                       IF ( p_in_chem .EQ. p_iole ) THEN
1297                          e_bio(i,j,p_iole-1)  = e_bio(i,j,p_iole-1)  + gas_emis*convert2
1298                       END IF
1299                       IF ( p_in_chem .EQ. p_lim ) THEN
1300                          ebio_lim(i,j)        = ebio_lim(i,j)       + gas_emis
1301                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1302                       END IF
1303                       IF ( p_in_chem .EQ. p_meoh ) THEN
1304                          e_bio(i,j,p_meoh-1)  = e_bio(i,j,p_meoh-1)  + gas_emis*convert2
1305                       END IF
1306                       IF ( p_in_chem .EQ. p_nh3 ) THEN
1307                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
1308                       END IF
1309                       IF ( p_in_chem .EQ. p_no ) THEN
1310                          ebio_no(i,j)        = ebio_no(i,j)       + gas_emis
1311                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)  + gas_emis*convert2
1312                       END IF
1313                       IF ( p_in_chem .EQ. p_oci ) THEN
1314                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1) + gas_emis*convert2
1315                       END IF
1316                       IF ( p_in_chem .EQ. p_ole ) THEN
1317                          e_bio(i,j,p_ole-1)   = e_bio(i,j,p_ole-1)  + gas_emis*convert2
1318                       END IF
1319                       IF ( p_in_chem .EQ. p_par ) THEN
1320                          e_bio(i,j,p_par-1)   = e_bio(i,j,p_par-1)  + gas_emis*convert2
1321                       END IF
1322                       IF ( p_in_chem .EQ. p_terp ) THEN
1323                          ebio_terp(i,j)        = ebio_terp(i,j)       + gas_emis
1324                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1325                       END IF
1326                       IF ( p_in_chem .EQ. p_ter ) THEN
1327                          ebio_terp(i,j)        = ebio_terp(i,j)       + gas_emis
1328                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1329                       END IF
1330                       IF ( p_in_chem .EQ. p_tol ) THEN
1331                          e_bio(i,j,p_tol-1)   = e_bio(i,j,p_tol-1)  + gas_emis*convert2
1332                       END IF
1334                    END IF !( p_in_chem > param_first_scalar )
1336                 END IF
1337              END DO
1339           CASE (CB05_SORG_VBS_AQ_KPP)
1341              DO icount = 1, n_megan2cb05vbs
1342                 IF ( p_of_cb05vbs (icount) .NE. non_react ) THEN
1343                    ! Get index to chem array for the corresponding CB05
1344                    ! species.
1345                    p_in_chem = p_of_cb05vbs(icount)
1347                    ! Check if the species is actually in the mechanism
1348                    ! (e.g., NH3 is in the mechanism only if aerosols
1349                    ! are simulated)
1350                    ! Check if the species is actually in the mechanism
1351                    IF ( p_in_chem >= param_first_scalar ) THEN
1353                       ! Emission rate for mechanism species in mol km-2 hr-1
1354                       gas_emis = cb05vbs_per_megan(icount) * E_megan2(p_of_megan2cb05vbs(icount))
1356                       ! Add emissions to diagnostic output variables.
1357                       ! ebio_xxx (mol km-2 hr-1) were originally used by the
1358                       ! BEIS3.11 biogenic emissions module.
1359                       ! I have also borrowed variable e_bio (ppm m min-1).
1360                       IF ( p_in_chem .EQ. p_isop ) THEN
1361                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1362                          e_bio(i,j,p_isop-1)   = e_bio(i,j,p_isop-1)  + gas_emis*convert2
1363                       END IF
1364                       IF ( p_in_chem .EQ. p_aacd ) THEN
1365                          e_bio(i,j,p_aacd-1)  = e_bio(i,j,p_aacd-1)  + gas_emis*convert2
1366                       END IF
1367                       IF ( p_in_chem .EQ. p_ald2 ) THEN
1368                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1369                          e_bio(i,j,p_ald2-1)  = e_bio(i,j,p_ald2-1)  + gas_emis*convert2
1370                       END IF
1371                       IF ( p_in_chem .EQ. p_aldx ) THEN
1372                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1373                          e_bio(i,j,p_aldx-1)  = e_bio(i,j,p_aldx-1)  + gas_emis*convert2
1374                       END IF
1375                       IF ( p_in_chem .EQ. p_apin ) THEN
1376                          ebio_api(i,j)        = ebio_api(i,j)       + gas_emis
1377                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1378                       END IF
1379                       IF ( p_in_chem .EQ. p_bpin ) THEN
1380                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1381                       END IF
1382                       IF ( p_in_chem .EQ. p_ch4 ) THEN
1383                          e_bio(i,j,p_ch4-1)   = e_bio(i,j,p_ch4-1)  + gas_emis*convert2
1384                       END IF
1385                       IF ( p_in_chem .EQ. p_co ) THEN
1386                          ebio_co(i,j)        = ebio_co(i,j)       + gas_emis
1387                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)  + gas_emis*convert2
1388                       END IF
1389                       IF ( p_in_chem .EQ. p_eth ) THEN
1390                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1391                       END IF
1392                       IF ( p_in_chem .EQ. p_etha ) THEN
1393                          e_bio(i,j,p_etha-1)  = e_bio(i,j,p_etha-1)  + gas_emis*convert2
1394                       END IF
1395                       IF ( p_in_chem .EQ. p_etoh ) THEN
1396                          e_bio(i,j,p_etoh-1)  = e_bio(i,j,p_etoh-1)  + gas_emis*convert2
1397                       END IF
1398                       IF ( p_in_chem .EQ. p_facd ) THEN
1399                          e_bio(i,j,p_facd-1)  = e_bio(i,j,p_facd-1)  + gas_emis*convert2
1400                       END IF
1401                       IF ( p_in_chem .EQ. p_form ) THEN
1402                          ebio_hcho(i,j)        = ebio_hcho(i,j)       + gas_emis
1403                          e_bio(i,j,p_form-1)  = e_bio(i,j,p_form-1)  + gas_emis*convert2
1404                       END IF
1405                       IF ( p_in_chem .EQ. p_hum ) THEN
1406                           e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1407                       END IF
1408                       IF ( p_in_chem .EQ. p_iole ) THEN
1409                          e_bio(i,j,p_iole-1)  = e_bio(i,j,p_iole-1)  + gas_emis*convert2
1410                       END IF
1411                       IF ( p_in_chem .EQ. p_lim ) THEN
1412                          ebio_lim(i,j)        = ebio_lim(i,j)       + gas_emis
1413                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1414                       END IF
1415                       IF ( p_in_chem .EQ. p_meoh ) THEN
1416                          e_bio(i,j,p_meoh-1)  = e_bio(i,j,p_meoh-1)  + gas_emis*convert2
1417                       END IF
1418                       IF ( p_in_chem .EQ. p_nh3 ) THEN
1419                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
1420                       END IF
1421                       IF ( p_in_chem .EQ. p_no ) THEN
1422                          ebio_no(i,j)        = ebio_no(i,j)       + gas_emis
1423                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)  + gas_emis*convert2
1424                       END IF
1425                       IF ( p_in_chem .EQ. p_oci ) THEN
1426                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1427                       END IF
1428                       IF ( p_in_chem .EQ. p_ole ) THEN
1429                          e_bio(i,j,p_ole-1)   = e_bio(i,j,p_ole-1)  + gas_emis*convert2
1430                       END IF
1431                       IF ( p_in_chem .EQ. p_par ) THEN
1432                          e_bio(i,j,p_par-1)   = e_bio(i,j,p_par-1)  + gas_emis*convert2
1433                       END IF
1434                       IF ( p_in_chem .EQ. p_terp ) THEN
1435                          ebio_terp(i,j)        = ebio_terp(i,j)       + gas_emis
1436                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1437                       END IF
1438                       IF ( p_in_chem .EQ. p_ter ) THEN
1439                          ebio_terp(i,j)        = ebio_terp(i,j)       + gas_emis
1440                          e_bio(i,j,p_terp-1)   = e_bio(i,j,p_terp-1)  + gas_emis*convert2
1441                       END IF
1442                       IF ( p_in_chem .EQ. p_tol ) THEN
1443                          e_bio(i,j,p_tol-1)   = e_bio(i,j,p_tol-1)  + gas_emis*convert2
1444                       END IF
1446                    END IF !( p_in_chem > param_first_scalar )
1448                 END IF
1449              END DO
1451           CASE (RACM_SOA_VBS_KPP,RACM_SOA_VBS_AQCHEM_KPP,RACM_SOA_VBS_HET_KPP)
1453           DO icount = 1, n_megan2racmSOA
1455                 IF ( p_of_racmSOA(icount) .NE. non_react ) THEN
1457                    ! Get index to chem array for the corresponding RACM-SOA-VBS-KPP
1458                    ! species.  
1459                    p_in_chem = p_of_racmSOA(icount)
1461                    ! Check if the species is actually in the mechanism
1462                    IF( p_in_chem >= param_first_scalar ) THEN
1464                       ! Emission rate of mechanism species in mol km-2 hr-1
1465                       gas_emis =  racmSOA_per_megan(icount) * E_megan2(p_of_megan2racmSOA(icount))
1467                       ! Add emissions to diagnostic output variables.
1468                       ! ebio_xxx (mol km-2 hr-1) were originally used by the 
1469                       ! BEIS3.11 biogenic emissions module. 
1470                       ! I have also borrowed variable e_bio (ppm m min-1).
1471                       IF ( p_in_chem .EQ. p_iso ) THEN
1472                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1473                          e_bio(i,j,p_iso-1)   = e_bio(i,j,p_iso-1)  + gas_emis*convert2
1474                       ELSE IF ( p_in_chem .EQ. p_oli) THEN
1475                          ebio_oli(i,j)        = ebio_oli(i,j)       + gas_emis
1476                          e_bio(i,j,p_oli-1)   = e_bio(i,j,p_oli-1)  + gas_emis*convert2
1477                       ELSE IF ( p_in_chem .EQ. p_api) THEN
1478                          ebio_api(i,j)        = ebio_api(i,j)       + gas_emis
1479                          e_bio(i,j,p_api-1)   = e_bio(i,j,p_api-1)  + gas_emis*convert2
1480                       ELSE IF ( p_in_chem .EQ. p_lim) THEN
1481                          ebio_lim(i,j)        = ebio_lim(i,j)       + gas_emis
1482                          e_bio(i,j,p_lim-1)   = e_bio(i,j,p_lim-1)  + gas_emis*convert2
1483                       ELSE IF ( p_in_chem .EQ. p_hc3) THEN
1484                          ebio_hc3(i,j)        = ebio_hc3(i,j)       + gas_emis
1485                          e_bio(i,j,p_hc3-1)   = e_bio(i,j,p_hc3-1)  + gas_emis*convert2
1486                       ELSE IF ( p_in_chem .EQ. p_ete) THEN
1487                          ebio_ete(i,j)        = ebio_ete(i,j)       + gas_emis
1488                          e_bio(i,j,p_ete-1)   = e_bio(i,j,p_ete-1)  + gas_emis*convert2
1489                       ELSE IF ( p_in_chem .EQ. p_olt) THEN
1490                          ebio_olt(i,j)        = ebio_olt(i,j)       + gas_emis
1491                          e_bio(i,j,p_olt-1)   = e_bio(i,j,p_olt-1)  + gas_emis*convert2
1492                       ELSE IF ( p_in_chem .EQ. p_ket) THEN
1493                          ebio_ket(i,j)        = ebio_ket(i,j)       + gas_emis
1494                          e_bio(i,j,p_ket-1)   = e_bio(i,j,p_ket-1)  + gas_emis*convert2
1495                       ELSE IF ( p_in_chem .EQ. p_ald) THEN
1496                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1497                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
1498                       ELSE IF ( p_in_chem .EQ. p_hcho) THEN
1499                          ebio_hcho(i,j)       = ebio_hcho(i,j)      + gas_emis
1500                          e_bio(i,j,p_hcho-1)  = e_bio(i,j,p_hcho-1) + gas_emis*convert2
1501                       ELSE IF ( p_in_chem .EQ. p_eth) THEN
1502                          ebio_eth(i,j)        = ebio_eth(i,j)       + gas_emis
1503                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1504                       ELSE IF ( p_in_chem .EQ. p_ora2) THEN
1505                          ebio_ora2(i,j)       = ebio_ora2(i,j)      + gas_emis
1506                          e_bio(i,j,p_ora2-1)  = e_bio(i,j,p_ora2-1) + gas_emis*convert2
1507                       ELSE IF ( p_in_chem .EQ. p_co) THEN
1508                          ebio_co(i,j)         = ebio_co(i,j)        + gas_emis
1509                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)   + gas_emis*convert2
1510                       ELSE IF ( p_in_chem .EQ. p_no) THEN
1511                          ebio_no(i,j)         = ebio_no(i,j)        + gas_emis
1512                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)   + gas_emis*convert2
1513                       ELSE IF ( p_in_chem .EQ. p_hc5) THEN
1514                           e_bio(i,j,p_hc5-1)  = e_bio(i,j,p_hc5-1)  + gas_emis*convert2
1515                       ELSE IF ( p_in_chem .EQ. p_hc8) THEN
1516                           e_bio(i,j,p_hc8-1)  = e_bio(i,j,p_hc8-1)  + gas_emis*convert2
1517                       ELSE IF ( p_in_chem .EQ. p_ora1) THEN
1518                           e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
1519                       ELSE IF ( p_in_chem .EQ. p_sesq) THEN
1520                           ebio_sesq(i,j)      = ebio_sesq(i,j)      + gas_emis
1521                           e_bio(i,j,p_sesq-1) = e_bio(i,j,p_sesq-1) + gas_emis*convert2
1522                       ELSE IF ( p_in_chem .EQ. p_mbo) THEN
1523                           ebio_mbo(i,j)       = ebio_mbo(i,j)        + gas_emis
1524                           e_bio(i,j,p_mbo-1)  = e_bio(i,j,p_mbo-1)   + gas_emis*convert2
1525                       END IF
1527                    END IF !( p_in_chem > param_first_scalar )
1530                 END IF !( p_of_racm(icount) .NE. non_react )
1532              END DO
1533           CASE (CBMZ, CBMZ_BB, CBMZ_BB_KPP, CBMZ_MOSAIC_KPP, &
1534                 CBMZ_MOSAIC_4BIN, &
1535                 CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ, &
1536                 CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, &
1537                 CBMZ_MOSAIC_DMS_4BIN_AQ,CBMZ_MOSAIC_DMS_8BIN_AQ,CBMZSORG, CBMZSORG_AQ, &
1538                 CBMZ_CAM_MAM3_NOAQ, CBMZ_CAM_MAM3_AQ, CBMZ_CAM_MAM7_NOAQ, CBMZ_CAM_MAM7_AQ)
1540              DO icount = 1, n_megan2cbmz
1542                 IF ( p_of_cbmz (icount) .NE. non_react ) THEN
1544                    ! Get index to chem array for the corresponding CBMZ
1545                    ! species.  
1546                    p_in_chem = p_of_cbmz(icount)
1548                    ! Check if the species is actually in the mechanism
1549                    ! (e.g., NH3 is in the mechanism only if aerosols
1550                    ! are simulated)
1551                    IF( p_in_chem >= param_first_scalar ) THEN
1553                       ! Emission rate of mechanism species in mol km-2 hr-1
1554                       gas_emis = cbmz_per_megan(icount) * E_megan2(p_of_megan2cbmz(icount))
1557                       ! Add emissions to diagnostic output variables.
1558                       ! ebio_xxx (mol km-2 hr-1) were originally used by the 
1559                       ! BEIS3.11 biogenic emissions module. 
1560                       ! I have also borrowed variable e_bio (ppm m min-1).
1561                       IF ( p_in_chem .EQ. p_iso ) THEN
1562                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1563                          e_bio(i,j,p_iso-1)   = e_bio(i,j,p_iso-1)  + gas_emis*convert2
1564                       ELSE IF ( p_in_chem .EQ. p_oli) THEN
1565                          ebio_oli(i,j)        = ebio_oli(i,j)       + gas_emis
1566                          e_bio(i,j,p_oli-1)   = e_bio(i,j,p_oli-1)  + gas_emis*convert2
1567                       ELSE IF ( p_in_chem .EQ. p_olt) THEN
1568                          ebio_olt(i,j)        = ebio_olt(i,j)       + gas_emis
1569                          e_bio(i,j,p_olt-1)   = e_bio(i,j,p_olt-1)  + gas_emis*convert2
1570                       ELSE IF ( p_in_chem .EQ. p_ket) THEN
1571                          ebio_ket(i,j)        = ebio_ket(i,j)       + gas_emis
1572                          e_bio(i,j,p_ket-1)   = e_bio(i,j,p_ket-1)  + gas_emis*convert2
1573                       ELSE IF ( p_in_chem .EQ. p_ald) THEN
1574                          ebio_ald(i,j)        = ebio_ald(i,j)       + gas_emis
1575                          e_bio(i,j,p_ald-1)   = e_bio(i,j,p_ald-1)  + gas_emis*convert2
1576                       ELSE IF ( p_in_chem .EQ. p_hcho) THEN
1577                          ebio_hcho(i,j)       = ebio_hcho(i,j)      + gas_emis
1578                          e_bio(i,j,p_hcho-1)  = e_bio(i,j,p_hcho-1) + gas_emis*convert2
1579                       ELSE IF ( p_in_chem .EQ. p_eth) THEN
1580                          ebio_eth(i,j)        = ebio_eth(i,j)       + gas_emis
1581                          e_bio(i,j,p_eth-1)   = e_bio(i,j,p_eth-1)  + gas_emis*convert2
1582                       ELSE IF ( p_in_chem .EQ. p_ora2) THEN
1583                          ebio_ora2(i,j)       = ebio_ora2(i,j)      + gas_emis
1584                          e_bio(i,j,p_ora2-1)  = e_bio(i,j,p_ora2-1) + gas_emis*convert2
1585                       ELSE IF ( p_in_chem .EQ. p_co) THEN
1586                          ebio_co(i,j)         = ebio_co(i,j)        + gas_emis
1587                          e_bio(i,j,p_co-1)    = e_bio(i,j,p_co-1)   + gas_emis*convert2
1588                       ELSE IF ( p_in_chem .EQ. p_no) THEN
1589                          ebio_no(i,j)         = ebio_no(i,j)        + gas_emis   
1590                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)   + gas_emis*convert2
1591                       ELSE IF ( p_in_chem .EQ. p_ol2) THEN
1592                           e_bio(i,j,p_ol2-1)  = e_bio(i,j,p_ol2-1)  + gas_emis*convert2
1593                       ELSE IF ( p_in_chem .EQ. p_ora1) THEN
1594                           e_bio(i,j,p_ora1-1) = e_bio(i,j,p_ora1-1) + gas_emis*convert2
1595                       
1596                       ! SAN, 08/11/13 - adding missing CBMZ species to be mapped: 
1597                       ! missing: p_par, p_ch3oh, p_c2h5oh, p_nh3, p_tol
1598                       ELSE IF ( p_in_chem .EQ. p_par) THEN 
1599                          !ebio_par(i,j)        = ebio_par(i,j)       + gas_emis
1600                          e_bio(i,j,p_par-1)   = e_bio(i,j,p_par-1)  + gas_emis*convert2
1601                       ELSE IF ( p_in_chem .EQ. p_ch3oh) THEN    
1602                          ebio_ch3oh(i,j)      = ebio_ch3oh(i,j)     + gas_emis
1603                          e_bio(i,j,p_ch3oh-1) = e_bio(i,j,p_ch3oh-1)+ gas_emis*convert2
1604                       ELSE IF ( p_in_chem .EQ. p_c2h5oh) THEN   
1605                          ebio_c2h5oh(i,j)     = ebio_c2h5oh(i,j)      + gas_emis
1606                          e_bio(i,j,p_c2h5oh-1)= e_bio(i,j,p_c2h5oh-1) + gas_emis*convert2
1607                       ELSE IF ( p_in_chem .EQ. p_nh3) THEN      
1608                          ebio_nh3(i,j)        = ebio_nh3(i,j)       + gas_emis
1609                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
1610                       ELSE IF ( p_in_chem .EQ. p_tol) THEN      
1611                          ebio_tol(i,j)       = ebio_tol(i,j)        + gas_emis
1612                          e_bio(i,j,p_tol-1)  = e_bio(i,j,p_tol-1)   + gas_emis*convert2
1614                       END IF
1617                    END IF !( p_in_chem > param_first_scalar )
1618                    
1619                    
1620                 END IF ! ( p_of_cbmz (icount) .NE. non_react )
1622              END DO
1623             
1624           CASE (SAPRC99_KPP,SAPRC99_MOSAIC_4BIN_VBS2_KPP, &
1625                SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP,SAPRC99_MOSAIC_8BIN_VBS2_KPP)!BSINGH(12/03/2013): Added SAPRC 8 bin and non-aq on (04/07/2014) ! FIX FOR SAPRC99 AND SAPRC07
1627              DO icount = 1, n_megan2saprcnov
1629                 IF ( p_of_saprcnov(icount) .NE. non_react ) THEN
1631                    ! Get index to chem array for the corresponding RADM2
1632                    ! species.
1633                    p_in_chem = p_of_saprcnov(icount)
1635                    ! Check if the species is actually in the mechanism
1636                    IF ( p_in_chem >= param_first_scalar ) THEN
1638                       ! Emission rate for mechanism species in mol km-2 hr-1
1639                       gas_emis = saprcnov_per_megan(icount) * E_megan2(p_of_megan2saprcnov(icount))
1641                       ! Add emissions to diagnostic output variables.
1642                       ! ebio_xxx (mol km-2 hr-1) were originally used by the
1643                       ! BEIS3.11 biogenic emissions module.
1644                       ! I have also borrowed variable e_bio (ppm m min-1).
1645                       IF ( p_in_chem .EQ. p_isoprene ) THEN
1646                          ebio_iso(i,j)        = ebio_iso(i,j)       + gas_emis
1647                          e_bio(i,j,p_isoprene-1)   = e_bio(i,j,p_isoprene-1)  + gas_emis*convert2
1648                       ELSE IF ( p_in_chem .EQ. p_terp) THEN
1649                          ebio_api(i,j)       = ebio_api(i,j)      + gas_emis
1650                          e_bio(i,j,p_terp-1)  = e_bio(i,j,p_terp-1) + gas_emis*convert2
1651                       ELSE IF ( p_in_chem .EQ. p_sesq) THEN
1652                          ebio_lim(i,j)         = ebio_lim(i,j)        + gas_emis
1653                          e_bio(i,j,p_sesq-1)    = e_bio(i,j,p_sesq-1)   + gas_emis*convert2
1654                       ELSE IF ( p_in_chem .EQ. p_no) THEN
1655                          ebio_no(i,j)         = ebio_no(i,j)        + gas_emis
1656                          e_bio(i,j,p_no-1)    = e_bio(i,j,p_no-1)   + gas_emis*convert2
1657 !jdf
1658                       ELSE IF ( p_in_chem .EQ. p_alk3) THEN
1659                          ebio_alk3(i,j)         = ebio_alk3(i,j)        + gas_emis
1660                          e_bio(i,j,p_alk3-1)    = e_bio(i,j,p_alk3-1)   + gas_emis*convert2
1661                       ELSE IF ( p_in_chem .EQ. p_alk4) THEN
1662                          ebio_alk4(i,j)         = ebio_alk4(i,j)        + gas_emis
1663                          e_bio(i,j,p_alk4-1)    = e_bio(i,j,p_alk4-1)   + gas_emis*convert2
1664                       ELSE IF ( p_in_chem .EQ. p_alk5) THEN
1665                          ebio_alk5(i,j)         = ebio_alk5(i,j)        + gas_emis
1666                          e_bio(i,j,p_alk5-1)    = e_bio(i,j,p_alk5-1)   + gas_emis*convert2
1667                       ELSE IF ( p_in_chem .EQ. p_ole1) THEN
1668                          ebio_ole1(i,j)         = ebio_ole1(i,j)        + gas_emis
1669                          e_bio(i,j,p_ole1-1)    = e_bio(i,j,p_ole1-1)   + gas_emis*convert2
1670                       ELSE IF ( p_in_chem .EQ. p_ole2) THEN
1671                          ebio_ole2(i,j)         = ebio_ole2(i,j)        + gas_emis
1672                          e_bio(i,j,p_ole2-1)    = e_bio(i,j,p_ole2-1)   + gas_emis*convert2
1673                       ELSE IF ( p_in_chem .EQ. p_aro1) THEN
1674                          ebio_aro1(i,j)         = ebio_aro1(i,j)        + gas_emis
1675                          e_bio(i,j,p_aro1-1)    = e_bio(i,j,p_aro1-1)   + gas_emis*convert2
1676                       ELSE IF ( p_in_chem .EQ. p_aro2) THEN
1677                          ebio_aro2(i,j)         = ebio_aro2(i,j)        + gas_emis
1678                          e_bio(i,j,p_aro2-1)    = e_bio(i,j,p_aro2-1)   + gas_emis*convert2
1679                       ELSE IF ( p_in_chem .EQ. p_acet) THEN
1680                          ebio_acet(i,j)         = ebio_acet(i,j)        + gas_emis
1681                          e_bio(i,j,p_acet-1)    = e_bio(i,j,p_acet-1)   + gas_emis*convert2
1682                       ELSE IF ( p_in_chem .EQ. p_hcho) THEN
1683                          ebio_hcho(i,j)         = ebio_hcho(i,j)        + gas_emis
1684                          e_bio(i,j,p_hcho-1)    = e_bio(i,j,p_hcho-1)   + gas_emis*convert2
1685                       ELSE IF ( p_in_chem .EQ. p_ccho) THEN
1686                          ebio_ccho(i,j)         = ebio_ccho(i,j)        + gas_emis
1687                          e_bio(i,j,p_ccho-1)    = e_bio(i,j,p_ccho-1)   + gas_emis*convert2
1688                       ELSE IF ( p_in_chem .EQ. p_mek) THEN
1689                          ebio_mek(i,j)         = ebio_mek(i,j)        + gas_emis
1690                          e_bio(i,j,p_mek-1)    = e_bio(i,j,p_mek-1)   + gas_emis*convert2
1691                       ELSE IF ( p_in_chem .EQ. p_c2h6) THEN
1692                          ebio_c2h6(i,j)         = ebio_c2h6(i,j)        + gas_emis
1693                          e_bio(i,j,p_c2h6-1)    = e_bio(i,j,p_c2h6-1)   + gas_emis*convert2
1694                       ELSE IF ( p_in_chem .EQ. p_c3h6) THEN
1695                          ebio_c3h6(i,j)         = ebio_c3h6(i,j)        + gas_emis
1696                          e_bio(i,j,p_c3h6-1)    = e_bio(i,j,p_c3h6-1)   + gas_emis*convert2
1697                       ELSE IF ( p_in_chem .EQ. p_c3h8) THEN
1698                          ebio_c3h8(i,j)         = ebio_c3h8(i,j)        + gas_emis
1699                          e_bio(i,j,p_c3h8-1)    = e_bio(i,j,p_c3h8-1)   + gas_emis*convert2
1700                       ELSE IF ( p_in_chem .EQ. p_ethene) THEN
1701                          ebio_ethene(i,j)         = ebio_ethene(i,j)        + gas_emis
1702                          e_bio(i,j,p_ethene-1)    = e_bio(i,j,p_ethene-1)   + gas_emis*convert2
1703                       ELSE IF ( p_in_chem .EQ. p_bald) THEN
1704                          ebio_bald(i,j)         = ebio_bald(i,j)        + gas_emis
1705                          e_bio(i,j,p_bald-1)    = e_bio(i,j,p_bald-1)   + gas_emis*convert2
1706                       ELSE IF ( p_in_chem .EQ. p_meoh) THEN
1707                          ebio_meoh(i,j)         = ebio_meoh(i,j)        + gas_emis
1708                          e_bio(i,j,p_meoh-1)    = e_bio(i,j,p_meoh-1)   + gas_emis*convert2
1709                       ELSE IF ( p_in_chem .EQ. p_hcooh) THEN
1710                          ebio_hcooh(i,j)         = ebio_hcooh(i,j)        + gas_emis
1711                          e_bio(i,j,p_hcooh-1)    = e_bio(i,j,p_hcooh-1)   + gas_emis*convert2
1712                       ELSE IF ( p_in_chem .EQ. p_rco_oh) THEN
1713                          ebio_rco_oh(i,j)         = ebio_rco_oh(i,j)        + gas_emis
1714                          e_bio(i,j,p_rco_oh-1)    = e_bio(i,j,p_rco_oh-1)   + gas_emis*convert2
1715                       ELSE IF ( p_in_chem .EQ. p_terp) THEN
1716                          ebio_terp(i,j)         = ebio_terp(i,j)        + gas_emis
1717                          ebio_api(i,j)         = ebio_api(i,j)        + gas_emis
1718                          e_bio(i,j,p_terp-1)    = e_bio(i,j,p_terp-1)   + gas_emis*convert2
1719                       ELSE IF ( p_in_chem .EQ. p_sesq) THEN
1720                          ebio_sesq(i,j)         = ebio_sesq(i,j)        + gas_emis
1721                          ebio_lim(i,j)         = ebio_lim(i,j)        + gas_emis
1722                          e_bio(i,j,p_sesq-1)    = e_bio(i,j,p_sesq-1)   + gas_emis*convert2
1723 !jdf
1724                       END IF
1726                    END IF !( p_in_chem > param_first_scalar )
1728                 END IF !( p_of_saprcnov(icount) .NE. non_react )
1730              END DO
1732           CASE ( CRIMECH_KPP, CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP )
1734              DO icount = 1, n_megan2crimech
1735                 IF ( p_of_crimech(icount) .NE. non_react ) THEN
1737                    ! Get index to chem array for the corresponding crimech
1738                    ! species.  
1739                    p_in_chem = p_of_crimech(icount)
1740                    
1741                    ! Check if the species is actually in the mechanism
1742                    IF( p_in_chem >= param_first_scalar ) THEN
1744                       ! Emission rate of mechanism species in mol km-2 hr-1
1745                       gas_emis =  crimech_per_megan(icount) * E_megan2(p_of_megan2crimech(icount))
1747                       ! Add emissions to diagnostic output variables.
1748                       ! ebio_xxx (mol km-2 hr-1) were originally used by the 
1749                       ! BEIS3.11 biogenic emissions module. 
1750                       ! I have also borrowed variable e_bio (ppm m min-1).
1751                       
1752                       IF ( p_in_chem == p_c5h8 ) THEN
1753                          ebio_c5h8(i,j) = ebio_c5h8(i,j) + gas_emis
1754                          e_bio(i,j,p_c5h8-1)   = e_bio(i,j,p_c5h8-1)  + gas_emis*convert2
1755                       ELSE IF ( p_in_chem == p_no ) THEN
1756                          ebio_no(i,j)  = ebio_no(i,j) + gas_emis
1757                          e_bio(i,j,p_no-1)   = e_bio(i,j,p_no-1)  + gas_emis*convert2
1758                       ELSE IF ( p_in_chem == p_no2 ) THEN
1759                          ebio_no2(i,j)  = ebio_no2(i,j) + gas_emis
1760                          e_bio(i,j,p_no2-1)   = e_bio(i,j,p_no2-1)  + gas_emis*convert2
1761                       ELSE IF ( p_in_chem == p_co ) THEN
1762                          ebio_co(i,j)  = ebio_co(i,j) + gas_emis
1763                          e_bio(i,j,p_co-1)   = e_bio(i,j,p_co-1)  + gas_emis*convert2
1764                       ELSE IF ( p_in_chem == p_hcho ) THEN
1765                          ebio_hcho(i,j) = ebio_hcho(i,j) + gas_emis
1766                          e_bio(i,j,p_hcho-1)   = e_bio(i,j,p_hcho-1)  + gas_emis*convert2
1767                       ELSE IF ( p_in_chem == p_ket ) THEN
1768                          ebio_ket(i,j) = ebio_ket(i,j) + gas_emis
1769                          e_bio(i,j,p_ket-1)   = e_bio(i,j,p_ket-1)  + gas_emis*convert2
1770                       ELSE IF ( p_in_chem == p_toluene ) THEN
1771                          ebio_toluene(i,j) = ebio_toluene(i,j) + gas_emis
1772                          e_bio(i,j,p_toluene-1)   = e_bio(i,j,p_toluene-1)  + gas_emis*convert2
1773                       ELSE IF ( p_in_chem == p_apinene ) THEN
1774                          ebio_apinene(i,j) = ebio_apinene(i,j) + gas_emis
1775                          e_bio(i,j,p_apinene-1)   = e_bio(i,j,p_apinene-1)  + gas_emis*convert2
1776                       ELSE IF ( p_in_chem == p_bpinene ) THEN
1777                          ebio_bpinene(i,j) = ebio_bpinene(i,j) + gas_emis
1778                          e_bio(i,j,p_bpinene-1)   = e_bio(i,j,p_bpinene-1)  + gas_emis*convert2                         
1779                       ELSE IF ( p_in_chem == p_so2 ) THEN
1780                          ebio_so2(i,j) = ebio_so2(i,j) + gas_emis
1781                          e_bio(i,j,p_so2-1)   = e_bio(i,j,p_so2-1)  + gas_emis*convert2
1782                       ELSE IF ( p_in_chem == p_dms ) THEN
1783                          ebio_dms(i,j) = ebio_dms(i,j) + gas_emis
1784                          e_bio(i,j,p_dms-1)   = e_bio(i,j,p_dms-1)  + gas_emis*convert2
1785                       ELSE IF ( p_in_chem == p_nc4h10 ) THEN
1786                          ebio_nc4h10(i,j) = ebio_nc4h10(i,j) + gas_emis
1787                          e_bio(i,j,p_nc4h10-1)   = e_bio(i,j,p_nc4h10-1)  + gas_emis*convert2
1788                       ELSE IF ( p_in_chem == p_tbut2ene ) THEN
1789                          ebio_tbut2ene(i,j) = ebio_tbut2ene(i,j) + gas_emis
1790                          e_bio(i,j,p_tbut2ene-1)   = e_bio(i,j,p_tbut2ene-1)  + gas_emis*convert2
1791                       ELSE IF ( p_in_chem == p_nh3 ) THEN
1792                          ebio_nh3(i,j) = ebio_nh3(i,j) + gas_emis
1793                          e_bio(i,j,p_nh3-1)   = e_bio(i,j,p_nh3-1)  + gas_emis*convert2
1794                       ELSE IF ( p_in_chem == p_ch3oh ) THEN
1795                          ebio_ch3oh(i,j) = ebio_ch3oh(i,j) + gas_emis
1796                          e_bio(i,j,p_ch3oh-1)   = e_bio(i,j,p_ch3oh-1)  + gas_emis*convert2
1797                       ELSE IF ( p_in_chem == p_c2h5oh ) THEN
1798                          ebio_c2h5oh(i,j) = ebio_c2h5oh(i,j) + gas_emis
1799                          e_bio(i,j,p_c2h5oh-1)   = e_bio(i,j,p_c2h5oh-1)  + gas_emis*convert2
1800                       ELSE IF ( p_in_chem == p_ch3co2h ) THEN
1801                          ebio_ch3co2h(i,j) = ebio_ch3co2h(i,j) + gas_emis
1802                          e_bio(i,j,p_ch3co2h-1)   = e_bio(i,j,p_ch3co2h-1)  + gas_emis*convert2
1803                       ELSE IF ( p_in_chem == p_mek ) THEN
1804                          ebio_mek(i,j) = ebio_mek(i,j) + gas_emis
1805                          e_bio(i,j,p_mek-1)   = e_bio(i,j,p_mek-1)  + gas_emis*convert2
1806                       ELSE IF ( p_in_chem == p_c2h4 ) THEN
1807                          ebio_c2h4(i,j) = ebio_c2h4(i,j) + gas_emis
1808                          e_bio(i,j,p_c2h4-1)   = e_bio(i,j,p_c2h4-1)  + gas_emis*convert2
1809                       ELSE IF ( p_in_chem == p_c2h6 ) THEN
1810                          ebio_c2h6(i,j) = ebio_c2h6(i,j) + gas_emis
1811                          e_bio(i,j,p_c2h6-1)   = e_bio(i,j,p_c2h6-1)  + gas_emis*convert2
1812                       ELSE IF ( p_in_chem == p_c3h6 ) THEN
1813                          ebio_c3h6(i,j) = ebio_c3h6(i,j) + gas_emis
1814                          e_bio(i,j,p_c3h6-1)   = e_bio(i,j,p_c3h6-1)  + gas_emis*convert2
1815                       ELSE IF ( p_in_chem == p_c3h8 ) THEN
1816                          ebio_c3h8(i,j) = ebio_c3h8(i,j) + gas_emis
1817                          e_bio(i,j,p_c3h8-1)   = e_bio(i,j,p_c3h8-1)  + gas_emis*convert2                         
1818                       ELSE IF ( p_in_chem == p_ch3cho ) THEN
1819                          ebio_ch3cho(i,j) = ebio_ch3cho(i,j) + gas_emis
1820                          e_bio(i,j,p_ch3cho-1)   = e_bio(i,j,p_ch3cho-1)  + gas_emis*convert2                        
1821                       ELSE IF ( p_in_chem == p_hcooh ) THEN
1822                          ebio_hcooh(i,j) = ebio_hcooh(i,j) + gas_emis
1823                          e_bio(i,j,p_hcooh-1)   = e_bio(i,j,p_hcooh-1)  + gas_emis*convert2                         
1824                       END IF
1826                    END IF !( p_in_chem > param_first_scalar )
1827                    
1829                 END IF !( p_of_crimech(icount) .NE. non_react )
1831              END DO
1835              CASE DEFAULT
1837                 CALL wrf_error_fatal('Species conversion table for MEGAN v2.04 not available. ')
1839              END SELECT GAS_MECH_SELECT
1843        END DO i_loop ! i = its, ite
1844     END DO j_loop    ! j = jts, jte
1847   CONTAINS
1849     ! -----------------------------------------------------------------
1850     !  SUBROUTINE GAMMA_TISOP returns the GAMMA_T value for isoprene
1851     !  Orginally from Tan's gamma_etc.F
1852     ! -----------------------------------------------------------------
1854     SUBROUTINE GAMMA_TISOP( TEMP, D_TEMP, gam_T )
1855       !
1856       !  Description :
1857       !
1858       !    MEGAN biogenic emissions adjustment factor for temperature
1859       !    for isoprene
1860       !
1861       !  Reference: 
1862       !
1863       !    Estimates of global terrestial isoprene emissions using MEGAN
1864       !    (Model of Emissions of Gases and Aerosols from Nature )
1865       !    A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, 
1866       !    P.I. Palmer, and C. Geron
1867       !    Atmospheric Chemistry and Physics, 6, 3181-3210, 2006      !
1868       !
1870       IMPLICIT NONE
1872       ! hourly surface air temperature (K)
1873       ! (here use instantaneous temperature
1874       REAL, INTENT(IN)  :: TEMP
1875       ! daily-mean surface airtemperature (K)
1876       ! (here use the previous month's monthly mean)
1877       REAL, INTENT(IN)  :: D_TEMP
1878       !temperature adjustment factor
1879       REAL, INTENT(OUT) :: gam_T
1881       ! Local parameters
1882       REAL :: Eopt, Topt, X
1883       REAL :: AAA, BBB
1884       REAL, PARAMETER :: CT1 = 80.0
1885       REAL, PARAMETER :: CT2 = 200.0
1886       
1887       ! End header ----------------------------------------------------
1889       ! Below Eqn (14) of Guenther et al. [2006]
1890       ! (assuming T_daily = D_TEMP)
1891       Eopt = 1.75 * EXP(0.08*(D_TEMP-297.0))
1893       ! Eqn (8) of Guenther et al. [2006]
1894       ! (assuming T_daily = D_TEMP)
1895       Topt = 313.0 + ( 0.6*(D_TEMP-297.0) )
1897       ! Eqn (5) of Guenther et al. [2006]
1898       X = ( (1.0/Topt)-(1.0/TEMP) ) / 0.00831
1899       AAA = Eopt*CT2*EXP(CT1*X)
1900       BBB = (  CT2-CT1*( 1.0-EXP(CT2*X) )  )
1901       gam_T = AAA/BBB
1903     END SUBROUTINE GAMMA_TISOP
1905     ! -----------------------------------------------------------------
1906     ! SUBROUITINE GAMMA_TNISP returns the GAMMA_T value for 
1907     ! non-isoprene species
1908     ! Originally from Tan's gamma_etc.F
1909     !------------------------------------------------------------------
1911     SUBROUTINE GAMMA_TNISP( SPCNUM, TEMP, gam_T )
1912       !
1913       !  Description :
1914       !
1915       !    MEGAN biogenic emissions adjustment factor for temperature
1916       !    for non-isoprene species.
1917       !
1918       !  Reference:
1919       !
1920       !    MEGAN v2.0 Documentation
1921       !
1922       ! Method:
1923       !
1924       !    GAMMA_T =  exp[BETA*(T-Ts)]
1925       !      where BETA   = temperature dependent parameter
1926       !            Ts     = standard temperature (normally 303K, 30C)
1927       !
1929       IMPLICIT NONE
1931       INTEGER, INTENT(IN) :: SPCNUM               ! Species number
1932       REAL, INTENT(IN)    :: TEMP
1933       REAL, INTENT(OUT)   :: gam_T
1934       REAL, PARAMETER     :: Ts = 303.0
1935       
1936       ! End header ----------------------------------------------------
1938       ! TDF_PRM is defined in module_data_megan2.F
1939       gam_T = EXP( TDF_PRM(SPCNUM)*(TEMP-Ts) )
1941     END SUBROUTINE GAMMA_TNISP
1944     ! --------------------------------------------------------------------
1945     ! SUBROUTINE GAMMA_LAI
1946     ! Originally from Tan's gamma_etc.F
1947     ! --------------------------------------------------------------------
1949     SUBROUTINE GAMMA_LAI(LAI, gam_L )
1950       !  Description :
1951       !
1952       !    MEGAN biogenic emissions adjustment factor for leaf area
1953       !    index
1954       !
1955       !  Reference: 
1956       !
1957       !    Estimates of global terrestial isoprene emissions using MEGAN
1958       !    (Model of Emissions of Gases and Aerosols from Nature )
1959       !    A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, 
1960       !    P.I. Palmer, and C. Geron
1961       !    Atmospheric Chemistry and Physics, 6, 3181-3210, 2006      !
1962       !
1963       ! Method:
1964       !                       0.49[LAI]
1965       !        GAMMA_LAI = ----------------    (dimensionless)
1966       !                    (1+0.2LAI^2)^0.5
1967       !
1969       IMPLICIT NONE
1970       REAL, INTENT(IN)  ::  LAI 
1971       REAL, INTENT(OUT) :: gam_L
1973       ! End header ----------------------------------------------------
1975       
1976       ! Eqn (15) of Guenther et al. [2006]
1977       gam_L = (0.49*LAI) / ( SQRT(1.0+0.2*(LAI**2)) )
1979       RETURN
1980     END SUBROUTINE GAMMA_LAI
1982     !-------------------------------------------------------------------
1983     ! SUBROUTINE GAMMA_P 
1984     ! Originally from Tan's gamma_etc.F
1985     !-------------------------------------------------------------------
1986     SUBROUTINE GAMMA_P(             &
1987          DOY_in, tmidh, LAT, LONG,  &                    
1988          PPFD, D_PPFD, gam_P        &
1989          )
1990       !
1991       !  Description :
1992       !
1993       !    MEGAN biogenic emissions adjustment factor for
1994       !    photosynthetic photon flux density (PPFD or PAR)
1995       !
1996       !  Reference: 
1997       !
1998       !    Estimates of global terrestial isoprene emissions using MEGAN
1999       !    (Model of Emissions of Gases and Aerosols from Nature )
2000       !    A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, 
2001       !    P.I. Palmer, and C. Geron
2002       !    Atmospheric Chemistry and Physics, 6, 3181-3210, 2006      
2003       !  
2004       !  Method:
2005       !
2006       !    GAMMA_P = 0.0         sin(a)<=0
2007       !
2008       !    GAMMA_P = sin(a)[2.46*0.9*PHI^3*(1+0.0005(Pdaily-400))]
2009       !                                  0<a<180
2010       !           where PHI = above canopy PPFD transmission (non-dimension)
2011       !           Pdaily    = daily average above canopy PPFD (umol/m2s)
2012       !              a      = solar angle (degree)
2013       !
2014       !         Note: AAA = 2.46*BBB*PHI-0.9*PHI^2
2015       !               BBB = (1+0.0005(Pdaily-400))
2016       !           GAMMA_P = sin(a)*AAA
2017       !
2018       !                       Pac
2019       !             PHI = -----------
2020       !                   sin(a)*Ptoa
2021       !
2022       !     where Pac  = above canopy PPFD (umol/m2s)
2023       !                 Ptoa = PPFD at the top of atmosphere (umol/m2s)
2024       !
2025       !             Pac =  SRAD * 4.766 mmmol/m2-s * 0.5
2026       !
2027       !             Ptoa = 3000 + 99*cos[2*3.14-( DOY-10)/365 )]
2028       !        where DOY = day of year (julian day)
2029       !
2030       ! NOTE: This code has been corrected. The gamma P equation as defined in the
2031       ! original Guenther et al., 2006 (equation 11b) is incorrect. This has the
2032       ! corrected algorithm. (CW, 08/16/2007)
2033       !-----------------------------------------------------------------
2035       IMPLICIT NONE
2037       INTEGER, INTENT(IN) :: DOY_in ! julian day at GMT
2039       ! GMT hour plus minutes (in fractaionl hour) of the middle
2040       ! of the current time step
2041       REAL, INTENT(IN)  :: tmidh
2042       REAL, INTENT(IN)  ::  LAT    ! Latitude [=] degrees
2043       REAL, INTENT(IN)  ::  LONG   ! Longitude [=] degrees
2044       REAL, INTENT(IN)  ::  PPFD   ! Photosynthetic Photon Flus Density
2045       REAL, INTENT(IN)  ::  D_PPFD ! Daily PPFD
2046       REAL, INTENT(OUT) ::  gam_P  ! GAMMA_P
2049       !...  Local scalars
2050       INTEGER :: DOY                 ! local julian day
2051       REAL :: HOUR                   ! solar hour
2052       REAL :: AAA, BBB
2053       REAL :: SIN_solarangle         ! sin(solar angle)
2054       REAL :: Ptoa, Pac, PHI
2056       ! End header ----------------------------------------------------
2058       ! Convert time of the middle of the current time step
2059       ! from GMT to solar hour (include minutes in decimals)
2060       DOY = DOY_in
2061       HOUR = tmidh + long/15.
2062       IF ( HOUR .LT. 0.0 ) THEN
2063          HOUR = HOUR + 24.0
2064          DOY  = DOY - 1
2065       ENDIF
2067       ! Above canopy photosynthetic photo flux density (PPFD)
2068       ! ( micromole/m2/s )
2069       Pac = PPFD
2071       ! Get sin of solar elevation angle
2072       CALL SOLARANGLE( DOY, HOUR, LAT, SIN_solarangle )
2074       ! Calculate gamma_p in Eqn (10) of Guenther et al. [2006]
2075       IF ( SIN_solarangle .LE. 0.0 ) THEN
2076          ! Eqn (11a) of Guenther et al. [2006]
2077          gam_P = 0.0
2078       ELSE
2079          ! PPFD at top of the atmosphere
2080          ! Eqn (13) of Guenther et al. [2006]
2081          ! ( micromole/m2/s )
2082          Ptoa = 3000.0 + 99.0 * COS( 2.*3.14*(DOY-10.)/365. )
2083          ! Above canopy PPFD transmission
2084          ! Eqn (12) of Guenther et al. [2006]
2085          ! (nondimensional)
2086          PHI = Pac/(SIN_solarangle * Ptoa)
2087          ! Eqn (11b) of Guenther et al. [2006]
2088          ! (Note: typo in the paper; correction made 08/06/2007)
2089          BBB = 1. + 0.0005*( D_PPFD-400. )
2090          AAA = 2.46 * BBB * PHI - 0.9 * (PHI**2)
2091          gam_P = SIN_solarangle * AAA
2093       ENDIF
2094       ! Screening the unforced errors
2095       ! IF solar elevation angle is less than 1 THEN
2096       ! gamma_p can not be greater than 0.1.
2097       IF (SIN_solarangle .LE. 0.0175 .AND. gam_P .GT. 0.1) THEN
2098          gam_P = 0.1
2099       ENDIF
2102     END SUBROUTINE GAMMA_P
2104     ! ----------------------------------------------------------------
2105     ! SUBROUTINE GAMMA_A returns GAMMA_A
2106     ! Originally from Tan's gamma_etc.F
2107     !------------------------------------------------------------------
2108     SUBROUTINE GAMMA_A( i_spc, LAIp, LAIc, TSTLEN, D_TEMP, gam_A )
2109       !  Description :
2110       !
2111       !    MEGAN biogenic emissions adjustment factor for leaf age
2112       !
2113       !  Reference: 
2114       !
2115       !    Estimates of global terrestial isoprene emissions using MEGAN
2116       !    (Model of Emissions of Gases and Aerosols from Nature )
2117       !    A. Guenther, T. Karl, P. Harley, C. Wiedinmyer, 
2118       !    P.I. Palmer, and C. Geron
2119       !    Atmospheric Chemistry and Physics, 6, 3181-3210, 2006
2120       !
2121       !    MEGAN v2.0 Documentation
2122       !
2123       !
2124       ! Method:
2125       !
2126       !     GAMMA_age = Fnew*Anew + Fgro*Agro + Fmat*Amat + Fold*Aold
2127       !      where Fnew = new foliage fraction
2128       !            Fgro = growing foliage fraction
2129       !                 Fmat = mature foliage fraction
2130       !                 Fold = old foliage fraction
2131       !                 Anew = relative emission activity for new foliage
2132       !                 Agro = relative emission activity for growing foliage
2133       !                 Amat = relative emission activity for mature foliage
2134       !                 Aold = relative emission activity for old foliage
2135       !
2136       !
2137       !             For foliage fraction
2138       !             Case 1) LAIc = LAIp
2139       !             Fnew = 0.0  , Fgro = 0.1  , Fmat = 0.8  , Fold = 0.1
2140       !
2141       !             Case 2) LAIp > LAIc
2142       !             Fnew = 0.0  , Fgro = 0.0
2143       !             Fmat = 1-Fold
2144       !             Fold = (LAIp-LAIc)/LAIp
2145       !
2146       !             Case 3) LAIp < LAIc
2147       !             Fnew = 1-(LAIp/LAIc)                       t <= ti
2148       !                  = (ti/t) * ( 1-(LAIp/LAIc) )          t >  ti
2149       !
2150       !             Fmat = LAIp/LAIc                           t <= tm
2151       !                  = (LAIp/LAIc) +
2152       !                      ( (t-tm)/t ) * ( 1-(LAIp/LAIc) )  t >  tm
2153       !
2154       !             Fgro = 1 - Fnew - Fmat
2155       !             Fold = 0.0
2156       !
2157       !           where
2158       !             ti = 5 + (0.7*(300-Tt))                   Tt <= 303
2159       !                = 2.9                                  Tt >  303
2160       !             tm = 2.3*ti
2161       !
2162       !             t  = length of the time step (days)
2163       !             ti = number of days between budbreak and the induction of
2164       !                  emission
2165       !             tm = number of days between budbreak and the initiation of
2166       !                  peak emissions rates
2167       !             Tt = average temperature (K) near top of the canopy during
2168       !                  current time period (daily ave temp for this case)
2169       !
2170       !
2171       !             For relative emission activity
2172       !             Case 1) Constant
2173       !             Anew = 1.0  , Agro = 1.0  , Amat = 1.0  , Aold = 1.0
2174       !
2175       !             Case 2) Monoterpenes
2176       !             Anew = 2.0  , Agro = 1.8  , Amat = 0.95 , Aold = 1.0
2177       !
2178       !             Case 3) Sesquiterpenes
2179       !             Anew = 0.4  , Agro = 0.6  , Amat = 1.075, Aold = 1.0
2180       !
2181       !             Case 4) Methanol
2182       !             Anew = 3.0  , Agro = 2.6  , Amat = 0.85 , Aold = 1.0
2183       !
2184       !             Case 5) Isoprene
2185       !             Anew = 0.05 , Agro = 0.6  , Amat = 1.125, Aold = 1.0
2188       IMPLICIT NONE
2190       ! SUBROUTINE arguments
2192       !..."Pointer" for class of species
2193       INTEGER, INTENT(IN) :: i_spc
2194       !...average temperature of the previous 24-hours
2195       REAL, INTENT(IN) :: D_TEMP
2196       !...leaf area index of the current and previous
2197       !...month
2198       REAL, INTENT(IN) :: LAIp, LAIc
2199       !...time step between LAIc and LAIp (days)
2200       REAL, INTENT(IN) ::     TSTLEN
2201       !...emissions adjustment factor accounting for leaf age
2202       REAL, INTENT(OUT) :: gam_A
2204       ! Local scalars
2206       !...leaf age fractions
2207       REAL :: Fnew, Fgro, Fmat, Fold
2208       !...relative emission activity index
2209       INTEGER ::  AINDX 
2210       !...time step between LAIC and LAIp (days)
2211       INTEGER :: t 
2212       !...number of days between budbreak and the induction emission
2213       REAL     ti
2214       !...number of days between budbreak  and the initiation of peak
2215       !...emissions rates
2216       REAL     tm
2217       !
2219       REAL     Tt                   ! average temperature (K)
2220       ! daily ave temp
2222       ! End header ----------------------------------------------------
2224       ! Choose relative emission activity class
2225       ! See Table 2 of MEGAN v2.0 Documentation
2226       !
2228       IF (      (i_spc==imgn_acto) .OR. (i_spc==imgn_acta) .OR. (i_spc==imgn_form)   &
2229            .OR. (i_spc==imgn_ch4)  .OR. (i_spc==imgn_no)   .OR. (i_spc==imgn_co)     &
2230            ) THEN
2231          AINDX = 1
2233       ELSE IF ( (i_spc==imgn_myrc) .OR. (i_spc==imgn_sabi) .OR. (i_spc==imgn_limo)   &
2234            .OR. (i_spc==imgn_3car) .OR. (i_spc==imgn_ocim) .OR. (i_spc==imgn_bpin)   &
2235            .OR. (i_spc==imgn_apin) .OR. ( i_spc==imgn_omtp)                          &
2236            ) THEN
2237          AINDX = 2
2239       ELSE IF ( (i_spc==imgn_afarn) .OR. (i_spc==imgn_bcar) .OR. (i_spc==imgn_osqt)  &
2240            ) THEN
2241          AINDX = 3
2243       ELSE IF (i_spc==imgn_meoh) THEN
2244          aindx = 4
2246       ELSE IF ( (i_spc==imgn_isop) .OR. (i_spc==imgn_mbo) ) THEN
2247          aindx = 5
2248       ELSE
2249          WRITE(mesg,fmt = '("Invalid i_spc in SUBROUTINE GAMMA_A; i_spc = ", I3)') i_spc
2250          CALL wrf_error_fatal(mesg)
2251       END IF
2255       ! Time step between LAIp and LAIc (days)
2256       t = TSTLEN
2257       ! Tt is the average ambient air temperature (K) of the preceeding time
2258       ! interval.  Here, use the monthly-mean surface air temperature
2259       Tt   = D_TEMP
2261       ! Calculate foliage fraction
2262       ! (section 3.2.2 of Guenther et al. [2006])
2263       IF (LAIp .EQ. LAIc) THEN
2264          Fnew = 0.0
2265          Fgro = 0.1
2266          Fmat = 0.8
2267          Fold = 0.1
2268       ELSEIF (LAIp .GT. LAIc) THEN
2269          Fnew = 0.0
2270          Fgro = 0.0
2271          Fold = ( LAIp-LAIc ) / LAIp
2272          Fmat = 1.0-Fold
2273       ELSE ! LAIp < LAIc
2274          ! Calculate ti, which is the number of days between budbreak and
2275          ! the induction of isoprene emission.
2276          IF (Tt .LE. 303.0) THEN
2277             ! Eqn (18a) of Guenther et al. [2006]
2278             ti = 5.0 + 0.7*(300.0-Tt)
2279          ELSE
2280             ! Eqn (18b) of Guenther et al. [2006]
2281             ti = 2.9
2282          ENDIF
2283          ! tm is the number of days between budbreak and the initiation
2284          ! of peak isoprene emissions rates.
2285          ! Eqn (19) of Guenther et al. [2006]
2286          tm = 2.3*ti
2288          ! Calculate Fnew and Fmat, then Fgro and Fold
2289          !  Fnew
2290          IF (t .LE. ti) THEN
2291             ! Eqn (17a) of Guenther et al. [2006]
2292             Fnew = 1.0 - (LAIp/LAIc)
2293          ELSE
2294             ! Eqn (17b) of Guenther et al. [2006]
2295             Fnew = (ti/t) * ( 1-(LAIp/LAIc) )
2296          ENDIF
2298          ! Fmat
2299          IF (t .LE. tm) THEN
2300             ! Eqn (17c) of Guenther et al. [2006]
2301             Fmat = LAIp/LAIc
2302          ELSE
2303             ! Eqn (17d) of Guenther et al. [2006]
2304             Fmat = (LAIp/LAIc) + ( (t-tm)/t ) * ( 1-(LAIp/LAIc) )
2305          ENDIF
2307          Fgro = 1.0 - Fnew - Fmat
2308          Fold = 0.0
2310       ENDIF
2312       !Calculate GAMMA_A
2313       ! Anew, Agro, Amat, Aold are defined in module_data_megan2.F
2314       gam_A = Fnew*Anew(AINDX) + Fgro*Agro(AINDX)    &
2315            + Fmat*Amat(AINDX) + Fold*Aold(AINDX)
2318     END SUBROUTINE GAMMA_A
2320     ! ----------------------------------------------------------------
2321     ! SUBROUTINE SOLARANGLE calculates the solar angle
2322     ! Originally from Tan's solarangle.F
2323     !------------------------------------------------------------------
2324     SUBROUTINE SOLARANGLE( DAY, SHOUR, LAT, SIN_solarangle )
2325       !
2326       !
2327       !   Input:
2328       !            1) Day of year
2329       !            2) Latitude
2330       !            3) Hour
2331       !
2332       !   Output: sin of solar angle
2333       !
2335       IMPLICIT NONE
2337       ! Arguments
2338       INTEGER, INTENT(IN) :: DAY                  ! DOY or julian day
2339       REAL, INTENT(IN)    :: SHOUR                ! Solar hour
2340       REAL, INTENT(IN)    :: LAT                  ! Latitude
2341       REAL, INTENT(OUT)   :: SIN_solarangle
2343       ! Local scalars
2344       REAL    :: sindelta, cosdelta, A, B
2346       ! End header -----------------------------------------------------
2348       sindelta = -SIN(0.40907) * COS( 6.28*(REAL(DAY,KIND(0.))+10.)/365. )
2349       cosdelta = SQRT(1.-sindelta**2.)
2351       A = SIN( LAT*D2RAD ) * sindelta
2352       B = COS( LAT*D2RAD ) * cosdelta
2354       SIN_solarangle = A + B * COS(2.*PI*(SHOUR-12.)/24.)
2357     END SUBROUTINE SOLARANGLE
2362   END SUBROUTINE bio_emissions_megan2
2364 END MODULE module_bioemi_megan2