Initial (broken) version of half_plume.py

metrics/half_plume.py

+import numpy as np
+import netCDF4 as nc
+import matplotlib.pyplot as plt
+
+# flatten the conc arrays TAKE ONLY FINAL SIMULATION TIME, ONLY ONE SPECIES ON DYNAMICO SIDE
+
+def compute_pressure(ap, bp, p0, ps):
+    lev, lon, lat = len(ap), ps.shape[1], ps.shape[0]
+    p = np.zeros((lev, lat, lon)) # lat changes with row, lon with column
+    for i in range (0, lev): # 0 is ground level and not corresponding to a value of ilev?
+        p[i, :, :] = p0*ap[i] + ps*bp[i] # should return a pressure matrix(lat, lon) at each level
+    return p
+
+def check_compute_pressure(ap, bp, p0, ps, p):
+    pp = compute_pressure(ap, bp, p0, ps)
+    dp = np.absolute(pp-p)#.abs()
+    print("Checking formula for hybrid pressure coordinate : max(p)=%f, max(pp)=%f, max(dp)%f" % (p.max(), pp.max(), dp.max()) )
+    return dp.max()>1e-4  # true if there is a discrepancy
+
+def max_volume_plot(variable_chim, variable_dyn, day_list):
+    GTRC1_array = []
+    for day in day_list:
+        chimere_nc = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.201106%s00_24_PUY60-tier2.nc' %day
+        dsf = nc.Dataset(chimere_nc)    
+        GTRC1_array.append(np.amax(dsf[variable_chim][:,:,:])) # each time is in a separate file
+
+    ds = nc.Dataset('/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc')
+    q_array = ds[variable_dyn][:,0,:,:,:]
+    t_for_plot = []
+    q_for_plot = []
+    for t in range(0, q_array.shape[0], 48):
+        t_for_plot.append(t/48.)
+        q_for_plot.append(np.amax(ds[variable_dyn][t,0,:,:,:]))
+
+    plt.plot(t_for_plot, q_for_plot, '-r+')
+    plt.plot(day_list, GTRC1_array, '-b+')
+    plt.ylabel('Concentration')
+    plt.xlabel('t')
+    plt.legend()#([q_for_plot, GTRC1_array],['q', 'gtrc'])
+    plt.savefig('metric_figures/max_concentrations_versus_time.png')
+
+
+def read_DYNAMICO(dynamico_nc, dynamico_hybrids, hour=427):
+    g, p0 = 9.81, 1e5 # gravity, reference pressure for hybrid coordinate
+
+    print('Reading DYNAMICO data ...')
+
+    ds = nc.Dataset(dynamico_nc)
+    print(ds)
+    q = ds['q'][hour, 0, :, :, :].flatten()
+    ps = ds['ps'][hour, :, :]
+    phi = ds['PHI'][hour, :, :, :] # geopotential
+    print('  geopotential last hour = ', phi[:,12,27])
+
+    ds1 = nc.Dataset(dynamico_hybrids)
+    ap = ds1['hyai'][:].flatten() # 1D array, with ilev
+    bp = ds1['hybi'][:].flatten() # same
+
+    print('  ap shape, bp shape, ps shape = ',ap.shape, bp.shape, ps.shape)
+
+    lev = len(ap) - 1
+    lat = ds['ps'][hour, :, :].shape[0]
+    lon = ds['ps'][hour, :, :].shape[1]
+
+    p = np.zeros((lev, lat, lon)) # lat changes with row, lon with column
+    print('  p shape, q shape = ', p.shape, ds['q'][hour, 0, :, :, :].shape)
+
+    pp = ds['P'][hour,:,12,27]
+    print("  Pressure read from file :", pp)
+
+    if check_compute_pressure(ds1['hyam'][:], ds1['hybm'][:], p0, ps, ds['P'][hour,:,:,:] ):
+        print("  Discrepancy in computed pressure. Stop.")
+        return True
+
+    p = compute_pressure(ap,bp,p0, ps)
+#    for i in range (1, lev): # 0 is ground level and not corresponding to a value of ilev?
+#        p[i, :, :] = 10**5*ap[i] + ps*bp[i] # should return a pressure matrix(lat, lon) at each level
+
+    mass = np.zeros((lev, lat, lon))  # array of air mass per cell
+    volume = np.zeros((lev, lat, lon))  # array of volume per cell
+    for i in range (1,lev-2):
+        mass[i, :, :] = (p[i, :, :] - p[i+1, :, :]) / g # remove g factor?
+        volume[i, :, :] = (phi[i+1, :, :] - phi[i, :, :]) / g
+
+    volume = volume.flatten()
+    print('volume array length = ', len(volume), volume.size)
+    print('volume array = ', volume)
+
+    air_density = mass.flatten() / volume.flatten() # can i use chimere value?
+    SO2_density = q * air_density
+    return volume, SO2_density
+    
+
+def read_CHIMERE(chimere_nc, hour=24):
+    #	double airm(Time, bottom_top, south_north, west_east) ;
+    #		airm:units = "molecule/cm**3" ;
+    #		airm:long_name = "Air density" ;
+    #	float GTRC1(Time, bottom_top, south_north, west_east) ;
+    #		GTRC1:units = "ppb vol" ;
+    #		GTRC1:long_name = "GTRC1 Concentration" ;
+    # FIXME
+    # * we assume linear variation of pressure with model level
+    # * we neglect horizontal variations of cell area
+    # => all cells have the same air mass
+    # * we work with molecule number rather than mass (constant factor = molar mass)
+    print('Reading CHIMERE data ...')
+
+    ds2 = nc.Dataset(chimere_nc)
+    GTRC1 = ds2['GTRC1'][hour,:,:,:].flatten() # SO2 molecular mixing ratio (ppb)
+    air_density = ds2['airm'][hour,:,:,:].flatten()   # air molecule number per unit volume
+    SO2_density = GTRC1 * air_density                # SO2 molecule number per unit volume
+
+# *(64.066)*(1.67377*10**(-27)) molar mass x mass of 1 AMU (constant needed?) x number of air molecules per cell.
+    cell_mass = 1. # see FIXME above
+    cell_volume = cell_mass / air_density
+
+    return cell_volume, SO2_density
+
+def volume_half_mass(cell_volume, SO2_density):
+    # sort SO2 density per cell from small to large
+    ind = np.argsort(SO2_density) # indices sorting SO2 density
+    sorted_mass = SO2_density[ind]
+    sorted_volume = cell_volume[ind]
+    cum_mass   = np.cumsum(sorted_mass)
+    half_mass = 0.5*cum_mass[-1]
+
+    j = bisection_method(half_mass, cum_mass)
+
+    volume = sorted_volume[j:-1].sum()
+
+    print('half-mass volume = ', volume)
+
+
+def bisection_method(half_mass, cum_mass):
+
+    print('Dichotomy for half-mass index ... ')
+
+    N = len(cum_mass)
+    a = 0   # lower_bound
+    b = N-1 # upper_bound
+    m = int((a+b)/2)
+
+    if (cum_mass[m] <= half_mass) and (half_mass <= cum_mass[m+1]):
+        return m
+
+    else:
+        while (cum_mass[m] > half_mass) or (half_mass > cum_mass[m+1]):
+            print('a, b, m = ', a, b, m)
+            print('m, m+1, half_mass = ', cum_mass[m], cum_mass[m+1], half_mass)
+            if (cum_mass[m]>half_mass):
+                b = m # new upper bound
+            else:
+                a = m # new lower bound
+            m = int((a+b)/2)
+
+        print('a, b, m = ', a, b, m)
+        print('m, m+1, half_mass = ', cum_mass[m], cum_mass[m+1], half_mass)
+        return m    
+
+
+# --------------------------------------------------------------
+
+def main():
+    max_volume_plot('GTRC1', 'q', ['04', '05', '06', '07', '08', '09', '10', '11'])
+
+    dynamico_nc = '/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc'
+    dynamico_hybrids = '/data/PLT8/pub/smailler/PUY60-dynamico/apbp.nc'
+    data2 = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
+    volume_half_mass_DYNAMICO= volume_half_mass(*data2)
+
+    chimere_nc = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.2011061100_24_PUY60-tier2.nc'
+    data = read_CHIMERE(chimere_nc)
+    volume_half_mass_CHIMERE = volume_half_mass(*data)
+
+
+    #plots(metrics_CHIMERE, metrics_DYNAMICO)
+
+main()
+