CMAQ-DDM-3D User's Guide - 6 January 2010 Sergey L. Napelenok U.S. EPA napelenok.sergey@epa.gov Daniel S. Cohan Rice University cohan@rice.edu TABLE OF CONTENTS 1. Overview of Direct Decoupled Method 2. Execution of CMAQ-DDM 3. Compatible CMAQ Options 4. Output of Data 5. Shortcomings and Unimplemented Features 6. Summary 7. References Appendices 1. OVERVIEW OF DECOUPLED DIRECT METHOD The decoupled direct method in 3 dimensions (DDM-3D) is a sensitivity analysis technique for computing sensitivity coefficients simultaneously while air pollutant concentrations are being computed (see References). The sensitivity coefficients represent the change in concentration, of any modeled species at any modeled time, associated with a change in a model input (e.g., an initial condition, boundary condition or emission rate) or a parameter (e.g., a reaction rate). CMAQ-DDM-3D generates concentration outputs that are essentially identical to normal CMAQ results, while simultaneously computing sensitivity coefficients for any species concentration to a change in initial conditions, boundary conditions, or emission rates. We have attempted to adhere to CMAQ programming conventions. Some changes have been made to the input and output of data to keep the sensitivity output files to a reasonable size and to facilitate the necessary passing of information from mother domains to daughter domains for nested runs (see Output of Data). Current DDM-3D implementation is available for version 4.7 of the Community Multiscale Air Quality (CMAQ) model. This version of the code has been the work of Sergey L. Napelenok, Daniel Cohan, Ted Russell, Yueh-Jiun Yang, Amir Hakami and James Boylan. This code can be compiled and executed in parallel mode identically to the base CMAQ model to take advantage of available multiprocessing capabilities. This guide is intended to assist users of our implementation of DDM sensitivity analysis into CMAQ. We acknowledge that our implementation is a work in progress and list cautionary notes (see Shortcomings and Unimplemented Features) that should be considered when using CMAQ-DDM. However, we have tested that for ozone chemistry and PM processes CMAQ-DDM gives results in good agreement with sensitivities calculated by differencing multiple brute-force runs of CMAQ, at a significant savings of computational time. 2. EXECUTION OF CMAQ-DDM The CMAQ code, as modified for sensitivity analysis, must be compiled with a modified bldit script, run using a modified runscript, and requires an additional input file specifying the parameters of interest. These required changes are as follows: bldit script You can modify the bldit.cctm.linux script as follows and build the DDM version of the CCTM: 43c43 < set APPL = e3a --- > set APPL = de3a 65c65 < set ModDriver = ( module ctm_yamo $Revision; ) --- > set ModDriver = ( module ctm_yamo_ddm3d $Revision; ) 74c74 < set ModInit = ( module init_yamo $Revision; ) --- > set ModInit = ( module init_yamo_ddm3d $Revision; ) 80c80 < set ModCpl = ( module gencoor $Revision; ) --- > set ModCpl = ( module gencoor_ddm3d $Revision; ) 84c84 < set ModHadv = ( module hyamo $Revision; ) --- > set ModHadv = ( module hyamo_ddm3d $Revision; ) 88c88 < set ModVadv = ( module vyamo $Revision; ) --- > set ModVadv = ( module vyamo_ddm3d $Revision; ) 91c91 < set ModHdiff = ( module multiscale $Revision; ) --- > set ModHdiff = ( module multiscale_ddm3d $Revision; ) 94c94 < #set ModVdiff = ( module eddy $Revision; ) --- > set ModVdiff = ( module eddy_ddm3d $Revision; ) 96c96 < set ModVdiff = ( module acm2_inline $Revision; ) --- > #set ModVdiff = ( module acm2_inline $Revision; ) 106c106 < set ModChem = ( module ebi_cb05cl_ae5 $Revision; ) --- > set ModChem = ( module ebi_cb05cl_ae5_ddm3d $Revision; ) 112c112 < set ModAero = ( module aero5 $Revision; ) --- > set ModAero = ( module aero5_ddm3d $Revision; ) 119c119 < set ModCloud = ( module cloud_acm_ae5 $Revision; ) --- > set ModCloud = ( module cloud_acm_ae5_ddm3d $Revision; ) 126a127 > set SensOpt # Uncomment to use DDM 225a227,232 > if ( $?SensOpt ) then > set DDM = -Dsens > else > set DDM = > endif > 278c285 < set text = "$quote$CPP_FLAGS $CV $PAR $STX1 $STX2 $STX3$quote;" --- > set text = "$quote$CPP_FLAGS $DDM $CV $PAR $STX1 $STX2 $STX3$quote;" Note: You can change the default scripts by using the Unix "patch" utility. Cut the indented section listed above into a file, say "mod." Then type, for example, "patch bldit.cctm.pgf mod." DDM-3D includes additional files beyond the base model and these will be retreived from the archive using the modified bldit script. The netcdf and ioapi libraries used to compile the executable must be able to accommodate the appropriate number of variables in CMAQ output files. This number is generally equal to NP*NSPCS, where NP is the number of sensitivity parameters specified and NSPCS is the number of species tracked by the chemical mechanism employed. run script The switch to calculate DDM-3D calculations must be set to true: #> ddm-3d calculations [ T | Y | F | N ] setenv CTM_DDM3D Y # note: executable must be compiled with "SensOpt" Output files locations must be specified: #> sensitivity output files if ( $?CTM_DDM3D ) then if ( $CTM_DDM3D == 'Y' || $CTM_DDM3D == 'T' ) then set SENSfile = $EXEC"_"$EPI.SENGRID.${APPL} # CTM_SENS_1 (sens equivalent of S_CGRID) set ASENSfile = $EXEC"_"$EPI.ASENS.${APPL} # A_SENS_1 set SWDfile = $EXEC"_"$EPI.SENWDEP.${APPL} # CTM_SWETDEP_1 set SDDfile = $EXEC"_"$EPI.SENDDEP.${APPL} # CTM_SDRYDEP_1 endif endif Other options and input file locations must be set: #> ddm-3d specific input files if ( $?CTM_DDM3D ) then if ( $CTM_DDM3D == 'Y' || $CTM_DDM3D == 'T' ) then set NPMAX = 1 # number of sensitivity parameters set SN_INpath = /vortex/home/nsu/CMAQ/CMAQv4_7_DDM/run_prb/sensinput set SN_INfile = sensinput.prb.saprc99_boun setenv DDM3D_HIGH N # allow higher-order sensitivity parameters [ T | Y | F | N ] #> (default is N/F) setenv DDM3D_RGN Y # allows for a file specifying regions [ T | Y | F | N ] #> (default is N/F) set it below if Y/T set REGNpath = /vortex/home/nsu/CMAQ/CMAQv4_7_DDM/run_prb/sensinput set REGNfile = prb_rgns.ncf setenv DDM3D_ES N # emissions split into categories [ T | Y | F | N ] #> (default is N/F) set them below if Y/T set EMISpath_S = ${EMISpath} set EMIS2file = #area set EMIS3file = #biog set EMIS4file = #mobi set EMIS5file = #poin set EMIS6file = #nonroad if ( $APPL == $FIRSTDAY ) then setenv DDM3D_RST Y # begins from sensitivities from a restart file [ T | Y | F | N ] #> (default is Y/T) the setting of N/F is typically for first simulation day only set SGC_ICpath = $OUTDIR1 set SGC_ICfile = CCTM_ddm_saprc99_prb.SENGRID.20031231 else setenv DDM3D_RST Y set SGC_ICpath = $OUTDIR1 set SGC_ICfile = $EXEC"_"$EPI.SENGRID.${yesterday} endif setenv DDM3D_BCS F # use sensitivity bc file [ T | Y | F | N ] #> (default is N/F) set SGC_BCpath = set SGC_BCfile = source $BASE/ddm3d.q endif endif The script also calls an external file 'ddm3d.q'. The contents of this file are provided in Appendix A. Cut, paste, and save these lines to a file named ddm3d.q into the directory with your runscripts. sensinput.dat You must specify which sensitivity parameters you want to examine by carefully listing the specifications strictly following formatting specifications. See Appendix B for instructions and an example of this input file. 3. COMPATIBLE CMAQ OPTIONS CMAQ model allows the use some flexibility in choosing various options, such as the chemical mechanism. At this time, DDM-3D has been implemented for some, but not all, of these options and care must be given to insure that the desired combination of run parameters has actually been implemented. At best, the code will not compile or the run will crash for incompatible options, but at worst, erroneous sensitivities will be computed. The following is the list of CMAQ science options released with model version 4.7.1 and the corresponding DDM-3D compatibility: Science Process CMAQ module DDM-3D ================================================== Advection: ppm NO yamo YES Vertical Diffusion: eddy YES acm2_inline NO (soon) Horizontal Diffusion: multiscale YES Chemical Mechanism: cb05cl YES saprc99 YES Chemical Solver: smvgear NO ros3 NO ebi YES Cloud Processing: cloud_acm_ae5 YES Aerosol Processing: aero4 NO aero5 YES Photolysis: phot_table YES phot_inline YES ================================================== 4. INPUT AND OUTPUT OF DATA CMAQ-DDM-3D requires the same input files as a normal CMAQ run. Additional input files may be required depending on the choice of calculated sensitivity parameters. The following table includes a list of all possible files specific to sensitivity calculations. Base model FILE Type Contains analog ============================================================================ sensinput.dat Input Sensitivity Parameter definitions (required). N/A ASENS Output Averaged hourly sensitivities. List defined by 'AVG_CONC_SPCS' variable in the run script. ACONC SENGRID Output Last hour's sensitivity fields to be used as initial conditions for the following time period. CGRID SENWDEP Output Sensitivities of wet deposited species. WETDEP1 SENDDEP Output Sensitivities of dry deposited species. DRYDEP REGIONS Input Regional definitions. N/A EMIS_A Input Area sector emissions for split emissions. N/A EMIS_B Input Biogenic sector emissions for split emissions. N/A EMIS_M Input Mobile sector emissions for split emissions. N/A EMIS_P Input Point sector emissions for split emissions. N/A EMIS_N Input Nonroad sector emissions for split emissions. N/A SBCON Input Sensitivity field boundary conditions. BCON SICON Input Sensitivity field initial conditions. ICON ============================================================================ 5. SHORTCOMINGS AND UNIMPLEMENTED FEATURES The implementation of CMAQ-DDM is a work in progress. Features not yet implemented in CMAQ-DDM include: * Tracer species * 2nd order aerosol sensitivities (in development at GaTech) * SMVGEAR and Rosenbrock solvers (This should be straightforward to do if desired, by adopting the same approach as we used in hrdriver.F. ) * Only sensitivity to emissions, initial conditions, boundary conditions, and reaction rate constants have been implemented so far; it would be straightforward to extend CMAQ-DDM to handle sensitivities to deposition velocities (in development at Rice) or wind speed (see Yang et al. 1997). 6. SUMMARY DDM has proven to be a very effective tool for air quality studies. This implementation in CMAQ has been done with the intent to provide flexibility and computational efficiency, and also maintain the base CMAQ code structure. CMAQ-DDM has been found to accurately simulate sensitivity of ozone and PM species to initial conditions, boundary conditions, and emissions of precursor species. However, CMAQ-DDM remains a work in progress with known shortcomings and its accuracy has not been tested for all conceivable applications. Any errors should be reported to the contacts provided on the cover page. 7. REFERENCES Cohan, D., Y. Hu, A. Hakami, M. T. Odman, A. Russell, 2002: Implementation of a direct sensitivity method into CMAQ. Models-3 User's Workshop, RTP, North Carolina, October 22, 2002. Available at: www.cmascenter.org/workshop/session5/cohan_abstract.pdf Cohan, D., Y. Hu, A. Hakami, A. Russell, 2003: Sensitivity Analysis of Ozone in the Southeast. Models-3 User's Workshop, RTP, North Carolina, October 27, 2003. Available at: www.cmascenter.org/2003_workshop/session2/cohan_abstract.pdf Cohan, D., Y. Hu, A. Hakami, A. Russell, 2005: Nonlinear response of ozone to emissions: source apportionment and sensitivity analysis. Environ. Sci. Technol., 39, 6739-6748. Cohan, D., D. Tian, Y. Hu, and A. Russell (2006). Control strategy optimization for attainment and exposure mitigation: Case study for ozone in Macon, Georgia. Environmental Management, 38, 451-462. Cohan, D., Y. Hu, and A. Russell (2006). Dependence of ozone sensitivity analysis on grid resolution. Atmospheric Environment, 40, 126-135. Dunker, A., 1981: Calculation of sensitivity coefficients for complex atmospheric models. Atmos. Environ., 15, 1155-1161. Dunker, A. 1984: The decoupled direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys., 81, 2385-2393. Dunker, A., G. Yarwood, J. Ortmann, and G. Wilson, 2002: The decoupled direct method for sensitivity analysis in a three-dimensional air quality model'Implementation, accuracy, and efficiency. Environ. Sci. Technol., 36, 2965-2976. Napelenok, S., D. Cohan, Y. Hu, A. Russell, 2006: Decoupled direct 3D sensitivity analysis for particulate matter (DDM-3D/PM). Atmos. Environ., 40, 6112-6121. Napelenok, S., D. Cohan, M.T. Odman, S. Tonse, 2008: Extension and evaluation of sensitivity analysis capabilities in a photochemical model. Environ. Model. & Software., 23(8), 994-999. Yang, Y-J, J. Wilkinson, and A. Russell, 1997: Fast, direct sensitivity analysis of multidimensional photochemical models. Environ. Sci. Technol, 31, 2859-2868. Appendix A --------------------------------------------------------------------- #ddm3d.q # sensitivity parameter definition (required) if ( $?SN_INpath ) then setenv SEN_INPUT $SN_INpath/$SN_INfile if (! (-e $SEN_INPUT) ) then echo " $SEN_INPUT not found " exit 1 endif endif # for regions file if ( $?DDM3D_RGN ) then if ( $DDM3D_RGN == 'Y' || $DDM3D_RGN == 'T' ) then setenv REGIONS_1 $REGNpath/$REGNfile if (! (-e $REGIONS_1) ) then echo " $REGIONS_1 not found " exit 1 endif endif endif # for split emissions files if ( $?DDM3D_ES ) then if ( $DDM3D_ES == 'Y' || $DDM3D_ES == 'T' ) then setenv EMIS_A $EMISpath_S/$EMIS2file setenv EMIS_B $EMISpath_S/$EMIS3file setenv EMIS_M $EMISpath_S/$EMIS4file setenv EMIS_P $EMISpath_S/$EMIS5file setenv EMIS_N $EMISpath_S/$EMIS6file set flist = ( $EMIS_A\ $EMIS_B\ $EMIS_M\ $EMIS_P\ $EMIS_N ) foreach file ( $flist ) if (! (-e $file) ) then echo " $file not found " exit 1 endif end endif endif # for sensitivity IC files (not always there, so no existance check - RISKY!!!) set SAE_ICpath = $SGC_ICpath set SNR_ICpath = $SGC_ICpath set SAE_ICfile = $SGC_ICfile set SNR_ICfile = $SGC_ICfile setenv INIT_GASC_S $SGC_ICpath/$SGC_ICfile setenv INIT_AERO_S $SAE_ICpath/$SAE_ICfile setenv INIT_NONR_S $SNR_ICpath/$SNR_ICfile # for sensitivity BC files (not always there, so no existance check - RISKY!!!) set SAE_BCpath = $SGC_BCpath set SNR_BCpath = $SGC_BCpath set SAE_BCfile = $SGC_BCfile set SNR_BCfile = $SGC_BCfile setenv BNDY_GASC_S $SGC_BCpath/$SGC_BCfile setenv BNDY_AERO_S $SAE_BCpath/$SAE_BCfile setenv BNDY_NONR_S $SNR_BCpath/$SNR_BCfile # sensitivity output files setenv CTM_SENS_1 "$OUTDIR1/$SENSfile -v" setenv A_SENS_1 "$OUTDIR1/$ASENSfile -v" setenv CTM_SWETDEP_1 "$OUTDIR1/$SWDfile -v" setenv CTM_SDRYDEP_1 "$OUTDIR1/$SDDfile -v" set flist = ( $A_SENS_1\ $CTM_SWETDEP_1\ $CTM_SDRYDEP_1 ) unalias rm foreach file ( $flist ) if ( $file != '-v' ) then if ( -e $file ) then echo " $file already exists" if ( $DISP == 'delete' ) then echo " Deleting $file" rm $file else if ( $DISP == 'update' ) then echo " Updating $file" else echo "*** RUN ABORTED ***" exit 1 endif endif endif end Appendix B --------------------------------------------------------------------- Preparing the input data file for Sensitivity: For each sensitivity: 1) (mandatory) The first line is the name of the sensitivity parameter; any 8-character name of the user's choosing, no leading spaces 2) (mandatory) The next line specifies the type of sensitivity (One leading space followed by 4 capitalized characters) EMIS: Emissions INIT: Initial Conditions BOUN: Boundary Conditions RATE: Reaction rate HIGH: Higher-order sensitivity. 3) (mandatory) For EMIS, INIT, or BOUN sensitivity: The term ' SPECIES' (all-cap, one leading space) must appear next. For RATE sensitivity: The term ' REACTION' (all-cap, one leading space must appear next. For HIGH sensitivity: The next 2 lines must each be one leading space followed by the name of the sensitivity to which we're taking higher order sensitivity. That name must have already been defined as the name of a sensitivity parameter. No further information should be defined for a higher-order sensitivity parameter. 4) (mandatory) For EMIS, INIT, or BOUN sensitivity: Specify one or more species. Names must have two leading spaces and then exactly match a species from $M3MODEL/include/release/saprc99_ae3_aq/GC_SPC.EXT. For RATE sensitivity: Specify one or more reactions. Names must have two leading spaces and then exactly match the _label_ from mech.def (also in RXDT.EXT). 5) (optional) The term ' AMOUNT' (all-cap, one leading space). This may be used only for emissions (EMIS). After ' AMOUNT', the next line may either be: a positive real number to set the hourly amount of emissions a negative number, followed by 24 numbers which specify the hourly emissions rates. NOTE: All numbers must be followed by two leading spaces. NOTE2: Amount is the emission rate from each gridcell of the desired region (as specified by GRIDCELL and/or REGION). NOTE3: If amount is specified, amount numbers must be given for each species listed under SPECIES. Either 1 amount number (for a continuous, constant amount) or 25 amount numbers (a negative followed by 24 hourly values) must be specified FOR EACH SPECIES. NOTE4: AMOUNT is assumed to be a ground-level emission (layer 1). The LAYER term (below) must be used to change this if the desired emitting gridcell or region is not layer 1. DDM is not currently set up to allow "AMOUNT" to be used with emissions into more than one layer. DEFAULT: sensitivity to emissions is assessed relative to 100% of the emissions inventory. (b) (optional) The term ' LAYER' (all-cap, 2 leading spaces). Can be used only following AMOUNT. If this term is used, next line must give an integer (3-leading spaces) of the layer in which emissions are released. Only one layer may be set. DEFAULT: AMOUNT is assumed to be a ground-level emission (layer 1). 6) (optional) The term ' DATE' (all-cap, one leading space). If this term is used, the next line must give an integer (2-leading spaces) saying how many dates are desired. The following line(s) must give the datenumbers that are desired, one date per line (yyyyddd format). DEFAULT: all dates will be used. 7) (optional) The term ' TIME' (all-cap, one leading space). If this term is used, the next line must contain two numbers (two leading spaces, hhmmss, space, hhmmss). The numbers specify the beginning and end time of the desired interval. To straddle midnight, the first number may be later than the second number. DEFAULT: all times will be used. 8) (optional) The following three terms are used to specify IREGION, which defines the region over which a sensitivity parameter applies. If more than one of the following choices are used, IREGION will combine the regions, gridcells and corner-regions specified. DEFAULT: If neither REGIONS nor GRIDCELLS nor CORNERS nor CIRCLES are specified, domainwide is assumed. NOTE: If multiple forms are used, they must be specified in the order REGIONS, GRIDCELLS, CORNERS, then CIRCLES. (a) The term ' REGIONS' (all-cap, one leading space). If this term is used, the next line must give an integer (2-leading spaces) number of regions to be defined. The following lines are used to specify regions from a regions netcdf file, and must exactly match the region name (2- leading spaces). (b) The term ' GRIDCELLS' (all-cap, one leading space). If this term is used, the next line must give an integer (2-leading spaces) number of gridcells to be defined. The following lines specify one gridcell per line with format __##_##_## i.e., 2 leading spaces, 2-digit column number, space, 2-digit row number, space, 2-digit level number. Use 99 for the level number to include all levels. (c) The term ' CORNERS' (all-cap, one leading space). If this term is used, the next line must give an integer (2-leading spaces) number of rectangular zones to be defined. The following line(s) must, for each zone, give 6 integers in the format ' C1 C2 R1 R2 L1 L2' (all 2 digit numbers, 2 leading spaces, 1 space between each number). The numbers represent a region bounded by these minimum and maximum columns numbers (C1,C2), row numbers (R1,R2), and layer numbers (L1,L2). (d) The term ' CIRCLES' (all-cap, one leading space). If this term is used, the next line must give an integer number of circular (or ring- shaped) zones to be defined. The following line must give the gridsize in km. The following line(s) must, for each circular/ring zone, give 4 numbers in a single line: COL_CENTER, ROW_CENTER, MINRAD(in km), MAXRAD(in km). Cells will be included in the ring iff the distance from that cell to the center cell is .GE. MINRAD, and .LE. MAXRAD. 9) A line 'END' (no leading spaces) signifies the end of the list. A line with no leading spaces and any other name is interpreted as a new sensitivity parameter (return to step 1). NOTE1: This list must be consistent with the max # of sens parameters (NPMAX) set in the runscript. NOTE2: For better understanding of how this file is read, or to modify/add features, look at sinput.F in the CMAQ-DDM code. ================================================================================ Example of sensinput.dat: EMISNOX EMIS SPECIES NO NO2 2ENOX HIGH EMISNOX EMISNOX RATE1 RATE REACTION 1 RATE8 RATE REACTION 8 RATEisOH RATE REACTION isOH RATEaOH RATE REACTION a1OH a2OH a3OH a4OH a5OH END