Compute correct volumes and masses

Bug is there in computation of DYNAMICO area.

metrics/half_plume.py

@@ -13,11 +13,31 @@ def compute_pressure(ap, bp, p0, ps):
         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):
+def check_compute_pressure(ap, bp, p0, ps, p, rel_tol=1e-6):
     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
+    return dp.max()>p0*rel_tol  # true if there is a discrepancy
+
+def compute_volume(cell_area, hlay):
+    cell_area = 1 # FIXME
+    lev = hlay.shape[0]
+    volume = np.zeros(hlay.shape)
+    for l in range(lev):
+        if l==0:
+            height = hlay[0,:,:]
+        else:
+            height = hlay[l,:,:]-hlay[l-1,:,:]
+        vol =  cell_area * height
+        volume[l,:,:] = vol
+    return volume
+
+def compute_mass(p, g):
+    ilev, nlat, nlon = p.shape
+    mass = np.zeros((ilev-1, nlat, nlon))  # array of air mass per cell
+    for i in range(ilev-1):
+        mass[i, :, :] = (p[i, :, :] - p[i+1, :, :]) / g
+    return mass
 
 def max_volume_plot(variable_chim, variable_dyn, day_list):
     GTRC1_array = []
@@ -44,57 +64,52 @@ def max_volume_plot(variable_chim, variable_dyn, day_list):
 
 def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
     hour = 2*( (day-1)*24+hour ) # assume output every 30 min  
-    g, p0 = 9.81, 1e5 # gravity, reference pressure for hybrid coordinate
+    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)
-    print(ds)
-    q = ds['q'][hour, 0, :, :, :].flatten()
-    ps = ds['ps'][hour, :, :]
-    phi = ds['PHI'][hour, :, :, :] # geopotential
+    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)
+
+    lon, lat = np.meshgrid(lon, lat)
+    print('Min(cos(lat) = ', np.cos(lat).min() )
+    area = (radius**2)*np.cos(lat)*(2*np.pi/nlon)*(np.pi/nlat)
+
+    volume = compute_volume(area/g, phi)
+#    print('volume array length = ', len(volume), volume.size)
+#    print('volume array = ', volume)
+
+
     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
+    ap, bp = ds1['hyai'][:], ds1['hybi'][:]    
 
     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)
+    p = np.zeros((lev, nlat, nlon)) # 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)
+    p = compute_pressure(ap, bp, p0, ps)
+    air_mass = compute_mass(p, g)
+    SO2_density = q * air_mass / volume 
 
-    air_density = mass.flatten() / volume.flatten() # can i use chimere value?
-    SO2_density = q * air_density
-    return volume, SO2_density
-    
+    return radius, area, volume.flatten(), SO2_density.flatten()    
 
-def read_CHIMERE(chimere_nc, day=8, hour=24):
+def read_CHIMERE(dx_dy, wrf_nc, chimere_nc, day=8, hour=24):
     #	double airm(Time, bottom_top, south_north, west_east) ;
     #		airm:units = "molecule/cm**3" ;
     #		airm:long_name = "Air density" ;
@@ -106,33 +121,36 @@ def read_CHIMERE(chimere_nc, day=8, hour=24):
     # * 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)
-    filename = chimere_nc%CHIMERE_days[day-1]
-    print('Reading CHIMERE data from %s ...'%filename)
-    ds2 = nc.Dataset(filename)
-    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
+    print('Reading WRF data from %s ...'%wrf_nc)
+    wrf = nc.Dataset(wrf_nc)
+    cell_area = dx_dy * wrf['MAPFAC_MX'][0,:,:] * wrf['MAPFAC_MY'][0,:,:]
 
-# *(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)
+    def getvar(var):
+        return ds[var][hour,:,:,:]
 
+    filename = chimere_nc%CHIMERE_days[day-1]
+    print('Reading CHIMERE data from %s ...'%filename)
+    ds = nc.Dataset(filename)
+    hlay        = getvar('hlay')    # altitude above ground of upper interface
+    GTRC1       = getvar('GTRC1')   # SO2 molecular mixing ratio (ppb)
+    air_density = getvar('airm')    # air molecule number per unit volume
+    cell_volume = compute_volume(cell_area, hlay)
+    SO2_density = GTRC1 * air_density              # SO2 molecule number per unit volume
+
+    return cell_area, cell_volume.flatten(), SO2_density.flatten()
+
+def volume_half_mass(volume, density, label):
+    amount = volume*density
+    # sort density per cell from small to large
+    ind = np.argsort(density) # indices sorting density
+    sorted_amount = amount[ind]
+    sorted_volume = volume[ind]
+    cum_amount    = np.cumsum(amount)
+    half_amount   = 0.5*cum_amount[-1]
+
+    j = bisection_method(half_amount, cum_amount)
     volume = sorted_volume[j:-1].sum()
-
-    print('half-mass volume = ', volume)
-
+    print('(%s) half-mass volume = %d '%(label, volume))
 
 def bisection_method(half_mass, cum_mass):
 
@@ -166,18 +184,22 @@ def bisection_method(half_mass, cum_mass):
 import argparse
 
 def main():
-    dynamico_nc = '/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc'
+    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)
+    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'
+    dx, dy = 6e4, 6e4 # nominal CHIMERE cell sizes in m
 
-    chimere_nc = '/data/PLT8/pub/smailler/CHIMOUT-PUY60-tier2/out.201106%s00_24_PUY60-tier2.nc'
-    data = read_CHIMERE(chimere_nc)
-    volume_half_mass_CHIMERE = volume_half_mass(*data)
+    radius, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
+    print('Total DYNAMICO domain area/theory =', area_DYNAMICO.sum() / (4*np.pi*radius**2) )
+    volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')
 
-    max_volume_plot('GTRC1', 'q', CHIMERE_days)
+    area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(dx*dy, wrf_nc, chimere_nc)    
+    print('Total CHIMERE domain area =', area_CHIMERE.sum() )
+    volume_half_mass_CHIMERE = volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE')
+
+    #max_volume_plot('GTRC1', 'q', CHIMERE_days)
 
     #plots(metrics_CHIMERE, metrics_DYNAMICO)
 
 main()
-