Commit 33db3c32 authored by Thomas Dubos's avatar Thomas Dubos
Browse files

Fix missing area factor in computation of DYNAMICO mass

parent 24c61f76
......@@ -14,7 +14,7 @@ def compute_pressure(ap, bp, p0, ps):
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()) )
# 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()>p0*rel_tol # true if there is a discrepancy
def compute_volume(cell_area, hlay):
......@@ -29,11 +29,11 @@ def compute_volume(cell_area, hlay):
volume[l,:,:] = vol
return volume
def compute_mass(p, g):
def compute_mass(cell_area, 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
mass[i, :, :] = cell_area * (p[i, :, :] - p[i+1, :, :]) / g
return mass
def max_volume_plot(variable_chim, variable_dyn, day_list):
......@@ -65,7 +65,7 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
hour = day*24+hour # assume output every hour, first day is day 0
radius, g, p0 = 6.371e6, 9.81, 1e5 # Earth radius, gravity, reference pressure for hybrid coordinate
print('Reading DYNAMICO data from %s ...'%dynamico_nc)
# print('Reading DYNAMICO data from %s ...'%dynamico_nc)
ds = nc.Dataset(dynamico_nc)
lon, lat = to_radian(ds['lon'][:]), to_radian(ds['lat'][:])
......@@ -101,9 +101,11 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
return True
p = compute_pressure(ap, bp, p0, ps)
air_mass = compute_mass(p, g)
SO2_density = q * air_mass / volume # (kg/kg) * (kg) * (m/3) = kg / m^3
air_mass = compute_mass(area, p, g)
SO2_density = q * air_mass / volume # (kg/kg) * (kg) / (m^3) = kg / m^3
# print("Total SO2 mass in DYNAMICO = ", (q*air_mass).sum())
# print("Max SO2 mixing ratio (kg/kg) in DYNAMICO = ", q.max())
return radius, area, volume.flatten(), SO2_density.flatten()
def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
......@@ -112,7 +114,7 @@ def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
day, hour = hour//24, hour%24
# NOTE :we work with molecule number rather than mass (constant factor = molar mass)
print('Reading WRF data from %s ...'%wrf_nc)
# print('Reading WRF data from %s ...'%wrf_nc)
wrf = nc.Dataset(wrf_nc)
cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:] # m^2
......@@ -120,7 +122,7 @@ def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
return ds[var][hour,:,:,:]
filename = chimere_nc%CHIMERE_days[day]
print('Reading CHIMERE data from %s ...'%filename)
# print('Reading CHIMERE data from %s ...'%filename)
ds = nc.Dataset(filename)
hlay = getvar('hlay') # altitude above ground of upper interface (m)
GTRC1 = getvar('GTRC1') # SO2 molecular mixing ratio (ppb)
......@@ -148,6 +150,7 @@ def volume_half_mass(volume, density, label):
j = bisection_method(half_amount, cum_amount)
volume = sorted_volume[j:-1].sum()
print('(%s) half-mass = %f '%(label, half_amount))
print('(%s) half-mass volume = %f '%(label, volume/total_volume))
def bisection_method(half_mass, cum_mass):
......@@ -188,8 +191,8 @@ import argparse
import matplotlib.pyplot as plt
getargs = argparse.ArgumentParser(description='Analyze 3D transport experiments.')
getargs.add_argument("--amount_emitted", action='store_true', help='Total emitted SO2 vs time')
getargs.add_argument("--max_volume", action='store_true', help='Total emitted SO2 vs time')
getargs.add_argument("--amount-emitted", action='store_true', help='Total emitted SO2 vs time')
getargs.add_argument("--half-mass-volume", action='store_true', help='Minimal volume containing half of emitted SO2 mass')
args = getargs.parse_args()
def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
......@@ -221,8 +224,10 @@ def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
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/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'
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'
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment