Fix map factor (CHIMERE) and degrees->radians (DYNAMICO) in area computation

metrics/half_plume.py

@@ -20,7 +20,6 @@ def check_compute_pressure(ap, bp, p0, ps, p, rel_tol=1e-6):
     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):
@@ -61,6 +60,8 @@ def max_volume_plot(variable_chim, variable_dyn, day_list):
     plt.legend()#([q_for_plot, GTRC1_array],['q', 'gtrc'])
     plt.savefig('metric_figures/max_concentrations_versus_time.png')
 
+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  
@@ -69,7 +70,7 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
     print('Reading DYNAMICO data from %s ...'%dynamico_nc)
 
     ds = nc.Dataset(dynamico_nc)
-    lon, lat = ds['lon'][:], ds['lat'][:]
+    lon, lat = to_radian(ds['lon'][:]), to_radian(ds['lat'][:])
     nlon, nlat = len(lon), len(lat)
     q   = ds['q'][hour, 0, :, :, :]
     ps  = ds['ps'][hour, :, :]
@@ -95,9 +96,6 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
     lev = len(ap) - 1
 
     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.")
@@ -109,21 +107,11 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=8, hour=22):
 
     return radius, area, volume.flatten(), SO2_density.flatten()    
 
-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" ;
-    #	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)
+def read_CHIMERE(wrf_nc, chimere_nc, day=8, 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 = dx_dy * wrf['MAPFAC_MX'][0,:,:] * wrf['MAPFAC_MY'][0,:,:]
+    cell_area = wrf.DX*wrf.DY / wrf['MAPFAC_MX'][0,:,:]/ wrf['MAPFAC_MY'][0,:,:]
 
     def getvar(var):
         return ds[var][hour,:,:,:]
@@ -183,19 +171,23 @@ def bisection_method(half_mass, cum_mass):
 
 import argparse
 
+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) ))
+ 
 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'
-    dx, dy = 6e4, 6e4 # nominal CHIMERE cell sizes in m
 
     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) )
+    check_area(area_DYNAMICO, radius, 'DYNAMICO')
+
     volume_half_mass_DYNAMICO = volume_half_mass(volume_DYNAMICO, density_DYNAMICO, 'DYNAMICO')
 
-    area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(dx*dy, wrf_nc, chimere_nc)    
-    print('Total CHIMERE domain area =', area_CHIMERE.sum() )
+    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)