Remove layer: python script

components/mpas-ocean/src/sea_level_analysis.py

@@ -7,14 +7,24 @@
 conserv_prefix = 'mpaso.hist.am.conservationCheck'
 ts_prefix = 'mpaso.hist.am.timeSeriesStatsMonthly'
 output_dir = '/lcrc/group/acme/ac.cbegeman/ssh_analysis'
+
 decimal_year = []
 total_ssh = []
 ssh_rate = []
 A = []
+V = []
+accumulated_snow_flux = []
+accumulated_rain_flux = []
+accumulated_evaporation_flux = []
+accumulated_sea_ice_flux = []
 accumulated_land_ice_flux = []
+accumulated_frazil_flux = []
 accumulated_iceberg_flux = []
+accumulated_river_runoff_flux = []
+accumulated_ice_runoff_flux = []
 accumulated_removed_river_runoff_flux = []
 accumulated_removed_ice_runoff_flux = []
+
 s_year = 3.154e7 
 s_month = s_year / 12.
 rho_sw = 1026.
@@ -26,31 +36,75 @@
        ds = xr.open_dataset(f'{run_path}/{file_prefix}{conserv_prefix}.{date_suffix}')
        ds_ts = xr.open_dataset(f'{run_path}/{file_prefix}{ts_prefix}.{date_suffix}')
        A.append(ds_ts.timeMonthly_avg_areaCellGlobal.values[0])
-       V = ds_ts.timeMonthly_avg_volumeCellGlobal.values
+       V.append(ds_ts.timeMonthly_avg_volumeCellGlobal.values[0])
+       accumulated_ice_runoff_flux.append(
+           ds.accumulatedIceRunoffFlux.values[0])
+       accumulated_river_runoff_flux.append(
+           ds.accumulatedRiverRunoffFlux.values[0])
        accumulated_removed_ice_runoff_flux.append(
-           ds.accumulatedRemovedIceRunoffFlux.values[0])
+           -ds.accumulatedRemovedIceRunoffFlux.values[0])
        accumulated_removed_river_runoff_flux.append(
-           ds.accumulatedRemovedRiverRunoffFlux.values[0])
+           -ds.accumulatedRemovedRiverRunoffFlux.values[0])
        accumulated_iceberg_flux.append(
            ds.accumulatedIcebergFlux.values[0])
        accumulated_land_ice_flux.append(
            ds.accumulatedLandIceFlux.values[0])
+       accumulated_rain_flux.append(
+           ds.accumulatedRainFlux.values[0])
+       accumulated_snow_flux.append(
+           ds.accumulatedSnowFlux.values[0])
+       accumulated_evaporation_flux.append(
+           ds.accumulatedEvaporationFlux.values[0])
+       accumulated_sea_ice_flux.append(
+           ds.accumulatedSeaIceFlux.values[0])
+       accumulated_frazil_flux.append(
+           ds.accumulatedFrazilFlux.values[0])
        #total_ssh_from_components.append(total_ssh_from_components[-1] + total_flux * s_month / A)
-       total_ssh.append(V / A[-1])
 # print(accumulated_land_ice_flux)
 A = np.array(A)
+V = np.array(V)
+total_ssh = V / A
 
 # SSH rates are given in mm/year
 factor = mass_flux_to_SLR * s_year / A
-ssh_rate_land_ice = np.array(accumulated_land_ice_flux) * factor
+ssh_rate_rain = np.array(accumulated_rain_flux) * factor
+ssh_rate_snow = np.array(accumulated_snow_flux) * factor
+ssh_rate_evaporation = np.array(accumulated_evaporation_flux) * factor
+ssh_rate_sea_ice = np.array(accumulated_sea_ice_flux) * factor
+ssh_rate_river_runoff = np.array(accumulated_river_runoff_flux) * factor
+ssh_rate_ice_runoff = np.array(accumulated_ice_runoff_flux) * factor
 ssh_rate_iceberg = np.array(accumulated_iceberg_flux) * factor
