Check mean model top

metrics/half_plume.py

@@ -106,7 +106,7 @@ def read_DYNAMICO(dynamico_nc, dynamico_hybrids, day=7, hour=22):
 
     return radius, area, volume.flatten(), SO2_density.flatten()    
 
-def read_CHIMERE(wrf_nc, chimere_nc, day=7, hour=24):
+def read_CHIMERE(wrf_nc, chimere_nc, day=6, hour=24):
 
     hour = day*24+hour # assume output every hour
     day, hour = hour//24, hour%24
@@ -137,6 +137,7 @@ def read_CHIMERE(wrf_nc, chimere_nc, day=7, hour=24):
     return cell_area, cell_volume.flatten(), SO2_density.flatten() 
 
 def volume_half_mass(volume, density, label):
+    total_volume=volume.sum()
     amount = volume*density
     # sort density per cell from small to large
     ind = np.argsort(density) # indices sorting density
@@ -147,11 +148,11 @@ def volume_half_mass(volume, density, label):
 
     j = bisection_method(half_amount, cum_amount)
     volume = sorted_volume[j:-1].sum()
-    print('(%s) half-mass volume = %d '%(label, volume))
+    print('(%s) half-mass volume = %f '%(label, volume/total_volume))
 
 def bisection_method(half_mass, cum_mass):
 
-    print('Dichotomy for half-mass index ... ')
+#    print('Dichotomy for half-mass index ... ')
 
     N = len(cum_mass)
     a = 0   # lower_bound
@@ -163,21 +164,23 @@ def bisection_method(half_mass, cum_mass):
 
     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)
+#            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)
+#        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 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 check_area(area, volume, radius, label):
+    area = area.sum()
+    print('Total %s domain area= %f = %d percent of Earth surface'%(label, area, 100*area/(4*np.pi*radius**2) )) 
+    volume = volume.sum()
+    print('Mean model top for %s is at %d meters'%(label, volume/area))
 
 # --------------------------------------------------------------
 
@@ -232,12 +235,11 @@ def main():
         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')
+        check_area(area_DYNAMICO, volume_DYNAMICO, radius, 'DYNAMICO')
+        area_CHIMERE, volume_CHIMERE, density_CHIMERE = read_CHIMERE(wrf_nc, chimere_nc)    
+        check_area(area_CHIMERE, volume_CHIMERE, radius, 'CHIMERE')
 
         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')
 
     #max_volume_plot('GTRC1', 'q', CHIMERE_days)