Read native DYNAMICO data

metrics/half_plume.py

@@ -61,35 +61,19 @@ def max_volume_plot(variable_chim, variable_dyn, day_list):
 def to_radian(angle):
     return angle*(np.pi/180)
 
-def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
-    hour = day*24+hour # assume output every hour, first day is day 0  
+def read_DYNAMICO_common(hybrids, hour, area, ps, pmid, phi, q):
     radius, g, p0 = 6.371e6, 9.81, 1e5 # Earth radius, gravity, reference pressure for hybrid coordinate
-
-#    print('Reading DYNAMICO data from %s ...'%dynamico_nc)
-
-    ds = nc.Dataset(dynamico_nc)
-    lon, lat = ds['lon'][:], ds['lat'][:]
-    nlon, nlat = len(lon), len(lat)
-    q   = ds['q'][hour, 0, :, :, :]
-    ps  = ds['ps'][hour, :, :]
-    phi = ds['PHI'][hour, 1:, :, :] # geopotential'
-#    print('shape(phi) = ', phi.shape)
-#    print("Max of DYNAMICO q at hour=%d : %f"%(hour, q.max()))
-
-    lon_deg, lat_deg = np.meshgrid(lon, lat)          # degrees
-    lon, lat = to_radian(lon_deg), to_radian(lat_deg) # radians
-    area = (radius**2)*np.cos(lat)*(2*np.pi/nlon)*(np.pi/nlat)
-
     volume = compute_volume(area/g, phi) # m^3
 
 #    print('  geopotential last hour = ', phi[:,12,27])
 
-    ds1 = nc.Dataset(dynamico_hybrids)
+    ds1 = nc.Dataset(hybrids)
     ap, bp = ds1['hyai'][:], ds1['hybi'][:]    
     lev = len(ap) - 1
-    p = np.zeros((lev, nlat, nlon)) # lat changes with row, lon with column
+    nlat, nlon = ps.shape
+    p = np.zeros((lev, nlat, nlon))
 
-    if check_compute_pressure(ds1['hyam'][:], ds1['hybm'][:], p0, ps, ds['P'][hour,:,:,:] ):
+    if check_compute_pressure(ds1['hyam'][:], ds1['hybm'][:], p0, ps, pmid ):
         print("  Discrepancy in computed pressure. Stop.")
         return True
 
@@ -99,8 +83,56 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
 
 #    print("Total SO2 mass in DYNAMICO = ", (q*air_mass).sum())
 #    print("Max SO2 mixing ratio (kg/kg) in DYNAMICO = ", q.max())
+    return radius, volume, SO2_density
+
+def read_DYNAMICO(dynamico_nc, hybrids, day=7, hour=22):
+    radius, g, p0 = 6.371e6, 9.81, 1e5 # Earth radius, gravity, reference pressure for hybrid coordinate
+    hour = day*24+hour # assume output every hour, first day is day 0  
+
+#    print('Reading DYNAMICO data from %s ...'%dynamico_nc)
+
+    ds = nc.Dataset(dynamico_nc)
+    lon, lat = ds['lon'][:], ds['lat'][:]
+    nlon, nlat = len(lon), len(lat)
+    q    = ds['q'][hour, 0, :, :, :]
+    ps   = ds['ps'][hour, :, :]
+    pmid = ds['P'][hour,:,:,:]       # pressure at layers, used to check hybrid pressure formula
+    phi  = ds['PHI'][hour, 1:, :, :] # geopotential'
+#    print('shape(phi) = ', phi.shape)
+#    print("Max of DYNAMICO q at hour=%d : %f"%(hour, q.max()))
+
+    lon_deg, lat_deg = np.meshgrid(lon, lat)          # degrees
+    lon, lat = to_radian(lon_deg), to_radian(lat_deg) # radians
+    area = (radius**2)*np.cos(lat)*(2*np.pi/nlon)*(np.pi/nlat)
+
+    radius, volume, SO2_density = read_DYNAMICO_common(hybrids, hour, area, ps, pmid, phi, q)
     return radius, lon_deg, lat_deg, area, volume, SO2_density
 
+def read_DYNAMICO_native(dynamico_nc, area_nc, hybrids, day=7, hour=22):
+    def extend(x):
+        return x.reshape(x.shape + (1,)) 
+    hour = day*24+hour # assume output every hour, first day is day 0  
+    ds_area = nc.Dataset(area_nc)
+
+    lon  = ds_area['lon_i'][:]
+    lat  = ds_area['lat_i'][:]
+    area = ds_area['Ai'][0,:]
+    ncell = len(area)
+
+    ds   = nc.Dataset(dynamico_nc)
+    ps   = ds['ps'] [hour, :]
+    pmid = ds['P']  [hour, :, :]  # pressure at layers
+    phi  = ds['PHI'][hour, 1:, :] # geopotential'
+    q    = ds['q']  [hour, 0, :, :]
+    
+    print("Reading DYNAMICO native output with %d cells"%ncell)
+    print(area.shape)
+
+    lon, lat, area, ps, pmid, phi, q = map(extend, (lon, lat, area, ps, pmid, phi, q))
+
+    radius, volume, SO2_density = read_DYNAMICO_common(hybrids, hour, area, ps, pmid, phi, q)
+    return radius, lon, lat, area, volume, SO2_density
+
 def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
 
     hour = day*24+hour # assume output every hour
@@ -246,22 +278,24 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
     plt.close()
         
 def main():
-#    dynamico_nc      = '/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc'
-#    dynamico_hybrids = '/data/PLT8/pub/smailler/PUY60-dynamico/apbp.nc'
     dynamico_nc      = '/data/PLT8/pub/smailler/multiscale/dynamico.nc'
     dynamico_hybrids = '/data/PLT8/pub/smailler/multiscale/apbp.nc'
+    native_nc        = '/data/WORKSPACE/dubos/Transport/NetCDF/dynamico_native.nc'
+    native_area      = '/data/WORKSPACE/dubos/Transport/NetCDF/Ai.nc'
     wrf_nc           = '/home/smailler/PUY60/geog_USGS_PUY60.nc'
     chimere_nc       = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.201106%s00_24_PUY60-tier2.nc'
 
-    def read_DYN(day,hour):
-        return read_DYNAMICO(dynamico_nc, dynamico_hybrids, day, hour)
-    def read_CHI(day,hour):
-        return read_CHIMERE(wrf_nc, chimere_nc, day, hour)
-
     if args.amount_emitted:
+        def read_DYN(day,hour):
+            return read_DYNAMICO(dynamico_nc, dynamico_hybrids, day, hour)
+        def read_NAT(day,hour):
+            return read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids, day, hour)
+        def read_CHI(day,hour):
+            return read_CHIMERE(wrf_nc, chimere_nc, day, hour)
         plot_amount_emitted(read_DYN, read_CHI)
     else:
-        radius, _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
+#        radius, _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
+        radius, _, _, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids)
         check_area(area_DYNAMICO, volume_DYNAMICO, radius, 'DYNAMICO')
         _, _, area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc)    
         check_area(area_CHIMERE, volume_CHIMERE, radius, 'CHIMERE')