WRF mode temperature threshold issue

[CUIyige]

1.Recently, I attempted to use the Polar WRF4.3.3 model to simulate the explosive warming event that occurred in Antarctica on March 18, 2022. The model was driven by ERA5 data, and during the simulation, I found something incomprehensible. I set the start times to March 12, March 17, and March 18 respectively, and the model seemed unable to capture this heatwave. The simulation with the start time of March 12 was initially very close to the observed values, but after March 17, there was a significant difference between the simulated values and the observed values. The simulations with the start time of March 17 and 18 only showed a close relationship with the observed values in the first three hours, while the simulation results showed a similar relationship after the fourth hour. The same result as on the 12th was obtained - unable to continue simulating normal temperatures, After discussing with the teacher, I am considering whether there is a mechanism in the WRF mode that prevents it from simulating excessively high temperatures. Once the mode detects a temperature threshold, it will automatically drop and not allow the mode to simulate temperatures exceeding this threshold. The attachments are my namelist-wps and namelist-INput,
2. There is another question. I tried to use Matlab to draw the 500hPa U-V wind generated by my pattern, but found that my wind direction and speed were not very suitable. Below is the code I have drawn. What's wrong with my code, or what's wrong with the wind simulation of the pattern?

 

[CUIyige]

3.namelist.wps:
&share
wrf_core = 'ARW',
max_dom = 1,
start_date = '2022-03-12_00:00:00','2022-03-12_00:00:00','2022-03-12_00:00:00',
end_date = '2022-03-26_00:00:00','2022-03-26_00:00:00','2022-03-26_00:00:00',
interval_seconds = 21600,
io_form_geogrid = 2,
debug_level = 1000
/

&geogrid
parent_id = 1, 1, 1,
parent_grid_ratio = 1, 3, 3,
i_parent_start = 1, 137, 128,
j_parent_start = 1, 67, 155,
s_we = 1, 1, 1,
e_we = 760,, 382,, 562,
s_sn = 1, 1, 1,
e_sn = 600,, 364,, 478,
geog_data_res = 'default','default',
dx = 9000,
dy = 9000,
map_proj = 'polar',
ref_lat = -90.0,
ref_lon = 0.0,
truelat1 = -90.0,
truelat2 = -90.0,
stand_lon = 0.0,
geog_data_path = '/g13/yinjf/geog39/',
/

&ungrib
out_format = 'WPS',
prefix = 'ERA_PRES',
/

&metgrid
fg_name = 'ERA_SFC', 'ERA_PRES',
io_form_metgrid = 2,
/

&mod_levs
press_pa = 201300 , 200100 , 100000 , 95000 , 90000 , 85000 , 80000 , 75000 , 70000 , 65000 , 60000 , 55000 , 50000 , 45000 , 40000 , 35000 , 30000 , 25000 , 20000 , 15000 , 10000 , 5000 , 1000,
/

 

[CUIyige]

4.namelist.input:
&time_control
run_days = 0,
run_hours = 0,
run_minutes = 0,
run_seconds = 0,
start_year = 2022, 2022, 2022, 2022,
start_month = 03, 03, 03, 03,
start_day = 18, 18, 18, 18,
start_hour = 00, 00, 00, 00,
start_minute = 00, 00, 00, 00,
start_second = 00, 00, 00, 00,
end_year = 2022, 2022, 2022, 2022,
end_month = 03, 03, 03, 08,
end_day = 19, 19, 19, 02,
end_hour = 00, 00, 00, 00,
end_minute = 00, 00, 00, 00,
end_second = 00, 00, 00, 00,
interval_seconds = 21600
input_from_file = .true., .true., .true., .true.,
history_interval = 60, 60, 60, 10,
frames_per_outfile = 1000, 1000, 1, 1000,
restart = .false.,
restart_interval = 50000,
io_form_history = 2
io_form_restart = 2
io_form_input = 2
fine_input_stream = 0, 0, 0, 0,
io_form_auxinput2 = 2,
io_form_boundary = 2
all_ic_times = .false.
auxinput11_interval = 1,
auxinput11_end_h = 12,
debug_level = 0
iofields_filename = "my_file_d01.txt", "my_file_d02.txt", "my_file_d03.txt",
auxinput4_inname=wrflowinp_d<domain>,
auxinput4_interval=360,360,
io_form_auxinput4=2,
/

