[3298] | 1 | #################################################################################### |
---|
| 2 | ### Python script to output the stratification data from the "startpem.nc" files ### |
---|
| 3 | #################################################################################### |
---|
| 4 | |
---|
| 5 | import netCDF4 as nc |
---|
| 6 | import numpy as np |
---|
| 7 | import matplotlib.pyplot as plt |
---|
| 8 | import sys |
---|
[3300] | 9 | import os.path |
---|
[3298] | 10 | |
---|
| 11 | ############################## |
---|
| 12 | ### Parameters to fill in |
---|
[3458] | 13 | filename = 'startpem.nc' #'starts/restartpem226.nc' # File name |
---|
[3298] | 14 | igrid = 1 # Grid point |
---|
| 15 | islope = 1 # Slope number |
---|
| 16 | istr = 4 # Stratum number |
---|
| 17 | ############################## |
---|
| 18 | |
---|
[3458] | 19 | |
---|
| 20 | #################################################################################### |
---|
[3298] | 21 | ### Open the NetCDF file |
---|
[3300] | 22 | if os.path.isfile(filename): |
---|
| 23 | nc_file = nc.Dataset(filename,'r') |
---|
| 24 | else: |
---|
| 25 | sys.exit('The file \"' + filename + '\" does not exist!') |
---|
[3298] | 26 | |
---|
| 27 | ### Get the dimensions |
---|
| 28 | Time = len(nc_file.dimensions['Time']) |
---|
| 29 | ngrid = len(nc_file.dimensions['physical_points']) |
---|
| 30 | nslope = len(nc_file.dimensions['nslope']) |
---|
| 31 | nb_str_max = len(nc_file.dimensions['nb_str_max']) |
---|
| 32 | if igrid > ngrid or igrid < 1: |
---|
[3300] | 33 | sys.exit('Asked grid point is not possible!') |
---|
[3298] | 34 | if islope > nslope or islope < 1: |
---|
[3300] | 35 | sys.exit('Asked slope number is not possible!') |
---|
[3298] | 36 | if istr > nb_str_max or istr < 1: |
---|
[3300] | 37 | sys.exit('Asked stratum number is not possible!') |
---|
[3298] | 38 | |
---|
| 39 | ### Get the stratification properties |
---|
| 40 | stratif_thickness = [] |
---|
| 41 | stratif_top_elevation = [] |
---|
| 42 | stratif_co2ice_volfrac = [] |
---|
| 43 | stratif_h2oice_volfrac = [] |
---|
| 44 | stratif_dust_volfrac = [] |
---|
| 45 | stratif_air_volfrac = [] |
---|
| 46 | for i in range(1,nslope + 1): |
---|
| 47 | stratif_thickness.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_thickness'][:]) |
---|
| 48 | stratif_top_elevation.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_top_elevation'][:]) |
---|
| 49 | stratif_co2ice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_co2ice_volfrac'][:]) |
---|
| 50 | stratif_h2oice_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_h2oice_volfrac'][:]) |
---|
| 51 | stratif_dust_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_dust_volfrac'][:]) |
---|
| 52 | stratif_air_volfrac.append(nc_file.variables['stratif_slope' + str(i).zfill(2) + '_air_volfrac'][:]) |
---|
| 53 | |
---|
| 54 | ### Display the data |
---|
| 55 | igrid = igrid - 1 |
---|
| 56 | islope = islope - 1 |
---|
| 57 | labels = 'CO2 ice', 'H2O ice', 'Dust', 'Air' |
---|
[3300] | 58 | contents = np.zeros([4,len(stratif_top_elevation[islope][0,:,igrid]) + 1]) |
---|
| 59 | height = np.zeros(len(stratif_top_elevation[islope][0,:,igrid]) + 1) |
---|
| 60 | height[0] = stratif_top_elevation[islope][0,0,:] - stratif_thickness[islope][0,0,:] |
---|
| 61 | height[1:] = stratif_top_elevation[islope][0,:,igrid] |
---|
[3298] | 62 | for i in range(len(stratif_top_elevation[islope][0,:,igrid])): |
---|
[3300] | 63 | contents[0,1 + i] = stratif_co2ice_volfrac[islope][0,i,igrid] |
---|
| 64 | contents[1,1 + i] = stratif_h2oice_volfrac[islope][0,i,igrid] |
---|
| 65 | contents[2,1 + i] = stratif_dust_volfrac[islope][0,i,igrid] |
---|
| 66 | contents[3,1 + i] = stratif_air_volfrac[islope][0,i,igrid] |
---|
| 67 | contents[:,0] = contents[:,1] |
---|
[3298] | 68 | |
---|
[3300] | 69 | # Simple subplots for a layering |
---|
| 70 | fig, (ax1, ax2, ax3, ax4) = plt.subplots(1,4,layout = 'constrained',sharey = True) |
---|
| 71 | fig.suptitle('Simple content profiles for the layering') |
---|
| 72 | ax1.step(contents[0,:],height,where = 'post') |
---|
| 73 | ax2.step(contents[1,:],height,where = 'post') |
---|
| 74 | ax3.step(contents[2,:],height,where = 'post') |
---|
| 75 | ax4.step(contents[3,:],height,where = 'post') |
---|
| 76 | ax1.set_ylabel('Elevation [m]') |
---|
| 77 | ax1.set_xlabel('Volume fraction [m3/m3]') |
---|
| 78 | ax2.set_xlabel('Volume fraction [m3/m3]') |
---|
| 79 | ax3.set_xlabel('Volume fraction [m3/m3]') |
---|
| 80 | ax4.set_xlabel('Volume fraction [m3/m3]') |
---|
[3298] | 81 | ax1.set_title(labels[0]) |
---|
| 82 | ax2.set_title(labels[1]) |
---|
| 83 | ax3.set_title(labels[2]) |
---|
| 84 | ax4.set_title(labels[3]) |
---|
[3300] | 85 | plt.savefig('layering_simpleprofiles.png') |
---|
[3298] | 86 | |
---|
[3458] | 87 | # Content profiles for a layering |
---|
| 88 | plt.figure() |
---|
| 89 | plt.step(contents[0,:],height,where = 'post',color = 'r',label = labels[0]) |
---|
| 90 | #plt.plot(contents[0,:],height,'o--',color = 'r',alpha = 0.3) |
---|
| 91 | plt.step(contents[1,:],height,where = 'post',color = 'b',label = labels[1]) |
---|
| 92 | #plt.plot(contents[1,:],height,'o--',color = 'b',alpha = 0.3) |
---|
| 93 | plt.step(contents[2,:],height,where = 'post',color = 'y',label = labels[2]) |
---|
| 94 | #plt.plot(contents[2,:],height,'o--',color = 'y',alpha = 0.3) |
---|
| 95 | plt.step(contents[3,:],height,where = 'post',color = 'g',label = labels[3]) |
---|
| 96 | #plt.plot(contents[3,:],height,'o--',color = 'g',alpha = 0.3) |
---|
| 97 | plt.grid(axis='x', color='0.95') |
---|
| 98 | plt.grid(axis='y', color='0.95') |
---|
| 99 | plt.xlim(0,1) |
---|
| 100 | plt.xlabel('Volume fraction [m3/m3]') |
---|
| 101 | plt.ylabel('Elevation [m]') |
---|
| 102 | plt.title('Content profiles for the layering') |
---|
| 103 | plt.savefig('layering_simpleprofiles.png') |
---|
| 104 | |
---|
| 105 | # Stack content profiles for a layering |
---|
[3298] | 106 | fig, ax = plt.subplots() |
---|
[3300] | 107 | ax.fill_betweenx(height,0,contents[0,:],label = labels[0],color = 'r',step = 'pre') |
---|
| 108 | ax.fill_betweenx(height,contents[0,:],sum(contents[0:2,:]),label = labels[1],color = 'b',step = 'pre') |
---|
| 109 | ax.fill_betweenx(height,sum(contents[0:2,:]),sum(contents[0:3,:]),label = labels[2],color = 'y',step = 'pre') |
---|
| 110 | ax.fill_betweenx(height,sum(contents[0:3,:]),sum(contents),label = labels[3],color = 'g',step = 'pre') |
---|
| 111 | plt.vlines(x = 0.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-') |
---|
| 112 | plt.vlines(x = 1.,ymin = height[0],ymax = height[len(height) - 1],color = 'k',linestyle = '-') |
---|
| 113 | for i in range(len(height)): |
---|
| 114 | plt.hlines(y = height[i],xmin = 0.0,xmax = 1.0,color = 'k',linestyle = '--') |
---|
| 115 | ax.set_title('Stack content profiles for the layering') |
---|
| 116 | plt.xlabel('Volume fraction [m3/m3]') |
---|
| 117 | plt.ylabel('Elevation [m]') |
---|
[3298] | 118 | ax.legend(loc = 'center left',bbox_to_anchor = (1,0.5)) |
---|
[3300] | 119 | plt.savefig('layering_stackprofiles.png') |
---|
[3298] | 120 | |
---|
| 121 | plt.show() |
---|
| 122 | |
---|
| 123 | ### Close the NetCDF file |
---|
| 124 | nc_file.close() |
---|