diff --git a/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O1.png b/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O1.png new file mode 100644 index 000000000..0ea434baf Binary files /dev/null and b/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O1.png differ diff --git a/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O2.png b/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O2.png new file mode 100644 index 000000000..66878d31e Binary files /dev/null and b/diagnostics/top_heaviness_metric/doc/Long_term_mean_of_O2.png differ diff --git a/diagnostics/top_heaviness_metric/doc/Q1&Q2_R.png b/diagnostics/top_heaviness_metric/doc/Q1&Q2_R.png new file mode 100644 index 000000000..405c4d5bc Binary files /dev/null and b/diagnostics/top_heaviness_metric/doc/Q1&Q2_R.png differ diff --git a/diagnostics/top_heaviness_metric/doc/R2_Between_Recon_Omega&Original.png b/diagnostics/top_heaviness_metric/doc/R2_Between_Recon_Omega&Original.png new file mode 100644 index 000000000..870568431 Binary files /dev/null and b/diagnostics/top_heaviness_metric/doc/R2_Between_Recon_Omega&Original.png differ diff --git a/diagnostics/top_heaviness_metric/doc/Top_Heaviness_Ratio.png b/diagnostics/top_heaviness_metric/doc/Top_Heaviness_Ratio.png new file mode 100644 index 000000000..a756e5636 Binary files /dev/null and b/diagnostics/top_heaviness_metric/doc/Top_Heaviness_Ratio.png differ diff --git a/diagnostics/top_heaviness_metric/doc/top_heaviness_metric.rst b/diagnostics/top_heaviness_metric/doc/top_heaviness_metric.rst new file mode 100644 index 000000000..c45b4d1cd --- /dev/null +++ b/diagnostics/top_heaviness_metric/doc/top_heaviness_metric.rst @@ -0,0 +1,91 @@ +Top-Heaviness Metric Diagnostic Documentation +================================ + +Last update: 5/30/2021 + +The vertical profiles of diabatic heating have important implications for large-scale dynamics, especially for the coupling between the large-scale atmospheric circulation and precipitation processes. We adopt an objective approach to examine the top-heaviness of vertical motion (Back et al. 2017), which is closely related to the heating profiles and a commonly available model output variable. The diagnostic metric can also be used to evaluate the diabatic heating profile. + +Version & Contact info +---------------------- + +.. '-' starts items in a bulleted list: + https://docutils.sourceforge.io/docs/user/rst/quickref.html#bullet-lists + +- Version/revision information: version 1.0 (6/28/2021) +- Developer/point of contact (Jiacheng Ye, jye18@illinois.edu, DAS UIUC; Zhuo Wang, zhuowang@illinois.edu, DAS UIUC) + +.. Underline with '^'s to make a third-level heading. + +Open source copyright agreement +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The MDTF framework is distributed under the LGPLv3 license (see LICENSE.txt). + +Functionality +------------- + +The current package consists of following functionalities: + +(1) Calculation of the fractional variance of vertical velocity at each grid point explained by two base functions, Q1 (~idealized deep convection profile) and Q2 (~idealized deep stratiform profile) + +(2) Calculation of the top-heaviness ratio (O2/O1) + +As a module of the MDTF code package, all scripts of this package can be found under +``mdtf/MDTF_$ver/diagnostics/top_heaviness_ratio`` + +Required programming language and libraries +------------------------------------------- + +Python3 packages: "netCDF4", "xarray", "numpy", "scipy", "matplotlib", "cartopy" + +Required model output variables +------------------------------- + +1) wap (plev x lat x lon) : Vertical Velocity [Pa/s], which can be either the monthly mean for a certain year or the long-term monthly/season mean. + +References +---------- + +.. : + +Back, L. E., Hansen, Z., & Handlos, Z. (2017). Estimating vertical motion profile top-heaviness: Reanalysis compared to satellite-based observations and stratiform rain fraction. Journal of the Atmospheric Sciences, 74(3), 855-864. https://doi.org/10.1175/JAS-D-16-0062.1 + +Jiacheng and Zhuo's paper is in preparation. + +More about this diagnostic +-------------------------- + +Q1 and Q2 (Figure 1a) are two prescribed base functions. Following Back et al. (2017), Q1 as a half sine function, and Q2 as a full sine function, which represent the idealized deep convection profile and the idealized stratiform profile, respectively. The vertical velocity can be approximated by Q1 and Q2: + + ω'(x,y,p) = O1(x,y) * Q1(p) + O2(x,y)*Q2(p) + + +where O1 and O2 are the coefficients of Q1 and Q2, respectively. Back et al. (2017) showed that Q1 and Q2 resemble the first two EOF modes of vertical velocity profile variability. Then O1 and O2 can be approximately regarded as the corresponding principal component time series. + +For O1>0, ω' transitions from a bottom-heavy profile to a top-heavy profile when the ratio of r=O2/O1 increases from -1 to 1 (Figure 1b). The ratio r is thus defined as the top-heaviness ratio. + +To assess how well ω' approximates ω, the fractional variance is calculated over each grid point. The fractional variance is defined as the square of the pearson correlation between ω' and ω. As shown in Figure 2, ω' explains more than 80% of the vertical variance over most tropical/subtropical oceanic grid points. + +The top-heaviness ratio (r) is presented in Figure 3. The Western Pacific is dominated by more top-heavy vertical profiles while the Eastern Pacific and Atlantic are characterized by more bottom-heavy profiles, exhibiting a great contrast. + + + +.. figure:: Q1&Q2_R.png + :align: center + :width: 100 % + + Figure 1. (a) Q1 and Q2; (b) Vertical velocity profiles constructed from the varying top-heaviness ratio (r; r=-1: dark blue, r=1: dark red). + + +.. figure:: R2_Between_Recon_Omega&Original.png + :align: center + :width: 100 % + + Figure 2. The fractional variance of ω explained by ω'. + + +.. figure:: Top_Heaviness_Ratio.png + :align: center + :width: 100 % + + Figure 3. Long-term mean Top-Heaviness Ratio in July (2000-2019). The ratio is presented only for the grid points with O1 no less than 0.01. diff --git a/diagnostics/top_heaviness_metric/settings.jsonc b/diagnostics/top_heaviness_metric/settings.jsonc new file mode 100644 index 000000000..cd8b86918 --- /dev/null +++ b/diagnostics/top_heaviness_metric/settings.jsonc @@ -0,0 +1,34 @@ +{ + "settings": { + "driver": "top_heaviness_metric.py", + "long_name": "Top Heaviness Metric Diagnostic", + "realm": "atmos", + "description": "The vertical profiles of diabatic heating have important implications for large-scale dynamics, especially for the coupling between the large-scale atmospheric circulation and precipitation processes. We adopt an objective approach to examine the top-heaviness of vertical motion, which is closely related to the heating profiles and a commonly available model output variable. The diagnostic/metric can also be used to evaluate the top-heaviness of diabatic heating.", + "runtime_requirements": { + "python3": ["netCDF4", "xarray", "numpy", "scipy", "matplotlib", "cartopy"] + } + }, + "dimensions": { + "lat": { + "standard_name": "latitude" + }, + "lon": { + "standard_name": "longitude" + }, + "lev": { + "standard_name": "air_pressure", + "units": "hPa", + "positive": "down", + "axis": "Z" + } + }, + "varlist": { + "omega": { + "standard_name": "lagrangian_tendency_of_air_pressure", + "path_variable": "PR_FILE", + "units": "Pa s-1", + "dimensions": ["lev", "lat", "lon"] + } + } +} + diff --git a/diagnostics/top_heaviness_metric/top_heaviness_metric.html b/diagnostics/top_heaviness_metric/top_heaviness_metric.html new file mode 100644 index 000000000..e0f343900 --- /dev/null +++ b/diagnostics/top_heaviness_metric/top_heaviness_metric.html @@ -0,0 +1,30 @@ +MDTF example diagnostic + + +

