Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 6f8d7126 authored by Thomas Dubos's avatar Thomas Dubos
Browse files

Plot total amount of SO2 vs time - units not right yet

parent 3e1c7438
......@@ -4,8 +4,6 @@ import matplotlib.pyplot as plt
CHIMERE_days = ['04', '05', '06', '07', '08', '09', '10', '11']
# 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
......@@ -63,8 +61,8 @@ 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=8, hour=22):
hour = 2*( (day-1)*24+hour ) # assume output every 30 min
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)
......@@ -75,23 +73,24 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
q = ds['q'][hour, 0, :, :, :]
ps = ds['ps'][hour, :, :]
phi = ds['PHI'][hour, 1:, :, :] # geopotential'
print('shape(phi) = ', phi.shape)
# print('shape(phi) = ', phi.shape)
# print("Max of DYNAMICO a at hour=%d : %f"%(hour, q.max()))
lon, lat = np.meshgrid(lon, lat)
print('Min(cos(lat) = ', np.cos(lat).min() )
# 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)
volume = compute_volume(area/g, phi) # m^3
# print('volume array length = ', len(volume), volume.size)
# print('volume array = ', volume)
print(' geopotential last hour = ', phi[:,12,27])
# print(' geopotential last hour = ', phi[:,12,27])
ds1 = nc.Dataset(dynamico_hybrids)
ap, bp = ds1['hyai'][:], ds1['hybi'][:]
print(' ap shape, bp shape, ps shape = ',ap.shape, bp.shape, ps.shape)
# print(' ap shape, bp shape, ps shape = ',ap.shape, bp.shape, ps.shape)
lev = len(ap) - 1
......@@ -103,29 +102,39 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
p = compute_pressure(ap, bp, p0, ps)
air_mass = compute_mass(p, g)
SO2_density = q * air_mass / volume
SO2_density = q * air_mass / volume # (kg/kg) * (kg) * (m/3) = kg / m^3
return radius, area, volume.flatten(), SO2_density.flatten()
def read_CHIMERE(wrf_nc, chimere_nc, day=8, hour=24):
def read_CHIMERE(wrf_nc, chimere_nc, day=7, hour=24):
hour = day*24+hour # assume output every hour
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)
wrf = nc.Dataset(wrf_nc)
cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:]
cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:] # m^2
def getvar(var):
return ds[var][hour,:,:,:]
filename = chimere_nc%CHIMERE_days[day-1]
filename = chimere_nc%CHIMERE_days[day]
print('Reading CHIMERE data from %s ...'%filename)
ds = nc.Dataset(filename)
hlay = getvar('hlay') # altitude above ground of upper interface
hlay = getvar('hlay') # altitude above ground of upper interface (m)
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
air_density = getvar('airm') # air molecule number per unit volume (cm^3)
cell_volume = compute_volume(cell_area, hlay) # m^3
return cell_area, cell_volume.flatten(), SO2_density.flatten()
SO2_molar_mass = 64.066/1e3 # kg/mol
avogadro = 6.02214076e23 # Avogadro number
air_density = air_density*1e6 # air molecule number per unit volume (m^3)
SO2_density = 1e-9*GTRC1*air_density # SO2 molecule number per unit volume (m^3)
SO2_density = SO2_density*SO2_molar_mass/avogadro # SO2 mass per unit volume (m^3)
return cell_area, cell_volume.flatten(), SO2_density.flatten()
def volume_half_mass(volume, density, label):
amount = volume*density
......@@ -166,29 +175,70 @@ def bisection_method(half_mass, cum_mass):
print('m, m+1, half_mass = ', cum_mass[m], cum_mass[m+1], half_mass)
return m
def check_area(area, radius, label):
tot = area.sum()
print('Total %s domain area= %f = %d percent of Earth surface'%(label, tot, 100*tot/(4*np.pi*radius**2) ))
# --------------------------------------------------------------
import argparse
import matplotlib.pyplot as plt
def check_area(area, radius, label):
tot = area.sum()
print('Total %s domain area= %f = %d percent of Earth surface'%(label, tot, 100*tot/(4*np.pi*radius**2) ))
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')
args = getargs.parse_args()
def plot_amount_emitted(read_DYN, read_CHI, ndays=4, hour_step=3):
n=24*ndays//hour_step
hours = [ i*hour_step for i in range(n)]
amount_DYN=np.zeros(n)
amount_CHI=np.zeros(n)
def amount(vol, dens): return (vol*dens).sum()
for i in range(n):
radius, area_DYN, volume_DYN, density_DYN = read_DYN(3, hours[i]) # we shift the DYNAMICO data by 3 days
area_CHI, volume_CHI, density_CHI = read_CHI(0, hours[i])
amount_DYN[i] = amount(volume_DYN, density_DYN)
amount_CHI[i] = amount(volume_CHI, density_CHI)
# print(amount_DYN)
# print(amount_CHI)
plt.figure();
plt.plot(hours, amount_DYN, label='DYNAMICO')
plt.plot(hours, amount_CHI, label='CHIMERE')
plt.yscale('log')
plt.xlabel('hour')
# plt.legend()
plt.title('Total amount of SO2 (should be kg)')
plt.savefig('amount_emitted.png')
def main():
dynamico_nc = '/data/PLT8/pub/smailler/PUY60-dynamico/dynamico.nc'
dynamico_hybrids = '/data/PLT8/pub/smailler/PUY60-dynamico/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'
radius, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
check_area(area_DYNAMICO, radius, 'DYNAMICO')
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:
plot_amount_emitted(read_DYN, read_CHI)
else:
radius, area_DYNAMICO, volume_DYNAMICO, density_DYNAMICO = read_DYNAMICO(dynamico_nc, dynamico_hybrids)
check_area(area_DYNAMICO, radius, 'DYNAMICO')
volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')
volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')
area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc)
check_area(area_CHIMERE, radius, 'CHIMERE')
volume_half_mass_CHIMERE = volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE')
area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc)
check_area(area_CHIMERE, radius, 'CHIMERE')
volume_half_mass_CHIMERE = volume_half_mass(volume_CHIMERE, density_CHIMERE, 'CHIMERE')
#max_volume_plot('GTRC1', 'q', CHIMERE_days)
......
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