Merge branch 'aerosol_profile' into 'main'

fixed aerosol profile

See merge request mmbees-git/wrf_chem_for_palm!22

dynamic/palm_dynamic.py

@@ -172,7 +172,6 @@ if grid_from_static:
     # close static driver nc file
     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)
@@ -254,15 +253,14 @@ if nested_domain:
     end_time = start_time
 
 # get wrf times and sort wrf files by time
-print('\nAnalyse WRF files dates:')
+#print('\nAnalyse WRF files dates:')
 file_times = []
 for wrf_file in wrf_file_list:
-    #print('\nWRF file',wrf_file)
     nc_wrf = netCDF4.Dataset(wrf_file, "r", format="NETCDF4")
     ta = nc_wrf.variables['Times'][:]
     t = ta.tobytes().decode("utf-8")
     td = datetime.strptime(t, '%Y-%m-%d_%H:%M:%S')
-    print(os.path.basename(wrf_file), ': ', td)
+    #print(os.path.basename(wrf_file), ': ', td)
     file_times.append((td,wrf_file))
     nc_wrf.close()
 
@@ -273,7 +271,7 @@ for tf in file_times:
     if end_time is None or tf[0] <= end_time:
         times.append(tf[0])
         wrf_files.append(tf[1])
-#print('\n'.join('{}'.format(t) for t in times))
+
 print('PALM output times:', ', '.join('{}'.format(t) for t in times))
 
 if not times.__contains__(start_time):
@@ -364,11 +362,8 @@ for wrf_file in wrf_files_proc:
                 # add chemical species if included
                 if len(wrfchem_spec)>0:
                     wrfchem_variables = wrfchem_dynamic + wrfchem_spec
-                    #wrfchem_variables.append('ALT')        # inverse density
 
                 # add aerosol species if included
-                #if aerosol_wrfchem:
-                #    wrfchem_variables = wrfchem_variables + wrfchem_aerosols
                 N_avr = 6.022e23 # Avargardo constant
                 inv_den = f_wrf.variables['ALT'][0] # inverse density
 
@@ -471,7 +466,6 @@ for wrf_file in wrf_files_proc:
                     v_out = f_out.createVariable('aerosol'+ aero_bin, 'f4', nbn.dimensions)
                     v_out[:] = regridder.regrid(nbn_val[...,regridder.ys,regridder.xs])
 
-
                 # U and V have special treatment (unstaggering)
                 v_out = f_out.createVariable('U', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east'))
                 v_out[:] = regridder_u.regrid(f_wrf.variables['U'][...,regridder_u.ys,regridder_u.xs])
@@ -597,4 +591,4 @@ palm_dynamic_output(wrf_files_proc, interp_files, dynamic_driver_file, times_sec
                     z_levels, z_levels_stag, ztop, z_soil_levels, dx, dy, cent_lon, cent_lat,
                     rad_times_proc, rad_values_proc, soil_moisture_adjust, nested_domain)
 
-print('Creation of dynamic driver finished.')
+print('\nCreation of dynamic driver finished.')

dynamic/palm_dynamic_aerosol.py

@@ -31,21 +31,18 @@ def upwind_location(zlev, u, v):
     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)
+    if  0 < wnd_avg <= 45 or 315 < wnd_avg <=0:
         prf_y = 0
+        prf_x = round(wnd_dir.shape[1]/2)
     elif 45 < wnd_avg <= 135:
-        prf_x = round(wnd_dir.shape[0]/2)
-        prf_y = round(wnd_dir.shape[1]/2)
+        prf_y = round(wnd_dir.shape[0]/2)
+        prf_x = round(wnd_dir.shape[1])
     elif 135 < wnd_avg <= 225:
-        prf_x = round(wnd_dir.shape[0]/2)
-        prf_y = wnd_dir.shape[1]
+        prf_y = wnd_dir.shape[0]
+        prf_x = round(wnd_dir.shape[1]/2)
     else:
-        prf_x = round(wnd_dir.shape[0]/2)
-        prf_y = round(wnd_dir.shape[1]/2)
+        prf_y = round(wnd_dir.shape[0]/2)
+        prf_x = 0
 
     return prf_x, prf_y
 

dynamic/palm_dynamic_output.py

@@ -84,7 +84,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
     #---------------------------------------------------------------------------
     # include aerosols, add dimensions and variables
     if aerosol_wrfchem:
-        print('Adding aerosol variables')
+        print('\tCreating aerosol variables in dynamic driver')
         # create dimensions
         outfile.createDimension('Dmid', sum(nbin))
 
@@ -169,7 +169,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
     #---------------------------------------------------------------------------
     # create dynamical & chemical variables in output
     def add_interpDim(dynam_chem_variables):
-        print('Adding dynamical variables')
+        print('\tCreating dynamical variables in dynamic driver')
         # surface pressure
         _val_surface_forcing_surface_pressure = outfile.createVariable('surface_forcing_surface_pressure', "f4",("time"))
         # geostrophic wind
@@ -270,7 +270,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
     
     # read interpolated files and write values for dynamical & chemical variables
     def add_interpValues(dynam_chem_variables):
-        print('Adding initializing variables to Dynamic Driver')
+        print('\tCreating initialisation variablesin dynamic driver')
         infile = netCDF4.Dataset(interp_files[0], "r", format="NETCDF4")
         outfile = netCDF4.Dataset(dynamic_driver_file, "r+", format="NETCDF4")
         # initialization variables
@@ -301,11 +301,11 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
 
         # time dependent variables - dynamical & chemical variables
         if not nested_domain:
-            print('\nAdding time dependent variables to Dynamic Driver')
+            print('\tCreating time dependent variables in dynamic driver')
             outfile = netCDF4.Dataset(dynamic_driver_file, "r+", format="NETCDF4")
             for ts in range(0, len(interp_files)):
                 # geostrophic wind
-                print('Open wrf file: '+wrf_files[ts])
+                #print('Open wrf file: '+wrf_files[ts])
                 nc_wrf = netCDF4.Dataset(wrf_files[ts], 'r')
                 ug, vg = palm_wrf_gw(nc_wrf, lon_center, lat_center, z_levels)
                 _val_ls_forcing_ug = outfile.variables['ls_forcing_ug']
@@ -314,7 +314,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
                 _val_ls_forcing_vg = vg
                 nc_wrf.close()
 
-                print("Processing interpolated file: ",interp_files[ts])
+                print("\tProcessing interpolated file: ",interp_files[ts])
                 infile = netCDF4.Dataset(interp_files[ts], "r", format="NETCDF4")
                 # surface pressure
                 surface_forcing_surface_pressure = infile.variables['surface_forcing_surface_pressure']
@@ -358,7 +358,7 @@ def palm_dynamic_output(wrf_files, interp_files, dynamic_driver_file, times_sec,
                         - (vynorth * areas_yb).sum()
                         - (wztop * areas_zb).sum())
                 mass_corr_v = mass_disbalance / area_boundaries
-                print('Mass disbalance: {0:8g} m3/s (avg = {1:8g} m/s)'.format(mass_disbalance, mass_corr_v))
+                #print('Mass disbalance: {0:8g} m3/s (avg = {1:8g} m/s)'.format(mass_disbalance, mass_corr_v))
                 uxleft  -= mass_corr_v
                 uxright += mass_corr_v
                 vysouth -= mass_corr_v

dynamic/palm_wrf_utils.py

@@ -255,7 +255,6 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag,
         target_terrain = ndimage.gaussian_filter(terrain, sigma=vinterp_terrain_smoothing, order=0)
     print('Morphing WRF terrain ({0} ~ {1}) to PALM terrain ({2} ~ {3})'.format(
         wrfterr.min(), wrfterr.max(), target_terrain.min(), target_terrain.max()))
-    #print_dstat('terrain shift', wrfterr - target_terrain[:,:])
 
     # Load original dry air column pressure
     mu = nc_infile.variables['MUB'][0,:,:] + nc_infile.variables['MU'][0,:,:]
@@ -283,9 +282,6 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag,
 
     # Report
     gpdelta = gpf2 - gpf
-    #print('GP deltas by level:')
-    #for k in range(gpf.shape[0]):
-    #    print_dstat(k, gpdelta[k])
 
     # Because we require levels below the lowest level from WRF, we will always
     # add one layer at zero level with repeated values from the lowest level.
@@ -327,7 +323,6 @@ def palm_wrf_vertical_interp(infile, outfile, wrffile, z_levels, z_levels_stag,
     init_atmosphere_v = interpolate_1d(z_levels, height, v_raw)
     
     vdata = nc_outfile.createVariable('init_atmosphere_v', "f4", ("Time", "z","south_north","west_east"))
-    #vdata.coordinates = "XLONG_V XLAT_V XTIME"
     vdata[0,:,:,:] = init_atmosphere_v
     
     w_raw = nc_infile.variables['W'][0]
@@ -525,7 +520,6 @@ def calcgw_gfs(v, lat, lon):
     j = np.searchsorted(lons[0,:], lon)
     if abs(lons[0,i+1] - lon) < abs(lons[0,i] - lon):
         j = j+1
-    #print('level', v.level, 'height', height[i,j], lats[i,j], lons[i,j])
 
     # Set up some constants based on our projection, including the Coriolis
     # parameter and grid spacing, converting lon/lat spacing to Cartesian
@@ -587,15 +581,11 @@ if __name__ == '__main__':
         phf, phh = calc_ph_sigma(f, mu)
         gp_calc = calc_gp(f, phf)
         delta = gp_calc - gp
-        #for lev in range(delta.shape[0]):
-        #    print_dstat(lev, delta[lev])
 
         print('\nUsing hybrid:')
         phf, phh = calc_ph_hybrid(f, mu)
         gp_calc = calc_gp(f, phf)
         delta = gp_calc - gp
-        #for lev in range(delta.shape[0]):
-        #    print_dstat(lev, delta[lev])
 
         f.close()
     if args.camx: