The impact of large-scale forcing scales on clouds.

# Import packages

from datetime import datetime
import numpy as np
import xarray as xr
import xwrf

import matplotlib.pyplot as plt

# The Large Scale Forcing Scale in Simulation 6 is 75 km
# The Large Scale Forcing Scale in Simulation 7 is 150 km
# The Large Scale Forcing Scale in Simulation 8 is 300 km

# Load Data
# To recap, wrfstat is 10 minutes-average values, wrfout is a snapshot at a given time.

path_shcu_root = "/data/project/ARM_Summer_School_2024_Data/lasso_tutorial/ShCu/untar"  # on Jupyter

case_date = datetime(2019, 4, 4)

sim_id2 = 6
sim_id3 = 7
sim_id4 = 8

ds_stat6 =  xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id2:04d}/raw_model/wrfstat_d01_{case_date:%Y-%m-%d_12:00:00}.nc")
ds_stat7 =  xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id3:04d}/raw_model/wrfstat_d01_{case_date:%Y-%m-%d_12:00:00}.nc")
ds_stat8 =  xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id4:04d}/raw_model/wrfstat_d01_{case_date:%Y-%m-%d_12:00:00}.nc")

ds_stat6 = ds_stat6.xwrf.postprocess()
ds_stat7 = ds_stat7.xwrf.postprocess()
ds_stat8 = ds_stat8.xwrf.postprocess()

ds_stat7

sim_id2 = 6
ds_cogs6 = xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id2:04d}/obs_model/sgplassocogsdiagobsmod{sim_id2}C1.m1.{case_date:%Y%m%d}.120000.nc")

sim_id3 = 7
ds_cogs7 = xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id3:04d}/obs_model/sgplassocogsdiagobsmod{sim_id3}C1.m1.{case_date:%Y%m%d}.120000.nc")

sim_id4 = 8
ds_cogs8 = xr.open_dataset(f"{path_shcu_root}/{case_date:%Y%m%d}/sim{sim_id4:04d}/obs_model/sgplassocogsdiagobsmod{sim_id4}C1.m1.{case_date:%Y%m%d}.120000.nc")

# Plot cloud core fractions

fig, ax = plt.subplots(ncols=1)

# source_type = 1 is simulations
# source_type = 0 is observations

ds_cogs6["low_cloud_fraction_cogs"].isel(source_type=1).plot(ax=ax, marker="1", label="Sim6-75")
ds_cogs7["low_cloud_fraction_cogs"].isel(source_type=1).plot(ax=ax, marker="1", label="Sim7-150")
ds_cogs8["low_cloud_fraction_cogs"].isel(source_type=1).plot(ax=ax, marker="1", label="Sim8-300")
# Observations
ds_cogs7["low_cloud_fraction_cogs"].isel(source_type=0).plot(ax=ax, marker="1", label="COG-OBS")

ax.legend()
ax.set_title("Comparisons of Simulations for WRF Cloud Fraction")
ax.set_xlabel("Time (UTC)")
ax.set_ylabel("Cloud Fraction")

plt.show()

# Fraction of cloud core grid points

fig, ax = plt.subplots(1, 3, figsize=(20, 5))

ds_stat6["CSP_A_CC"].plot(ax=ax[0], y ="z") #one plot
ds_stat7["CSP_A_CC"].plot(ax=ax[1], y = "z")
ds_stat8["CSP_A_CC"].plot(ax=ax[2], y = "z")

plt.show()

# simulation 6, 7, and 8 from the left to the right.

# Compare TKE with three simulations
fig, ax = plt.subplots(ncols=1)

ds_stat6["CST_TKE"].plot(ax=ax, marker="1", label="Sim6-75")
ds_stat7["CST_TKE"].plot(ax=ax, marker="1", label="Sim7-150")
ds_stat8["CST_TKE"].plot(ax=ax, marker="1", label="Sim8-300")

ax.legend()
ax.set_title("Comparisons of TKE")

plt.show()

# To further investigate the impact of large-scale forcing scales,
# I looked at each component of TKE and specific humidity (water vapor mixing ratio) variance


fig, ax = plt.subplots(4, 3, figsize=(20, 10))

