****************************************************************** * MEANING OF THE VARIABLES * * dictionary and species lists : * ------------------------------ * - ninorg : number of inorganic species * - nrec : number of species recorded * - nrectot : total number of record (not including inorganic) * - dict(j) : dictionary line (name + formula + functional * group info) of species number j * - namlst(j) : name (lco=6 characters) of the species already * used at position number j * - inorglst(j) : list of inorganic species - include the name + * formula + functional group info * * - name : used for a given element of namlst * - chem : used for a given formula * * data in stack * -------------- * - nhldvoc : number of (stable) VOC in the stack * - holdvoc(i) : list of the VOC in the stack (name[a6]+ * formula[a120]+stabl[i3]+level[i3]) * - nhldrad : number of radical in the stack * - holdrad(i) : list of the radicals in the stack(name[a6]+ * formula[a120]+stabl[i3]+level[i3]) * - stabl : gives the number of generations that were necessary * to produce the current species * - level : gives the number of intermediates (including * radical) that was necessary to produce the * current species * * species having a name fixed (not set by the generator) * --------------------------- * - nfn : total nb. of species having a fixed name * - namfn(i) : table of the fixed name (6 character) * - chemfn(i) : formula corresponding the ith species having a * fixed name * * species having a structure that can not be held by the generator * (cyclic molecule for example). * --------------------------- * - nspsp : total number of "special" species (e.g. furane) * - dictsp(i) : list of the "known" special species * * data for photolytic reaction (see photo.dat) * ---------------------------- * - jdat : nb. of photolytic reaction given as input in * the file photo.dat * - jlabel(i) : ID number of photolytic reaction i * - jchem(i) : chemical formula of the reactant for reaction i * - jprod(2,i) : the two main species arising from the photolytic * bond break. If a third product is formed, then * it must be stored in coprodj(i) * - coprodj(i) : additional product in reaction (i) - species * in coprod are not the organic fragments but * only minor inorganic coproduct * - jlab40(i) : label of the ith data for which J value has been * evaluated (at a zenith angle of 40) * - j40(i) : J value (at a zenith angle of 40) corresponding * to the reaction having the label jlab40(i) * * Data for the reaction with OH * ----------------------------- * ARRHENIUS PARAMETER (see oh_rate.dat) * - nohdat : total number of known OH rate constant * - ohchem(i) : chemical formula of ith species having a known OH * rate constant * - ohrate(i,3) : arrhenius parameter of the known OH rate constant * corresponding to the ith species * KNOWN REACTION (see oh_prod.dat) * - nkwoh : number of VOC having known products for the react * with OH * - kwohrct(i) : chemical formula of the VOC having a known * chemistry for the reaction with OH * - kwohpd(i,j) : chemical formula of the jth product for the * reaction of the ith VOC with OH * - kwohcopd(i,j,k): name (6 character) of the kth coproduct of the * jth product for the reaction of the ith VOC with OH * - kwohyld(i,j) : yield of the jth product for the reaction of the * ith VOC with OH * * Data for the reactions with NO3 (see no3_rate.dat) * ----------------------------- * ARRHENIUS PARAMETER (see no3_rate.dat) * - nno3dat : total number of known NO3 rate constant * - no3chem(i) : chemical formula of ith species having a known NO3 * rate constant * - no3rate(i,3) : arrhenius parameter of the known NO3 rate constant * corresponding to the ith species * KNOWN REACTION (see no3_prod.dat) * - nkwno3 : number of VOC having known products for the react * with OH * - kwno3rct(i) : chemical formula of the VOC having a known * chemistry for the reaction with OH * - kwno3pd(i,j) : chemical formula of the jth product for the * reaction of the ith VOC with OH * - kwno3copd(i,j,k): name (6 character) of the k thcoproduct of the * jth product for the reaction of the ith VOC with OH * - kwno3yld(i,j): yield of the jth product for the reaction of the * ith VOC with OH * * Data for the reactions with O3 (see o3_rate.dat) * ----------------------------- * - no3dat : total number of known O3 rate constant * - o3chem(i) : chemical formula of ith species having a known O3 * rate constant * - o3rate(i,3) : arrhenius parameter of the known O3 rate constant * corresponding to the ith species * * Data for the reactions that cannot be managed by the generator * ------------------------------------------------------------- * - noge : number of reaction given (oge=out of generator) * - ogertve(i,3): reactants for reaction i. * - ogeprod(i,j): products for the reaction i * - ogearh(i,3) : arrehnius coefficients (A, n, Ea) * - ogestoe(i,j): stochiometric coefficients for product j * - ogelab(i) : label for the reaction (if EXTRA or HV) * - ogeaux(i,7) : auxiliary information for reaction i (e.g. data for * fall off reaction * * Benson group data (see benson.dat) * ---------------------------- * - nbson : total number of benson group * - bsongrp(i) : table of benson group * - bsonval(i) : heat of formation corresponding to group i * * primary species given as input * ---------------------------- * - ninp : number of "primary" species * - input(i) : table of the primary species * * general parameter controlling how the schemes are generated * - cut_off : ratio under which a reactive channel is not taken * into account (e.g. if a reactive path represent less * than 1 % of the total paths) * * Data for Pvap estimates * ---------------------------- * - molmass(i) : molar mass of species in dict(i) * - Tb : boiling point * - HBN : Hydrogen bond number * - tau : number of torsional bond ******************************************************************* PROGRAM main IMPLICIT NONE INCLUDE 'general.h' INCLUDE 'organic.h' INCLUDE 'common.h' INCLUDE 'isomer.h' REAL etime, durat REAL(kind = 4) tdat(2) CHARACTER*(ldi) dict(mni) CHARACTER*(lcf) rdct,rdcttmp CHARACTER*(lfo) chem CHARACTER*(lco) name, namlst(mni) CHARACTER*(lco) tr_species(mni) INTEGER srch,cnum INTEGER nc, ind, ic REAL brch, dbrch(mni) INTEGER i, j,k,l,ntr_sp REAL cut_off CHARACTER*(lfl) fgrp * special species * INTEGER nspsp * CHARACTER*(ldi) dictsp(mfn) * LOGICAL lospsp(mfn) * stack data INTEGER level INTEGER stabl CHARACTER*(lst) holdvoc(mlv) INTEGER nhldvoc CHARACTER*(lst) holdrad(mra) INTEGER nhldrad * Benson group data INTEGER nbson REAL bsonval(mbg) CHARACTER*(lgb) bsongrp(mbg) * heats of formation REAL heat,dHeatf * Data for the conversion to hydrofuran INTEGER fdhf,fd INTEGER ncha,itest CHARACTER*(lco) chatab(10000) * Data for the reaction with OH INTEGER nohdat CHARACTER*(lfo) ohchem(mrd) REAL ohrate(mrd,3) INTEGER nkwoh,nkwohpd(mkr) CHARACTER*(lfo) kwohrct(mkr), kwohpd(mkr,mnr) CHARACTER*(lco) kwohcopd(mkr,mnr,mcp) REAL kwohyld(mkr,mnr) * Data for the reaction with NO3 INTEGER nno3dat CHARACTER*(lfo) no3chem(mrd) REAL no3rate(mrd,3) INTEGER nkwno3,nkwno3pd(mkr) CHARACTER*(lfo) kwno3rct(mkr), kwno3pd(mkr,mnr) CHARACTER*(lco) kwno3copd(mkr,mnr,mcp) REAL kwno3yld(mkr,mnr) * Data for the reaction with O3 CHARACTER*(lfo) o3chem(mrd) REAL o3rate(mrd,3) INTEGER no3dat INTEGER nkwo3,nkwo3pd(mkr) CHARACTER*(lfo) kwo3rct(mkr), kwo3pd(mkr,mnr) CHARACTER*(lco) kwo3copd(mkr,mnr,mcp) REAL kwo3yld(mkr,mnr) * Data for the reaction with RO CHARACTER*(lfo) rokrct(mkr,2),rokprd(mkr,3) INTEGER nrokr REAL roarrh(mkr,3),rocost(mkr,3) * Data for the reaction with RO2 CHARACTER*(lfo) ro2krct(mkr,2),ro2kprd(mkr,3) INTEGER nro2kr REAL ro2arrh(mkr,3),ro2cost(mkr,3) * Data for the reaction with RCO2--Added April 2015 Renee CHARACTER*(lfo) rcoo2krct(mkr,2),rcoo2kprd(mkr,3) INTEGER nrcoo2kr REAL rcoo2arrh(mkr,3),rcoo2cost(mkr,3) * Data for the reactions that cannot be managed by the generator CHARACTER*(lfo) ogertve(mog,3) CHARACTER*(lfo) ogeprod(mog,mnp) REAL ogearh(mog,3) REAL ogestoe(mog,mnp) REAL ogeaux(mog,7) INTEGER noge,ogelab(mog) * primary species INTEGER ninp CHARACTER*(lfo) input(mps) CHARACTER*(lgr) group(mca) INTEGER bond(mca,mca) LOGICAL chromophores INTEGER nfn CHARACTER*(lco) namfn(mfn) CHARACTER*(lfo) chemfn(mfn) REAL tbm CHARACTER*(110) line INTEGER dbflg, nring INTEGER nca CHARACTER*(lfo) jchem(mkr),jprod(2,mkr) CHARACTER*(lco) coprodj(mkr) INTEGER jlabel(mkr),jdat,jlab40(mkr) REAL j40(mkr) INTEGER nam_int(12) CHARACTER*(llin) filename INTEGER id,srh5 INTEGER ndepth LOGICAL lorad REAL tempcut_off CHARACTER*(lfo) tchem INTEGER rjg(mri,2) ! ring-join group pairs * data for autooxidation INTEGER oxdone * data for dimer INTEGER nalde,nalco,nooh,ncooh,ncoooh REAL malde,malco,mooh,mcooh,mcoooh LOGICAL lodimer * data for Pvap estimates REAL molmass(mni) REAL Tb, HBN, tau REAL log10Psat, mw REAL vdmol REAL Nangroup(219),logPvap,T,dB,Trb REAL bk(0:30,4),simpgroup(0:30) REAL dH_nan,dH_sim,dH_MY,psat ! data for known Henry constants INTEGER nhenrydat CHARACTER*(lfo) chemhenry(mrd) REAL henrydat(mrd,3),rf,cf ! data for known hydration constants INTEGER nhydratdat CHARACTER*(lfo) chemhydrat(mrd) REAL hydratdat(mrd) REAL keff * parameters for parallel running !$ INTEGER omp_get_num_threads !$ INTEGER omp_get_thread_num INTEGER ired REAL dummy CHARACTER(lsb) :: progname='*main* ' CHARACTER(ler) :: mesg * flags for molecular parameter writing ** now activated using "postflag" ** * (to save time, use only for last precursor iteration) ! INTEGER,PARAMETER :: wrt_dhf = 1 ! (note that 'A' and 'W' masses are only written if wrt_pvap = 1) ! INTEGER,PARAMETER :: wrt_pvap = 0 ! INTEGER,PARAMETER :: wrt_henry = 0 * time estimate INCLUDE 'timecommon.h' CHARACTER*10 date, begtime, endtime REAL elapsed, time_diff INTEGER irewind INTEGER ih, in, io, ir, is, ifl, ibr, icl INTEGER tc LOGICAL flg !------------------------------ ! DEFAULT USER-SET PARAMETERS & FLAGS !------------------------------ ! vapor pressure exponent below which species partition completely to aerosol. REAL :: critvp = -13. ! branching ratio below which reaction pathways not considered REAL :: cutoff_default = 0.05 REAL :: cutoff_OH = 0.00 REAL :: cutoff_O3 = 0.00 REAL :: cutoff_NO3 = 0.00 REAL :: cutoff_PAN = 0.00 REAL :: cutoff_HV = 0.00 REAL :: cutoff_RO = 0.00 REAL :: cutoff_RO2 = 0.00 REAL :: cutoff_RCOO2 = 0.00 REAL :: temp = 298. ! maximum number of generations allowed INTEGER :: maxgen = 20 ! number of species before rewind INTEGER :: nrewind = 100 !------------------------------ !Declare namelist to be used for user parameters file NAMELIST /userparams/ critvp, cutoff_default, cutoff_OH, & cutoff_O3, cutoff_NO3, cutoff_PAN, cutoff_HV, cutoff_RO, & cutoff_RO2, cutoff_RCOO2, temp, maxgen, nrewind, & z_conformer_prop !Declare namelist to be used for user flags NAMELIST /flags/ wtflag,wtopeflag,per3_typefg,high_NOxfg, & zero_NOXfg,iflag,pvap_fg,kdissfg,kisomfg,hoadd_c1_fg,flgdhf, & autoox_fg,dimer_fg,masstransfg, spinup_fg, spinup_mech, & aeropart_fg, wallpart_fg, criegee_fg, stab_criegee_fg, & ro2ho2_fg, isopsoa_fg, aerophot_fg, OFR_fg, VBS_fg, kill_fg, & ro2no2_fg, ro2cond_fg, ro2dep_fg NAMELIST /manage/ prevflag,postflag CALL date_and_time(date,begtime) telstd=0. * ---------- * initialize * ---------- irewind=0 ninorg = 0 inorglst(:) = ' ' line = ' ' ndepth=0 * data in the stack lorad=.false. nhldvoc = 0 nhldrad = 0 level = 0 stabl = 0 holdvoc(:) = ' ' holdrad(:) = ' ' group(:) = ' ' ninp=0 ! read user parameter in namelist userparams OPEN(73, file = 'userparams.input') !default values cutoff_PAN = cutoff_default cutoff_HV = cutoff_default cutoff_RO = cutoff_default cutoff_RO2 = cutoff_default cutoff_RCOO2 = cutoff_default ! temporary value for Z-E conformation of -CH=CH- alkenes z_conformer_prop = 0.5 READ(73, nml = userparams ) CLOSE(73) * other parameters WRITE(6,*) ' temp = ', temp tbm = 2.5e19 WRITE(6,*) ' air density = ', tbm brch = 1. ntr_sp = 0 dbrch = 0. diccri=0 tr_species(:) = ' ' ntr_sp=1 tr_species(1) = 'CH4' id = 0 chatab(:)=' ' ncha=0 nalde=0 nalco=0 nooh=0 ncooh=0 ncoooh=0 malde=0 malco=0 mooh=0 ! default flags spinup_fg = 0 ro2ho2_fg = 1 ro2no2_fg = 0 ro2cond_fg = 0 ro2dep_fg = 0 OFR_fg = 0 VBS_fg = 0 ! note for future development: kill_fg could potentially have different ! values for different severities of error kill_fg = 0 !Read user overwrites for the flags OPEN(73, file = 'userparams.input') READ(73, nml = flags ) CLOSE(73) OPEN(73, file = 'manage.input') READ(73, nml = manage ) CLOSE(73) CALL checkflags() * -------------------------------------------------- * WRITE MECHANISM PARAMETERS TO STDOUT, FOR LOGGING * -------------------------------------------------- OPEN(16,file='scheme.log') WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' GECKO VERSION = GECKO-A_September_2016 ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' maxgen = ', maxgen WRITE(16,*) ' nrewind = ', nrewind WRITE(16,*) ' temp = ', temp WRITE(16,*) ' air density = ', tbm WRITE(16,*) ' ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' --- Chem Pathways --- ' WRITE(16,*) ' ---------------------------- ' IF(flgdhf.EQ.0)THEN WRITE(16,*) ' no DHF formation ' ELSEIF (flgdhf.EQ.1)THEN WRITE(16,*) ' includes DHF formation ' ENDIF IF(autoox_fg.EQ.0)THEN WRITE(16,*) ' no auto-oxidation of peroxy radicals ' ELSEIF (autoox_fg.EQ.1)THEN WRITE(16,*) ' includes auto-oxidation of peroxy radicals' ENDIF IF(dimer_fg.EQ.0)THEN WRITE(16,*) ' no dimerisation ' ELSEIF (dimer_fg.EQ.1)THEN WRITE(16,*) ' includes dimerisation' ENDIF WRITE(16,*) ' ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' --------- SAR -------- ' WRITE(16,*) ' ---------------------------- ' IF(pvap_fg.EQ.1)THEN WRITE(16,*) ' vapor pressure scheme = M&Y ' ELSEIF (pvap_fg.EQ.2)THEN WRITE(16,*) ' vapor pressure scheme = Nannoolal ' ELSEIF (pvap_fg.EQ.3)THEN WRITE(16,*) ' vapor pressure scheme = SIMPOL-1 ' ENDIF WRITE(6,*) ' critical vapor pressure = ', critvp IF(kdissfg.EQ.1)THEN WRITE(16,*) ' alkoxy decomposition method = Atkinson 2007 ' ELSEIF (kdissfg.EQ.2)THEN WRITE(16,*) ' alkoxy decomposition method = Vereecken 2009 ' ENDIF IF(kisomfg.EQ.1)THEN WRITE(16,*) ' kisomfg = 1 ' ELSEIF (kisomfg.EQ.2)THEN WRITE(16,*) ' kisomfg = 2 ' ENDIF IF(hoadd_c1_fg.EQ.1)THEN WRITE(16,*) ' OH addition SAR = Peeters 1997 ' ELSEIF (hoadd_c1_fg.EQ.2)THEN WRITE(16,*) ' OH addition SAR = Ziemann 2009 ' ELSEIF (hoadd_c1_fg.EQ.3)THEN WRITE(16,*) ' OH addition SAR = Jenkin 2016 (MAGNIFY) ' ENDIF IF(masstransfg.EQ.1)THEN WRITE(16,*) ' mass transfer method = thermodynamic ' ELSEIF (masstransfg.EQ.2)THEN WRITE(16,*) ' mass transfer method = dynamic ' ENDIF IF (ro2ho2_fg .eq. 2) then WRITE(16,*) ' RO2+HO2 = Wennberg et al., 2018 ' ELSEIF (ro2ho2_fg .eq. 1) then WRITE(16,*) ' RO2+HO2 = Orlando 2013, Hasson 2013 ' WRITE(16,*) ' H-abstr group contribution for -OOH = 3.5' ENDIF IF (ro2no2_fg .eq. 1) then WRITE(16,*) ' RO2+NO2 = Orlando pers comm 2020 ' ELSEIF (ro2no2_fg .eq. 0) then WRITE(16,*) ' RO2+NO2 = disallowed ' ENDIF IF (ro2dep_fg .eq. 1) then WRITE(16,*) ' RO2 deposition enabled ' ENDIF IF (ro2cond_fg .eq. 1) then WRITE(16,*) ' RO2 condensation per R(OOH)enabled ' ENDIF WRITE(16,*) ' ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' ----- Reductions ----- ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' critical vapor pressure = ', critvp IF(iflag.EQ.1)THEN WRITE(16,*) ' isomerisation allowed: miso =',miso ELSE WRITE(16,*) ' no isomerisation ' ENDIF WRITE(16,*) ' high-NOx flag = ',high_NOxfg WRITE(16,*) ' zero-NOx flag = ',zero_NOxfg WRITE(16,*) ' per3_typeflag = ',per3_typefg WRITE(16,*) ' ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' ----- Flags ----- ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' cutoff_default = ', cutoff_default WRITE(16,*) ' cutoff_OH = ', cutoff_OH WRITE(16,*) ' cutoff_O3 = ', cutoff_O3 WRITE(16,*) ' cutoff_NO3 = ', cutoff_NO3 WRITE(16,*) ' cutoff_PAN = ', cutoff_PAN WRITE(16,*) ' cutoff_HV = ', cutoff_HV WRITE(16,*) ' cutoff_RO = ', cutoff_RO WRITE(16,*) ' cutoff_RO2 = ', cutoff_RO2 WRITE(16,*) ' cutoff_RCOO2 = ', cutoff_RCOO2 WRITE(16,*) ' ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' ----- Precursors ----- ' WRITE(16,*) ' ---------------------------- ' WRITE(16,*) ' ' * ------------------------------------------ * OPEN SOME OUTPUT FILES * ------------------------------------------ IF (wtflag.NE.0) OPEN(35,file='debug.out') IF (wtflag.NE.0) OPEN(26,file='info_ro.out') OPEN(14,file='lstprim.out') OPEN(99,file='warnings.out') * ------------------------------------------ * READ ALL DATA * ------------------------------------------ OPEN(10,file='../DATA/simpol.dat',form='FORMATTED', status='OLD') DO i=0,30 READ(10,'(f10.4,2x,f9.6,2x,f11.8,2x,f11.8)') & bk(i,1),bk(i,2),bk(i,3),bk(i,4) ENDDO CLOSE(10) * read species already in the dictionary (inorganic and C1 organic * species) and sort namlst. WRITE (6,*) ' ...reading species already in the dictionary' * IF THIS IS A STAND_ALONE RUN, OR THE FIRST IN A SEQUENCE ..... IF(prevflag.EQ.0)THEN CALL rddict(ninorg,nrec,nrectot, & dict,namlst,inorglst,filename) if (spinup_mech .eq. 1) then filename = '../DATA/spinup_dic.dat' CALL rddict(ninorg,nrec,nrectot, & dict,namlst,inorglst, & filename) endif * ... read inorganic reactions and CH4 chemistry and write out WRITE (6,*) ' ...reading inorganic reactions' IF(ofr_fg.EQ.1)THEN filename='../DATA/inorg_mch_pam.dat' WRITE(6,*)"USING OFR MECH" ELSE filename='../DATA/inorg_mch.dat' WRITE(6,*)"USING STANDARD MECH" ENDIF CALL rdfixmch(filename) WRITE (6,*) ' ...reading CH4 chemistry' filename='../DATA/singlec_mch.dat' CALL rdfixmch(filename) IF(spinup_mech.EQ.1) THEN WRITE(6,*) ' ...reading spinup mechanism' filename='../DATA/spinup_mch.dat' CALL rdfixmch(filename) ENDIF IF (isopsoa_fg.EQ.1) THEN WRITE(6,*) ' ...reading isoprene soa mechanism' filename='../DATA/isop_soa_mch.dat' call rdfixmch(filename) ENDIF IF (VBS_fg == 1) then filename='../DATA/vbs_mch.dat' call rdfixmch(filename) ENDIF * ... write CH3O2+counters reactions WRITE (6,*) ' ...writing CH3O2+counters reactions' CALL rxmeo2ro2 * IF WE NEED TO READ IN A PREVIOUSLY-CREATED DICTIONARY... * ... (Don't need to write single-C and inorg reactions) ELSE WRITE (6,*) ' ...reading previous dictionary' CALL prevdict(ninorg,nrec,nrectot, & dict,namlst,inorglst) ENDIF * END SECTION ABOUT PREVIOUSLY_CREATED DICTIONARY * sort list of 6-character names * sortbis is VERY SLOW - use unix system call first if possible CALL sortunix(nrec,namlst) CALL sortbis(nrec,namlst) write(6,*) 'nrectot=',nrectot * read known species having a "special" chemical structure that * may be produced by the generator WRITE (6,*) ' ...reading special species' CALL rdspsp(nspsp,lospsp,dictsp) * read fixed names for the most common species WRITE (6,*) ' ...reading fixed name species' CALL rdfixnam(nfn,namfn,chemfn) * read bsongrp group names and their heat of formation : WRITE (6,*) ' ...reading benson groups' CALL rdbenson(nbson,bsongrp,bsonval) * read known species for photolysis reactions WRITE (6,*) ' ...reading photolysis reactions' CALL rdhv(jdat,jlabel,jchem,jprod,coprodj,jlab40,j40) * read prescribed rate constant (OH, NO3, O3) WRITE (6,*) ' ...reading OH rate constants' filename='../DATA/oh_rate.dat' CALL rdrate(filename,nohdat,ohchem,ohrate) WRITE (6,*) ' ...reading NO3 rate constants' filename='../DATA/no3_rate.dat' CALL rdrate(filename,nno3dat,no3chem,no3rate) WRITE (6,*) ' ...reading O3 rate constants' filename='../DATA/o3_rate.dat' CALL rdrate(filename,no3dat,o3chem,o3rate) * read prescribed reactions (OH, NO3, O3) WRITE (6,*) ' ...reading known VOC+OH reactions ...' filename='../DATA/oh_prod.dat' CALL rdkwreact(filename, & nkwoh,nkwohpd,kwohrct,kwohpd,kwohcopd,kwohyld) WRITE (6,*) ' ...reading known VOC+NO3 reactions ...' filename='../DATA/no3_prod.dat' CALL rdkwreact(filename, & nkwno3,nkwno3pd,kwno3rct,kwno3pd,kwno3copd,kwno3yld) WRITE (6,*) ' ...reading known VOC+O3 reactions ...' filename='../DATA/o3_prod.dat' CALL rdkwreact(filename, & nkwo3,nkwo3pd,kwo3rct,kwo3pd,kwo3copd,kwo3yld) WRITE (6,*) ' ...reading known RO reactions ...' filename='../DATA/ro.dat' CALL rdkwrad(filename,nrokr,rokrct,rokprd,roarrh,rocost) WRITE (6,*) ' ...reading known RO2 reactions ...' filename='../DATA/ro2.dat' CALL rdkwrad(filename,nro2kr,ro2krct,ro2kprd,ro2arrh,ro2cost) WRITE (6,*) ' ...reading known RCOO2 reactions ...' filename='../DATA/rcoo2.dat' CALL rdkwrad(filename, & nrcoo2kr,rcoo2krct,rcoo2kprd,rcoo2arrh,rcoo2cost) * read reactions for the special species WRITE (6,*) ' ...reading reactions for "special species" ...' !filename='../DATA/mcm.3.3.1_mxyl.dat' ! self-generating mxyl !filename='../DATA/mcm.3.3.1_myrcene_ocimene.dat' ! mcm mxyl ! mcm mxyl + NCAR/Lindsay Hatch furans ! filename='../DATA/mcm.3.3.1_myrc_oc_furan.dat' ! as above, with Flocke/Orlando additions ! filename='../DATA/mcm.3.3.1_myrc_oc_furan.dat' ! as above, with Flocke/Orlando additions and JMLT mxyl edits ! filename='../DATA/mcm.3.3.1_plus_200514.dat' ! as above, with Flocke/Orlando additions and JMLT mxyl edits ! filename='../DATA/mcm.3.3.1_plus_200518.dat' ! as above, with 2015 MXYL no-NOx improvements !filename='../DATA/mcm.3.3.1_plus_200731.dat' ! as v 200518, with MXYL split off into separate files filename='../DATA/mcm.3.3.1_200814.dat' CALL rdoutgene(filename,1,ninorg,inorglst,nrec,dict,nspsp,dictsp, & noge,ogertve,ogeprod,ogearh,ogestoe,ogelab,ogeaux) * read reactions for m-xylene IF (zero_NOxfg.LT.1) THEN filename='../DATA/mcm.3.3.1_mxyl_orig.dat' ELSE filename='../DATA/mcm.3.3.1_mxyl_noNOx.dat' ENDIF PRINT*,filename CALL rdoutgene(filename,2,ninorg,inorglst,nrec,dict,nspsp,dictsp, & noge,ogertve,ogeprod,ogearh,ogestoe,ogelab,ogeaux) ! read known hydration constants WRITE (6,*) ' ...reading known hydration constants ...' filename='../DATA/hydrat_const.dat' CALL rdhydrat(filename,nhydratdat,chemhydrat,hydratdat) ! read known Henry's law constants WRITE (6,*) ' ...reading known Henrys law constants ...' filename='../DATA/henry_const.dat' CALL rdhenry(filename,nhenrydat,chemhenry,henrydat, & nhydratdat,hydratdat,chemhydrat) * ----------------------------------------------- * LOAD SPECIES FROM CHEMINPUT.DAT (INPUT SPECIES) * ----------------------------------------------- WRITE (6,*) ' ...reading the list of the primary species' filename='./cheminput.dat' CALL rdchemin(filename,ninorg,inorglst,dict,nrec, & nspsp,dictsp, & ninp,input) WRITE (6,*) ' ...done reading the list of the primary species' ! if it is only a spinup run, we won't generate any organic chemistry IF( spinup_fg .eq. 1) THEN ninp = 0 ! need to put a dummy species in the primary species list ! it will be initialised to 0 in the simulation ! lostcarbon is a safe choice: it should never be initialised ! or constrained in a boxmodel simulation WRITE(14,*) "xxlost" ENDIF ! open the dictionaries at the beginning because we need to write info ! about aerosol and wall species WRITE(21,'(a7)') 'SPECIES' ! create more detailed dictionnary. ! for easier extraction ! OPEN(!!22, file="outdat.gasdict") ! WRITE(22,*) "name, chem, nC, nH, nO, nN, molW, functions" * ----------------------------------------------- * LOAD INPUT SPECIES IN STACK - START LOOP * ----------------------------------------------- DO 421 k=1,ninp chem = input(k) WRITE(14,*) chem WRITE(16,*) chem * check the species, give it a name (if not already in the * dictionary), update the dictionary and put the species * at the beginning of the stack (i.e. in holdvoc(1)) CALL in1chm(chem,nrec,nrectot, & nfn,namfn,chemfn,nspsp,dictsp, & dbrch,dict,namlst, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl) ! For 'successive-mode' run, skip chemistry for species already in dictionary. ! (if re-starting after mid-species crash, disable this line.) IF(prevflag.GT.0.AND.level.LT.0) CYCLE * sort namlst: sortbis is VERY SLOW - use unix system call first if possible CALL sortunix(nrec,namlst) CALL sortbis(nrec,namlst) * --------------------------------------------------- * RETRY POINT TO GET THE NEXT SPECIES IN THE STACKS * --------------------------------------------------- * ired is used to sort the stack when species that belong to the * generation are all treated - required to play with isomers ired=1 1000 CONTINUE * Give the priority to the stack holding radical species (holdrad) * If this stack is empty, then take the next stable VOC in * the stack (holdvoc). If both stacks are empty then load next input * species (label 421) IF (holdrad(1).NE.' ') THEN READ (holdrad(1),'(a126,i3,i3)') rdct,stabl,level lorad=.true. ELSE lorad=.false. IF (holdvoc(1).EQ.' ') THEN WRITE(6,*) '..........................................' * If doing single-species successive runs, exit primary species loop GO TO 987 !GO TO 421 ENDIF READ(holdvoc(1),'(a126,i3,i3)') rdct,stabl,level * sort the stack if next generation is reached IF (stabl.eq.ired) THEN ired=ired+1 DO i=1,nhldvoc READ(holdvoc(i),'(a126,i3,i3)') rdct,stabl,level ind = srch(nrec,rdct(7:126),dict) WRITE(60,'(f20.10,2x,a132)') dbrch(ind), holdvoc(i) ENDDO CLOSE(60) CALL SYSTEM('sort -r -n fort.60 > sortlist.60') OPEN(60,file='sortlist.60') DO i=1,nhldvoc READ(60,'(f20.10,2x,a132)') dummy, holdvoc(i) ENDDO CLOSE(60) READ(holdvoc(1),'(a126,i3,i3)') rdct,stabl,level ENDIF ENDIF name = rdct(1:lco) chem = rdct(lco+1:lcf) nc = index(chem,' ') - 1 irewind=irewind+1 * screen writing decreases efficiency by about 30%. Write to file 35 * and clean every 100 species written. Uncomment screen writing * if debugging is required * every species ! DEBUG ! ! IF(nrectot.EQ.3373)THEN ! GOTO 987 ! ENDIF ! DEBUG! ! WRITE(6,'(i8,1x,i6,1x,i3,1x,a)') ! & nrectot,nhldvoc,stabl,chem(1:nc) ! IF(wtflag.GT.0)THEN ! WRITE(6,*) '..........................................' ! WRITE(6,'(i8,1x,i6,1x,i3,1x,a)') ! & nrectot,nhldvoc,stabl,chem(1:nc) ! ELSE ! IF(nhldvoc.EQ.1.AND.stabl.EQ.0)THEN !! WRITE(6,*) '..........................................' ! WRITE(6,'(i8,1x,i6,1x,i3,1x,a)') ! & nrectot,nhldvoc,stabl,chem(1:nc) ! ENDIF ! ENDIF IF (wtflag.NE.0) THEN WRITE(35,'(i8,1x,i6,1x,i3,1x,a)') & nrectot,nhldvoc,stabl,chem(1:nc) WRITE(6,*)'=================================================' WRITE(6,'(i8,1x,i6,1x,i3,1x,a)') & nrectot,nhldvoc,stabl,chem(1:nc) WRITE(6,*)'=================================================' ENDIF * every 'nrewind'th species IF (irewind.eq.nrewind) THEN irewind=0 REWIND(35) IF(wtflag.EQ.0)THEN WRITE(6,'(i8,1x,i6,1x,i3,1x,a)') & nrectot,nhldvoc,stabl,chem(1:nc) ENDIF ENDIF IF(nhldvoc.LT.0)STOP * write information about depth (species in stack vs depth) c IF (stabl.gt.ndepth) THEN c WRITE(82,*) ndepth+1,nrectot,nhldvoc c ndepth=stabl c ENDIF * If the species is a special species then don_t go to the * reaction subroutine. Chemistry must be read in the "oge..." tables. * special chemistry for aromatics IF (chem(1:1).eq.'#') THEN CALL spreact(rdct,brch, & dbrch,dict,namlst,nfn,namfn,chemfn, & noge,ogertve,ogeprod,ogearh,ogestoe,ogelab,ogeaux, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl) IF (INDEX(chem,'.').EQ.0)THEN CALL write_partitionning(chem,name,temp,chatab,ncha) ENDIF ! JMLT test: convert peroxy radicals to corresponding hydroperoxide IF (ro2cond_fg.eq.1.AND.INDEX(chem,'(OO.)').NE.0) THEN tchem = chem j = INDEX(tchem,'(OO.)') IF(j.NE.0) tchem(j:j+4)='(OOH)' CALL write_partitionning(tchem,name,temp,chatab,ncha) tchem = chem ENDIF holdvoc(1:nhldvoc)=holdvoc(2:nhldvoc+1) nhldvoc = nhldvoc-1 GOTO 1000 ENDIF * rebuild the info linked to the species (bond matrix, ...) and * check that the species is effectively known in the dictionary * (just in case). ic = cnum (chem,nc) CALL grbond(chem,nc,group,bond,dbflg,nring) ind = srch(nrec,chem,dict) IF (ind.LT.0) THEN mesg ='species loaded from stack is unknown in dictionary.' CALL errexit(progname,mesg,chem(1:nc)) DO i=1,nrec WRITE(18,'(a)') dict(i) ENDDO STOP ENDIF ************************************************************ ! last stable generation to react - get atom numbers for mass budget. ! warning, information stored in the common block (see common.h) * get the branching ratio for the species currently treated brch = dbrch(ind) IF (brch.LT.1E-20) brch=1E-20 iflost=0 !IF (.not.lorad) THEN IF (.NOT.lorad.OR. & (ro2cond_fg.EQ.1.AND.INDEX(chem,'(OO.)').ne.0)) THEN tc = 1 IF (chem(1:3).EQ.'#mm' ) THEN tchem=' ' tchem(1:)=chem(4:) ELSE IF (chem(1:3).EQ.'#bb' ) THEN tchem=' ' tchem(1:)=chem(4:) ELSE IF (chem(1:1).EQ.'#' ) THEN tchem=' ' tchem(1:)=chem(2:) ELSE tchem=' ' tchem=chem ENDIF nc = index(tchem,' ') - 1 CALL number(tchem,nc,tc,ih,in,io,ir,is,ifl,ibr,icl) IF ( stabl.GE.maxgen ) THEN iflost=1 CALL number(tchem,nc,xxc,xxh,xxn,xxo,xxr,xxs,xxfl,xxbr,xxcl) ENDIF ENDIF *********************************************************** * compute the vapor pressure - exit if below threshold ! IF (.not.lorad) THEN ! JMLT test: convert peroxy radicals to corresponding hydroperoxide IF (.NOT.lorad.OR. & (ro2cond_fg.EQ.1.AND.INDEX(chem,'(OO.)').ne.0)) THEN tchem = chem j = INDEX(tchem,'(OO.)') IF(j.NE.0) tchem(j:j+4)='(OOH)' ! end test section nc = index(chem,' ') - 1 !CALL grbond(chem,nc,group,bond,dbflg,nring) CALL grbond(tchem,nc,group,bond,dbflg,nring) IF (aeropart_fg.NE.0.OR.wallpart_fg.NE.0) THEN IF (pvap_fg.EQ.1) THEN CALL molarw(tchem,mw) CALL myrdaldata(tchem,mw,Tb,HBN,tau) log10Psat= - (86. + 0.4*tau + 1421*HBN) & *(Tb-temp)/(19.1*temp) & + (-90.0-2.1*tau)/19.1 & *((Tb-temp)/temp-log(Tb/temp)) ! JMLT: reset group identities for chemistry routines IF(j.NE.0) THEN tchem = chem CALL grbond(tchem,nc,group,bond,dbflg,nring) ENDIF ELSE IF (pvap_fg.EQ.2) THEN CALL nannoolal_tb(tchem,group,bond,nring,Tb,Nangroup) CALL nannoolal_pvap(tchem,temp,Tb,Nangroup,dB,log10Psat) Trb=temp/Tb log10Psat = (4.1012+dB)*((Trb-1.)/(Trb-(1./8.)))! write mass transfer equation !write aerosol and wall partitioning reactions CALL write_partitionning(tchem,name,temp,chatab,ncha) ! JMLT: reset group identities for chemistry routines IF(j.NE.0) THEN tchem = chem CALL grbond(chem,nc,group,bond,dbflg,nring) ENDIF ELSE IF (pvap_fg.EQ.3) THEN CALL simpol1(tchem,group,bond,nring,simpgroup) ! JMLT: reset group identities for chemistry routines IF(j.NE.0) THEN tchem = chem CALL grbond(chem,nc,group,bond,dbflg,nring) ENDIF log10Psat=0. DO i=0,30 IF (simpgroup(i).NE.0) THEN log10Psat = log10Psat + simpgroup(i)* & ((bk(i,1)/temp)+bk(i,2)+(bk(i,3)*temp)+(bk(i,4)*alog(temp))) ENDIF ENDDO CALL write_partitionning(chem,name,temp,chatab,ncha) ELSE WRITE(6,*) 'aerosol condensation allowed BUT' WRITE(6,*) 'no vapour pressure estimation method selected' WRITE(6,*) 'see main.f' STOP ENDIF ENDIF WRITE(56,*) chem(1:60),' ',exp(logPvap*log(10.)) IF (log10Psat.lt.critvp) THEN WRITE(55,*) chem(1:60),' ',log10Psat !WRITE(6,*) chem(1:60),' ',log10Psat !holdvoc(1:nhldvoc)=holdvoc(2:nhldvoc+1) !nhldvoc = nhldvoc-1 ! JMLT RO2 condensation requires updating the RADICAL list size IF(.NOT.lorad)THEN holdvoc(1:nhldvoc)=holdvoc(2:nhldvoc+1) nhldvoc = nhldvoc-1 ELSE holdrad(1:nhldrad)=holdrad(2:nhldrad+1) nhldrad = nhldrad-1 ENDIF ! END JMLT GO TO 1000 ENDIF ENDIF * ------------------------------------- * SEND SPECIES TO REACTION SUBROUTINES * ------------------------------------- * STABLE VOC * =========================== IF (.not.lorad) THEN !! DEBUG !! ! write(*,*)'Press Enter to continue' ! read(*,*) !! END DEBUG !! !$OMP PARALLEL SECTIONS !$OMP+ FIRSTPRIVATE(bond,group,nring) !$ IF(wtflag.NE.0) print*,'#thread',omp_get_thread_num(),'PAN decomp' * Conversion to dihydrofuran IF (flgdhf.EQ.1) THEN IF (((INDEX(chem, 'CHO').NE.0).OR.(INDEX(chem, 'CO').NE.0)) & .AND.(INDEX(chem, 'OH').NE.0)) THEN cut_off = cutoff_default IF (masstransfg.EQ.1) THEN CALL dhf(rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn,fdhf) ELSE CALL dhf_thf(rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn,fdhf, & ncha,chatab) ENDIF cc PRINT *, 'fin' c! test si hydrofuran trouve - exit si oui c IF (fdhf.ne.0) GOTO 55 ENDIF ENDIF * PAN decomposition IF (index(chem,'OONO2').ne.0) THEN cut_off = cutoff_PAN ! non-PAN peroxy-nitrate ONLY decomposes to RO2+NO2, even if the ! molecule has a PAN IF ((index(chem,'C(OONO2)').ne.0).OR. & (index(chem,'c(OONO2)').ne.0).OR. & (index(chem,'H(OONO2)').ne.0).OR. & (index(chem,'1(OONO2)').ne.0).OR. & (index(chem,'2(OONO2)').ne.0)) THEN CALL decomp_RO2NO2(rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) GOTO 55 ! normal PAN treatment in absence of non-PAN peroxy-nitrate ELSE CALL decomp(rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF ! no further reactions for non-PAN species IF (index(chem,'CO(OONO2)').eq.0) GOTO 55 ENDIF !$OMP SECTION !$ IF(wtflag.NE.0) print*,'#thread',omp_get_thread_num(),'OH reac' * reaction with OH IF ((INDEX(chem,'H').ne.0).OR.(INDEX(chem,'Cd').ne.0)) THEN c IF (stabl.EQ.0) THEN c cut_off=0.01 c ELSE cut_off=cutoff_OH c ENDIF CALL ho_rad2(rdct,bond,group,nring,brch,temp, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nohdat,ohchem,ohrate, & nfn,namfn,chemfn, & nkwoh,nkwohpd,kwohrct,kwohpd,kwohcopd,kwohyld) ENDIF c GOTO 987 !$OMP SECTION !$ IF(wtflag.NE.0) print*,'#thread',omp_get_thread_num(),'O3 reac' * reaction with O3 IF (INDEX(chem,'Cd').ne.0) THEN cut_off=cutoff_O3 CALL rozone2(chem,rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & no3dat,o3chem,o3rate,nfn,namfn,chemfn, & nkwo3,nkwo3pd,kwo3rct,kwo3pd,kwo3copd,kwo3yld) ENDIF !$OMP SECTION !$ IF(wtflag.NE.0) print*,'#thread',omp_get_thread_num(),'NO3 reac' * reaction with NO3 IF (zero_NOxfg.LT.1) THEN cut_off=cutoff_NO3 IF((INDEX(chem,'H').ne.0).OR.(INDEX(chem,'Cd').ne.0))THEN CALL no3rad(rdct,bond,group,nring,brch, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nno3dat,no3chem,no3rate,nfn,namfn,chemfn, & nkwno3,nkwno3pd,kwno3rct,kwno3pd,kwno3copd,kwno3yld) ENDIF ENDIF 210 CONTINUE !$OMP SECTION !$ IF(wtflag.NE.0) print*,'#thread',omp_get_thread_num(),'HV reac' * test for chromophores and perform reaction with HV rdcttmp=rdct IF (nring.gt.0) THEN CALL rjsrm(nring,rdct(lco+1:lcf),rjg) ENDIF chromophores = (INDEX(rdct(lco+1:lcf),'COC').NE.0) .OR. & (INDEX(rdct(lco+1:lcf),'CO(OOH)').NE.0) .OR. & (INDEX(rdct(lco+1:lcf),'CO(ONO2)').NE.0) .OR. & ( (INDEX(rdct(lco+1:lcf),'CO(').NE.0) .AND. & (INDEX(rdct(lco+1:lcf),'CO(O').EQ.0) ).OR. & (INDEX(rdct(lco+1:lcf),ketene).NE.0) .OR. & (INDEX(rdct(lco+1:lcf),aldehyde).NE.0) .OR. & (INDEX(rdct(lco+1:lcf),hydro_peroxide).NE.0) .OR. & (INDEX(rdct(lco+1:lcf),nitrate).NE.0).OR. & (INDEX(rdct(lco+1:lcf),pan).NE.0).OR. & (INDEX(rdct(lco+1:lcf),carboxylic_acid).NE.0) rdct=rdcttmp c IF ((chromophores).AND.(nring.EQ.0)) THEN IF (chromophores) THEN cut_off=cutoff_HV CALL hvdiss2(rdct,bond,group,nring,brch, & nbson,bsongrp,bsonval, & dbrch,dict,namlst, & cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & jchem,jprod, & coprodj,jlabel,jdat,j40,jlab40,nfn,namfn,chemfn) ENDIF !$OMP END PARALLEL SECTIONS 55 CONTINUE * remove the species from the stack holdvoc(1:nhldvoc)=holdvoc(2:nhldvoc+1) nhldvoc = nhldvoc-1 * load the next species in stacks (label 1000) GOTO 1000 ENDIF * RADICAL VOC * ============================= IF (lorad) THEN !! DEBUG !! ! write(*,*)'Press Enter to continue' ! read(*,*) !! END DEBUG !! ! Autooxidation of peroxy radicals IF (autoox_fg.EQ.1) THEN IF ((name(1:1).EQ.'2').AND.INDEX(chem,'CHO').ne.0) THEN CALL autoox(rdct,bond,group,nring,dict,namlst, & dbrch,nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn,oxdone) ! remove the species from the stack if the autooxidation was done IF (oxdone.GT.0) THEN holdrad(1:nhldrad)=holdrad(2:nhldrad+1) nhldrad = nhldrad-1 GOTO 1000 ! go to next species ENDIF ENDIF ENDIF * alkoxy radical IF (name(1:1).EQ.'1') THEN cut_off=cutoff_RO CALL ro(rdct,bond,group,nring,brch,temp,tbm, & nbson,bsongrp,bsonval, & dbrch,dict,namlst,cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn, & nrokr,rokrct,rokprd,roarrh,rocost) * peroxy radical ELSE IF (name(1:1).EQ.'2') THEN cut_off=cutoff_RO2 CALL ro2(rdct,bond,group,nring,brch,temp,tbm, & dbrch,dict,namlst,cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn, & nro2kr,ro2krct,ro2kprd,ro2arrh,ro2cost) * acyl peroxy radical ELSE IF (name(1:1).EQ.'3') THEN cut_off=cutoff_RCOO2 CALL rcoo2(rdct,bond,group,nring,brch, & dbrch,dict,namlst,cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn, & nrcoo2kr,rcoo2krct,rcoo2kprd,rcoo2arrh,rcoo2cost) ! criegee radical ELSE IF (name(1:1).EQ.'4') THEN ! cut_off=cutoff_criegee IF ( stab_criegee_fg .eq. 2) THEN CALL stabilized_criegee_cmv(rdct,bond,group,nring,brch, & dbrch,dict,namlst,cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ELSE CALL stabilized_criegee(rdct,bond,group,nring,brch, & dbrch,dict,namlst,cut_off, & nhldvoc,holdvoc,nhldrad,holdrad,level,stabl, & nfn,namfn,chemfn) ENDIF ELSE WRITE(6,*) 'problem in main. Cannot find the radical type' WRITE(6,*) 'for a following species loaded from the stack ' WRITE(6,*) 'name=', name WRITE(6,*) 'chem=', chem(1:nc) WRITE(99,*) 'main',chem ENDIF * remove the species from the stack and load next stack species (1000) holdrad(1:nhldrad)=holdrad(2:nhldrad+1) nhldrad = nhldrad - 1 GOTO 1000 ENDIF * load the next "primary" species in the input table 421 CONTINUE * ---------------- * END OF REACTIONS * ---------------- * this continue is useful for debugging (when dictionary writing is * needed 987 CONTINUE * write the total rate of production for each species in the mechanism * do it before sorting the name in dict WRITE(6,*)'dictionary writing ...' DO i=1,ninorg WRITE(7,'(a)') inorglst(i) ENDDO DO i=2,nrec WRITE(7,'(a)') dict(i) ENDDO WRITE(7,'(a)') '************' CLOSE(7) DO i=2,nrec IF (INDEX(dict(i),'.').eq.0) THEN WRITE(8,'(f15.12,2x,a)') dbrch(i),dict(i) ENDIF ENDDO WRITE(*,'(a)') '************' CLOSE(8) ! create more detailed dictionnary. ! for easier extraction ! OPEN(22, file="outdat.gasdict") ! WRITE(22,*) "name, chem, nC, nH, nO, nN, molW, functions" * --------------------------------------------------------------------- * output for the "chemical interpreter" * --------------------------------------------------------------------- DO i=1,ninorg IF ((INDEX(inorglst(i),'(M)').eq.0).and. & (INDEX(inorglst(i)(1:2),'M ').eq.0).and. & (INDEX(inorglst(i),'HV').eq.0) ) THEN WRITE(21,'(A1,A6,10X,A4)') 'G',inorglst(i),'/1./' ! WRITE(22,*) 'G',inorglst(i)(1:6),",", trim(inorglst(i)(10:129)), ! & ",NA,NA,NA,NA,NA,NA" ENDIF ENDDO ! DO i=2,nrec ! WRITE(21,'(A1,A6,10X,A1,F6.1,A1)') ! & 'G',dict(i),'/',molmass(i) ,'/' ! ENDDO * --------------------------------------------------------------------- * compute molecular parameters: * 1) molecular weight (default is 1) * 2) parameters for vapor pressure estimation, * 3) heats of formation of non-radicals * 4) parameters for Henry's Law estimation * --------------------------------------------------------------------- ! DO THIS AFTER GENERATION, IN SCRIPT gen_package_* ! WRITE(6,*)'...assign peroxy classes' ! CALL count_ro2(dict,nrec) WRITE(6,*)'...compute molecular weight' IF(postflag.GT.0) & WRITE(6,*)'...compute data for Psat & Henry evaluation' ! IF(wrt_pvap.GT.0) WRITE(6,*)'...compute data for Psat evaluation' ! IF(wrt_dhf.GT.0) WRITE(6,*)'...compute heats of formation' ! IF(wrt_henry.GT.0)WRITE(6,*)'...compute data for Henry parameters' * set default molecular mass molmass=1. * open file to store Pvap parameters IF(postflag.GT.0) THEN OPEN(76, file='difvol.dat') IF (pvap_fg.EQ.1) THEN OPEN(77, file='pmy.dat') WRITE(77,'(a6)') 'MYRDAL' ENDIF !(pvap_fg.EQ.1) THEN OPEN(78, file='pnan.dat') WRITE(78,'(a6)') 'NANNOO' OPEN(79, file='psim.dat') WRITE(79,'(a6)') 'SIMPOL' * open file to store dHf parameters IF(postflag.GT.0) OPEN(80, file='dHeatf.dat') ! IF(wrt_dhf.GT.0) OPEN(80, file='dHeatf.dat') * open file to store Henry parameters IF(postflag.GT.0) OPEN(81, file='henry.dat') ! IF(wrt_henry.GT.0) OPEN(81, file='henry.dat') ! open file to store dimerisation data IF(dimer_fg.EQ.1) THEN OPEN(83, file='dimerisation.log') OPEN(84, file='indat1.dim') OPEN(85, file='indat2.dim') OPEN(86, file='indat3.dim') OPEN(87, file='indat4.dim') ENDIF * loop through species DO i=2, nrec nc = 0 ic = 0 ih = 0 in = 0 io = 0 ir = 0 is = 0 ifl = 0 ibr = 0 icl = 0 chem=dict(i)(10:129) IF (chem(1:3).eq."#mm") THEN chem = chem(4:) ELSE IF (chem(1:3).eq."#bb") THEN chem = chem(4:) ELSE IF (chem(1:1).eq."#") THEN chem = chem(2:) ENDIF nc = index(chem,' ') - 1 fgrp=trim(dict(i)(132:)) CALL molarw_number_withexceptions(chem,molmass(i), flg, & nc,ic,ih,in,io,ir,is,ifl,ibr,icl) IF (ic.le.1 .or. flg) THEN WRITE(21,'(A1,A6,10X,A1,F9.4,A1)') & 'G',dict(i),'/',molmass(i) ,'/' ! write the C1 in the fort.21 file ! DO NOT WRITE PARTITIONING SPECIES FOR C1 ! IF (aeropart_fg .eq. 1) THEN ! WRITE(21,'(A1,A6,10X,A1,F9.4,A1)') ! & 'A',dict(i),'/',molmass(i) ,'/' ! write the C1 in the fort.21 file ! ENDIF ! IF (wallpart_fg .eq. 1) THEN ! WRITE(21,'(A1,A6,10X,A1,F9.4,A1)') ! & 'W',dict(i),'/',molmass(i) ,'/' ! write the C1 in the fort.21 file ! ENDIF CYCLE ! C1 ENDIF ! C1 *-1) molecular weight WRITE(74,*) chem(1:60),' ',molmass(i) WRITE(21,'(A1,A6,10X,A1,F6.1,A1)') & 'G',dict(i),'/',molmass(i) ,'/' ! if a radical or a C1, don't output partitioned-species mass !IF (INDEX(chem,'.').ne.0.OR.ic.LE.1) CYCLE ! no radicals or C1 IF (ic.LE.1) CYCLE ! no radicals or C1 !IF (INDEX(chem,'.').eq.0.OR.INDEX(chem,'(OO.)').NE.0) THEN IF (INDEX(chem,'.').eq.0.OR. & (ro2cond_fg.EQ.1.AND.INDEX(chem,'(OO.)').ne.0)) THEN IF (aeropart_fg .eq. 1) THEN WRITE(21,'(A1,A6,10X,A1,F6.1,A1)') & 'A',dict(i),'/',molmass(i) ,'/' ! write the C1 in the fort.21 file ENDIF IF (wallpart_fg .eq. 1) THEN WRITE(21,'(A1,A6,10X,A1,F6.1,A1)') & 'W',dict(i),'/',molmass(i) ,'/' ! write the C1 in the fort.21 file ENDIF ENDIF *-2a) Get parameters for Pvap estimates based on M&Y methods (boiling * points Tb, Hydrogen Bond Number HBN, effective number of torsional * bonds tau). ! IF(wrt_pvap.gt.0)THEN IF(postflag.GT.0)THEN !IF (INDEX(chem,'.').eq.0) THEN ! JMLT test: calculate also for peroxy radicals IF (INDEX(chem,'.').eq.0.OR. & (ro2cond_fg.EQ.1.AND.INDEX(chem,'(OO.)').ne.0)) THEN ! JMLT test: convert peroxy radicals to corresponding hydroperoxide tchem = chem j = INDEX(tchem,'(OO.)') IF(j.NE.0) tchem(j:j+4)='(OOH)' IF (pvap_fg.EQ.1) THEN CALL myrdaldata(tchem,molmass(i),Tb,HBN,tau) WRITE(77,'(A1,A6,2x,f6.1,2x,f8.4,2x,f7.4,2x,f6.1)') & 'G',dict(i)(1:6),Tb,HBN,tau ELSE IF (pvap_fg.GT.1) THEN *-2bc) Write BOTH Nannoolal and SIMPOL parameters out *-2b) Get parameters for Pvap estimates based on Nannoolal methods !ELSE IF (pvap_fg.EQ.2) THEN CALL grbond(tchem,nc,group,bond,dbflg,nring) CALL nannoolal_tb(tchem,group,bond,nring,Tb,Nangroup) CALL nannoolal_pvap(tchem,temp,Tb,Nangroup,dB,logPvap) WRITE(78,'(A1,A6,2x,f6.1,2x,f8.4)') & 'G',dict(i)(1:6),Tb,dB *-2c) Get parameters for Pvap estimates based on SIMPOL methods !ELSE IF (pvap_fg.EQ.3) THEN CALL grbond(tchem,nc,group,bond,dbflg,nring) CALL simpol1(tchem,group,bond,nring,simpgroup) WRITE(79,'(A1,A6,2x,31(f3.0))') & 'G',dict(i)(1:6),(simpgroup(j),j=0,30) *-2d) Write diffusion volumes to file on unit 76 CALL calc_difvol(tchem,name,vdmol) WRITE(76,'(A1A6,2x,f7.2)') 'G',dict(i)(1:6),vdmol ENDIF ! (pvap_fg) tchem = chem *-3) get heats of formation ! IF(wrt_dhf.gt.0)THEN dHeatf = HEAT(chem,nbson,bsongrp,bsonval) WRITE(80,'(A1,A6,2x,f10.2)') 'G',dict(i)(1:6),dHeatf ENDIF ! non-radicals *-4) get Henry parameters IF (INDEX(chem,'.').eq.0) THEN CALL gromhe(chem,keff,nhenrydat,henrydat,chemhenry, & nhydratdat,hydratdat,chemhydrat) rf=0. IF (INDEX(chem,'(OOH)').NE.0) rf=0.1 IF (INDEX(chem,'CO(OONO2)').NE.0) rf=0.1 cf= (molmass(i)/18.01)**0.5 Keff=10.**Keff ! Define upper limit for Henry's contant Keff = MIN(Keff,1.0E18) WRITE(81,'(E7.1,2x,E7.1,2x,E7.1,2x,a1,a6,a1)') & cf,Keff,rf,'G',dict(i)(1:6),' ' ELSE IF (INDEX(chem,'(OO.)').ne.0.AND.ro2dep_fg.EQ.1) THEN rf=0. cf= (molmass(i)/18.01)**0.5 Keff=1.0E+14 WRITE(81,'(E7.1,2x,E7.1,2x,E7.1,2x,a1,a6,a1)') & cf,Keff,rf,'G',dict(i)(1:6),' ' ENDIF ENDIF ! (postflag) IF (INDEX(chem,'.').NE.0) CYCLE !-5) count aldehyde and hydroperoxide for dimerisation name=dict(i)(1:lco) CALL compteur(name,chem, & nalde,nalco,nooh,ncooh,ncoooh, & molmass(i), & malde,malco,mooh,mcooh,mcoooh) !-------- end species loop ENDDO WRITE(6,*)'...done compteur' IF (dimer_fg.EQ.1) THEN malde=malde/nalde malco=malco/nalco mooh=mooh/nooh mcooh=mcooh/ncooh WRITE(83,*)' number of aldehyde =',nalde WRITE(83,*)' number of alcohol =',nalco WRITE(83,*)' number of hydroperoxide =',nooh WRITE(83,*)' number of carboxylic acid =',ncooh WRITE(83,*)' M aldehyde =',malde WRITE(83,*)' M alcohol =',malco WRITE(83,*)' M hydroperoxide =',mooh WRITE(83,*)' M carboxylic acid =',mcooh WRITE(84,'(a5)') '*****' CLOSE(84) WRITE(85,'(a5)') '*****' CLOSE(85) WRITE(86,'(a6)') '*****' CLOSE(86) WRITE(87,'(a6)') '*****' CLOSE(87) ENDIF WRITE(6,*)'...done molwt' IF(postflag.gt.0)THEN ! IF(wrt_pvap.gt.0)THEN WRITE(6,*)'...done Psat' DO i=76,79 WRITE(i,'(a4)') 'END ' CLOSE(i) ENDDO ENDIF ! ENDIF ! IF(wrt_dhf.gt.0)THEN WRITE(6,*)'...done dHeatf' WRITE(80,'(a4)') 'END ' CLOSE(80) ! ENDIF ! IF(wrt_henry.gt.0)THEN WRITE(6,*)'...done Henry' WRITE(81,'(a4)') 'END ' CLOSE(81) CALL SYSTEM('cat ../DATA/henry.dep.clean henry.dat > henry.dep') CALL SYSTEM('mv henry.dep henry.dat') ! ENDIF ENDIF ! (postflag) ! diagnostic output for partitioning species IF ((wtflag.NE.0).AND.(pvap_fg.EQ.2).AND.(masstransfg.EQ.2)) THEN OPEN(82, file='chaname.dat') DO i=1,ncha WRITE(82,'(A1,A6)') 'G',chatab(i) ENDDO WRITE(6,*)'...done CHA' CLOSE(82) ENDIF ! DIMERISATION IF (dimer_fg.EQ.1) THEN WRITE(6,*)'...compute dimerisation' ! loop through species DO i=2, nrec name=dict(i)(1:lco) chem=dict(i)(10:129) nc = index(chem,' ') - 1 IF (cnum(chem,nc).le.1) CYCLE ! C1 * remove species not provided by a chem formula (eg in special dict). tchem=' ' IF (chem(1:1).EQ.'C'.or.chem(1:1).EQ.'c'.or. & chem(1:2).EQ.'-O' ) THEN tchem=chem ELSE IF (chem(1:3).EQ.'#mm' ) THEN tchem(1:)=chem(4:) ELSE IF (chem(1:3).EQ.'#bb' ) THEN tchem(1:)=chem(4:) ELSE WRITE(6,*) chem,': cycling' CYCLE ENDIF ENDIF ! remove radicals IF (INDEX(chem,'.').ne.0) CYCLE ! remove radicals ! check if the current species has an aldehyde, if yes, make it react ! with en generic hydroperoxyde to form a dimer, and remove one to the ! hydroperoxide counter ! last argument is an id which gives the type of reaction : ! 1 for reaction with -CHO ! 2 for reaction with -OOH ! 3 for reaction with -OH ! 4 for reaction with -CO(OH) lodimer=.false. IF (INDEX(chem,'CHO').NE.0) THEN CALL dimerisation(17,21,name,tchem,nooh,molmass(i),2,lodimer) CALL dimerisation(17,21,name,tchem,nalco,molmass(i),3,lodimer) ELSE IF (INDEX(chem,'(OOH)').NE.0) THEN CALL dimerisation(17,21,name,tchem,nalde,molmass(i),1,lodimer) ! CALL dimerisation(17,21,name,tchem,nalco,molmass(i),3,lodimer) ELSE IF (INDEX(chem,'(OH)').NE.0) THEN CALL dimerisation(17,21,name,tchem,nalde,molmass(i),1,lodimer) ! CALL dimerisation(17,21,name,tchem,nooh,molmass(i),2,lodimer) CALL dimerisation(17,21,name,tchem,ncooh,molmass(i),4,lodimer) ENDIF IF (lodimer) THEN WRITE(21,'(A1,A6,10X,A1,F6.1,A1)') & 'D',name,'/',molmass(i) & ,'/' ! WRITE(22,*) 'D',dict(i)(1:6),",",trim(tchem),",",ic,",",ih,",",io, ! & ",", in, ",", molmass(i),",", fgrp ENDIF ENDDO WRITE(6,*)'...done dimerisation' ENDIF !(dimer_fg) CLOSE(83) WRITE(21,'(a3)') 'END' WRITE(21,'(a)') '!_______________________________________________' WRITE(21,'(A9,11X,A9,48X,A7)') 'REACTIONS','MOLECULES','KELVINS' WRITE(21,'(a)') '!_______________________________________________' CLOSE(21) CLOSE(15) CLOSE(16) WRITE(17,'(A3)')'END' IF (wtopeflag.EQ.1) WRITE(10,'(A3)')'END/' CLOSE(17) * ---------------- * WRITE CPU * ---------------- tdat(1)=0. tdat(2)=0. durat=etime(tdat) WRITE (6,'(a14, 1pe11.4)') 'duration=', durat CALL date_and_time(date,endtime) WRITE(6,*) 'begtime=',begtime WRITE(6,*) 'endtime=',endtime elapsed=time_diff(begtime,endtime) WRITE(6,*) 'telstd =',telstd WRITE(6,*) 'elapsed=',elapsed ! report on "no-sink" species listed in fort.50 ! WRITE(50,*) ' ' CLOSE (50) INQUIRE(FILE="fort.50",EXIST=flg) IF(flg) CALL system('cat fort.50 ') END