- Timestamp:
- Oct 15, 2024, 1:41:48 PM (2 months ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3458 r3460 445 445 == 11/10/2024 == JBC 446 446 Addition of a python script to visualize the layerings over time and on a reference grid for elevation which is specified by the user. 447 448 == 15/10/2024 == JBC 449 Update and correction for the visualization of the layerings over time. -
trunk/LMDZ.COMMON/libf/evolution/deftank/README
r3403 r3460 49 49 Bash script file to clean the folder after a PEM simulation. 50 50 51 # output_layering.py: 51 # multiple_exec.sh: 52 Bash script to execute multiple scripts in subdirectories, useful to launch multiple simulations at once. 53 54 # visu_layering.py: 52 55 Python script file to output the stratification data from the "startpem.nc" files. 56 57 # visu_evol_layering.py: 58 Python script file to output the stratification data over time from the "startpem.nc" files. 53 59 54 60 Note: -
trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py
r3458 r3460 13 13 folder_path = "starts" 14 14 base_name = "restartpem" 15 dz = 1 # Discrization step of the reference grid for the elevation [m]15 dz = 0.1 # Discrization step of the reference grid for the elevation [m] 16 16 infofile = 'info_PEM.txt' 17 17 ######################### … … 39 39 ngrid = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['physical_points'].size # ngrid is the same for all files 40 40 nslope = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['nslope'].size # nslope is the same for all files 41 longitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['longitude'][:] 42 latitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['latitude'][:] 41 43 for i in range(1,nfile + 1): 42 44 file_path = os.path.join(folder_path,base_name + str(i) + ".nc") … … 53 55 max_top_elevation = max(max_top_elevation,np.max(ds.variables[slope_var_name][:])) 54 56 55 print(f" 56 print(f" 57 print(f" 58 print(f" 59 print(f" 57 print(f"> number of files = {nfile}") 58 print(f"> ngrid = {ngrid}") 59 print(f"> nslope = {nslope}") 60 print(f"> max(nb_str_max) = {max_nb_str}") 61 print(f"> max(top_elevation) = {max_top_elevation}") 60 62 61 63 # Concatenate stratif variables with dimension 'nb_str_max' along the "Time" dimension 62 64 stratif_data = [] 63 65 stratif_heights = np.zeros((ngrid,nfile,nslope,max_nb_str)) 64 stratif_co2ice = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))65 stratif_h2oice = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))66 stratif_dust = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))67 stratif_air = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))66 stratif_co2ice = np.zeros((ngrid,nfile,nslope,max_nb_str)) 67 stratif_h2oice = np.zeros((ngrid,nfile,nslope,max_nb_str)) 68 stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str)) 69 stratif_air = np.zeros((ngrid,nfile,nslope,max_nb_str)) 68 70 for var_name in datasets[0].variables: 69 71 if 'top_elevation' in var_name: … … 108 110 109 111 stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_air] 110 return stratif_data, max_top_elevation 112 return stratif_data, max_top_elevation, longitude, latitude 111 113 112 114 ### Function to interpolate the stratification data on a reference grid … … 114 116 # Define the reference ref_grid 115 117 ref_grid = np.arange(0,max_top_elevation,dz) 116 print(f" > Number of ref_grid points = {len(ref_grid)}") 117 118 # Extrapolate the top_elevation of strata 119 for i in range(np.shape(stratif_data)[1]): 120 for j in range(np.shape(stratif_data)[2]): 121 for k in range(np.shape(stratif_data)[3]): 122 if len(np.where(stratif_data[0][i,j,k,:] == 0)[0]) > 1: # If there are trailing zeros in top_elevation 123 i0 = np.where(stratif_data[0][i,j,k,:] == 0)[0][1] 124 stratif_data[0][i,j,k,i0:] = np.linspace(stratif_data[0][i,j,k,i0 - 1],max_top_elevation,np.shape(stratif_data)[4] - i0) 118 print(f"> Number of ref_grid points = {len(ref_grid)}") 125 119 126 120 # Interpolate the strata properties on the ref_grid 127 gridded_stratif_data = np.zeros((np.shape(stratif_data)[0] - 1,np.shape(stratif_data)[1],np.shape(stratif_data)[2],np.shape(stratif_data)[3],len(ref_grid)))121 gridded_stratif_data = -1.*np.ones((np.shape(stratif_data)[0] - 1,np.shape(stratif_data)[1],np.shape(stratif_data)[2],np.shape(stratif_data)[3],len(ref_grid))) 128 122 for iprop in range(1,np.shape(stratif_data)[0]): 129 123 for i in range(np.shape(stratif_data)[1]): 130 124 for j in range(np.shape(stratif_data)[2]): 131 125 for k in range(np.shape(stratif_data)[3]): 132 f = interpolate.interp1d(stratif_data[0][i,j,k,:],stratif_data[iprop][i,j,k,:],kind = 'next',fill_value = "extrapolate") 133 gridded_stratif_data[iprop - 1,i,j,k,:] = f(ref_grid) 126 i_hmax = np.max(np.nonzero(stratif_data[0][i,j,k,:])) 127 hmax = stratif_data[0][i,j,k,i_hmax] 128 i_zmax = np.searchsorted(ref_grid,hmax,'left') 129 f = interpolate.interp1d(stratif_data[0][i,j,k,:i_hmax + 1],stratif_data[iprop][i,j,k,:i_hmax + 1],kind = 'next')#,fill_value = "extrapolate") 130 gridded_stratif_data[iprop - 1,i,j,k,:i_zmax] = f(ref_grid[:i_zmax]) 134 131 135 132 return ref_grid, gridded_stratif_data … … 153 150 154 151 ### Processing 155 stratif_data, max_top_elevation = process_files(folder_path,base_name)152 stratif_data, max_top_elevation, longitude, latitude = process_files(folder_path,base_name) 156 153 ref_grid, gridded_stratif_data = interpolate_data(stratif_data,max_top_elevation,dz) 157 154 date_time = read_infofile(infofile)
Note: See TracChangeset
for help on using the changeset viewer.