!************************************************************************ ! This computer software was prepared by Battelle Memorial Institute, ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY ! LIABILITY FOR THE USE OF THIS SOFTWARE. ! ! Chemistry Option: CBMZ (Carbon Bond Mechanism IV - Zaveri) ! * Primary investigator: Rahul A. Zaveri ! * Co-investigator: Richard C. Easter, William I. Gustafson Jr. ! Last update: February 2009 ! ! Contacts: ! Rahul A. Zaveri, PhD Jerome D. Fast, PhD ! Senior Research Scientist Staff Scientist ! Pacific Northwest National Laboratory Pacific Northwest National Laboratory ! P.O. Box 999, MSIN K9-30 P.O. Box 999, MSIN K9-30 ! Richland, WA 99352 Richland, WA, 99352 ! Phone: (509) 372-6159 Phone: (509) 372-6116 ! Email: Rahul.Zaveri@pnl.gov Email: Jerome.Fast@pnl.gov ! ! Please report any bugs or problems to Rahul Zaveri, the primary author ! of the code, or Jerome Fast, the WRF-chem implementation team leader. ! ! Terms of Use: ! 1) Users are requested to consult the primary author prior to ! modifying the CBMZ code or incorporating it or its submodules in ! another code. This is meant to ensure that the any linkages and/or ! assumptions will not adversely affect the operation of CBMZ. ! 2) The CBMZ source code is intended for research and educational ! purposes. Users are requested to contact the primary author ! regarding the use of CBMZ code for any commercial application. ! 3) Users preparing publications resulting from the usage of CBMZ are ! requested to cite one or more of the references below (depending on ! the application) for proper acknowledgement. ! ! References: ! 1) Zaveri R.A., and L.K. Peters (1999), A new lumped structure ! photochemical mechanism for large-scale applications, J. Geophys. ! Res., 104, 30387-30415. ! 2) Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. ! Barnard, E.G., Chapman, G.A. Grell, and S.E. Peckham (2005), ! Evolution of ozone, particulates, and aerosol direct radiative ! forcing in the vicinity of Houston using a fully-coupled ! meteorology- chemistry-aerosol model, J. Geophys. Res., 111, D21305, ! doi:10.1029/2005JD006721. ! ! Additional information: ! * www.pnl.gov/atmospheric/research/wrf-chem ! ! Support: ! Funding for developing and evaluating CBMZ was provided by the U.S. ! Department of Energy under the auspices of Atmospheric Science Program ! of the Office of Biological and Environmental Research the the PNNL ! Laboratory Research and Directed Research and Development program. !************************************************************************ module module_cbmz use module_peg_util contains !*********************************************************************** ! < 1.> subr cbmz_driver ! ! purpose: serves as an interface between subr. gas_chemistry and ! the actual solver subr such as lsodes, rodas, etc. ! ! grid : fixed i,j,k (box-model) ! ! author : Rahul A. Zaveri ! date : November 1998 ! !----------------------------------------------------------------------- subroutine cbmz_driver( & id, curr_secs, ktau, dtstep, ktauc, dtstepc, & config_flags, gmt, julday, t_phy, moist, p8w, t8w, & p_phy, chem, rho_phy, dz8w, z, z_at_w, vdrog3, & ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, & ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, & ph_ch3o2h, ph_n2o5, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_configure, only: grid_config_rec_type, num_moist, num_chem, & ! p_qv, p_so2, p_ho2, p_so4aj, p_corn, p_hcl, p_mtf !jdf p_qv, p_so2, p_ho2, p_so4aj, p_corn, p_hcl, p_mtf, & p_xyl, p_tol, p_csl, p_oli, p_olt, p_ho, p_o3, p_no !jdf USE module_data_sorgam, only: ldrog USE module_data_cbmz IMPLICIT NONE !----------------------------------------------------------------------- ! subr arguments INTEGER, INTENT(IN ) :: id, julday, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN ) :: ktau, ktauc REAL(KIND=8), INTENT(IN) :: curr_secs REAL, INTENT(IN ) :: dtstep, dtstepc, gmt ! ! advected moisture variables ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & INTENT(IN ) :: moist ! ! advected chemical tracers ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! ! arrays that hold photolysis rates ! REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT ) :: & ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, & ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, & ph_ch3o2h, ph_n2o5 ! ! on input from met model ! REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: & t_phy, & ! temperature rho_phy, & ! air density (kg/m3) p_phy, & ! NOT USED z, z_at_w, & ! NOT USED dz8w, & ! NOT USED t8w, p8w ! NOT USED ! ! for interaction with aerosols (really is output) ! REAL, DIMENSION( ims:ime , kms:kme-0 , jms:jme , ldrog ) , & INTENT(INOUT ) :: & vdrog3 ! NOT USED TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags !----------------------------------------------------------------------- ! local variables integer :: idum, iok,iprog integer :: iregime integer :: igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 integer :: igas_solver, iregime_forced integer :: i_boxtest_units_convert integer :: i_print_gasode_stats integer :: i_force_dump, mode_force_dump integer :: it, jt, kt integer :: jsolver integer :: lunerr, lunout, levdbg_err, levdbg_info integer :: mgaschem real(kind=8) :: trun real :: abs_error, rel_error real :: dtchem, tchem, tstart, tstop real :: airdenbox, pressbox, tempbox real :: cair_mlc real :: h2o, ch4, oxygen, nitrogen, hydrogen real :: cboxnew(ngas_z), cboxold(ngas_z),prodrog(ldrog) !MS real :: Aperox(nperox,nperox), Bperox(nperox,nperox) real :: rk_param(nperox), rk_photo(nphoto) real :: rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4) real factorprog integer, dimension(2,6), save :: inforodas=0 integer, dimension(6), save :: iodestatus_count=0, ioderegime_count=0 #ifdef CHEM_DBG_I !rcetestb diagnostics -------------------------------------------------- print 93010, ' ' print 93010, 'rcetestb diagnostics from cbmz_driver' print 93010, 'id, chem_opt, ktau, ktauc, julday ', & id, config_flags%chem_opt, ktau, ktauc, julday print 93020, 'dtstep, dtstepc, gmt ', & dtstep, dtstepc, gmt print 93010, 'ids/e, j, k', ids, ide, jds, jde, kds, kde print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte print 93010, 'num_moist, p_qv ', num_moist, p_qv print 93010, 'num_chem, p_so2, p_ho2 ', num_chem, p_so2, p_ho2 print 93010, 'p_so4aj, p_corn, p_hcl, p_mtf', p_so4aj, p_corn, p_hcl, p_mtf 93010 format( a, 8(1x,i6) ) 93020 format( a, 8(1p,e14.6) ) !rcetestb diagnostics -------------------------------------------------- #endif ! set some control variables to their "standard for wrf-chem" values igas_solver = 1 iregime_forced = -1 mgaschem = +1 i_boxtest_units_convert = 0 i_print_gasode_stats = 1 mode_force_dump = 0 lunerr = -1 lunout = -1 levdbg_err = 0 levdbg_info = 15 abs_error = 1.0e1 ! solver absolute tolerance (molecules/cm3) rel_error = 1.0e-3 ! solver relative tolerance ! set some control variables to non-standard values for testing ! force dumps for center column, every 3rd level ! mode_force_dump = +77 ! force dumps for center column, 1st level ! mode_force_dump = +7 ! do levdbg_info output always levdbg_info = 0 ! following call is for boxwrf testing only ! it must be commented out for actual wrf applications ! call boxtest_get_extra_args( & ! igas_solver, iregime_forced, & ! i_boxtest_units_convert, lunerr, lunout, & ! abs_error, rel_error, trun ) ! currently nothing is done with vdrog3 ! vdrog3(its:ite,kts:kte,jts:jte,:) = 0.0 !This is already set to zero in chem_driver. ! initialize rate constant arrays to 0.0 since not all elements are set ! to valid values before some were used in loops, wig 16-Nov-2007 rk_m1(:) = 0. rk_m2(:) = 0. rk_m3(:) = 0. rk_m4(:) = 0. ! determine which regimes are allowed ! based on which gas species are "active" call set_gaschem_allowed_regimes( lunerr, & igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 ) ! ! main loop -- do gas chemistry at each i,j,k ! do 2900 jt = jts, jte do 2900 kt = kts, kte do 2900 it = its, ite trun = curr_secs ! run time in s tchem = mod( real(gmt*3600.0,8)+trun, 86400._8 ) ! time from 00 UTC in s dtchem = dtstepc tstart = tchem ! s tstop = tstart + dtchem ! s ! skip integration for very small dtchem if ((tstop-tstart) .le. 1.0e-5) goto 2900 ! initial species mapping from host array call mapgas_tofrom_host( 0, & i_boxtest_units_convert, & it,jt,kt, ims,ime, jms,jme, kms,kme, & num_moist, num_chem, moist, chem, & t_phy, p_phy, rho_phy, & cboxold, tempbox, pressbox, airdenbox, & cair_mlc, & h2o, ch4, oxygen, nitrogen, hydrogen ) cboxnew(:) = cboxold(:) ! determine regime call selectgasregime( iregime, iregime_forced, cboxold, & igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 ) idum = iregime if ((idum .lt. 1) .or. (idum .ge. 6)) idum = 6 ioderegime_count(idum) = ioderegime_count(idum) + 1 iodestatus_count(6) = iodestatus_count(6) + 1 ! compute rate constants ! transfer/map incoming photolysis rate contants to local array call gasphotoconstants( rk_photo, & i_boxtest_units_convert, & it,jt,kt, ims,ime, jms,jme, kms,kme, & ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, & ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, & ph_ch3o2h, ph_n2o5 ) ! loads Aperox and Bperox call loadperoxyparameters( Aperox, Bperox ) ! calculate parameterized rate constants call peroxyrateconstants( tempbox, cboxold, & Aperox, Bperox, rk_param ) ! calculate thermal rate constants call gasrateconstants( iregime, tempbox, cair_mlc, & rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 ) ! mode_force_dump selects a detailed dump of gaschem at either ! first ijk grid, first ij column, all ijk, or no ijk i_force_dump = 0 if (mode_force_dump .eq. 1) then if ((it.eq.its) .and. (jt.eq.jts) & .and. (kt.eq.kts)) i_force_dump = 1 else if (mode_force_dump .eq. 10) then if ((it.eq.its) .and. (jt.eq.jts)) i_force_dump = 1 else if (mode_force_dump .eq. 100) then i_force_dump = 1 else if (mode_force_dump .eq. 7) then if ( (it .eq. (its+ite)/2) .and. & (jt .eq. (jts+jte)/2) .and. & (kt .eq. kts) ) i_force_dump = 1 else if (mode_force_dump .eq. 77) then if ( (it .eq. (its+ite)/2) .and. & (jt .eq. (jts+jte)/2) .and. & (mod(kt-kts,3) .eq. 0) ) i_force_dump = 1 end if ! rodas iok = 0 jsolver = 0 if (igas_solver .eq. 1) then jsolver = 1 call gasodesolver_rodas( tstart, tstop, iok, & it, jt, kt, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info, & i_force_dump, inforodas, iodestatus_count, & abs_error, rel_error, trun, & tempbox, pressbox, airdenbox, cboxnew, cboxold, & rk_m1, rk_m2, rk_m3, rk_m4 ) endif ! lsodes if (igas_solver.eq.2 .or. iok.le.0) then jsolver = 2 call gasodesolver_lsodes( tstart, tstop, iok, & it, jt, kt, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info, & i_force_dump, iodestatus_count, & abs_error, rel_error, trun, & tempbox, pressbox, airdenbox, cboxnew, cboxold, & rk_m1, rk_m2, rk_m3, rk_m4 ) endif ! final species mapping back to host array -- only when iok > 0 if (iok .gt. 0) then call mapgas_tofrom_host( 1, & i_boxtest_units_convert, & it,jt,kt, ims,ime, jms,jme, kms,kme, & num_moist, num_chem, moist, chem, & t_phy, p_phy, rho_phy, & cboxnew, tempbox, pressbox, airdenbox, & cair_mlc, & h2o, ch4, oxygen, nitrogen, hydrogen ) end if ! following call is for boxwrf testing only ! it must be commented out for actual wrf applications ! call boxtest_set_extra_args( iregime, it, jt, kt ) if (iok .gt. 0) then call update_vdrog3_cbmz(rk_m2,rk_m3,dtstepc,cboxnew, cboxold, & prodrog) !prodrog units are in molecules/cm3 ! chem units = (ppm) airdenbox = rho_phy(it,kt,jt)/28.966e3 if (i_boxtest_units_convert .eq. 10) then airdenbox = rho_phy(it,kt,jt) end if factorprog = airdenbox*1.0e-6 if (i_boxtest_units_convert .eq. 10) factorprog = airdenbox factorprog = factorprog*avognumkpp ! convert prodrog from units of molecules/cm3 to ppmv do iprog=1,ldrog prodrog(iprog)=prodrog(iprog)/factorprog enddo ! Now assign the prodrog values to vdrog3 array do iprog=1,ldrog vdrog3(it,kt,jt,iprog)=prodrog(iprog) vdrog3(it,kt,jt,iprog)=max(0.,vdrog3(it,kt,jt,iprog)) enddo endif ! iok .gt. 0 !MS 2900 continue if (i_print_gasode_stats .gt. 0) & call print_gasode_stats( lunout, levdbg_info, & inforodas, iodestatus_count, ioderegime_count ) return end subroutine cbmz_driver !*********************************************************************** ! < xx.> subr update_vdrog3_cbmz ! ! purpose: updates vdrog3 concentrations for MADE-SORGAM-CBMZ package ! Written by Manish Shrivastava on 12/20/2011 !----------------------------------------------------------------------- subroutine update_vdrog3_cbmz(rk_m2,rk_m3,dtstepc,cboxnew, cboxold, & prodrog) USE module_data_sorgam, only: ldrog USE module_data_cbmz implicit none real,intent(in) :: rk_m2(nrxn_m2),cboxnew(ngas_z), cboxold(ngas_z) real,intent(in) :: rk_m3(nrxn_m3) real,intent(in) :: dtstepc real,intent(inout) :: prodrog(ldrog) integer i real temp_par_new,temp_par_old,f_api,f_lim ! f_api is fraction of isoprene attributed to API ! f_lim is fraction of isoprene attributed to LIM f_api=0.292 f_lim=0.064 ! First initialize prodrog temp_par_new=0.0 temp_par_old=0.0 do i=1,ldrog prodrog(i)=0.0 enddo temp_par_new=temp_par_new+cboxnew(ipar_z)/7.9 !since hc8 emissions are multiplied by 7.9 to calculate par, we divide by this number temp_par_old=temp_par_old+cboxold(ipar_z)/7.9 !since hc8 emissions are multiplied by 7.9 to calculate par, we divide by this number ! The order of following reactions is from module_data_sorgam.F ! prodrog(1) is rxn of XYL+OH rk_m2(19) prodrog(1)=0.5*(cboxnew(ixyl_z)+cboxold(ixyl_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(19)*dtstepc !prodrog(2) is rxn of TOL+OH rk_m2(18) prodrog(2)=0.5*(cboxnew(itol_z)+cboxold(itol_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(18)*dtstepc !prodrog(3) is rxn of CRES+OH rk_m2(21) prodrog(3)=0.5*(cboxnew(icres_z)+cboxold(icres_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(21)*dtstepc !prodrog(4) is rxn of CRES+NO3 rk_m2(22) prodrog(4)=0.5*(cboxnew(icres_z)+cboxold(icres_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) & *rk_m2(22)*dtstepc !prodrog(5) is rxn of HC+OH rk_m2(1) prodrog(5)=0.5*(temp_par_new+temp_par_old)*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(1)*dtstepc !Set to 0.0 as PAR should be replaced by HC8 or higher alkanes !prodrog(6) is rxn of OLEI+OH rk_m2(15) prodrog(6)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(15)*dtstepc !prodrog(7) is rxn of OLEI+NO3 rk_m2(17) prodrog(7)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) & *rk_m2(17)*dtstepc !prodrog(8) is rxn of OLEI+O3 rk_m2(13) prodrog(8)=0.5*(cboxnew(iolei_z)+cboxold(iolei_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) & *rk_m2(13)*dtstepc !prodrog(9) is rxn of OLET+OH rk_m2(14) prodrog(9)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m2(14)*dtstepc !prodrog(10) is rxn of OLET+NO3 rk_m2(16) prodrog(10)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) & *rk_m2(16)*dtstepc !prodrog(11) is rxn of OLET+O3 rk_m2(12) prodrog(11)=0.5*(cboxnew(iolet_z)+cboxold(iolet_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) & *rk_m2(12)*dtstepc prodrog(12)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m3(1)*dtstepc*f_api prodrog(13)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) & *rk_m3(3)*dtstepc*f_api prodrog(14)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) & *rk_m3(2)*dtstepc*f_api prodrog(15)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ioh_z)+cboxold(ioh_z)) & *rk_m3(1)*dtstepc*f_lim prodrog(16)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(ino3_z)+cboxold(ino3_z)) & *rk_m3(3)*dtstepc*f_lim prodrog(17)=0.5*(cboxnew(iisop_z)+cboxold(iisop_z))*0.5*(cboxnew(io3_z)+cboxold(io3_z)) & *rk_m3(2)*dtstepc*f_lim ! prodrog(12) to prodrog(17) are for biogenic species API and LIM ! Since these species are not available in CBMZ prodrog(12) to prodrog(17) are not calculated return end subroutine update_vdrog3_cbmz !*********************************************************************** ! < xx.> subr print_gasode_stats ! ! purpose: writes some statistics on ode solver performance to unit lunout ! !----------------------------------------------------------------------- subroutine print_gasode_stats( lunout, levdbg, & inforodas, iodestatus_count, ioderegime_count ) implicit none ! subr arguments integer lunout, levdbg integer inforodas(2,6), iodestatus_count(6), ioderegime_count(6) ! local variables integer i, j character*80 msg msg = ' ' call peg_debugmsg( lunout, levdbg, msg ) msg = 'output from dump_cbmz_gasodeinfo' call peg_debugmsg( lunout, levdbg, msg ) write(msg,9100) 'oderegime(1-6)', (ioderegime_count(i), i=1,6) call peg_debugmsg( lunout, levdbg, msg ) write(msg,9100) 'odestatus(1-6)', (iodestatus_count(i), i=1,6) call peg_debugmsg( lunout, levdbg, msg ) write(msg,9200) & 'inforodas(1-3)', ((inforodas(j,i), j=1,2), i=1,3) call peg_debugmsg( lunout, levdbg, msg ) write(msg,9200) & 'inforodas(4-6)', ((inforodas(j,i), j=1,2), i=4,6) call peg_debugmsg( lunout, levdbg, msg ) 9100 format( a, 6i11 ) 9200 format( a, 3( i11, '--', i9.9 ) ) return end subroutine print_gasode_stats !*********************************************************************** ! < xx.> subr gasodesolver_rodas ! ! purpose: interfaces to rodas ode solver ! !----------------------------------------------------------------------- subroutine gasodesolver_rodas( tstart, tstop, iok, & isvode, jsvode, ksvode, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info, & i_force_dump, inforodas, iodestatus_count, & abs_error, rel_error, trun, & tempbox, pressbox, airdenbox, cboxnew, cboxold, & rk_m1, rk_m2, rk_m3, rk_m4 ) use module_data_cbmz use module_cbmz_rodas_prep, only: & cbmz_v02r01_mapconcs, cbmz_v02r01_maprates, cbmz_v02r01_torodas, & cbmz_v02r02_mapconcs, cbmz_v02r02_maprates, cbmz_v02r02_torodas, & cbmz_v02r03_mapconcs, cbmz_v02r03_maprates, cbmz_v02r03_torodas, & cbmz_v02r04_mapconcs, cbmz_v02r04_maprates, cbmz_v02r04_torodas, & cbmz_v02r05_mapconcs, cbmz_v02r05_maprates, cbmz_v02r05_torodas, & cbmz_v02r06_mapconcs, cbmz_v02r06_maprates, cbmz_v02r06_torodas implicit none ! subr arguments real(kind=8) :: trun integer iok, isvode, jsvode, ksvode, i_force_dump, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info integer inforodas(2,6), iodestatus_count(6) real tstart, tstop, abs_error, rel_error real tempbox, pressbox, airdenbox real cboxnew(ngas_z), cboxold(ngas_z) real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4) ! local variables integer :: ia, idum, idydt_sngldble, ig, l, ntot integer, save :: nrodas_failures = 0 integer, dimension(6) :: inforodas_cur real hmin, hstart, taa, tzz real atolvec(ngas_z), rtolvec(ngas_z), & stot(ngas_z), & yposlimit(ngas_z), yneglimit(ngas_z) real sfixedkpp(nfixed_kppmax), rconstkpp(nreact_kppmax) character*80 msg ! map reaction rate constants (pegasus --> kpp) ! map concentrations (cboxold --> stot) ! dump rates (for debugging) if (iregime .eq. 1) then call cbmz_v02r01_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r01_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) else if (iregime .eq. 2) then call cbmz_v02r02_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r02_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) else if (iregime .eq. 3) then call cbmz_v02r03_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r03_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) else if (iregime .eq. 4) then call cbmz_v02r04_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r04_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) else if (iregime .eq. 5) then call cbmz_v02r05_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r05_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) else call cbmz_v02r06_maprates( rk_m1, rk_m2, rk_m3, rk_m4, & rconstkpp ) call cbmz_v02r06_mapconcs( 0, ntot, stot, sfixedkpp, cboxold ) end if ! set parameters for rodas call do l = 1, ntot atolvec(l) = abs_error rtolvec(l) = rel_error yposlimit(l) = 1.0e20 yneglimit(l) = -1.0e8 end do taa = tstart tzz = tstop hmin = 1.0e-5 hstart = 60.0 idydt_sngldble = 1 ! call rodas integrator ! subr cbmz_v02r06_torodas( ! + ngas, taa, tzz, ! + stot, atol, rtol, yposlimit, yneglimit, ! + hmin, hstart, ! + inforodas_cur, iok, lunerr, idydt_sngldble ) if (iregime .eq. 1) then call cbmz_v02r01_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) else if (iregime .eq. 2) then call cbmz_v02r02_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) else if (iregime .eq. 3) then call cbmz_v02r03_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) else if (iregime .eq. 4) then call cbmz_v02r04_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) else if (iregime .eq. 5) then call cbmz_v02r05_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) else call cbmz_v02r06_torodas( & ntot, taa, tzz, & stot, atolvec, rtolvec, yposlimit, yneglimit, & sfixedkpp, rconstkpp, & hmin, hstart, & inforodas_cur, iok, lunerr, idydt_sngldble ) end if ! increment odeinfo counters if (iok .gt. 0) then if (inforodas_cur(6) .le. 0) then ia = 1 else ia = 2 end if else ia = 3 end if iodestatus_count(ia) = iodestatus_count(ia) + 1 ! do following to avoid overflow of the "inforodas" numbers ! inforodas(2,i) contains rightmost 9 digits of each inforodas number ! inforodas(1,i) contains any higher digits of each inforodas number do ia = 1, 6 idum = inforodas(2,ia) + inforodas_cur(ia) inforodas(1,ia) = inforodas(1,ia) + (idum/1000000000) inforodas(2,ia) = mod(idum, 1000000000) end do ! map concentrations (stot --> cboxnew) if (iregime .eq. 1) then call cbmz_v02r01_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) else if (iregime .eq. 2) then call cbmz_v02r02_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) else if (iregime .eq. 3) then call cbmz_v02r03_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) else if (iregime .eq. 4) then call cbmz_v02r04_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) else if (iregime .eq. 5) then call cbmz_v02r05_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) else call cbmz_v02r06_mapconcs( 1, ntot, stot, sfixedkpp, cboxnew ) end if ! diagnostic output if integration fails OR if i_force_dump > 0 if (iok .gt. 0) then if (i_force_dump .le. 0) goto 20000 else nrodas_failures = nrodas_failures + 1 if (nrodas_failures .gt. 100) goto 20000 end if msg = ' ' call peg_debugmsg( lunout, levdbg_err, msg ) if (iok .gt. 0) then msg = '*** gasodesolver_rodas forced dump' else write(msg,*) '*** gasodesolver_rodas failure no.', & nrodas_failures end if call peg_debugmsg( lunout, levdbg_err, msg ) msg = 'iregime, iok, i, j, k / t' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97010) iregime, iok, isvode, jsvode, ksvode call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) trun call peg_debugmsg( lunout, levdbg_err, msg ) msg = 'inforodas_cur(1-6) =' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97010) inforodas_cur call peg_debugmsg( lunout, levdbg_err, msg ) msg = & 'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza =' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) tstart, tstop, abs_error, rel_error call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) tempbox, pressbox, airdenbox, -99.0 call peg_debugmsg( lunout, levdbg_err, msg ) idum = 0 do ig = nreact_kppmax, 1, -1 if ((idum .eq. 0) .and. (rconstkpp(ig) .ne. 0.0)) idum = ig end do msg = 'ngas_z, nrconst_nonzero =' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97010) ngas_z, idum call peg_debugmsg( lunout, levdbg_err, msg ) msg = 'l, name, cboxold, cboxnew for l=1,ngas_z' call peg_debugmsg( lunout, levdbg_err, msg ) do l = 1, ngas_z write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l) call peg_debugmsg( lunout, levdbg_err, msg ) end do msg = 'rconst for i=1,nrconst_nonzero' call peg_debugmsg( lunout, levdbg_err, msg ) do ia = 1, idum, 4 write(msg,97020) ( rconstkpp(ig), ig = ia, min(ia+3,idum) ) call peg_debugmsg( lunout, levdbg_err, msg ) end do 97010 format( 6i12 ) 97020 format( 4(1pe18.10) ) 97030 format(( i3, 1x, a, 2(1pe18.10) )) ! force non-negative values 20000 do l = 1, ngas_z cboxnew(l) = max( cboxnew(l), 0.0 ) end do return end subroutine gasodesolver_rodas !*********************************************************************** ! < xx.> subr gasodesolver_lsodes ! ! purpose: interface to lsodes ode solver ! ! author : Rahul A. Zaveri ! date : May, 2000 ! !----------------------------------------------------------------------- subroutine gasodesolver_lsodes( tstart, tstop, iok, & isvode, jsvode, ksvode, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info, & i_force_dump, iodestatus_count, & abs_error, rel_error, trun, & tempbox, pressbox, airdenbox, cboxnew, cboxold, & rk_m1, rk_m2, rk_m3, rk_m4 ) use module_data_cbmz use module_cbmz_lsodes_solver, only: lsodes_solver, xsetf, & set_lsodes_common_vars implicit none ! subr arguments real(kind=8) :: trun integer i, iok, isvode, jsvode, ksvode, i_force_dump, iregime, & mgaschem, lunerr, lunout, levdbg_err, levdbg_info integer iodestatus_count(6) real tstart, tstop, abs_error, rel_error real tempbox, pressbox, airdenbox real cboxnew(ngas_z), cboxold(ngas_z) real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4) ! lsodes parameters and local variables integer itoler, itask, iopt, mf, lwm, nrdim, nidim integer nruserpar, niuserpar parameter( itoler = 1, itask = 1, iopt = 1, mf= 222 ) parameter( lwm = 3*ngas_tot*ngas_tot + 12*ngas_tot ) parameter( nrdim = 20 + 9*ngas_tot + lwm ) parameter( nidim = 31 + ngas_tot + ngas_tot*ngas_tot ) parameter( nruserpar = 5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4) parameter( niuserpar = ngas_z + 1 ) integer ia, idum, ig, ioffset, istate, iwork(nidim), l integer ntotvec(1), iuserpar(niuserpar) integer indx(ngas_z) integer, save :: iflagout = 0 integer, save :: nlsodes_failures = 0 real dtchem, rwork(nrdim), stot(ngas_tot) real atolvec(1), rtolvec(1), ruserpar(nruserpar) character*80 msg iok = 1 ! reset call set_lsodes_common_vars() ! sets gas species indices for iregime call setgasindices( iregime, indx ) ! map cboxold --> stot call mapgasspecies( cboxold, stot, 0, iregime, indx ) !---------------------------------------------------------------------- ! set number of species (ntot) for the selected regime for LSODES if (iregime .eq. 1) then ntotvec(1) = ngas_r1 else if (iregime .eq. 2) then ntotvec(1) = ngas_r2 else if (iregime .eq. 3) then ntotvec(1) = ngas_r3 else if (iregime .eq. 4) then ntotvec(1) = ngas_r4 else if (iregime .eq. 5) then ntotvec(1) = ngas_r5 else ntotvec(1) = ngas_r6 end if 100 continue ! set other LSODES parameters... iwork(6) = 1000 ! max iterations for a time step iwork(7) = 1 istate = 1 rwork(6) = dtchem if(iflagout.eq.0)then call xsetf(iflagout) endif atolvec(1) = abs_error rtolvec(1) = rel_error do ig = 1, 5 ruserpar(ig) = ig*7.0 end do ruserpar(1) = cboxold(ih2o_z) ruserpar(2) = cboxold(ich4_z) ruserpar(3) = cboxold(io2_z) ruserpar(4) = cboxold(in2_z) ruserpar(5) = cboxold(ih2_z) ioffset = 5 do ig = 1, nrxn_m1 ruserpar(ioffset+ig) = rk_m1(ig) end do ioffset = ioffset + nrxn_m1 do ig = 1, nrxn_m2 ruserpar(ioffset+ig) = rk_m2(ig) end do ioffset = ioffset + nrxn_m2 do ig = 1, nrxn_m3 ruserpar(ioffset+ig) = rk_m3(ig) end do ioffset = ioffset + nrxn_m3 do ig = 1, nrxn_m4 ruserpar(ioffset+ig) = rk_m4(ig) end do iuserpar(1) = iregime do ig = 1, ngas_z iuserpar(1+ig) = indx(ig) end do call lsodes_solver( & gasode_cbmz, ntotvec, stot, tstart, tstop, & itoler, rtolvec, atolvec, itask, istate, iopt, & rwork, nrdim, iwork, nidim, jcs, mf, & ruserpar, nruserpar, iuserpar, niuserpar ) if (istate .le. 0) iok = -1 ! increment odeinfo counters if (iok .gt. 0) then ia = 4 else ia = 5 end if iodestatus_count(ia) = iodestatus_count(ia) + 1 ! map stot --> cboxnew call mapgasspecies( cboxnew, stot, 1, iregime, indx ) ! do diagnostic output if integration fails OR if i_force_dump > 0 if (iok .gt. 0) then if (i_force_dump .le. 0) goto 20000 else nlsodes_failures = nlsodes_failures + 1 end if msg = ' ' call peg_debugmsg( lunout, levdbg_err, msg ) if (iok .gt. 0) then msg = '*** gasodesolver_lsodes forced dump' else write(msg,*) '*** gasodesolver_lsodes failure no.', & nlsodes_failures end if call peg_debugmsg( lunout, levdbg_err, msg ) msg = 'iregime, iok, i, j, k / t' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97010) iregime, iok, isvode, jsvode, ksvode call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) trun call peg_debugmsg( lunout, levdbg_err, msg ) if (nlsodes_failures .gt. 1000) then write(msg,*) '*** exceeded lsodes failure limit =', 1000 call peg_debugmsg( lunout, levdbg_err, msg ) call peg_error_fatal( lunerr, msg ) end if if (nlsodes_failures .gt. 100) goto 20000 write(msg,*) 'istate -', istate call peg_debugmsg( lunout, levdbg_err, msg ) msg = & 'tstart, tstop, abs_error, rel_error / temp, press, cair, cos_sza =' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) tstart, tstop, abs_error, rel_error call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97020) tempbox, pressbox, airdenbox, -99.0 call peg_debugmsg( lunout, levdbg_err, msg ) idum = nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4 msg = 'ngas_z, nrconst_m1+m2+m3+m4 =' call peg_debugmsg( lunout, levdbg_err, msg ) write(msg,97010) ngas_z, idum call peg_debugmsg( lunout, levdbg_err, msg ) msg = 'l, name, cboxold, cboxnew for l=1,ngas_z' call peg_debugmsg( lunout, levdbg_err, msg ) do l = 1, ngas_z write(msg,97030) l, name_z(l), cboxold(l), cboxnew(l) call peg_debugmsg( lunout, levdbg_err, msg ) end do msg = 'rconst for i=1,nrconst_nonzero' call peg_debugmsg( lunout, levdbg_err, msg ) do ia = 1, nrxn_m1, 4 write(msg,97020) ( rk_m1(ig), ig = ia, min(ia+3,nrxn_m1) ) call peg_debugmsg( lunout, levdbg_err, msg ) end do do ia = 1, nrxn_m2, 4 write(msg,97020) ( rk_m2(ig), ig = ia, min(ia+3,nrxn_m2) ) call peg_debugmsg( lunout, levdbg_err, msg ) end do do ia = 1, nrxn_m3, 4 write(msg,97020) ( rk_m3(ig), ig = ia, min(ia+3,nrxn_m3) ) call peg_debugmsg( lunout, levdbg_err, msg ) end do do ia = 1, nrxn_m4, 4 write(msg,97020) ( rk_m4(ig), ig = ia, min(ia+3,nrxn_m4) ) call peg_debugmsg( lunout, levdbg_err, msg ) end do 97010 format( 6i12 ) 97020 format( 4(1pe18.10) ) 97030 format(( i3, 1x, a, 2(1pe18.10) )) ! force non-negative values 20000 do l = 1, ngas_z cboxnew(l) = max( cboxnew(l), 0.0 ) end do return end subroutine gasodesolver_lsodes !*********************************************************************** ! < 2.> subr selectgasregime ! ! purpose: selects an optimum combination of gas-phase ! mechanisms based on sensitivity of [OH] ! concentrations to some lumped structure ! hydrocarbon groups concentrations and [DMS] ! concentration. ! ! input : cbox = full species concentrations array (mol/cc) ! ! output: iregime = 1 : com ! = 2 : com + urb ! = 3 : com + urb + bio ! = 4 : com + mar ! = 5 : com + urb + mar ! = 6 : com + urb + bio + mar ! ntot = number of gas-phase species in the selected mechanism ! ! author: Rahul A. Zaveri ! date : April 2000 ! !--------------------------------------------------------------------- subroutine selectgasregime( iregime, iregime_forced, cbox, & igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 ) use module_data_cbmz implicit none ! subr arguments integer iregime, iregime_forced integer igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 real cbox(ngas_z) ! local variables integer iwork(6) integer m_m1, m_m2, m_m3, m_m4 real cut_molecpcc cut_molecpcc = 5.e+6 ! molecules/cc ! initialize mechanism flags m_m1 = 1 ! 1 (always) m_m2 = 0 ! 0 or 1 m_m3 = 0 ! 0 or 2 m_m4 = 0 ! 0 or 3 if (igaschem_allowed_m2 .gt. 0) then if (cbox(ipar_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iaone_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(imgly_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ieth_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iolet_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iolei_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ixyl_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(icres_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ito2_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(icro_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iopen_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ionit_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(irooh_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iro2_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(iano2_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(inap_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ixo2_z ) .gt. cut_molecpcc) m_m2 = 1 if (cbox(ixpar_z ) .gt. cut_molecpcc) m_m2 = 1 end if if (igaschem_allowed_m3 .gt. 0) then if (cbox(iisop_z ) .gt. cut_molecpcc) m_m3 = 2 if (cbox(iisoprd_z ) .gt. cut_molecpcc) m_m3 = 2 if (cbox(iisopp_z ) .gt. cut_molecpcc) m_m3 = 2 if (cbox(iisopn_z ) .gt. cut_molecpcc) m_m3 = 2 if (cbox(iisopo2_z ) .gt. cut_molecpcc) m_m3 = 2 end if if (igaschem_allowed_m4 .gt. 0) then if (cbox(idms_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(imsa_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(idmso_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(idmso2_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3so2h_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3sch2oo_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3so2_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3so3_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3so2oo_z ) .gt. cut_molecpcc) m_m4 = 3 if (cbox(ich3so2ch2oo_z) .gt. cut_molecpcc) m_m4 = 3 if (cbox(imtf_z ) .gt. cut_molecpcc) m_m4 = 3 end if iregime = m_m1 + m_m2*((2-m_m3)/2) + m_m3 + m_m4 ! force iregime = iregime_forced if ((iregime_forced .ge. 1) .and. (iregime_forced .le. 6)) & iregime = iregime_forced return end subroutine selectgasregime !*********************************************************************** ! < 3.> subr setgasindices ! ! purpose: sets gas species indices ! ! author : Rahul A. Zaveri ! date : May, 2000 ! !----------------------------------------------------------------------- subroutine setgasindices( iregime, indx ) use module_data_cbmz implicit none ! subr arguments integer iregime, indx(ngas_z) ! local variables integer ilast ilast = 0 indx(:) = -999888777 goto (1,2,3,4,5,6), iregime 1 call setgasindex_m1( ilast, indx ) ! regime 1 return 2 call setgasindex_m1( ilast, indx ) ! regime 2 call setgasindex_m2( ilast, indx ) return 3 call setgasindex_m1( ilast, indx ) ! regime 3 call setgasindex_m2( ilast, indx ) call setgasindex_m3( ilast, indx ) return 4 call setgasindex_m1( ilast, indx ) ! regime 4 call setgasindex_m4( ilast, indx ) return 5 call setgasindex_m1( ilast, indx ) ! regime 5 call setgasindex_m2( ilast, indx ) call setgasindex_m4( ilast, indx ) return 6 call setgasindex_m1( ilast, indx ) ! regime 6 call setgasindex_m2( ilast, indx ) call setgasindex_m3( ilast, indx ) call setgasindex_m4( ilast, indx ) return end subroutine setgasindices !*********************************************************************** ! < 4.> subr gasrateconstants ! ! purpose: calls regime-dependent subrs for calculating ! gas-phase thermal reaction rate constants ! ! author : Rahul A. Zaveri ! date : May, 2000 ! !----------------------------------------------------------------------- subroutine gasrateconstants( iregime, tempbox, cair_mlc, & rk_photo, rk_param, rk_m1, rk_m2, rk_m3, rk_m4 ) use module_data_cbmz implicit none ! subr arguments integer iregime real tempbox, cair_mlc real rk_photo(nphoto), rk_param(nperox) real rk_m1(nrxn_m1), rk_m2(nrxn_m2), rk_m3(nrxn_m3), rk_m4(nrxn_m4) ! iregime=1 --> do m1 ! iregime=2 --> do m1, m2 ! iregime=3 --> do m1, m2, m3 ! iregime=4 --> do m1, --, --, m4 ! iregime=5 --> do m1, m2, --, m4 ! iregime=6 --> do m1, m2, m3, m4 call gasthermrk_m1( tempbox, cair_mlc, & rk_photo, rk_param, rk_m1, rk_m2 ) if ((iregime .eq. 2) .or. & (iregime .eq. 3) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasthermrk_m2( tempbox, cair_mlc, & rk_photo, rk_param, rk_m2 ) if ((iregime .eq. 3) .or. & (iregime .eq. 6)) & call gasthermrk_m3( tempbox, cair_mlc, & rk_photo, rk_param, rk_m3 ) if ((iregime .eq. 4) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasthermrk_m4( tempbox, cair_mlc, & rk_photo, rk_param, rk_m4 ) return end subroutine gasrateconstants !*********************************************************************** ! < 3.> subr setgasindex_m1 ! ! purpose: defines gas species indices for regime 1 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! !----------------------------------------------------------------------- subroutine setgasindex_m1( ilast, indx ) use module_data_cbmz implicit none ! subr arguments integer ilast, indx(ngas_z) indx(ino_z) = 1 indx(ino2_z) = 2 indx(ino3_z) = 3 indx(in2o5_z) = 4 indx(ihono_z) = 5 indx(ihno3_z) = 6 indx(ihno4_z) = 7 indx(io3_z) = 8 indx(io1d_z) = 9 indx(io3p_z) = 10 indx(ioh_z) = 11 indx(iho2_z) = 12 indx(ih2o2_z) = 13 indx(ico_z) = 14 indx(iso2_z) = 15 indx(ih2so4_z) = 16 indx(inh3_z) = 17 indx(ihcl_z) = 18 indx(ic2h6_z) = 19 indx(ich3o2_z) = 20 indx(iethp_z) = 21 indx(ihcho_z) = 22 indx(ich3oh_z) = 23 indx(ic2h5oh_z) = 24 indx(ich3ooh_z) = 25 indx(iethooh_z) = 26 indx(iald2_z) = 27 indx(ihcooh_z) = 28 indx(ircooh_z) = 29 indx(ic2o3_z) = 30 indx(ipan_z) = 31 ilast = indx(ipan_z) return end subroutine setgasindex_m1 !*********************************************************************** ! < 4.> subr setgasindex_m2 ! ! purpose: defines gas species indices for regime 2 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! !----------------------------------------------------------------------- subroutine setgasindex_m2( ilast, indx ) use module_data_cbmz implicit none ! subr arguments integer ilast, indx(ngas_z) indx(ipar_z) = ilast + 1 indx(iaone_z) = ilast + 2 indx(imgly_z) = ilast + 3 indx(ieth_z) = ilast + 4 indx(iolet_z) = ilast + 5 indx(iolei_z) = ilast + 6 indx(itol_z) = ilast + 7 indx(ixyl_z) = ilast + 8 indx(icres_z) = ilast + 9 indx(ito2_z) = ilast + 10 indx(icro_z) = ilast + 11 indx(iopen_z) = ilast + 12 indx(ionit_z) = ilast + 13 ! indx(ipan_z) = ilast + 14 ! indx(ircooh_z) = ilast + 15 indx(irooh_z) = ilast + 14 ! indx(ic2o3_z) = ilast + 17 indx(iro2_z) = ilast + 15 indx(iano2_z) = ilast + 16 indx(inap_z) = ilast + 17 indx(ixo2_z) = ilast + 18 indx(ixpar_z) = ilast + 19 ilast = indx(ixpar_z) return end subroutine setgasindex_m2 !*********************************************************************** ! < 5.> subr setgasindex_m3 ! ! purpose: defines gas species indices for regime 3 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! !----------------------------------------------------------------------- subroutine setgasindex_m3( ilast, indx ) use module_data_cbmz implicit none ! subr arguments integer ilast, indx(ngas_z) indx(iisop_z) = ilast + 1 indx(iisoprd_z) = ilast + 2 indx(iisopp_z) = ilast + 3 indx(iisopn_z) = ilast + 4 indx(iisopo2_z) = ilast + 5 ilast = indx(iisopo2_z) return end subroutine setgasindex_m3 !*********************************************************************** ! < 6.> subr setgasindex_m4 ! ! purpose: defines gas species indices for regime 4 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! !----------------------------------------------------------------------- subroutine setgasindex_m4( ilast, indx ) use module_data_cbmz implicit none ! subr arguments integer ilast, indx(ngas_z) ! indx(idms_z) = ilast + 1 indx(imsa_z) = ilast + 2 indx(idmso_z) = ilast + 3 indx(idmso2_z) = ilast + 4 indx(ich3so2h_z) = ilast + 5 indx(ich3sch2oo_z) = ilast + 6 indx(ich3so2_z) = ilast + 7 indx(ich3so3_z) = ilast + 8 indx(ich3so2oo_z) = ilast + 9 indx(ich3so2ch2oo_z) = ilast + 10 indx(imtf_z) = ilast + 11 ilast = indx(imtf_z) return end subroutine setgasindex_m4 !*********************************************************************** ! < 9.> subr mapgasspecies ! ! purpose: map gas species between stot and cbox arrays ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine mapgasspecies( cbox, stot, imap, iregime, indx ) use module_data_cbmz implicit none ! subr arguments integer imap, iregime, indx(ngas_z) real cbox(ngas_z) real stot(ngas_tot) ! ! goto (1,2,3,4,5,6), iregime ! ! 1 call mapgas_m1( cbox, stot, imap, indx ) return ! ! 2 call mapgas_m1( cbox, stot, imap, indx ) call mapgas_m2( cbox, stot, imap, indx ) return ! ! 3 call mapgas_m1( cbox, stot, imap, indx ) call mapgas_m2( cbox, stot, imap, indx ) call mapgas_m3( cbox, stot, imap, indx ) return ! ! 4 call mapgas_m1( cbox, stot, imap, indx ) call mapgas_m4( cbox, stot, imap, indx ) return ! ! 5 call mapgas_m1( cbox, stot, imap, indx ) call mapgas_m2( cbox, stot, imap, indx ) call mapgas_m4( cbox, stot, imap, indx ) return ! ! 6 call mapgas_m1( cbox, stot, imap, indx ) call mapgas_m2( cbox, stot, imap, indx ) call mapgas_m3( cbox, stot, imap, indx ) call mapgas_m4( cbox, stot, imap, indx ) return end subroutine mapgasspecies !*********************************************************************** ! <10.> subr mapgas_m1 ! ! purpose: maps gas species between stot and cbox arrays for mechanism 1 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine mapgas_m1( cbox, stot, imap, indx ) use module_data_cbmz implicit none ! subr arguments integer imap, indx(ngas_z) real cbox(ngas_z) real stot(ngas_tot) if(imap.eq.0)then ! map cbox --> stot (both molec/cc) stot(indx(ino_z)) = cbox(ino_z) stot(indx(ino2_z)) = cbox(ino2_z) stot(indx(ino3_z)) = cbox(ino3_z) stot(indx(in2o5_z)) = cbox(in2o5_z) stot(indx(ihono_z)) = cbox(ihono_z) stot(indx(ihno3_z)) = cbox(ihno3_z) stot(indx(ihno4_z)) = cbox(ihno4_z) stot(indx(io3_z)) = cbox(io3_z) stot(indx(io1d_z)) = cbox(io1d_z) stot(indx(io3p_z)) = cbox(io3p_z) stot(indx(ioh_z)) = cbox(ioh_z) stot(indx(iho2_z)) = cbox(iho2_z) stot(indx(ih2o2_z)) = cbox(ih2o2_z) stot(indx(ico_z)) = cbox(ico_z) stot(indx(iso2_z)) = cbox(iso2_z) stot(indx(ih2so4_z)) = cbox(ih2so4_z) stot(indx(inh3_z)) = cbox(inh3_z) stot(indx(ihcl_z)) = cbox(ihcl_z) stot(indx(ic2h6_z)) = cbox(ic2h6_z) stot(indx(ich3o2_z)) = cbox(ich3o2_z) stot(indx(iethp_z)) = cbox(iethp_z) stot(indx(ihcho_z)) = cbox(ihcho_z) stot(indx(ich3oh_z)) = cbox(ich3oh_z) stot(indx(ic2h5oh_z)) = cbox(ic2h5oh_z) stot(indx(ich3ooh_z)) = cbox(ich3ooh_z) stot(indx(iethooh_z)) = cbox(iethooh_z) stot(indx(iald2_z)) = cbox(iald2_z) stot(indx(ihcooh_z)) = cbox(ihcooh_z) stot(indx(ircooh_z)) = cbox(ircooh_z) stot(indx(ic2o3_z)) = cbox(ic2o3_z) stot(indx(ipan_z)) = cbox(ipan_z) ! else ! map stot --> cbox (both molec/cc) cbox(ino_z) = stot(indx(ino_z)) cbox(ino2_z) = stot(indx(ino2_z)) cbox(ino3_z) = stot(indx(ino3_z)) cbox(in2o5_z) = stot(indx(in2o5_z)) cbox(ihono_z) = stot(indx(ihono_z)) cbox(ihno3_z) = stot(indx(ihno3_z)) cbox(ihno4_z) = stot(indx(ihno4_z)) cbox(io3_z) = stot(indx(io3_z)) cbox(io1d_z) = stot(indx(io1d_z)) cbox(io3p_z) = stot(indx(io3p_z)) cbox(ioh_z) = stot(indx(ioh_z)) cbox(iho2_z) = stot(indx(iho2_z)) cbox(ih2o2_z) = stot(indx(ih2o2_z)) cbox(ico_z) = stot(indx(ico_z)) cbox(iso2_z) = stot(indx(iso2_z)) cbox(ih2so4_z) = stot(indx(ih2so4_z)) cbox(inh3_z) = stot(indx(inh3_z)) cbox(ihcl_z) = stot(indx(ihcl_z)) cbox(ic2h6_z) = stot(indx(ic2h6_z)) cbox(ich3o2_z) = stot(indx(ich3o2_z)) cbox(iethp_z) = stot(indx(iethp_z)) cbox(ihcho_z) = stot(indx(ihcho_z)) cbox(ich3oh_z) = stot(indx(ich3oh_z)) cbox(ic2h5oh_z) = stot(indx(ic2h5oh_z)) cbox(ich3ooh_z) = stot(indx(ich3ooh_z)) cbox(iethooh_z) = stot(indx(iethooh_z)) cbox(iald2_z) = stot(indx(iald2_z)) cbox(ihcooh_z) = stot(indx(ihcooh_z)) cbox(ircooh_z) = stot(indx(ircooh_z)) cbox(ic2o3_z) = stot(indx(ic2o3_z)) cbox(ipan_z) = stot(indx(ipan_z)) endif return end subroutine mapgas_m1 !*********************************************************************** ! <11.> subr mapgas_m2 ! ! purpose: maps gas species between stot and cbox arrays for mechanism 2 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine mapgas_m2( cbox, stot, imap, indx ) use module_data_cbmz implicit none ! subr arguments integer imap, indx(ngas_z) real cbox(ngas_z) real stot(ngas_tot) if(imap.eq.0)then ! map cbox --> stot (both molec/cc) stot(indx(ipar_z)) = cbox(ipar_z) stot(indx(iaone_z)) = cbox(iaone_z) stot(indx(imgly_z)) = cbox(imgly_z) stot(indx(ieth_z)) = cbox(ieth_z) stot(indx(iolet_z)) = cbox(iolet_z) stot(indx(iolei_z)) = cbox(iolei_z) stot(indx(itol_z)) = cbox(itol_z) stot(indx(ixyl_z)) = cbox(ixyl_z) stot(indx(icres_z)) = cbox(icres_z) stot(indx(ito2_z)) = cbox(ito2_z) stot(indx(icro_z)) = cbox(icro_z) stot(indx(iopen_z)) = cbox(iopen_z) stot(indx(ionit_z)) = cbox(ionit_z) ! stot(indx(ipan_z)) = cbox(ipan_z) ! stot(indx(ircooh_z)) = cbox(ircooh_z) stot(indx(irooh_z)) = cbox(irooh_z) ! stot(indx(ic2o3_z)) = cbox(ic2o3_z) stot(indx(iro2_z)) = cbox(iro2_z) stot(indx(iano2_z)) = cbox(iano2_z) stot(indx(inap_z)) = cbox(inap_z) stot(indx(ixo2_z)) = cbox(ixo2_z) stot(indx(ixpar_z)) = cbox(ixpar_z) ! else ! map stot --> cbox (both molec/cc) cbox(ipar_z) = stot(indx(ipar_z)) cbox(iaone_z) = stot(indx(iaone_z)) cbox(imgly_z) = stot(indx(imgly_z)) cbox(ieth_z) = stot(indx(ieth_z)) cbox(iolet_z) = stot(indx(iolet_z)) cbox(iolei_z) = stot(indx(iolei_z)) cbox(itol_z) = stot(indx(itol_z)) cbox(ixyl_z) = stot(indx(ixyl_z)) cbox(icres_z) = stot(indx(icres_z)) cbox(ito2_z) = stot(indx(ito2_z)) cbox(icro_z) = stot(indx(icro_z)) cbox(iopen_z) = stot(indx(iopen_z)) cbox(ionit_z) = stot(indx(ionit_z)) ! cbox(ipan_z) = stot(indx(ipan_z)) ! cbox(ircooh_z) = stot(indx(ircooh_z)) cbox(irooh_z) = stot(indx(irooh_z)) ! cbox(ic2o3_z) = stot(indx(ic2o3_z)) cbox(iro2_z) = stot(indx(iro2_z)) cbox(iano2_z) = stot(indx(iano2_z)) cbox(inap_z) = stot(indx(inap_z)) cbox(ixo2_z) = stot(indx(ixo2_z)) cbox(ixpar_z) = stot(indx(ixpar_z)) endif return end subroutine mapgas_m2 !*********************************************************************** ! <12.> subr mapgas_m3 ! ! purpose: maps gas species between stot and cbox arrays for mechanism 3 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine mapgas_m3( cbox, stot, imap, indx ) use module_data_cbmz implicit none ! subr arguments integer imap, indx(ngas_z) real cbox(ngas_z) real stot(ngas_tot) if(imap.eq.0)then ! map cbox --> stot (both molec/cc) stot(indx(iisop_z)) = cbox(iisop_z) stot(indx(iisoprd_z)) = cbox(iisoprd_z) stot(indx(iisopp_z)) = cbox(iisopp_z) stot(indx(iisopn_z)) = cbox(iisopn_z) stot(indx(iisopo2_z)) = cbox(iisopo2_z) ! else ! map stot --> cbox (both molec/cc) cbox(iisop_z) = stot(indx(iisop_z)) cbox(iisoprd_z) = stot(indx(iisoprd_z)) cbox(iisopp_z) = stot(indx(iisopp_z)) cbox(iisopn_z) = stot(indx(iisopn_z)) cbox(iisopo2_z) = stot(indx(iisopo2_z)) endif return end subroutine mapgas_m3 !*********************************************************************** ! <13.> subr mapgas_m4 ! ! purpose: maps gas species between stot and cbox arrays for mechanism 4 ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine mapgas_m4( cbox, stot, imap, indx ) use module_data_cbmz implicit none ! subr arguments integer imap, indx(ngas_z) real cbox(ngas_z) real stot(ngas_tot) if(imap.eq.0)then ! map cbox --> stot (both molec/cc) stot(indx(idms_z)) = cbox(idms_z) stot(indx(imsa_z)) = cbox(imsa_z) stot(indx(idmso_z)) = cbox(idmso_z) stot(indx(idmso2_z)) = cbox(idmso2_z) stot(indx(ich3so2h_z)) = cbox(ich3so2h_z) stot(indx(ich3sch2oo_z)) = cbox(ich3sch2oo_z) stot(indx(ich3so2_z)) = cbox(ich3so2_z) stot(indx(ich3so3_z)) = cbox(ich3so3_z) stot(indx(ich3so2oo_z)) = cbox(ich3so2oo_z) stot(indx(ich3so2ch2oo_z))= cbox(ich3so2ch2oo_z) stot(indx(imtf_z)) = cbox(imtf_z) ! else ! map stot --> cbox (both molec/cc) cbox(idms_z) = stot(indx(idms_z)) cbox(imsa_z) = stot(indx(imsa_z)) cbox(idmso_z) = stot(indx(idmso_z)) cbox(idmso2_z) = stot(indx(idmso2_z)) cbox(ich3so2h_z) = stot(indx(ich3so2h_z)) cbox(ich3sch2oo_z) = stot(indx(ich3sch2oo_z)) cbox(ich3so2_z) = stot(indx(ich3so2_z)) cbox(ich3so3_z) = stot(indx(ich3so3_z)) cbox(ich3so2oo_z) = stot(indx(ich3so2oo_z)) cbox(ich3so2ch2oo_z)= stot(indx(ich3so2ch2oo_z)) cbox(imtf_z) = stot(indx(imtf_z)) endif return end subroutine mapgas_m4 !*********************************************************************** ! subr check_userpar ! ! purpose: called by lsodes (external) ! computes time derivatives of species concentrations ds/dt. ! by calling subr. gasrate and ode ! ! author : Rahul A. Zaveri ! date : December 1998 ! !--------------------------------------------------------------------------- subroutine check_userpar( ruserpar, nruserpar, iuserpar, niuserpar ) use module_data_cbmz implicit none ! subr arguments integer nruserpar, niuserpar integer iuserpar(niuserpar) real ruserpar(nruserpar) ! local variables integer i real dum character*80 msg if (nruserpar .ne. (5 + nrxn_m1 + nrxn_m2 + nrxn_m3 + nrxn_m4)) then write(msg,9010) 'nruserpar', -1, nruserpar call wrf_error_fatal( msg ) end if if (niuserpar .ne. (ngas_z + 1)) then write(msg,9010) 'niuserpar', -1, niuserpar call wrf_error_fatal( msg ) end if 9010 format( '*** check_userpar error -- ', a, 1x, i8, 1x, i8 ) return end subroutine check_userpar !*********************************************************************** ! <14.> subr gasode_cbmz ! ! purpose: called by lsodes (external) ! computes time derivatives of species concentrations ds/dt. ! by calling subr. gasrate and ode ! ! author : Rahul A. Zaveri ! date : December 1998 ! !--------------------------------------------------------------------------- subroutine gasode_cbmz( ntot, tt, s, sdot, & ruserpar, nruserpar, iuserpar, niuserpar ) use module_data_cbmz implicit none ! subr arguments integer ntot, nruserpar, niuserpar integer iuserpar(niuserpar) real tt real s(ngas_tot), sdot(ngas_tot), ruserpar(nruserpar) ! local variables integer ig, ioffset, iregime, irxn integer indx(ngas_z) real h2o, ch4, oxygen, nitrogen, hydrogen real rk_m1(nrxn_m1), r_m1(nrxn_m1) real rk_m2(nrxn_m2), r_m2(nrxn_m2) real rk_m3(nrxn_m3), r_m3(nrxn_m3) real rk_m4(nrxn_m4), r_m4(nrxn_m4) real p_m1(ngas_tot), d_m1(ngas_tot) real p_m2(ngas_tot), d_m2(ngas_tot) real p_m3(ngas_tot), d_m3(ngas_tot) real p_m4(ngas_tot), d_m4(ngas_tot) ! test on userpar call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar ) iregime = iuserpar(1) do ig = 1, ngas_z indx(ig) = iuserpar(ig+1) end do h2o = ruserpar(1) ch4 = ruserpar(2) oxygen = ruserpar(3) nitrogen = ruserpar(4) hydrogen = ruserpar(5) ioffset = 5 do ig = 1, nrxn_m1 rk_m1(ig) = ruserpar(ioffset+ig) end do ioffset = ioffset + nrxn_m1 do ig = 1, nrxn_m2 rk_m2(ig) = ruserpar(ioffset+ig) end do ioffset = ioffset + nrxn_m2 do ig = 1, nrxn_m3 rk_m3(ig) = ruserpar(ioffset+ig) end do ioffset = ioffset + nrxn_m3 do ig = 1, nrxn_m4 rk_m4(ig) = ruserpar(ioffset+ig) end do ! initialize to zero do irxn=1,nrxn_m1 r_m1(irxn) = 0. enddo do irxn=1,nrxn_m2 r_m2(irxn) = 0. enddo do irxn=1,nrxn_m3 r_m3(irxn) = 0. enddo do irxn=1,nrxn_m4 r_m4(irxn) = 0. enddo ! ! ! initialize to zero do ig=1,ngas_tot p_m1(ig) = 0. p_m2(ig) = 0. p_m3(ig) = 0. p_m4(ig) = 0. ! d_m1(ig) = 0. d_m2(ig) = 0. d_m3(ig) = 0. d_m4(ig) = 0. enddo ! iregime=1 --> do m1 ! iregime=2 --> do m1, m2 ! iregime=3 --> do m1, m2, m3 ! iregime=4 --> do m1, --, --, m4 ! iregime=5 --> do m1, m2, --, m4 ! iregime=6 --> do m1, m2, m3, m4 call gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2, & h2o, ch4, oxygen, nitrogen, hydrogen ) if ((iregime .eq. 2) .or. & (iregime .eq. 3) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasrate_m2( indx, s, r_m2, rk_m2 ) if ((iregime .eq. 3) .or. & (iregime .eq. 6)) & call gasrate_m3( indx, s, r_m3, rk_m3 ) if ((iregime .eq. 4) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasrate_m4( indx, s, r_m4, rk_m4, oxygen ) call gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 ) if ((iregime .eq. 2) .or. & (iregime .eq. 3) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasode_m2( indx, r_m2, p_m2, d_m2 ) if ((iregime .eq. 3) .or. & (iregime .eq. 6)) & call gasode_m3( indx, r_m3, p_m3, d_m3 ) if ((iregime .eq. 4) .or. & (iregime .eq. 5) .or. & (iregime .eq. 6)) & call gasode_m4( indx, r_m4, p_m4, d_m4 ) if (iregime .eq. 1) then ! regime = 1 do ig = 1, ngas_r1 sdot(ig) = real( dble(p_m1(ig)) - & dble(d_m1(ig)) ) end do else if (iregime .eq. 2) then ! regime = 2 do ig = 1, ngas_r2 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)) - & dble(d_m1(ig)+d_m2(ig)) ) end do else if (iregime .eq. 3) then ! regime = 3 do ig = 1, ngas_r3 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)) - & dble(d_m1(ig)+d_m2(ig)+d_m3(ig)) ) end do else if (iregime .eq. 4) then ! regime = 4 do ig = 1, ngas_r4 sdot(ig) = real( dble(p_m1(ig)+p_m4(ig)) - & dble(d_m1(ig)+d_m4(ig)) ) end do else if (iregime .eq. 5) then ! regime = 5 do ig = 1, ngas_r5 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m4(ig)) - & dble(d_m1(ig)+d_m2(ig)+d_m4(ig)) ) end do else if (iregime .eq. 6) then ! regime = 6 do ig = 1, ngas_r6 sdot(ig) = real( dble(p_m1(ig)+p_m2(ig)+p_m3(ig)+p_m4(ig)) - & dble(d_m1(ig)+d_m2(ig)+d_m3(ig)+d_m4(ig)) ) end do end if return end subroutine gasode_cbmz !*********************************************************************** ! <15.> subr jcs ! ! purpose: external dummy jacobian evaluation for LSODES (when mf=222) ! !----------------------------------------------------------------------- subroutine jcs( ngas, tt, s, j, ian, jan, pdj, & ruserpar, nruserpar, iuserpar, niuserpar ) implicit none integer ngas, j, ian(*), jan(*), nruserpar, niuserpar integer iuserpar(niuserpar) real tt, s(*), pdj(*) real ruserpar(nruserpar) ! test on userpar call check_userpar( ruserpar, nruserpar, iuserpar, niuserpar ) return end subroutine jcs !*********************************************************************** ! <16.> subr gasode_m1 ! ! purpose: updates production and destruction rates for mechanism 1 ! background troposphere ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine gasode_m1( indx, r_m1, p_m1, d_m1, r_m2 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real r_m1(nrxn_m1), p_m1(ngas_tot), d_m1(ngas_tot) real r_m2(nrxn_m2) p_m1(indx(ino_z))= r_m1(1)+0.11*r_m1(2) & +r_m1(3)+r_m1(15)+r_m1(38) d_m1(indx(ino_z))= r_m1(17)+r_m1(18)+r_m1(23) & +r_m1(33)+r_m1(37)+r_m1(57) & +r_m1(58) & +r_m2(34) ! p_m1(indx(ino2_z))= 0.89*r_m1(2)+r_m1(4) & +r_m1(5)+r_m1(6)+r_m1(17) & +r_m1(18)+r_m1(25) & +r_m1(26)+r_m1(28) & +r_m1(33)+r_m1(36) & +r_m1(37)+r_m1(37) & +r_m1(38)+r_m1(40) & +r_m1(40)+.7*r_m1(41) & +r_m1(43)+r_m1(57) & +r_m1(58)+r_m1(59) & +r_m1(60) & +r_m2(32)+r_m2(34)+r_m2(39) d_m1(indx(ino2_z))= r_m1(1)+r_m1(15)+r_m1(16) & +r_m1(19)+r_m1(24) & +r_m1(34)+r_m1(35) & +r_m1(38)+r_m1(39) & +r_m2(31) ! p_m1(indx(ino3_z))= r_m1(6)+r_m1(16)+r_m1(19) & +r_m1(27)+r_m1(43) d_m1(indx(ino3_z))= r_m1(2)+r_m1(25)+r_m1(37) & +r_m1(38)+r_m1(39) & +r_m1(40)+r_m1(40) & +r_m1(41)+r_m1(52) & +r_m1(59)+r_m1(60) & +r_m2(4)+r_m2(39) ! p_m1(indx(in2o5_z))= r_m1(39) d_m1(indx(in2o5_z))= r_m1(6)+r_m1(42) & +r_m1(43) ! p_m1(indx(ihono_z))= r_m1(23)+r_m1(35) d_m1(indx(ihono_z))= r_m1(3)+r_m1(26) ! p_m1(indx(ihno3_z))= r_m1(24)+.3*r_m1(41) & +r_m1(42)+r_m1(42) & +r_m1(52) & +r_m2(4) d_m1(indx(ihno3_z))= r_m1(4)+r_m1(27) ! p_m1(indx(ihno4_z))= r_m1(34) d_m1(indx(ihno4_z))= r_m1(5)+r_m1(28) & +r_m1(36) ! p_m1(indx(io3_z))= r_m1(13) & +.4*r_m2(44) d_m1(indx(io3_z))= r_m1(7)+r_m1(8)+r_m1(14) & +r_m1(18)+r_m1(19)+r_m1(20) & +r_m1(21) ! p_m1(indx(io1d_z))= r_m1(8) d_m1(indx(io1d_z))= r_m1(10)+r_m1(11) & +r_m1(12) ! p_m1(indx(io3p_z))= r_m1(1)+0.89*r_m1(2) & +r_m1(7)+r_m1(10)+r_m1(11) d_m1(indx(io3p_z))= r_m1(13)+r_m1(14) & +r_m1(15)+r_m1(16) & +r_m1(17) ! p_m1(indx(ioh_z))= r_m1(3)+r_m1(4)+2*r_m1(9) & +2*r_m1(12)+r_m1(21) & +r_m1(33)+.7*r_m1(41) & +r_m1(53)+r_m1(54)+.3*r_m1(55) & +.5*r_m1(56) d_m1(indx(ioh_z))= r_m1(20)+r_m1(22)+r_m1(23) & +r_m1(24)+r_m1(25)+r_m1(26) & +r_m1(27)+r_m1(28)+r_m1(29) & +r_m1(30)+r_m1(44)+r_m1(45) & +r_m1(46)+r_m1(47)+r_m1(48) & +r_m1(51)+r_m1(55)+r_m1(56) & +r_m1(65) & +r_m2(3) ! p_m1(indx(iho2_z))= r_m1(5)+r_m1(20)+r_m1(22) & +r_m1(25)+r_m1(30) & +r_m1(36)+r_m1(44) & +r_m1(45)+r_m1(48) & +2*r_m1(49)+r_m1(51) & +r_m1(52)+r_m1(53) & +r_m1(54)+r_m1(57) & +r_m1(58)+r_m1(59) & +r_m1(60)+.32*r_m1(63) & +.6*r_m1(64)+r_m1(65) & +r_m2(2) d_m1(indx(iho2_z))= r_m1(21)+r_m1(29) & +r_m1(31)+r_m1(31) & +r_m1(32)+r_m1(32) & +r_m1(33)+r_m1(34) & +r_m1(35)+r_m1(41) & +r_m1(61)+r_m1(62) & +r_m2(44) ! p_m1(indx(ih2o2_z))= r_m1(31)+r_m1(32) d_m1(indx(ih2o2_z))= r_m1(9)+r_m1(30) ! p_m1(indx(ico_z))= r_m1(49)+r_m1(50)+r_m1(51) & +r_m1(52) & +r_m2(2) d_m1(indx(ico_z))= r_m1(44) ! p_m1(indx(iso2_z))= 0.0 d_m1(indx(iso2_z))= r_m1(45) ! p_m1(indx(ih2so4_z))= r_m1(45) d_m1(indx(ih2so4_z))= 0.0 ! p_m1(indx(inh3_z))= 0.0 d_m1(indx(inh3_z))= 0.0 ! p_m1(indx(ihcl_z))= 0.0 d_m1(indx(ihcl_z))= 0.0 ! p_m1(indx(ic2h6_z))= .2*r_m1(64) d_m1(indx(ic2h6_z))= r_m1(47) ! p_m1(indx(ich3o2_z))= r_m1(46)+.7*r_m1(55) & +r_m2(2)+r_m2(34)+r_m2(39)+r_m2(49) d_m1(indx(ich3o2_z))= r_m1(57)+r_m1(59) & +r_m1(61)+r_m1(63) ! p_m1(indx(iethp_z))= r_m1(47)+.5*r_m1(56) d_m1(indx(iethp_z))= r_m1(58)+r_m1(60) & +r_m1(62)+r_m1(64) ! p_m1(indx(ihcho_z))= r_m1(48)+r_m1(53) & +.3*r_m1(55)+r_m1(57) & +r_m1(59)+.66*r_m1(63) d_m1(indx(ihcho_z))= r_m1(49)+r_m1(50) & +r_m1(51)+r_m1(52) ! p_m1(indx(ich3oh_z))= .34*r_m1(63) d_m1(indx(ich3oh_z))= r_m1(48) ! p_m1(indx(ic2h5oh_z))= 0.0 d_m1(indx(ic2h5oh_z))= r_m1(65) ! p_m1(indx(ich3ooh_z))= r_m1(61) d_m1(indx(ich3ooh_z))= r_m1(53)+r_m1(55) ! p_m1(indx(iethooh_z))= r_m1(62) d_m1(indx(iethooh_z))= r_m1(54)+r_m1(56) ! p_m1(indx(iald2_z))= r_m1(54)+.5*r_m1(56) & +r_m1(58)+r_m1(60) & +.8*r_m1(64)+r_m1(65) d_m1(indx(iald2_z))= r_m2(2)+r_m2(3)+r_m2(4) ! p_m1(indx(ihcooh_z))= 0.0 d_m1(indx(ihcooh_z))= 0.0 ! p_m1(indx(ircooh_z))= .4*r_m2(44) d_m1(indx(ircooh_z))= 0.0 ! p_m1(indx(ic2o3_z))= r_m2(3)+r_m2(4)+r_m2(32) d_m1(indx(ic2o3_z))= r_m2(31)+r_m2(34)+r_m2(39)+r_m2(44)+r_m2(49) ! p_m1(indx(ipan_z))= r_m2(31) d_m1(indx(ipan_z))= r_m2(32) return end subroutine gasode_m1 !*********************************************************************** ! <17.> subr gasode_m2 ! ! purpose: updates production and destruction rates for mechanism 2 ! anthropogenic hydrocarbons ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine gasode_m2( indx, r_m2, p_m2, d_m2 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real r_m2(nrxn_m2), p_m2(ngas_tot), d_m2(ngas_tot) p_m2(indx(ino_z))= 0.0 d_m2(indx(ino_z))= r_m2(20)+r_m2(33) +r_m2(35) & +r_m2(36)+r_m2(37) ! p_m2(indx(ino2_z))= .95*r_m2(20)+r_m2(30) +.84*r_m2(33) & +r_m2(35)+1.5*r_m2(36)+r_m2(37) & +r_m2(38) +r_m2(40)+1.5*r_m2(41)+r_m2(42) & +.5*r_m2(51) d_m2(indx(ino2_z))= r_m2(23) ! p_m2(indx(ino3_z))= 0.0 d_m2(indx(ino3_z))= +r_m2(9)+r_m2(16)+r_m2(17)+r_m2(22) & +r_m2(38) +r_m2(40)+r_m2(41)+r_m2(42) ! p_m2(indx(ihno3_z))= +r_m2(9)+r_m2(22) d_m2(indx(ihno3_z))= 0.0 ! p_m2(indx(io3_z))= 0.0 d_m2(indx(io3_z))= r_m2(10)+r_m2(12)+r_m2(13)+r_m2(26) ! p_m2(indx(ioh_z))= .12*r_m2(10)+.33*r_m2(12)+.60*r_m2(13) & +.08*r_m2(26)+r_m2(27)+.23*r_m2(28) d_m2(indx(ioh_z))= r_m2(1) +r_m2(6)+r_m2(8)+r_m2(11) & +r_m2(14)+r_m2(15)+r_m2(18)+r_m2(19)+r_m2(21) & +r_m2(24)+r_m2(28)+r_m2(29) ! p_m2(indx(iho2_z))= +r_m2(7)+.22*r_m2(10)+r_m2(11) & +.26*r_m2(12)+.22*r_m2(13)+r_m2(14)+r_m2(15) & +.2*r_m2(18)+.55*r_m2(19)+.95*r_m2(20) & +.6*r_m2(21)+2*r_m2(24)+r_m2(25)+.76*r_m2(26) & +.9*r_m2(27)+.9*r_m2(30)+.76*r_m2(33)+.5*r_m2(36) & +.9*r_m2(38)+.5*r_m2(41)+.54*r_m2(48) d_m2(indx(iho2_z))= r_m2(43) +r_m2(45)+r_m2(46)+r_m2(47) ! p_m2(indx(ico_z))= +r_m2(7)+r_m2(9)+.24*r_m2(10) & +.31*r_m2(12)+.30*r_m2(13)+2*r_m2(24)+r_m2(25) & +.69*r_m2(26) d_m2(indx(ico_z))= 0.0 ! p_m2(indx(ipar_z))= 1.1*r_m2(19) d_m2(indx(ipar_z))= r_m2(1) + r_m2(53) ! p_m2(indx(ich3oh_z))= .03*r_m2(12)+.04*r_m2(13) d_m2(indx(ich3oh_z))= 0.0 ! p_m2(indx(ihcho_z))= r_m2(10)+1.56*r_m2(11)+.57*r_m2(12)+r_m2(14) & +r_m2(24)+.7*r_m2(26)+r_m2(35)+.5*r_m2(36) & +r_m2(40)+.5*r_m2(41)+.7*r_m2(50)+.5*r_m2(51) d_m2(indx(ihcho_z))= 0.0 ! p_m2(indx(iald2_z))= .22*r_m2(11)+.47*r_m2(12)+1.03*r_m2(13) & +r_m2(14)+1.77*r_m2(15)+.03*r_m2(26)+.3*r_m2(27) & +.04*r_m2(28)+.3*r_m2(30)+.25*r_m2(33)+.5*r_m2(36) & +.3*r_m2(38)+.5*r_m2(41)+.21*r_m2(48)+.5*r_m2(51) d_m2(indx(iald2_z))= 0.0 ! p_m2(indx(ihcooh_z))= .52*r_m2(10)+.22*r_m2(12) d_m2(indx(ihcooh_z))= 0.0 ! p_m2(indx(iaone_z))= .07*r_m2(13)+.23*r_m2(15)+.74*r_m2(27) & +.74*r_m2(30)+.62*r_m2(33)+.74*r_m2(38) & +.57*r_m2(48)+.15*r_m2(50) d_m2(indx(iaone_z))= r_m2(5)+r_m2(6) ! p_m2(indx(imgly_z))= .04*r_m2(12)+.07*r_m2(13)+.8*r_m2(19) & +.2*r_m2(26)+.19*r_m2(28)+.15*r_m2(50) d_m2(indx(imgly_z))= r_m2(7)+r_m2(8)+r_m2(9) ! p_m2(indx(ieth_z))= 0.0 d_m2(indx(ieth_z))= r_m2(10)+r_m2(11) ! p_m2(indx(iolet_z))= 0.0 d_m2(indx(iolet_z))= r_m2(12)+r_m2(14)+r_m2(16) ! p_m2(indx(iolei_z))= 0.0 d_m2(indx(iolei_z))= r_m2(13)+r_m2(15)+r_m2(17) ! p_m2(indx(itol_z))= 0.0 d_m2(indx(itol_z))= r_m2(18) ! p_m2(indx(ixyl_z))= 0.0 d_m2(indx(ixyl_z))= r_m2(19) ! p_m2(indx(icres_z))= .12*r_m2(18)+.05*r_m2(19) d_m2(indx(icres_z))= r_m2(21)+r_m2(22) ! p_m2(indx(ito2_z))= .8*r_m2(18)+.45*r_m2(19) d_m2(indx(ito2_z))= r_m2(20) ! p_m2(indx(icro_z))= .4*r_m2(21)+r_m2(22) d_m2(indx(icro_z))= r_m2(23) ! p_m2(indx(iopen_z))= .95*r_m2(20)+.3*r_m2(21) d_m2(indx(iopen_z))= r_m2(24)+r_m2(25)+r_m2(26) ! p_m2(indx(ionit_z))= .05*r_m2(20)+r_m2(23)+.16*r_m2(33) & +.5*r_m2(36)+.5*r_m2(41)+r_m2(46)+.5*r_m2(51) d_m2(indx(ionit_z))= r_m2(29)+r_m2(30) ! p_m2(indx(ipan_z))= 0.0 d_m2(indx(ipan_z))= 0.0 ! p_m2(indx(ircooh_z))= .09*r_m2(12)+.16*r_m2(13) d_m2(indx(ircooh_z))= 0.0 ! p_m2(indx(irooh_z))= r_m2(43)+r_m2(45) d_m2(indx(irooh_z))= r_m2(27)+r_m2(28) ! p_m2(indx(ich3o2_z))= +r_m2(5)+.07*r_m2(12)+.10*r_m2(13) d_m2(indx(ich3o2_z))= 0.0 ! p_m2(indx(iethp_z))= .06*r_m2(12)+.05*r_m2(13)+.1*r_m2(27) & +.1*r_m2(30)+.08*r_m2(33)+.1*r_m2(38)+.06*r_m2(48) d_m2(indx(iethp_z))= 0.0 ! p_m2(indx(ic2o3_z))= +r_m2(5)+r_m2(7)+r_m2(8) & +r_m2(9)+.13*r_m2(12)+.19*r_m2(13)+r_m2(24) & +r_m2(25)+.62*r_m2(26) +r_m2(35) & +r_m2(40)+.7*r_m2(50) d_m2(indx(ic2o3_z))= 0.0 ! p_m2(indx(iro2_z))= r_m2(1)+.03*r_m2(12)+.09*r_m2(13)+.77*r_m2(28) d_m2(indx(iro2_z))= r_m2(33)+r_m2(38)+r_m2(43)+r_m2(48) ! p_m2(indx(iano2_z))= r_m2(6)+.11*r_m2(13) d_m2(indx(iano2_z))= r_m2(35)+r_m2(40)+r_m2(45)+r_m2(50) ! p_m2(indx(inap_z))= r_m2(16)+r_m2(17)+r_m2(29) d_m2(indx(inap_z))= r_m2(36)+r_m2(41)+r_m2(46)+r_m2(51) ! p_m2(indx(ixo2_z))= r_m2(8)+r_m2(11)+r_m2(14)+r_m2(15) & +.08*r_m2(18)+.5*r_m2(19)+.6*r_m2(21) & +r_m2(24)+.03*r_m2(26)+.4*r_m2(27)+.4*r_m2(30) & +.34*r_m2(33)+.4*r_m2(38)+.24*r_m2(48) d_m2(indx(ixo2_z))= r_m2(37)+r_m2(42)+r_m2(47)+r_m2(52) p_m2(indx(ixpar_z))= 1.06*r_m2(12)+2.26*r_m2(13)+r_m2(14) & +2.23*r_m2(15)+1.98*r_m2(27)+.42*r_m2(28) & +1.98*r_m2(30)+1.68*r_m2(33)+r_m2(36) & +1.98*r_m2(38)+r_m2(41)+1.25*r_m2(48)+r_m2(51) d_m2(indx(ixpar_z))= r_m2(53) ! return end subroutine gasode_m2 !*********************************************************************** ! <18.> subr gasode_m3 ! ! purpose: updates production and destruction rates for mechanism 3 ! isoprene ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine gasode_m3( indx, r_m3, p_m3, d_m3 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real r_m3(nrxn_m3), p_m3(ngas_tot), d_m3(ngas_tot) p_m3(indx(ino_z))= 0.0 d_m3(indx(ino_z))= r_m3(8)+r_m3(9)+r_m3(10) ! p_m3(indx(ino2_z))= .91*r_m3(8)+1.2*r_m3(9)+r_m3(10) d_m3(indx(ino2_z))= 0.0 ! p_m3(indx(ino3_z))= 0.0 d_m3(indx(ino3_z))= r_m3(3)+r_m3(7) ! p_m3(indx(ihno3_z))= .07*r_m3(7) d_m3(indx(ihno3_z))= 0.0 ! p_m3(indx(io3_z))= 0.0 d_m3(indx(io3_z))= r_m3(2)+r_m3(6) ! p_m3(indx(ioh_z))= .27*r_m3(2)+.27*r_m3(6) d_m3(indx(ioh_z))= r_m3(1)+r_m3(5) ! p_m3(indx(iho2_z))= .07*r_m3(2)+.33*r_m3(4)+.1*r_m3(6)+.93*r_m3(7) & +.91*r_m3(8)+.8*r_m3(9)+r_m3(10) d_m3(indx(iho2_z))= r_m3(11)+r_m3(12)+r_m3(13) ! p_m3(indx(ico_z))= .07*r_m3(2)+.33*r_m3(4)+.16*r_m3(6)+.64*r_m3(7) & +.59*r_m3(10) d_m3(indx(ico_z))= 0.0 ! p_m3(indx(ipar_z))= 1.86*r_m3(7)+0.18*r_m3(8)+1.6*r_m3(9)+2*r_m3(12) & +2*r_m3(15) d_m3(indx(ipar_z))= 0.0 ! p_m3(indx(ihcho_z))= .6*r_m3(2)+.2*r_m3(4)+.15*r_m3(6)+.28*r_m3(7) & +.63*r_m3(8)+.25*r_m3(10) d_m3(indx(ihcho_z))= 0.0 ! p_m3(indx(iald2_z))= .15*r_m3(2)+.07*r_m3(4)+.02*r_m3(6)+.28*r_m3(7) & +.8*r_m3(9)+.55*r_m3(10)+r_m3(15)+.5*r_m3(16) d_m3(indx(iald2_z))= 0.0 ! p_m3(indx(iaone_z))= .03*r_m3(4)+.09*r_m3(6)+.63*r_m3(10)+.5*r_m3(16) d_m3(indx(iaone_z))= 0.0 ! p_m3(indx(imgly_z))= .85*r_m3(6)+.34*r_m3(10) d_m3(indx(imgly_z))= 0.0 ! p_m3(indx(ionit_z))= .93*r_m3(7)+.09*r_m3(8)+.8*r_m3(9)+r_m3(12) & +r_m3(15) d_m3(indx(ionit_z))= 0.0 ! p_m3(indx(ihcooh_z))= .39*r_m3(2)+.46*r_m3(6) d_m3(indx(ihcooh_z))= 0.0 ! p_m3(indx(irooh_z))= r_m3(11)+r_m3(13) d_m3(indx(irooh_z))= 0.0 ! p_m3(indx(ich3o2_z))= .7*r_m3(4)+.05*r_m3(6) d_m3(indx(ich3o2_z))= 0.0 ! p_m3(indx(ic2o3_z))= .2*r_m3(2)+.97*r_m3(4)+.5*r_m3(5)+.11*r_m3(6) & +.07*r_m3(7) d_m3(indx(ic2o3_z))= 0.0 ! p_m3(indx(ixo2_z))= .08*r_m3(1)+.2*r_m3(2)+.2*r_m3(5)+.07*r_m3(6) & +.93*r_m3(7) d_m3(indx(ixo2_z))= 0.0 ! p_m3(indx(iisop_z))= 0.0 d_m3(indx(iisop_z))= r_m3(1)+r_m3(2)+r_m3(3) ! p_m3(indx(iisoprd_z))= .65*r_m3(2)+.91*r_m3(8)+.2*r_m3(9)+r_m3(14) d_m3(indx(iisoprd_z))= r_m3(4)+r_m3(5)+r_m3(6)+r_m3(7) ! p_m3(indx(iisopp_z))= r_m3(1) d_m3(indx(iisopp_z))= r_m3(8)+r_m3(11)+r_m3(14) ! p_m3(indx(iisopn_z))= r_m3(3) d_m3(indx(iisopn_z))= r_m3(9)+r_m3(12)+r_m3(15) ! p_m3(indx(iisopo2_z))= .5*r_m3(5) d_m3(indx(iisopo2_z))= r_m3(10)+r_m3(13)+r_m3(16) ! return end subroutine gasode_m3 !*********************************************************************** ! <19.> subr gasode_m4 ! ! purpose: updates production and destruction rates for mechanism 4 ! dimethylsulfide ! ! author : Rahul A. Zaveri ! date : December, 1998 ! ! ---------------------------------------------------------------------- subroutine gasode_m4( indx, r_m4, p_m4, d_m4 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real r_m4(nrxn_m4), p_m4(ngas_tot), d_m4(ngas_tot) p_m4(indx(ino_z))= r_m4(19) d_m4(indx(ino_z))= r_m4(5)+r_m4(11)+r_m4(26)+r_m4(30) ! p_m4(indx(ino2_z))= r_m4(5)+r_m4(11)+r_m4(26) d_m4(indx(ino2_z))= r_m4(19)+r_m4(29) ! p_m4(indx(ino3_z))= 0.0 d_m4(indx(ino3_z))= r_m4(2)+r_m4(14) ! p_m4(indx(ihono_z))= r_m4(30) d_m4(indx(ihono_z))= 0.0 ! p_m4(indx(ihno3_z))= r_m4(2)+r_m4(14)+r_m4(29) d_m4(indx(ihno3_z))= 0.0 ! p_m4(indx(io3_z))= 0.0 d_m4(indx(io3_z))= r_m4(20) ! p_m4(indx(io3p_z))= 0.0 d_m4(indx(io3p_z))= r_m4(3) ! p_m4(indx(ioh_z))= r_m4(21) d_m4(indx(ioh_z))= r_m4(1)+r_m4(4)+r_m4(9)+r_m4(10)+r_m4(16)+r_m4(23) ! p_m4(indx(iho2_z))= .965*r_m4(4)+r_m4(6)+.27*r_m4(9)+r_m4(12)+r_m4(22) & +r_m4(27)+r_m4(32) d_m4(indx(iho2_z))= r_m4(13)+r_m4(21)+r_m4(31) ! p_m4(indx(ih2o2_z))= r_m4(13) d_m4(indx(ih2o2_z))= 0.0 ! p_m4(indx(ico_z))= r_m4(32) d_m4(indx(ico_z))= 0.0 ! p_m4(indx(iso2_z))= r_m4(18) d_m4(indx(iso2_z))= 0.0 ! p_m4(indx(ih2so4_z))= r_m4(28) d_m4(indx(ih2so4_z))= 0.0 ! p_m4(indx(ihcho_z))= r_m4(5)+2*r_m4(6)+r_m4(7)+r_m4(11)+2*r_m4(12) & +r_m4(22)+r_m4(27) d_m4(indx(ihcho_z))= r_m4(32) ! p_m4(indx(ich3o2_z))= r_m4(3)+.035*r_m4(4)+.73*r_m4(9)+r_m4(18)+r_m4(28) d_m4(indx(ich3o2_z))= r_m4(6)+r_m4(12)+r_m4(15)+r_m4(22)+r_m4(27) ! p_m4(indx(ich3ooh_z))= r_m4(15) d_m4(indx(ich3ooh_z))= 0.0 ! p_m4(indx(idms_z))= 0.0 d_m4(indx(idms_z))= r_m4(1)+r_m4(2)+r_m4(3)+r_m4(4) ! p_m4(indx(imsa_z))= r_m4(17)+r_m4(23)+r_m4(29)+r_m4(30)+r_m4(31)+r_m4(32) d_m4(indx(imsa_z))= 0.0 ! p_m4(indx(idmso_z))= .965*r_m4(4) d_m4(indx(idmso_z))= r_m4(9) ! p_m4(indx(idmso2_z))= .27*r_m4(9) d_m4(indx(idmso2_z))= r_m4(10) ! p_m4(indx(ich3so2h_z))= .73*r_m4(9) d_m4(indx(ich3so2h_z))= r_m4(13)+r_m4(14)+r_m4(15)+r_m4(16)+r_m4(17) ! p_m4(indx(ich3sch2oo_z))= r_m4(1)+r_m4(2) d_m4(indx(ich3sch2oo_z))= r_m4(5)+r_m4(6)+r_m4(7)+r_m4(8)+r_m4(8) ! p_m4(indx(ich3so2_z))= r_m4(3)+.035*r_m4(4)+r_m4(5)+r_m4(6)+r_m4(7) & +1.85*r_m4(8) & +r_m4(11)+r_m4(12)+r_m4(13)+r_m4(14)+r_m4(15) & +r_m4(16)+r_m4(17)+r_m4(25) d_m4(indx(ich3so2_z))= r_m4(7)+r_m4(18)+r_m4(19)+r_m4(20)+r_m4(21) & +r_m4(22)+r_m4(23)+r_m4(24) ! p_m4(indx(ich3so3_z))= r_m4(7)+r_m4(19)+r_m4(20)+r_m4(21)+r_m4(22) & +r_m4(26)+r_m4(27) d_m4(indx(ich3so3_z))= r_m4(17)+r_m4(28)+r_m4(29)+r_m4(30)+r_m4(31) & +r_m4(32) ! p_m4(indx(ich3so2oo_z))= r_m4(24) d_m4(indx(ich3so2oo_z))= r_m4(25)+r_m4(26)+r_m4(27) ! p_m4(indx(ich3so2ch2oo_z))= r_m4(10) d_m4(indx(ich3so2ch2oo_z))= r_m4(11)+r_m4(12) ! p_m4(indx(imtf_z))= .15*r_m4(8) d_m4(indx(imtf_z))= 0.0 ! return end subroutine gasode_m4 !*********************************************************************** ! <20.> subr gasrate_m1 ! ! purpose: computes reaction rates for mechanism 1 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- subroutine gasrate_m1( indx, s, r_m1, r_m2, rk_m1, rk_m2, & h2o, ch4, oxygen, nitrogen, hydrogen ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real s(ngas_tot), r_m1(nrxn_m1), r_m2(nrxn_m2) real rk_m1(nrxn_m1), rk_m2(nrxn_m2) real h2o, ch4, oxygen, nitrogen, hydrogen r_m1(1) = rk_m1(1)*s(indx(ino2_z)) r_m1(2) = rk_m1(2)*s(indx(ino3_z)) r_m1(3) = rk_m1(3)*s(indx(ihono_z)) r_m1(4) = rk_m1(4)*s(indx(ihno3_z)) r_m1(5) = rk_m1(5)*s(indx(ihno4_z)) r_m1(6) = rk_m1(6)*s(indx(in2o5_z)) r_m1(7) = rk_m1(7)*s(indx(io3_z)) r_m1(8) = rk_m1(8)*s(indx(io3_z)) r_m1(9) = rk_m1(9)*s(indx(ih2o2_z)) r_m1(10) = rk_m1(10)*s(indx(io1d_z))*oxygen r_m1(11) = rk_m1(11)*s(indx(io1d_z))*nitrogen r_m1(12) = rk_m1(12)*s(indx(io1d_z))*h2o r_m1(13) = rk_m1(13)*s(indx(io3p_z))*oxygen r_m1(14) = rk_m1(14)*s(indx(io3p_z))*s(indx(io3_z)) r_m1(15) = rk_m1(15)*s(indx(io3p_z))*s(indx(ino2_z)) r_m1(16) = rk_m1(16)*s(indx(io3p_z))*s(indx(ino2_z)) r_m1(17) = rk_m1(17)*s(indx(io3p_z))*s(indx(ino_z)) r_m1(18) = rk_m1(18)*s(indx(io3_z))*s(indx(ino_z)) r_m1(19) = rk_m1(19)*s(indx(io3_z))*s(indx(ino2_z)) r_m1(20) = rk_m1(20)*s(indx(io3_z))*s(indx(ioh_z)) r_m1(21) = rk_m1(21)*s(indx(io3_z))*s(indx(iho2_z)) r_m1(22) = rk_m1(22)*s(indx(ioh_z))*hydrogen r_m1(23) = rk_m1(23)*s(indx(ioh_z))*s(indx(ino_z)) r_m1(24) = rk_m1(24)*s(indx(ioh_z))*s(indx(ino2_z)) r_m1(25) = rk_m1(25)*s(indx(ioh_z))*s(indx(ino3_z)) r_m1(26) = rk_m1(26)*s(indx(ioh_z))*s(indx(ihono_z)) r_m1(27) = rk_m1(27)*s(indx(ioh_z))*s(indx(ihno3_z)) r_m1(28) = rk_m1(28)*s(indx(ioh_z))*s(indx(ihno4_z)) r_m1(29) = rk_m1(29)*s(indx(ioh_z))*s(indx(iho2_z)) r_m1(30) = rk_m1(30)*s(indx(ioh_z))*s(indx(ih2o2_z)) r_m1(31) = rk_m1(31)*s(indx(iho2_z))*s(indx(iho2_z)) r_m1(32) = rk_m1(32)*s(indx(iho2_z))*s(indx(iho2_z))*h2o r_m1(33) = rk_m1(33)*s(indx(iho2_z))*s(indx(ino_z)) r_m1(34) = rk_m1(34)*s(indx(iho2_z))*s(indx(ino2_z)) r_m1(35) = rk_m1(35)*s(indx(iho2_z))*s(indx(ino2_z)) r_m1(36) = rk_m1(36)*s(indx(ihno4_z)) r_m1(37) = rk_m1(37)*s(indx(ino3_z))*s(indx(ino_z)) r_m1(38) = rk_m1(38)*s(indx(ino3_z))*s(indx(ino2_z)) r_m1(39) = rk_m1(39)*s(indx(ino3_z))*s(indx(ino2_z)) r_m1(40) = rk_m1(40)*s(indx(ino3_z))*s(indx(ino3_z)) r_m1(41) = rk_m1(41)*s(indx(ino3_z))*s(indx(iho2_z)) r_m1(42) = rk_m1(42)*s(indx(in2o5_z))*h2o r_m1(43) = rk_m1(43)*s(indx(in2o5_z)) r_m1(44) = rk_m1(44)*s(indx(ico_z))*s(indx(ioh_z)) r_m1(45) = rk_m1(45)*s(indx(iso2_z))*s(indx(ioh_z)) r_m1(46) = rk_m1(46)*s(indx(ioh_z))*ch4 r_m1(47) = rk_m1(47)*s(indx(ic2h6_z))*s(indx(ioh_z)) r_m1(48) = rk_m1(48)*s(indx(ich3oh_z))*s(indx(ioh_z)) r_m1(49) = rk_m1(49)*s(indx(ihcho_z)) r_m1(50) = rk_m1(50)*s(indx(ihcho_z)) r_m1(51) = rk_m1(51)*s(indx(ihcho_z))*s(indx(ioh_z)) r_m1(52) = rk_m1(52)*s(indx(ihcho_z))*s(indx(ino3_z)) r_m1(53) = rk_m1(53)*s(indx(ich3ooh_z)) r_m1(54) = rk_m1(54)*s(indx(iethooh_z)) r_m1(55) = rk_m1(55)*s(indx(ich3ooh_z))*s(indx(ioh_z)) r_m1(56) = rk_m1(56)*s(indx(iethooh_z))*s(indx(ioh_z)) r_m1(57) = rk_m1(57)*s(indx(ich3o2_z))*s(indx(ino_z)) r_m1(58) = rk_m1(58)*s(indx(iethp_z))*s(indx(ino_z)) r_m1(59) = rk_m1(59)*s(indx(ich3o2_z))*s(indx(ino3_z)) r_m1(60) = rk_m1(60)*s(indx(iethp_z))*s(indx(ino3_z)) r_m1(61) = rk_m1(61)*s(indx(ich3o2_z))*s(indx(iho2_z)) r_m1(62) = rk_m1(62)*s(indx(iethp_z))*s(indx(iho2_z)) r_m1(63) = rk_m1(63)*s(indx(ich3o2_z)) r_m1(64) = rk_m1(64)*s(indx(iethp_z)) r_m1(65) = rk_m1(65)*s(indx(ic2h5oh_z))*s(indx(ioh_z)) r_m2(2) = rk_m2(2)*s(indx(iald2_z)) r_m2(3) = rk_m2(3)*s(indx(iald2_z))*s(indx(ioh_z)) r_m2(4) = rk_m2(4)*s(indx(iald2_z))*s(indx(ino3_z)) r_m2(31) = rk_m2(31)*s(indx(ic2o3_z))*s(indx(ino2_z)) r_m2(32) = rk_m2(32)*s(indx(ipan_z)) r_m2(34) = rk_m2(34)*s(indx(ic2o3_z))*s(indx(ino_z)) r_m2(39) = rk_m2(39)*s(indx(ic2o3_z))*s(indx(ino3_z)) r_m2(44) = rk_m2(44)*s(indx(ic2o3_z))*s(indx(iho2_z)) r_m2(49) = rk_m2(49)*s(indx(ic2o3_z)) return end subroutine gasrate_m1 !*********************************************************************** ! <21.> subr gasrate_m2 ! ! purpose: computes reaction rates for mechanism 2 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- ! subroutine gasrate_m2( indx, s, r_m2, rk_m2 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real s(ngas_tot), r_m2(nrxn_m2), rk_m2(nrxn_m2) r_m2(1) = rk_m2(1)*s(indx(ipar_z))*s(indx(ioh_z)) r_m2(5) = rk_m2(5)*s(indx(iaone_z)) r_m2(6) = rk_m2(6)*s(indx(iaone_z))*s(indx(ioh_z)) r_m2(7) = rk_m2(7)*s(indx(imgly_z)) r_m2(8) = rk_m2(8)*s(indx(imgly_z))*s(indx(ioh_z)) r_m2(9) = rk_m2(9)*s(indx(imgly_z))*s(indx(ino3_z)) r_m2(10) = rk_m2(10)*s(indx(ieth_z))*s(indx(io3_z)) r_m2(11) = rk_m2(11)*s(indx(ieth_z))*s(indx(ioh_z)) r_m2(12) = rk_m2(12)*s(indx(iolet_z))*s(indx(io3_z)) r_m2(13) = rk_m2(13)*s(indx(iolei_z))*s(indx(io3_z)) r_m2(14) = rk_m2(14)*s(indx(iolet_z))*s(indx(ioh_z)) r_m2(15) = rk_m2(15)*s(indx(iolei_z))*s(indx(ioh_z)) r_m2(16) = rk_m2(16)*s(indx(iolet_z))*s(indx(ino3_z)) r_m2(17) = rk_m2(17)*s(indx(iolei_z))*s(indx(ino3_z)) r_m2(18) = rk_m2(18)*s(indx(itol_z))*s(indx(ioh_z)) r_m2(19) = rk_m2(19)*s(indx(ixyl_z))*s(indx(ioh_z)) r_m2(20) = rk_m2(20)*s(indx(ito2_z))*s(indx(ino_z)) r_m2(21) = rk_m2(21)*s(indx(icres_z))*s(indx(ioh_z)) r_m2(22) = rk_m2(22)*s(indx(icres_z))*s(indx(ino3_z)) r_m2(23) = rk_m2(23)*s(indx(icro_z))*s(indx(ino2_z)) r_m2(24) = rk_m2(24)*s(indx(iopen_z))*s(indx(ioh_z)) r_m2(25) = rk_m2(25)*s(indx(iopen_z)) r_m2(26) = rk_m2(26)*s(indx(iopen_z))*s(indx(io3_z)) r_m2(27) = rk_m2(27)*s(indx(irooh_z)) r_m2(28) = rk_m2(28)*s(indx(irooh_z))*s(indx(ioh_z)) r_m2(29) = rk_m2(29)*s(indx(ionit_z))*s(indx(ioh_z)) r_m2(30) = rk_m2(30)*s(indx(ionit_z)) r_m2(33) = rk_m2(33)*s(indx(iro2_z))*s(indx(ino_z)) r_m2(35) = rk_m2(35)*s(indx(iano2_z))*s(indx(ino_z)) r_m2(36) = rk_m2(36)*s(indx(inap_z))*s(indx(ino_z)) r_m2(37) = rk_m2(37)*s(indx(ixo2_z))*s(indx(ino_z)) r_m2(38) = rk_m2(38)*s(indx(iro2_z))*s(indx(ino3_z)) r_m2(40) = rk_m2(40)*s(indx(iano2_z))*s(indx(ino3_z)) r_m2(41) = rk_m2(41)*s(indx(inap_z))*s(indx(ino3_z)) r_m2(42) = rk_m2(42)*s(indx(ixo2_z))*s(indx(ino3_z)) r_m2(43) = rk_m2(43)*s(indx(iro2_z))*s(indx(iho2_z)) r_m2(45) = rk_m2(45)*s(indx(iano2_z))*s(indx(iho2_z)) r_m2(46) = rk_m2(46)*s(indx(inap_z))*s(indx(iho2_z)) r_m2(47) = rk_m2(47)*s(indx(ixo2_z))*s(indx(iho2_z)) r_m2(48) = rk_m2(48)*s(indx(iro2_z)) r_m2(50) = rk_m2(50)*s(indx(iano2_z)) r_m2(51) = rk_m2(51)*s(indx(inap_z)) r_m2(52) = rk_m2(52)*s(indx(ixo2_z)) r_m2(53) = rk_m2(53)*s(indx(ipar_z))*s(indx(ixpar_z)) ! return end subroutine gasrate_m2 !*********************************************************************** ! <22.> subr gasrate_m3 ! ! purpose: computes reaction rates for mechanism 3 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- ! subroutine gasrate_m3( indx, s, r_m3, rk_m3 ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real s(ngas_tot), r_m3(nrxn_m3), rk_m3(nrxn_m3) r_m3(1) = rk_m3(1)*s(indx(iisop_z))*s(indx(ioh_z)) r_m3(2) = rk_m3(2)*s(indx(iisop_z))*s(indx(io3_z)) r_m3(3) = rk_m3(3)*s(indx(iisop_z))*s(indx(ino3_z)) r_m3(4) = rk_m3(4)*s(indx(iisoprd_z)) r_m3(5) = rk_m3(5)*s(indx(iisoprd_z))*s(indx(ioh_z)) r_m3(6) = rk_m3(6)*s(indx(iisoprd_z))*s(indx(io3_z)) r_m3(7) = rk_m3(7)*s(indx(iisoprd_z))*s(indx(ino3_z)) r_m3(8) = rk_m3(8)*s(indx(iisopp_z))*s(indx(ino_z)) r_m3(9) = rk_m3(9)*s(indx(iisopn_z))*s(indx(ino_z)) r_m3(10) = rk_m3(10)*s(indx(iisopo2_z))*s(indx(ino_z)) r_m3(11) = rk_m3(11)*s(indx(iisopp_z))*s(indx(iho2_z)) r_m3(12) = rk_m3(12)*s(indx(iisopn_z))*s(indx(iho2_z)) r_m3(13) = rk_m3(13)*s(indx(iisopo2_z))*s(indx(iho2_z)) r_m3(14) = rk_m3(14)*s(indx(iisopp_z)) r_m3(15) = rk_m3(15)*s(indx(iisopn_z)) r_m3(16) = rk_m3(16)*s(indx(iisopo2_z)) return end subroutine gasrate_m3 !*********************************************************************** ! <23.> subr gasrate_m4 ! ! purpose: computes reaction rates for mechanism 4 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- ! subroutine gasrate_m4( indx, s, r_m4, rk_m4, oxygen ) use module_data_cbmz implicit none ! subr arguments integer indx(ngas_z) real s(ngas_tot), r_m4(nrxn_m4), rk_m4(nrxn_m4) real oxygen r_m4(1) = rk_m4(1)*s(indx(idms_z))*s(indx(ioh_z)) r_m4(2) = rk_m4(2)*s(indx(idms_z))*s(indx(ino3_z)) r_m4(3) = rk_m4(3)*s(indx(idms_z))*s(indx(io3p_z)) r_m4(4) = rk_m4(4)*s(indx(idms_z))*s(indx(ioh_z)) r_m4(5) = rk_m4(5)*s(indx(ich3sch2oo_z))*s(indx(ino_z)) r_m4(6) = rk_m4(6)*s(indx(ich3sch2oo_z))*s(indx(ich3o2_z)) r_m4(7) = rk_m4(7)*s(indx(ich3sch2oo_z))*s(indx(ich3so2_z)) r_m4(8) = rk_m4(8)*s(indx(ich3sch2oo_z))*s(indx(ich3sch2oo_z)) r_m4(9) = rk_m4(9)*s(indx(idmso_z))*s(indx(ioh_z)) r_m4(10) = rk_m4(10)*s(indx(idmso2_z))*s(indx(ioh_z)) r_m4(11) = rk_m4(11)*s(indx(ich3so2ch2oo_z))*s(indx(ino_z)) r_m4(12) = rk_m4(12)*s(indx(ich3so2ch2oo_z))*s(indx(ich3o2_z)) r_m4(13) = rk_m4(13)*s(indx(ich3so2h_z))*s(indx(iho2_z)) r_m4(14) = rk_m4(14)*s(indx(ich3so2h_z))*s(indx(ino3_z)) r_m4(15) = rk_m4(15)*s(indx(ich3so2h_z))*s(indx(ich3o2_z)) r_m4(16) = rk_m4(16)*s(indx(ich3so2h_z))*s(indx(ioh_z)) r_m4(17) = rk_m4(17)*s(indx(ich3so2h_z))*s(indx(ich3so3_z)) r_m4(18) = rk_m4(18)*s(indx(ich3so2_z)) r_m4(19) = rk_m4(19)*s(indx(ich3so2_z))*s(indx(ino2_z)) r_m4(20) = rk_m4(20)*s(indx(ich3so2_z))*s(indx(io3_z)) r_m4(21) = rk_m4(21)*s(indx(ich3so2_z))*s(indx(iho2_z)) r_m4(22) = rk_m4(22)*s(indx(ich3so2_z))*s(indx(ich3o2_z)) r_m4(23) = rk_m4(23)*s(indx(ich3so2_z))*s(indx(ioh_z)) r_m4(24) = rk_m4(24)*s(indx(ich3so2_z))*oxygen r_m4(25) = rk_m4(25)*s(indx(ich3so2oo_z)) r_m4(26) = rk_m4(26)*s(indx(ich3so2oo_z))*s(indx(ino_z)) r_m4(27) = rk_m4(27)*s(indx(ich3so2oo_z))*s(indx(ich3o2_z)) r_m4(28) = rk_m4(28)*s(indx(ich3so3_z)) r_m4(29) = rk_m4(29)*s(indx(ich3so3_z))*s(indx(ino2_z)) r_m4(30) = rk_m4(30)*s(indx(ich3so3_z))*s(indx(ino_z)) r_m4(31) = rk_m4(31)*s(indx(ich3so3_z))*s(indx(iho2_z)) r_m4(32) = rk_m4(32)*s(indx(ich3so3_z))*s(indx(ihcho_z)) return end subroutine gasrate_m4 !************************************************************************** ! <24.> subr loadperoxyparameters ! ! purpose: loads thermal rate coefficients for peroxy-peroxy ! permutation reactions ! ! author : Rahul A. Zaveri ! date : June 1998 ! ! nomenclature: ! Aperox = Pre-exponential factor (molec-cc-s) ! Bperox = activation energy (-E/R) (K) ! !------------------------------------------------------------------------- subroutine loadperoxyparameters( Aperox, Bperox ) use module_data_cbmz implicit none ! subr arguments real Aperox(nperox,nperox), Bperox(nperox,nperox) ! local variables integer i, j Aperox(jch3o2,jch3o2) = 2.5e-13 Aperox(jethp,jethp) = 6.8e-14 Aperox(jc2o3,jc2o3) = 2.9e-12 Aperox(jano2,jano2) = 8.0e-12 Aperox(jnap,jnap) = 1.0e-12 Aperox(jro2,jro2) = 5.3e-16 Aperox(jisopp,jisopp) = 3.1e-14 Aperox(jisopn,jisopn) = 3.1e-14 Aperox(jisopo2,jisopo2) = 3.1e-14 Aperox(jxo2,jxo2) = 3.1e-14 Bperox(jch3o2,jch3o2) = 190. Bperox(jethp,jethp) = 0.0 Bperox(jc2o3,jc2o3) = 500. Bperox(jano2,jano2) = 0.0 Bperox(jnap,jnap) = 0.0 Bperox(jro2,jro2) = 1980. Bperox(jisopp,jisopp) = 1000. Bperox(jisopn,jisopn) = 1000. Bperox(jisopo2,jisopo2) = 1000. Bperox(jxo2,jxo2) = 1000. do i = 1, nperox do j = 1, nperox if(i.ne.j)then Aperox(i,j) = 2.0*sqrt(Aperox(i,i)*Aperox(j,j)) Bperox(i,j) = 0.5*(Bperox(i,i) + Bperox(j,j)) endif enddo enddo ! except for Aperox(jc2o3,jch3o2) = 1.3e-12 Aperox(jch3o2,jc2o3) = 1.3e-12 Bperox(jc2o3,jch3o2) = 640. Bperox(jch3o2,jc2o3) = 640. return end subroutine loadperoxyparameters !************************************************************************** ! <25.> subr peroxyrateconstants ! ! purpose: computes parameterized thermal rate coefficients ! for the alkylperoxy radical permutation reactions ! for the entire mechanism. ! ! author : Rahul A. Zaveri ! date : June 1998 ! ! nomenclature: ! rk_param = parameterized reaction rate constants (1/s) ! rk_perox = individual permutation reaction rate constants (molec-cc-s) ! te = ambient atmospheric temperature (K) ! !------------------------------------------------------------------------- subroutine peroxyrateconstants( tempbox, cbox, & Aperox, Bperox, rk_param ) use module_data_cbmz implicit none ! subr arguments real tempbox, cbox(ngas_z) real Aperox(nperox,nperox), Bperox(nperox,nperox), rk_param(nperox) ! local variables integer i, j real te real sperox(nperox), rk_perox(nperox,nperox) te = tempbox sperox(jch3o2) = cbox(ich3o2_z) sperox(jethp) = cbox(iethp_z) sperox(jro2) = cbox(iro2_z) sperox(jc2o3) = cbox(ic2o3_z) sperox(jano2) = cbox(iano2_z) sperox(jnap) = cbox(inap_z) sperox(jisopp) = cbox(iisopp_z) sperox(jisopn) = cbox(iisopn_z) sperox(jisopo2) = cbox(iisopo2_z) sperox(jxo2) = cbox(ixo2_z) ! ! initialize to zero do i = 1, nperox rk_param(i) = 0.0 enddo do i = 1, nperox do j = 1, nperox rk_perox(i,j) = arr( Aperox(i,j), Bperox(i,j), te ) rk_param(i) = rk_param(i) + rk_perox(i,j)*sperox(j) enddo enddo return end subroutine peroxyrateconstants !*********************************************************************** ! <26.> subr gasthermrk_m1 ! ! purpose: computes thermal reaction rate coefficients for ! mechanism 1 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- subroutine gasthermrk_m1( tempbox, cair_mlc, & rk_photo, rk_param, rk_m1, rk_m2 ) use module_data_cbmz implicit none ! subr arguments real tempbox, cair_mlc real rk_photo(nphoto), rk_param(nperox) real rk_m1(nrxn_m1), rk_m2(nrxn_m2) ! local variables integer i real rk0, rk2, rk3, rki, rko, rmm, rnn, te ! real arr, troe te = tempbox rk_m1(1) = rk_photo(jphoto_no2) rk_m1(2) = rk_photo(jphoto_no3) rk_m1(3) = rk_photo(jphoto_hono) rk_m1(4) = rk_photo(jphoto_hno3) rk_m1(5) = rk_photo(jphoto_hno4) rk_m1(6) = rk_photo(jphoto_n2o5) rk_m1(7) = rk_photo(jphoto_o3a) rk_m1(8) = rk_photo(jphoto_o3b) rk_m1(9) = rk_photo(jphoto_h2o2) rk_m1(10) = arr(3.2e-11, 70., te) rk_m1(11) = arr(1.8e-11, 110., te) rk_m1(12) = 2.2e-10 rk_m1(13) = cair_mlc*6.e-34*(te/300.)**(-2.3) rk_m1(14) = arr(8.0e-12, -2060., te) rk_m1(15) = arr(6.5e-12, -120., te) ! rk0 = 9.0e-32 rnn = 2.0 rki = 2.2e-11 rmm = 0.0 rk_m1(16) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk0 = 9.0e-32 rnn = 1.5 rki = 3.0e-11 rmm = 0.0 rk_m1(17) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m1(18) = arr(2.0e-12, -1400., te) rk_m1(19) = arr(1.2e-13, -2450., te) rk_m1(20) = arr(1.6e-12, -940., te) rk_m1(21) = arr(1.1e-14, -500., te) rk_m1(22) = arr(5.5e-12, -2000., te) ! rk0 = 7.0e-31 rnn = 2.6 rki = 3.6e-11 rmm = 0.1 rk_m1(23) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk0 = 2.5e-30 rnn = 4.4 rki = 1.6e-11 rmm = 1.7 rk_m1(24) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m1(25) = 2.2e-11 rk_m1(26) = arr(1.8e-11, -390., te) rko = 7.2e-15 * exp(785./te) rk2 = 4.1e-16 * exp(1440./te) rk3 = 1.9e-33 * exp(725./te)*cair_mlc rk_m1(27) = rko + rk3/(1.+rk3/rk2) rk_m1(28) = arr(1.3e-12, 380., te) rk_m1(29) = arr(4.8e-11, 250., te) rk_m1(30) = arr(2.9e-12, -160., te) rk_m1(31) = 2.3e-13 * exp(600./te) + & 1.7e-33 * exp(1000./te)*cair_mlc ! ho2 + ho2 --> h2o2 rk_m1(32) = rk_m1(31)*1.4e-21*exp(2200./te) ! ho2 + ho2 + h2o --> h2o2 rk_m1(33) = arr(3.5e-12, 250., te) ! rk0 = 1.8e-31 rnn = 3.2 rki = 4.7e-12 rmm = 1.4 rk_m1(34) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m1(35) = 5.0e-16 rk_m1(36) = rk_m1(34)*arr(4.8e26, -10900., te) rk_m1(37) = arr(1.5e-11, 170., te) rk_m1(38) = arr(4.5e-14, -1260., te) ! rk0 = 2.2e-30 rnn = 3.9 rki = 1.5e-12 rmm = 0.7 rk_m1(39) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m1(40) = arr(8.5e-13, -2450., te) rk_m1(41) = 3.5e-12 rk_m1(42) = 2.0e-21 rk_m1(43) = rk_m1(39)*arr(3.7e26, -11000., te) rk_m1(44) = 1.5e-13 * (1.+8.18e-23*te*cair_mlc) ! co + oh --> ho2 rk0 = 3.0e-31 rnn = 3.3 rki = 1.5e-12 rmm = 0.0 rk_m1(45) = Troe(cair_mlc,te,rk0,rnn,rki,rmm) rk_m1(46) = te**.667*arr(2.8e-14, -1575., te) rk_m1(47) = te**2*arr(1.5e-17, -492., te) rk_m1(48) = arr(6.7e-12, -600., te) rk_m1(49) = rk_photo(jphoto_hchoa) ! hcho + hv --> 2ho2 + co rk_m1(50) = rk_photo(jphoto_hchob) ! hcho + hv --> co rk_m1(51) = 1.0e-11 rk_m1(52) = arr(3.4e-13, -1900., te) rk_m1(53) = rk_photo(jphoto_ch3ooh) rk_m1(54) = rk_photo(jphoto_ethooh) rk_m1(55) = arr(3.8e-12, 200., te) rk_m1(56) = arr(3.8e-12, 200., te) rk_m1(57) = arr(3.0e-12, 280., te) rk_m1(58) = arr(2.6e-12, 365., te) rk_m1(59) = 1.1e-12 rk_m1(60) = 2.5e-12 rk_m1(61) = arr(3.8e-13, 800., te) rk_m1(62) = arr(7.5e-13, 700., te) rk_m1(63) = rk_param(jch3o2) rk_m1(64) = rk_param(jethp) rk_m1(65) = arr(7.0e-12, -235.,te) rk_m2(2) = rk_photo(jphoto_ald2) rk_m2(3) = arr(5.6e-12, 270., te) rk_m2(4) = arr(1.4e-12, -1900., te) ! rk0 = 9.7e-29 rnn = 5.6 rki = 9.3e-12 rmm = 1.5 rk_m2(31) = troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m2(32) = rk_m2(31)*arr(1.1e28, -14000., te) rk_m2(34) = arr(5.3e-12, 360., te) rk_m2(39) = 4.0e-12 rk_m2(44) = arr(4.5e-13, 1000., te) rk_m2(49) = rk_param(jc2o3) ! Heterogeneous reactions ! rk_m1(65) = rk_het(1) ! O3 --> ! rk_m1(66) = rk_het(2) ! HO2 --> 0.5H2O2 ! rk_m1(67) = rk_het(3) ! NO2 --> 0.5HONO + 0.5HNO3 ! rk_m1(68) = rk_het(4) ! N2O5 --> 2HNO3 ! rk_m1(69) = rk_het(5) ! HNO3 --> NO2 ! rk_m1(70) = rk_het(6) ! HNO3 --> NO ! rk_m1(71) = rk_het(7) ! NO3 --> NO + O2 !!$Turn off this check since initializing to 0's; wig 16-Nov-2007 !!$! all rate constants but be >= 0 !!$ do i = 1, nrxn_m1 !!$ rk_m1(i) = max( rk_m1(i), 0.0 ) !!$ end do !!$ do i = 1, nrxn_m2 !!$ rk_m2(i) = max( rk_m2(i), 0.0 ) !!$ end do return end subroutine gasthermrk_m1 !*********************************************************************** ! <27.> subr gasthermrk_m2 ! ! purpose: computes thermal reaction rate coefficients for ! mechanism 2 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- subroutine gasthermrk_m2( tempbox, cair_mlc, & rk_photo, rk_param, rk_m2 ) use module_data_cbmz implicit none ! subr arguments real tempbox, cair_mlc real rk_photo(nphoto), rk_param(nperox), rk_m2(nrxn_m2) ! local variables integer i real rk0, rki, rmm, rnn, te ! real arr, troe te = tempbox rk_m2(1) = 8.1e-13 rk_m2(5) = rk_photo(jphoto_aone) rk_m2(6) = te**2*arr(5.3e-18, -230., te) rk_m2(7) = rk_photo(jphoto_mgly) rk_m2(8) = 1.7e-11 rk_m2(9) = arr(1.4e-12, -1900., te) rk_m2(10) = arr(1.2e-14, -2630., te) ! rk0 = 1.0e-28 rnn = 0.8 rki = 8.8e-12 rmm = 0.0 rk_m2(11) = troe(cair_mlc,te,rk0,rnn,rki,rmm) ! rk_m2(12) = arr(4.2e-15, -1800., te) rk_m2(13) = arr(8.9e-16, -392., te) rk_m2(14) = arr(5.8e-12, 478., te) rk_m2(15) = arr(2.9e-11, 255., te) rk_m2(16) = arr(3.1e-13, -1010., te) rk_m2(17) = 2.5e-12 rk_m2(18) = arr(2.1e-12, 322., te) rk_m2(19) = arr(1.7e-11, 116., te) rk_m2(20) = 8.1e-12 rk_m2(21) = 4.1e-11 rk_m2(22) = 2.2e-11 rk_m2(23) = 1.4e-11 rk_m2(24) = 3.0e-11 rk_m2(25) = rk_photo(jphoto_open) rk_m2(26) = arr(5.4e-17, -500., te) rk_m2(27) = rk_photo(jphoto_rooh) rk_m2(28) = arr(3.8e-12, 200., te) rk_m2(29) = arr(1.6e-11, -540., te) rk_m2(30) = rk_photo(jphoto_onit) rk_m2(33) = 4.0e-12 rk_m2(35) = 4.0e-12 rk_m2(36) = 4.0e-12 rk_m2(37) = 4.0e-12 rk_m2(38) = 2.5e-12 rk_m2(40) = 1.2e-12 rk_m2(41) = 4.0e-12 rk_m2(42) = 2.5e-12 rk_m2(43) = arr(1.7e-13, 1300., te) rk_m2(45) = arr(1.2e-13, 1300., te) rk_m2(46) = arr(1.7e-13, 1300., te) rk_m2(47) = arr(1.7e-13, 1300., te) rk_m2(48) = rk_param(jro2) rk_m2(50) = rk_param(jano2) rk_m2(51) = rk_param(jnap) rk_m2(52) = rk_param(jxo2) rk_m2(53) = 1.0e-11 ! XPAR + PAR --> !!$Turn off this check since initializing to 0's; wig 16-Nov-2007 !!$! all rate constants but be >= 0 !!$ do i = 1, nrxn_m2 !!$ rk_m2(i) = max( rk_m2(i), 0.0 ) !!$ end do return end subroutine gasthermrk_m2 !*********************************************************************** ! <28.> subr gasthermrk_m3 ! ! purpose: computes thermal reaction rate coefficients for ! mechanism 3 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- subroutine gasthermrk_m3( tempbox, cair_mlc, & rk_photo, rk_param, rk_m3 ) use module_data_cbmz implicit none ! subr arguments real tempbox, cair_mlc real rk_photo(nphoto), rk_param(nperox), rk_m3(nrxn_m3) ! local variables integer i real te ! real arr te = tempbox ! rk_m3(1) = arr(2.6e-11, 409., te) rk_m3(2) = arr(1.2e-14, -2013., te) rk_m3(3) = arr(3.0e-12, -446., te) rk_m3(4) = rk_photo(jphoto_isoprd) rk_m3(5) = 3.3e-11 rk_m3(6) = 7.0e-18 rk_m3(7) = 1.0e-15 rk_m3(8) = 4.0e-12 rk_m3(9) = 4.0e-12 rk_m3(10) = 4.0e-12 rk_m3(11) = arr(1.7e-13, 1300., te) rk_m3(12) = arr(1.7e-13, 1300., te) rk_m3(13) = arr(1.7e-13, 1300., te) rk_m3(14) = rk_param(jisopp) rk_m3(15) = rk_param(jisopn) rk_m3(16) = rk_param(jisopo2) ! all rate constants but be >= 0 do i = 1, nrxn_m3 rk_m3(i) = max( rk_m3(i), 0.0 ) end do return end subroutine gasthermrk_m3 !*********************************************************************** ! <29.> subr gasthermrk_m4 ! ! purpose: computes thermal reaction rate coefficients for ! mechanism 4 ! ! author : Rahul A. Zaveri ! date : December 1998 ! !------------------------------------------------------------------------- subroutine gasthermrk_m4( tempbox, cair_mlc, & rk_photo, rk_param, rk_m4 ) use module_data_cbmz implicit none ! subr arguments real tempbox, cair_mlc real rk_photo(nphoto), rk_param(nperox), rk_m4(nrxn_m4) ! local variables integer i real B_abs, B_add, rk_tot, rk_tot_den, rk_tot_num, te ! real arr te = tempbox ! rk_m4(1) = arr(9.6e-12, -234., te) ! ch3sch3 + oh --> ch3sch2 rk_m4(2) = arr(1.4e-13, 500., te) rk_m4(3) = arr(1.3e-11, 409., te) ! Hynes et al. (1986) rk_tot_num = te * exp(-234./te) + & 8.46e-10 * exp(7230./te) + & 2.68e-10 * exp(7810./te) rk_tot_den = 1.04e+11 * te + 88.1 * exp(7460./te) rk_tot = rk_tot_num/rk_tot_den B_abs = rk_m4(1)/rk_tot B_add = 1. - B_abs rk_m4(4) = B_add*rk_tot ! ch3sch3 + oh --> ch3s(oh)ch3 rk_m4(5) = 8.0e-12 rk_m4(6) = 1.8e-13 rk_m4(7) = 2.5e-13 rk_m4(8) = 8.6e-14 rk_m4(9) = 5.8e-11 rk_m4(10) = 1.0e-14 rk_m4(11) = 5.0e-12 rk_m4(12) = 1.8e-13 rk_m4(13) = 1.0e-15 rk_m4(14) = 1.0e-13 rk_m4(15) = 1.0e-15 rk_m4(16) = 1.6e-11 rk_m4(17) = 1.0e-13 !Bug fixing, the coefficient changed from 2.5e-13 to 2.5e13 by Qing.Yang@pnl.gov rk_m4(18) = arr(2.5e13, -8686., te) rk_m4(19) = 1.0e-14 rk_m4(20) = 5.0e-15 rk_m4(21) = 2.5e-13 rk_m4(22) = 2.5e-13 rk_m4(23) = 5.0e-11 rk_m4(24) = 2.6e-18 rk_m4(25) = 3.3 rk_m4(26) = 1.0e-11 rk_m4(27) = 5.5e-12 rk_m4(28) = arr(2.0e17, -12626., te) rk_m4(29) = 3.0e-15 rk_m4(30) = 3.0e-15 rk_m4(31) = 5.0e-11 rk_m4(32) = 1.6e-15 ! all rate constants but be >= 0 do i = 1, nrxn_m4 rk_m4(i) = max( rk_m4(i), 0.0 ) end do return end subroutine gasthermrk_m4 !*********************************************************************** ! <26.> subr hetrateconstants ! ! purpose: computes heterogeneous reaction rate coefficients ! ! author : Rahul A. Zaveri ! date : May 2000 ! !------------------------------------------------------------------------- subroutine hetrateconstants implicit none return end subroutine hetrateconstants !*********************************************************************** ! <31.> func troe ! ! purpose: calculates Troe reaction rate coefficient ! ! author : Rahul A. Zaveri ! date : December 1998 !----------------------------------------------------------------------- real function troe( cairmlc, te, rk0, rnn, rki, rmm ) implicit none ! func parameters real cairmlc, te, rk0, rnn, rki, rmm ! local variables real expo rk0 = rk0*cairmlc*(te/300.)**(-rnn) rki = rki*(te/300.)**(-rmm) expo= 1./(1. + (ALOG10(rk0/rki))**2) troe = (rk0*rki/(rk0+rki))*.6**expo return end function troe !*********************************************************************** ! <32.> func arr ! ! purpose: calculates arrhenius rate coefficient ! ! author : Rahul A. Zaveri ! date : December 1998 !----------------------------------------------------------------------- real function arr( aa, bb, te ) implicit none ! func parameters real aa, bb, te arr = aa*exp(bb/te) return end function arr !*********************************************************************** ! subr mapgas_tofrom_host ! ! purpose: maps gas species between cboxold/new and host arrays ! ! author : R. C. Easter ! date : November, 2003 ! ! ---------------------------------------------------------------------- subroutine mapgas_tofrom_host( imap, & i_boxtest_units_convert, & it,jt,kt, ims,ime, jms,jme, kms,kme, & num_moist, num_chem, moist, chem, & t_phy, p_phy, rho_phy, & cbox, tempbox, pressbox, airdenbox, & cair_mlc, & h2o, ch4, oxygen, nitrogen, hydrogen ) use module_configure, only: & p_qv, & p_so2, p_sulf, p_no2, p_no, p_o3, & p_hno3, p_h2o2, p_ald, p_hcho, p_op1, & p_op2, p_paa, p_ora1, p_ora2, p_nh3, & p_n2o5, p_no3, p_pan, p_hc3, p_hc5, & p_hc8, p_eth, p_co, p_ol2, p_olt, & p_oli, p_tol, p_xyl, p_aco3, p_tpan, & p_hono, p_hno4, p_ket, p_gly, p_mgly, & p_dcb, p_onit, p_csl, p_iso, p_ho, & p_ho2, & p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh, & p_par, p_to2, p_cro, p_open, p_op3, & p_c2o3, p_ro2, p_ano2, p_nap, p_xo2, & p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2, & p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h, & p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, & p_mtf use module_data_cbmz implicit none ! subr arguments INTEGER, INTENT(IN) :: imap, it,jt,kt, ims,ime, jms,jme, kms,kme, & num_moist, num_chem, i_boxtest_units_convert REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & INTENT(IN) :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT) :: chem REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN) :: t_phy, & ! temperature p_phy, & ! air pressure (Pa) rho_phy ! air density (kg/m3) REAL, INTENT(INOUT) :: cbox(ngas_z) REAL, INTENT(INOUT) :: tempbox, pressbox, airdenbox REAL, INTENT(INOUT) :: cair_mlc REAL, INTENT(INOUT) :: h2o, ch4, oxygen, nitrogen, hydrogen ! local variables integer l real factoraa real, parameter :: eps=0.622 tempbox = t_phy(it,kt,jt) ! p_phy = (Pa); pressbox = (dynes/cm2) pressbox = p_phy(it,kt,jt)*10.0 ! rho_phy = (kg_air/m3); airdenbox = (mole_air/cm3) airdenbox = rho_phy(it,kt,jt)/28.966e3 if (i_boxtest_units_convert .eq. 10) then airdenbox = rho_phy(it,kt,jt) end if if (imap .gt. 0) goto 2000 ! ! imap==0 -- initial species mapping from host array to cboxold ! chem --> czz --> cbox ! ! note: do not map nh3, hcl ! cbox(:) = 0.0 ! cair_mlc = (molecules_air/cm3) cair_mlc = airdenbox*avognumkpp ! moist = (kg_h2o/kg_air); czz = (mole_h2o/cm3); h2o = (molecules_h2o/cm3) ! czz(ih2o_z) = (moist(it,kt,jt,p_qv)/eps)*airdenbox ! if (i_boxtest_units_convert .eq. 10) then ! czz(ih2o_z) = moist(it,kt,jt,p_qv)*airdenbox ! end if ! h2o = czz(ih2o_z)*avognumkpp h2o = (moist(it,kt,jt,p_qv)/eps)*airdenbox if (i_boxtest_units_convert .eq. 10) then h2o = moist(it,kt,jt,p_qv)*airdenbox end if h2o = h2o*avognumkpp ! czz(ich4_z) = 1.7e-6*airdenbox ! ch4 conc. in mol/cc ! ch4 = czz(ich4_z)*avognumkpp ! ch4 conc. in molec/cc ch4 = 1.7e-6*airdenbox*avognumkpp ! ch4 conc. in molec/cc oxygen = 0.21*cair_mlc ! o2 conc. in molec/cc nitrogen = 0.79*cair_mlc ! n2 conc. in molec/cc hydrogen = 0.58e-6*cair_mlc ! h2 conc. in molec/cc ! chem units = (ppm); czz units = (mole/cm3); cbox units = (molecules/cm3) factoraa = airdenbox*1.0e-6 if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox factoraa = factoraa*avognumkpp cbox(iso2_z) = chem(it,kt,jt,p_so2)*factoraa cbox(ih2so4_z) = chem(it,kt,jt,p_sulf)*factoraa cbox(ino2_z) = chem(it,kt,jt,p_no2)*factoraa cbox(ino_z) = chem(it,kt,jt,p_no)*factoraa cbox(io3_z) = chem(it,kt,jt,p_o3)*factoraa cbox(ihno3_z) = chem(it,kt,jt,p_hno3)*factoraa cbox(ih2o2_z) = chem(it,kt,jt,p_h2o2)*factoraa cbox(iald2_z) = chem(it,kt,jt,p_ald)*factoraa cbox(ihcho_z) = chem(it,kt,jt,p_hcho)*factoraa cbox(ich3ooh_z) = chem(it,kt,jt,p_op1)*factoraa cbox(iethooh_z) = chem(it,kt,jt,p_op2)*factoraa cbox(ihcooh_z) = chem(it,kt,jt,p_ora1)*factoraa cbox(ircooh_z) = chem(it,kt,jt,p_ora2)*factoraa cbox(inh3_z) = chem(it,kt,jt,p_nh3)*factoraa cbox(in2o5_z) = chem(it,kt,jt,p_n2o5)*factoraa cbox(ino3_z) = chem(it,kt,jt,p_no3)*factoraa cbox(ipan_z) = chem(it,kt,jt,p_pan)*factoraa cbox(ic2h6_z) = chem(it,kt,jt,p_eth)*factoraa cbox(ico_z) = chem(it,kt,jt,p_co)*factoraa cbox(ieth_z) = chem(it,kt,jt,p_ol2)*factoraa cbox(iolet_z) = chem(it,kt,jt,p_olt)*factoraa cbox(iolei_z) = chem(it,kt,jt,p_oli)*factoraa cbox(itol_z) = chem(it,kt,jt,p_tol)*factoraa cbox(ixyl_z) = chem(it,kt,jt,p_xyl)*factoraa cbox(ihono_z) = chem(it,kt,jt,p_hono)*factoraa cbox(ihno4_z) = chem(it,kt,jt,p_hno4)*factoraa cbox(iaone_z) = chem(it,kt,jt,p_ket)*factoraa cbox(imgly_z) = chem(it,kt,jt,p_mgly)*factoraa cbox(ionit_z) = chem(it,kt,jt,p_onit)*factoraa cbox(icres_z) = chem(it,kt,jt,p_csl)*factoraa cbox(iisop_z) = chem(it,kt,jt,p_iso)*factoraa cbox(ioh_z) = chem(it,kt,jt,p_ho)*factoraa cbox(iho2_z) = chem(it,kt,jt,p_ho2)*factoraa cbox(ihcl_z) = chem(it,kt,jt,p_hcl)*factoraa cbox(ich3o2_z) = chem(it,kt,jt,p_ch3o2)*factoraa cbox(iethp_z) = chem(it,kt,jt,p_ethp)*factoraa cbox(ich3oh_z) = chem(it,kt,jt,p_ch3oh)*factoraa cbox(ic2h5oh_z) = chem(it,kt,jt,p_c2h5oh)*factoraa cbox(ipar_z) = chem(it,kt,jt,p_par)*factoraa cbox(ito2_z) = chem(it,kt,jt,p_to2)*factoraa cbox(icro_z) = chem(it,kt,jt,p_cro)*factoraa cbox(iopen_z) = chem(it,kt,jt,p_open)*factoraa cbox(irooh_z) = chem(it,kt,jt,p_op3)*factoraa cbox(ic2o3_z) = chem(it,kt,jt,p_c2o3)*factoraa cbox(iro2_z) = chem(it,kt,jt,p_ro2)*factoraa cbox(iano2_z) = chem(it,kt,jt,p_ano2)*factoraa cbox(inap_z) = chem(it,kt,jt,p_nap)*factoraa cbox(ixo2_z) = chem(it,kt,jt,p_xo2)*factoraa cbox(ixpar_z) = chem(it,kt,jt,p_xpar)*factoraa cbox(iisoprd_z) = chem(it,kt,jt,p_isoprd)*factoraa cbox(iisopp_z) = chem(it,kt,jt,p_isopp)*factoraa cbox(iisopn_z) = chem(it,kt,jt,p_isopn)*factoraa cbox(iisopo2_z) = chem(it,kt,jt,p_isopo2)*factoraa cbox(idms_z) = chem(it,kt,jt,p_dms)*factoraa cbox(imsa_z) = chem(it,kt,jt,p_msa)*factoraa cbox(idmso_z) = chem(it,kt,jt,p_dmso)*factoraa cbox(idmso2_z) = chem(it,kt,jt,p_dmso2)*factoraa cbox(ich3so2h_z) = chem(it,kt,jt,p_ch3so2h)*factoraa cbox(ich3sch2oo_z) = chem(it,kt,jt,p_ch3sch2oo)*factoraa cbox(ich3so2_z) = chem(it,kt,jt,p_ch3so2)*factoraa cbox(ich3so3_z) = chem(it,kt,jt,p_ch3so3)*factoraa cbox(ich3so2oo_z) = chem(it,kt,jt,p_ch3so2oo)*factoraa cbox(ich3so2ch2oo_z) = chem(it,kt,jt,p_ch3so2ch2oo)*factoraa cbox(imtf_z) = chem(it,kt,jt,p_mtf)*factoraa cbox(ih2o_z) = h2o cbox(ich4_z) = ch4 cbox(io2_z) = oxygen cbox(in2_z) = nitrogen cbox(ih2_z) = hydrogen return ! ! imap==1 -- final species mapping from cbox back to host array ! cbox --> czz --> chem ! ! note1: do not map nh3, hcl, ch4 ! 2000 continue ! chem = (ppm); czz = (mole/cm3); cbox = (molecules/cm3) factoraa = airdenbox*1.0e-6 if (i_boxtest_units_convert .eq. 10) factoraa = airdenbox factoraa = factoraa*avognumkpp chem(it,kt,jt,p_so2) = cbox(iso2_z)/factoraa chem(it,kt,jt,p_sulf) = cbox(ih2so4_z)/factoraa chem(it,kt,jt,p_no2) = cbox(ino2_z)/factoraa chem(it,kt,jt,p_no) = cbox(ino_z)/factoraa chem(it,kt,jt,p_o3) = cbox(io3_z)/factoraa chem(it,kt,jt,p_hno3) = cbox(ihno3_z)/factoraa chem(it,kt,jt,p_h2o2) = cbox(ih2o2_z)/factoraa chem(it,kt,jt,p_ald) = cbox(iald2_z)/factoraa chem(it,kt,jt,p_hcho) = cbox(ihcho_z)/factoraa chem(it,kt,jt,p_op1) = cbox(ich3ooh_z)/factoraa chem(it,kt,jt,p_op2) = cbox(iethooh_z)/factoraa chem(it,kt,jt,p_ora1) = cbox(ihcooh_z)/factoraa chem(it,kt,jt,p_ora2) = cbox(ircooh_z)/factoraa chem(it,kt,jt,p_nh3) = cbox(inh3_z)/factoraa chem(it,kt,jt,p_n2o5) = cbox(in2o5_z)/factoraa chem(it,kt,jt,p_no3) = cbox(ino3_z)/factoraa chem(it,kt,jt,p_pan) = cbox(ipan_z)/factoraa chem(it,kt,jt,p_eth) = cbox(ic2h6_z)/factoraa chem(it,kt,jt,p_co) = cbox(ico_z)/factoraa chem(it,kt,jt,p_ol2) = cbox(ieth_z)/factoraa chem(it,kt,jt,p_olt) = cbox(iolet_z)/factoraa chem(it,kt,jt,p_oli) = cbox(iolei_z)/factoraa chem(it,kt,jt,p_tol) = cbox(itol_z)/factoraa chem(it,kt,jt,p_xyl) = cbox(ixyl_z)/factoraa chem(it,kt,jt,p_hono) = cbox(ihono_z)/factoraa chem(it,kt,jt,p_hno4) = cbox(ihno4_z)/factoraa chem(it,kt,jt,p_ket) = cbox(iaone_z)/factoraa chem(it,kt,jt,p_mgly) = cbox(imgly_z)/factoraa chem(it,kt,jt,p_onit) = cbox(ionit_z)/factoraa chem(it,kt,jt,p_csl) = cbox(icres_z)/factoraa chem(it,kt,jt,p_iso) = cbox(iisop_z)/factoraa chem(it,kt,jt,p_ho) = cbox(ioh_z)/factoraa chem(it,kt,jt,p_ho2) = cbox(iho2_z)/factoraa chem(it,kt,jt,p_hcl) = cbox(ihcl_z)/factoraa chem(it,kt,jt,p_ch3o2) = cbox(ich3o2_z)/factoraa chem(it,kt,jt,p_ethp) = cbox(iethp_z)/factoraa chem(it,kt,jt,p_ch3oh) = cbox(ich3oh_z)/factoraa chem(it,kt,jt,p_c2h5oh) = cbox(ic2h5oh_z)/factoraa chem(it,kt,jt,p_par) = cbox(ipar_z)/factoraa chem(it,kt,jt,p_to2) = cbox(ito2_z)/factoraa chem(it,kt,jt,p_cro) = cbox(icro_z)/factoraa chem(it,kt,jt,p_open) = cbox(iopen_z)/factoraa chem(it,kt,jt,p_op3) = cbox(irooh_z)/factoraa chem(it,kt,jt,p_c2o3) = cbox(ic2o3_z)/factoraa chem(it,kt,jt,p_ro2) = cbox(iro2_z)/factoraa chem(it,kt,jt,p_ano2) = cbox(iano2_z)/factoraa chem(it,kt,jt,p_nap) = cbox(inap_z)/factoraa chem(it,kt,jt,p_xo2) = cbox(ixo2_z)/factoraa chem(it,kt,jt,p_xpar) = cbox(ixpar_z)/factoraa chem(it,kt,jt,p_isoprd) = cbox(iisoprd_z)/factoraa chem(it,kt,jt,p_isopp) = cbox(iisopp_z)/factoraa chem(it,kt,jt,p_isopn) = cbox(iisopn_z)/factoraa chem(it,kt,jt,p_isopo2) = cbox(iisopo2_z)/factoraa chem(it,kt,jt,p_dms) = cbox(idms_z)/factoraa chem(it,kt,jt,p_msa) = cbox(imsa_z)/factoraa chem(it,kt,jt,p_dmso) = cbox(idmso_z)/factoraa chem(it,kt,jt,p_dmso2) = cbox(idmso2_z)/factoraa chem(it,kt,jt,p_ch3so2h) = cbox(ich3so2h_z)/factoraa chem(it,kt,jt,p_ch3sch2oo) = cbox(ich3sch2oo_z)/factoraa chem(it,kt,jt,p_ch3so2) = cbox(ich3so2_z)/factoraa chem(it,kt,jt,p_ch3so3) = cbox(ich3so3_z)/factoraa chem(it,kt,jt,p_ch3so2oo) = cbox(ich3so2oo_z)/factoraa chem(it,kt,jt,p_ch3so2ch2oo) = cbox(ich3so2ch2oo_z)/factoraa chem(it,kt,jt,p_mtf) = cbox(imtf_z)/factoraa return end subroutine mapgas_tofrom_host !*********************************************************************** ! subr set_gaschem_allowed_regimes ! ! purpose: determines which gas-phase chemistry regimes are allowed based ! on which species are active in the simulation ! ! author : ! date : ! ! ---------------------------------------------------------------------- subroutine set_gaschem_allowed_regimes( lunerr, & igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 ) ! ! determines which gas-phase chemistry regimes are allowed based ! on which species are active in the simulation ! use module_configure, only: & p_qv, & p_so2, p_sulf, p_no2, p_no, p_o3, & p_hno3, p_h2o2, p_ald, p_hcho, p_op1, & p_op2, p_paa, p_ora1, p_ora2, p_nh3, & p_n2o5, p_no3, p_pan, p_hc3, p_hc5, & p_hc8, p_eth, p_co, p_ol2, p_olt, & p_oli, p_tol, p_xyl, p_aco3, p_tpan, & p_hono, p_hno4, p_ket, p_gly, p_mgly, & p_dcb, p_onit, p_csl, p_iso, p_ho, & p_ho2, & p_hcl, p_ch3o2, p_ethp, p_ch3oh, p_c2h5oh, & p_par, p_to2, p_cro, p_open, p_op3, & p_c2o3, p_ro2, p_ano2, p_nap, p_xo2, & p_xpar, p_isoprd, p_isopp, p_isopn, p_isopo2, & p_dms, p_msa, p_dmso, p_dmso2, p_ch3so2h, & p_ch3sch2oo, p_ch3so2, p_ch3so3, p_ch3so2oo, p_ch3so2ch2oo, & p_mtf use module_state_description, only: param_first_scalar use module_data_cbmz implicit none ! subr arguments integer lunerr integer igaschem_allowed_m1, igaschem_allowed_m2, & igaschem_allowed_m3, igaschem_allowed_m4 ! local variables integer nactive, ndum, p1st character*80 msg ! index for first "active" scalar (= 2) p1st = param_first_scalar ! determine if regime 1 is allowed ! (note: p_xxx>1 if xxx is active, p_xxx=1 if inactive) if (p_qv .lt. p1st) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) msg = '*** water vapor IS NOT ACTIVE' call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if ! determine if regime 1 is allowed ! (note: p_xxx>1 if xxx is active, p_xxx=1 if inactive) nactive = 0 ndum = 27 if (p_no .ge. p1st) nactive = nactive + 1 if (p_no2 .ge. p1st) nactive = nactive + 1 if (p_no3 .ge. p1st) nactive = nactive + 1 if (p_n2o5 .ge. p1st) nactive = nactive + 1 if (p_hono .ge. p1st) nactive = nactive + 1 if (p_hno3 .ge. p1st) nactive = nactive + 1 if (p_hno4 .ge. p1st) nactive = nactive + 1 if (p_o3 .ge. p1st) nactive = nactive + 1 ! if (p_o1d .ge. p1st) nactive = nactive + 1 ! if (p_o3p .ge. p1st) nactive = nactive + 1 if (p_ho .ge. p1st) nactive = nactive + 1 if (p_ho2 .ge. p1st) nactive = nactive + 1 if (p_h2o2 .ge. p1st) nactive = nactive + 1 if (p_co .ge. p1st) nactive = nactive + 1 if (p_so2 .ge. p1st) nactive = nactive + 1 if (p_sulf .ge. p1st) nactive = nactive + 1 ! if (p_nh3 .ge. p1st) nactive = nactive + 1 ! if (p_hcl .ge. p1st) nactive = nactive + 1 if (p_eth .ge. p1st) nactive = nactive + 1 if (p_ch3o2 .ge. p1st) nactive = nactive + 1 if (p_ethp .ge. p1st) nactive = nactive + 1 if (p_hcho .ge. p1st) nactive = nactive + 1 if (p_ch3oh .ge. p1st) nactive = nactive + 1 if (p_c2h5oh .ge. p1st) nactive = nactive + 1 if (p_op1 .ge. p1st) nactive = nactive + 1 if (p_op2 .ge. p1st) nactive = nactive + 1 if (p_ald .ge. p1st) nactive = nactive + 1 if (p_ora1 .ge. p1st) nactive = nactive + 1 if (p_pan .ge. p1st) nactive = nactive + 1 if (p_ora2 .ge. p1st) nactive = nactive + 1 if (p_c2o3 .ge. p1st) nactive = nactive + 1 if (nactive .le. 0) then igaschem_allowed_m1 = 0 else if (nactive .eq. ndum) then igaschem_allowed_m1 = 1 else msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90200) 1, nactive, ndum call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if 90200 format( ' error for regime ', i1, ', nactive, nexpected = ', 2i5 ) ! determine if regime 2 is allowed nactive = 0 ndum = 19 if (p_par .ge. p1st) nactive = nactive + 1 if (p_ket .ge. p1st) nactive = nactive + 1 if (p_mgly .ge. p1st) nactive = nactive + 1 if (p_ol2 .ge. p1st) nactive = nactive + 1 if (p_olt .ge. p1st) nactive = nactive + 1 if (p_oli .ge. p1st) nactive = nactive + 1 if (p_tol .ge. p1st) nactive = nactive + 1 if (p_xyl .ge. p1st) nactive = nactive + 1 if (p_csl .ge. p1st) nactive = nactive + 1 if (p_to2 .ge. p1st) nactive = nactive + 1 if (p_cro .ge. p1st) nactive = nactive + 1 if (p_open .ge. p1st) nactive = nactive + 1 if (p_onit .ge. p1st) nactive = nactive + 1 if (p_op3 .ge. p1st) nactive = nactive + 1 if (p_ro2 .ge. p1st) nactive = nactive + 1 if (p_ano2 .ge. p1st) nactive = nactive + 1 if (p_nap .ge. p1st) nactive = nactive + 1 if (p_xo2 .ge. p1st) nactive = nactive + 1 if (p_xpar .ge. p1st) nactive = nactive + 1 if (nactive .le. 0) then igaschem_allowed_m2 = 0 else if (nactive .eq. ndum) then igaschem_allowed_m2 = 2 else msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90200) 2, nactive, ndum call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if ! determine if regime 3 is allowed nactive = 0 ndum = 5 if (p_iso .ge. p1st) nactive = nactive + 1 if (p_isoprd .ge. p1st) nactive = nactive + 1 if (p_isopp .ge. p1st) nactive = nactive + 1 if (p_isopn .ge. p1st) nactive = nactive + 1 if (p_isopo2 .ge. p1st) nactive = nactive + 1 if (nactive .le. 0) then igaschem_allowed_m3 = 0 else if (nactive .eq. ndum) then igaschem_allowed_m3 = 3 else msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90200) 3, nactive, ndum call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if ! determine if regime 4 is allowed nactive = 0 ndum = 11 if (p_dms .ge. p1st) nactive = nactive + 1 if (p_msa .ge. p1st) nactive = nactive + 1 if (p_dmso .ge. p1st) nactive = nactive + 1 if (p_dmso2 .ge. p1st) nactive = nactive + 1 if (p_ch3so2h .ge. p1st) nactive = nactive + 1 if (p_ch3sch2oo .ge. p1st) nactive = nactive + 1 if (p_ch3so2 .ge. p1st) nactive = nactive + 1 if (p_ch3so3 .ge. p1st) nactive = nactive + 1 if (p_ch3so2oo .ge. p1st) nactive = nactive + 1 if (p_ch3so2ch2oo .ge. p1st) nactive = nactive + 1 if (p_mtf .ge. p1st) nactive = nactive + 1 if (nactive .le. 0) then igaschem_allowed_m4 = 0 else if (nactive .eq. ndum) then igaschem_allowed_m4 = 4 else msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90200) 4, nactive, ndum call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if ! regime 1 must always be allowed if (igaschem_allowed_m1 .le. 0) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90300) 1 call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if 90300 format( ' regime ', i1, ' must always be allowed' ) ! if regime 2 is allowed, then regime 1 must be allowed if (igaschem_allowed_m2 .gt. 0) then if (igaschem_allowed_m1 .le. 0) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90400) 2, 1 call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if end if 90400 format( ' regime ', i1, ' allowed BUT regime ', i1, ' unallowed' ) ! if regime 3 is allowed, then regimes 1&2 must be allowed if (igaschem_allowed_m3 .gt. 0) then if (igaschem_allowed_m1 .le. 0) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90400) 3, 1 call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) else if (igaschem_allowed_m2 .le. 0) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90400) 3, 2 call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if end if ! if regime 4 is allowed, then regime 1 must be allowed if (igaschem_allowed_m4 .gt. 0) then if (igaschem_allowed_m1 .le. 0) then msg = '*** subr set_gaschem_allowed_regimes' call peg_message( lunerr, msg ) write(msg,90400) 4, 1 call peg_message( lunerr, msg ) call peg_error_fatal( lunerr, msg ) end if end if return end subroutine set_gaschem_allowed_regimes !*********************************************************************** ! subr gasphotoconstants ! ! purpose: copy photolytic rate constants from host arrays to local array ! !----------------------------------------------------------------------- subroutine gasphotoconstants( rk_photo, & i_boxtest_units_convert, & it,jt,kt, ims,ime, jms,jme, kms,kme, & ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, & ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, & ph_ch3o2h, ph_n2o5 ) ! ! copies photolytic rate constants from host arrays to local arrays ! note1: currently 8 rate constants are scaled to other rate constants ! as is done in zz06gasphotolysis.f ! note2: currently the n2o5 rate is set to zero ! use module_data_cbmz implicit none ! subr arguments integer it,jt,kt, ims,ime, jms,jme, kms,kme integer i_boxtest_units_convert real rk_photo(nphoto) real, dimension( ims:ime, kms:kme, jms:jme ) :: & ph_o31d, ph_o33p, ph_no2, ph_no3o2, ph_no3o, ph_hno2, & ph_hno3, ph_hno4, ph_h2o2, ph_ch2or, ph_ch2om, & ph_ch3o2h, ph_n2o5 ! local variables real ft rk_photo(:) = 0.0 ! these from wrf/madronnich rate constants rk_photo(jphoto_no2) = ph_no2(it,kt,jt) rk_photo(jphoto_no3) = ph_no3o(it,kt,jt) & + ph_no3o2(it,kt,jt) rk_photo(jphoto_o3a) = ph_o33p(it,kt,jt) rk_photo(jphoto_o3b) = ph_o31d(it,kt,jt) rk_photo(jphoto_hono) = ph_hno2(it,kt,jt) rk_photo(jphoto_hno3) = ph_hno3(it,kt,jt) rk_photo(jphoto_hno4) = ph_hno4(it,kt,jt) rk_photo(jphoto_h2o2) = ph_h2o2(it,kt,jt) rk_photo(jphoto_ch3ooh) = ph_ch3o2h(it,kt,jt) rk_photo(jphoto_hchoa) = ph_ch2or(it,kt,jt) rk_photo(jphoto_hchob) = ph_ch2om(it,kt,jt) rk_photo(jphoto_n2o5) = ph_n2o5(it,kt,jt) ! these scaled to other rate constants rk_photo(jphoto_ethooh) = 0.7 *rk_photo(jphoto_h2o2) rk_photo(jphoto_ald2) = 4.6e-4*rk_photo(jphoto_no2) rk_photo(jphoto_aone) = 7.8e-5*rk_photo(jphoto_no2) rk_photo(jphoto_mgly) = 9.64 *rk_photo(jphoto_hchoa) rk_photo(jphoto_open) = 9.04 *rk_photo(jphoto_hchoa) rk_photo(jphoto_rooh) = 0.7 *rk_photo(jphoto_h2o2) rk_photo(jphoto_onit) = 1.0e-4*rk_photo(jphoto_no2) rk_photo(jphoto_isoprd) = .025 *rk_photo(jphoto_hchob) ! convert from (1/min) to (1/s) ! (except when i_boxtest_units_convert = 10 or 20) ft = 60.0 if (i_boxtest_units_convert .eq. 10) ft = 1.0 if (i_boxtest_units_convert .eq. 20) ft = 1.0 if (ft .ne. 1.0) then rk_photo(:) = rk_photo(:)/ft end if return end subroutine gasphotoconstants !----------------------------------------------------------------------- end module module_cbmz