ds_stat6["CSP_U2"].isel(z= slice(0, 150)).plot(ax=ax[0,0], y ="z") #one plot
ds_stat7["CSP_U2"].isel(z= slice(0, 150)).plot(ax=ax[0,1], y = "z")
ds_stat8["CSP_U2"].isel(z= slice(0, 150)).plot(ax=ax[0,2], y = "z")

ds_stat6["CSP_V2"].isel(z= slice(0, 150)).plot(ax=ax[1,0], y ="z") #one plot
ds_stat7["CSP_V2"].isel(z= slice(0, 150)).plot(ax=ax[1,1], y = "z")
ds_stat8["CSP_V2"].isel(z= slice(0, 150)).plot(ax=ax[1,2], y = "z")

ds_stat6["CSP_W2"].isel(z_stag= slice(0, 150)).plot(ax=ax[2,0], y ="z_stag") #one plot
ds_stat7["CSP_W2"].isel(z_stag= slice(0, 150)).plot(ax=ax[2,1], y = "z_stag")
ds_stat8["CSP_W2"].isel(z_stag= slice(0, 150)).plot(ax=ax[2,2], y = "z_stag")

ds_stat6["CSP_QV2"].isel(z= slice(0, 150)).plot(ax=ax[3,0], y ="z") #one plot
ds_stat7["CSP_QV2"].isel(z= slice(0, 150)).plot(ax=ax[3,1], y = "z")
ds_stat8["CSP_QV2"].isel(z= slice(0, 150)).plot(ax=ax[3,2], y = "z")


plt.show()

# U wind variance, V wind variance, W wind variance, the total water vapor mixing ratio from top to bottom.
# Simulation 6, 7, and 8 from the left to the right.

# If you find something from the below figure, please let me know!
# lim
# What I can see is the variances rapidly increase in Simulation 6. 
# I guess that is the atmospheric boundary layer (ABL) growth. 
# That is, the ABL grows at the fastest rate in the smallest Large Scale Forcing Scale. 

# However, the water vapor mixing ratio (the bottom row) is the lowest in Simulation 6. 
# Is it because the strong turbulence in Simulation 6 makes a well-mixed layer??
# If you can interpret this, please let me know.

fig = plt.figure(figsize=(20, 10))

# Create a plot for temperature
ax = fig.add_subplot(4, 3, 1)

ds_stat6["CSP_U2"].isel(z= slice(0, 150)).plot( y ="z") #one plot

#ax.set_title('Temperature')
#ax.set_xlim(110, 130)
#ax.set_ylim(290, 315)

# Create a plot for dewpoint
ax2 = fig.add_subplot(4, 3, 2)
ds_stat7["CSP_U2"].isel(z= slice(0, 150)).plot( y = "z")

ax3 = fig.add_subplot(4, 3, 3)
ds_stat8["CSP_U2"].isel(z= slice(0, 150)).plot( y = "z")

ax4 = fig.add_subplot(4, 3, 4)
ds_stat6["CSP_V2"].isel(z= slice(0, 150)).plot( y ="z") #one plot

ax5 = fig.add_subplot(4, 3, 5)
ds_stat7["CSP_V2"].isel(z= slice(0, 150)).plot( y = "z")

ax6 = fig.add_subplot(4, 3, 6)
ds_stat8["CSP_V2"].isel(z= slice(0, 150)).plot( y = "z")

ax7 = fig.add_subplot(4, 3, 7)
ds_stat6["CSP_W2"].isel(z_stag= slice(0, 150)).plot( y ="z_stag") #one plot

ax8 = fig.add_subplot(4, 3, 8)
ds_stat7["CSP_W2"].isel(z_stag= slice(0, 150)).plot( y = "z_stag")

ax9 = fig.add_subplot(4, 3, 9)
ds_stat8["CSP_W2"].isel(z_stag= slice(0, 150)).plot( y = "z_stag")

ax10 = fig.add_subplot(4, 3, 10)
ds_stat6["CSP_QV2"].isel(z= slice(0, 150)).plot( y ="z") #one plot
#ax10.set_ylim(0, 4)

ax11 = fig.add_subplot(4, 3, 11)
ds_stat7["CSP_QV2"].isel(z= slice(0, 150)).plot( y = "z")
#ax11.set_ylim(0, 4)

ax12 = fig.add_subplot(4, 3, 12)
ds_stat8["CSP_QV2"].isel(z= slice(0, 150)).plot( y = "z")
#ax12.set_ylim(0, 4)