获取Piggy_Packages
Piggy_Packages 是⼀个在 Windows 操作系统环境下“开箱即⽤”的便捷式气象科研环境。⽆需复杂的安装和配置步骤，就能快速投⼊到实际的科研⼯作中。
还没有Piggy_Packages的同学，请参考这篇帖子获取：
Piggy_Packages V2026.1 帮助文档（一）开箱即用
四维变分同化试验
四维变分同化与三维变分同化的主要区别在于引入了时间维度和数值预报模式。它假设在一段同化窗口内，大气状态的演变遵循动力学模式约束。通过这种方式，即使是在不同时间点获得的观测资料，也能在模式物理规律的约束下被一致地同化。
四维变分同化通常需要消耗大量计算资源，因此很难在业务上开展。这一节，我们来体验一下四维变分同化试验的运行流程。
本节同化运行时间和前面讲的三维变分保持一致。
背景数据准备
由于WRFDA 4D-Var 目前不支持混合垂直坐标和湿位温，我们在预报时必须将其关闭。
&dynamics
hybrid_opt = 0,
use_theta_m = 0,...【其他选项】
/
首先我们回到Home目录，创建一个名称为Case7的目录，并进入。将WRF拷贝进来，进入run目录。
mkdir -p Case7 cd Case7 cp /pp_model/WRF-4.7.1/WRF-Real . cd WRF-Real/run将Case1中准备好的气象场资料链接过来，编辑namelist：
ln -s ~/Case1/WPS/met_em.* . open -e namelist.input粘贴以下内容：
&time_control
run_days = 1,
run_hours = 1,
run_minutes = 0,
run_seconds = 0,
start_year = 2026,
start_month = 03,
start_day = 29,
start_hour = 00,
end_year = 2026,
end_month = 03,
end_day = 31,
end_hour = 01,
interval_seconds = 10800
input_from_file = .true.,
history_interval = 60,
frames_per_outfile = 1,
restart = .false.,
restart_interval = 7200,
io_form_history = 2
io_form_restart = 2
io_form_input = 2
io_form_boundary = 2
/&domains
time_step = 120,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 50,
e_sn = 50,
e_vert = 48,
dzbot = 30.
dzstretch_s = 1.11
dzstretch_u = 1.10
p_top_requested = 5000,
num_metgrid_levels = 34,
num_metgrid_soil_levels = 4,
dx = 81000,
dy = 81000,
grid_id = 1,
parent_id = 0,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/&physics
mp_physics = 3, -1,
cu_physics = 1, -1,
ra_lw_physics = 4, -1,
ra_sw_physics = 4, -1,
bl_pbl_physics = 1, -1,
sf_sfclay_physics = 1, -1,
sf_surface_physics = 2, -1,
radt = 15, 15,
bldt = 0, 0,
cudt = 0, 0,
icloud = 1,
num_land_cat = 21,
sf_urban_physics = 0, 0,
fractional_seaice = 1,
/&fdda
/&dynamics
hybrid_opt = 0,
use_theta_m = 0,
w_damping = 0,
diff_opt = 2, 2,
km_opt = 4, 4,
diff_6th_opt = 0, 0,
diff_6th_factor = 0.12, 0.12,
base_temp = 290.
damp_opt = 3,
zdamp = 5000., 5000.,
dampcoef = 0.2, 0.2,
khdif = 0, 0,
kvdif = 0, 0,
non_hydrostatic = .true., .true.,
moist_adv_opt = 1, 1,
scalar_adv_opt = 1, 1,
gwd_opt = 1, 0,
/&bdy_control
spec_bdy_width = 5,
specified = .true.
/&grib2
/&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/
运行real.exe和wrf.exe
mpirun -np 4 ./real.exe mpirun -np 4 ./wrf.exe四维变分同化流程
准备好背景场之后，我们来运行同化。
运行四维变分同化需要消耗大量资源，因此有必要对观测进行稀疏化，稀疏化网格大小近似于模式网格。在namelist.input中设置如下：
&wrfvar4
thin_conv = .true.,
thin_mesh_conv = 80,
thin_conv_opt = 1,
/
业务运行通常设置6小时观测时间窗，本文作为演示目的，只设置1小时时间窗。
即本文以3月29日00时起报的24小时预报场作为背景场，同化30日00时-01时观测资料。
编辑namelist.input
open -e namelist.input粘贴以下内容：
&wrfvar1
var4d = .true.,
var4d_bin = 1200,
var4d_lbc = .true.,
print_detail_grad = .true.,
/&wrfvar2
analysis_accu = 900,
/&wrfvar3
ob_format = 1,
/&wrfvar4
thin_conv = .true.,
thin_mesh_conv = 80,
thin_conv_opt = 1,
/&wrfvar5
check_max_iv = .true.,
/&wrfvar6
max_ext_its = 1,
ntmax = 50,
orthonorm_gradient = .true.,
/&wrfvar7
cv_options = 3,
/&wrfvar11
/&wrfvar18
analysis_date = "2026-03-30_00:00:00.0000",
/&wrfvar21
time_window_min = "2026-03-30_00:00:00.0000",
/&wrfvar22
time_window_max = "2026-03-30_01:00:00.0000",
/
&time_control
run_hours = 1,
start_year = 2026,
start_month = 03,
start_day = 30,
start_hour = 00,
end_year = 2026,
end_month = 03,
end_day = 30,
end_hour = 01,
interval_seconds = 3600,
/&domains
time_step = 120,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 50,
e_sn = 50,
e_vert = 48,
dzbot = 30.
dzstretch_s = 1.11
dzstretch_u = 1.10
p_top_requested = 5000,
num_metgrid_levels = 34,
num_metgrid_soil_levels = 4,
dx = 81000,
dy = 81000,
grid_id = 1,
parent_id = 0,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/&physics
mp_physics = 3,
mp_physics_ad = 99,
cu_physics = 1,
ra_lw_physics = 4,
ra_sw_physics = 4,
bl_pbl_physics = 1,
sf_sfclay_physics = 1,
sf_surface_physics = 2,
num_land_cat = 21,
/&dynamics
hybrid_opt = 0,
use_theta_m = 0,
w_damping = 0,
diff_opt = 2, 2,
km_opt = 4, 4,
diff_6th_opt = 0, 0,
diff_6th_factor = 0.12, 0.12,
base_temp = 290.
damp_opt = 3,
zdamp = 5000., 5000.,
dampcoef = 0.2, 0.2,
khdif = 0, 0,
kvdif = 0, 0,
non_hydrostatic = .true., .true.,
moist_adv_opt = 1, 1,
scalar_adv_opt = 1, 1,
gwd_opt = 1, 0,
base_temp = 290.,
/&bdy_control
spec_bdy_width = 5,
specified = .true.,
/&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/&perturbation
trajectory_io = .true.,
enable_identity = .false.,
/
删除 wrfinput_d01，将 wrfout_d01_2026-03-30_00:00:00 文件链接为fg和wrfinput_d01，我们在运行WRF时预报到了30日01时，与同化时间窗对应，将 wrfout_d01_2026-03-30_01:00:00 文件链接为 fg02。
rm wrfinput_d01 ln -s wrfout_d01_2026-03-30_00:00:00 fg ln -s wrfout_d01_2026-03-30_01:00:00 fg02 ln -s wrfout_d01_2026-03-30_00:00:00 wrfinput_d01WRFDA默认使用双精度（REAL8）编译，需要将辐射数据替换为双精度版本：
rm RRTM_DATA RRTMG_LW_DATA RRTMG_SW_DATA mv RRTM_DATA_DBL RRTM_DATA mv RRTMG_LW_DATA_DBL RRTMG_LW_DATA mv RRTMG_SW_DATA_DBL RRTMG_SW_DATA拷贝背景误差协方差矩阵：
cp /pp_model/WRF-4.7.1/WRFDA/var/run/be.dat.cv3 be.dat拷贝观测资料。由于在三维变分中，我们已经体验过将prepbufr格式观测资料转换为little_r格式了，本节我们直接使用prepbufr格式观测资料。
cp ~/Case5/WRFDA/var/da/prepbufr.gdas.20260330.t00z.nr.48h ob.bufr编辑namelist.input
open -e namelist.input粘贴以下内容：
&wrfvar1
var4d = .true.,
var4d_bin = 1200,
var4d_lbc = .true.,
print_detail_grad = .true.,
/&wrfvar2
analysis_accu = 900,
/&wrfvar3
ob_format = 1,
/&wrfvar4
thin_conv = .true.,
thin_mesh_conv = 80,
thin_conv_opt = 1,
/&wrfvar5
check_max_iv = .true.,
/&wrfvar6
max_ext_its = 1,
ntmax = 50,
orthonorm_gradient = .true.,
/&wrfvar7
cv_options = 3,
/&wrfvar11
/&wrfvar18
analysis_date = "2026-03-30_00:00:00.0000",
/&wrfvar21
time_window_min = "2026-03-30_00:00:00.0000",
/&wrfvar22
time_window_max = "2026-03-30_01:00:00.0000",
/
&time_control
run_hours = 1,
start_year = 2026,
start_month = 03,
start_day = 30,
start_hour = 00,
end_year = 2026,
end_month = 03,
end_day = 30,
end_hour = 01,
interval_seconds = 3600,
/&domains
time_step = 120,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 50,
e_sn = 50,
e_vert = 48,
dzbot = 30.
dzstretch_s = 1.11
dzstretch_u = 1.10
p_top_requested = 5000,
num_metgrid_levels = 34,
num_metgrid_soil_levels = 4,
dx = 81000,
dy = 81000,
grid_id = 1,
parent_id = 0,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/&physics
mp_physics = 3,
mp_physics_ad = 99,
cu_physics = 1,
ra_lw_physics = 4,
ra_sw_physics = 4,
bl_pbl_physics = 1,
sf_sfclay_physics = 1,
sf_surface_physics = 2,
num_land_cat = 21,
/&dynamics
hybrid_opt = 0,
use_theta_m = 0,
w_damping = 0,
diff_opt = 2, 2,
km_opt = 4, 4,
diff_6th_opt = 0, 0,
diff_6th_factor = 0.12, 0.12,
base_temp = 290.
damp_opt = 3,
zdamp = 5000., 5000.,
dampcoef = 0.2, 0.2,
khdif = 0, 0,
kvdif = 0, 0,
non_hydrostatic = .true., .true.,
moist_adv_opt = 1, 1,
scalar_adv_opt = 1, 1,
gwd_opt = 1, 0,
base_temp = 290.,
/&bdy_control
spec_bdy_width = 5,
specified = .true.,
/&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/&perturbation
trajectory_io = .true.,
enable_identity = .false.,
/
运行WRFDA：
mpirun -np 4 /pp_model/WRF-4.7.1/WRFDA/var/build/da_wrfvar.exe运行需要较长时间，请耐心等待。
更新边界条件
编辑配置文件parame.in
touch parame.in open -e parame.in粘贴以下内容：
&control_param
da_file = './wrfvar_output'
wrf_bdy_file = './wrfbdy_d01'
da_file_02 = './ana02'
domain_id = 1
update_lateral_bdy = .true.
update_low_bdy = .false.
update_lsm = .false.
iswater = 16
var4d_lbc = .true.
/
更新边界条件：
/pp_model/WRF-4.7.1/WRFDA/var/build/da_update_bc.exe同化结果分析
计算同化分析增量：
cdo sub wrfvar_output fg analysis_increment.nc可视化查看：
Panoply analysis_increment.nc运行预报
将同化结果作为预报初始场：
ln -sf wrfvar_output wrfinput_d01编辑namelist.input
open -e namelist.input粘贴以下内容：
&time_control
run_days = 1,
run_hours = 0,
run_minutes = 0,
run_seconds = 0,
start_year = 2026,
start_month = 03,
start_day = 30,
start_hour = 00,
end_year = 2026,
end_month = 03,
end_day = 31,
end_hour = 00,
interval_seconds = 10800
input_from_file = .true.,
history_interval = 60,
frames_per_outfile = 1,
restart = .false.,
restart_interval = 7200,
io_form_history = 2
io_form_restart = 2
io_form_input = 2
io_form_boundary = 2
/&domains
time_step = 120,
time_step_fract_num = 0,
time_step_fract_den = 1,
max_dom = 1,
e_we = 50,
e_sn = 50,
e_vert = 48,
dzbot = 30.
dzstretch_s = 1.11
dzstretch_u = 1.10
p_top_requested = 5000,
num_metgrid_levels = 34,
num_metgrid_soil_levels = 4,
dx = 81000,
dy = 81000,
grid_id = 1,
parent_id = 0,
i_parent_start = 1,
j_parent_start = 1,
parent_grid_ratio = 1,
parent_time_step_ratio = 1,
feedback = 1,
smooth_option = 0
/&physics
mp_physics = 3, -1,
cu_physics = 1, -1,
ra_lw_physics = 4, -1,
ra_sw_physics = 4, -1,
bl_pbl_physics = 1, -1,
sf_sfclay_physics = 1, -1,
sf_surface_physics = 2, -1,
radt = 15, 15,
bldt = 0, 0,
cudt = 0, 0,
icloud = 1,
num_land_cat = 21,
sf_urban_physics = 0, 0,
fractional_seaice = 1,
/&fdda
/&dynamics
hybrid_opt = 0,
use_theta_m = 0,
w_damping = 0,
diff_opt = 2, 2,
km_opt = 4, 4,
diff_6th_opt = 0, 0,
diff_6th_factor = 0.12, 0.12,
base_temp = 290.
damp_opt = 3,
zdamp = 5000., 5000.,
dampcoef = 0.2, 0.2,
khdif = 0, 0,
kvdif = 0, 0,
non_hydrostatic = .true., .true.,
moist_adv_opt = 1, 1,
scalar_adv_opt = 1, 1,
gwd_opt = 1, 0,
/&bdy_control
spec_bdy_width = 5,
specified = .true.
/&grib2
/&namelist_quilt
nio_tasks_per_group = 0,
nio_groups = 1,
/
由于WRF默认是单精度（REAL4）编译的，需要将辐射文件替换回来：
rm RRTM_DATA RRTMG_SW_DATA RRTMG_LW_DATA cp /pp_model/WRF-4.7.1/WRF-Real/run/RRTM* .运行WRF：
mpirun -np 4 ./wrf.exe可视化查看预报结果：
Panoply wrfout_d01_2026-03-31_00:00:00