+ssh_rate_frazil = np.array(accumulated_frazil_flux) * factor
+ssh_rate_land_ice = np.array(accumulated_land_ice_flux) * factor
 ssh_rate_removed_river_runoff = np.array(accumulated_removed_river_runoff_flux) * factor
 ssh_rate_removed_ice_runoff = np.array(accumulated_removed_ice_runoff_flux) * factor
-ssh_rate_total = (ssh_rate_land_ice + ssh_rate_iceberg) - \
-                 (ssh_rate_removed_river_runoff + ssh_rate_removed_ice_runoff)
-ssh_anomaly = np.array(total_ssh) - total_ssh[0]
-dt = s_month
-ssh_rate_derived = 1000 * (ssh_anomaly[1:] - ssh_anomaly[:-1]) * s_year / dt # [m] [1/s] * [mm/m] * [s/year]
+
+ssh_rate_total = ssh_rate_rain + \
+                 ssh_rate_snow + \
+                 ssh_rate_evaporation + \
+                 ssh_rate_river_runoff + \
+                 ssh_rate_ice_runoff + \
+                 ssh_rate_sea_ice + \
+                 ssh_rate_iceberg + \
+                 ssh_rate_frazil + \
+                 ssh_rate_land_ice
+
+ssh_anomaly = (total_ssh - total_ssh[0]) * 1000. # convert from m to mm
+dt = 1 / 12 # in years
+ssh_rate_derived = (ssh_anomaly[1:] - ssh_anomaly[:-1]) / dt # [m] [1/s] * [mm/m] * [s/year]
+ssh_anomaly_cum = np.zeros((len(ssh_rate_total)))
+ssh_start = 0
+for i in range(len(ssh_anomaly_cum)):
+    ssh_anomaly_cum[i] = ssh_start + ssh_rate_total[i] / 12
+    ssh_start = ssh_anomaly_cum[i]
+print(ssh_anomaly_cum)
+
+#print(f'mean(A) = {np.mean(A)}')
+#print(f'mean(V) = {np.mean(V)}')
+#print(f'mean(total_ssh) = {np.mean(total_ssh)}')
+#print(f'mean(ssh_rate_total) = {np.mean(ssh_rate_total)}')
+#print(f'min,max(ssh_anomaly) = {np.min(ssh_anomaly)},{np.max(ssh_anomaly)}')
+#print(f'min,max(ssh_anomaly_cum) = {np.min(ssh_anomaly_cum)},{np.max(ssh_anomaly_cum)}')
 
 fig = plt.figure()
 plt.plot(decimal_year, accumulated_land_ice_flux, label='land ice')
@@ -61,8 +115,11 @@
 fig.savefig(f'{output_dir}/ssh_flux_components.png')
 
 fig = plt.figure()
-plt.plot(decimal_year[1:], ssh_rate_derived[:, 0], label='from volume')
-plt.plot(decimal_year, ssh_rate_total, label='total')
+plt.plot(decimal_year[1:], ssh_rate_derived[:], label='from volume')
+plt.plot(decimal_year, ssh_rate_total, label='from flux')
+fig.savefig(f'{output_dir}/ssh_rate_total.png')
+
+fig = plt.figure()
 plt.plot(decimal_year, ssh_rate_land_ice, label='land ice')
 plt.plot(decimal_year, ssh_rate_iceberg, label='iceberg')
 plt.plot(decimal_year, ssh_rate_removed_river_runoff, label='removed river runoff')
@@ -71,5 +128,7 @@
 fig.savefig(f'{output_dir}/ssh_rate_components.png')
 
 fig = plt.figure()
-plt.plot(decimal_year, ssh_anomaly)
+plt.plot(decimal_year, ssh_anomaly, label='volume-derived')
+plt.plot(decimal_year, ssh_anomaly_cum - ssh_anomaly_cum[0], label='flux-derived')
+plt.legend()
 fig.savefig(f'{output_dir}/ssh_time_series.png')