adding guides, fixing anthro_emis.inp

blueprints/leeds/emis_edgarhtap2-meic2015_mozmos.inp

@@ -4,7 +4,7 @@
  src_file_prefix = 'EDGARHTAP2_MEIC2015_'
  src_file_suffix = '_2010.0.1x0.1.nc'
  src_names = 'CO(28)','NOx(30)','SO2(64)','NH3(17)','BC(12)',
-             'OM(12)','OIN_PM2.5(1)','PM2.5_10(1)',
+             'POM(12)','OIN_PM2.5(1)','PM2.5_10(1)',
              'BIGALK(72)','BIGENE(56)','C2H4(28)','C2H5OH(46)','C2H6(30)','CH2O(30)',
              'CH3CHO(44)','CH3COCH3(58)','CH3OH(32)','MEK(72)','TOLUENE(92)','C3H6(42)','C3H8(44)',
              'BENZENE(78)','XYLENE(106)','CH4(16)'
@@ -73,16 +73,16 @@
             'NO2_SHP->0.2*NOx(emis_shp)','NO2_CDS->0.2*NOx(emis_cds)',
             'NO2_CRS->0.2*NOx(emis_crs)','NO2_LTO->0.2*NOx(emis_lto)',
             'NO2->0.2*NOx(emis_tot)',
-            'ORGI_TRA(a)->0.1*OM(emis_tra)','ORGI_IND(a)->0.1*OM(emis_ind)',
-            'ORGI_RES(a)->0.1*OM(emis_res)','ORGI_POW(a)->0.1*OM(emis_ene)',
-            'ORGI_SHP(a)->0.1*OM(emis_shp)','ORGI_CDS(a)->0.1*OM(emis_cds)',
-            'ORGI_CRS(a)->0.1*OM(emis_crs)','ORGI_LTO(a)->0.1*OM(emis_lto)',
-            'ORGI(a)->0.1*OM(emis_tot)',
-            'ORGJ_TRA(a)->0.9*OM(emis_tra)','ORGJ_IND(a)->0.9*OM(emis_ind)',
-            'ORGJ_RES(a)->0.9*OM(emis_res)','ORGJ_POW(a)->0.9*OM(emis_ene)',
-            'ORGJ_SHP(a)->0.9*OM(emis_shp)','ORGJ_CDS(a)->0.9*OM(emis_cds)',
-            'ORGJ_CRS(a)->0.9*OM(emis_crs)','ORGJ_LTO(a)->0.9*OM(emis_lto)',
-            'ORGJ(a)->0.9*OM(emis_tot)',
+            'ORGI_TRA(a)->0.1*POM(emis_tra)','ORGI_IND(a)->0.1*POM(emis_ind)',
+            'ORGI_RES(a)->0.1*POM(emis_res)','ORGI_POW(a)->0.1*POM(emis_ene)',
+            'ORGI_SHP(a)->0.1*POM(emis_shp)','ORGI_CDS(a)->0.1*POM(emis_cds)',
+            'ORGI_CRS(a)->0.1*POM(emis_crs)','ORGI_LTO(a)->0.1*POM(emis_lto)',
+            'ORGI(a)->0.1*POM(emis_tot)',
+            'ORGJ_TRA(a)->0.9*POM(emis_tra)','ORGJ_IND(a)->0.9*POM(emis_ind)',
+            'ORGJ_RES(a)->0.9*POM(emis_res)','ORGJ_POW(a)->0.9*POM(emis_ene)',
+            'ORGJ_SHP(a)->0.9*POM(emis_shp)','ORGJ_CDS(a)->0.9*POM(emis_cds)',
+            'ORGJ_CRS(a)->0.9*POM(emis_crs)','ORGJ_LTO(a)->0.9*POM(emis_lto)',
+            'ORGJ(a)->0.9*POM(emis_tot)',
             'PM_10_TRA(a)->PM2.5_10(emis_tra)','PM_10_IND(a)->PM2.5_10(emis_ind)',
             'PM_10_RES(a)->PM2.5_10(emis_res)','PM_10_POW(a)->PM2.5_10(emis_ene)',
             'PM_10_SHP(a)->PM2.5_10(emis_shp)','PM_10_CDS(a)->PM2.5_10(emis_cds)',

guides/add_bbinjectscheme.md