Top-Heaviness Metric Diagnostics

+

+The vertical profiles of diabatic heating have important implications for large-scale dynamics, +especially for the coupling between the large-scale atmospheric circulation and precipitation +processes. We adopt an objective approach to examine the top-heaviness of vertical motion, which +is closely related to the heating profiles and a commonly available model output variable. +The diagnostic/metric can also be used to evaluate the top-heaviness of diabatic heating. +

+ + + + + + +
Top-Heaviness Ratio +plot +plot +
O1 +plot +plot +
O2 +plot +plot +
R2 +plot +plot +
diff --git a/diagnostics/top_heaviness_metric/top_heaviness_metric.py b/diagnostics/top_heaviness_metric/top_heaviness_metric.py new file mode 100644 index 000000000..fd233f01c --- /dev/null +++ b/diagnostics/top_heaviness_metric/top_heaviness_metric.py @@ -0,0 +1,75 @@ +# 28 June top_heaviness_metric.py +# Top-Heaviness Metric +# +# ================================================================================ +# +# Last update: 20 May, 2021 +# Contributors: Jiacheng Ye (jye18@illinois.edu; UIUC), Zhuo Wang (zhuowang@illinois.edu; UIUC) +# +# Evaluate model performance for the representation of vertical motion (omega) profile; +# The diabatic heating profile is closely related to vertical motion profile. Thus, diagnosing omega vertical structure +# would help us to better understand the coupling between large-scale circulation and precipitation process +# +# Version and contact info +# +# - Version: 1.0 +# - Contact info: Jiacheng Ye (jye18@illinois.edu) and +# Zhuo Wang (zhuowang@illinois.edu) +# +# ================================================================================ +# Functionality +# 1) calculate the coefficient of Q1 and Q2 (Q1 ~= idealized deep convection profile; Q2 ~= idealized deep stratiform profile); +# 2) calculate top-heaviness ratio (defined as O2/O1) +# +# ================================================================================ +# +# All scripts of this package can be found under: /diagnostics/top_heaviness_metric +# & observational data under: /obs_data/top_heaviness_metric +# +# Monthly 3-D (time-lat-lon) vertical motion (wap) fields are required; +# +# Required programming language and libraries: Tested in the Python 3.7 environment; +# Required Python libraries: Numpy, Scipy +# +# ================================================================================ +# Reference: +# 1) Back, L. E., Hansen, Z., & Handlos, Z. (2017). Estimating vertical motion profile top-heaviness: +# Reanalysis compared to satellite-based observations and stratiform rain fraction. +# Journal of the Atmospheric Sciences, 74(3), 855-864. +# 2) Our paper which focuses on GEFS v12 diagnostics is under progress... +# + +# driver file +import os +import glob + +missing_file=0 +if len(glob.glob(os.environ["OMEGA_FILE"]))==0: + print("Required OMEGA data missing!") + missing_file=1 + +if missing_file==1: + print("Top-heaviness metric diagnostics Package will NOT be executed!") +else: + try: + os.system("python3 "+os.environ["POD_HOME"]+"/"+"top_heaviness_ratio_calculation.py") + except OSError as e: + print('WARNING',e.errno,e.strerror) + print("**************************************************") + print("Top-Heaviness Metric Diagnostics (top_heaviness_ratio_calculation.py) is NOT Executed as Expected!") + print("**************************************************") + # if the user only focuses on calculating top-heaviess ratio instead of applying some tests on + # ratio robustness, the user can choose not to run the following python file. + try: + os.system("python3 "+os.environ["POD_HOME"]+"/"+"top_heaviness_ratio_robustness_calc.py") + except OSError as e: + print('WARNING',e.errno,e.strerror) + print("**************************************************") + print("Top-Heaviness Metric Diagnostics (top_heaviness_ratio_robustness_calc.py) is NOT Executed as Expected!") + print("**************************************************") + + print("**************************************************") + print("Top-Heaviness Metric Diagnostics Executed!") + print("**************************************************") + + diff --git a/diagnostics/top_heaviness_metric/top_heaviness_ratio_calculation.py b/diagnostics/top_heaviness_metric/top_heaviness_ratio_calculation.py new file mode 100644 index 000000000..306a00bb7 --- /dev/null +++ b/diagnostics/top_heaviness_metric/top_heaviness_ratio_calculation.py @@ -0,0 +1,246 @@ +# 28 June top_heaviness_ratio_calculation.py +#Python packages/ modules imported for the diagnostic +# the sample monthly data is from ERA5 in July from 2000 to 2019 +import os +import xarray as xr +import numpy as np +from scipy import integrate +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +import cartopy.mpl.ticker as cticker + + +#Setting variables equal to environment variables set by the diagnostic +lat_coord = os.environ["lat_coord"] +lon_coord = os.environ["lon_coord"] +lev_coord = os.environ["lev_coord"] +omega_var = os.environ["omega_var"] +WK_DIR = os.environ["WK_DIR"] +OBS_DATA = os.environ["OBS_DATA"] +CASENAME = os.environ["CASENAME"] + +#====================== deriving model output ======================= +def top_heaviness_ratio_calculation_model(reanalysis_path, reanalysis_var): + # read omega and land-sea fraction data + ds= xr.open_dataset(reanalysis_path) + lev_model=ds[lev_coord].values + lat_model=ds[lat_coord].values + lon_model=ds[lon_coord].values + isort=np.argsort(lev_model)[::-1] # descending + mid_omega_model=ds[reanalysis_var].values # mon x lev x lat x lon; for the sample data (July over 2000-2019) + mid_omega_model=mid_omega_model[:,isort] + #======================deriving O1 and O2======================= + # construct Q1_model as half a sine wave and Q2_model as a full sine wave + # two base functions; Q1: idealized deep convection profile; Q2: Deep stratiform profile (Back et al. 2017) + Q1_model=np.zeros(len(lev_model)) + Q2_model=np.zeros(len(lev_model)) + dp=lev_model[-1]-lev_model[0] + for i in range(len(lev_model)): + Q1_model[i]=-np.sin(np.pi*(lev_model[i]-lev_model[0])/dp) + Q2_model[i]=np.sin(2*np.pi*(lev_model[i]-lev_model[0])/dp) + #Normalize + factor=integrate.trapz(Q1_model*Q1_model,lev_model)/dp + Q1_model=Q1_model/np.sqrt(factor) + factor=integrate.trapz(Q2_model*Q2_model,lev_model)/dp + Q2_model=Q2_model/np.sqrt(factor) + # deriving O1 and O2; O1 and O2 are coefs of Q1 and Q2 + mid_omega_model_ltm=np.nanmean(mid_omega_model,axis=0) + O1_model=integrate.trapz(mid_omega_model_ltm*Q1_model[:,None,None],lev_model,axis=0)/dp + O2_model=integrate.trapz(mid_omega_model_ltm*Q2_model[:,None,None],lev_model,axis=0)/dp + #====================== set up figures ======================= + #====================== O1 ======================= + fig, axes = plt.subplots(figsize=(8,4)) + ilat=np.argsort(lat_model) + ilon=np.argsort(lon_model) + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.06,0.07,0.01) + im0 = axes.contourf(lon_model, lat_model, O1_model, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_model[ilon][0],lon_model[ilon][-1],60) + lat_grid = np.arange(lat_model[ilat][0],lat_model[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('O1 [Pa/s]',loc='center',fontsize=16) + fig.tight_layout() + fig.savefig(WK_DIR+"/model/"+CASENAME+"_O1.png", format='png',bbox_inches='tight') + #====================== O2 ======================= + fig, axes = plt.subplots(figsize=(8,4)) + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.06,0.07,0.01) + im0 = axes.contourf(lon_model, lat_model, O2_model, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_model[ilon][0],lon_model[ilon][-1],60) + lat_grid = np.arange(lat_model[ilat][0],lat_model[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('O2 [Pa/s]',loc='center',fontsize=16) + fig.tight_layout() + fig.savefig(WK_DIR+"/model/"+CASENAME+"_O2.png", format='png',bbox_inches='tight') + #====================== O2/O1 top-heaviness ratio ======================= + fig, axes = plt.subplots(figsize=(8,4)) + mmid1=O2_model/O1_model + midi=O1_model<0.01 # We only investigate areas with O1 larger than zero + mmid1[midi]=np.nan + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.6,0.7,0.1) + im0 = axes.contourf(lon_model, lat_model, mmid1, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_model[ilon][0],lon_model[ilon][-1],60) + lat_grid = np.arange(lat_model[ilat][0],lat_model[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('Top-heaviness Ratio (O2/O1)',loc='center',fontsize=18) + fig.tight_layout() + fig.savefig(WK_DIR+"/model/"+CASENAME+"_Top_Heaviness_Ratio.png", format='png',bbox_inches='tight') + print("Plotting Completed") + + +top_heaviness_ratio_calculation_model(os.environ["OMEGA_FILE"],os.environ["omega_var"]) + + + +#====================== deriving obs output ======================= +# run obs data (ERA5 2000-2019 July only) + +def top_heaviness_ratio_calculation_obs(obs_data_full_dir): + # read omega + ds= xr.open_dataset(obs_data_full_dir) + lev_obs=ds['lev'].values + lat_obs=ds['lat'].values + lon_obs=ds['lon'].values + isort=np.argsort(lev_obs)[::-1] # descending + mid_omega_obs=ds['omega'].values # mon x lev x lat x lon; for the sample data (July over 2000-2019) + mid_omega_obs=mid_omega_obs[:,isort] + #======================deriving O1 and O2======================= + # construct Q1_obs as half a sine wave and Q2_obs as a full sine wave + # two base functions; Q1: idealized deep convection profile; Q2: Deep stratiform profile (Back et al. 2017) + Q1_obs=np.zeros(len(lev_obs)) + Q2_obs=np.zeros(len(lev_obs)) + dp=lev_obs[-1]-lev_obs[0] + for i in range(len(lev_obs)): + Q1_obs[i]=-np.sin(np.pi*(lev_obs[i]-lev_obs[0])/dp) + Q2_obs[i]=np.sin(2*np.pi*(lev_obs[i]-lev_obs[0])/dp) + #Normalize + factor=integrate.trapz(Q1_obs*Q1_obs,lev_obs)/dp + Q1_obs=Q1_obs/np.sqrt(factor) + factor=integrate.trapz(Q2_obs*Q2_obs,lev_obs)/dp + Q2_obs=Q2_obs/np.sqrt(factor) + # deriving O1 and O2; O1 and O2 are coefs of Q1 and Q2 + mid_omega_obs_ltm=np.nanmean(mid_omega_obs,axis=0) + O1_obs=integrate.trapz(mid_omega_obs_ltm*Q1_obs[:,None,None],lev_obs,axis=0)/dp + O2_obs=integrate.trapz(mid_omega_obs_ltm*Q2_obs[:,None,None],lev_obs,axis=0)/dp + #====================== set up figures ======================= + #====================== O1 ======================= + fig, axes = plt.subplots(figsize=(8,4)) + ilat=np.argsort(lat_obs) + ilon=np.argsort(lon_obs) + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.06,0.07,0.01) + im0 = axes.contourf(lon_obs, lat_obs, O1_obs, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_obs[ilon][0],lon_obs[ilon][-1],60) + lat_grid = np.arange(lat_obs[ilat][0],lat_obs[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('O1 [Pa/s]',loc='center',fontsize=16) + fig.tight_layout() + fig.savefig(WK_DIR+"/obs/ERA5_O1_2000_2019_July.png", format='png',bbox_inches='tight') + #====================== O2 ======================= + fig, axes = plt.subplots(figsize=(8,4)) + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.06,0.07,0.01) + im0 = axes.contourf(lon_obs, lat_obs, O2_obs, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_obs[ilon][0],lon_obs[ilon][-1],60) + lat_grid = np.arange(lat_obs[ilat][0],lat_obs[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('O2 [Pa/s]',loc='center',fontsize=16) + fig.tight_layout() + fig.savefig(WK_DIR+"/obs/ERA5_O2_2000_2019_July.png", format='png',bbox_inches='tight') + #====================== O2/O1 top-heaviness ratio ======================= + fig, axes = plt.subplots(figsize=(8,4)) + mmid1=O2_obs/O1_obs + midi=O1_obs<0.01 # We only investigate areas with O1 larger than zero + mmid1[midi]=np.nan + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(-0.6,0.7,0.1) + im0 = axes.contourf(lon_obs, lat_obs, mmid1, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_obs[ilon][0],lon_obs[ilon][-1],60) + lat_grid = np.arange(lat_obs[ilat][0],lat_obs[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('Top-heaviness Ratio (O2/O1)',loc='center',fontsize=18) + fig.tight_layout() + fig.savefig(WK_DIR+"/obs/ERA5_Top_Heaviness_Ratio_2000_2019_July.png", format='png',bbox_inches='tight') + print("Plotting Completed") + +top_heaviness_ratio_calculation_obs(OBS_DATA+'/ERA5_omega_mon_2000_2019_July.nc') + diff --git a/diagnostics/top_heaviness_metric/top_heaviness_ratio_robustness_calc.py b/diagnostics/top_heaviness_metric/top_heaviness_ratio_robustness_calc.py new file mode 100644 index 000000000..a8e12e783 --- /dev/null +++ b/diagnostics/top_heaviness_metric/top_heaviness_ratio_robustness_calc.py @@ -0,0 +1,198 @@ +# Last update: 28 June top_heaviness_ratio_robustness_calc.py +import os +import xarray as xr +import numpy as np +import scipy +from scipy import interpolate +from scipy import integrate +import matplotlib.pyplot as plt +import cartopy.crs as ccrs +import cartopy.mpl.ticker as cticker + + +#Setting variables equal to environment variables set by the diagnostic +time_coord = os.environ["time_coord"] +lat_coord = os.environ["lat_coord"] +lon_coord = os.environ["lon_coord"] +lev_coord = os.environ["lev_coord"] +omega_var = os.environ["omega_var"] +WK_DIR = os.environ["WK_DIR"] +OBS_DATA = os.environ["OBS_DATA"] +CASENAME = os.environ["CASENAME"] + +#from numba import jit +#@jit(nopython=True) +def corr2d(a1,a2): + ij=np.shape(a1[:]) + mid_corr=np.zeros(ij[1:]) + for i in np.arange(ij[1]): + for j in np.arange(ij[2]): + mid_corr[i,j]=np.corrcoef(a1[:,i,j],a2[:,i,j])[0,1] + return mid_corr + +#====================== deriving model output ======================= +def top_heaviness_ratio_robustness_calc_model(reanalysis_path, reanalysis_var): + # read omega and land-sea fraction data + ds= xr.open_dataset(reanalysis_path) + lev_model=ds[lev_coord].values + lat_model=ds[lat_coord].values + lon_model=ds[lon_coord].values + isort=np.argsort(lev_model)[::-1] # descending + mid_omega_model=ds[reanalysis_var].values # mon x lev x lat x lon; for the sample data (JAS over 2000-2019) + mid_omega_model=mid_omega_model[:,isort] + mid_omega_model_ltm=np.nanmean(mid_omega_model,axis=0) # lev x lat x lon + #======================Interpolation======================= + dp=lev_model[-1]-lev_model[0] + levs_interp=np.linspace(lev_model[0], lev_model[-1], num=len(lev_model)) + f=interpolate.interp1d(lev_model, mid_omega_model_ltm, kind='cubic', axis=0) # you can choose linear which consumes less time + mid_omega_model_ltm=f(levs_interp) + #======================deriving O1 and O2======================= + # construct Q1_model as half a sine wave and Q2_model as a full sine wave + # two base functions; Q1: idealized deep convection profile; Q2: Deep stratiform profile (Back et al. 2017) + Q1_model=np.zeros(len(levs_interp)) + Q2_model=np.zeros(len(levs_interp)) + for i in range(len(levs_interp)): + Q1_model[i]=-np.sin(np.pi*(levs_interp[i]-levs_interp[0])/dp) + Q2_model[i]=np.sin(2*np.pi*(levs_interp[i]-levs_interp[0])/dp) + #Normalize + factor=integrate.trapz(Q1_model*Q1_model,lev_model)/dp + Q1_model=Q1_model/np.sqrt(factor) + factor=integrate.trapz(Q2_model*Q2_model,lev_model)/dp + Q2_model=Q2_model/np.sqrt(factor) + #======================calculate explained variance over the globe======================= + # Often times, the pres level is not equally distributed. People tend to increase density in the + # ... upper and lower atmopshere, leaving a less dense distribution in the mid atmosphere. + # ... Thus, it is important to weight Q1 and Q2 before we calculate R2 + # Remember, we already take weights into account when calculating O1 and O2 + mid_Q1_model=Q1_model/np.sqrt(np.sum(Q1_model**2)) # unitize Q1 + mid_Q2_model=Q2_model/np.sqrt(np.sum(Q2_model**2)) # unitize Q2 + # calc ltm O1 and O2, because calculating correlation is not a linear operation + OQ1_model_ltm=np.nansum(mid_Q1_model[:,None,None]*mid_omega_model_ltm,axis=0) + OQ2_model_ltm=np.nansum(mid_Q2_model[:,None,None]*mid_omega_model_ltm,axis=0) + OQQ1_model=OQ1_model_ltm[None,:,:]*mid_Q1_model[:,None,None] # reconstruct Q1 field + OQQ2_model=OQ2_model_ltm[None,:,:]*mid_Q2_model[:,None,None] # reconstruct Q2 field + OQQ_sum=OQQ1_model+OQQ2_model + corr2d_model=corr2d(mid_omega_model_ltm,OQQ_sum) + R2_model=corr2d_model**2 + #====================== setting up figures ======================= + ilat=np.argsort(lat_model) + ilon=np.argsort(lon_model) + #====================== R2 ======================= + # R2 measures the proportion of ltm omega profile explained by Q1 and Q2 + fig, axes = plt.subplots(figsize=(8,4)) + mmid=R2_model + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(0,1.,0.1) + im0 = axes.contourf(lon_model, lat_model, mmid, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_model[ilon][0],lon_model[ilon][-1],60) + lat_grid = np.arange(lat_model[ilat][0],lat_model[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('$R^{2}$ Between Recon. Omega & Original',loc='center',fontsize=18) + fig.savefig(WK_DIR+"/model/"+CASENAME+"_R2.png", format='png',bbox_inches='tight') + + print("Plotting Completed") + + +top_heaviness_ratio_robustness_calc_model(os.environ["OMEGA_FILE"],os.environ["omega_var"]) + +#====================== deriving obs output ======================= +def top_heaviness_ratio_robustness_calc_obs(obs_data_full_dir): + # read omega + ds= xr.open_dataset(obs_data_full_dir) + lev_obs=ds['lev'].values + lat_obs=ds['lat'].values + lon_obs=ds['lon'].values + isort=np.argsort(lev_obs)[::-1] # descending + mid_omega_obs=ds['omega'].values # mon x lev x lat x lon; for the sample data (July over 2000-2019) + mid_omega_obs=mid_omega_obs[:,isort] + mid_omega_obs_ltm=np.nanmean(mid_omega_obs,axis=0) # lev x lat x lon + #======================Interpolation======================= + dp=lev_obs[-1]-lev_obs[0] + levs_interp=np.linspace(lev_obs[0], lev_obs[-1], num=len(lev_obs)) + f=interpolate.interp1d(lev_obs, mid_omega_obs_ltm, kind='cubic', axis=0) # you can choose linear which consumes less time + mid_omega_obs_ltm=f(levs_interp) + #======================deriving O1 and O2======================= + # construct Q1_obs as half a sine wave and Q2_obs as a full sine wave + # two base functions; Q1: idealized deep convection profile; Q2: Deep stratiform profile (Back et al. 2017) + Q1_obs=np.zeros(len(levs_interp)) + Q2_obs=np.zeros(len(levs_interp)) + for i in range(len(levs_interp)): + Q1_obs[i]=-np.sin(np.pi*(levs_interp[i]-levs_interp[0])/dp) + Q2_obs[i]=np.sin(2*np.pi*(levs_interp[i]-levs_interp[0])/dp) + #Normalize + factor=scipy.integrate.trapz(Q1_obs*Q1_obs,lev_obs)/dp + Q1_obs=Q1_obs/np.sqrt(factor) + factor=scipy.integrate.trapz(Q2_obs*Q2_obs,lev_obs)/dp + Q2_obs=Q2_obs/np.sqrt(factor) + #======================calculate explained variance over the globe======================= + # Often times, the pres level is not equally distributed. People tend to increase density in the + # ... upper and lower atmopshere, leaving a less dense distribution in the mid atmosphere. + # ... Thus, it is important to weight Q1 and Q2 before we calculate explained variance + # Remember, we already take weights into account when calculating O1 and O2 + mid_Q1_obs=Q1_obs/np.sqrt(np.sum(Q1_obs**2)) # unitize Q1 + mid_Q2_obs=Q2_obs/np.sqrt(np.sum(Q2_obs**2)) # unitize Q2 + # calc ltm O1 and O2, because calculating correlation is not a linear operation + OQ1_obs_ltm=np.nansum(mid_Q1_obs[:,None,None]*mid_omega_obs_ltm,axis=0) + OQ2_obs_ltm=np.nansum(mid_Q2_obs[:,None,None]*mid_omega_obs_ltm,axis=0) + OQQ1_obs=OQ1_obs_ltm[None,:,:]*mid_Q1_obs[:,None,None] # reconstruct Q1 field + OQQ2_obs=OQ2_obs_ltm[None,:,:]*mid_Q2_obs[:,None,None] # reconstruct Q2 field + OQQ_sum=OQQ1_obs+OQQ2_obs + corr2d_obs=corr2d(mid_omega_obs_ltm,OQQ_sum) + R2_obs=corr2d_obs**2 + #====================== setting up figures ======================= + ilat=np.argsort(lat_obs) + ilon=np.argsort(lon_obs) + #====================== R2 ======================= + # R2 measures the proportion of ltm omega profile explained by Q1 and Q2 + fig, axes = plt.subplots(figsize=(8,4)) + mmid=R2_obs + axes = plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) + clevs=np.arange(0,1.,0.1) + im0 = axes.contourf(lon_obs, lat_obs, mmid, clevs,cmap = plt.get_cmap('RdBu_r'),extend='both', + transform=ccrs.PlateCarree()) + axes.coastlines() + lon_grid = np.arange(lon_obs[ilon][0],lon_obs[ilon][-1],60) + lat_grid = np.arange(lat_obs[ilat][0],lat_obs[ilat][-1],30) + # set x labels + axes.set_xticks(lon_grid, crs=ccrs.PlateCarree()) + axes.set_xticklabels(lon_grid, rotation=0, fontsize=14) + lon_formatter = cticker.LongitudeFormatter() + axes.xaxis.set_major_formatter(lon_formatter) + # set y labels + axes.set_yticks(lat_grid, crs=ccrs.PlateCarree()) + axes.set_yticklabels(lat_grid, rotation=0, fontsize=14) + lat_formatter = cticker.LatitudeFormatter() + axes.yaxis.set_major_formatter(lat_formatter) + # colorbar + fig.colorbar(im0, ax=axes, orientation="horizontal", pad=0.15,shrink=.9,aspect=45) + axes.set_title('$R^{2}$ Between Recon. Omega & Original',loc='center',fontsize=18) + fig.tight_layout() + fig.savefig(WK_DIR+"/obs/ERA5_R2_Between_Recon_Omega&Original_2000_2019_July.png", format='png',bbox_inches='tight') + + print("Plotting Completed") + + +top_heaviness_ratio_robustness_calc_obs(OBS_DATA+'/ERA5_omega_mon_2000_2019_July.nc') + + + + + + + + + +