chem/module_cbmz.F

   1 !************************************************************************
   2 ! This computer software was prepared by Battelle Memorial Institute,
   3 ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
   4 ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
   5 ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
   6 ! LIABILITY FOR THE USE OF THIS SOFTWARE.
   7 !
   8 ! Chemistry Option:  CBMZ (Carbon Bond Mechanism IV - Zaveri)
   9 ! * Primary investigator: Rahul A. Zaveri
  10 ! * Co-investigator: Richard C. Easter, William I. Gustafson Jr.
  11 ! Last update: February 2009
  12 !
  13 ! Contacts:
  14 ! Rahul A. Zaveri, PhD                    Jerome D. Fast, PhD
  15 ! Senior Research Scientist               Staff Scientist
  16 ! Pacific Northwest National Laboratory   Pacific Northwest National Laboratory
  17 ! P.O. Box 999, MSIN K9-30                P.O. Box 999, MSIN K9-30
  18 ! Richland, WA 99352                      Richland, WA, 99352
  19 ! Phone: (509) 372-6159                   Phone: (509) 372-6116
  20 ! Email: Rahul.Zaveri@pnl.gov             Email: Jerome.Fast@pnl.gov
  21 !
  22 ! Please report any bugs or problems to Rahul Zaveri, the primary author
  23 ! of the code, or Jerome Fast, the WRF-chem implementation team leader.
  24 !
  25 ! Terms of Use:
  26 !  1) Users are requested to consult the primary author prior to
  27 !     modifying the CBMZ code or incorporating it or its submodules in
  28 !     another code. This is meant to ensure that the any linkages and/or
  29 !     assumptions will not adversely affect the operation of CBMZ.
  30 !  2) The CBMZ source code is intended for research and educational
  31 !     purposes. Users are requested to contact the primary author
  32 !     regarding the use of CBMZ code for any commercial application.
  33 !  3) Users preparing publications resulting from the usage of CBMZ are
  34 !     requested to cite one or more of the references below (depending on
  35 !     the application) for proper acknowledgement.
  36 !
  37 ! References:
  38 ! 1) Zaveri R.A., and L.K. Peters (1999), A new lumped structure
  39 !    photochemical mechanism for large-scale applications, J. Geophys.
  40 !    Res., 104, 30387-30415.
  41 ! 2) Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C.
  42 !    Barnard, E.G., Chapman, G.A. Grell, and S.E. Peckham (2005),
  43 !    Evolution of ozone, particulates, and aerosol direct radiative
  44 !    forcing in the vicinity of Houston using a fully-coupled
  45 !    meteorology- chemistry-aerosol model, J. Geophys. Res., 111, D21305,
  46 !    doi:10.1029/2005JD006721.
  47 !
  48 ! Additional information:
  49 ! * www.pnl.gov/atmospheric/research/wrf-chem
  50 !
  51 ! Support:
  52 ! Funding for developing and evaluating CBMZ was provided by the U.S.
  53 ! Department of Energy under the auspices of Atmospheric Science Program
  54 ! of the Office of Biological and Environmental Research the the PNNL
  55 ! Laboratory Research and Directed Research and Development program.
  56 !************************************************************************
  57       module module_cbmz
  58
  59
  60
  61       use module_peg_util
  62
  63       contains
  64
  65
  66 !***********************************************************************
  67 ! < 1.> subr cbmz_driver
  68 !
  69 ! purpose: serves as an interface between subr. gas_chemistry and
  70 !          the actual solver subr such as lsodes, rodas, etc.
  71 !
  72 ! grid   : fixed i,j,k  (box-model)
  73 !
  74 ! author : Rahul A. Zaveri
  75 ! date   : November 1998
  76 !
  77 !-----------------------------------------------------------------------
  78
  79       subroutine cbmz_driver( &
  80                id, curr_secs, ktau, dtstep, ktauc, dtstepc, &
  81                config_flags, gmt, julday, t_phy, moist, p8w, t8w, &
  82                p_phy, chem, rho_phy, dz8w, z, z_at_w, vdrog3, &
  83                ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
  84                ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
  85                ph_ch3o2h, ph_n2o5, &
  86                ids,ide, jds,jde, kds,kde, &
  87                ims,ime, jms,jme, kms,kme, &
  88                its,ite, jts,jte, kts,kte  )
  89
  90    USE module_configure, only:  grid_config_rec_type, num_moist, num_chem,  &
  91 !       p_qv, p_so2, p_ho2, p_so4aj, p_corn, p_hcl, p_mtf
  92 !jdf
  93         p_qv, p_so2, p_ho2, p_so4aj, p_corn, p_hcl, p_mtf, &
  94         p_xyl, p_tol, p_csl, p_oli, p_olt, p_ho, p_o3, p_no
  95 !jdf
  96    USE module_data_sorgam, only:  ldrog
  97    USE module_data_cbmz
  98    IMPLICIT NONE
  99
 100
 101 !-----------------------------------------------------------------------
 102 ! subr arguments
 103
 104    INTEGER, INTENT(IN ) :: id, julday, &
 105                            ids,ide, jds,jde, kds,kde, &
 106                            ims,ime, jms,jme, kms,kme, &
 107                            its,ite, jts,jte, kts,kte
 108
 109    INTEGER, INTENT(IN ) :: ktau, ktauc
 110
 111    REAL(KIND=8), INTENT(IN) :: curr_secs
 112
 113    REAL, INTENT(IN ) :: dtstep, dtstepc, gmt
 114 !
 115 ! advected moisture variables
 116 !
 117    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
 118          INTENT(IN ) :: moist
 119 !
 120 ! advected chemical tracers
 121 !
 122    REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
 123          INTENT(INOUT ) :: chem
 124 !
 125 ! arrays that hold photolysis rates
 126 !
 127    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
 128          INTENT(INOUT ) :: &
 129            ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
 130            ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
 131            ph_ch3o2h, ph_n2o5
 132 !
 133 ! on input from met model
 134 !
 135    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
 136           INTENT(IN ) :: &
 137                          t_phy, &       ! temperature
 138                          rho_phy, &     ! air density (kg/m3)
 139                          p_phy, &       ! NOT USED
 140                          z, z_at_w, &   ! NOT USED
 141                          dz8w, &        ! NOT USED
 142                          t8w, p8w       ! NOT USED
 143 !
 144 ! for interaction with aerosols (really is output)
 145 !
 146    REAL, DIMENSION( ims:ime , kms:kme-0 , jms:jme , ldrog ) , &
 147           INTENT(INOUT ) :: &
 148                          vdrog3         ! NOT USED
 149
 150    TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
 151
 152
 153 !-----------------------------------------------------------------------
 154
 155
 156 !   local variables
 157         integer :: idum, iok,iprog
 158         integer :: iregime
 159         integer :: igaschem_allowed_m1, igaschem_allowed_m2,   &
 160                    igaschem_allowed_m3, igaschem_allowed_m4
 161         integer :: igas_solver, iregime_forced
 162         integer :: i_boxtest_units_convert
 163         integer :: i_print_gasode_stats
 164         integer :: i_force_dump, mode_force_dump
 165         integer :: it, jt, kt
 166         integer :: jsolver
 167         integer :: lunerr, lunout, levdbg_err, levdbg_info
 168         integer :: mgaschem
 169
 170     real(kind=8) :: trun
 171         real :: abs_error, rel_error
 172         real :: dtchem, tchem, tstart, tstop
 173         real :: airdenbox, pressbox, tempbox
 174         real :: cair_mlc
 175         real :: h2o, ch4, oxygen, nitrogen, hydrogen
 176         real :: cboxnew(ngas_z), cboxold(ngas_z),prodrog(ldrog) !MS
 177         real :: Aperox(nperox,nperox), Bperox(nperox,nperox)
 178         real :: rk_param(nperox), rk_photo(nphoto)
 179         real :: rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
 180         real factorprog
 181         integer, dimension(2,6), save :: inforodas=0
 182         integer, dimension(6), save :: iodestatus_count=0, ioderegime_count=0
 183
 184 #ifdef CHEM_DBG_I
 185 !rcetestb diagnostics --------------------------------------------------
 186         print 93010, ' '
 187         print 93010, 'rcetestb diagnostics from cbmz_driver'
 188         print 93010, 'id, chem_opt, ktau, ktauc, julday    ',   &
 189                       id, config_flags%chem_opt, ktau, ktauc, julday
 190         print 93020, 'dtstep, dtstepc, gmt                 ',   &
 191                       dtstep, dtstepc, gmt
 192         print 93010, 'ids/e, j, k', ids, ide, jds, jde, kds, kde
 193         print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme
 194         print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte
 195         print 93010, 'num_moist, p_qv              ', num_moist, p_qv
 196         print 93010, 'num_chem, p_so2, p_ho2       ', num_chem, p_so2, p_ho2
 197         print 93010, 'p_so4aj, p_corn, p_hcl, p_mtf', p_so4aj, p_corn, p_hcl, p_mtf
 198 93010   format( a, 8(1x,i6) )
 199 93020   format( a, 8(1p,e14.6) )
 200 !rcetestb diagnostics --------------------------------------------------
 201 #endif
 202
 203
 204 !   set some control variables to their "standard for wrf-chem" values
 205         igas_solver = 1
 206         iregime_forced = -1
 207         mgaschem = +1
 208         i_boxtest_units_convert = 0
 209
 210         i_print_gasode_stats = 1
 211         mode_force_dump = 0
 212         lunerr = -1
 213         lunout = -1
 214         levdbg_err = 0
 215         levdbg_info = 15
 216
 217         abs_error = 1.0e1       ! solver absolute tolerance (molecules/cm3)
 218         rel_error = 1.0e-3      ! solver relative tolerance
 219
 220 !   set some control variables to non-standard values for testing
 221 !   force dumps for center column, every 3rd level
 222 !       mode_force_dump = +77
 223 !   force dumps for center column, 1st level
 224 !       mode_force_dump = +7
 225 !   do levdbg_info output always
 226         levdbg_info = 0
 227
 228
 229 !   following call is for boxwrf testing only
 230 !   it must be commented out for actual wrf applications
 231 !       call boxtest_get_extra_args( &
 232 !           igas_solver, iregime_forced, &
 233 !           i_boxtest_units_convert, lunerr, lunout, &
 234 !           abs_error, rel_error, trun   )
 235
 236
 237 !   currently nothing is done with vdrog3
 238 !       vdrog3(its:ite,kts:kte,jts:jte,:) = 0.0  !This is already set to zero in chem_driver.
 239
 240 !   initialize rate constant arrays to 0.0 since not all elements are set
 241 !   to valid values before some were used in loops, wig 16-Nov-2007
 242     rk_m1(:) = 0.
 243     rk_m2(:) = 0.
 244     rk_m3(:) = 0.
 245     rk_m4(:) = 0.
 246
 247 !   determine which regimes are allowed
 248 !   based on which gas species are "active"
 249         call set_gaschem_allowed_regimes( lunerr,   &
 250                 igaschem_allowed_m1, igaschem_allowed_m2,   &
 251                 igaschem_allowed_m3, igaschem_allowed_m4 )
 252
 253 !
 254 !   main loop -- do gas chemistry at each i,j,k
 255 !
 256         do 2900 jt = jts, jte
 257         do 2900 kt = kts, kte
 258         do 2900 it = its, ite
 259
 260         trun   = curr_secs              ! run time in s
 261         tchem  = mod( real(gmt*3600.0,8)+trun, 86400._8 ) ! time from 00 UTC in s
 262         dtchem = dtstepc
 263         tstart = tchem                  ! s
 264         tstop  = tstart + dtchem        ! s
 265
 266 !   skip integration for very small dtchem
 267         if ((tstop-tstart) .le. 1.0e-5) goto 2900
 268
 269
 270 !   initial species mapping from host array
 271         call mapgas_tofrom_host( 0,                    &
 272                 i_boxtest_units_convert,               &
 273                 it,jt,kt, ims,ime, jms,jme, kms,kme,   &
 274                 num_moist, num_chem, moist, chem,      &
 275                 t_phy, p_phy, rho_phy,                 &
 276                 cboxold, tempbox, pressbox, airdenbox, &
 277                 cair_mlc,                              &
 278                 h2o, ch4, oxygen, nitrogen, hydrogen   )
 279         cboxnew(:) = cboxold(:)
 280
 281 !   determine regime
 282         call selectgasregime( iregime, iregime_forced, cboxold,   &
 283                 igaschem_allowed_m1, igaschem_allowed_m2,   &
 284                 igaschem_allowed_m3, igaschem_allowed_m4 )
 285         idum = iregime
 286         if ((idum .lt. 1) .or. (idum .ge. 6)) idum = 6
 287         ioderegime_count(idum) = ioderegime_count(idum) + 1
 288         iodestatus_count(6) = iodestatus_count(6) + 1
 289
 290 !   compute rate constants
 291 !   transfer/map incoming photolysis rate contants to local array
 292       call gasphotoconstants( rk_photo,   &
 293             i_boxtest_units_convert,               &
 294             it,jt,kt, ims,ime, jms,jme, kms,kme,   &
 295             ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
 296             ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
 297             ph_ch3o2h, ph_n2o5 )
 298 !   loads Aperox and Bperox
 299         call loadperoxyparameters( Aperox, Bperox )
 300 !   calculate parameterized rate constants
 301         call peroxyrateconstants( tempbox, cboxold,   &
 302                  Aperox, Bperox, rk_param )
 303 !   calculate thermal rate constants
 304         call gasrateconstants( iregime, tempbox, cair_mlc,   &
 305                 rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 )
 306
 307 !   mode_force_dump selects a detailed dump of gaschem at either
 308 !   first ijk grid, first ij column, all ijk, or no ijk
 309         i_force_dump = 0
 310         if (mode_force_dump .eq. 1) then
 311             if ((it.eq.its) .and. (jt.eq.jts)   &
 312                             .and. (kt.eq.kts)) i_force_dump = 1
 313         else if (mode_force_dump .eq. 10) then
 314             if ((it.eq.its) .and. (jt.eq.jts)) i_force_dump = 1
 315         else if (mode_force_dump .eq. 100) then
 316             i_force_dump = 1
 317         else if (mode_force_dump .eq. 7) then
 318             if ( (it .eq.  (its+ite)/2) .and.   &
 319                  (jt .eq.  (jts+jte)/2) .and.   &
 320                  (kt .eq.          kts) ) i_force_dump = 1
 321         else if (mode_force_dump .eq. 77) then
 322             if ( (it .eq.  (its+ite)/2) .and.   &
 323                  (jt .eq.  (jts+jte)/2) .and.   &
 324                  (mod(kt-kts,3) .eq. 0) ) i_force_dump = 1
 325         end if
 326
 327
 328 !   rodas
 329         iok = 0
 330         jsolver = 0
 331         if (igas_solver .eq. 1) then
 332             jsolver = 1
 333             call gasodesolver_rodas( tstart, tstop, iok,   &
 334                 it, jt, kt, iregime,   &
 335                 mgaschem, lunerr, lunout, levdbg_err, levdbg_info,   &
 336                 i_force_dump, inforodas, iodestatus_count,   &
 337                 abs_error, rel_error, trun,   &
 338                 tempbox, pressbox, airdenbox, cboxnew, cboxold,   &
 339                 rk_m1, rk_m2, rk_m3, rk_m4 )
 340         endif
 341
 342 !   lsodes
 343         if (igas_solver.eq.2 .or. iok.le.0) then
 344             jsolver = 2
 345             call gasodesolver_lsodes( tstart, tstop, iok,   &
 346                 it, jt, kt, iregime,   &
 347                 mgaschem, lunerr, lunout, levdbg_err, levdbg_info,   &
 348                 i_force_dump, iodestatus_count,   &
 349                 abs_error, rel_error, trun,   &
 350                 tempbox, pressbox, airdenbox, cboxnew, cboxold,   &
 351                 rk_m1, rk_m2, rk_m3, rk_m4 )
 352         endif
 353
 354 !   final species mapping back to host array -- only when iok > 0
 355         if (iok .gt. 0) then
 356             call mapgas_tofrom_host( 1,                &
 357                 i_boxtest_units_convert,               &
 358                 it,jt,kt, ims,ime, jms,jme, kms,kme,   &
 359                 num_moist, num_chem, moist, chem,      &
 360                 t_phy, p_phy, rho_phy,                 &
 361                 cboxnew, tempbox, pressbox, airdenbox, &
 362                 cair_mlc,                              &
 363                 h2o, ch4, oxygen, nitrogen, hydrogen   )
 364         end if
 365
 366 !   following call is for boxwrf testing only
 367 !   it must be commented out for actual wrf applications
 368 !       call boxtest_set_extra_args( iregime, it, jt, kt )
 369
 370          if (iok .gt. 0) then
 371         call update_vdrog3_cbmz(rk_m2,rk_m3,dtstepc,cboxnew, cboxold, &
 372                                 prodrog)
 373 !prodrog units are in molecules/cm3
 374 !   chem units = (ppm)
 375         airdenbox = rho_phy(it,kt,jt)/28.966e3
 376         if (i_boxtest_units_convert .eq. 10) then
 377             airdenbox = rho_phy(it,kt,jt)
 378         end if
 379
 380
 381         factorprog = airdenbox*1.0e-6
 382         if (i_boxtest_units_convert .eq. 10) factorprog = airdenbox
 383         factorprog = factorprog*avognumkpp
 384
 385
 386 !  convert prodrog from units of molecules/cm3 to ppmv
 387          do iprog=1,ldrog
 388            prodrog(iprog)=prodrog(iprog)/factorprog
 389           enddo
 390
 391
 392 !  Now assign the prodrog values to vdrog3 array
 393           do iprog=1,ldrog
 394           vdrog3(it,kt,jt,iprog)=prodrog(iprog)
 395           vdrog3(it,kt,jt,iprog)=max(0.,vdrog3(it,kt,jt,iprog))
 396           enddo
 397
 398         endif ! iok .gt. 0
 399 !MS
 400 2900    continue
 401
 402         if (i_print_gasode_stats .gt. 0)   &
 403            call print_gasode_stats( lunout, levdbg_info,   &
 404                 inforodas, iodestatus_count, ioderegime_count )
 405         return
 406         end subroutine cbmz_driver
 407
 408
 409 !***********************************************************************
 410 ! < xx.> subr update_vdrog3_cbmz
 411 !
 412 ! purpose: updates vdrog3 concentrations for MADE-SORGAM-CBMZ package
 413 ! Written by Manish Shrivastava on 12/20/2011
 414 !-----------------------------------------------------------------------
 415
 416       subroutine update_vdrog3_cbmz(rk_m2,rk_m3,dtstepc,cboxnew, cboxold, &
 417                                 prodrog)
 418    USE module_data_sorgam, only:  ldrog
 419    USE module_data_cbmz
 420
 421       implicit none
 422
 423       real,intent(in) :: rk_m2(nrxn_m2),cboxnew(ngas_z), cboxold(ngas_z)
 424       real,intent(in) :: rk_m3(nrxn_m3)
 425       real,intent(in) :: dtstepc
 426       real,intent(inout) :: prodrog(ldrog)
 427       integer i
 428       real temp_par_new,temp_par_old,f_api,f_lim
 429
 430 !  f_api is fraction of isoprene attributed to API
 431 !  f_lim is fraction of isoprene attributed to LIM
 432
 433       f_api=0.292
 434       f_lim=0.064
 435
 436 ! First initialize prodrog
 437
 438         temp_par_new=0.0
 439         temp_par_old=0.0
 440         do i=1,ldrog
 441         prodrog(i)=0.0
 442         enddo
 443
 444         temp_par_new=temp_par_new+cboxnew(ipar_z)/7.9
 445 !since hc8 emissions are multiplied by 7.9 to calculate par, we divide by this number
 446
 447        temp_par_old=temp_par_old+cboxold(ipar_z)/7.9
 448 !since hc8 emissions are multiplied by 7.9 to calculate par, we divide by this number
 449
 450
 451
 452 ! The order of following reactions is from module_data_sorgam.F
 453 ! prodrog(1) is rxn of XYL+OH rk_m2(19)
 454       prodrog(1)=0.5*(cboxnew(ixyl_z)+cboxold(ixyl_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 455                     *rk_m2(19)*dtstepc
 456 !prodrog(2) is rxn of TOL+OH rk_m2(18)
 457      prodrog(2)=0.5*(cboxnew(itol_z)+cboxold(itol_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 458                     *rk_m2(18)*dtstepc
 459 !prodrog(3) is rxn of CRES+OH rk_m2(21)
 460      prodrog(3)=0.5*(cboxnew(icres_z)+cboxold(icres_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 461                     *rk_m2(21)*dtstepc
 462 !prodrog(4) is rxn of CRES+NO3 rk_m2(22)
 463      prodrog(4)=0.5*(cboxnew(icres_z)+cboxold(icres_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) &
 464                     *rk_m2(22)*dtstepc
 465 !prodrog(5) is rxn of HC+OH rk_m2(1)
 466      prodrog(5)=0.5*(temp_par_new+temp_par_old)*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 467                     *rk_m2(1)*dtstepc !Set to 0.0 as PAR should be replaced by HC8 or higher alkanes
 468 !prodrog(6) is rxn of OLEI+OH rk_m2(15)
 469      prodrog(6)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 470                     *rk_m2(15)*dtstepc
 471 !prodrog(7) is rxn of OLEI+NO3 rk_m2(17)
 472      prodrog(7)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) &
 473                     *rk_m2(17)*dtstepc
 474 !prodrog(8) is rxn of OLEI+O3 rk_m2(13)
 475      prodrog(8)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) &
 476                     *rk_m2(13)*dtstepc
 477 !prodrog(9) is rxn of OLET+OH rk_m2(14)
 478      prodrog(9)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 479                     *rk_m2(14)*dtstepc
 480 !prodrog(10) is rxn of OLET+NO3 rk_m2(16)
 481      prodrog(10)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) &
 482                     *rk_m2(16)*dtstepc
 483 !prodrog(11) is rxn of OLET+O3 rk_m2(12)
 484      prodrog(11)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) &
 485                     *rk_m2(12)*dtstepc
 486
 487      prodrog(12)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 488                     *rk_m3(1)*dtstepc*f_api
 489      prodrog(13)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) &
 490                     *rk_m3(3)*dtstepc*f_api
 491      prodrog(14)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) &
 492                     *rk_m3(2)*dtstepc*f_api
 493      prodrog(15)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) &
 494                     *rk_m3(1)*dtstepc*f_lim
 495      prodrog(16)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) &
 496                     *rk_m3(3)*dtstepc*f_lim
 497      prodrog(17)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) &
 498                     *rk_m3(2)*dtstepc*f_lim
 499
 500
 501
 502 ! prodrog(12) to prodrog(17) are for biogenic species API and LIM
 503 ! Since these species are not available in CBMZ prodrog(12) to prodrog(17) are not calculated
 504
 505         return
 506         end subroutine update_vdrog3_cbmz
 507
 508 !***********************************************************************
 509 ! < xx.> subr print_gasode_stats
 510 !
 511 ! purpose: writes some statistics on ode solver performance to unit lunout
 512 !
 513 !-----------------------------------------------------------------------
 514
 515         subroutine print_gasode_stats( lunout, levdbg,   &
 516                 inforodas, iodestatus_count, ioderegime_count )
 517
 518         implicit none
 519
 520 !   subr arguments
 521         integer lunout, levdbg
 522         integer inforodas(2,6), iodestatus_count(6), ioderegime_count(6)
 523
 524 !   local variables
 525         integer i, j
 526         character*80 msg
 527
 528
 529         msg = ' '
 530         call peg_debugmsg( lunout, levdbg, msg )
 531         msg = 'output from dump_cbmz_gasodeinfo'
 532         call peg_debugmsg( lunout, levdbg, msg )
 533         write(msg,9100) 'oderegime(1-6)', (ioderegime_count(i), i=1,6)
 534         call peg_debugmsg( lunout, levdbg, msg )
 535         write(msg,9100) 'odestatus(1-6)', (iodestatus_count(i), i=1,6)
 536         call peg_debugmsg( lunout, levdbg, msg )
 537
 538         write(msg,9200)   &
 539                 'inforodas(1-3)', ((inforodas(j,i), j=1,2), i=1,3)
 540         call peg_debugmsg( lunout, levdbg, msg )
 541         write(msg,9200)   &
 542                 'inforodas(4-6)', ((inforodas(j,i), j=1,2), i=4,6)
 543         call peg_debugmsg( lunout, levdbg, msg )
 544
 545 9100    format( a, 6i11 )
 546 9200    format( a, 3( i11, '--', i9.9 ) )
 547
 548         return
 549         end subroutine print_gasode_stats
 550
 551
 552
 553 !***********************************************************************
 554 ! < xx.> subr gasodesolver_rodas
 555 !
 556 ! purpose: interfaces to rodas ode solver
 557 !
 558 !-----------------------------------------------------------------------
 559
 560         subroutine gasodesolver_rodas( tstart, tstop, iok,   &
 561                 isvode, jsvode, ksvode, iregime,   &
 562                 mgaschem, lunerr, lunout, levdbg_err, levdbg_info,   &
 563                 i_force_dump, inforodas, iodestatus_count, &
 564                 abs_error, rel_error, trun,   &
 565                 tempbox, pressbox, airdenbox, cboxnew, cboxold,   &
 566                 rk_m1, rk_m2, rk_m3, rk_m4 )
 567
 568         use module_data_cbmz
 569         use module_cbmz_rodas_prep, only:                                    &
 570             cbmz_v02r01_mapconcs, cbmz_v02r01_maprates, cbmz_v02r01_torodas, &
 571             cbmz_v02r02_mapconcs, cbmz_v02r02_maprates, cbmz_v02r02_torodas, &
 572             cbmz_v02r03_mapconcs, cbmz_v02r03_maprates, cbmz_v02r03_torodas, &
 573             cbmz_v02r04_mapconcs, cbmz_v02r04_maprates, cbmz_v02r04_torodas, &
 574             cbmz_v02r05_mapconcs, cbmz_v02r05_maprates, cbmz_v02r05_torodas, &
 575             cbmz_v02r06_mapconcs, cbmz_v02r06_maprates, cbmz_v02r06_torodas
 576
 577         implicit none
 578
 579 !   subr arguments
 580     real(kind=8) :: trun
 581         integer iok, isvode, jsvode, ksvode, i_force_dump, iregime,   &
 582               mgaschem, lunerr, lunout, levdbg_err, levdbg_info
 583         integer inforodas(2,6), iodestatus_count(6)
 584         real tstart, tstop, abs_error, rel_error
 585         real tempbox, pressbox, airdenbox
 586         real cboxnew(ngas_z), cboxold(ngas_z)
 587         real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
 588
 589 !   local variables
 590         integer               :: ia, idum, idydt_sngldble, ig, l, ntot
 591         integer, save         :: nrodas_failures = 0
 592         integer, dimension(6) :: inforodas_cur
 593
 594         real hmin, hstart, taa, tzz
 595         real atolvec(ngas_z), rtolvec(ngas_z),   &
 596                 stot(ngas_z),   &
 597                 yposlimit(ngas_z), yneglimit(ngas_z)
 598         real sfixedkpp(nfixed_kppmax), rconstkpp(nreact_kppmax)
 599
 600         character*80 msg
 601
 602 !   map reaction rate constants (pegasus --> kpp)
 603 !   map concentrations (cboxold --> stot)
 604 !   dump rates (for debugging)
 605         if (iregime .eq. 1) then
 606             call cbmz_v02r01_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 607                 rconstkpp )
 608             call cbmz_v02r01_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 609
 610         else if (iregime .eq. 2) then
 611             call cbmz_v02r02_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 612                 rconstkpp )
 613             call cbmz_v02r02_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 614
 615         else if (iregime .eq. 3) then
 616             call cbmz_v02r03_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 617                 rconstkpp )
 618             call cbmz_v02r03_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 619
 620         else if (iregime .eq. 4) then
 621             call cbmz_v02r04_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 622                 rconstkpp )
 623             call cbmz_v02r04_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 624
 625         else if (iregime .eq. 5) then
 626             call cbmz_v02r05_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 627                 rconstkpp )
 628             call cbmz_v02r05_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 629
 630         else
 631             call cbmz_v02r06_maprates( rk_m1, rk_m2, rk_m3, rk_m4,   &
 632                 rconstkpp )
 633             call cbmz_v02r06_mapconcs( 0, ntot, stot, sfixedkpp, cboxold )
 634         end if
 635
 636 !   set parameters for rodas call
 637         do l = 1, ntot
 638             atolvec(l) = abs_error
 639             rtolvec(l) = rel_error
 640             yposlimit(l) = 1.0e20
 641             yneglimit(l) = -1.0e8
 642         end do
 643
 644         taa = tstart
 645         tzz = tstop
 646         hmin = 1.0e-5
 647         hstart = 60.0
 648         idydt_sngldble = 1
 649
 650 !   call rodas integrator
 651 !       subr cbmz_v02r06_torodas(
 652 !    +    ngas, taa, tzz,
 653 !    +    stot, atol, rtol, yposlimit, yneglimit,
 654 !    +    hmin, hstart,
 655 !    +    inforodas_cur, iok, lunerr, idydt_sngldble )
 656
 657         if (iregime .eq. 1) then
 658             call cbmz_v02r01_torodas(   &
 659                 ntot, taa, tzz,   &
 660                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 661                 sfixedkpp, rconstkpp,   &
 662                 hmin, hstart,   &
 663                 inforodas_cur, iok, lunerr, idydt_sngldble )
 664
 665         else if (iregime .eq. 2) then
 666             call cbmz_v02r02_torodas(   &
 667                 ntot, taa, tzz,   &
 668                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 669                 sfixedkpp, rconstkpp,   &
 670                 hmin, hstart,   &
 671                 inforodas_cur, iok, lunerr, idydt_sngldble )
 672
 673         else if (iregime .eq. 3) then
 674             call cbmz_v02r03_torodas(   &
 675                 ntot, taa, tzz,   &
 676                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 677                 sfixedkpp, rconstkpp,   &
 678                 hmin, hstart,   &
 679                 inforodas_cur, iok, lunerr, idydt_sngldble )
 680
 681         else if (iregime .eq. 4) then
 682             call cbmz_v02r04_torodas(   &
 683                 ntot, taa, tzz,   &
 684                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 685                 sfixedkpp, rconstkpp,   &
 686                 hmin, hstart,   &
 687                 inforodas_cur, iok, lunerr, idydt_sngldble )
 688
 689         else if (iregime .eq. 5) then
 690             call cbmz_v02r05_torodas(   &
 691                 ntot, taa, tzz,   &
 692                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 693                 sfixedkpp, rconstkpp,   &
 694                 hmin, hstart,   &
 695                 inforodas_cur, iok, lunerr, idydt_sngldble )
 696
 697         else
 698             call cbmz_v02r06_torodas(   &
 699                 ntot, taa, tzz,   &
 700                 stot, atolvec, rtolvec, yposlimit, yneglimit,   &
 701                 sfixedkpp, rconstkpp,   &
 702                 hmin, hstart,   &
 703                 inforodas_cur, iok, lunerr, idydt_sngldble )
 704         end if
 705
 706
 707 !   increment odeinfo counters
 708         if (iok .gt. 0) then
 709             if (inforodas_cur(6) .le. 0) then
 710                 ia = 1
 711             else
 712                 ia = 2
 713             end if
 714         else
 715             ia = 3
 716         end if
 717         iodestatus_count(ia) = iodestatus_count(ia) + 1
 718 !   do following to avoid overflow of the "inforodas" numbers
 719 !       inforodas(2,i) contains rightmost 9 digits of each inforodas number
 720 !       inforodas(1,i) contains any higher  digits of each inforodas number
 721         do ia = 1, 6
 722             idum = inforodas(2,ia) + inforodas_cur(ia)
 723             inforodas(1,ia) = inforodas(1,ia) + (idum/1000000000)
 724             inforodas(2,ia) = mod(idum, 1000000000)
 725         end do
 726
 727
 728 !   map concentrations (stot --> cboxnew)
 729         if (iregime .eq. 1) then
 730             call cbmz_v02r01_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 731         else if (iregime .eq. 2) then
 732             call cbmz_v02r02_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 733         else if (iregime .eq. 3) then
 734             call cbmz_v02r03_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 735         else if (iregime .eq. 4) then
 736             call cbmz_v02r04_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 737         else if (iregime .eq. 5) then
 738             call cbmz_v02r05_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 739         else
 740             call cbmz_v02r06_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew )
 741         end if
 742
 743
 744 !   diagnostic output if integration fails OR if i_force_dump > 0
 745         if (iok .gt. 0) then
 746             if (i_force_dump .le. 0) goto 20000
 747         else
 748             nrodas_failures = nrodas_failures + 1
 749             if (nrodas_failures .gt. 100) goto 20000
 750         end if
 751
 752         msg = ' '
 753         call peg_debugmsg( lunout, levdbg_err, msg )
 754         if (iok .gt. 0) then
 755             msg = '*** gasodesolver_rodas forced dump'
 756         else
 757             write(msg,*) '*** gasodesolver_rodas failure no.',   &
 758                 nrodas_failures
 759         end if
 760         call peg_debugmsg( lunout, levdbg_err, msg )
 761         msg = 'iregime, iok, i, j, k / t'
 762         call peg_debugmsg( lunout, levdbg_err, msg )
 763         write(msg,97010) iregime, iok, isvode, jsvode, ksvode
 764         call peg_debugmsg( lunout, levdbg_err, msg )
 765         write(msg,97020) trun
 766         call peg_debugmsg( lunout, levdbg_err, msg )
 767         msg = 'inforodas_cur(1-6) ='
 768         call peg_debugmsg( lunout, levdbg_err, msg )
 769         write(msg,97010) inforodas_cur
 770         call peg_debugmsg( lunout, levdbg_err, msg )
 771         msg =   &
 772         'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza ='
 773         call peg_debugmsg( lunout, levdbg_err, msg )
 774         write(msg,97020) tstart, tstop, abs_error, rel_error
 775         call peg_debugmsg( lunout, levdbg_err, msg )
 776         write(msg,97020) tempbox, pressbox, airdenbox, -99.0
 777         call peg_debugmsg( lunout, levdbg_err, msg )
 778
 779         idum = 0
 780         do ig = nreact_kppmax, 1, -1
 781             if ((idum .eq. 0) .and. (rconstkpp(ig) .ne. 0.0)) idum = ig
 782         end do
 783         msg = 'ngas_z, nrconst_nonzero ='
 784         call peg_debugmsg( lunout, levdbg_err, msg )
 785         write(msg,97010) ngas_z, idum
 786         call peg_debugmsg( lunout, levdbg_err, msg )
 787         msg = 'l, name, cboxold, cboxnew for l=1,ngas_z'
 788         call peg_debugmsg( lunout, levdbg_err, msg )
 789         do l = 1, ngas_z
 790             write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l)
 791             call peg_debugmsg( lunout, levdbg_err, msg )
 792         end do
 793         msg = 'rconst for i=1,nrconst_nonzero'
 794         call peg_debugmsg( lunout, levdbg_err, msg )
 795         do ia = 1, idum, 4
 796             write(msg,97020) ( rconstkpp(ig), ig = ia, min(ia+3,idum) )
 797             call peg_debugmsg( lunout, levdbg_err, msg )
 798         end do
 799
 800 97010   format( 6i12 )
 801 97020   format( 4(1pe18.10) )
 802 97030   format(( i3, 1x, a, 2(1pe18.10) ))
 803
 804
 805 !   force non-negative values
 806 20000   do l = 1, ngas_z
 807             cboxnew(l) = max( cboxnew(l), 0.0 )
 808         end do
 809
 810         return
 811         end subroutine gasodesolver_rodas
 812
 813
 814
 815 !***********************************************************************
 816 ! < xx.> subr gasodesolver_lsodes
 817 !
 818 ! purpose: interface to lsodes ode solver
 819 !
 820 ! author : Rahul A. Zaveri
 821 ! date   : May, 2000
 822 !
 823 !-----------------------------------------------------------------------
 824
 825       subroutine gasodesolver_lsodes( tstart, tstop, iok,   &
 826                 isvode, jsvode, ksvode, iregime,   &
 827                 mgaschem, lunerr, lunout, levdbg_err, levdbg_info,   &
 828                 i_force_dump, iodestatus_count, &
 829                 abs_error, rel_error, trun,   &
 830                 tempbox, pressbox, airdenbox, cboxnew, cboxold,   &
 831                 rk_m1, rk_m2, rk_m3, rk_m4 )
 832
 833       use module_data_cbmz
 834       use module_cbmz_lsodes_solver, only:  lsodes_solver, xsetf,   &
 835                                             set_lsodes_common_vars
 836       implicit none
 837
 838 !   subr arguments
 839       real(kind=8) :: trun
 840       integer i, iok, isvode, jsvode, ksvode, i_force_dump, iregime,   &
 841               mgaschem, lunerr, lunout, levdbg_err, levdbg_info
 842       integer iodestatus_count(6)
 843       real tstart, tstop, abs_error, rel_error
 844       real tempbox, pressbox, airdenbox
 845       real cboxnew(ngas_z), cboxold(ngas_z)
 846       real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
 847
 848 ! lsodes parameters and local variables
 849       integer itoler, itask, iopt, mf, lwm, nrdim, nidim
 850       integer nruserpar, niuserpar
 851       parameter( itoler = 1, itask = 1, iopt = 1, mf= 222 )
 852       parameter( lwm = 3*ngas_tot*ngas_tot + 12*ngas_tot )
 853       parameter( nrdim = 20 + 9*ngas_tot + lwm )
 854       parameter( nidim = 31 + ngas_tot + ngas_tot*ngas_tot )
 855       parameter( nruserpar = 5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4)
 856       parameter( niuserpar = ngas_z + 1 )
 857
 858       integer ia, idum, ig, ioffset, istate, iwork(nidim), l
 859       integer ntotvec(1), iuserpar(niuserpar)
 860       integer indx(ngas_z)
 861       integer, save :: iflagout = 0
 862       integer, save :: nlsodes_failures = 0
 863
 864       real dtchem, rwork(nrdim), stot(ngas_tot)
 865       real atolvec(1), rtolvec(1), ruserpar(nruserpar)
 866
 867       character*80 msg
 868
 869
 870
 871       iok = 1                           ! reset
 872
 873       call set_lsodes_common_vars()
 874
 875 ! sets gas species indices for iregime
 876       call setgasindices( iregime, indx )
 877
 878 ! map cboxold --> stot
 879       call mapgasspecies( cboxold, stot, 0, iregime, indx )
 880
 881 !----------------------------------------------------------------------
 882 ! set number of species (ntot) for the selected regime for LSODES
 883       if      (iregime .eq. 1) then
 884           ntotvec(1) = ngas_r1
 885       else if (iregime .eq. 2) then
 886           ntotvec(1) = ngas_r2
 887       else if (iregime .eq. 3) then
 888           ntotvec(1) = ngas_r3
 889       else if (iregime .eq. 4) then
 890           ntotvec(1) = ngas_r4
 891       else if (iregime .eq. 5) then
 892           ntotvec(1) = ngas_r5
 893       else
 894           ntotvec(1) = ngas_r6
 895       end if
 896
 897 100   continue
 898
 899 ! set other LSODES parameters...
 900       iwork(6) = 1000           ! max iterations for a time step
 901       iwork(7) = 1
 902       istate   = 1
 903       rwork(6) = dtchem
 904       if(iflagout.eq.0)then
 905           call xsetf(iflagout)
 906       endif
 907
 908       atolvec(1) = abs_error
 909       rtolvec(1) = rel_error
 910
 911       do ig = 1, 5
 912           ruserpar(ig) = ig*7.0
 913       end do
 914       ruserpar(1) = cboxold(ih2o_z)
 915       ruserpar(2) = cboxold(ich4_z)
 916       ruserpar(3) = cboxold(io2_z)
 917       ruserpar(4) = cboxold(in2_z)
 918       ruserpar(5) = cboxold(ih2_z)
 919       ioffset = 5
 920       do ig = 1, nrxn_m1
 921           ruserpar(ioffset+ig) = rk_m1(ig)
 922       end do
 923       ioffset = ioffset + nrxn_m1
 924       do ig = 1, nrxn_m2
 925           ruserpar(ioffset+ig) = rk_m2(ig)
 926       end do
 927       ioffset = ioffset + nrxn_m2
 928       do ig = 1, nrxn_m3
 929           ruserpar(ioffset+ig) = rk_m3(ig)
 930       end do
 931       ioffset = ioffset + nrxn_m3
 932       do ig = 1, nrxn_m4
 933           ruserpar(ioffset+ig) = rk_m4(ig)
 934       end do
 935
 936       iuserpar(1) = iregime
 937       do ig = 1, ngas_z
 938           iuserpar(1+ig) = indx(ig)
 939       end do
 940
 941       call lsodes_solver(   &
 942                 gasode_cbmz, ntotvec, stot, tstart, tstop,   &
 943                 itoler, rtolvec, atolvec, itask, istate, iopt,   &
 944                 rwork, nrdim, iwork, nidim, jcs, mf,   &
 945                 ruserpar, nruserpar, iuserpar, niuserpar )
 946
 947       if (istate .le. 0) iok = -1
 948
 949
 950 !   increment odeinfo counters
 951       if (iok .gt. 0) then
 952           ia = 4
 953       else
 954           ia = 5
 955       end if
 956       iodestatus_count(ia) = iodestatus_count(ia) + 1
 957
 958
 959 ! map stot --> cboxnew
 960         call mapgasspecies( cboxnew, stot, 1, iregime, indx )
 961
 962
 963 ! do diagnostic output if integration fails OR if i_force_dump > 0
 964         if (iok .gt. 0) then
 965             if (i_force_dump .le. 0) goto 20000
 966         else
 967             nlsodes_failures = nlsodes_failures + 1
 968         end if
 969
 970         msg = ' '
 971         call peg_debugmsg( lunout, levdbg_err, msg )
 972         if (iok .gt. 0) then
 973             msg = '*** gasodesolver_lsodes forced dump'
 974         else
 975             write(msg,*) '*** gasodesolver_lsodes failure no.',   &
 976                 nlsodes_failures
 977         end if
 978         call peg_debugmsg( lunout, levdbg_err, msg )
 979         msg = 'iregime, iok, i, j, k / t'
 980         call peg_debugmsg( lunout, levdbg_err, msg )
 981         write(msg,97010) iregime, iok, isvode, jsvode, ksvode
 982         call peg_debugmsg( lunout, levdbg_err, msg )
 983         write(msg,97020) trun
 984         call peg_debugmsg( lunout, levdbg_err, msg )
 985         if (nlsodes_failures .gt. 1000) then
 986             write(msg,*) '*** exceeded lsodes failure limit =', 1000
 987             call peg_debugmsg( lunout, levdbg_err, msg )
 988             call peg_error_fatal( lunerr, msg )
 989         end if
 990         if (nlsodes_failures .gt. 100) goto 20000
 991
 992         write(msg,*) 'istate -', istate
 993         call peg_debugmsg( lunout, levdbg_err, msg )
 994         msg =   &
 995         'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza ='
 996         call peg_debugmsg( lunout, levdbg_err, msg )
 997         write(msg,97020) tstart, tstop, abs_error, rel_error
 998         call peg_debugmsg( lunout, levdbg_err, msg )
 999         write(msg,97020) tempbox, pressbox, airdenbox, -99.0
