Changeset 3883

Timestamp:

Aug 8, 2025, 7:50:27 PM (19 hours ago)

Author:

jbclement

Message:

PEM:

• Replacing the Bash script "concat_diagpem.sh" by a Python script "concat_pem.py" which is much faster.
• Corrections and updates for "visu_evol_layering.py".

JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• deftank/README (modified) (1 diff)
• deftank/concat_diagpem.sh (deleted)
• deftank/concat_pem.py (added)
• deftank/visu_evol_layering.py (modified) (14 diffs)

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

@@ -766 +766 @@
 == 28/07/2025 == JBC
 Bug correction to detect the job time limit with PBS/TORQUE. Making it more robust and automatic for the launching script and the Fortran code.
+
+== 08/08/2025 == JBC
+- Replacing the Bash script "concat_diagpem.sh" by a Python script "concat_pem.py" which is much faster.
+- Corrections and updates for "visu_evol_layering.py".

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

@@ -103 +103 @@
   Bash script file to set the orbital parameters of a file "startfi.nc" from Laskar's data contained in "obl_ecc_lsp.asc" according to the initial date 'year_earth_bp_ini' defined in "run_PEM.def". See also "modify_startfi_orbit.sh".
 
-# concat_diagpem.sh:
-  Bash script file to concatenate along the variable 'Time' all the "diagpem.nc" files of the PEM into one NetCDF file. 'Time' is re-indexed with the numbering of Martian years simulated by the PEM run.
+# concat_pem.py:
+  Python script file to concatenate the NetCDF files of the PEM along the dimension 'Time' into one NetCDF file. 'Time' is re-indexed by increment.
 
 # clean.sh:

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

@@ -509 +509 @@
                 output_folder, f"layering_evolution_ig{ig+1}_is{isl+1}.png"
             )
-            fig.savefig(fname, dpi=150)
+            fig.savefig(fname, dpi=1200, bbox_inches='tight')
 
 
@@ -577 +577 @@
             rgb = np.ones((nz, ntime, 3), dtype=float)
 
-            frac_all = np.zeros((nz, ntime, 3), dtype=float)  # store fH2O, fCO2, fDust
+            frac_all = np.full((nz, ntime, 3), np.nan, dtype=float)  # store fH2O, fCO2, fDust
             for t in range(ntime):
                 depth = ti[ig, t, isl]
@@ -598 +598 @@
             display_frac = frac_all[elevation_mask, :, :]
 
-            display_rgb = rgb[elevation_mask, :, :]
-
             # Compute edges for pcolormesh
             dt = date_time[1] - date_time[0] if len(date_time) > 1 else 1
@@ -609 +607 @@
             # Create figure with legend
             fig, (ax_main, ax_leg) = plt.subplots(
-                1, 2, figsize=(8, 4), dpi=200,
+                1, 2, figsize=(8, 4), dpi=150,
                 gridspec_kw={'width_ratios': [5, 1]}
             )
@@ -630 +628 @@
                 i = np.searchsorted(x_edges, x) - 1
                 j = np.searchsorted(y_edges, y) - 1
-                i = np.clip(i, 0, display_rgb.shape[1]-1)
-                j = np.clip(j, 0, display_rgb.shape[0]-1)
+                i = np.clip(i, 0, display_rgb.shape[1] - 1)
+                j = np.clip(j, 0, display_rgb.shape[0] - 1)
                 # get fractions
                 fH2O, fCO2, fDust = display_frac[j, i]
-                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.2f}, CO2={fCO2:.2f}, Dust={fDust:.2f}"
+                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.4f}, CO2={fCO2:.4f}, Dust={fDust:.4f}"
             ax_main.format_coord = main_format
             ax_main.set_facecolor('white')
@@ -684 +682 @@
             plt.tight_layout()
             fname = os.path.join(output_folder, f"layering_rgb_evolution_ig{ig+1}_is{isl+1}.png")
-            fig.savefig(fname, dpi=150, bbox_inches='tight')
+            fig.savefig(fname, dpi=1200, bbox_inches='tight')
 
 
@@ -722 +720 @@
     vmin, vmax = -2, 1
     epsilon = 1e-6
-    
+
     # Loop over grids and slopes
     for ig in range(ngrid):
         for isl in range(nslope):
             ti = top_index[ig, :, isl].copy().astype(int)
-            for t in range(1, ntime):
-                if ti[t] <= 0:
-                    ti[t] = ti[t-1]
 
             # Compute log10(dust/ice) profile at each time step
             log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
             for t in range(ntime):
