diff --git a/README.md b/README.md
index c4b5d35..6f73b53 100644
--- a/README.md
+++ b/README.md
@@ -1,8 +1,6 @@
 EMIPS：EMission Inventory Processing System
 ============================================
 
-[![Join the chat at https://gitter.im/meteoinfo/community](https://badges.gitter.im/meteoinfo/community/meteoinfo.svg)](https://gitter.im/meteoinfo/community)
-
 Installation
 ------------
 
@@ -34,6 +32,8 @@ Learn more about MeteoInfo and EMIPS in its official documentation at http://met
 Version
 -------
 
+EMIPS-1.1 was released (2023-06-24).   Added support for CMAQ and bug fixes.
+
 EMIPS-1.0 was released (2023-01-26).
 
 Author
@@ -61,7 +61,7 @@ Sci. Total Environ. 870, 161909. https://doi.org/10.1016/j.scitotenv.2023.161909
 License
 -------
 
-Copyright 2019-2022, EMIPS Developers
+Copyright 2019-2023, EMIPS Developers
 
 Licensed under the LGPL License, Version 3.0 (the "License");
 you may not use this file except in compliance with the License.
diff --git a/emips/__init__$py.class b/emips/__init__$py.class
index 65e9f72..8e12bef 100644
Binary files a/emips/__init__$py.class and b/emips/__init__$py.class differ
diff --git a/emips/__init__.py b/emips/__init__.py
index 7685688..0c2c5ad 100644
--- a/emips/__init__.py
+++ b/emips/__init__.py
@@ -4,10 +4,10 @@
 import chem_spec
 import vertical_alloc
 
-__version__ = '1.0'
+__version__ = '1.1'
 
 import os
 ge_data_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'ge_data')
 
-version = 'EMIPS-1.0 (2022.08.04)'
+version = 'EMIPS-1.1 (2023.06.24)'
 __all__ = ['ge_data_dir', 'utils', 'temp_alloc', 'spatial_alloc', 'chem_spec', 'vertical_alloc', 'version']
diff --git a/emips/chem_spec/chemical_speciation$py.class b/emips/chem_spec/chemical_speciation$py.class
index a43bf84..bf8b988 100644
Binary files a/emips/chem_spec/chemical_speciation$py.class and b/emips/chem_spec/chemical_speciation$py.class differ
diff --git a/emips/chem_spec/chemical_speciation.py b/emips/chem_spec/chemical_speciation.py
index 63d5cb0..4afdc68 100644
--- a/emips/chem_spec/chemical_speciation.py
+++ b/emips/chem_spec/chemical_speciation.py
@@ -4,7 +4,7 @@
 from .pollutant import Pollutant
 import mipylib.numeric as np
 
-__all__ = ['get_pollutant_profile', 'read_file', 'speciation', 'get_model_species_wrf']
+__all__ = ['get_pollutant_profile', 'read_file', 'speciation', 'get_model_species_WRFChem', 'get_model_species_CMAQ']
 
 
 def get_pollutant_profile(poll_profiles, pollutant):
@@ -89,9 +89,9 @@ def speciation(data, pollutant_profile):
     return s_data
 
 
-def get_model_species_wrf(mechanism_name):
+def get_model_species_WRFChem(mechanism_name):
     """
-    Get species in wrf-chem in different chemical mechanism.
+    Get species under different chemical mechanisms in WRF-Chem.
     (E_PM10: saprc99, cb05 E_PM_10: radm2, mozart)
     :param mechanism_name: (*string*) The name of chemical mechanism.
     :return: (*list*) All species(out_species) and aerosols(out_species_aer) under chemical mechanisms.
@@ -99,7 +99,7 @@ def get_model_species_wrf(mechanism_name):
     mechanism_name = mechanism_name.lower()
     if mechanism_name == 'cb05':
         ##########################################
-        # ------CB05, emiss_opt=15, (52, 15)------#
+        #------CB05, emiss_opt=15, (52, 15)------#
         ##########################################
         out_species = ['E_ACET', 'E_PAR', 'E_ALK3', 'E_ALK4', 'E_ALK5', 'E_TOL', 'E_XYL', 'E_BALD',
                        'E_ALD2', 'E_CCOOH', 'E_CO', 'E_CRES', 'E_ETH', 'E_ETHA', 'E_GLY', 'E_FORM',
@@ -113,7 +113,7 @@ def get_model_species_wrf(mechanism_name):
 
     elif mechanism_name == 'radm2':
         ##########################################
-        # ------RADM2, emiss_opt=3, (42, 19)------#
+        #------RADM2, emiss_opt=3, (42, 19)------#
         ##########################################
         out_species = ['E_ISO', 'E_SO2', 'E_NO', 'E_NO2', 'E_CO', 'E_CH4', 'E_ETH', 'E_HC3',
                        'E_HC5', 'E_HC8', 'E_XYL', 'E_OL2', 'E_OLT', 'E_OLI', 'E_TOL', 'E_CSL',
@@ -126,7 +126,7 @@ def get_model_species_wrf(mechanism_name):
                            'E_ORGJ_BB', 'E_CLI', 'E_CLJ']
     elif mechanism_name == 'saprc99':
         ###########################################
-        # -----SAPRC99, emiss_opt=13, (55, 15)-----#
+        #-----SAPRC99, emiss_opt=13, (55, 15)-----#
         ###########################################
         out_species = ['E_SO2', 'E_C2H6', 'E_C3H8', 'E_C2H2', 'E_ALK3', 'E_ALK4', 'E_ALK5', 'E_ETHENE',
                        'E_C3H6', 'E_OLE1', 'E_OLE2', 'E_ARO1', 'E_ARO2', 'E_HCHO', 'E_CCHO', 'E_RCHO',
@@ -140,7 +140,7 @@ def get_model_species_wrf(mechanism_name):
 
     elif mechanism_name == 'mozart':
         ##########################################
-        # -----MOZART, emiss_opt=10, (53, 23)-----#
+        #-----MOZART, emiss_opt=10, (53, 23)-----#
         ##########################################
         out_species = ['E_CO', 'E_NO', 'E_NO2', 'E_BIGALK', 'E_BIGENE', 'E_C2H4', 'E_C2H5OH', 'E_C2H6',
                        'E_C3H6', 'E_C3H8', 'E_CH2O', 'E_CH3CHO', 'E_CH3COCH3', 'E_CH3OH', 'E_MEK', 'E_SO2',
@@ -154,3 +154,25 @@ def get_model_species_wrf(mechanism_name):
                            'E_CO_A', 'E_ORGI_A', 'E_ORGJ_A', 'E_CO_BB', 'E_ORGI_BB', 'E_ORGJ_BB', 'E_PM_10']
 
     return out_species, out_species_aer
+    
+
+def get_model_species_CMAQ(mechanism_name):
+	"""
+    Get species under different chemical mechanisms in CMAQ.
+    :param mechanism_name: (*string*) The name of chemical mechanism.
+    :return: (*list*) All species(out_species) and aerosols(out_species_aer) under chemical mechanisms.
+    """
+	mechanism_name = mechanism_name.lower()
+	if mechanism_name == 'cb05':
+		########################################
+		#-----CB05, emiss_opt=10, (42, 19)-----#
+		########################################
+		out_species = ['ALD2', 'CO', 'ETH', 'FORM', 'ISOP', 'NH3', 'NO', 'NO2', 'UNR',
+					   'OLE', 'PAR', 'PEC', 'PMC', 'PMOTHR', 'PNO3', 'POC', 'PSO4', 'PCL', 
+					   'PNH4', 'PNA', 'PMG', 'PK', 'PCA', 'PNCOM', 'PFE', 'PAL', 'PSI', 'PTI', 
+					   'PMN', 'PH2O', 'SO2', 'SULF', 'TERP', 'TOL', 'XYL', 'MEOH', 'ETOH', 
+					   'ETHA', 'ALDX', 'IOLE', 'CH4', 'AACD']
+		out_species_aer = ['PEC', 'PMC', 'PMOTHR', 'PNO3', 'POC', 'PSO4', 'PCL', 'PNH4', 'PNA', 
+						   'PMG', 'PK', 'PCA', 'PNCOM', 'PFE', 'PAL', 'PSI', 'PTI', 'PMN', 'PH2O']
+
+	return out_species, out_species_aer
diff --git a/emips/chem_spec/species_profile$py.class b/emips/chem_spec/species_profile$py.class
index e34cc08..b666a1e 100644
Binary files a/emips/chem_spec/species_profile$py.class and b/emips/chem_spec/species_profile$py.class differ
diff --git a/emips/chem_spec/species_profile.py b/emips/chem_spec/species_profile.py
index 844b283..36c6a40 100644
--- a/emips/chem_spec/species_profile.py
+++ b/emips/chem_spec/species_profile.py
@@ -50,7 +50,12 @@ def read_string(cls, line, mechanism=None):
         data = line.split()
         pollutant = Pollutant(data[1])
         if mechanism is None:
-            species = Species(data[2])
+			if data[2] == 'NO':
+				species = Species(data[2], molar_mass=30)
+			elif data[2] == 'NO2':
+				species = Species(data[2], molar_mass=46)
+			else:
+				species = Species(data[2])
         else:
             species = mechanism.species(data[2])
         return SpeciesProfile(pollutant, species, float(data[3]), float(data[4]), float(data[5]))
diff --git a/emips/loadApp.py b/emips/loadApp.py
index 4edf49a..ec79639 100644
--- a/emips/loadApp.py
+++ b/emips/loadApp.py
@@ -10,7 +10,7 @@ class LoadApp(PluginBase):
 
     def __init__(self):
         self.setName("EMIPS")
-        self.setAuthor("Yaqiang Wang & Wenchong Chen")
+        self.setAuthor("Yaqiang Wang & Wencong Chen")
         self.setVersion("1.0")
         self.setDescription("EMission Inventory Processing System")
         self.app_menu_item = None
diff --git a/emips/run/run_with_script/for_CMAQ/README.txt b/emips/run/run_with_script/for_CMAQ/README.txt
new file mode 100644
index 0000000..325214f
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/README.txt
@@ -0,0 +1,8 @@
+Read before running
+----------------------
+
+1, Set parameters(time, directory and dimension etc.) in "for_CMAQ.py".
+
+
+
+
diff --git a/emips/run/run_with_script/for_CMAQ/__init__.py b/emips/run/run_with_script/for_CMAQ/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/emips/run/run_with_script/for_CMAQ/for_CMAQ.py b/emips/run/run_with_script/for_CMAQ/for_CMAQ.py
new file mode 100644
index 0000000..88750c4
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/for_CMAQ.py
@@ -0,0 +1,137 @@
+"""
+# Author: Wencong Chen
+# Date: 2023-06-21
+# Purpose: Convert model-ready emission file to meet CMAQ's input requirements.
+           (height allocation, merge files of different sectors, convert pollutants and projection)
+"""
+
+import time 
+start_time = time.clock()
+#Set current working directory
+from inspect import getsourcefile
+dir_run = os.path.dirname(os.path.abspath(getsourcefile(lambda:0)))
+if not dir_run in sys.path:
+    sys.path.append(dir_run)
+
+import os
+from emips.spatial_alloc import GridDesc
+from emips import ge_data_dir
+from emips.utils import SectorEnum
+from emips.chem_spec import get_model_species_CMAQ
+from collections import OrderedDict
+import datetime
+import height
+import merge
+import transform
+import proj
+
+#Choose Chemical mechanism
+mechanism_name = 'cb05'
+print('--------CHEMICAL MECHANISM: {}--------'.format(mechanism_name.upper()))
+
+#Set year and month
+year = 2019
+month = 1
+
+#Set input&output file path
+dir_inter = os.path.join(r'G:\CMAQ', mechanism_name, r'merge\{0:}\{0:}{1:>02d}'.format(year, month))
+
+#Set original grid
+original_proj = geolib.projinfo()
+original_grid = GridDesc(original_proj, x_orig=64., x_cell=0.25, x_num=324,
+    y_orig=15., y_cell=0.25, y_num=180)
+    
+#Set target grid
+ncols = 164             #the number of x dimension
+nrows = 97              #the number of y dimension
+nlays = 7               #the number of z dimension
+xcell = 36000.0         #Grid spacing in x dimensions
+ycell = 36000.0         #Grid spacing in y dimensions
+proj_name = 'lcc'       #Name of the projection ('lcc' is Lambert Conformal Conic projection)
+xcent = 110.0           #Longitude of projection center
+ycent = 34.0            #Latitude of projection center
+p_alp = 25.0            #First standard parallel
+p_bet = 40.0            #Second standard parallel
+p_gam = 110.0
+radius = 6370000.0      #Radius of the sphere(unit: meter)
+
+if ncols % 2 == 0:
+    xorig = (ncols-2) / 2.0 * xcell + xcell / 2.0
+else:
+    xorig = (ncols-1) / 2.0 * xcell
+if nrows % 2 == 0:
+    yorig = (nrows-2) / 2.0 * ycell + ycell / 2.0
+else:
+    yorig = (nrows-1) / 2.0 * ycell
+target_proj = geolib.projinfo(proj=proj_name, lon_0=xcent, lat_0=ycent, lat_1=p_alp, lat_2=p_bet, a=radius, b=radius)
+target_grid = GridDesc(target_proj, x_orig=-xorig, x_cell=xcell, x_num=ncols,
+        y_orig=-yorig, y_cell=ycell, y_num=nrows)
+
+#Set path of the vertical allocate file
+vert_prof_file = os.path.join(ge_data_dir, 'height.txt')
+
+#Set sectors that need to be processed
+sectors = [SectorEnum.INDUSTRY, SectorEnum.AGRICULTURE, SectorEnum.ENERGY,
+        SectorEnum.RESIDENTIAL, SectorEnum.TRANSPORT] 
+
+#Acquisition of model species (gases and aerosol) and aerosol names based on chemical mechanisms
+out_species, out_species_aer = get_model_species_CMAQ(mechanism_name)
+varlist = ''
+for i in out_species:
+    varlist = varlist + '{:<16s}'.format(i)
+
+#Set global attributes
+gattrs = OrderedDict()
+#gattrs['Conventions'] = 'CF-1.6'
+gattrs['TITLE'] = "Created using MeteoInfo, mechanism: {}".format(mechanism_name.upper())
+gattrs['IOAPI_VERSION'] = "{:<80s}".format('$Id: @(#) ioapi library version 3.1 $')
+gattrs['EXEC_ID'] = "{:<80s}".format('????????????????')
+gattrs['FTYPE'] = 1
+gattrs['CDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j")  #test
+gattrs['CTIME'] = 81255                                                          #test
+gattrs['WDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j")  #test
+gattrs['WTIME'] = 81255                                                          #test
+gattrs['SDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j")
+gattrs['STIME'] = 0
+gattrs['TSTEP'] = 10000
+gattrs['NTHIK'] = 1
+gattrs['NCOLS'] = ncols
+gattrs['NROWS'] = nrows
+gattrs['NLAYS'] = nlays
+gattrs['NVARS'] = len(out_species)
+gattrs['GDTYP'] = 2
+gattrs['P_ALP'] = p_alp
+gattrs['P_BET'] = p_bet
+gattrs['P_GAM'] = p_gam
+gattrs['XCENT'] = xcent
+gattrs['YCENT'] = ycent
+gattrs['XORIG'] = -xorig
+gattrs['YORIG'] = -yorig
+gattrs['XCELL'] = xcell
+gattrs['YCELL'] = ycell
+gattrs['VGTYP'] = 2                #test
+gattrs['VGTOP'] = 10000.0
+gattrs['VGLVLS'] = 1.0, 0.995, 0.988, 0.98, 0.97, 0.956, 0.938, 0.893, 0.839, 0.777, 0.702, 0.582, 0.4, 0.2, 0.0
+gattrs['GDNAM'] = "{:<16s}".format('M_36_08CHINA')
+gattrs['UPNAM'] = ""
+gattrs['VAR-LIST'] = varlist
+gattrs['FILEDESC'] = ""
+gattrs['HISTORY'] = ""
+
+#Run all scripts
+print('Allocate according to height...')
+height.run(year, month, dir_inter, original_grid, sectors, nlays, vert_prof_file)
+print('Merge data from different sectors...')
+merge.run(year, month, dir_inter, original_grid, sectors, nlays)
+print('Distribution of particulate matter and change unit...')
+transform.run(year, month, dir_inter, original_grid, out_species, out_species_aer, nlays)
+print('Conversion projection...')
+proj.run(year, month, dir_inter, original_grid, target_grid, out_species, out_species_aer, gattrs, mechanism_name, nlays)
+print('-----------------------------')
+print('--------All finished!--------')
+print('-----------------------------')
+
+#Calculate running time
+end_time = time.clock()
+time = (end_time - start_time) / 60.0
+print('Time: {:.2f}min'.format(time))
\ No newline at end of file
diff --git a/emips/run/run_with_script/for_CMAQ/height.py b/emips/run/run_with_script/for_CMAQ/height.py
new file mode 100644
index 0000000..1331145
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/height.py
@@ -0,0 +1,93 @@
+from emips.spatial_alloc import GridDesc
+from collections import OrderedDict
+from mipylib import dataset
+import mipylib.numeric as np
+import os
+from emips.utils import emis_util
+import emips
+
+out_species_unit = ['PEC', 'POA', 'PMFINE', 'PNO3', 'PSO4', 'PMC']
+
+
+def run(year, month, dir_inter, model_grid, sectors, nlays, z_file):
+    """
+    Allocate data to different heights.
+
+    :param year: (*int*) Year.
+    :param month: (*int*) Month.
+    :param dir_inter: (*string*) The directory where data is stored.
+    :param model_grid: (*GridDesc*) Model data grid describe.
+    :param sectors: (*GridDesc*) The sectors need to be processed.
+    :param nlays: (*int*) The z-dimension of the output data.
+    :param z_file: (*string*) The path of the vertical allocate file.
+    """
+    print('Define dimension and global attributes...')
+    tdim = np.dimension(np.arange(24), 'HOUR')
+    ydim = np.dimension(model_grid.y_coord, 'ROW', 'Y')
+    xdim = np.dimension(model_grid.x_coord, 'COL', 'X')
+    zdim = np.dimension(np.arange(nlays), 'LAY')
+    dims = [tdim, zdim, ydim, xdim]
+
+    gattrs = OrderedDict()
+    gattrs['Conventions'] = 'CF-1.6'
+    gattrs['Tools'] = 'Created using MeteoInfo'
+
+    for sector in sectors:
+        fn = dir_inter + '\emis_{}_{}_{}_hour.nc'.format(sector.name, year, month)
+        print('File input: {}'.format(fn))
+        dimvars = []
+        if os.path.exists(fn):
+            f = dataset.addfile(fn)
+            for var in f.varnames:
+                if var == 'lat' or var == 'lon':
+                    continue
+                else:
+                    dimvar = dataset.DimVariable()
+                    dimvar.name = var
+                    dimvar.dtype = np.dtype.float
+                    dimvar.dims = dims
+                    dimvar.addattr('description', "EMISSION_{}".format(var))
+                    if var in out_species_unit:
+                        dimvar.addattr('units', 'g/m2/s')
+                    else:
+                        dimvar.addattr('units', 'mole/m2/s')
+                    dimvars.append(dimvar)
+
+            out_fn = dir_inter + '\emis_{}_{}_{}_hour_height.nc'.format(sector.name, year, month)
+            print('Create output file:{}'.format(out_fn))
+            ncfile = dataset.addfile(out_fn, 'c', largefile=True)
+            ncfile.nc_define(dims, gattrs, dimvars)
+
+            data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+            dd = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+            # get vertical profiles
+            scc = emis_util.get_scc(sector)
+            vertical_pro = emips.vertical_alloc.read_file(z_file, scc)
+            if nlays != len(vertical_pro.get_ratios()):
+                print('----------Warning-----------\nThe z-dimensions in the vertical profile file and the output file do not match.\n----------Warning-----------')
+            # read, merge and output
+            if round(vertical_pro.get_ratios()[0], 2) != 1.0:
+                print('Allocating: {}'.format(sector.name))
+            else:
+                print('Do not need to be allocated: {}'.format(sector.name))
+            print('Write data to file...')
+            for var in f.varnames:
+                if var == 'lat' or var == 'lon':
+                    continue
+                else:
+                    print(var)
+                    dd[:, 0, :, :] = f[var][:]
+                    if round(vertical_pro.get_ratios()[0], 2) == 1.0:
+                        data[:, 0, :, :] = dd[:, 0, :, :]
+                    else:
+                        for lay in np.arange(len(vertical_pro.get_ratios())):
+                            data[:, lay, :, :] = dd[:, 0, :, :] * vertical_pro.get_ratios()[lay]
+                    # Turn nan to zero
+                    data[data == np.nan] = 0
+                ncfile.write(var, data)
+            ncfile.close()
+            f.close()
+        else:
+            print('File not exist: {}'.format(fn))
+            continue
+    print('Allocate of height finished!')
diff --git a/emips/run/run_with_script/for_CMAQ/merge.py b/emips/run/run_with_script/for_CMAQ/merge.py
new file mode 100644
index 0000000..19cead3
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/merge.py
@@ -0,0 +1,92 @@
+from emips.spatial_alloc import GridDesc
+from collections import OrderedDict
+from mipylib import dataset
+import mipylib.numeric as np
+import os
+
+out_species_aer = ['PEC', 'POA', 'PMFINE', 'PNO3', 'PSO4', 'PMC']
+
+def run(year, month, dir_inter, model_grid, sectors, nlays):
+    """
+    Combine all sectors into one file
+
+    :param year: (*int*) Year.
+    :param month: (*int*) Month.
+    :param dir_inter: (*string*) The directory where data is stored.
+    :param model_grid: (*GridDesc*) Model data grid describe.
+    :param sectors: (*list*) Sectors that needs to be merged.
+    :param nlays: (*int*) The z-dimension of the output data.
+    """
+    print('Define dimension and global attributes...')
+    tdim = np.dimension(np.arange(24), 'HOUR')
+    zdim = np.dimension(np.arange(nlays), 'LAY')
+    ydim = np.dimension(model_grid.y_coord, 'ROW', 'Y')
+    xdim = np.dimension(model_grid.x_coord, 'COL', 'X')
+    dims = [tdim, zdim, ydim, xdim]
+    #set the definition of the output variable
+    print('Define variables...')
+    dimvars = []
+    count = []
+    dict_spec = {}
+    print('Add files...')
+    for sector in sectors:
+        fn = dir_inter + '\emis_{}_{}_{}_hour_height.nc'.format(sector.name, year, month)
+        if os.path.exists(fn):
+            print(fn)
+            f = dataset.addfile(fn)
+            for var in f.variables:
+                if var.ndim == 4:
+                    if dict_spec.has_key(var.name):
+                        dict_spec[var.name].append(fn)
+                    else:
+                        dict_spec[var.name] = [fn]
+            for var in f.varnames:
+                if var == 'ROW' or var == 'COL':
+                        continue
+                if var in count:
+                    continue
+                else:
+                    dimvar = dataset.DimVariable()
+                    dimvar.name = var
+                    dimvar.dtype = np.dtype.float
+                    dimvar.dims = dims
+                    dimvar.addattr('description', "EMISSION_{}".format(var))
+                    if var in out_species_aer:
+                        dimvar.addattr('units', 'g/m2/s')
+                    else:
+                        dimvar.addattr('units', 'mole/m2/s')
+                count.append(var)
+                dimvars.append(dimvar)
+            f.close()
+        else:
+            print('File not exist: {}'.format(fn))
+            continue
+                 
+    #set dimension and define ncfile 
+    out_fn = dir_inter + '\emis_{}_{}_hour.nc'.format(year, month)
+    gattrs = OrderedDict()
+    gattrs['Conventions'] = 'CF-1.6'
+    gattrs['Tools'] = 'Created using MeteoInfo'
+    print('Create output data file:{}'.format(out_fn))
+    ncfile = dataset.addfile(out_fn, 'c', largefile=True)
+    ncfile.nc_define(dims, gattrs, dimvars)
+    #read, merge and output
+    print('Write variables data...')
+    for sname, fns in dict_spec.iteritems():
+        print(sname)
+        spec_data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+        dd = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+        for fn in fns:
+            f = dataset.addfile(fn)
+            dd = f[sname][:]
+            #turn nan to zero
+            dd[dd==np.nan] = 0.0
+            if spec_data.sum() == 0.0:
+                spec_data = dd
+            else:
+                spec_data = spec_data + dd
+            f.close()
+        ncfile.write(sname, spec_data)
+    ncfile.close()
+    print('Merge data finished!')
+    
\ No newline at end of file
diff --git a/emips/run/run_with_script/for_CMAQ/proj.py b/emips/run/run_with_script/for_CMAQ/proj.py
new file mode 100644
index 0000000..db48243
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/proj.py
@@ -0,0 +1,117 @@
+from emips.spatial_alloc import GridDesc, transform
+from collections import OrderedDict
+from mipylib import dataset
+import mipylib.numeric as np
+import datetime
+
+def run(year, month, dir_inter, model_grid, target_grid, out_species, out_species_aer, global_attributes, mechanism_name, nlays):
+    """
+    Write Times variable, add global attributes, convert data's projection.
+
+    :param year: (*int*) Year.
+    :param month: (*int*) Month.
+    :param dir_inter: (*string*) The directory where data is stored.
+    :param model_grid: (*GridDesc*) Model data grid describe.
+    :param target_grid: (*GridDesc*) Target data grid describe.
+    :param out_species: (*list*) The name of the output species (gases and aerosol).
+    :param out_species_aer: (*list*) The name of the output species (aerosol).
+    :param global_attributes: (*OrderedDict*) The global attributes of the output file.
+    :param mechanism_name: (*string*) The name of the chemical mechanism.
+    :param nlays: (*int*) The z-dimension of the output data.
+    """
+    print('Add input file...')
+    fn_in = dir_inter + '\emis_{}_{}_hour_transform.nc'.format(year, month)
+    print(fn_in)
+    f_in = dataset.addfile(fn_in)
+    
+    #set dimension 
+    tdim = np.dimension(np.arange(25), 'TSTEP')
+    ddim = np.dimension(np.arange(2), 'DATE-TIME')
+    ydim = np.dimension(target_grid.y_coord, 'ROW', 'Y')
+    xdim = np.dimension(target_grid.x_coord, 'COL', 'X')
+    zdim = np.dimension(np.arange(nlays), 'LAY')
+    vdim = np.dimension(np.arange(len(out_species)), 'VAR')
+    vardims = [tdim, zdim, ydim, xdim]
+    all_dims = [tdim, ddim, zdim, vdim, ydim, xdim]
+    
+    #Set variables
+    dimvars = []
+    dimvar = dataset.DimVariable()
+    dimvar.name = 'TFLAG'
+    dimvar.dtype = np.dtype.int
+    dimvar.dims = [tdim, vdim, ddim]
+    dimvar.addattr('units', '<YYYYDDD,HHMMSS>')
+    dimvar.addattr('long_name', '{:<16s}'.format('TFLAG'))
+    dimvar.addattr('var_desc', '{:<80s}'.format('Timestep-valid flags:  (1) YYYYDDD or (2) HHMMSS'))
+    dimvars.append(dimvar)
+    for out_specie in out_species:
+        dimvar = dataset.DimVariable()
+        dimvar.name = out_specie
+        dimvar.dtype = np.dtype.float
+        dimvar.dims = vardims
+        dimvar.addattr('long_name', "{:<16s}".format(out_specie))
+        if out_specie in out_species_aer:
+            dimvar.addattr('units', "{:<16s}".format('g/s'))
+        else:
+            dimvar.addattr('units', "{:<16s}".format('moles/s'))
+        dimvar.addattr('var_desc', "{:<80s}".format("Hourly emission"))
+        dimvars.append(dimvar)
+
+    #Output file
+    fn_out = dir_inter + '\GR_EMIS_{}_{}{:0>2d}'.format(mechanism_name, year, month)
+    print('Create output data file...')
+    print(fn_out)
+    ncfile = dataset.addfile(fn_out, 'c', largefile=True)
+    print('Define dimensions, global attributes and variables...')
+    ncfile.nc_define(all_dims, global_attributes, dimvars, write_dimvars=False)
+    
+    #Times
+    print('Write Times variable...')
+    daytime = datetime.date(year=year, month=month, day=1)
+    daytime1 = daytime.strftime("%Y%j")
+    daytime2 = (daytime + datetime.timedelta(days=1)).strftime("%Y%j")
+    t_out = np.zeros((tdim.length, vdim.length, ddim.length))
+    for k in range(0, ddim.length):
+        for i in range(0, tdim.length):
+            if k == 0:
+                if i != tdim.length - 1:
+                    for j in range(0, vdim.length):
+                        t_out[i, j, k] = int(daytime1)  
+                else:
+                    for j in range(0, vdim.length):
+                        t_out[i, j, k] = int(daytime2)
+            else:
+                if i != tdim.length - 1:
+                    for j in range(0, vdim.length):
+                        t_out[i, j, k] = int(i * 10000)
+                else:
+                    for j in range(0, vdim.length):
+                        t_out[i, j, k] = int(0)
+    ncfile.write('TFLAG', t_out)
+    
+    #
+    print('Write variable data except times...')
+    for out_specie in out_species:
+        data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+        dd = np.zeros((tdim.length, zdim.length, model_grid.y_num, model_grid.x_num))
+        if out_specie in f_in.varnames:
+            print(out_specie)
+            for i in range(tdim.length):
+                j = i
+                if i == tdim.length - 1:
+                    j = 0
+                dd[i, :, :, :] = f_in[out_specie][j, :, :, :]
+            #Conversion
+            dd = transform(dd, model_grid, target_grid)
+            #Set the fourth dimension 
+            #dd = dd.reshape(12, 1, ydim.length, xdim.length)
+            #Set default values
+            dd[dd==np.nan] = 0
+            data[:, :, :, :] = dd
+        else:
+            print('{} no data!'.format(out_specie))
+        ncfile.write(out_specie, data)
+    ncfile.close()
+    f_in.close()
+    print('Convert projection finished!')
+    
\ No newline at end of file
diff --git a/emips/run/run_with_script/for_CMAQ/transform.py b/emips/run/run_with_script/for_CMAQ/transform.py
new file mode 100644
index 0000000..61fbb5c
--- /dev/null
+++ b/emips/run/run_with_script/for_CMAQ/transform.py
@@ -0,0 +1,90 @@
+from emips.spatial_alloc import GridDesc
+from collections import OrderedDict
+from mipylib import dataset
+import mipylib.numeric as np
+
+def run(year, month, dir_inter, model_grid, out_species, out_species_aer, nlays):
+    """
+    Change unit and match variable name.
+
+    :param year: (*int*) Year.
+    :param month: (*int*) Month.
+    :param dir_inter: (*string*) The directory where data is stored.
+    :param model_grid: (*GridDesc*) Model data grid describe.
+    :param out_species: (*list*) The name of the output species (gases and aerosol).
+    :param out_species_aer: (*list*) The name of the output species (aerosol).
+    :param nlays: (*int*) The z-dimension of the output data.
+    """
+    print('Add input file...')
+    fn_in = dir_inter + '\emis_{}_{}_hour.nc'.format(year, month)
+    f_in = dataset.addfile(fn_in)
+    
+    #Set dimension
+    print('Define dimensions and global attributes...')
+    tdim = np.dimension(np.arange(24), 'TSTEP')
+    zdim = np.dimension(np.arange(nlays), 'LAY')
+    ydim = np.dimension(model_grid.y_coord, 'ROW', 'Y')
+    xdim = np.dimension(model_grid.x_coord, 'COL', 'X')
+    dims = [tdim, zdim, ydim, xdim]
+    gattrs = OrderedDict()
+    gattrs['Conventions'] = 'CF-1.6'
+    gattrs['Tools'] = 'Created using MeteoInfo'
+    
+    #Set the definition of the output variable and ncfile
+    fn_out = dir_inter + '\emis_{}_{}_hour_transform.nc'.format(year, month)
+    print('Define variables and output file...')
+    dimvars = []
+    for out_specie in out_species:
+        dimvar = dataset.DimVariable()
+        dimvar.name = out_specie
+        dimvar.dtype = np.dtype.float
+        dimvar.dims = dims
+        dimvar.addattr('long_name', "{:<16s}".format(out_specie))
+        if out_specie in out_species_aer:
+            #g/m2/s to g/grid/s
+            dimvar.addattr('units', "{:<16s}".format('g/s'))
+        else:
+            #mole/m2/s to moles/grid/s
+            dimvar.addattr('units', "{:<16s}".format('moles/s'))
+        dimvar.addattr('var_desc', "{:<80s}".format("Hourly emission"))
+        dimvars.append(dimvar)
+    print('Create output data file:{}'.format(fn_out))
+    ncfile = dataset.addfile(fn_out, 'c', largefile=True)
+    ncfile.nc_define(dims, gattrs, dimvars)    
+
+    #calculate grid areas (meter)
+    grid_areas = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+    for i in range(tdim.length):
+        for j in range(zdim.length):
+            grid_areas[i, j, :, :] = model_grid.grid_areas()
+    
+    #add data to ncfile
+    print('Process data and write to file...')
+    for name in out_species:
+        data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+        sname = name
+        print(sname)
+        if sname in f_in.varnames:
+            data = f_in[sname][:, :, :, :]
+#            data = data * 1e6
+            data = data * grid_areas
+            ncfile.write(name, data)
+        elif sname == 'UNR':
+            data = f_in['NR'][:, :, :]
+            data = data * grid_areas
+            ncfile.write(name, data)
+        elif sname == 'PMOTHR':
+            data = f_in['PMFINE'][:, :, :]
+            data = data * grid_areas
+            ncfile.write(name, data)
+        #####test#####
+        elif sname == 'POC':
+            data = f_in['POA'][:, :, :]
+            data = data * grid_areas
+            ncfile.write(name, data)
+        #####test#####
+        else:
+            ncfile.write(name, data)
+    f_in.close()
+    ncfile.close()      
+    print('Distribution of particulate matter and change unit finised!')
diff --git a/emips/run/run_with_script/for_CUACE/convert_grads$py.class b/emips/run/run_with_script/for_CUACE/convert_grads$py.class
deleted file mode 100644
index b2244c6..0000000
Binary files a/emips/run/run_with_script/for_CUACE/convert_grads$py.class and /dev/null differ
diff --git a/emips/run/run_with_script/for_CUACE/for_CUACE$py.class b/emips/run/run_with_script/for_CUACE/for_CUACE$py.class
deleted file mode 100644
index 2c1354f..0000000
Binary files a/emips/run/run_with_script/for_CUACE/for_CUACE$py.class and /dev/null differ
diff --git a/emips/run/run_with_script/for_CUACE/for_CUACE.py b/emips/run/run_with_script/for_CUACE/for_CUACE.py
index 76ac600..7c9954e 100644
--- a/emips/run/run_with_script/for_CUACE/for_CUACE.py
+++ b/emips/run/run_with_script/for_CUACE/for_CUACE.py
@@ -4,40 +4,39 @@
 # Purpose: Convert netcdf model-ready emission file to GrADS data format for CUACE model 
            and write the description file of the output binary data files.
 """
-import convert_grads
-import write_ctl
 import time
+time_start = time.time()
+#Set current working directory
+from inspect import getsourcefile
+dir_run = os.path.dirname(os.path.abspath(getsourcefile(lambda:0)))
+if not dir_run in sys.path:
+    sys.path.append(dir_run)
 
-def run(run_config):
-    """
-    To CUACE model ready emission data files.
-
-    :param run_config: (*RunConfigure*) The run configure.
-    """
-    time_start = time.time()
-
-    year = run_config.emission_year
-    months = [run_config.emission_month]
-    dir_in = run_config.run_output_dir
-    dir_out = run_config.run_output_dir
-    model_grid = run_config.spatial_model_grid
-    xn = model_grid.x_num
-    yn = model_grid.y_num
-    xmin = model_grid.x_orig
-    ymin = model_grid.y_orig
-    xdelta = model_grid.x_cell
-    ydelta = model_grid.y_cell
+#set year and months
+year = 2017
+months = [1]
+##Set directory for input and output files
+dir_in = r'G:\test'
+dir_out = r'G:\test'
+#Set dimension length of output file.
+xn = 324
+yn = 180
+xmin = 64.0
+ymin = 15.0
+xdelta = 0.25
+ydelta = 0.25
 
-    print('Convert to grads...')
-    convert_grads.run(year, months, dir_in, dir_out, xn, yn)
-
-    print('Write .ctl files...')
-    write_ctl.run(year, months, dir_out, xn, yn, xmin, ymin, xdelta, ydelta)
+print('Convert to grads...')
+import convert_grads
+convert_grads.run(year, months, dir_in, dir_out, xn, yn)
 
-    print('-------------------')
-    print('---All finished!---')
-    print('-------------------')
+print('Write .ctl files...')
+import write_ctl
+write_ctl.run(year, months, dir_out, xn, yn, xmin, ymin, xdelta, ydelta)
 
-    time_end = time.time()
-    time_count = (time_end - time_start) / 60
-    print('Time: {:.2f}min'.format(time_count))
+print('-------------------')
+print('---All finished!---')
+print('-------------------')
+time_end = time.time()
+time = (time_end - time_start) / 60
+print('Time: {:.2f}min'.format(time))
diff --git a/emips/run/run_with_script/for_CUACE/for_CUACE_bak.py b/emips/run/run_with_script/for_CUACE/for_CUACE_bak.py
deleted file mode 100644
index 7c9954e..0000000
--- a/emips/run/run_with_script/for_CUACE/for_CUACE_bak.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-# Author: Wencong Chen
-# Date: 2022-11-27
-# Purpose: Convert netcdf model-ready emission file to GrADS data format for CUACE model 
-           and write the description file of the output binary data files.
-"""
-import time
-time_start = time.time()
-#Set current working directory
-from inspect import getsourcefile
-dir_run = os.path.dirname(os.path.abspath(getsourcefile(lambda:0)))
-if not dir_run in sys.path:
-    sys.path.append(dir_run)
-
-#set year and months
-year = 2017
-months = [1]
-##Set directory for input and output files
-dir_in = r'G:\test'
-dir_out = r'G:\test'
-#Set dimension length of output file.
-xn = 324
-yn = 180
-xmin = 64.0
-ymin = 15.0
-xdelta = 0.25
-ydelta = 0.25
-
-print('Convert to grads...')
-import convert_grads
-convert_grads.run(year, months, dir_in, dir_out, xn, yn)
-
-print('Write .ctl files...')
-import write_ctl
-write_ctl.run(year, months, dir_out, xn, yn, xmin, ymin, xdelta, ydelta)
-
-print('-------------------')
-print('---All finished!---')
-print('-------------------')
-time_end = time.time()
-time = (time_end - time_start) / 60
-print('Time: {:.2f}min'.format(time))
diff --git a/emips/run/run_with_script/for_CUACE/write_ctl$py.class b/emips/run/run_with_script/for_CUACE/write_ctl$py.class
deleted file mode 100644
index 5e37e0d..0000000
Binary files a/emips/run/run_with_script/for_CUACE/write_ctl$py.class and /dev/null differ
diff --git a/emips/run/run_with_script/for_WRFChem/for_WRFChem.py b/emips/run/run_with_script/for_WRFChem/for_WRFChem.py
index 516f9a2..57fe678 100644
--- a/emips/run/run_with_script/for_WRFChem/for_WRFChem.py
+++ b/emips/run/run_with_script/for_WRFChem/for_WRFChem.py
@@ -17,7 +17,7 @@
 from emips.spatial_alloc import GridDesc
 from emips import ge_data_dir
 from emips.utils import SectorEnum
-from emips.chem_spec import get_model_species_wrf
+from emips.chem_spec import get_model_species_WRFChem
 from collections import OrderedDict
 import height
 import merge
@@ -26,29 +26,52 @@
 #Generate one(proj.py) or two(proj_2_file.py) emission files
 import proj_2_file as proj
 
+#Choose Chemical mechanism
 mechanism_name = 'radm2'
 print('--------CHEMICAL MECHANISM: {}--------'.format(mechanism_name.upper()))
-#set year and month
+
+#Set year and month
 year = 2017
 month = 1
-#set file path
+
+#Set input&output file path
 dir_inter = os.path.join(r'G:\test', mechanism_name, r'merge\{0:}\{0:}{1:>02d}'.format(year, month))
 
-#set  model grid
-model_proj = geolib.projinfo()
-model_grid = GridDesc(model_proj, x_orig=64., x_cell=0.25, x_num=324,
+#Set original grid
+original_proj = geolib.projinfo()
+original_grid = GridDesc(original_proj, x_orig=64., x_cell=0.25, x_num=324,
     y_orig=15., y_cell=0.25, y_num=180)
+    
 #set target grid
-target_proj = geolib.projinfo(proj='lcc', lon_0=103.5, lat_0=36.500008, lat_1=30.0, lat_2=60.0, a=6370000, b=6370000)
-target_grid = GridDesc(target_proj, x_orig=-2497499.597352108, x_cell=15000.0, x_num=334,
-        y_orig=-2047499.8096037393, y_cell=15000.0, y_num=274)
+e_we = 335              #the number of x dimension
+e_sn = 275              #the number of y dimension
+zdim = 7                #the number of z dimension
+dx = 15000.0            #Grid spacing in x dimensions
+dy = 15000.0            #Grid spacing in y dimensions
+proj_name = 'lcc'       #Name of the projection ('lcc' is Lambert Conformal Conic projection)
+ref_lat = 36.500008     #Latitude of projection center
+ref_lon = 103.5         #Longitude of projection center
+truelat1 = 30.0         #First standard parallel
+truelat2 = 60.0         #Second standard parallel
+stand_lon = 103.5
+radius = 6370000.0      #Radius of the sphere(unit: meter)
+
+if (e_we-1) % 2 == 0:
+    xmin = (e_we-3) / 2.0 * dx + dx / 2.0
+else:
+    xmin = (e_we-2) / 2.0 * dx
+if (e_sn-1) % 2 == 0:
+    ymin = (e_sn-3) / 2.0 * dy + dy / 2.0
+else:
+    ymin = (e_sn-2) / 2.0 * dy 
+target_proj = geolib.projinfo(proj=proj_name, lon_0=ref_lon, lat_0=ref_lat, lat_1=truelat1, lat_2=truelat2, a=radius, b=radius)
+target_grid = GridDesc(target_proj, x_orig=-xmin, x_cell=dx, x_num=e_we-1,
+        y_orig=-ymin, y_cell=dy, y_num=e_sn-1)
 
-#set the number of z dimension
-zdim = 7
-#set path of the vertical allocate file
-z_file = os.path.join(ge_data_dir, 'height.txt')    
+#Set path of the vertical allocate file
+vert_prof_file = os.path.join(ge_data_dir, 'height.txt')
 
-#set sectors that need to be processed
+#Set sectors that need to be processed
 sectors = [SectorEnum.INDUSTRY, SectorEnum.AGRICULTURE, SectorEnum.ENERGY,
         SectorEnum.RESIDENTIAL, SectorEnum.TRANSPORT] 
 
@@ -57,29 +80,29 @@
 #gattrs['Conventions'] = 'CF-1.6'
 gattrs['TITLE'] = 'Created using MeteoInfo, mechanism: {}'.format(mechanism_name.upper())
 gattrs['START_DATE'] = "{}-{:0>2d}-01_00:00:00".format(year, month)
-gattrs['WEST-EAST_GRID_DIMENSION'] = 335
-gattrs['SOUTH-NORTH_GRID_DIMENSION'] = 275
-gattrs['BOTTOM-TOP_GRID_DIMENSION'] = 8
-gattrs['DX'] = 15000.0
-gattrs['DY'] = 15000.0
-gattrs['WEST-EAST_PATCH_START_UNSTAG '] = 1
-gattrs['WEST-EAST_PATCH_END_UNSTAG'] = 334
+gattrs['WEST-EAST_GRID_DIMENSION'] = e_we
+gattrs['SOUTH-NORTH_GRID_DIMENSION'] = e_sn
+gattrs['BOTTOM-TOP_GRID_DIMENSION'] = zdim + 1
+gattrs['DX'] = dx
+gattrs['DY'] = dy
+gattrs['WEST-EAST_PATCH_START_UNSTAG'] = 1
+gattrs['WEST-EAST_PATCH_END_UNSTAG'] = e_we - 1
 gattrs['WEST-EAST_PATCH_START_STAG'] = 1
-gattrs['WEST-EAST_PATCH_END_STAG'] = 335
+gattrs['WEST-EAST_PATCH_END_STAG'] = e_we
 gattrs['SOUTH-NORTH_PATCH_START_UNSTAG'] = 1
-gattrs['SOUTH-NORTH_PATCH_END_UNSTAG'] = 274
+gattrs['SOUTH-NORTH_PATCH_END_UNSTAG'] = e_sn - 1
 gattrs['SOUTH-NORTH_PATCH_START_STAG'] = 1
-gattrs['SOUTH-NORTH_PATCH_END_STAG'] = 275
+gattrs['SOUTH-NORTH_PATCH_END_STAG'] = e_sn
 gattrs['BOTTOM-TOP_PATCH_START_UNSTAG'] = 1
-gattrs['BOTTOM-TOP_PATCH_END_UNSTAG'] = 7
+gattrs['BOTTOM-TOP_PATCH_END_UNSTAG'] = zdim
 gattrs['BOTTOM-TOP_PATCH_START_STAG'] = 1
-gattrs['BOTTOM-TOP_PATCH_END_STAG'] = 8
-gattrs['CEN_LAT'] = 36.500008
-gattrs['CEN_LON'] = 103.5
-gattrs['TRUELAT1'] = 30.0
-gattrs['TRUELAT2'] = 60.0
-gattrs['MOAD_CEN_LAT'] = 36.500008
-gattrs['STAND_LON'] = 103.5
+gattrs['BOTTOM-TOP_PATCH_END_STAG'] = zdim + 1
+gattrs['CEN_LAT'] = ref_lat
+gattrs['CEN_LON'] = ref_lon
+gattrs['TRUELAT1'] = truelat1
+gattrs['TRUELAT2'] = truelat2
+gattrs['MOAD_CEN_LAT'] = ref_lat
+gattrs['STAND_LON'] = stand_lon
 gattrs['POLE_LAT'] = 90.0
 gattrs['POLE_LON'] = 0.0
 gattrs['GMT'] = 0.0
@@ -93,25 +116,24 @@
 gattrs['ISURBAN'] = 13
 gattrs['ISOILWATER'] = 14
 
-######set out species(gases and aerosol) and out species(aerosol)
-#################################################################################################
+##Acquisition of model species (gases and aerosol) and aerosol names based on chemical mechanisms
 ###------CB05,emiss_opt=15; RADM2,emiss_opt=3; SAPRC99,emiss_opt=13; MOZART,emiss_opt=10------###
-#################################################################################################
-out_species, out_species_aer = get_model_species_wrf(mechanism_name)
+out_species, out_species_aer = get_model_species_WRFChem(mechanism_name)
 
-#run all scripts
+#Run all scripts
 print('Allocate according to height...')
-height.run(year, month, dir_inter, model_grid, sectors, zdim, z_file)
+height.run(year, month, dir_inter, original_grid, sectors, zdim, vert_prof_file)
 print('Merge data from different sectors...')
-merge.run(year, month, dir_inter, model_grid, sectors, zdim)
+merge.run(year, month, dir_inter, original_grid, sectors, zdim)
 print('Distribution of particulate matter and change unit...')
-transform.run(year, month, dir_inter, model_grid, out_species, out_species_aer, zdim)
+transform.run(year, month, dir_inter, original_grid, out_species, out_species_aer, zdim)
 print('Conversion projection...')
-proj.run(year, month, dir_inter, model_grid, target_grid, out_species, out_species_aer, gattrs, mechanism_name, zdim)
+proj.run(year, month, dir_inter, original_grid, target_grid, out_species, out_species_aer, gattrs, mechanism_name, zdim)
 print('-----------------------------')
 print('--------All finished!--------')
 print('-----------------------------')
 
+#Calculate running time
 end_time = time.clock()
-time = (end_time - start_time)/60
+time = (end_time - start_time) / 60.0
 print('Time: {:.2f}min'.format(time))
\ No newline at end of file
diff --git a/emips/run/run_with_script/for_WRFChem/height.py b/emips/run/run_with_script/for_WRFChem/height.py
index 06fe943..80c017c 100644
--- a/emips/run/run_with_script/for_WRFChem/height.py
+++ b/emips/run/run_with_script/for_WRFChem/height.py
@@ -54,15 +54,17 @@ def run(year, month, dir_inter, model_grid, sectors, z, z_file):
                     dimvars.append(dimvar)
 
             out_fn = dir_inter + '\emis_{}_{}_{}_hour_height.nc'.format(sector.name, year, month)
-            print('Create output data file:{}'.format(out_fn))
+            print('Create output file:{}'.format(out_fn))
             ncfile = dataset.addfile(out_fn, 'c', largefile=True)
             ncfile.nc_define(dims, gattrs, dimvars)
 
-            data = np.zeros((tdim.length, z, ydim.length, xdim.length))
-            dd = np.zeros((tdim.length, z, ydim.length, xdim.length))
+            data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+            dd = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
             # get vertical profiles
             scc = emis_util.get_scc(sector)
             vertical_pro = emips.vertical_alloc.read_file(z_file, scc)
+            if z != len(vertical_pro.get_ratios()):
+                print('----------Warning-----------\nThe z-dimensions in the vertical profile file and the output file do not match.\n----------Warning-----------')
             # read, merge and output
             if round(vertical_pro.get_ratios()[0], 2) != 1.0:
                 print('Allocating: {}'.format(sector.name))
diff --git a/emips/run/run_with_script/for_WRFChem/merge.py b/emips/run/run_with_script/for_WRFChem/merge.py
index 9e2cb6a..4e48669 100644
--- a/emips/run/run_with_script/for_WRFChem/merge.py
+++ b/emips/run/run_with_script/for_WRFChem/merge.py
@@ -74,8 +74,8 @@ def run(year, month, dir_inter, model_grid, sectors, z):
     print('Write variables data...')
     for sname, fns in dict_spec.iteritems():
         print(sname)
-        spec_data = np.zeros((tdim.length, z, ydim.length, xdim.length))
-        dd = np.zeros((tdim.length, z, ydim.length, xdim.length))
+        spec_data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
+        dd = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
         for fn in fns:
             f = dataset.addfile(fn)
             dd = f[sname][:]
diff --git a/emips/run/run_with_script/for_WRFChem/proj_2_file.py b/emips/run/run_with_script/for_WRFChem/proj_2_file.py
index babbea4..e478e68 100644
--- a/emips/run/run_with_script/for_WRFChem/proj_2_file.py
+++ b/emips/run/run_with_script/for_WRFChem/proj_2_file.py
@@ -90,6 +90,7 @@ def run(year, month, dir_inter, model_grid, target_grid, out_species, out_specie
                 #dd = dd.reshape(12, 1, ydim.length, xdim.length)
                 #Set default values
                 dd[dd==np.nan] = 0
+                data[:, :, :, :] = dd
             else:
                 print('{} no data!'.format(out_specie))
                 '''
@@ -101,7 +102,6 @@ def run(year, month, dir_inter, model_grid, target_grid, out_species, out_specie
                 #Set default values
                 dd[dd==np.nan] = 0
                 '''
-            data[:, :, :, :] = dd
             ncfile.write(out_specie, data)
         ncfile.close()
     f_in.close()
diff --git a/emips/run/run_with_script/for_WRFChem/transform.py b/emips/run/run_with_script/for_WRFChem/transform.py
index 9f4db04..cc58263 100644
--- a/emips/run/run_with_script/for_WRFChem/transform.py
+++ b/emips/run/run_with_script/for_WRFChem/transform.py
@@ -59,7 +59,7 @@ def run(year, month, dir_inter, model_grid, out_species, out_species_aer, z):
     #add data to ncfile
     print('Process data and write to file...')
     for name in out_species:
-        data = np.zeros((tdim.length, z, ydim.length, xdim.length))
+        data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length))
         sname = name[2:]
         print(sname)
         if sname in f_in.varnames:
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/emission$py.class b/emips/run/run_with_script/run_meic_cams_htap/meic/emission$py.class
deleted file mode 100644
index b13e055..0000000
Binary files a/emips/run/run_with_script/run_meic_cams_htap/meic/emission$py.class and /dev/null differ
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/emission.py b/emips/run/run_with_script/run_meic_cams_htap/meic/emission.py
deleted file mode 100644
index 347d25f..0000000
--- a/emips/run/run_with_script/run_meic_cams_htap/meic/emission.py
+++ /dev/null
@@ -1,56 +0,0 @@
-from emips.utils import EmissionReader, SectorEnum
-from emips.chem_spec import PollutantEnum
-from emips.spatial_alloc import GridDesc
-import os
-from mipylib import geolib
-from mipylib.dataset import addfile
-
-
-class MyEmissionReader(EmissionReader):
-
-    def get_emis_fn(self, sector, pollutant, year, month):
-        sector_name = sector.name.lower()
-        if sector == SectorEnum.ENERGY:
-            sector_name = 'power'
-        elif sector == SectorEnum.TRANSPORT:
-            sector_name = 'transportation'
-        pollutant_name = pollutant.name.upper()
-        if pollutant == PollutantEnum.PM2_5:
-            pollutant_name = 'PM25'
-        elif pollutant == PollutantEnum.NMVOC:
-            pollutant_name = 'VOC'
-        fn = 'meic-2017_%i_%s-%s.nc' % (month, sector_name, pollutant_name)
-        return os.path.join(self.dir_emission, fn)
-
-    def read_emis(self, sector, pollutant, year, month):
-        fn = get_emis_fn(sector, pollutant, year, month)
-        f = addfile(fn)
-        data = f['z'][:]
-        nx = 800
-        ny = 500
-        data = data.reshape(ny, nx)
-        data = data[::-1,:]
-        return data
-
-    def get_emis_grid(self):
-        return GridDesc(geolib.projinfo(), x_orig=70.05, x_cell=0.1, x_num=800,
-                        y_orig=10.05, y_cell=0.1, y_num=500)
-
-
-_emis_reader = MyEmissionReader(dir_emission=r'D:\KeyData\Emission\MEIC\2017\nc')
-
-
-def get_emis_fn(sector, pollutant, year, month):
-    return _emis_reader.get_emis_fn(sector, pollutant, year, month)
-
-
-def read_emis(sector, pollutant, year, month):
-    print("sector: {}".format(sector))
-    return _emis_reader.read_emis(sector, pollutant, year, month)
-
-
-def get_emis_grid():
-    return _emis_reader.get_emis_grid()
-
-
-grid_areas = _emis_reader.get_emis_grid().grid_areas()
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/emission_meic_2017.py b/emips/run/run_with_script/run_meic_cams_htap/meic/emission_meic_2017.py
index 4ca3fa8..442bb38 100644
--- a/emips/run/run_with_script/run_meic_cams_htap/meic/emission_meic_2017.py
+++ b/emips/run/run_with_script/run_meic_cams_htap/meic/emission_meic_2017.py
@@ -40,7 +40,7 @@ def get_emis_fn(sector, pollutant, year, month):
         pollutant_name = 'PM25'
     elif pollutant == PollutantEnum.NMVOC:
         pollutant_name = 'VOC'
-    fn = '2017_{:0>2d}_{}_{}.asc'.format(month, sector_name, pollutant_name)
+    fn = '{}_{:0>2d}_{}_{}.asc'.format(year, month, sector_name, pollutant_name)
     return os.path.join(dir_emission, fn)
 
 
@@ -54,7 +54,7 @@ def read_emis(sector, pollutant, year, month):
     :param month: (*int*) The month.
     :returns: (*array*) Emission data array.
     """
-    fn = get_emis_fn(sector, pollutant, month)
+    fn = get_emis_fn(sector, pollutant, year, month)
     print('File_in:{}'.format(fn))
     f = addfile_ascii_grid(fn)
     data = f['var'][:]
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec$py.class b/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec$py.class
deleted file mode 100644
index 4aa5eff..0000000
Binary files a/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec$py.class and /dev/null differ
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec.py b/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec.py
deleted file mode 100644
index a861c3d..0000000
--- a/emips/run/run_with_script/run_meic_cams_htap/meic/grid_spec.py
+++ /dev/null
@@ -1,80 +0,0 @@
-from emips.chem_spec import GridSpecReader
-from emips.utils import SectorEnum
-from emips.spatial_alloc import GridDesc
-import os
-from mipylib import dataset
-from mipylib import numeric as np
-
-
-def get_sector_str(sector):
-    """
-    Get sector string.
-    :param sector: (*Sector*) The sector.
-    :return: (*str*) Sector string.
-    """
-    if sector == SectorEnum.INDUSTRY:
-        return "inc"
-    elif sector == SectorEnum.AGRICULTURE:
-        return "agr"
-    elif sector == SectorEnum.ENERGY:
-        return "pow"
-    elif sector == SectorEnum.RESIDENTIAL:
-        return "res"
-    elif sector in [SectorEnum.TRANSPORT, SectorEnum.AIR, SectorEnum.SHIPS]:
-        return "tra"
-    else:
-        return None
-
-
-class RetroGridSpecReader(GridSpecReader):
-
-    def get_spec_fn(self, sector):
-        self.sector = sector
-        sector_str = get_sector_str(sector)
-        fn = os.path.join(self.dir_grid, 'retro_nmvoc_ratio_{}_2000_0.1deg.nc'.format(sector_str))
-        return fn
-
-    def get_spec_vars(self, sector, dims):
-        fn = self.get_spec_fn(sector)
-        self.f = dataset.addfile(fn)
-        spec_vars = []
-        for var in self.f.variables:
-            if var.ndim == 2:
-                spec_var = dataset.DimVariable()
-                spec_var.name = var.name
-                spec_var.dtype = np.dtype.float
-                spec_var.dims = dims
-                spec_var.addattr('units', 'g/m2/s')
-                spec_vars.append(spec_var)
-
-        return spec_vars
-
-    def read_spec(self, sector, spec_var):
-        if self.sector != sector or self.f is None:
-            fn = self.get_spec_fn(sector)
-            self.f = dataset.addfile(fn)
-
-        return self.f[spec_var.name][:]
-
-    def get_spec_grid(self):
-        return GridDesc(x_orig=0.05, x_cell=0.1, x_num=3600,
-                        y_orig=-89.95, y_cell=0.1, y_num=1800)
-
-
-_grid_spec_reader = RetroGridSpecReader(dir_grid="Z:\chen\Research\EMIPS\Grid_speciation_data(VOC)")
-
-
-def get_spec_fn(sector):
-    return _grid_spec_reader.get_spec_fn(sector)
-
-
-def get_spec_vars(sector, dims):
-    return _grid_spec_reader.get_spec_vars(sector, dims)
-
-
-def read_spec(sector, spec_var):
-    return _grid_spec_reader.read_spec(sector, spec_var)
-
-
-def get_spec_grid():
-    return _grid_spec_reader.get_spec_grid()
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/run_VOC.py b/emips/run/run_with_script/run_meic_cams_htap/meic/run_VOC.py
index 52b2dc8..24056ca 100644
--- a/emips/run/run_with_script/run_meic_cams_htap/meic/run_VOC.py
+++ b/emips/run/run_with_script/run_meic_cams_htap/meic/run_VOC.py
@@ -50,7 +50,7 @@ def run(year, month, dir_inter, emission, model_grid):
         scc = emis_util.get_scc(sector)
 
         print('Read emission data...')
-        emis_data = emission.read_emis(sector, pollutant, month)
+        emis_data = emission.read_emis(sector, pollutant, year, month)
 
         #### Spatial allocation  
         print('Spatial allocation...')
diff --git a/emips/run/run_with_script/run_meic_cams_htap/meic/run_pollutants.py b/emips/run/run_with_script/run_meic_cams_htap/meic/run_pollutants.py
index 61faa64..1675fce 100644
--- a/emips/run/run_with_script/run_meic_cams_htap/meic/run_pollutants.py
+++ b/emips/run/run_with_script/run_meic_cams_htap/meic/run_pollutants.py
@@ -61,15 +61,15 @@ def run(year, month, dir_inter, emission, model_grid):
             print(pollutant)
 
             print('Read emission data...')
-            emis_data = emission.read_emis(sector, pollutant, month)
+            emis_data = emission.read_emis(sector, pollutant, year, month)
 
             # Remove PM2.5 included in PM10, Remove BC and OC included in PM2.5
             if pollutant == PollutantEnum.PM10:
-                emis_data_pm25 = emission.read_emis(sector, PollutantEnum.PM2_5, month)
+                emis_data_pm25 = emission.read_emis(sector, PollutantEnum.PM2_5, year, month)
                 emis_data = emis_data - emis_data_pm25
             if pollutant == PollutantEnum.PM2_5:
-                emis_data_bc = emission.read_emis(sector, PollutantEnum.BC, month)
-                emis_data_oc = emission.read_emis(sector, PollutantEnum.OC, month)
+                emis_data_bc = emission.read_emis(sector, PollutantEnum.BC, year, month)
+                emis_data_oc = emission.read_emis(sector, PollutantEnum.OC, year, month)
                 emis_data = emis_data - emis_data_bc - emis_data_oc
 
             # Spatial allocation