Changeset 3300 for trunk/LMDZ.COMMON/libf
- Timestamp:
- Apr 15, 2024, 9:42:47 AM (8 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3298 r3300 274 274 == 10/04/2024 == JBC 275 275 Addition of a Python script "output_layering.py" in the deftank which allows to output layering data nicely from the "startpem.nc" files. 276 277 == 15/04/2024 == JBC 278 Update of "output_layering.py" in the deftank. -
trunk/LMDZ.COMMON/libf/evolution/deftank/output_layering.py
r3298 r3300 7 7 import matplotlib.pyplot as plt 8 8 import sys 9 import os.path 9 10 10 11 ############################## 11 12 ### Parameters to fill in 12 filename = 'startpem 3.nc' # File name13 filename = 'startpem9.nc' # File name 13 14 igrid = 1 # Grid point 14 15 islope = 1 # Slope number … … 17 18 18 19 ### Open the NetCDF file 19 nc_file = nc.Dataset(filename,'r') 20 if os.path.isfile(filename): 21 nc_file = nc.Dataset(filename,'r') 22 else: 23 sys.exit('The file \"' + filename + '\" does not exist!') 20 24 21 25 ### Get the dimensions … … 25 29 nb_str_max = len(nc_file.dimensions['nb_str_max']) 26 30 if igrid > ngrid or igrid < 1: 27 sys.exit( "Asked grid point is not possible!.")31 sys.exit('Asked grid point is not possible!') 28 32 if islope > nslope or islope < 1: 29 sys.exit( "Asked slope number is not possible!")33 sys.exit('Asked slope number is not possible!') 30 34 if istr > nb_str_max or istr < 1: 31 sys.exit("Asked stratum number is not possible!")35 sys.exit('Asked stratum number is not possible!') 32 36 33 37 ### Get the stratification properties … … 49 53 igrid = igrid - 1 50 54 islope = islope - 1 51 istr = istr - 152 55 labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air' 53 contents = stratif_co2ice_volfrac[islope][0,:,igrid], stratif_h2oice_volfrac[islope][0,:,igrid], stratif_dust_volfrac[islope][0,:,igrid], stratif_air_volfrac[islope][0,:,igrid] 54 x = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1]) 55 y = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1) 56 y[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:] 57 y[1:] = stratif_top_elevation[islope][0,:,igrid] 56 contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1]) 57 height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1) 58 height[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:] 59 height[1:] = stratif_top_elevation[islope][0,:,igrid] 58 60 for i in range(len(stratif_top_elevation[islope][0,:,igrid])): 59 x[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid]60 x[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid]61 x[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid]62 x[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid]63 x[:,0] = x[:,1]61 contents[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid] 62 contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid] 63 contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid] 64 contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid] 65 contents[:,0] = contents[:,1] 64 66 65 # Simple multiple subplots for a layering 66 fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4) 67 fig.suptitle('Volume fractions [m3/m3] in the layering') 68 ax1.plot(contents[0],stratif_top_elevation[islope][0,:,igrid]) 69 ax2.plot(contents[1],stratif_top_elevation[islope][0,:,igrid]) 70 ax3.plot(contents[2],stratif_top_elevation[islope][0,:,igrid]) 71 ax4.plot(contents[3],stratif_top_elevation[islope][0,:,igrid]) 67 # Simple subplots for a layering 68 fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,layout = 'constrained',sharey = True) 69 fig.suptitle('Simple content profiles for the layering') 70 ax1.step(contents[0,:],height,where = 'post') 71 ax2.step(contents[1,:],height,where = 'post') 72 ax3.step(contents[2,:],height,where = 'post') 73 ax4.step(contents[3,:],height,where = 'post') 74 ax1.set_ylabel('Elevation [m]') 75 ax1.set_xlabel('Volume fraction [m3/m3]') 76 ax2.set_xlabel('Volume fraction [m3/m3]') 77 ax3.set_xlabel('Volume fraction [m3/m3]') 78 ax4.set_xlabel('Volume fraction [m3/m3]') 72 79 ax1.set_title(labels[0]) 73 80 ax2.set_title(labels[1]) 74 81 ax3.set_title(labels[2]) 75 82 ax4.set_title(labels[3]) 76 77 # Pie chart for a stratum 78 fig, ax = plt.subplots(figsize = (6, 3),subplot_kw = dict(aspect = "equal")) 79 def func(pct,allvals): 80 absolute = int(np.round(pct/100.*np.sum(allvals))) 81 return f"{pct:.1f}%\n({absolute:d} m3/m3)" 82 wedges, texts, autotexts = ax.pie(x[:,istr + 1],autopct = lambda pct: func(pct,x[:,istr + 1]),textprops = dict(color = "w")) 83 ax.legend(wedges,labels,title = "Content",loc = "center left",bbox_to_anchor = (1, 0, 0.5, 1)) 84 plt.setp(autotexts,size = 8,weight = "bold") 85 ax.set_title("Content of the stratum " + str(istr + 1)) 83 plt.savefig('layering_simpleprofiles.png') 86 84 87 85 # Stackplot for a layering 88 86 fig, ax = plt.subplots() 89 ax.fill_betweenx( y,0,x[0,:],label = labels[0],step = 'pre')90 ax.fill_betweenx( y,x[0,:],sum(x[0:2,:]),label = labels[1],step = 'pre')91 ax.fill_betweenx( y,sum(x[0:2,:]),sum(x[0:3,:]),label = labels[2],step = 'pre')92 ax.fill_betweenx( y,sum(x[0:3,:]),sum(x),label = labels[3],step = 'pre')93 plt.vlines(x = 0.,ymin = y[0],ymax = y[len(y) - 1],color = 'k',linestyle = '-')94 plt.vlines(x = 1.,ymin = y[0],ymax = y[len(y) - 1],color = 'k',linestyle = '-')95 for i in range(len( y)):96 plt.hlines(y = y[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--')97 ax.set_title( "Layering")98 plt.xlabel( "Volume fraction [m3/m3]")99 plt.ylabel( "Elevation [m]")87 ax.fill_betweenx(height,0,contents[0,:],label = labels[0],color = 'r',step = 'pre') 88 ax.fill_betweenx(height,contents[0,:],sum(contents[0:2,:]),label = labels[1],color = 'b',step = 'pre') 89 ax.fill_betweenx(height,sum(contents[0:2,:]),sum(contents[0:3,:]),label = labels[2],color = 'y',step = 'pre') 90 ax.fill_betweenx(height,sum(contents[0:3,:]),sum(contents),label = labels[3],color = 'g',step = 'pre') 91 plt.vlines(x = 0.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-') 92 plt.vlines(x = 1.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-') 93 for i in range(len(height)): 94 plt.hlines(y = height[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--') 95 ax.set_title('Stack content profiles for the layering') 96 plt.xlabel('Volume fraction [m3/m3]') 97 plt.ylabel('Elevation [m]') 100 98 ax.legend(loc = 'center left',bbox_to_anchor = (1,0.5)) 99 plt.savefig('layering_stackprofiles.png') 101 100 102 101 plt.show()
Note: See TracChangeset
for help on using the changeset viewer.