+### How to add the `bbinjectscheme` options to WRFChem
+
+#### Steps
+1. `Registry/registry.chem`  
+    a. Include `bbinjectscheme` flag within `registry.chem` to enable injection scheme switching using this chemistry namelist flag. Add the following code within the chemistry scheme declaration section with the other namelist flags, previously used line 3825:  
+    ```fortran
+    rconfig   integer     bbinjectscheme      namelist,chem          max_domains    4       rh    "bbinjectscheme"      ""      ""
+    ```
+
+2. `chem/chem_driver.F`  
+    a. Include `grid%kpbl` within the `emissions_driver` call to indicate the model level containing the top of boundary layer. Add the following code to line 879, underneath `grid%dust_flux, grid%seas_flux,`:  
+    ```fortran
+             ! bbinjectscheme requires record of layer which PBL goes up to
+                  grid%kpbl,                                                                   & 
+    ```
+
+3. `chem/emissions_driver.F`  
+    a. Insert `lpbl` to the interface for `emissions_driver` in the same position as `grid%kpbl` is in `chem_driver.F`. Add the following code to line 74, underneath `dust_flux, seas_flux,`:  
+    ```fortran
+         ! bbinjectscheme requires record of layer which PBL goes up to
+         lpbl,                                                             &
+    ```
+
+    b. Declare `lpbl` as an integer. Add the following code to line 336, underneath `end stuff for lightning NOx`:  
+    ```fortran
+    ! bbinjectscheme requires record of boundary layer maximum level  
+        INTEGER,  DIMENSION( ims:ime , jms:jme ) , INTENT(in) :: lpbl  
+    ```
+
+    c. Insert `lpbl` to `do_plumerisefire`. Add the following code to line 672, underneath `emis_ant,z_at_w,z,config_flags%scale_fire_emiss,`:  
+    ```fortran
+               lpbl,                                                           &
+    ```
+
+4. `chem/module_plumerise1.F`  
+    a. Insert `lpbl` to the `plumerise_driver` in the same position as it is in `emissions_driver.F`. Add the following code to line 29, underneath `emis_ant,z_at_w,z,scale_fire_emiss,`:  
+    ```fortran
+               lpbl,                                              &           
+    ```  
+
+    b. Declare `lpbl` as an integer. Add the following code to line 53, underneath `its,ite, jts,jte, kts,kte`:  
+    ```fortran
+    ! bbinjectscheme requires record of boundary layer maximum level  
+       INTEGER,  DIMENSION( ims:ime , jms:jme ) , INTENT(in) :: lpbl  
+    ```
+
+    c. Declare real variables. Add the following code to line 96, underneath `,z_lev`:  
+    ```fortran
+          real, dimension (num_ebu) :: ebu_sum 
+	      real :: pbl_depth, flame_frac, emiss_frac
+    ```
+
+    d. Insert the conditionals which control the different injection options. Add the following code to line 251, underneath `if( maxval( eburn_in(:) ) == 0. ) cycle`:  
+    ```fortran
+    !++ SAN, 2015-04-08: adding namelist option to turn on/off plumerise calculations
+               if(config_flags%bbinjectscheme > 0) then
+
+                 ! If we are not using the plumerise routine then we need to make sure
+                 ! we deal with all emissions, not just the smouldering emissions at 
+                 ! level 1.
+
+                 !!!!!! scale_fire_emiss = .false. !!!!!!!!
+                 ! 3BEM emissions contain only the smouldering emissions, so we need to
+                 ! determine the amount of flaming emissions, and distribute that as directed.
+
+                 !!!!!! scale_fire_emiss = .true.  !!!!!!!!!
+                 ! FINN (and other emissions?) contain the whole set of emissions, so
+                 ! we already have the total emissions input. What we will do here is 
+                 ! calculate what the flaming and smouldering emissions parts of this are,
+                 ! separate them out, and then process the flaming emissions in the same
+                 ! manner as above.
+                 !
+                 ! This will mean that we don't need to use the code after the plumerise
+                 ! call which rescales the emissions once they have been generated.
+                 if(scale_fire_emiss == .false.) then
+
+                    ! calculate what the flaming emissions should be
+                    do nv=1,num_ebu
+                        ebu_sum(nv) = ebu(i,kts,j,nv) *     &
+                                        (mean_fct(1) + mean_fct(2) + mean_fct(3) + mean_fct(4))
+                    enddo
+
+                 else  !! need to account for scaling of fire emissions
+
+                    !! calculate the fraction of full emissions that should be flaming emissions
+                    !! emiss_frac = (mean_fct(1) + mean_fct(2) + mean_fct(3) + mean_fct(4)) / & !*** CHANGE ****
+                    !!						(1 + mean_fct(1) + mean_fct(2) + mean_fct(3) + mean_fct(4)) !*** CHANGE ****
+                    emiss_frac = 1.0  !*** CHANGE ****
+                    do nv=1,num_ebu
+                        !! calculate the amount of flaming emissions, to use in calculations below
+                        ebu_sum(nv) = ebu(i,kts,j,nv) * emiss_frac
+                        !! replace emission amount at ground level with proper smouldering emissions
+                        ebu(i,kts,j,nv) = ebu(i,kts,j,nv) * (1.0 - emiss_frac)
+                    end do
+
+                 end if ! check for scale_fire_emiss
+
+
+                    !! DL (23/7/2015) Modify the biomass burning injection heights
+                    !! DL (13/5/2016)  adding options 3 & 4
+                    !!   bbinjectscheme:      
+                    !!			1 = all emissions go in at ground level
+                    !!			2 = flaming emissions are distributed evenly through the boundary layer
+                    !!          3 = flaming emissions are injected at the top of the boundary layer
+                    !!          4 = flaming emissions are injected at a predetermined height
+
+
+                    !!! option 1
+                    if(config_flags%bbinjectscheme == 1)then ! everything goes into the bottom level 
+                        do nv=1,num_ebu
+                            ebu(i,kts,j,nv) = ebu(i,kts,j,nv) + ebu_sum(nv)
+                        enddo					
+
+                    !!! option 2
+                    else if(config_flags%bbinjectscheme == 2 .and. lpbl(i,j) == kts )then ! everything goes into the bottom level (as PBL very low)
+                        do nv=1,num_ebu
+                            ebu(i,kts,j,nv) = ebu(i,kts,j,nv) + ebu_sum(nv) !*** CHANGE ****
+                                                    !!ebu(i,kts,j,nv) = ebu(i,kts,j,nv)                !*** CHANGE ****
+                        enddo					
+                    else if(config_flags%bbinjectscheme == 2 .and. lpbl(i,j) > kts )then ! distribute flaming emissions through the PBL
+                        ! determine division of emissions through PBL
+                        !   1) get full depth of PBL (minus the depth of lowest model layer)
+                        !              z_at_w is the height of the lower boundary of a model layer,
+                        !              so to get top of the PBL layer we have to use the bottom of the layer above
+                        !! pbl_depth = z_at_w(i,lpbl(i,j)+1,j) - z_at_w(i,kts+1,j) !*** CHANGE ****
+                                            pbl_depth = z_at_w(i,lpbl(i,j)+1,j) - z_at_w(i,kts,j)    !*** CHANGE **** To get full depth including depth of lowest layer, not minus lowest layer
+
+                        ! loop through layers of the PBL
+                        !! do k=kts+1,lpbl(i,j) !*** CHANGE ****
+                                            do k=kts,lpbl(i,j)   !*** CHANGE **** include loop from lowest layer
+                            !   2) calculate the fraction of this height that each layer is
+                            flame_frac = (z_at_w(i,k+1,j) - z_at_w(i,k,j)) / pbl_depth
+
+
+                            !   3) multiple flaming emission total by this fraction to get injection into layer
+                            do nv=1,num_ebu
+                                ebu(i,k,j,nv) = ebu_sum(nv) * flame_frac !*** CHANGE ****
+                                                            !!ebu(i,k,j,nv) = ebu(i,k,j,nv) * flame_frac  !*** CHANGE ****
+                            end do				
+                        end do
+
+                    !!! option 3
+                    else if(config_flags%bbinjectscheme == 3 )then ! insert flaming emissions in the PBL model level (doesn't matter where this is)
+
+                        do nv=1,num_ebu
+                            ebu(i,lpbl(i,j),j,nv) = ebu(i,lpbl(i,j),j,nv) + ebu_sum(nv)
+                        enddo					
+
+
+                    !!! option 4
+    !				else if(config_flags%bbinjectscheme == 4 )then ! insert flaming emissions at a predetermined height
+    !
+    !					do nv=1,num_ebu
+    !						ebu(i,lpbl(i,j),j,nv) = ebu(i,lpbl(i,j),j,nv) + ebu_sum(nv)
+    !					enddo					
+
+
+                    else
+                        call wrf_error_fatal("bbinjectscheme setting is unsupported - check your config file")
+                    end if
+
+
+
+               else 
+    !-- SAN. Only do following calculations if plumerise is on:
+    ```
+
+    e. End the conditional. Add the following code to line 544, underneath `end if has_total_emissions`:  
+    ```fortran
+    !++ SAN, 2015-04-08
+                endif ! plumerise_off
+    !-- SAN
+    ```
+
+5. `WRFotron/namelist.chem/blueprint`  
+    a. Define the `bbinjectscheme` setting in `namelist.chem.blueprint`. Add the following line to the bottom of the namelist:  
+    ```fortran
+     bbinjectscheme                      = 2, 2, 2,                   ! 0 = use plume rise under biomass_burn_opt, 1 = all ground level, 2 = flaming evenly in BL (recommened), 3 = flaming injected top of BL, 4 = flaming injected at specific height
+    ```
+
+6. [Recompile](https://github.com/wrfchem-leeds/WRFotron/blob/master/WRFotron_user_guide.md#compile)  
+

guides/changing_domain.md

+# Changing the WRF-Chem domain
+
+### Viewing the current domain
+WRFotron includes a script which makes a plot of the domain, `plotgrids.ncl`.
+'ncl' stands for NCAR Command Language and is another programming language like python,
+but hopefully you will never need to learn it.
+1. On ARC4, load ncl by running `module load ncl`
+2. Run the script: `ncl plotgrids.ncl`. This will create a file, `wps_show_dom.pdf`
+3. Open the file. You can open files on linux with the default program by running
+`xdg-open`. Try   
+`xdg-open wps_show_dom.pdf`. You should see the default domain which includes
+all of China.
+
+### Changing the domains
+To change the domain, we need to edit several variables in two of the namelist files,
+`namelist.wps.blueprint` and `namelist.wrf.blueprint`.
+
+`namelist.wps.blueprint` contains the following lines which define the domain:
+```bash
+&geogrid                                              ! geogrid
+ parent_id         =   1, 1,                            ! domain number of the nest’s parent
+ parent_grid_ratio =   1, 3,                            ! nesting ratio relative to the domain’s parent
+ i_parent_start    =   1, 35,                           ! x coordinate of the lower-left corner
+ j_parent_start    =   1, 10,                           ! y coordinate of the lower-left corner
+ e_we              =  170, 220,                         ! westeast dimension
+ e_sn              =  170, 120,                         ! southnorth dimension
+ geog_data_res     = 'default', 'default',              ! resolution, or list of resolutions separated by + symbols, of source data to be used when interpolating static terrestrial data to the nest’s grid
+ dx = 30000,                                            ! westeast resolution, in metres for 'polar', 'lambert', and 'mercator'
+ dy = 30000,                                            ! westeast resolution, in metres for 'polar', 'lambert', and 'mercator'
+ map_proj = 'lambert',                                  ! map projection
+ ref_lat   =  33.0,                                     ! center latitude
+ ref_lon   =  103.0,                                    ! centre longitude
+ truelat1  =  23.0,                                     ! the first true latitude for the Lambert conformal projection, or the only true latitude for the Mercator and polar stereographic projections
+ truelat2  =  43.0,                                     ! the second true latitude for the Lambert conformal conic projection. For all other projections, truelat2 is ignored
+ stand_lon =  103.0,                                    ! the longitude that is parallel with the y-axis in the Lambert conformal and polar stereographic projections
+ geog_data_path = '__geogDir__'                         ! path to the static geography files
+/
+```
+
+The following variables should be modified when changing domain:  
+1. Change `ref_lat` and `ref_lon` to move the centre of the domain.
+2. Change the size of the domain by changing `dx` and `dy` as well as `e_we` and `e_sn`. `dx`
+and `dy` change the size of grid cells in the x (longitude) and y (latitude) directions.
+Changing them will change the model resolution, but will not affect the run time,
+just how far apart the grid cells are spaced. `e_we` and `e_sn` change the number of grid cells
+across the longitude and latitude of the domain, respectively.
+Adding or removing gridcells can greatly increase/decrease the run time.
+3. Change `truelat1` and `truelat2`. These should be equidistant to the centre
+and top/bottom of the domain, respectively.
+`stand_lon` should be the same as `ref_lon`
+4. Finally, rerun `plotgrids.ncl` and check that the domain has been adjusted correctly.
+
+`namelist.wrf.blueprint` also contains the variables
+`e_we`, `e_sn`, `dx` and `dy`. The values need to be changed so they are the same as
+in  `namelist.wps.blueprint`.

guides/postprocessing_guide.md

+---
+
+# WRF-Chem post-processing instructions
+
+
+### File tree
+Here is a diagram showing the layout of the files that is created in your nobackup when you run WRF-Chem
+```bash
+./simulation_WRFChem4.2_test # can be changed for each project
+└── output
+│   └── base # where the post processed output files are
+└── restart
+│   └── base # where the restart file should be moved to
+└── run
+│   └── base
+│   │   └── 2015-10-11_18:00:00-2015-10-13_00:00:00 # where the run happens
+│   │   └── staging # where wrfout files are copied before they are post-processed
+```
+
+### Checking the output
+Your test simulation should now have finished. First we will check that all the 'wrfout' files were created. On ARC4, go to the run directory of the test simulation (`simulation_WRFChem4.2_test`).
+```bash
+cd /nobackup/${USER}/simulation_WRFChem4.2_test/run/base/2015-10-11_18:00:00-2015-10-13_00:00:00
+```
+
+Then you can use `ls` to list the wrfout files, and you can count the wrfout files using `ls | wc`
+```bash
+ls wrfout_d01_2015-10-1*
+ls wrfout_d01_2015-10-1* | wc
+```
+
+The test run should generate 31 output files. This is because the meteorology spinup was 6 hours, and the run was 24, and an extra hour is added so that a restart file for the next day (`wrfout_d01_2015-10-13_00:00:00`) can be generated. This would allow us to restart the run from the next day if necessary.
+
+Then we need to check whether `post.bash` has run successfully. `post.bash` compresses the files to reduce their size. At the end of `main.bash`, the wrfout files are moved to the staging directory (`/simulation_WRFChem4.2_test/run/base/staging`). If you navigate there, you should see the wrfout files there (but not the ones generated during the spinup).
+
+Then, `post.bash` outputs postprocessed files to `/simulation_WRFChem4.2_test/output/base`. There should be 23 wrfout files there (It should be 24 but there is a bug in the code, I have sent a [pull request](https://github.com/wrfchem-leeds/WRFotron/pull/46)). If you have 23 wrfout files in your `/output/base` directory, move onto the next step
+
+### Regridding and concatenating the output
+
+This process used to be a bit of a pain but now we have Luke's `pp_concat_regrid.py` script. This script does two things:
+1. It concatenates the hourly wrfout files into a single file, i.e. the 2d arrays (lat, lon) are joined into a 3d array (time, lat, lon).
+2. It regrids the arrays from a curvilinear grid to a rectilinear/regular grid. This makes the data easier to analyse and compare with other datasets.
+
+First, please make sure you have these lines in your `~/.bashrc` file:
+```bash
+if [ -r /nobackup/cemac/cemac.sh ] ; then
+  . /nobackup/cemac/cemac.sh
+fi
+```
+
+You will see that `pp_concat_regrid.py` is already in the `/output/base` directory. There is also the `pp_concat_regrid.bash` script, which is used to submit `pp_concat_regrid.py` to the supercomputer queue and load the python modules needed to run the script.
+
+Before submitting `pp_concat_regrid.py` you need to change a couple of variables in the file. Open `pp_concat_regrid.py` using `vim` or `gedit`, and you will see some variables that you can define at the top of the file, you will need to change the `year` variable as the test run was in 2015 not 2016
+
+```python
+# setup - ensure meets your requirements
+year = "2016" # change this to 2015
+month = "10"
+domains = ["1"]
+variables = ["PM2_5_DRY", "o3", "AOD550_sfc"]
+aerosols = ["bc", "oc", "nh4", "so4", "no3", "asoaX", "bsoaX", "oin"] 
+variables.extend(aerosols)
+surface_only = True
+res = 0.25
+regrid = True
+convert_aerosol_units_from_ugkg_to_ugm3 = True
+```
+
+Now submit the regridding script `qsub pp_concat_regrid.bash`. The script may take around an hour to run. If it finishes quickly, there was probably an error. If this happens, the error will be in the `pp_concat_regrid.bash.eXXXXXX` file. If the script is working, it will create these files:
+
+```bash
+wrfout_d01_global_0.25deg_2015-10_AOD550_sfc.nc  wrfout_d01_global_0.25deg_2015-10_nh4_2p5.nc  wrfout_d01_global_0.25deg_2015-10_oin_2p5.nc
+wrfout_d01_global_0.25deg_2015-10_asoaX_2p5.nc   wrfout_d01_global_0.25deg_2015-10_no3_2p5.nc  wrfout_d01_global_0.25deg_2015-10_PM2_5_DRY.nc
+wrfout_d01_global_0.25deg_2015-10_bc_2p5.nc      wrfout_d01_global_0.25deg_2015-10_o3.nc       wrfout_d01_global_0.25deg_2015-10_so4_2p5.nc
+wrfout_d01_global_0.25deg_2015-10_bsoaX_2p5.nc   wrfout_d01_global_0.25deg_2015-10_oc_2p5.nc
+```