| 1 | #------------------------------------------------------------------------------# |
|---|
| 2 | # This file is part of the PALM model system. |
|---|
| 3 | # |
|---|
| 4 | # PALM is free software: you can redistribute it and/or modify it under the |
|---|
| 5 | # terms of the GNU General Public License as published by the Free Software |
|---|
| 6 | # Foundation, either version 3 of the License, or (at your option) any later |
|---|
| 7 | # version. |
|---|
| 8 | # |
|---|
| 9 | # PALM is distributed in the hope that it will be useful, but WITHOUT ANY |
|---|
| 10 | # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR |
|---|
| 11 | # A PARTICULAR PURPOSE. See the GNU General Public License for more details. |
|---|
| 12 | # |
|---|
| 13 | # You should have received a copy of the GNU General Public License along with |
|---|
| 14 | # PALM. If not, see <http://www.gnu.org/licenses/>. |
|---|
| 15 | # |
|---|
| 16 | # Copyright 2019-2019 University of Helsinki |
|---|
| 17 | #------------------------------------------------------------------------------# |
|---|
| 18 | |
|---|
| 19 | import math |
|---|
| 20 | from netCDF4 import Dataset |
|---|
| 21 | import numpy as np |
|---|
| 22 | import datetime |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | class SalsaDriver: |
|---|
| 26 | """ This is an example script to generate salsa drivers for PALM. |
|---|
| 27 | |
|---|
| 28 | You can use it as a starting point for creating your setup specific |
|---|
| 29 | driver. |
|---|
| 30 | """ |
|---|
| 31 | |
|---|
| 32 | |
|---|
| 33 | def __init__(self): |
|---|
| 34 | """ Open the salsa driver as NetCDF4 file. Here, you have to give the |
|---|
| 35 | full path to the static driver that shall be created. Existing file |
|---|
| 36 | with same name is deleted. |
|---|
| 37 | """ |
|---|
| 38 | print('Opening file...') |
|---|
| 39 | self.nc_file = Dataset('salsa_driver.nc', 'w', format='NETCDF4') |
|---|
| 40 | |
|---|
| 41 | |
|---|
| 42 | def write_global_attributes(self): |
|---|
| 43 | """ Write global attributes to static driver. """ |
|---|
| 44 | print("Writing global attributes...") |
|---|
| 45 | |
|---|
| 46 | # optional global attributes |
|---|
| 47 | self.nc_file.origin_lon = 55.0 |
|---|
| 48 | self.nc_file.origin_lat = 0.0 |
|---|
| 49 | self.nc_file.origin_time = '2000-06-21 12:00:00 +00' |
|---|
| 50 | self.nc_file.origin_x = 308124 |
|---|
| 51 | self.nc_file.origin_y = 6098908 |
|---|
| 52 | self.nc_file.origin_z = 0.0 |
|---|
| 53 | self.nc_file.rotation_angle = 0.0 |
|---|
| 54 | self.nc_file.author = 'Your Name' |
|---|
| 55 | self.nc_file.comment = 'Miscellaneous information about the data ' \ |
|---|
| 56 | 'or methods to produce it.' |
|---|
| 57 | self.nc_file.creation_date = str(datetime.datetime.now()) |
|---|
| 58 | self.nc_file.institution = 'INAR/Physics, University of Helsinki' |
|---|
| 59 | self.nc_file.history = '' |
|---|
| 60 | self.nc_file.palm_revision = '' |
|---|
| 61 | self.nc_file.title = 'Salsa driver for some arbitrary PALM setup' |
|---|
| 62 | |
|---|
| 63 | |
|---|
| 64 | def define_dimensions(self): |
|---|
| 65 | """ Set dimensions on which variables are defined. """ |
|---|
| 66 | print("Writing dimensions...") |
|---|
| 67 | |
|---|
| 68 | # General grid parameters: |
|---|
| 69 | |
|---|
| 70 | self.nx = 19 |
|---|
| 71 | self.ny = 19 |
|---|
| 72 | self.nz = 20 |
|---|
| 73 | dx = 2 |
|---|
| 74 | dy = 2 |
|---|
| 75 | dz = 2 |
|---|
| 76 | |
|---|
| 77 | # time steps in the emission data: |
|---|
| 78 | self.times = np.arange( 0.0, 7201.0, 3600.0 ) |
|---|
| 79 | |
|---|
| 80 | |
|---|
| 81 | # Aerosol emissions: |
|---|
| 82 | |
|---|
| 83 | # number of emission categories |
|---|
| 84 | self.nncat = 3 |
|---|
| 85 | |
|---|
| 86 | # number of chemical components |
|---|
| 87 | self.ncomposition_index = 7 |
|---|
| 88 | |
|---|
| 89 | # level of detail of aerosol emissions: |
|---|
| 90 | self.lod_aerosol_emission = 2 # 1: yearly PM emissions per emission category |
|---|
| 91 | # 2: per emission category at given points in time and space |
|---|
| 92 | |
|---|
| 93 | # level of detail for scaling lod=1 emissions |
|---|
| 94 | self.lod_emission_time_factor = 1 # 1: scaling according to month-day-hour, |
|---|
| 95 | # 2: scaling according to hour of year |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | # Aerosol size distribution: |
|---|
| 99 | |
|---|
| 100 | # number of aerosol size bins in the subrange 1 and 2 |
|---|
| 101 | nbin = [1, 7] |
|---|
| 102 | |
|---|
| 103 | # subrange diameter limit (e.g. subrange 1: 3-10 nm, subrange 2: 10 nm - 2.5 um |
|---|
| 104 | reglim = [3.0e-9, 10.0e-9, 2.5e-6] |
|---|
| 105 | |
|---|
| 106 | |
|---|
| 107 | # Mandatory dimensions |
|---|
| 108 | self.nmax_string_length = 25 |
|---|
| 109 | self.nnhoursyear = np.arange( 1, 8760+1, 1 ) |
|---|
| 110 | self.nnmonthdayhour = np.arange( 1, 91+1, 1 ) |
|---|
| 111 | |
|---|
| 112 | self.ntime = len( self.times ) |
|---|
| 113 | |
|---|
| 114 | self.nc_file.createDimension('x' ,self.nx+1) |
|---|
| 115 | self.x = self.nc_file.createVariable('x', 'f8', ('x',)) |
|---|
| 116 | self.x.units = 'm' |
|---|
| 117 | self.x.standard_name = 'x coordinate of cell centers' |
|---|
| 118 | self.x[:] = np.arange( 0.5*dx, ( self.nx+1 ) * dx, dx ) |
|---|
| 119 | |
|---|
| 120 | self.nc_file.createDimension('y', self.ny+1) |
|---|
| 121 | self.y = self.nc_file.createVariable('y', 'f8', ('y',)) |
|---|
| 122 | self.y.units = 'm' |
|---|
| 123 | self.y.standard_name = 'y coordinate of cell centers' |
|---|
| 124 | self.y[:] = np.arange( 0.5*dy, ( self.ny+1 ) * dy, dy ) |
|---|
| 125 | |
|---|
| 126 | self.nc_file.createDimension('max_string_length', self.nmax_string_length ) |
|---|
| 127 | self.max_string_length = self.nc_file.createVariable('max_string_length', 'i4', |
|---|
| 128 | ('max_string_length',)) |
|---|
| 129 | self.max_string_length.units = '' |
|---|
| 130 | self.max_string_length[:] = np.linspace( 1, self.nmax_string_length, self.nmax_string_length ) |
|---|
| 131 | |
|---|
| 132 | self.nc_file.createDimension('ncat', self.nncat ) |
|---|
| 133 | self.ncat = self.nc_file.createVariable('ncat', 'i4', ('ncat',)) |
|---|
| 134 | self.ncat.units = '' |
|---|
| 135 | self.ncat.standard_name = 'number of emission categories' |
|---|
| 136 | self.ncat[:] = np.linspace( 1, self.nncat, self.nncat ) |
|---|
| 137 | |
|---|
| 138 | self.nc_file.createDimension('composition_index', self.ncomposition_index ) |
|---|
| 139 | self.composition_index = self.nc_file.createVariable('composition_index', 'i4', |
|---|
| 140 | ('composition_index',)) |
|---|
| 141 | self.composition_index.units = '' |
|---|
| 142 | self.composition_index.long_name = 'composition index' |
|---|
| 143 | self.composition_index[:] = np.linspace( 1, self.ncomposition_index, self.ncomposition_index ) |
|---|
| 144 | |
|---|
| 145 | |
|---|
| 146 | if self.lod_aerosol_emission==1: |
|---|
| 147 | |
|---|
| 148 | if self.lod_emission_time_factor==1: |
|---|
| 149 | self.nc_file.createDimension('nmonthdayhour', len( self.nnmonthdayhour ) ) |
|---|
| 150 | self.nmonthdayhour = self.nc_file.createVariable('nmonthdayhour', 'i4', ('nmonthdayhour',)) |
|---|
| 151 | self.nmonthdayhour.units = '' |
|---|
| 152 | self.nmonthdayhour.standard_name = 'number of required input values for emission time '\ |
|---|
| 153 | 'rescaling factors' |
|---|
| 154 | self.nmonthdayhour[:] = self.nnmonthdayhour |
|---|
| 155 | |
|---|
| 156 | elif self.lod_emission_time_factor==2: |
|---|
| 157 | self.nc_file.createDimension('nhoursyear', len( self.nnhoursyear ) ) |
|---|
| 158 | self.nhoursyear = self.nc_file.createVariable('nhoursyear', 'i4', ('nhoursyear',)) |
|---|
| 159 | self.nhoursyear.units = '' |
|---|
| 160 | self.nhoursyear.standard_name = 'number of hours per year' |
|---|
| 161 | self.nhoursyear[:] = self.nnhoursyear |
|---|
| 162 | |
|---|
| 163 | elif self.lod_aerosol_emission==2: |
|---|
| 164 | self.nc_file.createDimension('time' , self.ntime ) |
|---|
| 165 | self.time = self.nc_file.createVariable('time', 'f8', ('time',)) |
|---|
| 166 | self.time.units = 's' |
|---|
| 167 | self.time.standard_name = 'time since utc init' |
|---|
| 168 | self.time[:] = self.times |
|---|
| 169 | |
|---|
| 170 | self.ndmid, self.bin_limits = define_bins( nbin, reglim ) |
|---|
| 171 | self.nc_file.createDimension('Dmid' , np.sum( nbin ) ) |
|---|
| 172 | self.Dmid = self.nc_file.createVariable('Dmid', 'f8', ('Dmid',)) |
|---|
| 173 | self.Dmid.units = 'm' |
|---|
| 174 | self.Dmid.standard_name = 'mean diamater per size bin' |
|---|
| 175 | self.Dmid[:] = self.ndmid |
|---|
| 176 | |
|---|
| 177 | |
|---|
| 178 | def add_variables(self): |
|---|
| 179 | """ Uncomment variables below as you like. |
|---|
| 180 | |
|---|
| 181 | Be aware that some variables depend on others. For a description of |
|---|
| 182 | each variable, please have a look at the documentation at |
|---|
| 183 | https://palm.muk.uni-hannover.de/trac/wiki/doc/app/iofiles/pids/aerosol |
|---|
| 184 | |
|---|
| 185 | An example of how you modify the variables is given below: |
|---|
| 186 | |
|---|
| 187 | building_2d_array = np.ones((self.ny+1,self.nx+1)) * -9999.0 |
|---|
| 188 | south_wall, north_wall, left_wall, right_wall = 20, 25, 20, 25 |
|---|
| 189 | building_2d_array[south_wall:north_wall,left_wall:right_wall] = 50 |
|---|
| 190 | nc_buildings_2d = self.nc_file.createVariable( |
|---|
| 191 | 'buildings_2d', 'f4', ('y','x'),fill_value=-9999.0) |
|---|
| 192 | nc_buildings_2d.lod = 1 |
|---|
| 193 | nc_buildings_2d[:,:] = building_2d_array |
|---|
| 194 | """ |
|---|
| 195 | print("Writing variables...") |
|---|
| 196 | |
|---|
| 197 | emission_category_name_list = ['traffic exhaust ', |
|---|
| 198 | 'road dust ', |
|---|
| 199 | 'wood combustion '] |
|---|
| 200 | composition_name_list = ['H2SO4 ','OC ', |
|---|
| 201 | 'BC ','DU ', |
|---|
| 202 | 'SS ','HNO3 ', |
|---|
| 203 | 'NH3 '] |
|---|
| 204 | |
|---|
| 205 | mass_fracs_array = np.zeros( [ self.nncat, self.ncomposition_index ] ) |
|---|
| 206 | mass_fracs_array[0,:] = np.array( [0.04, 0.48, 0.48, 0.0, 0.0, 0.0, 0.0] ) # traffic |
|---|
| 207 | mass_fracs_array[1,:] = np.array( [0.0, 0.05, 0.0, 0.85, 0.0, 0.05, 0.05] ) # road dust |
|---|
| 208 | mass_fracs_array[2,:] = np.array( [0.0, 0.5, 0.5, 0.0, 0.0, 0.0, 0.0] ) # wood combustion |
|---|
| 209 | |
|---|
| 210 | esh = np.ones( ( self.ny+1, self.nx+1 ) ) * -9999.0 # emission stack height (e.g. chimneys) |
|---|
| 211 | esh[1:3,1:3] = 10.0 |
|---|
| 212 | esh[15:18,1:3] = 20.0 |
|---|
| 213 | |
|---|
| 214 | road = np.zeros( ( self.ny+1, self.nx+1 ) ) |
|---|
| 215 | road[:,int(self.nx/2-2):int(self.nx/2+4)] = 1 |
|---|
| 216 | |
|---|
| 217 | if self.lod_aerosol_emission==1: |
|---|
| 218 | if self.lod_emission_time_factor==1: |
|---|
| 219 | # emission time factors: |
|---|
| 220 | etf = np.zeros( [self.nncat, len( self.nnmonthdayhour )] ) |
|---|
| 221 | # traffic exhaust and dust |
|---|
| 222 | etf[0:2,0:12] = 1.0/12 # months |
|---|
| 223 | etf[0:2,12:17] = 0.17; etf[0:2,17:19] = 0.075 # days of week |
|---|
| 224 | etf[0:2,19:43] = 0.03; etf[0:2,19+7:19+10] += 0.04; etf[0:2,19+15:19+19] += 0.04 # working day hours |
|---|
| 225 | etf[0:2,43:67] = 1.0/24 # hours of a saturday |
|---|
| 226 | etf[0:2,67:91] = 1.0/24 # hours of a sunday |
|---|
| 227 | # wood combustion |
|---|
| 228 | etf[2,0:3] = 0.2; etf[2,10:12] = 0.2 # months |
|---|
| 229 | etf[2,12:19] = 1.0/7 # days of week |
|---|
| 230 | etf[2,19+18:19+23] = 0.2 # hours of a working day |
|---|
| 231 | etf[2,43+18:43+23] = 0.2 # hours of a saturday |
|---|
| 232 | etf[2,67+18:67+23] = 0.2 # hours of a saturday |
|---|
| 233 | |
|---|
| 234 | else: |
|---|
| 235 | etf = np.zeros( [self.nncat, len( self.nnhoursyear )] ) |
|---|
| 236 | # For now, set everything to an equal value |
|---|
| 237 | etf[:,:] = 1.0 / len( self.nnhoursyear ) |
|---|
| 238 | |
|---|
| 239 | aerosol_emission_values = np.ones( ( self.ny+1, self.nx+1, self.nncat ) ) * -9999.0 |
|---|
| 240 | i = emission_category_name_list.index('traffic exhaust ') |
|---|
| 241 | aerosol_emission_values[road>0,i] = 8.0 # g/m2/year |
|---|
| 242 | i = emission_category_name_list.index('road dust ') |
|---|
| 243 | aerosol_emission_values[road>0,i] = 4.0 # g/m2/year |
|---|
| 244 | i = emission_category_name_list.index('wood combustion ') |
|---|
| 245 | aerosol_emission_values[esh>0,i] = 2.0 # g/m2/year |
|---|
| 246 | |
|---|
| 247 | elif self.lod_aerosol_emission==2: |
|---|
| 248 | # sectional size distribution |
|---|
| 249 | nsect = np.zeros( [self.nncat, len( self.ndmid ) ], dtype=float ) |
|---|
| 250 | relative_per_bin = np.copy( nsect ) # sum over all bins = 1 |
|---|
| 251 | for n in range( len( self.ncat ) ): |
|---|
| 252 | |
|---|
| 253 | if 'traffic exhaust' in emission_category_name_list[n]: |
|---|
| 254 | dpg = np.array([20.3, 72.0]) |
|---|
| 255 | n_lognorm = np.array([18960.0, 13750.0]) |
|---|
| 256 | sigmag = np.array([1.7, 1.6]) |
|---|
| 257 | |
|---|
| 258 | elif 'road dust' in emission_category_name_list[n]: |
|---|
| 259 | dpg = np.array([1400.0]) |
|---|
| 260 | n_lognorm = np.array([4.0]) |
|---|
| 261 | sigmag = np.array([1.4]) |
|---|
| 262 | |
|---|
| 263 | elif 'wood combustion' in emission_category_name_list[n]: |
|---|
| 264 | dpg = np.array([540.0]) |
|---|
| 265 | n_lognorm = np.array([5000.0]) |
|---|
| 266 | sigmag = np.array([1.7]) |
|---|
| 267 | |
|---|
| 268 | for l in range( 1, len( self.ndmid )+1 ): |
|---|
| 269 | d1 = self.bin_limits[l-1] |
|---|
| 270 | d2 = self.bin_limits[l] |
|---|
| 271 | delta_d = ( d2 - d1 ) / 10.0 |
|---|
| 272 | for ib in range( 1, len( self.ndmid )+1 ): |
|---|
| 273 | d1 = self.bin_limits[l-1] + ( ib - 1 ) * delta_d |
|---|
| 274 | d2 = d1 + delta_d |
|---|
| 275 | dmidi = ( d1 + d2 ) / 2.0 |
|---|
| 276 | deltadp = np.log10( d2 / d1 ) |
|---|
| 277 | nsect[n,l-1] = nsect[n,l-1] + np.sum( n_lognorm * 1.0E6 * deltadp / |
|---|
| 278 | ( np.sqrt( 2.0 * np.pi ) * np.log10( sigmag ) ) * |
|---|
| 279 | np.exp( -np.log10( dmidi / ( 1.0E-9 * dpg ) )**2.0 / |
|---|
| 280 | ( 2.0 * np.log10( sigmag ) ** 2.0 ) ) ) |
|---|
| 281 | |
|---|
| 282 | relative_per_bin[n,:] = nsect[n,:] / np.sum( nsect[n,:] ) |
|---|
| 283 | |
|---|
| 284 | aerosol_emission_values = np.ones( ( self.ntime, self.ny+1, self.nx+1, self.nncat ) ) * -9999.0 |
|---|
| 285 | for t in range( self.ntime ): |
|---|
| 286 | i = emission_category_name_list.index('traffic exhaust ') |
|---|
| 287 | aerosol_emission_values[t,road>0,i] = 1.0E+9 * (t+1) # units #/m2/s |
|---|
| 288 | i = emission_category_name_list.index('road dust ') |
|---|
| 289 | aerosol_emission_values[t,road>0,i] = 1.0E+3 * (t+1) # units #/m2/s |
|---|
| 290 | i = emission_category_name_list.index('wood combustion ') |
|---|
| 291 | aerosol_emission_values[t,esh>0,i] = 5.0E+6 # units #/m2/s |
|---|
| 292 | |
|---|
| 293 | |
|---|
| 294 | # Save into the file |
|---|
| 295 | nc_emission_category_index = self.nc_file.createVariable( 'emission_category_index','i1', |
|---|
| 296 | ('ncat',) ) |
|---|
| 297 | nc_emission_category_index[:] = np.linspace( 1, self.nncat, self.nncat ) |
|---|
| 298 | nc_emission_category_index.units = '' |
|---|
| 299 | nc_emission_category_index.long_name = 'emission_category_index' |
|---|
| 300 | |
|---|
| 301 | nc_emission_category_name = self.nc_file.createVariable( 'emission_category_name', 'S1', |
|---|
| 302 | ('ncat','max_string_length',)) |
|---|
| 303 | nc_emission_category_name[:] = list( map( lambda x : list(x), emission_category_name_list ) ) |
|---|
| 304 | nc_emission_category_name.long_name = 'emission category name' |
|---|
| 305 | nc_emission_category_name.standard_name = 'emission_category_name' |
|---|
| 306 | |
|---|
| 307 | nc_composition_name = self.nc_file.createVariable( 'composition_name', 'S1', |
|---|
| 308 | ('composition_index','max_string_length',) ) |
|---|
| 309 | nc_composition_name[:] = list( map( lambda x : list(x), composition_name_list ) ) |
|---|
| 310 | nc_composition_name.long_name = 'aerosol composition name' |
|---|
| 311 | nc_composition_name.standard_name = 'composition_name' |
|---|
| 312 | |
|---|
| 313 | nc_emission_mass_fracs = self.nc_file.createVariable( 'emission_mass_fracs', 'f4', |
|---|
| 314 | ('ncat','composition_index',), |
|---|
| 315 | fill_value=-9999.0 ) |
|---|
| 316 | nc_emission_mass_fracs.units = '' |
|---|
| 317 | nc_emission_mass_fracs[:] = mass_fracs_array |
|---|
| 318 | nc_emission_mass_fracs.long_name = 'mass fractions of chemical components in aerosol emissions' |
|---|
| 319 | nc_emission_mass_fracs.standard_name = 'emission_mass_fractions' |
|---|
| 320 | |
|---|
| 321 | if self.lod_aerosol_emission==1: |
|---|
| 322 | if self.lod_emission_time_factor==1: |
|---|
| 323 | nc_emission_time_factors = self.nc_file.createVariable( 'emission_time_factors', 'f4', |
|---|
| 324 | ('ncat','nmonthdayhour',), |
|---|
| 325 | fill_value=-9999.0 ) |
|---|
| 326 | nc_emission_time_factors.lod = 1 |
|---|
| 327 | else: |
|---|
| 328 | nc_emission_time_factors = self.nc_file.createVariable( 'emission_time_factors', 'f4', |
|---|
| 329 | ('ncat','nhoursyear',), |
|---|
| 330 | fill_value=-9999.0 ) |
|---|
| 331 | nc_emission_time_factors.lod = 2 |
|---|
| 332 | nc_emission_time_factors[:] = etf |
|---|
| 333 | nc_emission_time_factors.long_name = 'emission time scaling factors' |
|---|
| 334 | nc_emission_time_factors.standard_name = 'emission_time_factors' |
|---|
| 335 | |
|---|
| 336 | nc_emission_stack_height = self.nc_file.createVariable( 'emission_stack_height', 'f4', |
|---|
| 337 | ('y','x',), fill_value=-9999.0 ) |
|---|
| 338 | nc_emission_stack_height[:] = esh |
|---|
| 339 | nc_emission_stack_height.units = 'm' |
|---|
| 340 | nc_emission_stack_height.long_name = 'emission stack height' |
|---|
| 341 | nc_emission_stack_height.standard_name = 'emission_stack_height' |
|---|
| 342 | |
|---|
| 343 | nc_aerosol_emission_values = self.nc_file.createVariable( 'aerosol_emission_values', 'f4', |
|---|
| 344 | ('y','x','ncat',), |
|---|
| 345 | fill_value=-9999.0 ) |
|---|
| 346 | nc_aerosol_emission_values[:] = aerosol_emission_values |
|---|
| 347 | nc_aerosol_emission_values.units = 'g/m2/yr' |
|---|
| 348 | nc_aerosol_emission_values.long_name = 'aerosol emission values' |
|---|
| 349 | nc_aerosol_emission_values.standard_name = 'aerosol_emission_values' |
|---|
| 350 | nc_aerosol_emission_values.lod = 1 |
|---|
| 351 | else: |
|---|
| 352 | nc_emission_number_fracs = self.nc_file.createVariable( 'emission_number_fracs', 'f4', |
|---|
| 353 | ('ncat','Dmid',), fill_value=-9999.0 ) |
|---|
| 354 | nc_emission_number_fracs.units = '' |
|---|
| 355 | nc_emission_number_fracs[:] = relative_per_bin |
|---|
| 356 | nc_emission_number_fracs.long_name = 'number fractions of aerosol size bins in aerosol emissions' |
|---|
| 357 | nc_emission_number_fracs.standard_name = 'emission_number_fractions' |
|---|
| 358 | |
|---|
| 359 | nc_aerosol_emission_values = self.nc_file.createVariable( 'aerosol_emission_values', 'f4', |
|---|
| 360 | ('time','y','x','ncat',), |
|---|
| 361 | fill_value=-9999.0 ) |
|---|
| 362 | nc_aerosol_emission_values[:] = aerosol_emission_values |
|---|
| 363 | nc_aerosol_emission_values.units = '#/m2/s' |
|---|
| 364 | nc_aerosol_emission_values.long_name = 'aerosol emission values' |
|---|
| 365 | nc_aerosol_emission_values.standard_name = 'aerosol_emission_values' |
|---|
| 366 | nc_aerosol_emission_values.lod = 2 |
|---|
| 367 | |
|---|
| 368 | |
|---|
| 369 | def finalize(self): |
|---|
| 370 | """ Close file """ |
|---|
| 371 | print("Closing file...") |
|---|
| 372 | |
|---|
| 373 | self.nc_file.close() |
|---|
| 374 | |
|---|
| 375 | def define_bins( nbin, reglim ): |
|---|
| 376 | """ This function defines the sectional bin limits based on number of bins |
|---|
| 377 | (nbin) and diameter limits (reglim) per subrange |
|---|
| 378 | """ |
|---|
| 379 | |
|---|
| 380 | nbins = np.sum( nbin ) # = subrange 1 + subrange 2 |
|---|
| 381 | |
|---|
| 382 | # Log-normal to sectional |
|---|
| 383 | |
|---|
| 384 | vlolim = np.zeros( nbins ) |
|---|
| 385 | vhilim = np.zeros( nbins ) |
|---|
| 386 | dmid = np.zeros( nbins ) |
|---|
| 387 | bin_limits = np.zeros( nbins ) |
|---|
| 388 | |
|---|
| 389 | # Sectional bin limits |
|---|
| 390 | ratio_d = reglim[1] / reglim[0] |
|---|
| 391 | for b in range( nbin[0] ): |
|---|
| 392 | vlolim[b] = np.pi / 6.0 * ( reglim[0] * ratio_d **( float(b) / nbin[0] ) )**3 |
|---|
| 393 | vhilim[b] = np.pi / 6.0 * ( reglim[0] * ratio_d **( float(b+1) / nbin[0] ) )**3 |
|---|
| 394 | dmid[b] = np.sqrt( ( 6.0 * vhilim[b] / np.pi )**0.33333333 * ( 6.0 * vlolim[b] / np.pi |
|---|
| 395 | )**0.33333333 ) |
|---|
| 396 | |
|---|
| 397 | ratio_d = reglim[2] / reglim[1] |
|---|
| 398 | for b in np.arange( nbin[0], np.sum( nbin ),1 ): |
|---|
| 399 | c = b-nbin[0] |
|---|
| 400 | vlolim[b] = np.pi / 6.0 * ( reglim[1] * ratio_d ** ( float(c) / nbin[1] ) )**3 |
|---|
| 401 | vhilim[b] = np.pi / 6.0 * ( reglim[1] * ratio_d ** ( float(c+1) / nbin[1] ) ) ** 3 |
|---|
| 402 | dmid[b] = np.sqrt( ( 6.0 * vhilim[b] / np.pi )**0.33333333 * ( 6.0 * vlolim[b] / np.pi |
|---|
| 403 | )**0.33333333 ) |
|---|
| 404 | |
|---|
| 405 | bin_limits = ( 6.0 * vlolim / np.pi )**0.33333333 |
|---|
| 406 | bin_limits = np.append( bin_limits, reglim[-1] ) |
|---|
| 407 | |
|---|
| 408 | return dmid, bin_limits |
|---|
| 409 | |
|---|
| 410 | if __name__ == '__main__': |
|---|
| 411 | driver = SalsaDriver() |
|---|
| 412 | driver.write_global_attributes() |
|---|
| 413 | driver.define_dimensions() |
|---|
| 414 | driver.add_variables() |
|---|
| 415 | driver.finalize() |
|---|
| 416 | |
|---|