&domains
time_step = 30,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_vert = 58, 58, 58, 58,
p_top_requested = 5000,
interp_type = 1
lagrange_order = 1
force_sfc_in_vinterp = 6
zap_close_levels = 500
adjust_heights = .false.
num_metgrid_levels = 38,
num_metgrid_soil_levels = 4,
eta_levels = 1.000, 0.998, 0.997, 0.992, 0.985, 0.978, 0.969, 0.960, 0.950,
0.938, 0.925, 0.910, 0.894, 0.876, 0.857, 0.835, 0.812,
0.787, 0.760, 0.731, 0.700, 0.668, 0.635, 0.600, 0.565,
0.530, 0.494, 0.458, 0.423, 0.388, 0.355, 0.323, 0.293,
0.264, 0.237, 0.212, 0.188, 0.167, 0.147, 0.130, 0.114,
0.099, 0.086, 0.074, 0.064, 0.054, 0.046, 0.039, 0.032,
0.027, 0.022, 0.017, 0.013, 0.010, 0.007, 0.004, 0.002,
0.000
dx = 9000, 3000, 1000,
dy = 9000, 3000, 1000,
grid_id = 1, 2, 3, 4,
parent_id = 0, 1, 2, 3,
parent_grid_ratio = 1, 3, 3, 3,
i_parent_start = 1, 310, 290, 220,
j_parent_start = 1, 75, 530, 85,
e_we = 760, 973, 1231,
e_sn = 600, 1231, 1591,
parent_time_step_ratio = 1, 3, 3, 3
feedback = 0,
smooth_option = 0
interp_type = 2,
smooth_cg_topo = .false.,
force_sfc_in_vinterp = 1,
lagrange_order = 1,
lowest_lev_from_sfc = .true.,
rh2qv_wrt_liquid = .true.,
rh2qv_method = 2,
/

&physics
mp_physics =53, 53, 53, 53,
ra_lw_physics = 4, 4, 4, 4,
ra_sw_physics = 4, 4, 4, 4,
radt = 5, 5, 5, 5,
sf_sfclay_physics = 1, 1, 1, 1,
sf_surface_physics = 2, 2, 2, 2,
bl_pbl_physics = 1, 1, 1, 1,
bldt = 0, 0, 0, 0,
cu_physics = 1, 1, 0, 0,
kfeta_trigger = 1,
cudt = 0, 0, 0, 0,
isfflx = 1,
ifsnow = 1,
icloud = 1,
surface_input_source = 1,
num_soil_layers = 4,
sf_urban_physics = 0, 0, 0, 0,
mp_zero_out = 2,
mp_zero_out_thresh = 1.0e-8,
iz0tlnd = 1,
surface_input_source = 1,
NUM_LAND_CAT = 21,
fractional_seaice = 1,
tice2tsk_if2cold = .true.，
sst_update =1,
/