+                if ti[t] <= 0:
+                    ti[t] = ti[t-1]
                 zmax = ti[t]
                 if zmax <= 0:
                     continue
 
-                h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None)
-                dust_profile = np.clip(dust[ig, t, isl, :zmax], 0, None)
+                cH2O = np.clip(h2o[ig, t, isl, :zmax], 0, None)
+                cCO2 = np.clip(co2[ig, t, isl, :zmax], 0, None)
+                cDust = np.clip(dust[ig, t, isl, :zmax], 0, None)
 
                 with np.errstate(divide='ignore', invalid='ignore'):
-                    ratio_profile = np.where(
-                        h2o_profile > 0,
-                        dust_profile / h2o_profile,
+                    ratio = np.where(
+                        cH2O > 0,
+                        cDust / cH2O,
                         10**(vmax + 1)
                     )
-                    log_ratio = np.log10(ratio_profile + epsilon)
+                    log_ratio = np.log10(ratio + epsilon)
                     log_ratio = np.clip(log_ratio, vmin, vmax)
 
                 log_ratio_array[:zmax, t] = log_ratio
 
-            # Convert back to linear ratio and apply elevation mask
+            # Mask and prepare for display
             ratio_array = 10**log_ratio_array
             ratio_display = ratio_array[elevation_mask, :]
+            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1] - date_time[-2])]])
+            y_edges = np.concatenate([elev, [elev[-1] + (elev[-1] - elev[-2])]])
 
             # Plot
@@ -766 +766 @@
                 norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
             )
-            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]])
-            attach_format_coord(ax, ratio_display, x_edges, np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]]), is_pcolormesh=True)
+            attach_format_coord(ax, ratio_display, x_edges, y_edges, is_pcolormesh=True)
             ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
             ax.set_xlabel('Time (Mars years)')
             ax.set_ylabel('Elevation (m)')
 
-            # Add colorbar with simplified ratio labels
-            cbar = fig.colorbar(im, ax=ax, orientation='vertical')
+            # Add colorbar
+            cbar = fig.colorbar(im, ax=ax, orientation='vertical', pad=0.15)
             cbar.set_label('Dust / H₂O ice (ratio)')
-
-            # Define custom ticks and labels
-            ticks = [1e-2, 1e-1, 1, 1e1]
-            labels = ['1:100', '1:10', '1:1', '10:1']
-            cbar.set_ticks(ticks)
-            cbar.set_ticklabels(labels)
+            cbar.set_ticks([1e-2, 1e-1, 1, 1e1])
+            cbar.set_ticklabels(['1:100', '1:10', '1:1', '10:1'])
 
             # Save figure
             plt.tight_layout()
