Corrections to the Community Multiscale Air Quality (CMAQ) model version 4.7.1 Sanford Sillman University of Michigan 21 November, 2010 The Community Multiscale Air Quality (CMAQ) model, developed by the U.S. Environmental Protection Agency, is widely used to simulate the formation, transport, photochemical transportation and ultimate removal of atmospheric pollutants. The most recent version (4.7.1) was released in June, 2010. This site reports some minor errors in the CMAQ code and provides corrected versions. These are not 'official' CMAQ releases, but are provided as a convenience to CMAQ users. The corrections are all intended to address situations in which the current version of CMAQ is likely to crash or to 'hang up' in an infinite loop. None of the changes are expected to affect model results. The corrections all involve CMAQ applications with meteorological input from the WRF model (and may result from properties of WRF). CMAQ is distributed by the Community Modeling and Analysis System (CMAS) Center at the University of North Carolina at Chapel Hill (http://www.cmascenter.org/ and http://www.cmaq-model.org/). The corrected files here are based on the original CMAQ 4.7.1 files. For convenience, both the original files (with suffix 'cmaq') and corrected versions (with suffix 'sillman') are made available here. To use the corrected files, rename each file to end with "...F,v" and copy them to the appropriate CMAQ subdirectory in place of the original files. The corrections: 1. Numerical problems in vertical advection/yamo: The CMAQ calculation of vertical advection (.../models/CCTM/src/vadv/yamo) includes a numerical iteration that in rare situations can lead to an infinite loop. When this happens the model 'hangs up' and stops making further progress. The potential infinite loop occurs in file vppm.F at index labels 66 and 77 (involving the iteration counter ICNT). In the modified version an exit has been added from the iterative loop when the counter ICNT exceeds 100. The maximum value of the iteration counter is passed to the parent subroutine (in zadvyppm.F) which identifies the date, time and location of the problem. In addition, the maximum number of iterations in zadvyppm.F (for a different iterative solution) has been increased from 31 (the default value) to 100 to reduce the frequency of crashes due to nonconvergence. The modified version also does not crash when there is a nonconvergence error in vertical advection. Instead, it prints out a warning message ('MAJOR ERROR') and continues until the cumulative number of nonconvergence errors exceeds a given threshhold (currently set at 50). Until that point the model skips the vertical advection calculation for the individual grid column and time step associated with the nonconvergence. I have found that nonconvergence errors occur very rarely (approximately once per million calculations) and do not persist for extended time periods. This have appeared only for locations in the Atlantic Ocean. Skipping vertical advection during these rare occurrences was adopted here as an interim solution. The iterative loop in vppm.F at labels 66 and 77 both appear to represent solutions to a cubic equation (FM(I)=FDN and FM(I)=FUP, where FM(I) is effectively a cubic function of the variable X). If this is the case a more efficient solution should be possible. Corrected files: vppm.F.sillman and zadvyppm.F.sillman CMAQ directory: .../models/CCTM/src/vadv/yamo 2. Numerical problems in horizontal advection/yamo: The CMAQ calculation of horizontal advection (.../models/CCTM/src/hadv/yamo) also includes a numerical iteration that can lead to an infinite loop. The potential infinite loop occurs in file hadvyppm.F at index label 101 (involving the iteration counter ICNT). In addition, I have found program crashes due to negative values for species miing ratios originating from the horizontal advection algorithm. These negative values occurred because the advection calculation was performed with Courant number greater than 1 (meaning that the calculation time step was too large). CMAQ determines the calculation time step to insure that the Courant number is less than one (in the program advstep.F in .../models/CCTM/src/driver/yamo). However, the horizontal advection algorithm includes an adjustment factor for wind velocity (ADJFAC, to insure conservation of mass). Failures occurred when the adjustmentor was large, pushing the Courant number above one. These tended to occur in locations with high wind shear. In the modified version the Courant number problem has been corrected by strengthening the criteria for setting the time step in advstep.F. The wind shear criteria parameter (HDIV_LIM) was lowered from 1.0 to 0.75. In addition, wind shear was incorporated in the maximum wind array (WIND_IDX) which is used to set the Courant number. Some other related modifications were made. In the modified version an exit has been added from the loop at statement 100 when the counter ICNT exceeds 100. This prevents the possible infinite loop in this subroutine. Corrected files: zadvyppm.F.sillman (in .../hadv/yamo) and advstep.F (in .../driver/yamo). CMAQ directory: .../models/CCTM/src/hadv/yamo and ..../driver/yamo 3. Dry deposition calculation with zero leaf area index (LAI): The CMAQ dry deposition calculation uses the leaf area index (LAI) to calculate one layer of resistance, and also the closely associated vegetation index. This calculation occurs in inline versions of CMAQ (...models/CCTM/src/vdiff/acm2_inline_txhgsim). A similar calculation also appears in the MCIP preprocessor. Both cases are associated with the file m3dry.F. Numerical problems can occur in these algorithms because the WRF model input has LAI=0 for grids located over the ocean or in ice-bound regions. The m3dry.F algorithm includes divide-by-LAI without protection for LAI being zero. This can cause the program to generate Inf or NaN values. The most recent version of MCIP (3.6) includes a correction for this problem: the MCIP file m3dry.F includes separate calculations for land and water, and the solution for water is used whenever vegetation (xveg) is zero. This appears to cover all occurrences of zero LAI. However the CMAQ algorithm still uses solutions that cannot handle zero LAI and are not bypassed. The inline m3dry.F for bidirectional Hg also includes a call to a numerical solver (ATX) for soil concentrations with control statements that call for the solution for land grids only. (A separate solution is provided for grids representing water.) The solution for soil concentrations can fail to converge or generate negative values if leaf area index and other soil parameters are zero or near-zero. Normally, LAI is 0.1 or higher for land grids. However the WRF output includes a small number of grids that are flagged as 'land' grids even though LAI = 0 and the soil type category is 'water'. (These have been found for grids at the southern end of Hudsons Bay in Canada.) In order to avoid these numerical crashes the following changes have been made: CMAQ inline version (m3dry.F and HGSIM.F): When leaf area index (LAI) equals zero, it is replaced with a small nonzero value to prevent divide-by-zero errors. In addition, zero LAI and soil type category equal to water, ice or bedrock both cause the model to invoke the 'water' calculation of dry deposition rather than the 'soil' calculation. (Originally the 'water' calculation was used only for grids that were flagged as water with the flag lwmask.) (MCIP version: In MCIP3.6 various protections against zero leaf area index, or xlai, are included. However resistcalc.f90 may still have a divide by zero xlai. This has not yet been included here.) MCIP version 3.4.1: Here, protection against divide-by-zero in cases with zero leaf area index has been added (for example, see 'rstom'). This occurs in m3dry.F and resistcalc.F. In addition, a bypass was added to set dry deposition equal to zero and continue (rather than crashing) when certain errors occur (usually associated with very low leaf area index). (In m3dry.F) It was found that the parameter 'badval3' in m3dry.F was occasionally negative. To prevent a crash, average values for xradyn and xrstom were used in these cases. (In m3dry.F) Corrected files: m3dry.F.hgsim_sillman and HGSIM.F.sillman CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim MCIP corrected files: version 3.4 only! These are differnt in version 3.6 and may not be needed. m3dry.F.mcip34_sillman resistcalc.F.mcip34_sillman 4. Soil type classification: mismatch between WRF/USGS and CMAQ This problem occurs only in the bidirectional soil Hg version of the model (file HGSIM.F in /vdiff/acm2_inline_txhgsim). The file HGSIM.F assigns soil parameters based on input soil type, that input represents one of 11 soil type categories (sltype from 1 to 11). However, the output from the WRF model, passed to CMAQ through MCIP, uses a soil categorization with 16 soil types (adapted from the USGS). This can result in erroneous parameter values because the 16 WRF soil types do not match the 11 soil types in HGSIM.F. In addition, grids with input soil types 12 or higher (including ocean grids) will cause CMAQ to crash at this point. The corrected version converts from the 16 soil types used by WRF to the 11 types requested by CMAQ. This correction should be used only when the meteorology is from WRF. Documentation: USGS soil classification used by WRF is identified in Table 3 at http://www.mmm.ucar.edu/wrf/users/docs/user_guide_V3.1/users_guide_chap3.htm# _Land_Use_and Corrected files: HGSIM.F.sillman CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim 5. Timing roundoff error for written output - bidirectional Hg only This error occurs in subroutine WRASX_MEDIA (in file HGSIM.F), which is called by the bidirectional Hg model to write output. The time identified in the subroutine needs to be exactly equal to the time for output, or else the program will crash. The time identified in this subroutine is occasionally different from the time for output by 1 second (for example, 065959 versus 070000). A temporary correction has been added to round the time setting to the nearest hour if the identified time is within 1 minute of the hour. (This solves the problem if the time step is 1 hour or another exact hour interval. A more complete fix would round off the time setting to the nearest timestep if within 60 seconds.) Corrected files: HGSIM.F.sillman CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim 6. Error in initial soil content for mercury - bidirectional Hg only This error occurs in subroutine OPASX_MEDIA (in file HGSIM.F), which is called by the bidirectional Hg model to set initial soil concentrations (CMEDIA) at the start of a model run. The intent of the algorithm is to set initial soil concentrations equal to 1E-30 if there is no previous-day value. However the algorithm erroneously fails to enter these values into the CMEDIA array. The CMEDIA array is assigned values inside an IF loop that is restricted to cases where previous-day input is available (see lines 990, 997, 1012, 1052 and 1060 in HGSIM.F from cmaq version 4.7.1). This is corrected by adding a loop to assign values to CMEDIA at the end of OPASX_MEDIA, after the end of the IF loop instead of before the end. Corrected files: HGSIM.F.sillman CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim 7. Toxic emissions - modification to allow missing species When multipollutant versions of CMAQ are used, including toxic pollutants, trace metals and/or mercury, the files that process inline emissions contain hard-wired species names. These files (rdemis_aetox.F and rdemis_pm_hg.F) will exit if any of the listed species are not included in the emission inventory. This is a problem if anyone wants to do a simulation that includes some toxic species but not all. To avoid this problem, modifications have been made to the rdemis files. If emissions for any of the listed species are not included in the inventory, the model now assumes that emissions are zero for the species and continues rather than crashing. Users should also note: the arrays EMIS_SPECIES and AE_TOX_SUR can be modified to match the exact species names in the input emission files. These arrays are currently consistent with the species names in the available chemistry (located at /models/include/release), but the species names can be changed by changing the array AE_TOX_SPC. Modified files: rdemis_aetox.F.txhgsim_sillman, rdemis_pm_hg.F.txhgsim_sillman nd PMEM_DEFN.F.sillman CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim 7. Another correction relevant for CMAQ 4.7.1 has been posted at: https://rc.usf.edu/trac/doc/wiki/CMAQ-4.5.1 =========================================== Summary of corrected files and their directories: To use these: (1) copy the corrected files ("...F.sillman") to the appropriate directory (2) rename the corrected file to its "use name" with the for "...F,v" for CMAQ files and "...F" for MCIP 3.4 files (3) Compile and run CMAQ using the settings for bidirectional Hg and using advection choice "yamo". (The same corrections made in these files may be applicable to other CMAQ settings, if applied to the equivalent files in parallel directories.) CMAQ directory: .../models/CCTM/src/vadv/yamo Original file Corrected file Working file name vppm.F.cmaq vppm.F.sillman vppm.F,v zadvyppm.F.cmaq zadvyppm.F.sillman zadvyppm.F,v CMAQ directory: .../models/CCTM/src/vdiff/acm2_inline_txhgsim Original file Corrected file Working file name m3dry.F.cmaq m3dry.F.hgsim_sillman m3dry.F,v HGSIM.F.cmaq HGSIM.F.sillman HGSIM.F,v rdemis_pm_hg.F.cmaq rdemis_aetox.F.txhgsim_sillman rdemis_aetox.F,v rdemis_aetox.F.cmaq rdemis_aetox.F.txhgsim_sillman rdemis_aetox.F,v PMEM_DEFN.F.cmaq PMEM_DEFN.F.sillman PMEM_DEFN.F,v MCIP directory: .../mcip3/BLD Original file Corrected file Working file name m3dry.F.mcip34 m3dry.F.mcip34_sillman m3dry.F resistcalc.F.mcip34 resistcalc.F.mcip34_sillman resistcalc.F =========================================== END