&fdda
grid_fdda = 1, 1, 1,
gfdda_inname = "wrffdda_d<domain>",
gfdda_end_h = 336, 336, 336,
gfdda_interval_m = 360, 360, 360,
fgdt = 0, 0, 0,
fgdtzero = 0, 0, 0,
if_no_pbl_nudging_uv = 0, 0, 0,
if_no_pbl_nudging_t = 0, 0, 0,
if_no_pbl_nudging_q = 0, 0, 0,
if_no_pbl_nudging_ph = 0, 0, 0,
if_zfac_uv = 0, 0, 0,
k_zfac_uv = 10, 10, 10,
if_zfac_t = 0, 0, 0,
k_zfac_t = 10, 10, 10,
if_zfac_q = 0, 0, 0,
k_zfac_q = 10, 10, 10
if_zfac_ph = 0, 0, 0,
k_zfac_ph = 10, 10, 10,
dk_zfac_uv = 1, 1, 1,
dk_zfac_t = 1, 1, 1,
dk_zfac_ph = 1, 1, 1,
guv = 5.0E-5, 5.0E-5, 5.0E-5,
gt = 5.0E-5, 5.0E-5, 5.0E-5,
gq = 5.0E-6, 5.0E-6, 5.0E-6,
gph = 0.0, 0.0, 0.0,
xwavenum = 3, 3, 3,
ywavenum = 3, 3, 3,
if_ramping = 1,
dtramp_min = 60.0,
io_form_gfdda = 2,
grid_sfdda = 0, 0, 0,
sgfdda_inname = "wrfsfdda_d<domain>",
sgfdda_end_h = 96, 96, 96,
sgfdda_interval_m = 60, 60, 60,
io_form_sgfdda = 2,
pxlsm_soil_nudge = 1, 1, 1,
guv_sfc = 8.3E-4, 8.3E-4, 8.3E-4,
gt_sfc = 8.3E-4, 8.3E-4, 8.3E-4,
gq_sfc = 8.3E-4, 8.3E-4, 8.3E-4,
rinblw = 60, 40, 40,
obs_nudge_opt = 0, 0, 0,
max_obs = 1500000,
fdda_start = 0, 0, 0,
fdda_end = 5760, 5760, 5760,
obs_nudge_wind = 1, 1, 1,
obs_coef_wind = 1.6e-3, 1.6e-3, 1.6e-3,
obs_nudge_temp = 1, 1, 1,
obs_coef_temp = 1.6e-3, 1.6e-3, 1.6e-3,
obs_nudge_mois = 1, 1, 1,
obs_coef_mois = 1.6e-3, 1.6e-3, 1.6e-3,
obs_nudge_pstr = 0, 0, 0,
obs_coef_pstr = 0.0, 0.0, 0.0,
obs_rinxy = 24, 16, 8,
obs_rinsig = 0.1,
obs_twindo = 0.6666667, 0.6666667, 0.6666667,
obs_npfi = 10,
obs_ionf = 1, 1, 1,
obs_idynin = 0,
obs_dtramp = 60.0,
obs_scl_neg_qv_innov = 1,
obs_ipf_nudob = .true.
/

&dynamics
hybrid_opt = 2
rk_ord = 3
w_damping = 1,
diff_opt = 1,
km_opt = 4,
diff_6th_opt = 2, 2, 2, 2,
diff_6th_factor = 0.10, 0.10, 0.10, 0.10,
base_temp = 290.
damp_opt = 3,
zdamp = 5000., 5000., 5000., 5000.,
dampcoef = 0.05, 0.05, 0.05, 0.05,
khdif = 0, 0, 0, 0,
kvdif = 0, 0, 0, 0,
smdiv = 0.1 0.1, 0.1, 0.1,
emdiv = 0.01, 0.01, 0.01, 0.01,
epssm = 0.5 0.5, 0.5, 0.5,
h_mom_adv_order = 5, 5, 5, 5,
v_mom_adv_order = 3, 3, 3, 3,
h_sca_adv_order = 5, 5, 5, 5,
v_sca_adv_order = 3, 3, 3, 3,
non_hydrostatic = .true., .true., .true., .true.,
moist_adv_opt = 1, 1, 1, 1,
scalar_adv_opt = 1, 1, 1, 1,
tke_adv_opt = 1, 1, 1, 1,
gwd_opt = 0,
use_theta_m = 0,
use_q_diabatic = 0,
/

&bdy_control
spec_bdy_width = 5,
spec_zone = 1,
relax_zone = 4,
specified = .true., .false.,.false., .false.,
nested = .false., .true., .true., .true.,
/

&grib2
/

&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/

 

[CUIyige]

6.MATLAB:
% ncdisp('wrfout_d01_2022-03-26_00_00_00');

addpath \wrf_post



f='wrfout_d01_2022-03-26_00_00_00';

lat=ncread(f,'XLAT');

lon=ncread(f,'XLONG');

latu=ncread(f,'XLAT_U');

lonu=ncread(f,'XLONG_U');



slp = GetVar(f,'height',double(lon),double(lat),500,'all',0,'p');

U = GetVar(f,'U',double(lonu),double(latu),500,'all',0,'p');

V = GetVar(f,'V',double(lonu),double(latu),500,'all',0,'p');





flag=15;



figure(1);

set(gcf,'position',[10,-20,700,700]);

m_proj('stereographic','lat',-90,'long',0,'radius',35);

[cs,h]=m_contourf(lon,lat,slp/10,100,'linestyle','none');hold on;

m_contour(lon,lat,slp/10,[480:4:540],'color',[.9 .9 .9],'linewidth',2,'ShowText','on');