-            fname = os.path.join(
+            outname = os.path.join(
                 output_folder,
                 f"dust_to_ice_ratio_grid{ig+1}_slope{isl+1}.png"
             )
-            fig.savefig(fname, dpi=150)
+            fig.savefig(outname, dpi=1200, bbox_inches='tight')
 
 
@@ -1060 +1055 @@
         elev = ref_grid
 
-    # Custom colormap: blue (ice) to orange (dust)
-    blue = np.array([0, 0, 255]) / 255
-    orange = np.array([255, 165, 0]) / 255
+    # Define custom blue-to-orange colormap
+    blue = np.array([0, 0, 255], dtype=float) / 255
+    orange = np.array([255, 165, 0], dtype=float) / 255
     custom_cmap = LinearSegmentedColormap.from_list('BlueOrange', [blue, orange], N=256)
 
@@ -1072 +1067 @@
     for ig in range(ngrid):
         for isl in range(nslope):
+            # Prepare fraction and ratio arrays
             ti = top_index[ig, :, isl].copy().astype(int)
-            frac_all = np.zeros((nz, ntime, 3), dtype=float)  # store fH2O, fCO2, fDust
-            for t in range(1, ntime):
+            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
+            frac_all = np.full((nz, ntime, 3), np.nan, dtype=float)  # store fH2O, fCO2, fDust
+            for t in range(ntime):
                 if ti[t] <= 0:
                     ti[t] = ti[t-1]
-
-            # Compute log10(dust/ice) profile at each time step
-            log_ratio_array = np.full((nz, ntime), np.nan, dtype=np.float32)
-            for t in range(ntime):
                 zmax = ti[t]
                 if zmax <= 0:
                     continue
+
                 cH2O = np.clip(h2o[ig, t, isl, :zmax], 0, None)
                 cCO2 = np.clip(co2[ig, t, isl, :zmax], 0, None)
@@ -1094 +1088 @@
                 frac_all[:zmax, t, :] = np.stack([fH2O, fCO2, fDust], axis=1)
 
-                h2o_profile = np.clip(h2o[ig, t, isl, :zmax], 0, None)
-                dust_profile = np.clip(dust[ig, t, isl, :zmax], 0, None)
-
                 with np.errstate(divide='ignore', invalid='ignore'):
-                    ratio_profile = np.where(
-                        h2o_profile > 0,
-                        dust_profile / h2o_profile,
+                    ratio = np.where(
+                        cH2O > 0,
+                        cDust / cH2O,
                         10**(vmax + 1)
                     )
-                    log_ratio = np.log10(ratio_profile + epsilon)
+                    log_ratio = np.log10(ratio + epsilon)
                     log_ratio = np.clip(log_ratio, vmin, vmax)
 
                 log_ratio_array[:zmax, t] = log_ratio
 
-            # Convert back to linear ratio and mask
+            # Mask elevation
             ratio_array = 10**log_ratio_array
             ratio_display = ratio_array[elevation_mask, :]
             display_frac = frac_all[elevation_mask, :, :]
 
+            # Compute edges for pcolormesh
+            dt = date_time[1] - date_time[0] if len(date_time) > 1 else 1
+            x_edges = np.concatenate([date_time, [date_time[-1] + dt]])
+            d_e = np.diff(elev)
+            last_e = elev[-1] + (d_e[-1] if len(d_e)>0 else 1)
+            y_edges = np.concatenate([elev, [last_e]])
+
             # Plot
             fig, ax = plt.subplots(figsize=(8, 6), dpi=150)
             im = ax.pcolormesh(
-                date_time,
-                elev,
+                x_edges,
+                y_edges,
                 ratio_display,
                 shading='auto',
@@ -1123 +1121 @@
                 norm=LogNorm(vmin=10**vmin, vmax=10**vmax),
             )
-            x_edges = np.concatenate([date_time, [date_time[-1] + (date_time[-1]-date_time[-2])]])
-            y_edges = np.concatenate([elev, [elev[-1] + (elev[-1]-elev[-2])]])
-            def format_coord_all(x, y):
+
+            def format_coord_custom(x_input, y_input):
+                # map onto the main axis
+                if plt.gca() is ax2:
+                    x_pix, y_pix = ax2.transData.transform((x_input, y_input))
+                    x, y = ax.transData.inverted().transform((x_pix, y_pix))
+                else:
+                    x, y = x_input, y_input
                 # check bounds
                 if x < x_edges[0] or x > x_edges[-1] or y < y_edges[0] or y > y_edges[-1]:
@@ -1132 +1135 @@
                 i = np.searchsorted(x_edges, x) - 1
                 j = np.searchsorted(y_edges, y) - 1
-                i = np.clip(i, 0, display_frac.shape[1]-1)
-                j = np.clip(j, 0, display_frac.shape[0]-1)
-                # get fractions
-                fH2O  = display_frac[j, i, 0]
-                fDust = display_frac[j, i, 2]
+                i = np.clip(i, 0, ratio_display.shape[1] - 1)
+                j = np.clip(j, 0, ratio_display.shape[0] - 1)
+                # get fractions and obliquity
+                fH2O, fCO2, fDust = display_frac[j, i]
                 obl   = np.interp(x, date_time, obliquity)
-                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.2f}, Dust={fDust:.2f}, Obl={obl:.2f}°"
-
-            ax.format_coord = format_coord_all
+                return f"Time={x:.2f}, Elev={y:.2f}, H2O={fH2O:.4f}, Dust={fDust:.4f}, Obl={obl:.2f}°"
+
             ax.set_title(f"Dust-to-Ice ratio over time (Grid point {ig+1}, Slope {isl+1})", fontweight='bold')
             ax.set_xlabel('Time (Mars years)')
@@ -1161 +1162 @@
                     linewidth=1.5
                 )
-            ax2.format_coord = format_coord_all
+            ax2.format_coord = format_coord_custom
             ax2.set_ylabel('Obliquity (°)')
             ax2.tick_params(axis='y')
             ax2.grid(False)
 
-            # Save
-            os.makedirs(output_folder, exist_ok=True)
+            # Save figure
+            plt.tight_layout()
             outname = os.path.join(
                 output_folder,
                 f'dust_ice_obliquity_grid{ig+1}_slope{isl+1}.png'
             )
-            plt.tight_layout()
-            fig.savefig(outname, dpi=150)
+            fig.savefig(outname, dpi=1200, bbox_inches='tight')