Changeset 3460

Timestamp:

Oct 15, 2024, 1:41:48 PM (6 weeks ago)

Author:

jbclement

Message:

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

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/README (modified) (1 diff)
• deftank/visu_evol_layering.py (modified) (6 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -445 +445 @@
 == 11/10/2024 == JBC
 Addition of a python script to visualize the layerings over time and on a reference grid for elevation which is specified by the user.
+
+== 15/10/2024 == JBC
+Update and correction for the visualization of the layerings over time.

trunk/LMDZ.COMMON/libf/evolution/deftank/README

@@ -49 +49 @@
   Bash script file to clean the folder after a PEM simulation.
 
-# output_layering.py:
+# multiple_exec.sh:
+  Bash script to execute multiple scripts in subdirectories, useful to launch multiple simulations at once.
+
+# visu_layering.py:
   Python script file to output the stratification data from the "startpem.nc" files.
+
+# visu_evol_layering.py:
+  Python script file to output the stratification data over time from the "startpem.nc" files.
 
 Note:

trunk/LMDZ.COMMON/libf/evolution/deftank/visu_evol_layering.py

@@ -13 +13 @@
 folder_path = "starts"
 base_name = "restartpem"
-dz = 1 # Discrization step of the reference grid for the elevation [m]
+dz = 0.1 # Discrization step of the reference grid for the elevation [m]
 infofile = 'info_PEM.txt'
 #########################
@@ -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
     nslope = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').dimensions['nslope'].size # nslope is the same for all files
+    longitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['longitude'][:]
+    latitude = Dataset(os.path.join(folder_path,base_name + "1.nc"),'r').variables['latitude'][:]
     for i in range(1,nfile + 1):
         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][:]))
 
-    print(f" > number of files    = {nfile}")
-    print(f" > ngrid              = {ngrid}")
-    print(f" > nslope             = {nslope}")
-    print(f" > max(nb_str_max)    = {max_nb_str}")
-    print(f" > max(top_elevation) = {max_top_elevation}")
+    print(f"> number of files    = {nfile}")
+    print(f"> ngrid              = {ngrid}")
+    print(f"> nslope             = {nslope}")
+    print(f"> max(nb_str_max)    = {max_nb_str}")
+    print(f"> max(top_elevation) = {max_top_elevation}")
 
     # Concatenate stratif variables with dimension 'nb_str_max' along the "Time" dimension
     stratif_data = []
     stratif_heights = np.zeros((ngrid,nfile,nslope,max_nb_str))
-    stratif_co2ice = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))
-    stratif_h2oice = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))
-    stratif_dust = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))
-    stratif_air = -1.*np.ones((ngrid,nfile,nslope,max_nb_str))
+    stratif_co2ice = np.zeros((ngrid,nfile,nslope,max_nb_str))
+    stratif_h2oice = np.zeros((ngrid,nfile,nslope,max_nb_str))
+    stratif_dust = np.zeros((ngrid,nfile,nslope,max_nb_str))
+    stratif_air = np.zeros((ngrid,nfile,nslope,max_nb_str))
     for var_name in datasets[0].variables:
         if 'top_elevation' in var_name:
@@ -108 +110 @@
 
     stratif_data = [stratif_heights,stratif_co2ice,stratif_h2oice,stratif_dust,stratif_air]
-    return stratif_data, max_top_elevation
+    return stratif_data, max_top_elevation, longitude, latitude
 
 ### Function to interpolate the stratification data on a reference grid
@@ -114 +116 @@
     # Define the reference ref_grid
     ref_grid = np.arange(0,max_top_elevation,dz)
-    print(f" > Number of ref_grid points = {len(ref_grid)}")
-
-    # Extrapolate the top_elevation of strata
-    for i in range(np.shape(stratif_data)[1]):
-        for j in range(np.shape(stratif_data)[2]):
-            for k in range(np.shape(stratif_data)[3]):
-                if len(np.where(stratif_data[0][i,j,k,:] == 0)[0]) > 1: # If there are trailing zeros in top_elevation
-                    i0 = np.where(stratif_data[0][i,j,k,:] == 0)[0][1]
-                    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)
+    print(f"> Number of ref_grid points = {len(ref_grid)}")
 
     # Interpolate the strata properties on the ref_grid
-    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)))
+    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)))
     for iprop in range(1,np.shape(stratif_data)[0]):
         for i in range(np.shape(stratif_data)[1]):
             for j in range(np.shape(stratif_data)[2]):
                 for k in range(np.shape(stratif_data)[3]):
-                    f = interpolate.interp1d(stratif_data[0][i,j,k,:],stratif_data[iprop][i,j,k,:],kind = 'next',fill_value = "extrapolate")
-                    gridded_stratif_data[iprop - 1,i,j,k,:] = f(ref_grid)
+                    i_hmax = np.max(np.nonzero(stratif_data[0][i,j,k,:]))
+                    hmax = stratif_data[0][i,j,k,i_hmax]
+                    i_zmax = np.searchsorted(ref_grid,hmax,'left')
+                    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")
+                    gridded_stratif_data[iprop - 1,i,j,k,:i_zmax] = f(ref_grid[:i_zmax])
 
     return ref_grid, gridded_stratif_data
@@ -153 +150 @@
 
 ### Processing
-stratif_data, max_top_elevation = process_files(folder_path,base_name)
+stratif_data, max_top_elevation, longitude, latitude = process_files(folder_path,base_name)
 ref_grid, gridded_stratif_data = interpolate_data(stratif_data,max_top_elevation,dz)
 date_time = read_infofile(infofile)