m_quiver(lonu(1:flag:759,1:10:599),latu(1:flag:759,1:10:599),V(1:flag:759,1:10:599),U(1:flag:759,1:10:599),3,'k-');

m_coast('line','color','k','linewidth',2);

m_grid('xtick',6,'xticklabels',{},'tickdir','in','ytick',-[80 70 60 50],'linest',':');

h = colorbar;

set(get(h,'title'),'string','dagpm')

caxis([480,520]);

 

[Ming Chen]

1.Recently, I attempted to use the Polar WRF4.3.3 model to simulate the explosive warming event that occurred in Antarctica on March 18, 2022. The model was driven by ERA5 data, and during the simulation, I found something incomprehensible. I set the start times to March 12, March 17, and March 18 respectively, and the model seemed unable to capture this heatwave. The simulation with the start time of March 12 was initially very close to the observed values, but after March 17, there was a significant difference between the simulated values and the observed values. The simulations with the start time of March 17 and 18 only showed a close relationship with the observed values in the first three hours, while the simulation results showed a similar relationship after the fourth hour. The same result as on the 12th was obtained - unable to continue simulating normal temperatures, After discussing with the teacher, I am considering whether there is a mechanism in the WRF mode that prevents it from simulating excessively high temperatures. Once the mode detects a temperature threshold, it will automatically drop and not allow the mode to simulate temperatures exceeding this threshold. The attachments are my namelist-wps and namelist-INput,

Your namlelis files look fine. However, I would suggest that you turn off grid nudging, i.e.,

grid_fdda = 0,

This is because grid nudfing will force WRF output to be closer to the large scale forcing data, which often miss detailed meso- and micro-scale features due to its coarse resolution. As a result, WRF simulations could be 'distorted' by nugding.
Also, I believe that Polar WRF has been modified to make it more suitable for studies in polar region. Thus it is different to standard WRF released by NCAR. Please talk to people in Polar WRF group and see whether they have some insights about the issues you have.

2. There is another question. I tried to use Matlab to draw the 500hPa U-V wind generated by my pattern, but found that my wind direction and speed were not very suitable. Below is the code I have drawn. What's wrong with my code, or what's wrong with the wind simulation of the pattern?

I am not familiar with Matlab. Are you able to process wrfout using NCL or Python?

 

[CUIyige]

您的 namlelis 文件看起来不错。但是，我建议您关闭网格微移，即
[代码]grid_fdda = 0，[/代码]
这是因为网格微移将迫使 WRF 输出更接近大尺度强制数据，而大尺度强制数据由于其粗略的分辨率而经常错过详细的中尺度和微观尺度特征。因此，WRF模拟可能会被微波“扭曲”。
此外，我相信 Polar WRF 已经过修改，使其更适合极地地区的研究。因此，它与NCAR发布的标准WRF不同。请与 Polar WRF 小组的人交谈，看看他们是否对您遇到的问题有一些见解。


我不熟悉Matlab。你能使用NCL或Python处理wrfout吗？

Thank you for your reply.
1. I tried to use Python to draw, but it seems that the result drawn by my code has obvious projection problems. There are strange polygons in the middle, and it can be seen that the wind direction is also incorrect and there is a mismatch between the circulation. Do you think this is a problem with my code or the pattern simulation? Attached is my drawing code
2.The second figure is a line graph of my simulated 2-meter temperature values. It can be observed that from the 17th and 18th simulations, the simulated values decreased without warning after 3 hours and converged with the simulation results on the 12th. This is illogical. What do you think could be the reason for this

 

[CUIyige]

Thank you for your reply.
1. I tried to use Python to draw, but it seems that the result drawn by my code has obvious projection problems. There are strange polygons in the middle, and it can be seen that the wind direction is also incorrect and there is a mismatch between the circulation. Do you think this is a problem with my code or the pattern simulation? Attached is my drawing code
2.The second figure is a line graph of my simulated 2-meter temperature values. It can be observed that from the 17th and 18th simulations, the simulated values decreased without warning after 3 hours and converged with the simulation results on the 12th. This is illogical. What do you think could be the reason for this

View attachment 14861View attachment 14862

