Ignore:
Timestamp:
Oct 15, 2024, 1:41:48 PM (2 months ago)
Author:
jbclement
Message:

PEM:
Update and correction for the visualization of the layerings over time.
JBC

Location:
trunk/LMDZ.COMMON/libf/evolution
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.COMMON/libf/evolution/changelog.txt

    r3458 r3460  
    445445== 11/10/2024 == JBC
    446446Addition 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
     449Update and correction for the visualization of the layerings over time.
  • trunk/LMDZ.COMMON/libf/evolution/deftank/README

    r3403 r3460  
    4949  Bash script file to clean the folder after a PEM simulation.
    5050
    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:
    5255  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.
    5359
    5460Note:
  • trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

    r3458 r3460  
    1313folder_path = "starts"
    1414base_name = "restartpem"
    15 dz = 1 # Discrization step of the reference grid for the elevation [m]
     15dz = 0.1 # Discrization step of the reference grid for the elevation [m]
    1616infofile = 'info_PEM.txt'
    1717#########################
     
    3939    ngrid = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['physical_points'].size # ngrid is the same for all files
    4040    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'][:]
    4143    for i in range(1,nfile + 1):
    4244        file_path = os.path.join(folder_path,base_name + str(i) + ".nc")
     
    5355            max_top_elevation = max(max_top_elevation,np.max(ds.variables[slope_var_name][:]))
    5456
    55     print(f" > number of files    = {nfile}")
    56     print(f" > ngrid              = {ngrid}")
    57     print(f" > nslope             = {nslope}")
    58     print(f" > max(nb_str_max)    = {max_nb_str}")
    59     print(f" > max(top_elevation) = {max_top_elevation}")
     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}")
    6062
    6163    # Concatenate stratif variables with dimension 'nb_str_max' along the "Time" dimension
    6264    stratif_data = []
    6365    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))
    6870    for var_name in datasets[0].variables:
    6971        if 'top_elevation' in var_name:
     
    108110
    109111    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
    111113
    112114### Function to interpolate the stratification data on a reference grid
     
    114116    # Define the reference ref_grid
    115117    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)}")
    125119
    126120    # 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)))
    128122    for iprop in range(1,np.shape(stratif_data)[0]):
    129123        for i in range(np.shape(stratif_data)[1]):
    130124            for j in range(np.shape(stratif_data)[2]):
    131125                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])
    134131
    135132    return ref_grid, gridded_stratif_data
     
    153150
    154151### Processing
    155 stratif_data, max_top_elevation = process_files(folder_path,base_name)
     152stratif_data, max_top_elevation, longitude, latitude = process_files(folder_path,base_name)
    156153ref_grid, gridded_stratif_data = interpolate_data(stratif_data,max_top_elevation,dz)
    157154date_time = read_infofile(infofile)
Note: See TracChangeset for help on using the changeset viewer.