diff --git a/configurations/augsburg_validation_summer_10.conf b/configurations/augsburg_validation_summer_10.conf index 26d0b71c601c27f17dc4328347550f60716d93e3..721faae30bb23e3b7b66153a21e6fc80d0bab6d4 100644 --- a/configurations/augsburg_validation_summer_10.conf +++ b/configurations/augsburg_validation_summer_10.conf @@ -8,11 +8,11 @@ scenario = 'validation_summer' nested_domain = False # dynamic driver output -dynamic_driver_file = "/cfs/home/d/u/dupreeda/MBEES/PALM/palm_model_system-v21.10/JOBS/augs11/INPUT/augs11_dynamic" +dynamic_driver_file = "/cfs/home/d/u/dupreeda/MBEES/PALM/palm_model_system-v21.10/JOBS/augs12/INPUT/augs12_dynamic" # import grid parameters for dynamic driver from static driver grid_from_static = True # static driver input -static_driver_file = "/cfs/home/d/u/dupreeda/MBEES/PALM/palm_model_system-v21.10/JOBS/augs11/INPUT/augs11_static" +static_driver_file = "/cfs/home/d/u/dupreeda/MBEES/PALM/palm_model_system-v21.10/JOBS/augs12/INPUT/augs12_static" # reference coordinate system of PALM simulation proj_palm = "EPSG:32633" # projection lon-lat @@ -38,7 +38,16 @@ interp_dir_name = '/cfs/home/d/u/dupreeda/MBEES/PALM/wrf_chem_data/interp' # WRF-chem species need for chemical reaction wrfchem_spec = ['no','no2','o3'] -# possible species: no, no2, no3, pm10,PM2_5_DRY, o3, co, hno3, ho, h2o2 +# possible species: no, no2, no3, pm10,PM2_5_DRY, o3, co, hno3, ho, h2o2, nh3 + +# Aerosol profile - SALSA parameters +aerosol_wrfchem = True +wrfchem_bin_limits = [3.9e-8, 1.56e-7, 6.25e-7, 2.5e-6, 1.0e-5] +listspec = ['SO4', 'NH', 'NO'] +# possibilities for listspec = ['SO4', 'OC', 'BC', 'DU', 'SS', 'NH', 'NO'] +nbin = [1,7] +reglim = [3.9e-8, 5.0e-8, 2.5e-6] +nf2a = 1.0 # insoluble not supported in PALM # radiation radiation_from_wrf = False diff --git a/dynamic/palm_dynamic.py b/dynamic/palm_dynamic.py index 83e1799951da45eeda248e6e804ca31d16bfed7d..326d2d752993742b5223e4fbc8dce1be84e9fcac 100644 --- a/dynamic/palm_dynamic.py +++ b/dynamic/palm_dynamic.py @@ -32,9 +32,9 @@ import glob import netCDF4 import palm_dynamic_config +import palm_dynamic_aerosol ################################## # read configuration from command -# line and parse it def print_help(): print() print(__doc__) @@ -73,8 +73,16 @@ if configname == '': palm_dynamic_config.configure(configname) # Import values (including loaded) from config module into globals from palm_dynamic_config import * -# Other modules are imported *AFTER* the config module has been filled with -# values (so that they can correctly import * from palm_dynamic_config). +# Load all aerosol chemical species +if aerosol_wrfchem: + wrfchem_aerosols = [] + for aero in listspec: + _aero = palm_dynamic_aerosol.translate_aerosol_species(aero).split(",") + for _aeros in _aero: + wrfchem_aerosols.append(_aeros+ '_a01') + wrfchem_aerosols.append(_aeros+ '_a02') + wrfchem_aerosols.append(_aeros+ '_a03') + wrfchem_aerosols.append(_aeros+ '_a04') import palm_wrf_utils from palm_dynamic_output import palm_dynamic_output @@ -88,6 +96,10 @@ print('WRF-CHEM file path:', wrf_dir_name) print('WRF-CHEM: dynamics file mask:', wrf_file_mask) print('Simulation start time:', origin_time) print('Simulation hours:', simulation_hours) +if aerosol_wrfchem: + print('\nNumber of aerosol size bins in subranges (nbin): ', nbin) + print('Aerosol diameter limits of subrnages (reglim): ', reglim) + print('Aerosol chemical components (listspec): ', listspec) ########################### # domain parameters print('\nDomain parameters:') @@ -173,7 +185,6 @@ if grid_from_static: ncs.close() else: #TODO check all necessary parameters are set and are correct - # # initialize necessary arrays try: terrain = np.zeros(shape=(ny, nx), dtype=float) @@ -359,6 +370,10 @@ for wrf_file in wrf_files_proc: # copied vars wrfchem_dynamic = ['PH', 'PHB', 'HGT', 'T', 'W', 'TSLB', 'SMOIS', 'MU', 'MUB','P', 'PB', 'PSFC'] wrfchem_variables = wrfchem_dynamic + wrfchem_spec + if aerosol_wrfchem: + wrfchem_variables = wrfchem_variables + wrfchem_aerosols + wrfchem_variables.append('ALT') + for varname in wrfchem_variables: v_wrf = f_wrf.variables[varname] v_out = f_out.createVariable(varname, 'f4', v_wrf.dimensions) @@ -384,7 +399,7 @@ for wrf_file in wrf_files_proc: print('Vertical interpolation...') palm_wrf_utils.palm_wrf_vertical_interp(hinterp_file, vinterp_file, wrf_file, z_levels, z_levels_stag, z_soil_levels, origin_z, terrain, - wrf_hybrid_levs, vinterp_terrain_smoothing) + wrf_hybrid_levs, vinterp_terrain_smoothing, nz, ny, nx) if radiation_from_wrf: print('Start processing of radiation inputs from the WRF radiation files.') diff --git a/dynamic/palm_dynamic_aerosol.py b/dynamic/palm_dynamic_aerosol.py new file mode 100644 index 0000000000000000000000000000000000000000..2995ccf331f5d58121ea00c6e94efd45be0a05d1 --- /dev/null +++ b/dynamic/palm_dynamic_aerosol.py @@ -0,0 +1,193 @@ +#------------------------------------------------------------------------------ +# +# Script for processing of WRF-CHEM files to PALM dynamic driver. +# +#------------------------------------------------------------------------------ +import os +import netCDF4 +import numpy as np +from palm_dynamic_config import * + +def define_bins(nbin, reglim): + # degine Dmid based on nbin and reglim + nbins = np.sum(nbin) # = subrange 1 + subrange 2 + # Log-normal to sectional + vlolim = np.zeros(nbins) + vhilim = np.zeros(nbins) + dmid = np.zeros(nbins) + bin_limits = np.zeros(nbins) + # Sectional bin limits + ratio_d = reglim[1] / reglim[0] + for b in range(nbin[0]): + vlolim[b] = np.pi / 6.0 * (reglim[0] * ratio_d **(float(b) / nbin[0]))**3 + vhilim[b] = np.pi / 6.0 * (reglim[0] * ratio_d **(float(b+1) / nbin[0]))**3 + dmid[b] = np.sqrt((6.0 * vhilim[b] / np.pi )**0.33333333 * (6.0 * vlolim[b] / np.pi)**0.33333333) + ratio_d = reglim[2] / reglim[1] + + for b in np.arange(nbin[0], np.sum(nbin),1): + c = b-nbin[0] + vlolim[b] = np.pi / 6.0 * (reglim[1] * ratio_d **(float(c) / nbin[1]))**3 + vhilim[b] = np.pi / 6.0 * (reglim[1] * ratio_d **(float(c+1) / nbin[1])) ** 3 + dmid[b] = np.sqrt((6.0 * vhilim[b] / np.pi)**0.33333333 * ( 6.0 * vlolim[b] / np.pi)**0.33333333) + bin_limits = (6.0 * vlolim / np.pi )**0.33333333 + bin_limits = np.append(bin_limits, reglim[-1]) + return dmid, bin_limits + +def translate_aerosol_species(name): + testname = name + translation_table = { + "SO4": "so4", + "NO" : "no3", + "NH" : "nh4", + "BC" : "bc", + "OC" : "smpa,smpbb,glysoa_sfc,biog1_c,biog1_o", + "SS" : "cl,na", + "DU" : "co3,ca,oin" + } + if testname in translation_table: + return translation_table[testname] + else: + return None + +def range_overlap(range1, range2): + x1, x2 = range1.start, range1.stop + y1, y2 = range2.start, range2.stop + return x1 <= y2 and y1 <= x2 + +def aerosol_binoverlap(palm_binlim, wrfchem_binlim): + overlap_ratio = np.zeros( (len(palm_binlim)-1, len(wrfchem_binlim)-1) ) + aerobin_open = [] + for pbin in range(0,len(palm_binlim)-1): + palm_range = range(int(palm_binlim[pbin]* 1e+9), int(palm_binlim[pbin+1]* 1e+9)) + for wbin in range(0, len(wrfchem_binlim)-1): + wrfchem_range = range(int(wrfchem_binlim[wbin]* 1e+9)+1, int(wrfchem_binlim[wbin+1]* 1e+9)) + + if range_overlap(palm_range, wrfchem_range): + aerobin_open.append('_a0'+ str(wbin+1)) + overlap = len(set(palm_range) & set(wrfchem_range)) + overlap_ratio[pbin, wbin] = overlap/len(wrfchem_range) + return aerobin_open, overlap_ratio + +def aerosolProfile(interp_files, listspec, nbin, reglim, wrfchem_bin_limits): + infile = netCDF4.Dataset(interp_files[0], "r", format="NETCDF4") + alt = infile.variables['init_atmosphere_alt'][0] + u = infile.variables['init_atmosphere_u'] + v = infile.variables['init_atmosphere_v'] + z_level = infile.variables['z'] + # upwind location & mass frac + aero_massfrac_a = np.zeros((z_level.size, len(listspec))) + for zlev in range(0, z_level.size): + u_wnd = u[0,zlev,:,:] + v_wnd = v[0,zlev,:,:] + wnd_dir = np.mod(180 + np.rad2deg(np.arctan2(u_wnd, v_wnd)),360) + wnd_avg = np.mean(wnd_dir) + if 0 < wnd_avg <= 45: + prf_x = round(wnd_dir.shape[0]/2) + prf_y = 0 + elif 315 < wnd_avg <=0: + prf_x = round(wnd_dir.shape[0]/2) + prf_y = 0 + elif 45 < wnd_avg <= 135: + prf_x = round(wnd_dir.shape[0]/2) + prf_y = round(wnd_dir.shape[1]/2) + elif 135 < wnd_avg <= 225: + prf_x = round(wnd_dir.shape[0]/2) + prf_y = wnd_dir.shape[1] + else: + prf_x = round(wnd_dir.shape[0]/2) + prf_y = round(wnd_dir.shape[1]/2) + # mass at each z-level (z,composition_index) + aero_mass = np.zeros((len(listspec))) + for naero in range(0, len(listspec)): + spec_mass = infile.variables['init_atmosphere_'+ listspec[naero]][0] + aero_mass[naero] = spec_mass[zlev, prf_y, prf_x] + total_mass = np.sum(aero_mass) + # mass fraction values + for naero in range(0, len(listspec)): + aero_massfrac_a[zlev, naero] = aero_mass[naero]/total_mass + # aerosol concen# + # ug/kg-dryair to aerosol# concen.(# m-3):inverse density (m3 kg-1) + # define palm bins + bin_dmid, bin_lims = define_bins(nbin, reglim) + # define overlap of palm & wrfchem bins + open_bins, overlap_ratio = aerosol_binoverlap(bin_lims, wrfchem_bin_limits) + open_bins = sorted(set(open_bins), key=open_bins.index) + # convert units, weight Dmid by size-bin use overlap ratio + aero_con = np.zeros((z_level.size, bin_dmid.size)) + for n_dmid in range(0, bin_dmid.size): + outval = 0 + for abin in range(0, len(open_bins)): + inval = infile.variables['init_aerosol'+ open_bins[abin]][0] + inval = inval[zlev, prf_y, prf_x] * alt[zlev, prf_y, prf_x] + inval = inval * overlap_ratio[n_dmid, abin] + outval = outval + inval + aero_con[zlev, n_dmid] = outval + infile.close() + return aero_massfrac_a, aero_con + +def aerosolMassWrfchemBoundary(dimensions, _dim1, _dim2, interp_files, listspec, side): + _dim1_size = dimensions[_dim1+'dim'] + _dim2_size = dimensions[_dim2+'dim'] + # mass_frac_a: val_side[time, z, y, composition_index] + val_side = np.zeros( (len(interp_files), _dim1_size, _dim2_size, len(listspec)) ) + for ts in range(0, len(interp_files)): + infile = netCDF4.Dataset(interp_files[ts], "r", format="NETCDF4") + val_spec = np.zeros( (_dim1_size, _dim2_size, len(listspec)) ) + for n_spec in range(0, len(listspec)): + if (side == 'left'): + val_spec[:,:,n_spec] = infile.variables['init_atmosphere_'+ listspec[n_spec]][0, :, :, 0] + elif (side == 'right'): + val_spec[:,:,n_spec] = infile.variables['init_atmosphere_'+ listspec[n_spec]][0, :, :, dimensions['xdim'] - 1] + elif (side =='south'): + val_spec[:,:,n_spec] = infile.variables['init_atmosphere_'+ listspec[n_spec]][0, :, 0, :] + elif (side == 'north'): + val_spec[:,:,n_spec] = infile.variables['init_atmosphere_'+ listspec[n_spec]][0, :, dimensions['ydim'] - 1, :] + elif (side == 'top'): + val_spec[:,:,n_spec] = infile.variables['init_atmosphere_'+ listspec[n_spec]][0, dimensions['zdim']-1, :, :] + val_sum = np.sum(val_spec,axis = 2) + # write mass frac values + for n_spec in range(0, len(listspec)): + val_side[ts, :, :, n_spec] = val_spec[:,:,n_spec]/val_sum + infile.close() + return val_side + +def aerosolConWrfchemBoundary(dimensions, _dim1, _dim2, interp_files, side, nbin, reglim, wrfchem_bin_limits): + # define palm bins + bin_dmid, bin_lims = define_bins(nbin, reglim) + # define overlap of palm & wrfchem bins + open_bins, overlap_ratio = aerosol_binoverlap(bin_lims, wrfchem_bin_limits) + open_bins = sorted(set(open_bins), key=open_bins.index) + # adjust dimensions + _dim1_size = dimensions[_dim1+'dim'] + _dim2_size = dimensions[_dim2+'dim'] + # aerosol concentration number [time, z, y, Dmid] + aero_conc = np.zeros((len(interp_files), _dim1_size, _dim2_size, bin_dmid.size)) + for ts in range(0, len(interp_files)): + infile = netCDF4.Dataset(interp_files[ts], "r", format="NETCDF4") + outval = np.zeros( (_dim1_size, _dim2_size) ) + for n_dmid in range(0, bin_dmid.size): + + for abin in range(0, len(open_bins)): + if (side == 'left'): + alt = infile.variables['init_atmosphere_alt'][0, :, :, 0] + inval = infile.variables['init_aerosol'+ open_bins[abin]][0, :, :, 0] + elif (side == 'right'): + alt = infile.variables['init_atmosphere_alt'][0, :, :, dimensions['xdim'] - 1] + inval = infile.variables['init_aerosol'+ open_bins[abin]][0, :, :, dimensions['xdim'] - 1] + elif (side =='south'): + alt = infile.variables['init_atmosphere_alt'][0, :, 0, :] + inval = infile.variables['init_aerosol'+ open_bins[abin]][0, :, 0, :] + elif (side == 'north'): + alt = infile.variables['init_atmosphere_alt'][0, :, dimensions['ydim'] - 1, :] + inval = infile.variables['init_aerosol'+ open_bins[abin]][0, :, dimensions['ydim'] - 1, :] + elif (side == 'top'): + alt = infile.variables['init_atmosphere_alt'][0, dimensions['zdim']-1, :, :] + inval = infile.variables['init_aerosol'+ open_bins[abin]][0, dimensions['zdim']-1, :, :] + # convert and factor for size-bin + inval = inval * alt + inval = inval * overlap_ratio[n_dmid, abin] + outval = outval + inval + + aero_conc[ts, :, :, n_dmid] = outval + infile.close() + return aero_conc diff --git a/dynamic/palm_dynamic_output.py b/dynamic/palm_dynamic_output.py index c9fdfbe48feedae643dd87b554a44d042dd4eb17..8c088185c0572a787d4eda164e94cdb34db4307c 100644 --- a/dynamic/palm_dynamic_output.py +++ b/dynamic/palm_dynamic_output.py @@ -15,12 +15,12 @@ import numpy as np import netCDF4 from palm_wrf_utils import palm_wrf_gw from palm_dynamic_config import * +import palm_dynamic_aerosol def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, dimensions, z_levels, z_levels_stag, ztop, z_soil_levels, dx, dy, lon_center, lat_center, rad_times_proc, rad_values_proc, sma, nested_domain): - print('\nProcessing interpolated files to dynamic driver') # dimension of the time coordinate dimtimes = len(times_sec) @@ -39,7 +39,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, boundary = ['left','right','south', 'north','top'] # atmospheric variables atmos_var = ['pt','qv','u','v','w','soil_m','soil_t'] - all_variables = atmos_var + wrfchem_spec + dynam_chem_variables = atmos_var + wrfchem_spec # prepare influx/outflux area sizes zstag_all = np.r_[0., z_levels_stag, ztop] @@ -73,17 +73,83 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, _val_y[:] = y[:] _val_x = outfile.createVariable('x',"f4", ("x")) _val_x[:] = x[:] - - # create all other variables in outout file - def add_interpDim(all_variables): + #--------------------------------------------------------------------------- + # include aerosols, add dimensions and variables + if aerosol_wrfchem: + # upwind values for aerosol concen and mass frac + aero_massfrac_a, aero_con = palm_dynamic_aerosol.aerosolProfile(interp_files, listspec, nbin, reglim, wrfchem_bin_limits) + + # create dimensions + outfile.createDimension('Dmid', sum(nbin)) + outfile.createDimension('composition_index', len(listspec)) + outfile.createDimension('max_string_length', 25) + # 2D vertical profile variables + _val_composition_name = outfile.createVariable('composition_name',"S1", ("composition_index", "max_string_length")) + _val_composition_name.setncattr('long_name', "aerosol composition name") + composition_name = np.array([listspec],dtype = 'S25') + _val_composition_name[:] = netCDF4.stringtochar(composition_name) + + _val_aerosol_mass_a = outfile.createVariable('init_atmosphere_mass_fracs_a', "f4", ("z", "composition_index"),fill_value=fillvalue_float) + _val_aerosol_mass_a.setncattr('long_name',"initial mass fraction profile: a bins") + _val_aerosol_mass_a[:,:] = aero_massfrac_a + + _val_aerosol_mass_b = outfile.createVariable('init_atmosphere_mass_fracs_b', "f4", ("z", "composition_index"),fill_value=fillvalue_float) + _val_aerosol_mass_b.setncattr('long_name',"initial mass fraction profile: b bins") + if nf2a == 1.0: + _val_aerosol_mass_b[:,:] = 0 + elif nf2a <= 1.0: + _val_aerosol_mass_b[:,:] = 1 - aero_massfrac_a + + _val_aerosol_con = outfile.createVariable('init_atmosphere_aerosol', "f4", ("z","Dmid"),fill_value=fillvalue_float) + _val_aerosol_con.setncattr('lod', 1) + _val_aerosol_con.setncattr('long_name',"initial vertical profile of aerosol concentration") + _val_aerosol_con.setncattr('units',"#/m3") + _val_aerosol_con = aero_con + + # time dependent variables - aerosols + for side in boundary: + if (side == 'left' or side == 'right'): + _dim1 = "z" + _dim2 = "y" + elif (side == 'south' or side == 'north'): + _dim1 = "z" + _dim2 = "x" + elif side =='top': + _dim1 = "y" + _dim2 = "x" + # mass fractions + _var_a = outfile.createVariable('ls_forcing_'+ side + '_mass_fracs_a', "f4", ("time", _dim1, _dim2, "composition_index"), + fill_value=fillvalue_float) + _var_a.setncattr('lod', 1) + _var_a.setncattr('long_name', "boundary conditions of mass fraction profile: a bins") + _var_a = palm_dynamic_aerosol.aerosolMassWrfchemBoundary(dimensions, _dim1, _dim2, interp_files, listspec, side) + + _var_b = outfile.createVariable('ls_forcing_'+ side + '_mass_fracs_b', "f4", ("time", _dim1, _dim2, "composition_index"), + fill_value=fillvalue_float) + _var_b.setncattr('lod', 1) + _var_b.setncattr('long_name', "boundary conditions of mass fraction profile: b bins") + if (nf2a >= 1.0): + _var_b = 0 + elif (nf2a < 1.0): + _var_b = 1-_var_a + # aerosol concen.# + _var_con = outfile.createVariable('ls_forcing_'+ side + '_aerosol', "f4", ("time", _dim1, _dim2, "Dmid"), + fill_value=fillvalue_float) + _var_con.setncattr('lod', 1) + _var_con.setncattr('long_name',"boundary condition of aerosol concentration") + _var_con = palm_dynamic_aerosol.aerosolConWrfchemBoundary(dimensions, _dim1, _dim2, interp_files, side, nbin, reglim, wrfchem_bin_limits) + + #--------------------------------------------------------------------------- + # create dynamical & chemical variables in output + def add_interpDim(dynam_chem_variables): print('Adding dimensions to Dynamic Driver') # surface pressure _val_surface_forcing_surface_pressure = outfile.createVariable('surface_forcing_surface_pressure', "f4",("time")) # geostrophic wind _val_ls_forcing_ug = outfile.createVariable('ls_forcing_ug', "f4", ("time", "z"),fill_value=fillvalue_float) _val_ls_forcing_vg = outfile.createVariable('ls_forcing_vg', "f4", ("time", "z"),fill_value=fillvalue_float) - # create variables in outfile for all variables - for var in all_variables: + # create variables other dynamic & chemical variables + for var in dynam_chem_variables: if (var == 'pt' or var == 'qv'): _val_init_var = outfile.createVariable('init_atmosphere_'+ var, "f4", ("z", "y", "x"),fill_value=fillvalue_float) elif (var == 'soil_m' or var == 'soil_t'): @@ -159,16 +225,16 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, else: _val_ls_forcing.setncattr('units', 'ppm') outfile.close() - # create all other variables in outout file - add_interpDim(all_variables) + # function - create dynamical & chemical variables + add_interpDim(dynam_chem_variables) - # read interpolated files and write values to outfile - def add_interpValues(all_variables): - print('\nAdding initializing variable values to Dynamic Driver') + # read interpolated files and write values for dynamical & chemical variables + def add_interpValues(dynam_chem_variables): + print('Adding initializing variable values to Dynamic Driver') infile = netCDF4.Dataset(interp_files[0], "r", format="NETCDF4") outfile = netCDF4.Dataset(dynamic_driver_file, "r+", format="NETCDF4") # initialization variables - for var in all_variables: + for var in dynam_chem_variables: if (var == 'soil_m' or var == 'soil_t'): init_var = infile.variables['init_'+var] _val_init_var = outfile.variables['init_'+var] @@ -190,7 +256,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, infile.close() outfile.close() - # time dependent variables (all time steps) + # time dependent variables - dynamical & chemical variables if not nested_domain: print('\nAdding time dependent variables value to Dynamic Driver') outfile = netCDF4.Dataset(dynamic_driver_file, "r+", format="NETCDF4") @@ -211,9 +277,8 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, surface_forcing_surface_pressure = infile.variables['surface_forcing_surface_pressure'] _val_surface_forcing_surface_pressure = outfile.variables['surface_forcing_surface_pressure'] _val_surface_forcing_surface_pressure[ts] = np.average(surface_forcing_surface_pressure[:,:,:], axis = (1,2))[0] - # all other variables except soil_t, soli_m, u, v, w - #['pt', 'qv', 'u', 'v', 'w', 'soil_m', 'soil_t', 'no', 'no2', 'o3', 'PM10', 'PM2_5_DRY'] - for var in all_variables: + # other dynamical & chemical variables except: soil_t, soli_m, u, v, w + for var in dynam_chem_variables: if (var == 'soil_m' or var == 'soil_t' or var == 'u' or var == 'v' or var == 'w'): continue else: @@ -231,7 +296,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, elif side == 'top': _val_ls_forcing_var[ts, :, :] = init_var[0, dimensions['zdim'] - 1, :, :] - # variables:u, v, w - mass balancing + # dynamical variables:u, v, w - mass balancing init_atmosphere_u = infile.variables['init_atmosphere_u'] init_atmosphere_v = infile.variables['init_atmosphere_v'] init_atmosphere_w = infile.variables['init_atmosphere_w'] @@ -292,9 +357,8 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, infile.close() outfile.close() - # read interpolated files and write values to outfile - add_interpValues(all_variables) + add_interpValues(dynam_chem_variables) # add radiation input if len(rad_times_proc) > 0: @@ -318,17 +382,16 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec, var[:] = rad_values_proc[2][:] outfile.close() - # correct chemical species names to PALM output outfile = netCDF4.Dataset(dynamic_driver_file, "a", format="NETCDF4") - if 'PM2_5_DRY' in all_variables: + if 'PM2_5_DRY' in dynam_chem_variables: outfile.renameVariable('init_atmosphere_PM2_5_DRY','init_atmosphere_pm2_5') outfile.renameVariable('ls_forcing_left_PM2_5_DRY', 'ls_forcing_left_pm2_5') outfile.renameVariable('ls_forcing_right_PM2_5_DRY','ls_forcing_right_pm2_5') outfile.renameVariable('ls_forcing_south_PM2_5_DRY','ls_forcing_south_pm2_5') outfile.renameVariable('ls_forcing_north_PM2_5_DRY','ls_forcing_north_pm2_5') outfile.renameVariable('ls_forcing_top_PM2_5_DRY','ls_forcing_top_pm2_5') - if 'PM10' in all_variables: + if 'PM10' in dynam_chem_variables: outfile.renameVariable('init_atmosphere_PM10','init_atmosphere_pm10') outfile.renameVariable('ls_forcing_left_PM10', 'ls_forcing_left_pm10') outfile.renameVariable('ls_forcing_right_PM10','ls_forcing_right_pm10') diff --git a/dynamic/palm_wrf_utils.py b/dynamic/palm_wrf_utils.py index 96c95c71e3947e3254f905cea36d016e772eb3ab..fb3f6569e0f7fb1231c2f62fd0ea1d7968ecfd92 100644 --- a/dynamic/palm_wrf_utils.py +++ b/dynamic/palm_wrf_utils.py @@ -20,6 +20,7 @@ from metpy.interpolate import interpolate_1d, log_interpolate_1d from metpy.units import units from palm_dynamic_config import * +import palm_dynamic_aerosol # Constants directly equivalent to WRF code radius = 6370000.0 @@ -29,13 +30,16 @@ rd_cp = 2./7. #from WRF v4 technote (R_d / c_p) wrf_base_temp = 300. #NOT wrfout T00 _ax = np.newaxis +if aerosol_wrfchem: + # Calculate aerosol bin overlaps + open_bin, overlap_ratio = palm_dynamic_aerosol.aerosol_binoverlap(reglim, wrfchem_bin_limits) + open_bin = sorted(set(open_bin), key=open_bin.index) class WRFCoordTransform(object): 'Coordinate transformer for WRFOUT files' def __init__(self, ncf): attr = lambda a: getattr(ncf, a) - # Define grids # see http://www.pkrc.net/wrf-lambert.html #latlon_wgs84 = pyproj.Proj(proj='latlong', @@ -215,7 +219,7 @@ def calc_gp(f, ph): return np.array(gp) def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag, - z_soil_levels, origin_z, terrain, wrf_hybrid_levs, vinterp_terrain_smoothing): + z_soil_levels, origin_z, terrain, wrf_hybrid_levs, vinterp_terrain_smoothing,nz,ny,nx): zdim = len(z_levels) zwdim = len(z_levels_stag) @@ -279,7 +283,7 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag, # Report gpdelta = gpf2 - gpf - print('GP deltas by level:') + #print('GP deltas by level:') # Changed here #for k in range(gpf.shape[0]): #print_dstat(k, gpdelta[k]) @@ -363,9 +367,7 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag, # zsoil is taken from wrf - not need to define it # ======================== CHEMISTRY =================================================== - # convert wrf-chem units to PALM units # convert ppmv to ppm for PALM except PM10 & PM2_5_DRY: micrograms m-3 to kg/m3 - def chem_from_wrfchem(wrfchem_spec): for spec in wrfchem_spec: chem_data_raw = nc_infile.variables[spec][0] @@ -377,11 +379,65 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag, vdata = nc_outfile.createVariable('init_atmosphere_'+spec,"f4",("Time", "z","south_north","west_east")) vdata[0,:,:,:] = chem_data - nc_infile.close() - nc_wrf.close() - nc_outfile.close() - chem_from_wrfchem(wrfchem_spec) + # ======================== AEROSOLS ===================================== + # aerosol_mass_fraction + def aerosol_mass_wrfchem(listspec,open_bin, overlap_ratio): + var_aero = palm_dynamic_aerosol.translate_aerosol_species(listspec[0]).split(",") + var_size = nc_infile.variables[var_aero[0]+'_a01'].shape[0] + # by species & sub-species + for spec in listspec: + val_specout = np.zeros((var_size)) + spec_wrf = palm_dynamic_aerosol.translate_aerosol_species(spec).split(",") + # by aerosol-size bin + naero = 0 + for aero_bin in open_bin: + val_specint = np.zeros((var_size)) + # interpolate sub-species & sum + for nsp in spec_wrf: + inval_spec = nc_infile.variables[nsp + aero_bin][0] + inval_spec = np.r_[inval_spec[0:1], inval_spec] + inval_spec = interpolate_1d(z_levels, height, inval_spec) + val_specint = val_specint + inval_spec + # factor for size-bin size + val_specint = val_specint*(sum(overlap_ratio[:,naero])) + val_specout = val_specout + val_specint + naero = naero + 1 + # write interp(sum) to file + vdata = nc_outfile.createVariable('init_atmosphere_'+spec,"f4",("Time", "z","south_north","west_east")) + vdata[0,:,:,:] = val_specout + + def aerosol_concen_wrfchem(listspec): + alt_data_raw = nc_infile.variables['ALT'][0] + alt_data_raw = np.r_[alt_data_raw[0:1], alt_data_raw] + alt_data = interpolate_1d(z_levels, height, alt_data_raw) + vdata = nc_outfile.createVariable('init_atmosphere_alt', "f4",("Time", "z", "south_north", "west_east")) + vdata[0,:,:,:] = alt_data + # by size bin + for ab in range(1,5): + # by species & sub-species + spec_wrf = [] + for spec in listspec: + _wrf = palm_dynamic_aerosol.translate_aerosol_species(spec).split(",") + spec_wrf.extend(_wrf) + # interpolate & sum + val_specout = np.zeros((nz,ny,nx)) + for wspec in spec_wrf: + inval_spec = nc_infile.variables[wspec +'_a0'+ str(ab)][0] + inval_spec = np.r_[inval_spec[0:1], inval_spec] + inval_specint = interpolate_1d(z_levels, height, inval_spec) + val_specout = val_specout + inval_specint + # write to .interp file + vdata = nc_outfile.createVariable('init_aerosol_a0'+ str(ab), "f4",("Time", "z", "south_north", "west_east")) + vdata[0,:,:,:] = val_specout + # include aerosols + if aerosol_wrfchem: + aerosol_mass_wrfchem(listspec, open_bin, overlap_ratio) + aerosol_concen_wrfchem(listspec) + + nc_infile.close() + nc_wrf.close() + nc_outfile.close() def palm_wrf_gw(f, lon, lat, levels, tidx=0): '''Calculate geostrophic wind from WRF using metpy'''