import numpy as np
from netCDF4 import Dataset
import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker
import cartopy.crs as ccrs
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cartopy.crs as crs
try:
import pykdtree.kdtree
_IS_PYKDTREE = True
except ImportError:
import scipy.spatial
_IS_PYKDTREE = False
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
from cartopy.feature import NaturalEarthFeature
from wrf import to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords, vertcross, smooth2d, CoordPair, GeoBounds,interpline
import warnings
warnings.filterwarnings('ignore')
wrf_file="E:/Polar WRF/wrfout_d01_2022-03-17_00_00_00"
# Open the NetCDF file
ncfile = Dataset(wrf_file)
times = getvar(ncfile, "Times")
dataset = nc.Dataset(wrf_file)
# 假设我们要选择第一个时间步（索引为0）
time_idx = 0
p = getvar(ncfile, "pressure", timeidx=time_idx)
z = getvar(ncfile, "z", units="dm", timeidx=time_idx)
ua = getvar(ncfile, "ua", units="kt", timeidx=time_idx)
va = getvar(ncfile, "va", units="kt", timeidx=time_idx)
wspd = getvar(ncfile, "wspd_wdir", units="kts", timeidx=time_idx)[0,:]

# Interpolate geopotential height, u, and v winds to 500 hPa
ht_500 = interplevel(z, p, 500)
u_500 = interplevel(ua, p, 500)
v_500 = interplevel(va, p, 500)
wspd_500 = interplevel(wspd, p, 500)

# Get the lat/lon coordinates
lats, lons = latlon_coords(ht_500)
# Get the map projection information
cart_proj = get_cartopy(ht_500)

# 设置地图的边界（如果你需要的话）
fig = plt.figure(figsize=(12,10))
lon_max = 180 #int(np.max(lons))
lon_min = -180 #int(np.min(lons))
lat_max = -60 #int(np.max(lats))
lat_min = -90 #int(np.min(lats))
ax= plt.axes(projection=cart_proj)
projection=cart_proj
ax.set_extent([lon_min, lon_max, lat_min, lat_max])
ax = fig.add_subplot(1, 1, 1, projection=ccrs.SouthPolarStereo())
ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())
# 添加海岸线等特征
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),linewidth=2.0, edgecolor="black",zorder=10)
# 添加网格线
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=0.5, color='grey', linestyle='--')
gl.xlines = list(np.arange(lon_min, lon_max, 1.))
gl.ylines = list(np.arange(lat_min, lat_max, 1.))
ax.set_title('Sea level pressure', fontsize=9)


# 调整布局并显示图形
plt.tight_layout()

# Add the 500 hPa geopotential height contours
levels = np.arange(520., 580., 6.)
contours = plt.contour(to_np(lons), to_np(lats), to_np(ht_500),
levels=levels, colors="black",
transform=crs.PlateCarree())
plt.clabel(contours, inline=1, fontsize=10, fmt="%i")
# Add the wind speed contours
levels = [25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 110, 120]
wspd_contours = plt.contourf(to_np(lons), to_np(lats), to_np(wspd_500),
levels=levels,
cmap=get_cmap("rainbow"),
transform=crs.PlateCarree())
plt.colorbar(wspd_contours, ax=ax, orientation="horizontal", pad=.05)

# Add the 500 hPa wind barbs, only plotting every 125th data point.
h1 = ax.quiver(to_np(lons[::30,::30]), to_np(lats[::30,::30]),
to_np(u_500[::30, ::30]), to_np(v_500[::30, ::30]),
color='black',
scale=1300,
transform=crs.PlateCarree())


# Set the map bounds
ax.set_xlim(cartopy_xlim(ht_500))
ax.set_ylim(cartopy_ylim(ht_500))
ax.gridlines()
plt.title("500 MB Height (dm), Wind Speed (kt), Barbs (kt)")

# 添加一个点（例如，Dome C）
lon, lat = 123.35, -75.11 # 东经123.35，南纬75.11（注意：这里的坐标可能不在你的数据范围内）
ax.plot(lon, lat, 'o', color='green', markersize=13, transform=ccrs.PlateCarree())
ax.text(lon, lat, 'Dome C', ha='center', va='center', fontsize=15,
transform=ccrs.PlateCarree())
# must call draw in order to get the axis boundary used to add ticks:
fig.canvas.draw()
plt.show()