1000         call peg_debugmsg( lunout, levdbg_err, msg )
1001
1002         idum = nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4
1003         msg = 'ngas_z, nrconst_m1+m2+m3+m4 ='
1004         call peg_debugmsg( lunout, levdbg_err, msg )
1005         write(msg,97010) ngas_z, idum
1006         call peg_debugmsg( lunout, levdbg_err, msg )
1007         msg = 'l, name, cboxold, cboxnew for l=1,ngas_z'
1008         call peg_debugmsg( lunout, levdbg_err, msg )
1009         do l = 1, ngas_z
1010             write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l)
1011             call peg_debugmsg( lunout, levdbg_err, msg )
1012         end do
1013         msg = 'rconst for i=1,nrconst_nonzero'
1014         call peg_debugmsg( lunout, levdbg_err, msg )
1015         do ia = 1, nrxn_m1, 4
1016             write(msg,97020) ( rk_m1(ig), ig = ia, min(ia+3,nrxn_m1) )
1017             call peg_debugmsg( lunout, levdbg_err, msg )
1018         end do
1019         do ia = 1, nrxn_m2, 4
1020             write(msg,97020) ( rk_m2(ig), ig = ia, min(ia+3,nrxn_m2) )
1021             call peg_debugmsg( lunout, levdbg_err, msg )
1022         end do
1023         do ia = 1, nrxn_m3, 4
1024             write(msg,97020) ( rk_m3(ig), ig = ia, min(ia+3,nrxn_m3) )
1025             call peg_debugmsg( lunout, levdbg_err, msg )
1026         end do
1027         do ia = 1, nrxn_m4, 4
1028             write(msg,97020) ( rk_m4(ig), ig = ia, min(ia+3,nrxn_m4) )
1029             call peg_debugmsg( lunout, levdbg_err, msg )
1030         end do
1031
1032 97010   format( 6i12 )
1033 97020   format( 4(1pe18.10) )
1034 97030   format(( i3, 1x, a, 2(1pe18.10) ))
1035
1036
1037 !   force non-negative values
1038 20000   do l = 1, ngas_z
1039             cboxnew(l) = max( cboxnew(l), 0.0 )
1040         end do
1041
1042       return
1043       end subroutine gasodesolver_lsodes
1044
1045
1046
1047 !***********************************************************************
1048 ! < 2.> subr selectgasregime
1049 !
1050 ! purpose: selects an optimum combination of gas-phase
1051 !          mechanisms based on sensitivity of [OH]
1052 !          concentrations to some lumped structure
1053 !          hydrocarbon groups concentrations and [DMS]
1054 !          concentration.
1055 !
1056 ! input : cbox      = full species concentrations array (mol/cc)
1057 !
1058 ! output: iregime   = 1     : com
1059 !                   = 2     : com + urb
1060 !                   = 3     : com + urb + bio
1061 !                   = 4     : com + mar
1062 !                   = 5     : com + urb + mar
1063 !                   = 6     : com + urb + bio + mar
1064 !         ntot      = number of gas-phase species in the selected mechanism
1065 !
1066 ! author: Rahul A. Zaveri
1067 ! date  : April 2000
1068 !
1069 !---------------------------------------------------------------------
1070
1071       subroutine selectgasregime( iregime, iregime_forced, cbox,   &
1072               igaschem_allowed_m1, igaschem_allowed_m2,   &
1073               igaschem_allowed_m3, igaschem_allowed_m4 )
1074
1075       use module_data_cbmz
1076       implicit none
1077
1078 !   subr arguments
1079       integer iregime, iregime_forced
1080       integer igaschem_allowed_m1, igaschem_allowed_m2,   &
1081               igaschem_allowed_m3, igaschem_allowed_m4
1082       real cbox(ngas_z)
1083
1084 !   local variables
1085       integer iwork(6)
1086       integer m_m1, m_m2, m_m3, m_m4
1087       real cut_molecpcc
1088
1089
1090       cut_molecpcc = 5.e+6              ! molecules/cc
1091
1092 ! initialize mechanism flags
1093       m_m1 = 1  ! 1 (always)
1094       m_m2 = 0  ! 0 or 1
1095       m_m3 = 0  ! 0 or 2
1096       m_m4 = 0  ! 0 or 3
1097
1098       if (igaschem_allowed_m2 .gt. 0) then
1099           if (cbox(ipar_z     ) .gt. cut_molecpcc) m_m2 = 1
1100           if (cbox(iaone_z    ) .gt. cut_molecpcc) m_m2 = 1
1101           if (cbox(imgly_z    ) .gt. cut_molecpcc) m_m2 = 1
1102           if (cbox(ieth_z     ) .gt. cut_molecpcc) m_m2 = 1
1103           if (cbox(iolet_z    ) .gt. cut_molecpcc) m_m2 = 1
1104           if (cbox(iolei_z    ) .gt. cut_molecpcc) m_m2 = 1
1105           if (cbox(ixyl_z     ) .gt. cut_molecpcc) m_m2 = 1
1106           if (cbox(icres_z    ) .gt. cut_molecpcc) m_m2 = 1
1107           if (cbox(ito2_z     ) .gt. cut_molecpcc) m_m2 = 1
1108           if (cbox(icro_z     ) .gt. cut_molecpcc) m_m2 = 1
1109           if (cbox(iopen_z    ) .gt. cut_molecpcc) m_m2 = 1
1110           if (cbox(ionit_z    ) .gt. cut_molecpcc) m_m2 = 1
1111           if (cbox(irooh_z    ) .gt. cut_molecpcc) m_m2 = 1
1112           if (cbox(iro2_z     ) .gt. cut_molecpcc) m_m2 = 1
1113           if (cbox(iano2_z    ) .gt. cut_molecpcc) m_m2 = 1
1114           if (cbox(inap_z     ) .gt. cut_molecpcc) m_m2 = 1
1115           if (cbox(ixo2_z     ) .gt. cut_molecpcc) m_m2 = 1
1116           if (cbox(ixpar_z    ) .gt. cut_molecpcc) m_m2 = 1
1117       end if
1118
1119       if (igaschem_allowed_m3 .gt. 0) then
1120           if (cbox(iisop_z    ) .gt. cut_molecpcc) m_m3 = 2
1121           if (cbox(iisoprd_z  ) .gt. cut_molecpcc) m_m3 = 2
1122           if (cbox(iisopp_z   ) .gt. cut_molecpcc) m_m3 = 2
1123           if (cbox(iisopn_z   ) .gt. cut_molecpcc) m_m3 = 2
1124           if (cbox(iisopo2_z  ) .gt. cut_molecpcc) m_m3 = 2
1125       end if
1126
1127       if (igaschem_allowed_m4 .gt. 0) then
1128           if (cbox(idms_z        ) .gt. cut_molecpcc) m_m4 = 3
1129           if (cbox(imsa_z        ) .gt. cut_molecpcc) m_m4 = 3
1130           if (cbox(idmso_z       ) .gt. cut_molecpcc) m_m4 = 3
1131           if (cbox(idmso2_z      ) .gt. cut_molecpcc) m_m4 = 3
1132           if (cbox(ich3so2h_z    ) .gt. cut_molecpcc) m_m4 = 3
1133           if (cbox(ich3sch2oo_z  ) .gt. cut_molecpcc) m_m4 = 3
1134           if (cbox(ich3so2_z     ) .gt. cut_molecpcc) m_m4 = 3
1135           if (cbox(ich3so3_z     ) .gt. cut_molecpcc) m_m4 = 3
1136           if (cbox(ich3so2oo_z   ) .gt. cut_molecpcc) m_m4 = 3
1137           if (cbox(ich3so2ch2oo_z) .gt. cut_molecpcc) m_m4 = 3
1138           if (cbox(imtf_z        ) .gt. cut_molecpcc) m_m4 = 3
1139       end if
1140
1141       iregime = m_m1 + m_m2*((2-m_m3)/2) + m_m3 + m_m4
1142
1143 !   force iregime = iregime_forced
1144       if ((iregime_forced .ge. 1) .and. (iregime_forced .le. 6))   &
1145           iregime = iregime_forced
1146
1147       return
1148       end subroutine selectgasregime
1149
1150
1151
1152 !***********************************************************************
1153 ! < 3.> subr setgasindices
1154 !
1155 ! purpose: sets gas species indices
1156 !
1157 ! author : Rahul A. Zaveri
1158 ! date   : May, 2000
1159 !
1160 !-----------------------------------------------------------------------
1161
1162       subroutine setgasindices( iregime, indx )
1163
1164       use module_data_cbmz
1165       implicit none
1166
1167 !   subr arguments
1168       integer iregime, indx(ngas_z)
1169
1170 !   local variables
1171       integer ilast
1172
1173       ilast = 0
1174       indx(:) = -999888777
1175
1176       goto (1,2,3,4,5,6), iregime
1177
1178 1     call setgasindex_m1( ilast, indx )                ! regime 1
1179       return
1180
1181 2     call setgasindex_m1( ilast, indx )                ! regime 2
1182       call setgasindex_m2( ilast, indx )
1183       return
1184
1185 3     call setgasindex_m1( ilast, indx )                ! regime 3
1186       call setgasindex_m2( ilast, indx )
1187       call setgasindex_m3( ilast, indx )
1188       return
1189
1190 4     call setgasindex_m1( ilast, indx )                ! regime 4
1191       call setgasindex_m4( ilast, indx )
1192       return
1193
1194 5     call setgasindex_m1( ilast, indx )                ! regime 5
1195       call setgasindex_m2( ilast, indx )
1196       call setgasindex_m4( ilast, indx )
1197       return
1198
1199 6     call setgasindex_m1( ilast, indx )                ! regime 6
1200       call setgasindex_m2( ilast, indx )
1201       call setgasindex_m3( ilast, indx )
1202       call setgasindex_m4( ilast, indx )
1203       return
1204
1205       end subroutine setgasindices
1206
1207
1208
1209
1210 !***********************************************************************
1211 ! < 4.> subr gasrateconstants
1212 !
1213 ! purpose: calls regime-dependent subrs for calculating
1214 !          gas-phase thermal reaction rate constants
1215 !
1216 ! author : Rahul A. Zaveri
1217 ! date   : May, 2000
1218 !
1219 !-----------------------------------------------------------------------
1220
1221       subroutine gasrateconstants( iregime, tempbox, cair_mlc,   &
1222                 rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 )
1223
1224       use module_data_cbmz
1225       implicit none
1226
1227 !   subr arguments
1228       integer iregime
1229       real tempbox, cair_mlc
1230       real rk_photo(nphoto), rk_param(nperox)
1231       real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4)
1232
1233
1234 !   iregime=1 --> do m1
1235 !   iregime=2 --> do m1, m2
1236 !   iregime=3 --> do m1, m2, m3
1237 !   iregime=4 --> do m1, --, --, m4
1238 !   iregime=5 --> do m1, m2, --, m4
1239 !   iregime=6 --> do m1, m2, m3, m4
1240
1241       call gasthermrk_m1( tempbox, cair_mlc,   &
1242                           rk_photo, rk_param, rk_m1, rk_m2 )
1243
1244       if ((iregime .eq. 2) .or.   &
1245           (iregime .eq. 3) .or.   &
1246           (iregime .eq. 5) .or.   &
1247           (iregime .eq. 6))       &
1248           call gasthermrk_m2( tempbox, cair_mlc,   &
1249                               rk_photo, rk_param, rk_m2 )
1250
1251       if ((iregime .eq. 3) .or.   &
1252           (iregime .eq. 6))       &
1253           call gasthermrk_m3( tempbox, cair_mlc,   &
1254                               rk_photo, rk_param, rk_m3 )
1255
1256       if ((iregime .eq. 4) .or.   &
1257           (iregime .eq. 5) .or.   &
1258           (iregime .eq. 6))       &
1259           call gasthermrk_m4( tempbox, cair_mlc,   &
1260                               rk_photo, rk_param, rk_m4 )
1261
1262       return
1263       end subroutine gasrateconstants
1264
1265
1266
1267 !***********************************************************************
1268 ! < 3.> subr setgasindex_m1
1269 !
1270 ! purpose: defines gas species indices for regime 1
1271 !
1272 ! author : Rahul A. Zaveri
1273 ! date   : December, 1998
1274 !
1275 !-----------------------------------------------------------------------
1276
1277       subroutine setgasindex_m1( ilast, indx )
1278
1279       use module_data_cbmz
1280       implicit none
1281
1282 !   subr arguments
1283       integer ilast, indx(ngas_z)
1284
1285
1286       indx(ino_z)       = 1
1287       indx(ino2_z)      = 2
1288       indx(ino3_z)      = 3
1289       indx(in2o5_z)     = 4
1290       indx(ihono_z)     = 5
1291       indx(ihno3_z)     = 6
1292       indx(ihno4_z)     = 7
1293       indx(io3_z)       = 8
1294       indx(io1d_z)      = 9
1295       indx(io3p_z)      = 10
1296       indx(ioh_z)       = 11
1297       indx(iho2_z)      = 12
1298       indx(ih2o2_z)     = 13
1299       indx(ico_z)       = 14
1300       indx(iso2_z)      = 15
1301       indx(ih2so4_z)    = 16
1302       indx(inh3_z)      = 17
1303       indx(ihcl_z)      = 18
1304       indx(ic2h6_z)     = 19
1305       indx(ich3o2_z)    = 20
1306       indx(iethp_z)     = 21
1307       indx(ihcho_z)     = 22
1308       indx(ich3oh_z)    = 23
1309       indx(ic2h5oh_z)   = 24
1310       indx(ich3ooh_z)   = 25
1311       indx(iethooh_z)   = 26
1312       indx(iald2_z)     = 27
1313       indx(ihcooh_z)    = 28
1314       indx(ircooh_z)    = 29
1315       indx(ic2o3_z)     = 30
1316       indx(ipan_z)      = 31
1317
1318       ilast = indx(ipan_z)
1319
1320       return
1321       end subroutine setgasindex_m1
1322
1323
1324
1325 !***********************************************************************
1326 ! < 4.> subr setgasindex_m2
1327 !
1328 ! purpose: defines gas species indices for regime 2
1329 !
1330 ! author : Rahul A. Zaveri
1331 ! date   : December, 1998
1332 !
1333 !-----------------------------------------------------------------------
1334       subroutine setgasindex_m2( ilast, indx )
1335
1336       use module_data_cbmz
1337       implicit none
1338
1339 !   subr arguments
1340       integer ilast, indx(ngas_z)
1341
1342
1343       indx(ipar_z)      = ilast + 1
1344       indx(iaone_z)     = ilast + 2
1345       indx(imgly_z)     = ilast + 3
1346       indx(ieth_z)      = ilast + 4
1347       indx(iolet_z)     = ilast + 5
1348       indx(iolei_z)     = ilast + 6
1349       indx(itol_z)      = ilast + 7
1350       indx(ixyl_z)      = ilast + 8
1351       indx(icres_z)     = ilast + 9
1352       indx(ito2_z)      = ilast + 10
1353       indx(icro_z)      = ilast + 11
1354       indx(iopen_z)     = ilast + 12
1355       indx(ionit_z)     = ilast + 13
1356 !     indx(ipan_z)      = ilast + 14
1357 !     indx(ircooh_z)    = ilast + 15
1358       indx(irooh_z)     = ilast + 14
1359 !     indx(ic2o3_z)     = ilast + 17
1360       indx(iro2_z)      = ilast + 15
1361       indx(iano2_z)     = ilast + 16
1362       indx(inap_z)      = ilast + 17
1363       indx(ixo2_z)      = ilast + 18
1364       indx(ixpar_z)     = ilast + 19
1365
1366       ilast = indx(ixpar_z)
1367
1368       return
1369       end subroutine setgasindex_m2
1370
1371
1372
1373 !***********************************************************************
1374 ! < 5.> subr setgasindex_m3
1375 !
1376 ! purpose: defines gas species indices for regime 3
1377 !
1378 ! author : Rahul A. Zaveri
1379 ! date   : December, 1998
1380 !
1381 !-----------------------------------------------------------------------
1382       subroutine setgasindex_m3( ilast, indx )
1383
1384       use module_data_cbmz
1385       implicit none
1386
1387 !   subr arguments
1388       integer ilast, indx(ngas_z)
1389
1390
1391       indx(iisop_z)     = ilast + 1
1392       indx(iisoprd_z)   = ilast + 2
1393       indx(iisopp_z)    = ilast + 3
1394       indx(iisopn_z)    = ilast + 4
1395       indx(iisopo2_z)   = ilast + 5
1396
1397       ilast = indx(iisopo2_z)
1398
1399       return
1400       end subroutine setgasindex_m3
1401
1402
1403
1404 !***********************************************************************
1405 ! < 6.> subr setgasindex_m4
1406 !
1407 ! purpose: defines gas species indices for regime 4
1408 !
1409 ! author : Rahul A. Zaveri
1410 ! date   : December, 1998
1411 !
1412 !-----------------------------------------------------------------------
1413
1414       subroutine setgasindex_m4( ilast, indx )
1415
1416       use module_data_cbmz
1417       implicit none
1418
1419 !   subr arguments
1420       integer ilast, indx(ngas_z)
1421
1422
1423 !
1424       indx(idms_z)              = ilast + 1
1425       indx(imsa_z)              = ilast + 2
1426       indx(idmso_z)             = ilast + 3
1427       indx(idmso2_z)            = ilast + 4
1428       indx(ich3so2h_z)          = ilast + 5
1429       indx(ich3sch2oo_z)        = ilast + 6
1430       indx(ich3so2_z)           = ilast + 7
1431       indx(ich3so3_z)           = ilast + 8
1432       indx(ich3so2oo_z) = ilast + 9
1433       indx(ich3so2ch2oo_z)      = ilast + 10
1434       indx(imtf_z)              = ilast + 11
1435
1436       ilast = indx(imtf_z)
1437
1438       return
1439       end subroutine setgasindex_m4
1440
1441
1442
1443 !***********************************************************************
1444 ! < 9.> subr mapgasspecies
1445 !
1446 ! purpose: map gas species between stot and cbox arrays
1447 !
1448 ! author : Rahul A. Zaveri
1449 ! date   : December, 1998
1450 !
1451 ! ----------------------------------------------------------------------
1452
1453       subroutine mapgasspecies( cbox, stot, imap, iregime, indx )
1454
1455       use module_data_cbmz
1456       implicit none
1457
1458 !   subr arguments
1459       integer imap, iregime, indx(ngas_z)
1460       real cbox(ngas_z)
1461       real stot(ngas_tot)
1462 !
1463 !
1464       goto (1,2,3,4,5,6), iregime
1465 !
1466 !
1467 1     call mapgas_m1( cbox, stot, imap, indx )
1468       return
1469 !
1470 !
1471 2     call mapgas_m1( cbox, stot, imap, indx )
1472       call mapgas_m2( cbox, stot, imap, indx )
1473       return
1474 !
1475 !
1476 3     call mapgas_m1( cbox, stot, imap, indx )
1477       call mapgas_m2( cbox, stot, imap, indx )
1478       call mapgas_m3( cbox, stot, imap, indx )
1479       return
1480 !
1481 !
1482 4     call mapgas_m1( cbox, stot, imap, indx )
1483       call mapgas_m4( cbox, stot, imap, indx )
1484       return
1485 !
1486 !
1487 5     call mapgas_m1( cbox, stot, imap, indx )
1488       call mapgas_m2( cbox, stot, imap, indx )
1489       call mapgas_m4( cbox, stot, imap, indx )
1490       return
1491 !
1492 !
1493 6     call mapgas_m1( cbox, stot, imap, indx )
1494       call mapgas_m2( cbox, stot, imap, indx )
1495       call mapgas_m3( cbox, stot, imap, indx )
1496       call mapgas_m4( cbox, stot, imap, indx )
1497       return
1498
1499       end subroutine mapgasspecies
1500
1501
1502
1503 !***********************************************************************
1504 ! <10.> subr mapgas_m1
1505 !
1506 ! purpose: maps gas species between stot and cbox arrays for mechanism 1
1507 !
1508 ! author : Rahul A. Zaveri
1509 ! date   : December, 1998
1510 !
1511 ! ----------------------------------------------------------------------
1512
1513       subroutine mapgas_m1( cbox, stot, imap, indx )
1514
1515       use module_data_cbmz
1516       implicit none
1517
1518 !   subr arguments
1519       integer imap, indx(ngas_z)
1520       real cbox(ngas_z)
1521       real stot(ngas_tot)
1522
1523       if(imap.eq.0)then    ! map cbox --> stot (both molec/cc)
1524       stot(indx(ino_z)) = cbox(ino_z)
1525       stot(indx(ino2_z))        = cbox(ino2_z)
1526       stot(indx(ino3_z))        = cbox(ino3_z)
1527       stot(indx(in2o5_z))       = cbox(in2o5_z)
1528       stot(indx(ihono_z))       = cbox(ihono_z)
1529       stot(indx(ihno3_z))       = cbox(ihno3_z)
1530       stot(indx(ihno4_z))       = cbox(ihno4_z)
1531       stot(indx(io3_z)) = cbox(io3_z)
1532       stot(indx(io1d_z))        = cbox(io1d_z)
1533       stot(indx(io3p_z))        = cbox(io3p_z)
1534       stot(indx(ioh_z)) = cbox(ioh_z)
1535       stot(indx(iho2_z))        = cbox(iho2_z)
1536       stot(indx(ih2o2_z))       = cbox(ih2o2_z)
1537       stot(indx(ico_z)) = cbox(ico_z)
1538       stot(indx(iso2_z))        = cbox(iso2_z)
1539       stot(indx(ih2so4_z))      = cbox(ih2so4_z)
1540       stot(indx(inh3_z))        = cbox(inh3_z)
1541       stot(indx(ihcl_z))        = cbox(ihcl_z)
1542       stot(indx(ic2h6_z))       = cbox(ic2h6_z)
1543       stot(indx(ich3o2_z))      = cbox(ich3o2_z)
1544       stot(indx(iethp_z))       = cbox(iethp_z)
1545       stot(indx(ihcho_z))       = cbox(ihcho_z)
1546       stot(indx(ich3oh_z))      = cbox(ich3oh_z)
1547       stot(indx(ic2h5oh_z))     = cbox(ic2h5oh_z)
1548       stot(indx(ich3ooh_z))     = cbox(ich3ooh_z)
1549       stot(indx(iethooh_z))     = cbox(iethooh_z)
1550       stot(indx(iald2_z))       = cbox(iald2_z)
1551       stot(indx(ihcooh_z))      = cbox(ihcooh_z)
1552       stot(indx(ircooh_z))      = cbox(ircooh_z)
1553       stot(indx(ic2o3_z))       = cbox(ic2o3_z)
1554       stot(indx(ipan_z))        = cbox(ipan_z)
1555 !
1556       else                 ! map stot --> cbox (both molec/cc)
1557       cbox(ino_z)       = stot(indx(ino_z))
1558       cbox(ino2_z)      = stot(indx(ino2_z))
1559       cbox(ino3_z)      = stot(indx(ino3_z))
1560       cbox(in2o5_z)     = stot(indx(in2o5_z))
1561       cbox(ihono_z)     = stot(indx(ihono_z))
1562       cbox(ihno3_z)     = stot(indx(ihno3_z))
1563       cbox(ihno4_z)     = stot(indx(ihno4_z))
1564       cbox(io3_z)       = stot(indx(io3_z))
1565       cbox(io1d_z)      = stot(indx(io1d_z))
1566       cbox(io3p_z)      = stot(indx(io3p_z))
1567       cbox(ioh_z)       = stot(indx(ioh_z))
1568       cbox(iho2_z)      = stot(indx(iho2_z))
1569       cbox(ih2o2_z)     = stot(indx(ih2o2_z))
1570       cbox(ico_z)       = stot(indx(ico_z))
1571       cbox(iso2_z)      = stot(indx(iso2_z))
1572       cbox(ih2so4_z)    = stot(indx(ih2so4_z))
1573       cbox(inh3_z)      = stot(indx(inh3_z))
1574       cbox(ihcl_z)      = stot(indx(ihcl_z))
1575       cbox(ic2h6_z)     = stot(indx(ic2h6_z))
1576       cbox(ich3o2_z)    = stot(indx(ich3o2_z))
1577       cbox(iethp_z)     = stot(indx(iethp_z))
1578       cbox(ihcho_z)     = stot(indx(ihcho_z))
1579       cbox(ich3oh_z)    = stot(indx(ich3oh_z))
1580       cbox(ic2h5oh_z)   = stot(indx(ic2h5oh_z))
1581       cbox(ich3ooh_z)   = stot(indx(ich3ooh_z))
1582       cbox(iethooh_z)   = stot(indx(iethooh_z))
1583       cbox(iald2_z)     = stot(indx(iald2_z))
1584       cbox(ihcooh_z)    = stot(indx(ihcooh_z))
1585       cbox(ircooh_z)    = stot(indx(ircooh_z))
1586       cbox(ic2o3_z)     = stot(indx(ic2o3_z))
1587       cbox(ipan_z)      = stot(indx(ipan_z))
1588       endif
1589
1590       return
1591       end subroutine mapgas_m1
1592
1593
1594
1595 !***********************************************************************
1596 ! <11.> subr mapgas_m2
1597 !
1598 ! purpose: maps gas species between stot and cbox arrays for mechanism 2
1599 !
1600 ! author : Rahul A. Zaveri
1601 ! date   : December, 1998
1602 !
1603 ! ----------------------------------------------------------------------
1604
1605       subroutine mapgas_m2( cbox, stot, imap, indx )
1606
1607       use module_data_cbmz
1608       implicit none
1609
1610 !   subr arguments
1611       integer imap, indx(ngas_z)
1612       real cbox(ngas_z)
1613       real stot(ngas_tot)
1614
1615       if(imap.eq.0)then    ! map cbox --> stot (both molec/cc)
1616       stot(indx(ipar_z))        = cbox(ipar_z)
1617       stot(indx(iaone_z))       = cbox(iaone_z)
1618       stot(indx(imgly_z))       = cbox(imgly_z)
1619       stot(indx(ieth_z))        = cbox(ieth_z)
1620       stot(indx(iolet_z))       = cbox(iolet_z)
1621       stot(indx(iolei_z))       = cbox(iolei_z)
1622       stot(indx(itol_z))        = cbox(itol_z)
1623       stot(indx(ixyl_z))        = cbox(ixyl_z)
1624       stot(indx(icres_z))       = cbox(icres_z)
1625       stot(indx(ito2_z))        = cbox(ito2_z)
1626       stot(indx(icro_z))        = cbox(icro_z)
1627       stot(indx(iopen_z))       = cbox(iopen_z)
1628       stot(indx(ionit_z))       = cbox(ionit_z)
1629 !     stot(indx(ipan_z))        = cbox(ipan_z)
1630 !     stot(indx(ircooh_z))      = cbox(ircooh_z)
1631       stot(indx(irooh_z))       = cbox(irooh_z)
1632 !     stot(indx(ic2o3_z))       = cbox(ic2o3_z)
1633       stot(indx(iro2_z))        = cbox(iro2_z)
1634       stot(indx(iano2_z))       = cbox(iano2_z)
1635       stot(indx(inap_z))        = cbox(inap_z)
1636       stot(indx(ixo2_z))        = cbox(ixo2_z)
1637       stot(indx(ixpar_z))       = cbox(ixpar_z)
1638 !
1639       else                 ! map stot --> cbox (both molec/cc)
1640       cbox(ipar_z)      = stot(indx(ipar_z))
1641       cbox(iaone_z)     = stot(indx(iaone_z))
1642       cbox(imgly_z)     = stot(indx(imgly_z))
1643       cbox(ieth_z)      = stot(indx(ieth_z))
1644       cbox(iolet_z)     = stot(indx(iolet_z))
1645       cbox(iolei_z)     = stot(indx(iolei_z))
1646       cbox(itol_z)      = stot(indx(itol_z))
1647       cbox(ixyl_z)      = stot(indx(ixyl_z))
1648       cbox(icres_z)     = stot(indx(icres_z))
1649       cbox(ito2_z)      = stot(indx(ito2_z))
1650       cbox(icro_z)      = stot(indx(icro_z))
1651       cbox(iopen_z)     = stot(indx(iopen_z))
1652       cbox(ionit_z)     = stot(indx(ionit_z))
1653 !     cbox(ipan_z)      = stot(indx(ipan_z))
1654 !     cbox(ircooh_z)    = stot(indx(ircooh_z))
1655       cbox(irooh_z)     = stot(indx(irooh_z))
1656 !     cbox(ic2o3_z)     = stot(indx(ic2o3_z))
1657       cbox(iro2_z)      = stot(indx(iro2_z))
1658       cbox(iano2_z)     = stot(indx(iano2_z))
1659       cbox(inap_z)      = stot(indx(inap_z))
1660       cbox(ixo2_z)      = stot(indx(ixo2_z))
1661       cbox(ixpar_z)     = stot(indx(ixpar_z))
1662       endif
1663
1664       return
1665       end subroutine mapgas_m2
1666
1667
1668
1669 !***********************************************************************
1670 ! <12.> subr mapgas_m3
1671 !
1672 ! purpose: maps gas species between stot and cbox arrays for mechanism 3
1673 !
1674 ! author : Rahul A. Zaveri
1675 ! date   : December, 1998
1676 !
1677 ! ----------------------------------------------------------------------
1678
1679       subroutine mapgas_m3( cbox, stot, imap, indx )
1680
1681       use module_data_cbmz
1682       implicit none
1683
1684 !   subr arguments
1685       integer imap, indx(ngas_z)
1686       real cbox(ngas_z)
1687       real stot(ngas_tot)
1688
1689       if(imap.eq.0)then    ! map cbox --> stot (both molec/cc)
1690       stot(indx(iisop_z))       = cbox(iisop_z)
1691       stot(indx(iisoprd_z))     = cbox(iisoprd_z)
1692       stot(indx(iisopp_z))      = cbox(iisopp_z)
1693       stot(indx(iisopn_z))      = cbox(iisopn_z)
1694       stot(indx(iisopo2_z))     = cbox(iisopo2_z)
1695 !
1696       else                 ! map stot --> cbox (both molec/cc)
1697       cbox(iisop_z)     = stot(indx(iisop_z))
1698       cbox(iisoprd_z)   = stot(indx(iisoprd_z))
1699       cbox(iisopp_z)    = stot(indx(iisopp_z))
1700       cbox(iisopn_z)    = stot(indx(iisopn_z))
1701       cbox(iisopo2_z)   = stot(indx(iisopo2_z))
1702       endif
1703
1704       return
1705       end subroutine mapgas_m3
1706
1707
1708
1709 !***********************************************************************
1710 ! <13.> subr mapgas_m4
1711 !
1712 ! purpose: maps gas species between stot and cbox arrays for mechanism 4
1713 !
1714 ! author : Rahul A. Zaveri
1715 ! date   : December, 1998
1716 !
1717 ! ----------------------------------------------------------------------
1718
1719       subroutine mapgas_m4( cbox, stot, imap, indx )
1720
1721       use module_data_cbmz
1722       implicit none
1723
1724 !   subr arguments
1725       integer imap, indx(ngas_z)
1726       real cbox(ngas_z)
1727       real stot(ngas_tot)
1728
1729       if(imap.eq.0)then    ! map cbox --> stot (both molec/cc)
1730       stot(indx(idms_z))        = cbox(idms_z)
1731       stot(indx(imsa_z))        = cbox(imsa_z)
1732       stot(indx(idmso_z))       = cbox(idmso_z)
1733       stot(indx(idmso2_z))      = cbox(idmso2_z)
1734       stot(indx(ich3so2h_z))    = cbox(ich3so2h_z)
1735       stot(indx(ich3sch2oo_z))  = cbox(ich3sch2oo_z)
1736       stot(indx(ich3so2_z))     = cbox(ich3so2_z)
1737       stot(indx(ich3so3_z))     = cbox(ich3so3_z)
1738       stot(indx(ich3so2oo_z))   = cbox(ich3so2oo_z)
1739       stot(indx(ich3so2ch2oo_z))= cbox(ich3so2ch2oo_z)
1740       stot(indx(imtf_z))        = cbox(imtf_z)
1741 !
1742       else                 ! map stot --> cbox (both molec/cc)
1743       cbox(idms_z)         = stot(indx(idms_z))
1744       cbox(imsa_z)         = stot(indx(imsa_z))
1745       cbox(idmso_z)        = stot(indx(idmso_z))
1746       cbox(idmso2_z)       = stot(indx(idmso2_z))
1747       cbox(ich3so2h_z)     = stot(indx(ich3so2h_z))
1748       cbox(ich3sch2oo_z)  = stot(indx(ich3sch2oo_z))
1749       cbox(ich3so2_z)      = stot(indx(ich3so2_z))
1750       cbox(ich3so3_z)      = stot(indx(ich3so3_z))
1751       cbox(ich3so2oo_z)   = stot(indx(ich3so2oo_z))
1752       cbox(ich3so2ch2oo_z)= stot(indx(ich3so2ch2oo_z))
1753       cbox(imtf_z)         = stot(indx(imtf_z))
1754       endif
1755
1756       return
1757       end subroutine mapgas_m4
1758
1759
1760
1761 !***********************************************************************
1762 ! <xx.> subr check_userpar
1763 !
1764 ! purpose: called by lsodes (external)
1765 !          computes time derivatives of species concentrations ds/dt.
1766 !          by calling subr. gasrate and ode
1767 !
1768 ! author : Rahul A. Zaveri
1769 ! date   : December 1998
1770 !
1771 !---------------------------------------------------------------------------
1772
1773       subroutine check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
1774
1775       use module_data_cbmz
1776       implicit none
1777
1778 !   subr arguments
1779       integer nruserpar, niuserpar
1780       integer iuserpar(niuserpar)
1781       real ruserpar(nruserpar)
1782
1783 !   local variables
1784       integer i
1785       real dum
1786       character*80 msg
1787
1788       if (nruserpar .ne. (5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4)) then
1789           write(msg,9010) 'nruserpar', -1, nruserpar
1790           call wrf_error_fatal( msg )
1791       end if
1792
1793       if (niuserpar .ne. (ngas_z + 1)) then
1794           write(msg,9010) 'niuserpar', -1, niuserpar
1795           call wrf_error_fatal( msg )
1796       end if
1797
1798 9010  format( '*** check_userpar error -- ', a, 1x, i8, 1x, i8 )
1799
1800
1801       return
1802       end subroutine check_userpar
1803
1804
1805
1806 !***********************************************************************
1807 ! <14.> subr gasode_cbmz
1808 !
1809 ! purpose: called by lsodes (external)
1810 !          computes time derivatives of species concentrations ds/dt.
1811 !          by calling subr. gasrate and ode
1812 !
1813 ! author : Rahul A. Zaveri
1814 ! date   : December 1998
1815 !
1816 !---------------------------------------------------------------------------
1817
1818       subroutine gasode_cbmz( ntot, tt, s, sdot,   &
1819           ruserpar, nruserpar, iuserpar, niuserpar )
1820
1821       use module_data_cbmz
1822       implicit none
1823
1824 !   subr arguments
1825       integer ntot, nruserpar, niuserpar
1826       integer iuserpar(niuserpar)
1827       real tt
1828       real s(ngas_tot), sdot(ngas_tot), ruserpar(nruserpar)
1829
1830 !   local variables
1831       integer ig, ioffset, iregime, irxn
1832       integer indx(ngas_z)
1833       real h2o, ch4, oxygen, nitrogen, hydrogen
1834       real rk_m1(nrxn_m1), r_m1(nrxn_m1)
1835       real rk_m2(nrxn_m2), r_m2(nrxn_m2)
1836       real rk_m3(nrxn_m3), r_m3(nrxn_m3)
1837       real rk_m4(nrxn_m4), r_m4(nrxn_m4)
1838       real p_m1(ngas_tot), d_m1(ngas_tot)
1839       real p_m2(ngas_tot), d_m2(ngas_tot)
1840       real p_m3(ngas_tot), d_m3(ngas_tot)
1841       real p_m4(ngas_tot), d_m4(ngas_tot)
1842
1843
1844 ! test on userpar
1845       call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
1846
1847       iregime = iuserpar(1)
1848       do ig = 1, ngas_z
1849           indx(ig) = iuserpar(ig+1)
1850       end do
1851
1852       h2o =      ruserpar(1)
1853       ch4 =      ruserpar(2)
1854       oxygen   = ruserpar(3)
1855       nitrogen = ruserpar(4)
1856       hydrogen = ruserpar(5)
1857       ioffset = 5
1858       do ig = 1, nrxn_m1
1859           rk_m1(ig) = ruserpar(ioffset+ig)
1860       end do
1861       ioffset = ioffset + nrxn_m1
1862       do ig = 1, nrxn_m2
1863           rk_m2(ig) = ruserpar(ioffset+ig)
1864       end do
1865       ioffset = ioffset + nrxn_m2
1866       do ig = 1, nrxn_m3
1867           rk_m3(ig) = ruserpar(ioffset+ig)
1868       end do
1869       ioffset = ioffset + nrxn_m3
1870       do ig = 1, nrxn_m4
1871           rk_m4(ig) = ruserpar(ioffset+ig)
1872       end do
1873
1874
1875 ! initialize to zero
1876       do irxn=1,nrxn_m1
1877       r_m1(irxn) = 0.
1878       enddo
1879
1880       do irxn=1,nrxn_m2
1881       r_m2(irxn) = 0.
1882       enddo
1883
1884       do irxn=1,nrxn_m3
1885       r_m3(irxn) = 0.
1886       enddo
1887
1888       do irxn=1,nrxn_m4
1889       r_m4(irxn) = 0.
1890       enddo
1891 !
1892 !
1893 ! initialize to zero
1894       do ig=1,ngas_tot
1895       p_m1(ig) = 0.
1896       p_m2(ig) = 0.
1897       p_m3(ig) = 0.
1898       p_m4(ig) = 0.
1899 !
1900       d_m1(ig) = 0.
1901       d_m2(ig) = 0.
1902       d_m3(ig) = 0.
1903       d_m4(ig) = 0.
1904       enddo
1905
1906
1907 !   iregime=1 --> do m1
1908 !   iregime=2 --> do m1, m2
1909 !   iregime=3 --> do m1, m2, m3
1910 !   iregime=4 --> do m1, --, --, m4
1911 !   iregime=5 --> do m1, m2, --, m4
1912 !   iregime=6 --> do m1, m2, m3, m4
1913
1914       call gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2, &
1915                            h2o, ch4, oxygen, nitrogen, hydrogen )
1916
1917       if ((iregime .eq. 2) .or.   &
1918           (iregime .eq. 3) .or.   &
1919           (iregime .eq. 5) .or.   &
1920           (iregime .eq. 6))       &
1921           call gasrate_m2( indx, s, r_m2, rk_m2 )
1922
1923       if ((iregime .eq. 3) .or.   &
1924           (iregime .eq. 6))       &
1925           call gasrate_m3( indx, s, r_m3, rk_m3 )
1926
1927       if ((iregime .eq. 4) .or.   &
1928           (iregime .eq. 5) .or.   &
1929           (iregime .eq. 6))       &
1930           call gasrate_m4( indx, s, r_m4, rk_m4, oxygen )
1931
1932       call gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 )
1933
1934       if ((iregime .eq. 2) .or.   &
1935           (iregime .eq. 3) .or.   &
1936           (iregime .eq. 5) .or.   &
1937           (iregime .eq. 6))       &
1938           call gasode_m2( indx, r_m2, p_m2, d_m2 )
1939
1940       if ((iregime .eq. 3) .or.   &
1941           (iregime .eq. 6))       &
1942           call gasode_m3( indx, r_m3, p_m3, d_m3 )
1943
1944       if ((iregime .eq. 4) .or.   &
1945           (iregime .eq. 5) .or.   &
1946           (iregime .eq. 6))       &
1947           call gasode_m4( indx, r_m4, p_m4, d_m4 )
1948
1949
1950       if (iregime .eq. 1) then
1951 ! regime = 1
1952       do ig = 1, ngas_r1
1953       sdot(ig) = real( dble(p_m1(ig)) -   &
1954                        dble(d_m1(ig)) )
1955       end do
1956
1957       else if (iregime .eq. 2) then
1958 ! regime = 2
1959       do ig = 1, ngas_r2
1960       sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)) -   &
1961                        dble(d_m1(ig)+d_m2(ig)) )
1962       end do
1963
1964       else if (iregime .eq. 3) then
1965 ! regime = 3
1966       do ig = 1, ngas_r3
1967       sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)) -   &
1968                        dble(d_m1(ig)+d_m2(ig)+d_m3(ig)) )
1969       end do
1970
1971       else if (iregime .eq. 4) then
1972 ! regime = 4
1973       do ig = 1, ngas_r4
1974       sdot(ig) = real( dble(p_m1(ig)+p_m4(ig)) -   &
1975                        dble(d_m1(ig)+d_m4(ig)) )
1976       end do
1977
1978       else if (iregime .eq. 5) then
1979 ! regime = 5
1980       do ig = 1, ngas_r5
1981       sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m4(ig)) -   &
1982                        dble(d_m1(ig)+d_m2(ig)+d_m4(ig)) )
1983       end do
1984
1985       else if (iregime .eq. 6) then
1986 ! regime = 6
1987       do ig = 1, ngas_r6
1988       sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)+p_m4(ig)) -   &
1989                        dble(d_m1(ig)+d_m2(ig)+d_m3(ig)+d_m4(ig)) )
1990       end do
1991
1992       end if
1993
1994       return
1995       end subroutine gasode_cbmz
1996
1997
1998
1999 !***********************************************************************
2000 ! <15.> subr jcs
2001 !
2002 ! purpose: external dummy jacobian evaluation for LSODES (when mf=222)
2003 !
2004 !-----------------------------------------------------------------------
2005
2006       subroutine jcs( ngas, tt, s, j, ian, jan, pdj,   &
2007           ruserpar, nruserpar, iuserpar, niuserpar )
2008
2009       implicit none
2010       integer ngas, j, ian(*), jan(*), nruserpar, niuserpar
2011       integer iuserpar(niuserpar)
2012       real tt, s(*), pdj(*)
2013       real ruserpar(nruserpar)
2014
2015 ! test on userpar
2016       call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar )
2017
2018       return
2019       end subroutine jcs
2020
2021
2022
2023 !***********************************************************************
2024 ! <16.> subr gasode_m1
2025 !
2026 ! purpose: updates production and destruction rates for mechanism 1
2027 !          background troposphere
2028 !
2029 ! author : Rahul A. Zaveri
2030 ! date   : December, 1998
2031 !
2032 ! ----------------------------------------------------------------------
2033
2034       subroutine gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 )
2035
2036       use module_data_cbmz
2037       implicit none
2038
2039 !   subr arguments
2040       integer indx(ngas_z)
2041       real r_m1(nrxn_m1), p_m1(ngas_tot), d_m1(ngas_tot)
2042       real r_m2(nrxn_m2)
2043
2044
2045       p_m1(indx(ino_z))= r_m1(1)+0.11*r_m1(2)   &
2046              +r_m1(3)+r_m1(15)+r_m1(38)
2047       d_m1(indx(ino_z))= r_m1(17)+r_m1(18)+r_m1(23)   &
2048              +r_m1(33)+r_m1(37)+r_m1(57)   &
2049              +r_m1(58)   &
2050              +r_m2(34)
2051 !
2052       p_m1(indx(ino2_z))= 0.89*r_m1(2)+r_m1(4)   &
2053               +r_m1(5)+r_m1(6)+r_m1(17)   &
2054               +r_m1(18)+r_m1(25)   &
2055               +r_m1(26)+r_m1(28)   &
2056               +r_m1(33)+r_m1(36)   &
2057               +r_m1(37)+r_m1(37)   &
2058               +r_m1(38)+r_m1(40)   &
2059               +r_m1(40)+.7*r_m1(41)   &
2060               +r_m1(43)+r_m1(57)   &
2061               +r_m1(58)+r_m1(59)   &
2062               +r_m1(60)   &
2063               +r_m2(32)+r_m2(34)+r_m2(39)
2064       d_m1(indx(ino2_z))= r_m1(1)+r_m1(15)+r_m1(16)   &
2065               +r_m1(19)+r_m1(24)   &
2066               +r_m1(34)+r_m1(35)   &
2067               +r_m1(38)+r_m1(39)   &
2068               +r_m2(31)
2069 !
2070       p_m1(indx(ino3_z))= r_m1(6)+r_m1(16)+r_m1(19)   &
2071               +r_m1(27)+r_m1(43)
2072       d_m1(indx(ino3_z))= r_m1(2)+r_m1(25)+r_m1(37)   &
2073               +r_m1(38)+r_m1(39)   &
2074               +r_m1(40)+r_m1(40)   &
2075               +r_m1(41)+r_m1(52)   &
2076               +r_m1(59)+r_m1(60)   &
2077               +r_m2(4)+r_m2(39)
2078 !
2079       p_m1(indx(in2o5_z))= r_m1(39)
2080       d_m1(indx(in2o5_z))= r_m1(6)+r_m1(42)   &
2081                +r_m1(43)
2082 !
2083       p_m1(indx(ihono_z))= r_m1(23)+r_m1(35)
2084       d_m1(indx(ihono_z))= r_m1(3)+r_m1(26)
2085 !
2086       p_m1(indx(ihno3_z))= r_m1(24)+.3*r_m1(41)   &
2087                +r_m1(42)+r_m1(42)   &
2088                +r_m1(52)   &
2089                +r_m2(4)
2090       d_m1(indx(ihno3_z))= r_m1(4)+r_m1(27)
2091 !
2092       p_m1(indx(ihno4_z))= r_m1(34)
2093       d_m1(indx(ihno4_z))= r_m1(5)+r_m1(28)   &
2094                +r_m1(36)
2095 !
2096       p_m1(indx(io3_z))= r_m1(13)   &
2097                +.4*r_m2(44)
2098       d_m1(indx(io3_z))= r_m1(7)+r_m1(8)+r_m1(14)   &
2099              +r_m1(18)+r_m1(19)+r_m1(20)   &
2100              +r_m1(21)
2101 !
2102       p_m1(indx(io1d_z))= r_m1(8)
2103       d_m1(indx(io1d_z))= r_m1(10)+r_m1(11)   &
2104               +r_m1(12)
2105 !
2106       p_m1(indx(io3p_z))= r_m1(1)+0.89*r_m1(2)   &
2107               +r_m1(7)+r_m1(10)+r_m1(11)
2108       d_m1(indx(io3p_z))= r_m1(13)+r_m1(14)   &
2109               +r_m1(15)+r_m1(16)   &
2110               +r_m1(17)
2111 !
2112       p_m1(indx(ioh_z))= r_m1(3)+r_m1(4)+2*r_m1(9)   &
2113              +2*r_m1(12)+r_m1(21)   &
2114              +r_m1(33)+.7*r_m1(41)   &
2115              +r_m1(53)+r_m1(54)+.3*r_m1(55)   &
2116              +.5*r_m1(56)
2117       d_m1(indx(ioh_z))= r_m1(20)+r_m1(22)+r_m1(23)   &
2118              +r_m1(24)+r_m1(25)+r_m1(26)   &
2119              +r_m1(27)+r_m1(28)+r_m1(29)   &
2120              +r_m1(30)+r_m1(44)+r_m1(45)   &
2121              +r_m1(46)+r_m1(47)+r_m1(48)   &
2122              +r_m1(51)+r_m1(55)+r_m1(56)   &
2123              +r_m1(65)   &
2124              +r_m2(3)
2125 !
2126       p_m1(indx(iho2_z))= r_m1(5)+r_m1(20)+r_m1(22)   &
2127               +r_m1(25)+r_m1(30)   &
2128               +r_m1(36)+r_m1(44)   &
2129               +r_m1(45)+r_m1(48)   &
2130               +2*r_m1(49)+r_m1(51)   &
2131               +r_m1(52)+r_m1(53)   &
2132               +r_m1(54)+r_m1(57)   &
2133               +r_m1(58)+r_m1(59)   &
2134               +r_m1(60)+.32*r_m1(63)   &
2135               +.6*r_m1(64)+r_m1(65)   &
2136               +r_m2(2)
2137       d_m1(indx(iho2_z))= r_m1(21)+r_m1(29)   &
2138               +r_m1(31)+r_m1(31)   &
2139               +r_m1(32)+r_m1(32)   &
2140               +r_m1(33)+r_m1(34)   &
2141               +r_m1(35)+r_m1(41)   &
2142               +r_m1(61)+r_m1(62)   &
2143               +r_m2(44)
2144 !
2145       p_m1(indx(ih2o2_z))= r_m1(31)+r_m1(32)
2146       d_m1(indx(ih2o2_z))= r_m1(9)+r_m1(30)
2147 !
2148       p_m1(indx(ico_z))= r_m1(49)+r_m1(50)+r_m1(51)   &
2149              +r_m1(52)   &
2150              +r_m2(2)
2151       d_m1(indx(ico_z))= r_m1(44)
2152 !
2153       p_m1(indx(iso2_z))= 0.0
2154       d_m1(indx(iso2_z))= r_m1(45)
2155 !
2156       p_m1(indx(ih2so4_z))= r_m1(45)
2157       d_m1(indx(ih2so4_z))= 0.0
2158 !
2159       p_m1(indx(inh3_z))= 0.0
2160       d_m1(indx(inh3_z))= 0.0
2161 !
2162       p_m1(indx(ihcl_z))= 0.0
2163       d_m1(indx(ihcl_z))= 0.0
2164 !
2165       p_m1(indx(ic2h6_z))= .2*r_m1(64)
2166       d_m1(indx(ic2h6_z))= r_m1(47)
2167 !
2168       p_m1(indx(ich3o2_z))= r_m1(46)+.7*r_m1(55)   &
2169                 +r_m2(2)+r_m2(34)+r_m2(39)+r_m2(49)
2170       d_m1(indx(ich3o2_z))= r_m1(57)+r_m1(59)   &
2171                 +r_m1(61)+r_m1(63)
2172 !
2173       p_m1(indx(iethp_z))= r_m1(47)+.5*r_m1(56)
2174       d_m1(indx(iethp_z))= r_m1(58)+r_m1(60)   &
2175                +r_m1(62)+r_m1(64)
2176 !
2177       p_m1(indx(ihcho_z))= r_m1(48)+r_m1(53)   &
2178                +.3*r_m1(55)+r_m1(57)   &
2179                +r_m1(59)+.66*r_m1(63)
2180       d_m1(indx(ihcho_z))= r_m1(49)+r_m1(50)   &
2181                +r_m1(51)+r_m1(52)
2182 !
2183       p_m1(indx(ich3oh_z))= .34*r_m1(63)
2184       d_m1(indx(ich3oh_z))= r_m1(48)
2185 !
2186       p_m1(indx(ic2h5oh_z))= 0.0
2187       d_m1(indx(ic2h5oh_z))= r_m1(65)
2188 !
2189       p_m1(indx(ich3ooh_z))= r_m1(61)
2190       d_m1(indx(ich3ooh_z))= r_m1(53)+r_m1(55)
2191 !
2192       p_m1(indx(iethooh_z))= r_m1(62)
2193       d_m1(indx(iethooh_z))= r_m1(54)+r_m1(56)
2194 !
2195       p_m1(indx(iald2_z))= r_m1(54)+.5*r_m1(56)   &
2196                +r_m1(58)+r_m1(60)   &
2197                +.8*r_m1(64)+r_m1(65)
2198       d_m1(indx(iald2_z))= r_m2(2)+r_m2(3)+r_m2(4)
2199 !
2200       p_m1(indx(ihcooh_z))= 0.0
2201       d_m1(indx(ihcooh_z))= 0.0
2202 !
2203       p_m1(indx(ircooh_z))= .4*r_m2(44)
2204       d_m1(indx(ircooh_z))= 0.0
2205 !
2206       p_m1(indx(ic2o3_z))= r_m2(3)+r_m2(4)+r_m2(32)
2207       d_m1(indx(ic2o3_z))= r_m2(31)+r_m2(34)+r_m2(39)+r_m2(44)+r_m2(49)
2208 !
2209       p_m1(indx(ipan_z))= r_m2(31)
2210       d_m1(indx(ipan_z))= r_m2(32)
2211
2212       return
2213       end subroutine gasode_m1
2214
2215
2216
2217 !***********************************************************************
2218 ! <17.> subr gasode_m2
2219 !
2220 ! purpose: updates production and destruction rates for mechanism 2
2221 !          anthropogenic hydrocarbons
2222 !
2223 ! author : Rahul A. Zaveri
2224 ! date   : December, 1998
2225 !
2226 ! ----------------------------------------------------------------------
2227
2228       subroutine gasode_m2( indx, r_m2, p_m2, d_m2 )
2229
2230       use module_data_cbmz
2231       implicit none
2232
2233 !   subr arguments
2234       integer indx(ngas_z)
2235       real r_m2(nrxn_m2), p_m2(ngas_tot), d_m2(ngas_tot)
2236
2237
2238       p_m2(indx(ino_z))= 0.0
2239       d_m2(indx(ino_z))= r_m2(20)+r_m2(33)         +r_m2(35)   &
2240               +r_m2(36)+r_m2(37)
2241 !
2242       p_m2(indx(ino2_z))= .95*r_m2(20)+r_m2(30)         +.84*r_m2(33)   &
2243                        +r_m2(35)+1.5*r_m2(36)+r_m2(37)   &
2244               +r_m2(38)         +r_m2(40)+1.5*r_m2(41)+r_m2(42)   &
2245               +.5*r_m2(51)
2246       d_m2(indx(ino2_z))= r_m2(23)
2247 !
2248       p_m2(indx(ino3_z))= 0.0
2249       d_m2(indx(ino3_z))=        +r_m2(9)+r_m2(16)+r_m2(17)+r_m2(22)   &
2250               +r_m2(38)         +r_m2(40)+r_m2(41)+r_m2(42)
2251 !
2252       p_m2(indx(ihno3_z))=        +r_m2(9)+r_m2(22)
2253       d_m2(indx(ihno3_z))= 0.0
2254 !
2255       p_m2(indx(io3_z))= 0.0
2256       d_m2(indx(io3_z))= r_m2(10)+r_m2(12)+r_m2(13)+r_m2(26)
2257 !
2258       p_m2(indx(ioh_z))= .12*r_m2(10)+.33*r_m2(12)+.60*r_m2(13)   &
2259              +.08*r_m2(26)+r_m2(27)+.23*r_m2(28)
2260       d_m2(indx(ioh_z))= r_m2(1)        +r_m2(6)+r_m2(8)+r_m2(11)   &
2261              +r_m2(14)+r_m2(15)+r_m2(18)+r_m2(19)+r_m2(21)   &
2262              +r_m2(24)+r_m2(28)+r_m2(29)
2263 !
2264       p_m2(indx(iho2_z))=        +r_m2(7)+.22*r_m2(10)+r_m2(11)   &
2265              +.26*r_m2(12)+.22*r_m2(13)+r_m2(14)+r_m2(15)   &
2266              +.2*r_m2(18)+.55*r_m2(19)+.95*r_m2(20)   &
2267              +.6*r_m2(21)+2*r_m2(24)+r_m2(25)+.76*r_m2(26)   &
2268              +.9*r_m2(27)+.9*r_m2(30)+.76*r_m2(33)+.5*r_m2(36)   &
2269              +.9*r_m2(38)+.5*r_m2(41)+.54*r_m2(48)
2270       d_m2(indx(iho2_z))= r_m2(43)         +r_m2(45)+r_m2(46)+r_m2(47)
2271 !
2272       p_m2(indx(ico_z))=        +r_m2(7)+r_m2(9)+.24*r_m2(10)   &
2273              +.31*r_m2(12)+.30*r_m2(13)+2*r_m2(24)+r_m2(25)   &
2274              +.69*r_m2(26)
2275       d_m2(indx(ico_z))= 0.0
2276 !
2277       p_m2(indx(ipar_z))= 1.1*r_m2(19)
2278       d_m2(indx(ipar_z))= r_m2(1) + r_m2(53)
2279 !
2280       p_m2(indx(ich3oh_z))= .03*r_m2(12)+.04*r_m2(13)
2281       d_m2(indx(ich3oh_z))= 0.0
2282 !
2283       p_m2(indx(ihcho_z))= r_m2(10)+1.56*r_m2(11)+.57*r_m2(12)+r_m2(14)   &
2284                +r_m2(24)+.7*r_m2(26)+r_m2(35)+.5*r_m2(36)   &
2285                +r_m2(40)+.5*r_m2(41)+.7*r_m2(50)+.5*r_m2(51)
2286       d_m2(indx(ihcho_z))= 0.0
2287 !
2288       p_m2(indx(iald2_z))= .22*r_m2(11)+.47*r_m2(12)+1.03*r_m2(13)   &
2289                +r_m2(14)+1.77*r_m2(15)+.03*r_m2(26)+.3*r_m2(27)   &
2290                +.04*r_m2(28)+.3*r_m2(30)+.25*r_m2(33)+.5*r_m2(36)   &
2291                +.3*r_m2(38)+.5*r_m2(41)+.21*r_m2(48)+.5*r_m2(51)
2292       d_m2(indx(iald2_z))= 0.0
2293 !
2294       p_m2(indx(ihcooh_z))= .52*r_m2(10)+.22*r_m2(12)
2295       d_m2(indx(ihcooh_z))= 0.0
2296 !
2297       p_m2(indx(iaone_z))= .07*r_m2(13)+.23*r_m2(15)+.74*r_m2(27)   &
2298                +.74*r_m2(30)+.62*r_m2(33)+.74*r_m2(38)   &
2299                +.57*r_m2(48)+.15*r_m2(50)
2300       d_m2(indx(iaone_z))= r_m2(5)+r_m2(6)
2301 !
2302       p_m2(indx(imgly_z))= .04*r_m2(12)+.07*r_m2(13)+.8*r_m2(19)   &
2303                +.2*r_m2(26)+.19*r_m2(28)+.15*r_m2(50)
2304       d_m2(indx(imgly_z))= r_m2(7)+r_m2(8)+r_m2(9)
2305 !
2306       p_m2(indx(ieth_z))= 0.0
2307       d_m2(indx(ieth_z))= r_m2(10)+r_m2(11)
2308 !
2309       p_m2(indx(iolet_z))= 0.0
2310       d_m2(indx(iolet_z))= r_m2(12)+r_m2(14)+r_m2(16)
2311 !
2312       p_m2(indx(iolei_z))= 0.0
2313       d_m2(indx(iolei_z))= r_m2(13)+r_m2(15)+r_m2(17)
2314 !
2315       p_m2(indx(itol_z))= 0.0
2316       d_m2(indx(itol_z))= r_m2(18)
2317 !
2318       p_m2(indx(ixyl_z))= 0.0
2319       d_m2(indx(ixyl_z))= r_m2(19)
2320 !
2321       p_m2(indx(icres_z))= .12*r_m2(18)+.05*r_m2(19)
2322       d_m2(indx(icres_z))= r_m2(21)+r_m2(22)
2323 !
2324       p_m2(indx(ito2_z))= .8*r_m2(18)+.45*r_m2(19)
2325       d_m2(indx(ito2_z))= r_m2(20)
2326 !
2327       p_m2(indx(icro_z))= .4*r_m2(21)+r_m2(22)
2328       d_m2(indx(icro_z))= r_m2(23)
2329 !
2330       p_m2(indx(iopen_z))= .95*r_m2(20)+.3*r_m2(21)
2331       d_m2(indx(iopen_z))= r_m2(24)+r_m2(25)+r_m2(26)
2332 !
2333       p_m2(indx(ionit_z))= .05*r_m2(20)+r_m2(23)+.16*r_m2(33)   &
2334                +.5*r_m2(36)+.5*r_m2(41)+r_m2(46)+.5*r_m2(51)
2335       d_m2(indx(ionit_z))= r_m2(29)+r_m2(30)
2336 !
2337       p_m2(indx(ipan_z))= 0.0
2338       d_m2(indx(ipan_z))= 0.0
2339 !
2340       p_m2(indx(ircooh_z))= .09*r_m2(12)+.16*r_m2(13)
2341       d_m2(indx(ircooh_z))= 0.0
2342 !
2343       p_m2(indx(irooh_z))= r_m2(43)+r_m2(45)
2344       d_m2(indx(irooh_z))= r_m2(27)+r_m2(28)
2345 !
2346       p_m2(indx(ich3o2_z))=        +r_m2(5)+.07*r_m2(12)+.10*r_m2(13)
2347       d_m2(indx(ich3o2_z))= 0.0
2348 !
2349       p_m2(indx(iethp_z))= .06*r_m2(12)+.05*r_m2(13)+.1*r_m2(27)   &
2350                +.1*r_m2(30)+.08*r_m2(33)+.1*r_m2(38)+.06*r_m2(48)
2351       d_m2(indx(iethp_z))= 0.0
2352 !
2353       p_m2(indx(ic2o3_z))=                +r_m2(5)+r_m2(7)+r_m2(8)   &
2354                +r_m2(9)+.13*r_m2(12)+.19*r_m2(13)+r_m2(24)   &
2355                +r_m2(25)+.62*r_m2(26)         +r_m2(35)   &
2356                +r_m2(40)+.7*r_m2(50)
2357       d_m2(indx(ic2o3_z))= 0.0
2358 !
2359       p_m2(indx(iro2_z))= r_m2(1)+.03*r_m2(12)+.09*r_m2(13)+.77*r_m2(28)
2360       d_m2(indx(iro2_z))= r_m2(33)+r_m2(38)+r_m2(43)+r_m2(48)
2361 !
2362       p_m2(indx(iano2_z))= r_m2(6)+.11*r_m2(13)
2363       d_m2(indx(iano2_z))= r_m2(35)+r_m2(40)+r_m2(45)+r_m2(50)
2364 !
2365       p_m2(indx(inap_z))= r_m2(16)+r_m2(17)+r_m2(29)
2366       d_m2(indx(inap_z))= r_m2(36)+r_m2(41)+r_m2(46)+r_m2(51)
2367 !
2368       p_m2(indx(ixo2_z))= r_m2(8)+r_m2(11)+r_m2(14)+r_m2(15)   &
2369               +.08*r_m2(18)+.5*r_m2(19)+.6*r_m2(21)   &
2370               +r_m2(24)+.03*r_m2(26)+.4*r_m2(27)+.4*r_m2(30)   &
2371               +.34*r_m2(33)+.4*r_m2(38)+.24*r_m2(48)
2372       d_m2(indx(ixo2_z))= r_m2(37)+r_m2(42)+r_m2(47)+r_m2(52)
2373
2374       p_m2(indx(ixpar_z))= 1.06*r_m2(12)+2.26*r_m2(13)+r_m2(14)   &
2375               +2.23*r_m2(15)+1.98*r_m2(27)+.42*r_m2(28)   &
2376               +1.98*r_m2(30)+1.68*r_m2(33)+r_m2(36)   &
2377               +1.98*r_m2(38)+r_m2(41)+1.25*r_m2(48)+r_m2(51)
2378       d_m2(indx(ixpar_z))= r_m2(53)
2379 !
2380       return
2381       end subroutine gasode_m2
2382
2383
2384
2385 !***********************************************************************
2386 ! <18.> subr gasode_m3
2387 !
2388 ! purpose: updates production and destruction rates for mechanism 3
2389 !          isoprene
2390 !
2391 ! author : Rahul A. Zaveri
2392 ! date   : December, 1998
2393 !
2394 ! ----------------------------------------------------------------------
2395
2396       subroutine gasode_m3( indx, r_m3, p_m3, d_m3 )
2397
2398       use module_data_cbmz
2399       implicit none
2400
2401 !   subr arguments
2402       integer indx(ngas_z)
2403       real r_m3(nrxn_m3), p_m3(ngas_tot), d_m3(ngas_tot)
2404
2405
2406       p_m3(indx(ino_z))= 0.0
2407       d_m3(indx(ino_z))= r_m3(8)+r_m3(9)+r_m3(10)
2408 !
2409       p_m3(indx(ino2_z))= .91*r_m3(8)+1.2*r_m3(9)+r_m3(10)
2410       d_m3(indx(ino2_z))= 0.0
2411 !
2412       p_m3(indx(ino3_z))= 0.0
2413       d_m3(indx(ino3_z))= r_m3(3)+r_m3(7)
2414 !
2415       p_m3(indx(ihno3_z))= .07*r_m3(7)
2416       d_m3(indx(ihno3_z))= 0.0
2417 !
2418       p_m3(indx(io3_z))= 0.0
2419       d_m3(indx(io3_z))= r_m3(2)+r_m3(6)
2420 !
2421       p_m3(indx(ioh_z))= .27*r_m3(2)+.27*r_m3(6)
2422       d_m3(indx(ioh_z))= r_m3(1)+r_m3(5)
2423 !
2424       p_m3(indx(iho2_z))= .07*r_m3(2)+.33*r_m3(4)+.1*r_m3(6)+.93*r_m3(7)   &
2425               +.91*r_m3(8)+.8*r_m3(9)+r_m3(10)
2426       d_m3(indx(iho2_z))= r_m3(11)+r_m3(12)+r_m3(13)
2427 !
2428       p_m3(indx(ico_z))= .07*r_m3(2)+.33*r_m3(4)+.16*r_m3(6)+.64*r_m3(7)   &
2429              +.59*r_m3(10)
2430       d_m3(indx(ico_z))= 0.0
2431 !
2432       p_m3(indx(ipar_z))= 1.86*r_m3(7)+0.18*r_m3(8)+1.6*r_m3(9)+2*r_m3(12)   &
2433               +2*r_m3(15)
2434       d_m3(indx(ipar_z))= 0.0
2435 !
2436       p_m3(indx(ihcho_z))= .6*r_m3(2)+.2*r_m3(4)+.15*r_m3(6)+.28*r_m3(7)   &
2437                +.63*r_m3(8)+.25*r_m3(10)
2438       d_m3(indx(ihcho_z))= 0.0
2439 !
2440       p_m3(indx(iald2_z))= .15*r_m3(2)+.07*r_m3(4)+.02*r_m3(6)+.28*r_m3(7)   &
2441                +.8*r_m3(9)+.55*r_m3(10)+r_m3(15)+.5*r_m3(16)
2442       d_m3(indx(iald2_z))= 0.0
2443 !
2444       p_m3(indx(iaone_z))= .03*r_m3(4)+.09*r_m3(6)+.63*r_m3(10)+.5*r_m3(16)
2445       d_m3(indx(iaone_z))= 0.0
2446 !
2447       p_m3(indx(imgly_z))= .85*r_m3(6)+.34*r_m3(10)
2448       d_m3(indx(imgly_z))= 0.0
2449 !
2450       p_m3(indx(ionit_z))= .93*r_m3(7)+.09*r_m3(8)+.8*r_m3(9)+r_m3(12)   &
2451                +r_m3(15)
2452       d_m3(indx(ionit_z))= 0.0
2453 !
2454       p_m3(indx(ihcooh_z))= .39*r_m3(2)+.46*r_m3(6)
2455       d_m3(indx(ihcooh_z))= 0.0
2456 !
2457       p_m3(indx(irooh_z))= r_m3(11)+r_m3(13)
2458       d_m3(indx(irooh_z))= 0.0
2459 !
2460       p_m3(indx(ich3o2_z))= .7*r_m3(4)+.05*r_m3(6)
2461       d_m3(indx(ich3o2_z))= 0.0
2462 !
2463       p_m3(indx(ic2o3_z))= .2*r_m3(2)+.97*r_m3(4)+.5*r_m3(5)+.11*r_m3(6)   &
2464                +.07*r_m3(7)
2465       d_m3(indx(ic2o3_z))= 0.0
2466 !
2467       p_m3(indx(ixo2_z))= .08*r_m3(1)+.2*r_m3(2)+.2*r_m3(5)+.07*r_m3(6)   &
2468               +.93*r_m3(7)
2469       d_m3(indx(ixo2_z))= 0.0
2470 !
2471       p_m3(indx(iisop_z))= 0.0
2472       d_m3(indx(iisop_z))= r_m3(1)+r_m3(2)+r_m3(3)
2473 !
2474       p_m3(indx(iisoprd_z))= .65*r_m3(2)+.91*r_m3(8)+.2*r_m3(9)+r_m3(14)
2475       d_m3(indx(iisoprd_z))= r_m3(4)+r_m3(5)+r_m3(6)+r_m3(7)
2476 !
2477       p_m3(indx(iisopp_z))= r_m3(1)
2478       d_m3(indx(iisopp_z))= r_m3(8)+r_m3(11)+r_m3(14)
2479 !
2480       p_m3(indx(iisopn_z))= r_m3(3)
2481       d_m3(indx(iisopn_z))= r_m3(9)+r_m3(12)+r_m3(15)
2482 !
2483       p_m3(indx(iisopo2_z))= .5*r_m3(5)
2484       d_m3(indx(iisopo2_z))= r_m3(10)+r_m3(13)+r_m3(16)
2485 !
2486       return
2487       end subroutine gasode_m3
2488
2489
2490
2491 !***********************************************************************
2492 ! <19.> subr gasode_m4
2493 !
2494 ! purpose: updates production and destruction rates for mechanism 4
2495 !          dimethylsulfide
2496 !
2497 ! author : Rahul A. Zaveri
2498 ! date   : December, 1998
2499 !
2500 ! ----------------------------------------------------------------------
2501
2502       subroutine gasode_m4( indx, r_m4, p_m4, d_m4 )
2503
2504       use module_data_cbmz
2505       implicit none
2506
2507 !   subr arguments
2508       integer indx(ngas_z)
2509       real r_m4(nrxn_m4), p_m4(ngas_tot), d_m4(ngas_tot)
2510
2511
2512       p_m4(indx(ino_z))= r_m4(19)
2513       d_m4(indx(ino_z))= r_m4(5)+r_m4(11)+r_m4(26)+r_m4(30)
2514 !
2515       p_m4(indx(ino2_z))= r_m4(5)+r_m4(11)+r_m4(26)
2516       d_m4(indx(ino2_z))= r_m4(19)+r_m4(29)
2517 !
2518       p_m4(indx(ino3_z))= 0.0
2519       d_m4(indx(ino3_z))= r_m4(2)+r_m4(14)
2520 !
2521       p_m4(indx(ihono_z))= r_m4(30)
2522       d_m4(indx(ihono_z))= 0.0
2523 !
2524       p_m4(indx(ihno3_z))= r_m4(2)+r_m4(14)+r_m4(29)
2525       d_m4(indx(ihno3_z))= 0.0
2526 !
2527       p_m4(indx(io3_z))= 0.0
2528       d_m4(indx(io3_z))= r_m4(20)
2529 !
2530       p_m4(indx(io3p_z))= 0.0
2531       d_m4(indx(io3p_z))= r_m4(3)
2532 !
2533       p_m4(indx(ioh_z))= r_m4(21)
2534       d_m4(indx(ioh_z))= r_m4(1)+r_m4(4)+r_m4(9)+r_m4(10)+r_m4(16)+r_m4(23)
2535 !
2536       p_m4(indx(iho2_z))= .965*r_m4(4)+r_m4(6)+.27*r_m4(9)+r_m4(12)+r_m4(22)   &
2537               +r_m4(27)+r_m4(32)
2538       d_m4(indx(iho2_z))= r_m4(13)+r_m4(21)+r_m4(31)
2539 !
2540       p_m4(indx(ih2o2_z))= r_m4(13)
2541       d_m4(indx(ih2o2_z))= 0.0
2542 !
2543       p_m4(indx(ico_z))= r_m4(32)
2544       d_m4(indx(ico_z))= 0.0
2545 !
2546       p_m4(indx(iso2_z))= r_m4(18)
2547       d_m4(indx(iso2_z))= 0.0
2548 !
2549       p_m4(indx(ih2so4_z))= r_m4(28)
2550       d_m4(indx(ih2so4_z))= 0.0
2551 !
2552       p_m4(indx(ihcho_z))= r_m4(5)+2*r_m4(6)+r_m4(7)+r_m4(11)+2*r_m4(12)   &
2553                +r_m4(22)+r_m4(27)
2554       d_m4(indx(ihcho_z))= r_m4(32)
2555 !
2556       p_m4(indx(ich3o2_z))= r_m4(3)+.035*r_m4(4)+.73*r_m4(9)+r_m4(18)+r_m4(28)
2557       d_m4(indx(ich3o2_z))= r_m4(6)+r_m4(12)+r_m4(15)+r_m4(22)+r_m4(27)
2558 !
2559       p_m4(indx(ich3ooh_z))= r_m4(15)
2560       d_m4(indx(ich3ooh_z))= 0.0
2561 !
2562       p_m4(indx(idms_z))= 0.0
2563       d_m4(indx(idms_z))= r_m4(1)+r_m4(2)+r_m4(3)+r_m4(4)
2564 !
2565       p_m4(indx(imsa_z))= r_m4(17)+r_m4(23)+r_m4(29)+r_m4(30)+r_m4(31)+r_m4(32)
2566       d_m4(indx(imsa_z))= 0.0
2567 !
2568       p_m4(indx(idmso_z))= .965*r_m4(4)
2569       d_m4(indx(idmso_z))= r_m4(9)
2570 !
2571       p_m4(indx(idmso2_z))= .27*r_m4(9)
2572       d_m4(indx(idmso2_z))= r_m4(10)
2573 !
2574       p_m4(indx(ich3so2h_z))= .73*r_m4(9)
2575       d_m4(indx(ich3so2h_z))= r_m4(13)+r_m4(14)+r_m4(15)+r_m4(16)+r_m4(17)
2576 !
2577       p_m4(indx(ich3sch2oo_z))= r_m4(1)+r_m4(2)
2578       d_m4(indx(ich3sch2oo_z))= r_m4(5)+r_m4(6)+r_m4(7)+r_m4(8)+r_m4(8)
2579 !
2580       p_m4(indx(ich3so2_z))= r_m4(3)+.035*r_m4(4)+r_m4(5)+r_m4(6)+r_m4(7)   &
2581                  +1.85*r_m4(8)   &
2582                  +r_m4(11)+r_m4(12)+r_m4(13)+r_m4(14)+r_m4(15)   &
2583                  +r_m4(16)+r_m4(17)+r_m4(25)
2584       d_m4(indx(ich3so2_z))= r_m4(7)+r_m4(18)+r_m4(19)+r_m4(20)+r_m4(21)   &
2585                  +r_m4(22)+r_m4(23)+r_m4(24)
2586 !
2587       p_m4(indx(ich3so3_z))= r_m4(7)+r_m4(19)+r_m4(20)+r_m4(21)+r_m4(22)   &
2588                  +r_m4(26)+r_m4(27)
2589       d_m4(indx(ich3so3_z))= r_m4(17)+r_m4(28)+r_m4(29)+r_m4(30)+r_m4(31)   &
2590                  +r_m4(32)
2591 !
2592       p_m4(indx(ich3so2oo_z))= r_m4(24)
2593       d_m4(indx(ich3so2oo_z))= r_m4(25)+r_m4(26)+r_m4(27)
2594 !
2595       p_m4(indx(ich3so2ch2oo_z))= r_m4(10)
2596       d_m4(indx(ich3so2ch2oo_z))= r_m4(11)+r_m4(12)
2597 !
2598       p_m4(indx(imtf_z))= .15*r_m4(8)
2599       d_m4(indx(imtf_z))= 0.0
2600 !
2601       return
2602       end subroutine gasode_m4
2603
2604
2605
2606 !***********************************************************************
2607 ! <20.> subr gasrate_m1
2608 !
2609 ! purpose: computes reaction rates for mechanism 1
2610 !
2611 ! author : Rahul A. Zaveri
2612 ! date   : December 1998
2613 !
2614 !-------------------------------------------------------------------------
2615
2616       subroutine gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2,   &
2617                              h2o, ch4, oxygen, nitrogen, hydrogen )
2618
2619       use module_data_cbmz
2620       implicit none
2621
2622 !   subr arguments
2623       integer indx(ngas_z)
2624       real s(ngas_tot), r_m1(nrxn_m1), r_m2(nrxn_m2)
2625       real rk_m1(nrxn_m1), rk_m2(nrxn_m2)
2626       real h2o, ch4, oxygen, nitrogen, hydrogen
2627
2628       r_m1(1)  = rk_m1(1)*s(indx(ino2_z))
2629       r_m1(2)  = rk_m1(2)*s(indx(ino3_z))
2630       r_m1(3)  = rk_m1(3)*s(indx(ihono_z))
2631       r_m1(4)  = rk_m1(4)*s(indx(ihno3_z))
2632       r_m1(5)  = rk_m1(5)*s(indx(ihno4_z))
2633       r_m1(6)  = rk_m1(6)*s(indx(in2o5_z))
2634       r_m1(7)  = rk_m1(7)*s(indx(io3_z))
2635       r_m1(8)  = rk_m1(8)*s(indx(io3_z))
2636       r_m1(9)  = rk_m1(9)*s(indx(ih2o2_z))
2637       r_m1(10) = rk_m1(10)*s(indx(io1d_z))*oxygen
2638       r_m1(11) = rk_m1(11)*s(indx(io1d_z))*nitrogen
2639       r_m1(12) = rk_m1(12)*s(indx(io1d_z))*h2o
2640       r_m1(13) = rk_m1(13)*s(indx(io3p_z))*oxygen
2641       r_m1(14) = rk_m1(14)*s(indx(io3p_z))*s(indx(io3_z))
2642       r_m1(15) = rk_m1(15)*s(indx(io3p_z))*s(indx(ino2_z))
2643       r_m1(16) = rk_m1(16)*s(indx(io3p_z))*s(indx(ino2_z))
2644       r_m1(17) = rk_m1(17)*s(indx(io3p_z))*s(indx(ino_z))
2645       r_m1(18) = rk_m1(18)*s(indx(io3_z))*s(indx(ino_z))
2646       r_m1(19) = rk_m1(19)*s(indx(io3_z))*s(indx(ino2_z))
2647       r_m1(20) = rk_m1(20)*s(indx(io3_z))*s(indx(ioh_z))
2648       r_m1(21) = rk_m1(21)*s(indx(io3_z))*s(indx(iho2_z))
2649       r_m1(22) = rk_m1(22)*s(indx(ioh_z))*hydrogen
2650       r_m1(23) = rk_m1(23)*s(indx(ioh_z))*s(indx(ino_z))
2651       r_m1(24) = rk_m1(24)*s(indx(ioh_z))*s(indx(ino2_z))
2652       r_m1(25) = rk_m1(25)*s(indx(ioh_z))*s(indx(ino3_z))
2653       r_m1(26) = rk_m1(26)*s(indx(ioh_z))*s(indx(ihono_z))
2654       r_m1(27) = rk_m1(27)*s(indx(ioh_z))*s(indx(ihno3_z))
2655       r_m1(28) = rk_m1(28)*s(indx(ioh_z))*s(indx(ihno4_z))
2656       r_m1(29) = rk_m1(29)*s(indx(ioh_z))*s(indx(iho2_z))
2657       r_m1(30) = rk_m1(30)*s(indx(ioh_z))*s(indx(ih2o2_z))
2658       r_m1(31) = rk_m1(31)*s(indx(iho2_z))*s(indx(iho2_z))
2659       r_m1(32) = rk_m1(32)*s(indx(iho2_z))*s(indx(iho2_z))*h2o
2660       r_m1(33) = rk_m1(33)*s(indx(iho2_z))*s(indx(ino_z))
2661       r_m1(34) = rk_m1(34)*s(indx(iho2_z))*s(indx(ino2_z))
2662       r_m1(35) = rk_m1(35)*s(indx(iho2_z))*s(indx(ino2_z))
2663       r_m1(36) = rk_m1(36)*s(indx(ihno4_z))
2664       r_m1(37) = rk_m1(37)*s(indx(ino3_z))*s(indx(ino_z))
2665       r_m1(38) = rk_m1(38)*s(indx(ino3_z))*s(indx(ino2_z))
2666       r_m1(39) = rk_m1(39)*s(indx(ino3_z))*s(indx(ino2_z))
2667       r_m1(40) = rk_m1(40)*s(indx(ino3_z))*s(indx(ino3_z))
2668       r_m1(41) = rk_m1(41)*s(indx(ino3_z))*s(indx(iho2_z))
2669       r_m1(42) = rk_m1(42)*s(indx(in2o5_z))*h2o
2670       r_m1(43) = rk_m1(43)*s(indx(in2o5_z))
2671       r_m1(44) = rk_m1(44)*s(indx(ico_z))*s(indx(ioh_z))
2672       r_m1(45) = rk_m1(45)*s(indx(iso2_z))*s(indx(ioh_z))
2673       r_m1(46) = rk_m1(46)*s(indx(ioh_z))*ch4
2674       r_m1(47) = rk_m1(47)*s(indx(ic2h6_z))*s(indx(ioh_z))
2675       r_m1(48) = rk_m1(48)*s(indx(ich3oh_z))*s(indx(ioh_z))
2676       r_m1(49) = rk_m1(49)*s(indx(ihcho_z))
2677       r_m1(50) = rk_m1(50)*s(indx(ihcho_z))
2678       r_m1(51) = rk_m1(51)*s(indx(ihcho_z))*s(indx(ioh_z))
2679       r_m1(52) = rk_m1(52)*s(indx(ihcho_z))*s(indx(ino3_z))
2680       r_m1(53) = rk_m1(53)*s(indx(ich3ooh_z))
2681       r_m1(54) = rk_m1(54)*s(indx(iethooh_z))
2682       r_m1(55) = rk_m1(55)*s(indx(ich3ooh_z))*s(indx(ioh_z))
2683       r_m1(56) = rk_m1(56)*s(indx(iethooh_z))*s(indx(ioh_z))
2684       r_m1(57) = rk_m1(57)*s(indx(ich3o2_z))*s(indx(ino_z))
2685       r_m1(58) = rk_m1(58)*s(indx(iethp_z))*s(indx(ino_z))
2686       r_m1(59) = rk_m1(59)*s(indx(ich3o2_z))*s(indx(ino3_z))
2687       r_m1(60) = rk_m1(60)*s(indx(iethp_z))*s(indx(ino3_z))
2688       r_m1(61) = rk_m1(61)*s(indx(ich3o2_z))*s(indx(iho2_z))
2689       r_m1(62) = rk_m1(62)*s(indx(iethp_z))*s(indx(iho2_z))
2690       r_m1(63) = rk_m1(63)*s(indx(ich3o2_z))
2691       r_m1(64) = rk_m1(64)*s(indx(iethp_z))
2692       r_m1(65) = rk_m1(65)*s(indx(ic2h5oh_z))*s(indx(ioh_z))
2693
2694       r_m2(2)  = rk_m2(2)*s(indx(iald2_z))
2695       r_m2(3)  = rk_m2(3)*s(indx(iald2_z))*s(indx(ioh_z))
2696       r_m2(4)  = rk_m2(4)*s(indx(iald2_z))*s(indx(ino3_z))
2697       r_m2(31) = rk_m2(31)*s(indx(ic2o3_z))*s(indx(ino2_z))
2698       r_m2(32) = rk_m2(32)*s(indx(ipan_z))
2699       r_m2(34) = rk_m2(34)*s(indx(ic2o3_z))*s(indx(ino_z))
2700       r_m2(39) = rk_m2(39)*s(indx(ic2o3_z))*s(indx(ino3_z))
2701       r_m2(44) = rk_m2(44)*s(indx(ic2o3_z))*s(indx(iho2_z))
2702       r_m2(49) = rk_m2(49)*s(indx(ic2o3_z))
2703
2704       return
2705       end subroutine gasrate_m1
2706
2707
2708
2709 !***********************************************************************
2710 ! <21.> subr gasrate_m2
2711 !
2712 ! purpose: computes reaction rates for mechanism 2
2713 !
2714 ! author : Rahul A. Zaveri
2715 ! date   : December 1998
2716 !
2717 !-------------------------------------------------------------------------
2718 !
2719       subroutine gasrate_m2( indx, s, r_m2, rk_m2 )
2720
2721       use module_data_cbmz
2722       implicit none
2723
2724 !   subr arguments
2725       integer indx(ngas_z)
2726       real s(ngas_tot), r_m2(nrxn_m2), rk_m2(nrxn_m2)
2727
2728       r_m2(1)  = rk_m2(1)*s(indx(ipar_z))*s(indx(ioh_z))
2729
2730       r_m2(5)  = rk_m2(5)*s(indx(iaone_z))
2731       r_m2(6)  = rk_m2(6)*s(indx(iaone_z))*s(indx(ioh_z))
2732       r_m2(7)  = rk_m2(7)*s(indx(imgly_z))
2733       r_m2(8)  = rk_m2(8)*s(indx(imgly_z))*s(indx(ioh_z))
2734       r_m2(9)  = rk_m2(9)*s(indx(imgly_z))*s(indx(ino3_z))
2735       r_m2(10) = rk_m2(10)*s(indx(ieth_z))*s(indx(io3_z))
2736       r_m2(11) = rk_m2(11)*s(indx(ieth_z))*s(indx(ioh_z))
2737       r_m2(12) = rk_m2(12)*s(indx(iolet_z))*s(indx(io3_z))
2738       r_m2(13) = rk_m2(13)*s(indx(iolei_z))*s(indx(io3_z))
2739       r_m2(14) = rk_m2(14)*s(indx(iolet_z))*s(indx(ioh_z))
2740       r_m2(15) = rk_m2(15)*s(indx(iolei_z))*s(indx(ioh_z))
2741       r_m2(16) = rk_m2(16)*s(indx(iolet_z))*s(indx(ino3_z))
2742       r_m2(17) = rk_m2(17)*s(indx(iolei_z))*s(indx(ino3_z))
2743       r_m2(18) = rk_m2(18)*s(indx(itol_z))*s(indx(ioh_z))
2744       r_m2(19) = rk_m2(19)*s(indx(ixyl_z))*s(indx(ioh_z))
2745       r_m2(20) = rk_m2(20)*s(indx(ito2_z))*s(indx(ino_z))
2746       r_m2(21) = rk_m2(21)*s(indx(icres_z))*s(indx(ioh_z))
2747       r_m2(22) = rk_m2(22)*s(indx(icres_z))*s(indx(ino3_z))
2748       r_m2(23) = rk_m2(23)*s(indx(icro_z))*s(indx(ino2_z))
2749       r_m2(24) = rk_m2(24)*s(indx(iopen_z))*s(indx(ioh_z))
2750       r_m2(25) = rk_m2(25)*s(indx(iopen_z))
2751       r_m2(26) = rk_m2(26)*s(indx(iopen_z))*s(indx(io3_z))
2752       r_m2(27) = rk_m2(27)*s(indx(irooh_z))
2753       r_m2(28) = rk_m2(28)*s(indx(irooh_z))*s(indx(ioh_z))
2754       r_m2(29) = rk_m2(29)*s(indx(ionit_z))*s(indx(ioh_z))
2755       r_m2(30) = rk_m2(30)*s(indx(ionit_z))
2756
2757       r_m2(33) = rk_m2(33)*s(indx(iro2_z))*s(indx(ino_z))
2758
2759       r_m2(35) = rk_m2(35)*s(indx(iano2_z))*s(indx(ino_z))
2760       r_m2(36) = rk_m2(36)*s(indx(inap_z))*s(indx(ino_z))
2761       r_m2(37) = rk_m2(37)*s(indx(ixo2_z))*s(indx(ino_z))
2762       r_m2(38) = rk_m2(38)*s(indx(iro2_z))*s(indx(ino3_z))
2763
2764       r_m2(40) = rk_m2(40)*s(indx(iano2_z))*s(indx(ino3_z))
2765       r_m2(41) = rk_m2(41)*s(indx(inap_z))*s(indx(ino3_z))
2766       r_m2(42) = rk_m2(42)*s(indx(ixo2_z))*s(indx(ino3_z))
2767       r_m2(43) = rk_m2(43)*s(indx(iro2_z))*s(indx(iho2_z))
2768
2769       r_m2(45) = rk_m2(45)*s(indx(iano2_z))*s(indx(iho2_z))
2770       r_m2(46) = rk_m2(46)*s(indx(inap_z))*s(indx(iho2_z))
2771       r_m2(47) = rk_m2(47)*s(indx(ixo2_z))*s(indx(iho2_z))
2772       r_m2(48) = rk_m2(48)*s(indx(iro2_z))
2773
2774       r_m2(50) = rk_m2(50)*s(indx(iano2_z))
2775       r_m2(51) = rk_m2(51)*s(indx(inap_z))
2776       r_m2(52) = rk_m2(52)*s(indx(ixo2_z))
2777       r_m2(53) = rk_m2(53)*s(indx(ipar_z))*s(indx(ixpar_z))
2778 !
2779       return
2780       end subroutine gasrate_m2
2781
2782
2783
2784 !***********************************************************************
2785 ! <22.> subr gasrate_m3
2786 !
2787 ! purpose: computes reaction rates for mechanism 3
2788 !
2789 ! author : Rahul A. Zaveri
2790 ! date   : December 1998
2791 !
2792 !-------------------------------------------------------------------------
2793 !
2794       subroutine gasrate_m3( indx, s, r_m3, rk_m3 )
2795
2796       use module_data_cbmz
2797       implicit none
2798
2799 !   subr arguments
2800       integer indx(ngas_z)
2801       real s(ngas_tot), r_m3(nrxn_m3), rk_m3(nrxn_m3)
2802
2803       r_m3(1)  = rk_m3(1)*s(indx(iisop_z))*s(indx(ioh_z))
2804       r_m3(2)  = rk_m3(2)*s(indx(iisop_z))*s(indx(io3_z))
2805       r_m3(3)  = rk_m3(3)*s(indx(iisop_z))*s(indx(ino3_z))
2806       r_m3(4)  = rk_m3(4)*s(indx(iisoprd_z))
2807       r_m3(5)  = rk_m3(5)*s(indx(iisoprd_z))*s(indx(ioh_z))
2808       r_m3(6)  = rk_m3(6)*s(indx(iisoprd_z))*s(indx(io3_z))
2809       r_m3(7)  = rk_m3(7)*s(indx(iisoprd_z))*s(indx(ino3_z))
2810       r_m3(8)  = rk_m3(8)*s(indx(iisopp_z))*s(indx(ino_z))
2811       r_m3(9)  = rk_m3(9)*s(indx(iisopn_z))*s(indx(ino_z))
2812       r_m3(10) = rk_m3(10)*s(indx(iisopo2_z))*s(indx(ino_z))
2813       r_m3(11) = rk_m3(11)*s(indx(iisopp_z))*s(indx(iho2_z))
2814       r_m3(12) = rk_m3(12)*s(indx(iisopn_z))*s(indx(iho2_z))
2815       r_m3(13) = rk_m3(13)*s(indx(iisopo2_z))*s(indx(iho2_z))
2816       r_m3(14) = rk_m3(14)*s(indx(iisopp_z))
2817       r_m3(15) = rk_m3(15)*s(indx(iisopn_z))
2818       r_m3(16) = rk_m3(16)*s(indx(iisopo2_z))
2819
2820       return
2821       end subroutine gasrate_m3
2822
2823
2824
2825 !***********************************************************************
2826 ! <23.> subr gasrate_m4
2827 !
2828 ! purpose: computes reaction rates for mechanism 4
2829 !
2830 ! author : Rahul A. Zaveri
2831 ! date   : December 1998
2832 !
2833 !-------------------------------------------------------------------------
2834 !
2835       subroutine gasrate_m4( indx, s, r_m4, rk_m4, oxygen )
2836
2837       use module_data_cbmz
2838       implicit none
2839
2840 !   subr arguments
2841       integer indx(ngas_z)
2842       real s(ngas_tot), r_m4(nrxn_m4), rk_m4(nrxn_m4)
2843       real oxygen
2844
2845       r_m4(1)  = rk_m4(1)*s(indx(idms_z))*s(indx(ioh_z))
2846       r_m4(2)  = rk_m4(2)*s(indx(idms_z))*s(indx(ino3_z))
2847       r_m4(3)  = rk_m4(3)*s(indx(idms_z))*s(indx(io3p_z))
2848       r_m4(4)  = rk_m4(4)*s(indx(idms_z))*s(indx(ioh_z))
2849       r_m4(5)  = rk_m4(5)*s(indx(ich3sch2oo_z))*s(indx(ino_z))
2850       r_m4(6)  = rk_m4(6)*s(indx(ich3sch2oo_z))*s(indx(ich3o2_z))
2851       r_m4(7)  = rk_m4(7)*s(indx(ich3sch2oo_z))*s(indx(ich3so2_z))
2852       r_m4(8)  = rk_m4(8)*s(indx(ich3sch2oo_z))*s(indx(ich3sch2oo_z))
2853       r_m4(9)  = rk_m4(9)*s(indx(idmso_z))*s(indx(ioh_z))
2854       r_m4(10) = rk_m4(10)*s(indx(idmso2_z))*s(indx(ioh_z))
2855       r_m4(11) = rk_m4(11)*s(indx(ich3so2ch2oo_z))*s(indx(ino_z))
2856       r_m4(12) = rk_m4(12)*s(indx(ich3so2ch2oo_z))*s(indx(ich3o2_z))
2857       r_m4(13) = rk_m4(13)*s(indx(ich3so2h_z))*s(indx(iho2_z))
2858       r_m4(14) = rk_m4(14)*s(indx(ich3so2h_z))*s(indx(ino3_z))
2859       r_m4(15) = rk_m4(15)*s(indx(ich3so2h_z))*s(indx(ich3o2_z))
2860       r_m4(16) = rk_m4(16)*s(indx(ich3so2h_z))*s(indx(ioh_z))
2861       r_m4(17) = rk_m4(17)*s(indx(ich3so2h_z))*s(indx(ich3so3_z))
2862       r_m4(18) = rk_m4(18)*s(indx(ich3so2_z))
2863       r_m4(19) = rk_m4(19)*s(indx(ich3so2_z))*s(indx(ino2_z))
2864       r_m4(20) = rk_m4(20)*s(indx(ich3so2_z))*s(indx(io3_z))
2865       r_m4(21) = rk_m4(21)*s(indx(ich3so2_z))*s(indx(iho2_z))
2866       r_m4(22) = rk_m4(22)*s(indx(ich3so2_z))*s(indx(ich3o2_z))
2867       r_m4(23) = rk_m4(23)*s(indx(ich3so2_z))*s(indx(ioh_z))
2868       r_m4(24) = rk_m4(24)*s(indx(ich3so2_z))*oxygen
2869       r_m4(25) = rk_m4(25)*s(indx(ich3so2oo_z))
2870       r_m4(26) = rk_m4(26)*s(indx(ich3so2oo_z))*s(indx(ino_z))
2871       r_m4(27) = rk_m4(27)*s(indx(ich3so2oo_z))*s(indx(ich3o2_z))
2872       r_m4(28) = rk_m4(28)*s(indx(ich3so3_z))
2873       r_m4(29) = rk_m4(29)*s(indx(ich3so3_z))*s(indx(ino2_z))
2874       r_m4(30) = rk_m4(30)*s(indx(ich3so3_z))*s(indx(ino_z))
2875       r_m4(31) = rk_m4(31)*s(indx(ich3so3_z))*s(indx(iho2_z))
2876       r_m4(32) = rk_m4(32)*s(indx(ich3so3_z))*s(indx(ihcho_z))
2877
2878       return
2879       end subroutine gasrate_m4
2880
2881
2882
2883 !**************************************************************************
2884 ! <24.> subr loadperoxyparameters
2885 !
2886 ! purpose: loads thermal rate coefficients for peroxy-peroxy
2887 !          permutation reactions
2888 !
2889 ! author : Rahul A. Zaveri
2890 ! date   : June 1998
2891 !
2892 ! nomenclature:
2893 ! Aperox  = Pre-exponential factor (molec-cc-s)
2894 ! Bperox  = activation energy (-E/R)  (K)
2895 !
2896 !-------------------------------------------------------------------------
2897       subroutine loadperoxyparameters( Aperox, Bperox )
2898
2899       use module_data_cbmz
2900       implicit none
2901
2902 !   subr arguments
2903       real Aperox(nperox,nperox), Bperox(nperox,nperox)
2904
2905 !   local variables
2906       integer i, j
2907
2908       Aperox(jch3o2,jch3o2)   = 2.5e-13
2909       Aperox(jethp,jethp)     = 6.8e-14
2910       Aperox(jc2o3,jc2o3)     = 2.9e-12
2911       Aperox(jano2,jano2)     = 8.0e-12
2912       Aperox(jnap,jnap)       = 1.0e-12
2913       Aperox(jro2,jro2)       = 5.3e-16
2914       Aperox(jisopp,jisopp)   = 3.1e-14
2915       Aperox(jisopn,jisopn)   = 3.1e-14
2916       Aperox(jisopo2,jisopo2) = 3.1e-14
2917       Aperox(jxo2,jxo2)       = 3.1e-14
2918
2919       Bperox(jch3o2,jch3o2)   = 190.
2920       Bperox(jethp,jethp)     = 0.0
2921       Bperox(jc2o3,jc2o3)     = 500.
2922       Bperox(jano2,jano2)     = 0.0
2923       Bperox(jnap,jnap)       = 0.0
2924       Bperox(jro2,jro2)       = 1980.
2925       Bperox(jisopp,jisopp)   = 1000.
2926       Bperox(jisopn,jisopn)   = 1000.
2927       Bperox(jisopo2,jisopo2) = 1000.
2928       Bperox(jxo2,jxo2)       = 1000.
2929
2930       do i = 1, nperox
2931       do j = 1, nperox
2932       if(i.ne.j)then
2933       Aperox(i,j) = 2.0*sqrt(Aperox(i,i)*Aperox(j,j))
2934       Bperox(i,j) = 0.5*(Bperox(i,i) + Bperox(j,j))
2935       endif
2936       enddo
2937       enddo
2938
2939 ! except for
2940       Aperox(jc2o3,jch3o2) = 1.3e-12
2941       Aperox(jch3o2,jc2o3) = 1.3e-12
2942       Bperox(jc2o3,jch3o2) = 640.
2943       Bperox(jch3o2,jc2o3) = 640.
2944
2945       return
2946       end subroutine loadperoxyparameters
2947
2948
2949
2950
2951 !**************************************************************************
2952 ! <25.> subr peroxyrateconstants
2953 !
2954 ! purpose: computes parameterized thermal rate coefficients
2955 !          for the alkylperoxy radical permutation reactions
2956 !          for the entire mechanism.
2957 !
2958 ! author : Rahul A. Zaveri
2959 ! date   : June 1998
2960 !
2961 ! nomenclature:
2962 ! rk_param  = parameterized reaction rate constants (1/s)
2963 ! rk_perox  = individual permutation reaction rate constants (molec-cc-s)
2964 ! te        = ambient atmospheric temperature (K)
2965 !
2966 !-------------------------------------------------------------------------
2967       subroutine peroxyrateconstants( tempbox, cbox,   &
2968                 Aperox, Bperox, rk_param )
2969
2970       use module_data_cbmz
2971       implicit none
2972
2973 !   subr arguments
2974       real tempbox, cbox(ngas_z)
2975       real Aperox(nperox,nperox), Bperox(nperox,nperox), rk_param(nperox)
2976
2977 !   local variables
2978       integer i, j
2979       real te
2980       real sperox(nperox), rk_perox(nperox,nperox)
2981
2982
2983       te = tempbox
2984
2985       sperox(jch3o2)  = cbox(ich3o2_z)
2986       sperox(jethp)   = cbox(iethp_z)
2987       sperox(jro2)    = cbox(iro2_z)
2988       sperox(jc2o3)   = cbox(ic2o3_z)
2989       sperox(jano2)   = cbox(iano2_z)
2990       sperox(jnap)    = cbox(inap_z)
2991       sperox(jisopp)  = cbox(iisopp_z)
2992       sperox(jisopn)  = cbox(iisopn_z)
2993       sperox(jisopo2) = cbox(iisopo2_z)
2994       sperox(jxo2)    = cbox(ixo2_z)
2995
2996 !
2997 ! initialize to zero
2998       do i = 1, nperox
2999       rk_param(i) = 0.0
3000       enddo
3001
3002       do i = 1, nperox
3003       do j = 1, nperox
3004       rk_perox(i,j) = arr( Aperox(i,j), Bperox(i,j), te )
3005       rk_param(i) = rk_param(i) + rk_perox(i,j)*sperox(j)
3006       enddo
3007       enddo
3008
3009       return
3010       end subroutine peroxyrateconstants
3011
3012
3013
3014 !***********************************************************************
3015 ! <26.> subr gasthermrk_m1
3016 !
3017 ! purpose: computes thermal reaction rate coefficients for
3018 !          mechanism 1
3019 !
3020 ! author : Rahul A. Zaveri
3021 ! date   : December 1998
3022 !
3023 !-------------------------------------------------------------------------
3024
3025       subroutine gasthermrk_m1( tempbox, cair_mlc,   &
3026                  rk_photo, rk_param, rk_m1, rk_m2 )
3027
3028       use module_data_cbmz
3029       implicit none
3030
3031 !   subr arguments
3032       real tempbox, cair_mlc
3033       real rk_photo(nphoto), rk_param(nperox)
3034       real rk_m1(nrxn_m1), rk_m2(nrxn_m2)
3035 !   local variables
3036       integer i
3037       real rk0, rk2, rk3, rki, rko, rmm, rnn, te
3038 !     real arr, troe
3039
3040
3041       te = tempbox
3042
3043       rk_m1(1)  = rk_photo(jphoto_no2)
3044       rk_m1(2)  = rk_photo(jphoto_no3)
3045       rk_m1(3)  = rk_photo(jphoto_hono)
3046       rk_m1(4)  = rk_photo(jphoto_hno3)
3047       rk_m1(5)  = rk_photo(jphoto_hno4)
3048       rk_m1(6)  = rk_photo(jphoto_n2o5)
3049       rk_m1(7)  = rk_photo(jphoto_o3a)
3050       rk_m1(8)  = rk_photo(jphoto_o3b)
3051       rk_m1(9)  = rk_photo(jphoto_h2o2)
3052       rk_m1(10) = arr(3.2e-11, 70., te)
3053       rk_m1(11) = arr(1.8e-11, 110., te)
3054       rk_m1(12) = 2.2e-10
3055       rk_m1(13) = cair_mlc*6.e-34*(te/300.)**(-2.3)
3056       rk_m1(14) = arr(8.0e-12, -2060., te)
3057       rk_m1(15) = arr(6.5e-12, -120., te)
3058 !
3059       rk0 = 9.0e-32
3060       rnn = 2.0
3061       rki = 2.2e-11
3062       rmm = 0.0
3063       rk_m1(16) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3064 !
3065       rk0 = 9.0e-32
3066       rnn = 1.5
3067       rki = 3.0e-11
3068       rmm = 0.0
3069       rk_m1(17) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3070 !
3071       rk_m1(18) = arr(2.0e-12, -1400., te)
3072       rk_m1(19) = arr(1.2e-13, -2450., te)
3073       rk_m1(20) = arr(1.6e-12, -940., te)
3074       rk_m1(21) = arr(1.1e-14, -500., te)
3075       rk_m1(22) = arr(5.5e-12, -2000., te)
3076 !
3077       rk0 = 7.0e-31
3078       rnn = 2.6
3079       rki = 3.6e-11
3080       rmm = 0.1
3081       rk_m1(23) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3082 !
3083       rk0 = 2.5e-30
3084       rnn = 4.4
3085       rki = 1.6e-11
3086       rmm = 1.7
3087       rk_m1(24) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3088 !
3089       rk_m1(25) = 2.2e-11
3090       rk_m1(26) = arr(1.8e-11, -390., te)
3091             rko = 7.2e-15 * exp(785./te)
3092             rk2 = 4.1e-16 * exp(1440./te)
3093             rk3 = 1.9e-33 * exp(725./te)*cair_mlc
3094       rk_m1(27) = rko + rk3/(1.+rk3/rk2)
3095       rk_m1(28) = arr(1.3e-12, 380., te)
3096       rk_m1(29) = arr(4.8e-11, 250., te)
3097       rk_m1(30) = arr(2.9e-12, -160., te)
3098       rk_m1(31) = 2.3e-13 * exp(600./te) +       &
3099                   1.7e-33 * exp(1000./te)*cair_mlc  ! ho2 + ho2 --> h2o2
3100       rk_m1(32) = rk_m1(31)*1.4e-21*exp(2200./te)   ! ho2 + ho2 + h2o --> h2o2
3101       rk_m1(33) = arr(3.5e-12, 250., te)
3102 !
3103       rk0 = 1.8e-31
3104       rnn = 3.2
3105       rki = 4.7e-12
3106       rmm = 1.4
3107       rk_m1(34) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3108 !
3109       rk_m1(35) = 5.0e-16
3110       rk_m1(36) = rk_m1(34)*arr(4.8e26, -10900., te)
3111       rk_m1(37) = arr(1.5e-11, 170., te)
3112       rk_m1(38) = arr(4.5e-14, -1260., te)
3113 !
3114       rk0 = 2.2e-30
3115       rnn = 3.9
3116       rki = 1.5e-12
3117       rmm = 0.7
3118       rk_m1(39) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3119 !
3120       rk_m1(40) = arr(8.5e-13, -2450., te)
3121       rk_m1(41) = 3.5e-12
3122       rk_m1(42) = 2.0e-21
3123       rk_m1(43) = rk_m1(39)*arr(3.7e26, -11000., te)
3124       rk_m1(44) = 1.5e-13 * (1.+8.18e-23*te*cair_mlc) ! co + oh --> ho2
3125
3126       rk0 = 3.0e-31
3127       rnn = 3.3
3128       rki = 1.5e-12
3129       rmm = 0.0
3130       rk_m1(45) = Troe(cair_mlc,te,rk0,rnn,rki,rmm)
3131
3132       rk_m1(46) = te**.667*arr(2.8e-14, -1575., te)
3133       rk_m1(47) = te**2*arr(1.5e-17, -492., te)
3134       rk_m1(48) = arr(6.7e-12, -600., te)
3135       rk_m1(49) = rk_photo(jphoto_hchoa)        ! hcho + hv --> 2ho2 + co
3136       rk_m1(50) = rk_photo(jphoto_hchob)        ! hcho + hv --> co
3137       rk_m1(51) = 1.0e-11
3138       rk_m1(52) = arr(3.4e-13, -1900., te)
3139       rk_m1(53) = rk_photo(jphoto_ch3ooh)
3140       rk_m1(54) = rk_photo(jphoto_ethooh)
3141       rk_m1(55) = arr(3.8e-12, 200., te)
3142       rk_m1(56) = arr(3.8e-12, 200., te)
3143       rk_m1(57) = arr(3.0e-12, 280., te)
3144       rk_m1(58) = arr(2.6e-12, 365., te)
3145       rk_m1(59) = 1.1e-12
3146       rk_m1(60) = 2.5e-12
3147       rk_m1(61) = arr(3.8e-13, 800., te)
3148       rk_m1(62) = arr(7.5e-13, 700., te)
3149       rk_m1(63) = rk_param(jch3o2)
3150       rk_m1(64) = rk_param(jethp)
3151       rk_m1(65) = arr(7.0e-12, -235.,te)
3152
3153       rk_m2(2)  = rk_photo(jphoto_ald2)
3154       rk_m2(3)  = arr(5.6e-12, 270., te)
3155       rk_m2(4)  = arr(1.4e-12, -1900., te)
3156 !
3157       rk0 = 9.7e-29
3158       rnn = 5.6
3159       rki = 9.3e-12
3160       rmm = 1.5
3161       rk_m2(31) = troe(cair_mlc,te,rk0,rnn,rki,rmm)
3162 !
3163       rk_m2(32) = rk_m2(31)*arr(1.1e28, -14000., te)
3164       rk_m2(34) = arr(5.3e-12, 360., te)
3165       rk_m2(39) = 4.0e-12
3166       rk_m2(44) = arr(4.5e-13, 1000., te)
3167       rk_m2(49) = rk_param(jc2o3)
3168
3169 ! Heterogeneous reactions
3170 !      rk_m1(65) = rk_het(1)    ! O3 -->
3171 !      rk_m1(66) = rk_het(2)    ! HO2 --> 0.5H2O2
3172 !      rk_m1(67) = rk_het(3)    ! NO2 --> 0.5HONO + 0.5HNO3
3173 !      rk_m1(68) = rk_het(4)    ! N2O5 --> 2HNO3
3174 !      rk_m1(69) = rk_het(5)    ! HNO3 --> NO2
3175 !      rk_m1(70) = rk_het(6)    ! HNO3 --> NO
3176 !      rk_m1(71) = rk_het(7)    ! NO3 --> NO + O2
3177
3178 !!$Turn off this check since initializing to 0's; wig 16-Nov-2007
3179 !!$! all rate constants but be >= 0
3180 !!$      do i = 1, nrxn_m1
3181 !!$          rk_m1(i) = max( rk_m1(i), 0.0 )
3182 !!$      end do
3183 !!$      do i = 1, nrxn_m2
3184 !!$          rk_m2(i) = max( rk_m2(i), 0.0 )
3185 !!$      end do
3186
3187       return
3188       end subroutine gasthermrk_m1
3189
3190
3191
3192 !***********************************************************************
3193 ! <27.> subr gasthermrk_m2
3194 !
3195 ! purpose: computes thermal reaction rate coefficients for
3196 !          mechanism 2
3197 !
3198 ! author : Rahul A. Zaveri
3199 ! date   : December 1998
3200 !
3201 !-------------------------------------------------------------------------
3202
3203       subroutine gasthermrk_m2( tempbox, cair_mlc,   &
3204                  rk_photo, rk_param, rk_m2 )
3205
3206       use module_data_cbmz
3207       implicit none
3208
3209 !   subr arguments
3210       real tempbox, cair_mlc
3211       real rk_photo(nphoto), rk_param(nperox), rk_m2(nrxn_m2)
3212 !   local variables
3213       integer i
3214       real rk0, rki, rmm, rnn, te
3215 !     real arr, troe
3216
3217
3218       te = tempbox
3219
3220       rk_m2(1)  = 8.1e-13
3221
3222       rk_m2(5)  = rk_photo(jphoto_aone)
3223       rk_m2(6)  = te**2*arr(5.3e-18, -230., te)
3224       rk_m2(7)  = rk_photo(jphoto_mgly)
3225       rk_m2(8)  = 1.7e-11
3226       rk_m2(9)  = arr(1.4e-12, -1900., te)
3227       rk_m2(10) = arr(1.2e-14, -2630., te)
3228 !
3229       rk0 = 1.0e-28
3230       rnn = 0.8
3231       rki = 8.8e-12
3232       rmm = 0.0
3233       rk_m2(11) = troe(cair_mlc,te,rk0,rnn,rki,rmm)
3234 !
3235       rk_m2(12) = arr(4.2e-15, -1800., te)
3236       rk_m2(13) = arr(8.9e-16, -392., te)
3237       rk_m2(14) = arr(5.8e-12, 478., te)
3238       rk_m2(15) = arr(2.9e-11, 255., te)
3239       rk_m2(16) = arr(3.1e-13, -1010., te)
3240       rk_m2(17) = 2.5e-12
3241       rk_m2(18) = arr(2.1e-12, 322., te)
3242       rk_m2(19) = arr(1.7e-11, 116., te)
3243       rk_m2(20) = 8.1e-12
3244       rk_m2(21) = 4.1e-11
3245       rk_m2(22) = 2.2e-11
3246       rk_m2(23) = 1.4e-11
3247       rk_m2(24) = 3.0e-11
3248       rk_m2(25) = rk_photo(jphoto_open)
3249       rk_m2(26) = arr(5.4e-17, -500., te)
3250       rk_m2(27) = rk_photo(jphoto_rooh)
3251       rk_m2(28) = arr(3.8e-12, 200., te)
3252       rk_m2(29) = arr(1.6e-11, -540., te)
3253       rk_m2(30) = rk_photo(jphoto_onit)
3254
3255       rk_m2(33) = 4.0e-12
3256
3257       rk_m2(35) = 4.0e-12
3258       rk_m2(36) = 4.0e-12
3259       rk_m2(37) = 4.0e-12
3260       rk_m2(38) = 2.5e-12
3261
3262       rk_m2(40) = 1.2e-12
3263       rk_m2(41) = 4.0e-12
3264       rk_m2(42) = 2.5e-12
3265       rk_m2(43) = arr(1.7e-13, 1300., te)
3266
3267       rk_m2(45) = arr(1.2e-13, 1300., te)
3268       rk_m2(46) = arr(1.7e-13, 1300., te)
3269       rk_m2(47) = arr(1.7e-13, 1300., te)
3270       rk_m2(48) = rk_param(jro2)
3271
3272       rk_m2(50) = rk_param(jano2)
3273       rk_m2(51) = rk_param(jnap)
3274       rk_m2(52) = rk_param(jxo2)
3275       rk_m2(53) = 1.0e-11               ! XPAR + PAR -->
3276
3277 !!$Turn off this check since initializing to 0's; wig 16-Nov-2007
3278 !!$! all rate constants but be >= 0
3279 !!$      do i = 1, nrxn_m2
3280 !!$          rk_m2(i) = max( rk_m2(i), 0.0 )
3281 !!$      end do
3282
3283       return
3284       end subroutine gasthermrk_m2
3285
3286
3287
3288 !***********************************************************************
3289 ! <28.> subr gasthermrk_m3
3290 !
3291 ! purpose: computes thermal reaction rate coefficients for
3292 !          mechanism 3
3293 !
3294 ! author : Rahul A. Zaveri
3295 ! date   : December 1998
3296 !
3297 !-------------------------------------------------------------------------
3298
3299       subroutine gasthermrk_m3( tempbox, cair_mlc,   &
3300                  rk_photo, rk_param, rk_m3 )
3301
3302       use module_data_cbmz
3303       implicit none
3304
3305 !   subr arguments
3306       real tempbox, cair_mlc
3307       real rk_photo(nphoto), rk_param(nperox), rk_m3(nrxn_m3)
3308 !   local variables
3309       integer i
3310       real te
3311 !     real arr
3312
3313
3314       te = tempbox
3315 !
3316       rk_m3(1)  = arr(2.6e-11, 409., te)
3317       rk_m3(2)  = arr(1.2e-14, -2013., te)
3318       rk_m3(3)  = arr(3.0e-12, -446., te)
3319       rk_m3(4)  = rk_photo(jphoto_isoprd)
3320       rk_m3(5)  = 3.3e-11
3321       rk_m3(6)  = 7.0e-18
3322       rk_m3(7)  = 1.0e-15
3323       rk_m3(8)  = 4.0e-12
3324       rk_m3(9)  = 4.0e-12
3325       rk_m3(10) = 4.0e-12
3326       rk_m3(11) = arr(1.7e-13, 1300., te)
3327       rk_m3(12) = arr(1.7e-13, 1300., te)
3328       rk_m3(13) = arr(1.7e-13, 1300., te)
3329       rk_m3(14) = rk_param(jisopp)
3330       rk_m3(15) = rk_param(jisopn)
3331       rk_m3(16) = rk_param(jisopo2)
3332
3333 ! all rate constants but be >= 0
3334       do i = 1, nrxn_m3
3335           rk_m3(i) = max( rk_m3(i), 0.0 )
3336       end do
3337
3338       return
3339       end subroutine gasthermrk_m3
3340
3341
3342
3343 !***********************************************************************
3344 ! <29.> subr gasthermrk_m4
3345 !
3346 ! purpose: computes thermal reaction rate coefficients for
3347 !          mechanism 4
3348 !
3349 ! author : Rahul A. Zaveri
3350 ! date   : December 1998
3351 !
3352 !-------------------------------------------------------------------------
3353
3354       subroutine gasthermrk_m4( tempbox, cair_mlc,   &
3355                  rk_photo, rk_param, rk_m4 )
3356
3357       use module_data_cbmz
3358       implicit none
3359
3360 !   subr arguments
3361       real tempbox, cair_mlc
3362       real rk_photo(nphoto), rk_param(nperox), rk_m4(nrxn_m4)
3363 !   local variables
3364       integer i
3365       real B_abs, B_add, rk_tot, rk_tot_den, rk_tot_num, te
3366 !     real arr
3367
3368
3369       te = tempbox
3370 !
3371       rk_m4(1) = arr(9.6e-12, -234., te)        ! ch3sch3 + oh --> ch3sch2
3372       rk_m4(2) = arr(1.4e-13, 500., te)
3373       rk_m4(3) = arr(1.3e-11, 409., te)
3374
3375 ! Hynes et al. (1986)
3376       rk_tot_num =       te * exp(-234./te) +   &
3377                    8.46e-10 * exp(7230./te) +   &
3378                    2.68e-10 * exp(7810./te)
3379       rk_tot_den = 1.04e+11 * te + 88.1 * exp(7460./te)
3380       rk_tot     = rk_tot_num/rk_tot_den
3381       B_abs      = rk_m4(1)/rk_tot
3382       B_add      = 1. - B_abs
3383
3384       rk_m4(4)  = B_add*rk_tot                  ! ch3sch3 + oh --> ch3s(oh)ch3
3385       rk_m4(5)  = 8.0e-12
3386       rk_m4(6)  = 1.8e-13
3387       rk_m4(7)  = 2.5e-13
3388       rk_m4(8)  = 8.6e-14
3389       rk_m4(9)  = 5.8e-11
3390       rk_m4(10) = 1.0e-14
3391       rk_m4(11) = 5.0e-12
3392       rk_m4(12) = 1.8e-13
3393       rk_m4(13) = 1.0e-15
3394       rk_m4(14) = 1.0e-13
3395       rk_m4(15) = 1.0e-15
3396       rk_m4(16) = 1.6e-11
3397       rk_m4(17) = 1.0e-13
3398 !Bug fixing, the coefficient changed from 2.5e-13 to 2.5e13 by Qing.Yang@pnl.gov
3399       rk_m4(18) = arr(2.5e13, -8686., te)
3400       rk_m4(19) = 1.0e-14
3401       rk_m4(20) = 5.0e-15
3402       rk_m4(21) = 2.5e-13
3403       rk_m4(22) = 2.5e-13
3404       rk_m4(23) = 5.0e-11
3405       rk_m4(24) = 2.6e-18
3406       rk_m4(25) = 3.3
3407       rk_m4(26) = 1.0e-11
3408       rk_m4(27) = 5.5e-12
3409       rk_m4(28) = arr(2.0e17, -12626., te)
3410       rk_m4(29) = 3.0e-15
3411       rk_m4(30) = 3.0e-15
3412       rk_m4(31) = 5.0e-11
3413       rk_m4(32) = 1.6e-15
3414
3415 ! all rate constants but be >= 0
3416       do i = 1, nrxn_m4
3417           rk_m4(i) = max( rk_m4(i), 0.0 )
3418       end do
3419
3420       return
3421       end subroutine gasthermrk_m4
3422
3423
3424
3425
3426 !***********************************************************************
3427 ! <26.> subr hetrateconstants
3428 !
3429 ! purpose: computes heterogeneous reaction rate coefficients
3430 !
3431 ! author : Rahul A. Zaveri
3432 ! date   : May 2000
3433 !
3434 !-------------------------------------------------------------------------
3435
3436       subroutine hetrateconstants
3437       implicit none
3438
3439       return
3440       end subroutine hetrateconstants
3441
3442
3443
3444 !***********************************************************************
3445 ! <31.> func troe
3446 !
3447 ! purpose: calculates Troe reaction rate coefficient
3448 !
3449 ! author : Rahul A. Zaveri
3450 ! date   : December 1998
3451 !-----------------------------------------------------------------------
3452
3453       real function troe( cairmlc, te, rk0, rnn, rki, rmm )
3454       implicit none
3455 !   func parameters
3456       real cairmlc, te, rk0, rnn, rki, rmm
3457 !   local variables
3458       real expo
3459
3460       rk0 = rk0*cairmlc*(te/300.)**(-rnn)
3461       rki = rki*(te/300.)**(-rmm)
3462       expo= 1./(1. + (ALOG10(rk0/rki))**2)
3463       troe  = (rk0*rki/(rk0+rki))*.6**expo
3464       return
3465       end function troe
3466
3467
3468
3469 !***********************************************************************
3470 ! <32.> func arr
3471 !
3472 ! purpose: calculates arrhenius rate coefficient
3473 !
3474 ! author : Rahul A. Zaveri
3475 ! date   : December 1998
3476 !-----------------------------------------------------------------------
3477
3478       real function arr( aa, bb, te )
3479       implicit none
3480 !   func parameters
3481       real aa, bb, te
3482
3483       arr = aa*exp(bb/te)
3484       return
3485       end function arr
3486
3487
3488
3489 !***********************************************************************
3490 ! <xx.> subr mapgas_tofrom_host
3491 !
3492 ! purpose: maps gas species between cboxold/new and host arrays
3493 !
3494 ! author : R. C. Easter
3495 ! date   : November, 2003
3496 !
3497 ! ----------------------------------------------------------------------
3498
3499         subroutine mapgas_tofrom_host( imap,          &
3500                 i_boxtest_units_convert,              &
3501                 it,jt,kt, ims,ime, jms,jme, kms,kme,  &
3502                 num_moist, num_chem, moist, chem,     &
3503                 t_phy, p_phy, rho_phy,                &
3504                 cbox, tempbox, pressbox, airdenbox,   &
3505                 cair_mlc,                             &
3506                 h2o, ch4, oxygen, nitrogen, hydrogen  )
3507
3508         use module_configure, only:                             &
3509                 p_qv,                                           &
3510                 p_so2, p_sulf, p_no2, p_no, p_o3,               &
3511                 p_hno3, p_h2o2, p_ald, p_hcho, p_op1,           &
3512                 p_op2, p_paa, p_ora1, p_ora2, p_nh3,            &
3513                 p_n2o5, p_no3, p_pan, p_hc3, p_hc5,             &
3514                 p_hc8, p_eth, p_co, p_ol2, p_olt,               &
3515                 p_oli, p_tol, p_xyl, p_aco3, p_tpan,            &
3516                 p_hono, p_hno4, p_ket, p_gly, p_mgly,           &
3517                 p_dcb, p_onit, p_csl, p_iso, p_ho,              &
3518                 p_ho2,                                          &
3519                 p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh,      &
3520                 p_par, p_to2, p_cro, p_open, p_op3,             &
3521                 p_c2o3, p_ro2, p_ano2, p_nap, p_xo2,            &
3522                 p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2,   &
3523                 p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h,       &
3524                 p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, &
3525                 p_mtf
3526         use module_data_cbmz
3527         implicit none
3528
3529 !   subr arguments
3530         INTEGER, INTENT(IN) :: imap, it,jt,kt, ims,ime, jms,jme, kms,kme, &
3531                 num_moist, num_chem, i_boxtest_units_convert
3532         REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
3533             INTENT(IN) :: moist
3534         REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
3535             INTENT(INOUT) :: chem
3536         REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3537             INTENT(IN) :: t_phy, &      ! temperature
3538                           p_phy, &      ! air pressure (Pa)
3539                           rho_phy       ! air density (kg/m3)
3540         REAL, INTENT(INOUT) :: cbox(ngas_z)
3541         REAL, INTENT(INOUT) :: tempbox, pressbox, airdenbox
3542         REAL, INTENT(INOUT) :: cair_mlc
3543         REAL, INTENT(INOUT) :: h2o, ch4, oxygen, nitrogen, hydrogen
3544
3545 !   local variables
3546         integer l
3547         real factoraa
3548         real, parameter :: eps=0.622
3549
3550
3551         tempbox = t_phy(it,kt,jt)
3552 !   p_phy = (Pa);  pressbox = (dynes/cm2)
3553         pressbox = p_phy(it,kt,jt)*10.0
3554 !   rho_phy = (kg_air/m3);  airdenbox = (mole_air/cm3)
3555         airdenbox = rho_phy(it,kt,jt)/28.966e3
3556         if (i_boxtest_units_convert .eq. 10) then
3557             airdenbox = rho_phy(it,kt,jt)
3558         end if
3559
3560         if (imap .gt. 0) goto 2000
3561
3562 !
3563 !   imap==0 -- initial species mapping from host array to cboxold
3564 !              chem --> czz --> cbox
3565 !
3566 !   note:  do not map nh3, hcl
3567 !
3568         cbox(:) = 0.0
3569
3570 !   cair_mlc = (molecules_air/cm3)
3571         cair_mlc = airdenbox*avognumkpp
3572
3573 !  moist = (kg_h2o/kg_air);  czz = (mole_h2o/cm3);  h2o = (molecules_h2o/cm3)
3574 !       czz(ih2o_z) = (moist(it,kt,jt,p_qv)/eps)*airdenbox
3575 !       if (i_boxtest_units_convert .eq. 10) then
3576 !           czz(ih2o_z) = moist(it,kt,jt,p_qv)*airdenbox
3577 !       end if
3578 !       h2o      = czz(ih2o_z)*avognumkpp
3579         h2o = (moist(it,kt,jt,p_qv)/eps)*airdenbox
3580         if (i_boxtest_units_convert .eq. 10) then
3581             h2o = moist(it,kt,jt,p_qv)*airdenbox
3582         end if
3583         h2o = h2o*avognumkpp
3584
3585 !       czz(ich4_z) = 1.7e-6*airdenbox           ! ch4 conc. in mol/cc
3586 !       ch4      = czz(ich4_z)*avognumkpp        ! ch4 conc. in molec/cc
3587         ch4      = 1.7e-6*airdenbox*avognumkpp   ! ch4 conc. in molec/cc
3588
3589         oxygen   = 0.21*cair_mlc                 ! o2 conc. in molec/cc
3590         nitrogen = 0.79*cair_mlc                 ! n2 conc. in molec/cc
3591         hydrogen = 0.58e-6*cair_mlc              ! h2 conc. in molec/cc
3592
3593 !   chem units = (ppm);  czz units = (mole/cm3);  cbox units = (molecules/cm3)
3594         factoraa = airdenbox*1.0e-6
3595         if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox
3596         factoraa = factoraa*avognumkpp
3597
3598         cbox(iso2_z)         = chem(it,kt,jt,p_so2)*factoraa
3599         cbox(ih2so4_z)       = chem(it,kt,jt,p_sulf)*factoraa
3600         cbox(ino2_z)         = chem(it,kt,jt,p_no2)*factoraa
3601         cbox(ino_z)          = chem(it,kt,jt,p_no)*factoraa
3602         cbox(io3_z)          = chem(it,kt,jt,p_o3)*factoraa
3603         cbox(ihno3_z)        = chem(it,kt,jt,p_hno3)*factoraa
3604         cbox(ih2o2_z)        = chem(it,kt,jt,p_h2o2)*factoraa
3605         cbox(iald2_z)        = chem(it,kt,jt,p_ald)*factoraa
3606         cbox(ihcho_z)        = chem(it,kt,jt,p_hcho)*factoraa
3607         cbox(ich3ooh_z)      = chem(it,kt,jt,p_op1)*factoraa
3608         cbox(iethooh_z)      = chem(it,kt,jt,p_op2)*factoraa
3609         cbox(ihcooh_z)       = chem(it,kt,jt,p_ora1)*factoraa
3610         cbox(ircooh_z)       = chem(it,kt,jt,p_ora2)*factoraa
3611         cbox(inh3_z)         = chem(it,kt,jt,p_nh3)*factoraa
3612         cbox(in2o5_z)        = chem(it,kt,jt,p_n2o5)*factoraa
3613         cbox(ino3_z)         = chem(it,kt,jt,p_no3)*factoraa
3614         cbox(ipan_z)         = chem(it,kt,jt,p_pan)*factoraa
3615         cbox(ic2h6_z)        = chem(it,kt,jt,p_eth)*factoraa
3616         cbox(ico_z)          = chem(it,kt,jt,p_co)*factoraa
3617         cbox(ieth_z)         = chem(it,kt,jt,p_ol2)*factoraa
3618         cbox(iolet_z)        = chem(it,kt,jt,p_olt)*factoraa
3619         cbox(iolei_z)        = chem(it,kt,jt,p_oli)*factoraa
3620         cbox(itol_z)         = chem(it,kt,jt,p_tol)*factoraa
3621         cbox(ixyl_z)         = chem(it,kt,jt,p_xyl)*factoraa
3622         cbox(ihono_z)        = chem(it,kt,jt,p_hono)*factoraa
3623         cbox(ihno4_z)        = chem(it,kt,jt,p_hno4)*factoraa
3624         cbox(iaone_z)        = chem(it,kt,jt,p_ket)*factoraa
3625         cbox(imgly_z)        = chem(it,kt,jt,p_mgly)*factoraa
3626         cbox(ionit_z)        = chem(it,kt,jt,p_onit)*factoraa
3627         cbox(icres_z)        = chem(it,kt,jt,p_csl)*factoraa
3628         cbox(iisop_z)        = chem(it,kt,jt,p_iso)*factoraa
3629         cbox(ioh_z)          = chem(it,kt,jt,p_ho)*factoraa
3630         cbox(iho2_z)         = chem(it,kt,jt,p_ho2)*factoraa
3631
3632         cbox(ihcl_z)         = chem(it,kt,jt,p_hcl)*factoraa
3633         cbox(ich3o2_z)       = chem(it,kt,jt,p_ch3o2)*factoraa
3634         cbox(iethp_z)        = chem(it,kt,jt,p_ethp)*factoraa
3635         cbox(ich3oh_z)       = chem(it,kt,jt,p_ch3oh)*factoraa
3636         cbox(ic2h5oh_z)      = chem(it,kt,jt,p_c2h5oh)*factoraa
3637         cbox(ipar_z)         = chem(it,kt,jt,p_par)*factoraa
3638         cbox(ito2_z)         = chem(it,kt,jt,p_to2)*factoraa
3639         cbox(icro_z)         = chem(it,kt,jt,p_cro)*factoraa
3640         cbox(iopen_z)        = chem(it,kt,jt,p_open)*factoraa
3641         cbox(irooh_z)        = chem(it,kt,jt,p_op3)*factoraa
3642         cbox(ic2o3_z)        = chem(it,kt,jt,p_c2o3)*factoraa
3643         cbox(iro2_z)         = chem(it,kt,jt,p_ro2)*factoraa
3644         cbox(iano2_z)        = chem(it,kt,jt,p_ano2)*factoraa
3645         cbox(inap_z)         = chem(it,kt,jt,p_nap)*factoraa
3646         cbox(ixo2_z)         = chem(it,kt,jt,p_xo2)*factoraa
3647         cbox(ixpar_z)        = chem(it,kt,jt,p_xpar)*factoraa
3648         cbox(iisoprd_z)      = chem(it,kt,jt,p_isoprd)*factoraa
3649         cbox(iisopp_z)       = chem(it,kt,jt,p_isopp)*factoraa
3650         cbox(iisopn_z)       = chem(it,kt,jt,p_isopn)*factoraa
3651         cbox(iisopo2_z)      = chem(it,kt,jt,p_isopo2)*factoraa
3652         cbox(idms_z)         = chem(it,kt,jt,p_dms)*factoraa
3653         cbox(imsa_z)         = chem(it,kt,jt,p_msa)*factoraa
3654         cbox(idmso_z)        = chem(it,kt,jt,p_dmso)*factoraa
3655         cbox(idmso2_z)       = chem(it,kt,jt,p_dmso2)*factoraa
3656         cbox(ich3so2h_z)     = chem(it,kt,jt,p_ch3so2h)*factoraa
3657         cbox(ich3sch2oo_z)   = chem(it,kt,jt,p_ch3sch2oo)*factoraa
3658         cbox(ich3so2_z)      = chem(it,kt,jt,p_ch3so2)*factoraa
3659         cbox(ich3so3_z)      = chem(it,kt,jt,p_ch3so3)*factoraa
3660         cbox(ich3so2oo_z)    = chem(it,kt,jt,p_ch3so2oo)*factoraa
3661         cbox(ich3so2ch2oo_z) = chem(it,kt,jt,p_ch3so2ch2oo)*factoraa
3662         cbox(imtf_z)         = chem(it,kt,jt,p_mtf)*factoraa
3663
3664         cbox(ih2o_z)         = h2o
3665         cbox(ich4_z)         = ch4
3666         cbox(io2_z)          = oxygen
3667         cbox(in2_z)          = nitrogen
3668         cbox(ih2_z)          = hydrogen
3669
3670         return
3671
3672 !
3673 !   imap==1 -- final species mapping from cbox back to host array
3674 !              cbox --> czz --> chem
3675 !
3676 !   note1:  do not map nh3, hcl, ch4
3677 !
3678 2000    continue
3679 !   chem = (ppm);  czz = (mole/cm3);  cbox = (molecules/cm3)
3680         factoraa = airdenbox*1.0e-6
3681         if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox
3682         factoraa = factoraa*avognumkpp
3683
3684         chem(it,kt,jt,p_so2)         = cbox(iso2_z)/factoraa
3685         chem(it,kt,jt,p_sulf)        = cbox(ih2so4_z)/factoraa
3686         chem(it,kt,jt,p_no2)         = cbox(ino2_z)/factoraa
3687         chem(it,kt,jt,p_no)          = cbox(ino_z)/factoraa
3688         chem(it,kt,jt,p_o3)          = cbox(io3_z)/factoraa
3689         chem(it,kt,jt,p_hno3)        = cbox(ihno3_z)/factoraa
3690         chem(it,kt,jt,p_h2o2)        = cbox(ih2o2_z)/factoraa
3691         chem(it,kt,jt,p_ald)         = cbox(iald2_z)/factoraa
3692         chem(it,kt,jt,p_hcho)        = cbox(ihcho_z)/factoraa
3693         chem(it,kt,jt,p_op1)         = cbox(ich3ooh_z)/factoraa
3694         chem(it,kt,jt,p_op2)         = cbox(iethooh_z)/factoraa
3695         chem(it,kt,jt,p_ora1)        = cbox(ihcooh_z)/factoraa
3696         chem(it,kt,jt,p_ora2)        = cbox(ircooh_z)/factoraa
3697         chem(it,kt,jt,p_nh3)         = cbox(inh3_z)/factoraa
3698         chem(it,kt,jt,p_n2o5)        = cbox(in2o5_z)/factoraa
3699         chem(it,kt,jt,p_no3)         = cbox(ino3_z)/factoraa
3700         chem(it,kt,jt,p_pan)         = cbox(ipan_z)/factoraa
3701         chem(it,kt,jt,p_eth)         = cbox(ic2h6_z)/factoraa
3702         chem(it,kt,jt,p_co)          = cbox(ico_z)/factoraa
3703         chem(it,kt,jt,p_ol2)         = cbox(ieth_z)/factoraa
3704         chem(it,kt,jt,p_olt)         = cbox(iolet_z)/factoraa
3705         chem(it,kt,jt,p_oli)         = cbox(iolei_z)/factoraa
3706         chem(it,kt,jt,p_tol)         = cbox(itol_z)/factoraa
3707         chem(it,kt,jt,p_xyl)         = cbox(ixyl_z)/factoraa
3708         chem(it,kt,jt,p_hono)        = cbox(ihono_z)/factoraa
3709         chem(it,kt,jt,p_hno4)        = cbox(ihno4_z)/factoraa
3710         chem(it,kt,jt,p_ket)         = cbox(iaone_z)/factoraa
3711         chem(it,kt,jt,p_mgly)        = cbox(imgly_z)/factoraa
3712         chem(it,kt,jt,p_onit)        = cbox(ionit_z)/factoraa
3713         chem(it,kt,jt,p_csl)         = cbox(icres_z)/factoraa
3714         chem(it,kt,jt,p_iso)         = cbox(iisop_z)/factoraa
3715         chem(it,kt,jt,p_ho)          = cbox(ioh_z)/factoraa
3716         chem(it,kt,jt,p_ho2)         = cbox(iho2_z)/factoraa
3717
3718         chem(it,kt,jt,p_hcl)         = cbox(ihcl_z)/factoraa
3719         chem(it,kt,jt,p_ch3o2)       = cbox(ich3o2_z)/factoraa
3720         chem(it,kt,jt,p_ethp)        = cbox(iethp_z)/factoraa
3721         chem(it,kt,jt,p_ch3oh)       = cbox(ich3oh_z)/factoraa
3722         chem(it,kt,jt,p_c2h5oh)      = cbox(ic2h5oh_z)/factoraa
3723         chem(it,kt,jt,p_par)         = cbox(ipar_z)/factoraa
3724         chem(it,kt,jt,p_to2)         = cbox(ito2_z)/factoraa
3725         chem(it,kt,jt,p_cro)         = cbox(icro_z)/factoraa
3726         chem(it,kt,jt,p_open)        = cbox(iopen_z)/factoraa
3727         chem(it,kt,jt,p_op3)         = cbox(irooh_z)/factoraa
3728         chem(it,kt,jt,p_c2o3)        = cbox(ic2o3_z)/factoraa
3729         chem(it,kt,jt,p_ro2)         = cbox(iro2_z)/factoraa
3730         chem(it,kt,jt,p_ano2)        = cbox(iano2_z)/factoraa
3731         chem(it,kt,jt,p_nap)         = cbox(inap_z)/factoraa
3732         chem(it,kt,jt,p_xo2)         = cbox(ixo2_z)/factoraa
3733         chem(it,kt,jt,p_xpar)        = cbox(ixpar_z)/factoraa
3734         chem(it,kt,jt,p_isoprd)      = cbox(iisoprd_z)/factoraa
3735         chem(it,kt,jt,p_isopp)       = cbox(iisopp_z)/factoraa
3736         chem(it,kt,jt,p_isopn)       = cbox(iisopn_z)/factoraa
3737         chem(it,kt,jt,p_isopo2)      = cbox(iisopo2_z)/factoraa
3738         chem(it,kt,jt,p_dms)         = cbox(idms_z)/factoraa
3739         chem(it,kt,jt,p_msa)         = cbox(imsa_z)/factoraa
3740         chem(it,kt,jt,p_dmso)        = cbox(idmso_z)/factoraa
3741         chem(it,kt,jt,p_dmso2)       = cbox(idmso2_z)/factoraa
3742         chem(it,kt,jt,p_ch3so2h)     = cbox(ich3so2h_z)/factoraa
3743         chem(it,kt,jt,p_ch3sch2oo)   = cbox(ich3sch2oo_z)/factoraa
3744         chem(it,kt,jt,p_ch3so2)      = cbox(ich3so2_z)/factoraa
3745         chem(it,kt,jt,p_ch3so3)      = cbox(ich3so3_z)/factoraa
3746         chem(it,kt,jt,p_ch3so2oo)    = cbox(ich3so2oo_z)/factoraa
3747         chem(it,kt,jt,p_ch3so2ch2oo) = cbox(ich3so2ch2oo_z)/factoraa
3748         chem(it,kt,jt,p_mtf)         = cbox(imtf_z)/factoraa
3749
3750         return
3751         end subroutine mapgas_tofrom_host
3752
3753
3754
3755 !***********************************************************************
3756 ! <xx.> subr set_gaschem_allowed_regimes
3757 !
3758 ! purpose: determines which gas-phase chemistry regimes are allowed based
3759 !          on which species are active in the simulation
3760 !
3761 ! author :
3762 ! date   :
3763 !
3764 ! ----------------------------------------------------------------------
3765
3766         subroutine set_gaschem_allowed_regimes( lunerr,   &
3767                 igaschem_allowed_m1, igaschem_allowed_m2,   &
3768                 igaschem_allowed_m3, igaschem_allowed_m4 )
3769 !
3770 !   determines which gas-phase chemistry regimes are allowed based
3771 !   on which species are active in the simulation
3772 !
3773         use module_configure, only:                             &
3774                 p_qv,                                           &
3775                 p_so2, p_sulf, p_no2, p_no, p_o3,               &
3776                 p_hno3, p_h2o2, p_ald, p_hcho, p_op1,           &
3777                 p_op2, p_paa, p_ora1, p_ora2, p_nh3,            &
3778                 p_n2o5, p_no3, p_pan, p_hc3, p_hc5,             &
3779                 p_hc8, p_eth, p_co, p_ol2, p_olt,               &
3780                 p_oli, p_tol, p_xyl, p_aco3, p_tpan,            &
3781                 p_hono, p_hno4, p_ket, p_gly, p_mgly,           &
3782                 p_dcb, p_onit, p_csl, p_iso, p_ho,              &
3783                 p_ho2,                                          &
3784                 p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh,      &
3785                 p_par, p_to2, p_cro, p_open, p_op3,             &
3786                 p_c2o3, p_ro2, p_ano2, p_nap, p_xo2,            &
3787                 p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2,   &
3788                 p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h,       &
3789                 p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, &
3790                 p_mtf
3791         use module_state_description, only:  param_first_scalar
3792         use module_data_cbmz
3793         implicit none
3794
3795 !   subr arguments
3796         integer lunerr
3797         integer igaschem_allowed_m1, igaschem_allowed_m2,   &
3798                 igaschem_allowed_m3, igaschem_allowed_m4
3799
3800 !   local variables
3801         integer nactive, ndum, p1st
3802         character*80 msg
3803
3804
3805 !   index for first "active" scalar (= 2)
3806         p1st = param_first_scalar
3807
3808 !   determine if regime 1 is allowed
3809 !   (note:  p_xxx>1 if xxx is active, p_xxx=1 if inactive)
3810         if (p_qv .lt. p1st) then
3811             msg = '*** subr set_gaschem_allowed_regimes'
3812             call peg_message( lunerr, msg )
3813             msg = '*** water vapor IS NOT ACTIVE'
3814             call peg_message( lunerr, msg )
3815             call peg_error_fatal( lunerr, msg )
3816         end if
3817
3818 !   determine if regime 1 is allowed
3819 !   (note:  p_xxx>1 if xxx is active, p_xxx=1 if inactive)
3820         nactive = 0
3821         ndum = 27
3822         if (p_no          .ge. p1st) nactive = nactive + 1
3823         if (p_no2         .ge. p1st) nactive = nactive + 1
3824         if (p_no3         .ge. p1st) nactive = nactive + 1
3825         if (p_n2o5        .ge. p1st) nactive = nactive + 1
3826         if (p_hono        .ge. p1st) nactive = nactive + 1
3827         if (p_hno3        .ge. p1st) nactive = nactive + 1
3828         if (p_hno4        .ge. p1st) nactive = nactive + 1
3829         if (p_o3          .ge. p1st) nactive = nactive + 1
3830 !       if (p_o1d         .ge. p1st) nactive = nactive + 1
3831 !       if (p_o3p         .ge. p1st) nactive = nactive + 1
3832         if (p_ho          .ge. p1st) nactive = nactive + 1
3833         if (p_ho2         .ge. p1st) nactive = nactive + 1
3834         if (p_h2o2        .ge. p1st) nactive = nactive + 1
3835         if (p_co          .ge. p1st) nactive = nactive + 1
3836         if (p_so2         .ge. p1st) nactive = nactive + 1
3837         if (p_sulf        .ge. p1st) nactive = nactive + 1
3838 !       if (p_nh3         .ge. p1st) nactive = nactive + 1
3839 !       if (p_hcl         .ge. p1st) nactive = nactive + 1
3840         if (p_eth         .ge. p1st) nactive = nactive + 1
3841         if (p_ch3o2       .ge. p1st) nactive = nactive + 1
3842         if (p_ethp        .ge. p1st) nactive = nactive + 1
3843         if (p_hcho        .ge. p1st) nactive = nactive + 1
3844         if (p_ch3oh       .ge. p1st) nactive = nactive + 1
3845         if (p_c2h5oh      .ge. p1st) nactive = nactive + 1
3846         if (p_op1         .ge. p1st) nactive = nactive + 1
3847         if (p_op2         .ge. p1st) nactive = nactive + 1
3848         if (p_ald         .ge. p1st) nactive = nactive + 1
3849         if (p_ora1        .ge. p1st) nactive = nactive + 1
3850         if (p_pan         .ge. p1st) nactive = nactive + 1
3851         if (p_ora2        .ge. p1st) nactive = nactive + 1
3852         if (p_c2o3        .ge. p1st) nactive = nactive + 1
3853
3854         if (nactive .le. 0) then
3855             igaschem_allowed_m1 = 0
3856         else if (nactive .eq. ndum) then
3857             igaschem_allowed_m1 = 1
3858         else
3859             msg = '*** subr set_gaschem_allowed_regimes'
3860             call peg_message( lunerr, msg )
3861             write(msg,90200) 1, nactive, ndum
3862             call peg_message( lunerr, msg )
3863             call peg_error_fatal( lunerr, msg )
3864         end if
3865 90200   format( '    error for regime ', i1, ', nactive, nexpected = ', 2i5 )
3866
3867 !   determine if regime 2 is allowed
3868         nactive = 0
3869         ndum = 19
3870         if (p_par         .ge. p1st) nactive = nactive + 1
3871         if (p_ket         .ge. p1st) nactive = nactive + 1
3872         if (p_mgly        .ge. p1st) nactive = nactive + 1
3873         if (p_ol2         .ge. p1st) nactive = nactive + 1
3874         if (p_olt         .ge. p1st) nactive = nactive + 1
3875         if (p_oli         .ge. p1st) nactive = nactive + 1
3876         if (p_tol         .ge. p1st) nactive = nactive + 1
3877         if (p_xyl         .ge. p1st) nactive = nactive + 1
3878         if (p_csl         .ge. p1st) nactive = nactive + 1
3879         if (p_to2         .ge. p1st) nactive = nactive + 1
3880         if (p_cro         .ge. p1st) nactive = nactive + 1
3881         if (p_open        .ge. p1st) nactive = nactive + 1
3882         if (p_onit        .ge. p1st) nactive = nactive + 1
3883         if (p_op3         .ge. p1st) nactive = nactive + 1
3884         if (p_ro2         .ge. p1st) nactive = nactive + 1
3885         if (p_ano2        .ge. p1st) nactive = nactive + 1
3886         if (p_nap         .ge. p1st) nactive = nactive + 1
3887         if (p_xo2         .ge. p1st) nactive = nactive + 1
3888         if (p_xpar        .ge. p1st) nactive = nactive + 1
3889         if (nactive .le. 0) then
3890             igaschem_allowed_m2 = 0
3891         else if (nactive .eq. ndum) then
3892             igaschem_allowed_m2 = 2
3893         else
3894             msg = '*** subr set_gaschem_allowed_regimes'
3895             call peg_message( lunerr, msg )
3896             write(msg,90200) 2, nactive, ndum
3897             call peg_message( lunerr, msg )
3898             call peg_error_fatal( lunerr, msg )
3899         end if
3900
3901 !   determine if regime 3 is allowed
3902         nactive = 0
3903         ndum = 5
3904         if (p_iso         .ge. p1st) nactive = nactive + 1
3905         if (p_isoprd      .ge. p1st) nactive = nactive + 1
3906         if (p_isopp       .ge. p1st) nactive = nactive + 1
3907         if (p_isopn       .ge. p1st) nactive = nactive + 1
3908         if (p_isopo2      .ge. p1st) nactive = nactive + 1
3909         if (nactive .le. 0) then
3910             igaschem_allowed_m3 = 0
3911         else if (nactive .eq. ndum) then
3912             igaschem_allowed_m3 = 3
3913         else
3914             msg = '*** subr set_gaschem_allowed_regimes'
3915             call peg_message( lunerr, msg )
3916             write(msg,90200) 3, nactive, ndum
3917             call peg_message( lunerr, msg )
3918             call peg_error_fatal( lunerr, msg )
3919         end if
3920
3921 !   determine if regime 4 is allowed
3922         nactive = 0
3923         ndum = 11
3924         if (p_dms         .ge. p1st) nactive = nactive + 1
3925         if (p_msa         .ge. p1st) nactive = nactive + 1
3926         if (p_dmso        .ge. p1st) nactive = nactive + 1
3927         if (p_dmso2       .ge. p1st) nactive = nactive + 1
3928         if (p_ch3so2h     .ge. p1st) nactive = nactive + 1
3929         if (p_ch3sch2oo   .ge. p1st) nactive = nactive + 1
3930         if (p_ch3so2      .ge. p1st) nactive = nactive + 1
3931         if (p_ch3so3      .ge. p1st) nactive = nactive + 1
3932         if (p_ch3so2oo    .ge. p1st) nactive = nactive + 1
3933         if (p_ch3so2ch2oo .ge. p1st) nactive = nactive + 1
3934         if (p_mtf         .ge. p1st) nactive = nactive + 1
3935         if (nactive .le. 0) then
3936             igaschem_allowed_m4 = 0
3937         else if (nactive .eq. ndum) then
3938             igaschem_allowed_m4 = 4
3939         else
3940             msg = '*** subr set_gaschem_allowed_regimes'
3941             call peg_message( lunerr, msg )
3942             write(msg,90200) 4, nactive, ndum
3943             call peg_message( lunerr, msg )
3944             call peg_error_fatal( lunerr, msg )
3945         end if
3946
3947 !   regime 1 must always be allowed
3948         if (igaschem_allowed_m1 .le. 0) then
3949             msg = '*** subr set_gaschem_allowed_regimes'
3950             call peg_message( lunerr, msg )
3951             write(msg,90300) 1
3952             call peg_message( lunerr, msg )
3953             call peg_error_fatal( lunerr, msg )
3954         end if
3955 90300   format( '    regime ', i1, ' must always be allowed' )
3956
3957 !   if regime 2 is allowed, then regime 1 must be allowed
3958         if (igaschem_allowed_m2 .gt. 0) then
3959             if (igaschem_allowed_m1 .le. 0) then
3960                 msg = '*** subr set_gaschem_allowed_regimes'
3961                 call peg_message( lunerr, msg )
3962                 write(msg,90400) 2, 1
3963                 call peg_message( lunerr, msg )
3964                 call peg_error_fatal( lunerr, msg )
3965             end if
3966         end if
3967 90400   format( '    regime ', i1, ' allowed BUT regime ', i1, ' unallowed' )
3968
3969 !   if regime 3 is allowed, then regimes 1&2 must be allowed
3970         if (igaschem_allowed_m3 .gt. 0) then
3971             if (igaschem_allowed_m1 .le. 0) then
3972                 msg = '*** subr set_gaschem_allowed_regimes'
3973                 call peg_message( lunerr, msg )
3974                 write(msg,90400) 3, 1
3975                 call peg_message( lunerr, msg )
3976                 call peg_error_fatal( lunerr, msg )
3977             else if (igaschem_allowed_m2 .le. 0) then
3978                 msg = '*** subr set_gaschem_allowed_regimes'
3979                 call peg_message( lunerr, msg )
3980                 write(msg,90400) 3, 2
3981                 call peg_message( lunerr, msg )
3982                 call peg_error_fatal( lunerr, msg )
3983             end if
3984         end if
3985
3986 !   if regime 4 is allowed, then regime 1 must be allowed
3987         if (igaschem_allowed_m4 .gt. 0) then
3988             if (igaschem_allowed_m1 .le. 0) then
3989                 msg = '*** subr set_gaschem_allowed_regimes'
3990                 call peg_message( lunerr, msg )
3991                 write(msg,90400) 4, 1
3992                 call peg_message( lunerr, msg )
3993                 call peg_error_fatal( lunerr, msg )
3994             end if
3995         end if
3996
3997         return
3998         end subroutine set_gaschem_allowed_regimes
3999
4000
4001
4002 !***********************************************************************
4003 ! <xx.> subr gasphotoconstants
4004 !
4005 ! purpose: copy photolytic rate constants from host arrays to local array
4006 !
4007 !-----------------------------------------------------------------------
4008         subroutine gasphotoconstants( rk_photo,   &
4009             i_boxtest_units_convert,               &
4010             it,jt,kt, ims,ime, jms,jme, kms,kme,   &
4011             ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
4012             ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
4013             ph_ch3o2h, ph_n2o5 )
4014 !
4015 !   copies photolytic rate constants from host arrays to local arrays
4016 !   note1:   currently 8 rate constants are scaled to other rate constants
4017 !            as is done in zz06gasphotolysis.f
4018 !   note2:   currently the n2o5 rate is set to zero
4019 !
4020         use module_data_cbmz
4021         implicit none
4022
4023 !   subr arguments
4024         integer it,jt,kt, ims,ime, jms,jme, kms,kme
4025         integer i_boxtest_units_convert
4026         real rk_photo(nphoto)
4027         real, dimension( ims:ime, kms:kme, jms:jme ) :: &
4028                ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, &
4029                ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, &
4030                ph_ch3o2h, ph_n2o5
4031
4032 !   local variables
4033         real ft
4034
4035
4036         rk_photo(:) = 0.0
4037
4038 !   these from wrf/madronnich rate constants
4039         rk_photo(jphoto_no2)    = ph_no2(it,kt,jt)
4040         rk_photo(jphoto_no3)    = ph_no3o(it,kt,jt)   &
4041                                 + ph_no3o2(it,kt,jt)
4042         rk_photo(jphoto_o3a)    = ph_o33p(it,kt,jt)
4043         rk_photo(jphoto_o3b)    = ph_o31d(it,kt,jt)
4044         rk_photo(jphoto_hono)   = ph_hno2(it,kt,jt)
4045         rk_photo(jphoto_hno3)   = ph_hno3(it,kt,jt)
4046         rk_photo(jphoto_hno4)   = ph_hno4(it,kt,jt)
4047         rk_photo(jphoto_h2o2)   = ph_h2o2(it,kt,jt)
4048         rk_photo(jphoto_ch3ooh) = ph_ch3o2h(it,kt,jt)
4049         rk_photo(jphoto_hchoa)  = ph_ch2or(it,kt,jt)
4050         rk_photo(jphoto_hchob)  = ph_ch2om(it,kt,jt)
4051         rk_photo(jphoto_n2o5)   = ph_n2o5(it,kt,jt)
4052
4053 !   these scaled to other rate constants
4054         rk_photo(jphoto_ethooh) = 0.7   *rk_photo(jphoto_h2o2)
4055         rk_photo(jphoto_ald2)   = 4.6e-4*rk_photo(jphoto_no2)
4056         rk_photo(jphoto_aone)   = 7.8e-5*rk_photo(jphoto_no2)
4057         rk_photo(jphoto_mgly)   = 9.64  *rk_photo(jphoto_hchoa)
4058         rk_photo(jphoto_open)   = 9.04  *rk_photo(jphoto_hchoa)
4059         rk_photo(jphoto_rooh)   = 0.7   *rk_photo(jphoto_h2o2)
4060         rk_photo(jphoto_onit)   = 1.0e-4*rk_photo(jphoto_no2)
4061         rk_photo(jphoto_isoprd) = .025  *rk_photo(jphoto_hchob)
4062
4063 !   convert from (1/min) to (1/s)
4064 !   (except when i_boxtest_units_convert = 10 or 20)
4065         ft = 60.0
4066         if (i_boxtest_units_convert .eq. 10) ft = 1.0
4067         if (i_boxtest_units_convert .eq. 20) ft = 1.0
4068         if (ft .ne. 1.0) then
4069             rk_photo(:) = rk_photo(:)/ft
4070         end if
4071
4072
4073         return
4074         end subroutine gasphotoconstants
4075
4076
4077
4078 !-----------------------------------------------------------------------
4079         end module module_cbmz