使用python处理wrfout文件

发表于 更新于 分类于

在前文已经成功驱动WRF模式的背景下，记录如何使用python处理wrfout文件。

参考网址：

wrf-python 1.3.4.1 documentation

首先在anaconda环境中安装wrf-python包：

1
conda install -c conda-forge wrf-python

读取文件，获取变量

这是wrf-python中的一些命令与ncl中的一些命令的对应关系。

NCL wrf-python 作用
wrf_user_get_var wrf.getvar 提取和计算变量，变量名参考wrf_user_getvar
wrf_user_interp_level wrf.vinterp或wrf.interplevel 将3D数据场插值到某个/些垂直层上
wrf_user_vert_cross wrf.vertcross 将3D数据场插值到某个水平线上
wrf_user_interp_line wrf.interpline 将2D数据场插值到某个线上
wrf_user_ll_to_xy wrf.ll_to_xy Lat/Lon <-> XY Routines

获取文件中的诊断变量：

1
2
3
4
5
6
7
8
9
10
from __future__ import print_function

from netCDF4 import Dataset
from wrf import getvar

ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")

p = getvar(ncfile, "P")

print(p)

常用的一些诊断变量：

varname Description
lat latitude
lon longitude
p/pres Full model pressure [Pa]
pressure Full model pressure [hPa]
z/height Full model height [m]
geopt/geopotential
pvo potential vorticity [PVU]
slp Sea level pressure [hPa]
ter Model terrain height [m]
tc Temperature [C]
tk Temperature [K]
th/theta Potential temperature [K]
times Times in file (return strings - recommended)
Times Times in file (return characters)
ua U component of wind on mass points
va V component of wind on mass points
wa W component of wind on mass points
uvmet10 10m U and V components of wind rotated to earth coordinates
uvmet U and V components of wind rotated to earth coordinates

以数组的方式获取变量：disable xarray

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from __future__ import print_function

from netCDF4 import Dataset
from wrf import getvar, disable_xarray

ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00")

# Disable xarray completely
disable_xarray()
p_no_meta = getvar(ncfile, "P")
print (type(p_no_meta))
enable_xarray()

# Disable by using the meta parameter
p_no_meta = getvar(ncfile, "P", meta=False)
print (type(p_no_meta))

读取多时次单文件

如果一个wrfout文件中存储了多个时次，读取的时候如何获得完成的数据：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from netCDF4 import Dataset
from wrf import getvar,disable_xarray,ALL_TIMES


ncfile = Dataset(file_dir+"wrfout_d01_2013-07-15_00_00_00")
disable_xarray()

# 获取数据坐标
lon = getvar(ncfile,"longitude")
lat = getvar(ncfile,"latitude")

# 获得文件中第三个时刻的数据
p = getvar(ncfile,"pressure",timeidx=2)

# 获得文件中所有时刻的数据
p = getvar(ncfile,"pressure",ALL_TIMES)

# Get the Sea Level Pressure
slp = getvar(ncfile, "slp")
z = getvar(ncfile, "z",ALL_TIMES)
# t = getvar(ncfile,'tc')
t2 = getvar(ncfile,'T2',ALL_TIMES)-273.15
hours = getvar(ncfile, "times",ALL_TIMES)
# print(hours)

批量读取单时次文件并拼接

用Python批处理指定数据-以WRF输出结果为例演示按照指定维度合并(附示例代码)-腾讯云开发者社区-腾讯云

1

画图

画全球图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
proj = ccrs.PlateCarree(central_longitude=0)
fig = plt.figure(figsize=[5.4,2.5],dpi=200)
ax = plt.axes(projection=proj)
ax.set_position([0.06,0.1,0.87,0.87])#有colorbar绘图（下端）
ax.coastlines(linewidth=0.25)
cycle_var, cycle_lon = add_cyclic_point(var, coord=lon)

ax.set_extent([-180,180,-90,90], crs=ccrs.PlateCarree())#一定要有！
ax.set_xticks([0,60,120,180,-180,-120,-60])#指定要显示的经纬度                      
ax.set_yticks([-90,-60,-30,0,30,60,90])  
# ax.set_yticks([-80,-70,-50,-40,-20,-10],minor = True)

ax.xaxis.set_major_formatter(LongitudeFormatter())#刻度格式转换为经纬度样式                       
ax.yaxis.set_major_formatter(LatitudeFormatter())                        
ax.tick_params(axis='both',which='major',labelsize=10,direction='out'
                ,length=3,width=0.3,pad=3,top=True,right=True)
yminorLocator = MultipleLocator(10)
# import matplotlib
# ax.yaxis.set_minor_locator(matplotlib.ticker.FixedLocator(yminorLocator))
ax.yaxis.set_minor_locator(yminorLocator)
xminorLocator = MultipleLocator(10)
ax.xaxis.set_minor_locator(xminorLocator)

ax.tick_params(top='on',right='on',which='both') 
ax.tick_params(axis='both',which='minor',direction='out'
                ,length=1.5,width=0.3,top=True,right=True)
ax.spines['geo'].set_linewidth(0.2)#调节边框粗细


levels = np.around(np.arange(-1.6,2.1,0.4),2) #bsf_diff
cs = ax.contourf(cycle_lon,lat, cycle_var[:,:],
                  levels=levels,cmap=nclcmaps.cmap('GMT_panoply'),
                  #MPL_YlOrRd,NCV_blue_red,BlueYellowRed,MPL_RdYlBu
                  #BlueDarkRed18,BlueWhiteOrangeRed
                  # ts:NCV_bright
                  # psi:ncl_default
                  # ss: NCV_bright
                  # bsf:GMT_panoply
                  extend='both')

v = np.around(np.arange(-1.6,2.1,0.4),2) #bsf_diff


cax = fig.add_axes([0.92,0.1,0.025,0.87])
bar = plt.colorbar(cs,cax = cax,orientation='vertical',
                    ticks=v,drawedges=False)

bar.ax.tick_params(labelcolor='k',labelsize=8,pad=0.5)
bar.outline.set_linewidth(0.3)
bar.ax.set_yticklabels(v,fontsize=8,fontname="Times New Roman")

fig.text(0.935,0.03,"(Sv)",fontsize=10) # bsf


plt.savefig(fp_0+'_yearmean89_'+var_name+"_diff2.png",dpi=800)

画wrf domain

Cartopy 系列：为 Lambert 投影地图添加刻度 - 炸鸡人博客

wrf模拟的domain图绘制 - xiaofeifeixd - 博客园

Python画wrf模式的模拟区域

Cartopy 系列：探索 shapefile - 炸鸡人博客

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
from __future__ import print_function
import xarray as xr
import numpy as np
from netCDF4 import Dataset
import netCDF4 as nc
from wrf import getvar,disable_xarray,ALL_TIMES
from datetime import datetime, timedelta
import wrf
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy.ma as ma
import matplotlib
from global_land_mask import globe
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
from matplotlib.pyplot import MultipleLocator
import nclcmaps
# from datetime import datetime, timedelta
import os
import mod_lambert
import shapely.geometry as sgeom

plt.rcParams['font.sans-serif']=['Times New Roman']
matplotlib.rcParams.update({'font.size': 8})
file_dir = 'D:\\LHwork\\wrfout\\'

ncfile = Dataset(file_dir+"met_em.d01.2022-08-04_00_00_00.nc")

disable_xarray()

hgt = getvar(ncfile,"HGT_M")
landmask = (getvar(ncfile,"LANDMASK"))
b = landmask==0

hgt = ma.masked_array(hgt,mask = b)
lon = getvar(ncfile,"longitude")
lat = getvar(ncfile,"latitude")

proj = ccrs.PlateCarree()  # 创建坐标系
ref_lat = 38
ref_lon = 107
true_lat1 = 30
true_lat2 = 60
false_easting = (216-1)/2*30000
false_northing = (179.5-1)/2*30000

proj_lambert = ccrs.LambertConformal(
                central_longitude=ref_lon,
                central_latitude=ref_lat,
                standard_parallels=(true_lat1,true_lat2),
                cutoff=-30,
                false_easting=false_easting,
                false_northing=false_northing,
            )
## 创建坐标系
fig = plt.figure(figsize=(4.5,3.2), dpi=200)  # 创建页面
ax = fig.add_axes([0.04,0.07,0.85,0.88], projection=proj_lambert)  
ax.set_extent([0, false_easting*2, 0, false_northing*2], crs=proj_lambert)
ax.coastlines(linewidth=0.8,resolution='50m')
ax.add_feature(cfeature.LAKES.with_scale('50m'),linewidth=0.5,edgecolor='black',facecolor='none')

# 添加shp文件
import cartopy.io.shapereader as shpreader
filepath ='D:\LHwork\shpfile\中国GIS地图\国家基础地理数据/hyd1_4l.shp'
crs = ccrs.PlateCarree()
reader = shpreader.Reader(filepath)
geoms = reader.geometries()
ax.add_geometries(geoms, crs, lw=0.5, fc='none')
reader.close()

## --设置网格属性, 不画默认的标签
gl=ax.gridlines(draw_labels=True,linestyle=":",linewidth=0.6 ,x_inline=False, y_inline=False,color='k')

## 关闭上面和右边的经纬度显示
gl.top_labels=False #关闭上部经纬标签                                  
# gl.bottom_labels = False
# gl.left_labels = False
import matplotlib.ticker as mticker
gl.right_labels=False
gl.xformatter = LONGITUDE_FORMATTER  #使横坐标转化为经纬度格式            
gl.yformatter = LATITUDE_FORMATTER        

gl.xlocator=mticker.FixedLocator(np.arange(60,160,10))                              
gl.ylocator=mticker.FixedLocator(np.arange(10,60,10)) 
gl.rotate_labels = False
gl.xlabel_style={'size':10,}
                  # 'rotation':'horizontal','rotation_mode':'default',
                  # 'horizontalalignment':'center','verticalalignment':'center'}#修改经纬度字体大小                             
gl.ylabel_style={'size':10,}#'rotation':'horizontal'}
ax.spines['geo'].set_linewidth(0.6)#调节边框粗细

plt.plot([110,122], [28,28], 'r--',linewidth=1,transform=ccrs.PlateCarree(),zorder=5)
plt.plot([110,122], [32,32], 'r--',linewidth=1,transform=ccrs.PlateCarree(),zorder=5)
plt.plot([110,110], [28,32], 'r--',linewidth=1,transform=ccrs.PlateCarree(),zorder=5)
plt.plot([122,122], [28,32], 'r--',linewidth=1,transform=ccrs.PlateCarree(),zorder=5)

levels = np.around(np.arange(0,5601,400),2) #bsf_diff
cs = ax.contourf(lon,lat, hgt,transform=ccrs.PlateCarree(),
                  levels=levels,cmap= nclcmaps.cmap('OceanLakeLandSnow'),#'terrain', 'OceanLakeLandSnow'
                  extend='both')
over_t = (224/255,224/255,224/255)
under_t = (0, 71/255 , 71/255 )
under_t = (1,1,1)
cs.cmap.set_over(over_t) #'fuchsia'
cs.cmap.set_under(under_t) #21
cs.changed()

v = np.around(np.arange(0,5601,400),2) #bsf_diff

cax = fig.add_axes([0.88,0.07,0.03,0.88])
bar = plt.colorbar(cs,cax = cax,orientation='vertical',
                    ticks=v,drawedges=False)

bar.ax.tick_params(labelcolor='k',labelsize=8,pad=0.5)
bar.outline.set_linewidth(0.3)
bar.ax.set_yticklabels(v,fontsize=9,fontname="Times New Roman")

fig.text(0.92,0.05,"(m)",fontsize=10) # bsf

plt.savefig("hgt.png",dpi=800)

画多子图

Cartopy地图投影绘制 | ZSYXY Meteorological workshop (yxy-biubiubiu.github.io)

运算

插值到规则经纬网格

将WRF模拟结果与站点观测或者再分析资料进行对比之前，我们需要对WRF输出的网格资料进行插值（或regrid）

利用xesmf将WRF数据插值到站点/网格

1
2
3
4
5
6
7
8
9
10
import salem
import xarray as xr 
import xesmf as xe
import numpy as np

lons = [100 , 120]
lats = [30  ,  40]
step = 1
ds_out_2D = xe.util.grid_2d(lons[0]-step/2, lons[1]+step/2, step, 
                        lats[0]-step/2, lats[1]+step/2, step)

python—站点数据插值到规则网格&EOF分解

1
2
3
4
5
6
from scipy.interpolate import griddata


z[i] = griddata((lon_real,lat_real),
                                    data_all[:,i],
                                    (x_t,y_t),method='cubic')

注意，使用以上这些插值的效果，特别是scipy中的插值效果都不如ncl的rcm2rgrid_wrap命令，所以考虑找到这个命令在python中的替代，发现了geocat项目。

Geocat

Installation — GeoCAT-f2py = “2023.03.0” %} documentation

Installation (geocat-comp.readthedocs.io)

不知道哪个可以用，就都安装试一下，但是貌似只能在linux下使用，在win下是不行的。

使用了，goecat.comp以及被淘汰，目前新的在使用的是geocat.f2py

wrfout文件解读（按照开头字母）

https://bbs.06climate.com/forum.php?mod=viewthread&tid=13277

float LU_INDEX(Time, south_north, west_east) ;
LU_INDEX:description = “LAND USE CATEGORY” ;
LU_INDEX:units = “” ;
土地利用类型（如城市、植被、湖泊等）

float ZNU(Time, bottom_top) ;
ZNU:description = “eta values on half (mass) levels” ;
ZNU:units = “” ;
eta半层（质量点）坐标值

float ZNW(Time, bottom_top_stag) ;
ZNW:description = “eta values on full (w) levels” ;
ZNW:units = “” ;
eta整层（w点）坐标值

float ZS(Time, soil_layers_stag) ;
ZS:description = “DEPTHS OF CENTERS OF SOIL LAYERS” ;
ZS:units = “m” ;
土壤层各层中间的深度

float DZS(Time, soil_layers_stag) ;
DZS:description = “THICKNESSES OF SOIL LAYERS” ;
DZS:units = “m” ;
土壤层厚度

float U(Time, bottom_top, south_north, west_east_stag) ;
U:description = “x-wind component” ;
U:units = “m s-1” ;
x方向风分量

float V(Time, bottom_top, south_north_stag, west_east) ;
V:description = “y-wind component” ;
V:units = “m s-1” ;
y方向风分量

float W(Time, bottom_top_stag, south_north, west_east) ;
W:description = “z-wind component” ;
W:units = “m s-1” ;
垂直风分量

float PH(Time, bottom_top_stag, south_north, west_east) ;
PH:description = “perturbation geopotential” ;
PH:units = “m2 s-2” ;
扰动位势高度

float PHB(Time, bottom_top_stag, south_north, west_east) ;
PHB:description = “base-state geopotential” ;
PHB:units = “m2 s-2” ;
基准态位势高度

float T(Time, bottom_top, south_north, west_east) ;
T:description = “perturbation potential temperature (theta-t0)” ;
T:units = “K” ;
扰动位温（theta-t0）

float MU(Time, south_north, west_east) ;
MU:description = “perturbation dry air mass in column” ;
MU:units = “Pa” ;
柱内扰动干空气质量

float MUB(Time, south_north, west_east) ;
MUB:description = “base state dry air mass in column” ;
MUB:units = “Pa” ;
柱内基准态干空气质量

float NEST_POS(Time, south_north, west_east) ;
NEST_POS:description = “-“ ;
NEST_POS:units = “-“ ;
作多层嵌套时粗（母）网格位置

float P(Time, bottom_top, south_north, west_east) ;
P:description = “perturbation pressure” ;
P:units = “Pa” ;
扰动气压

float PB(Time, bottom_top, south_north, west_east) ;
PB:description = “BASE STATE PRESSURE” ;
PB:units = “Pa” ;
基准态气压

float SR(Time, south_north, west_east) ;
SR:description = “fraction of frozen precipitation” ;
SR:units = “-“ ;
固体降水比例

float POTEVP(Time, south_north, west_east) ;
POTEVP:description = “accumulated potential evaporation” ;
POTEVP:units = “W m-2” ;
累计的潜在蒸发

float SNOPCX(Time, south_north, west_east) ;
SNOPCX:description = “snow phase change heat flux” ;
SNOPCX:units = “W m-2” ;
雪相态改变的热通量

float SOILTB(Time, south_north, west_east) ;
SOILTB:description = “bottom soil temperature” ;
SOILTB:units = “K” ;
土壤底部温度

float FNM(Time, bottom_top) ;
FNM:description = “upper weight for vertical stretching” ;
FNM:units = “” ;
作垂直方向展开时上层权重

float FNP(Time, bottom_top) ;
FNP:description = “lower weight for vertical stretching” ;
FNP:units = “” ; 作垂直方向展开时下层权重

float RDNW(Time, bottom_top) ;
RDNW:description = “inverse d(eta) values between full (w) levels” ;
RDNW:units = “” ;
1除以整层（w层）间eta值之差

float RDN(Time, bottom_top) ;
RDN:description = “inverse d(eta) values between half (mass) levels” ;
RDN:units = “” ;
1除以半层（质量层）间eta值之差

float DNW(Time, bottom_top) ;
DNW:description = “d(eta) values between full (w) levels” ;
DNW:units = “” ;
整层（w层）间eta值之差

float DN(Time, bottom_top) ;
DN:description = “d(eta) values between half (mass) levels” ;
DN:units = “” ;
半层（质量层）间eta值之差

float CFN(Time) ;
CFN:description = “extrapolation constant” ;
CFN:units = “” ;
外推常数

float CFN1(Time) ;
CFN1:description = “extrapolation constant” ;
CFN1:units = “” ;
外推常数1

float Q2(Time, south_north, west_east) ;
Q2:description = “QV at 2 M” ;
Q2:units = “kg kg-1” ;
地面上2m高度的比湿

float T2(Time, south_north, west_east) ;
T2:description = “TEMP at 2 M” ;
T2:units = “K” ;
地面上2m高度的温度

float TH2(Time, south_north, west_east) ;
TH2:description = “POT TEMP at 2 M” ;
TH2:units = “K” ;
地面上2m高度的位温

float PSFC(Time, south_north, west_east) ;
PSFC:description = “SFC PRESSURE” ;
PSFC:units = “Pa” ; 地面气压

float U10(Time, south_north, west_east) ;
U10:description = “U at 10 M” ;
U10:units = “m s-1” ;
地面上10m风场的纬向分量

float V10(Time, south_north, west_east) ;
V10:description = “V at 10 M” ;
V10:units = “m s-1” ;
地面上10m风场的经向分量

float RDX(Time) ;
RDX:description = “INVERSE X GRID LENGTH” ;
RDX:units = “” ;
1除以x方向网格距

float RDY(Time) ;
RDY:description = “INVERSE Y GRID LENGTH” ;
RDY:units = “” ;
1除以y方向网格距

float RESM(Time) ;
RESM:description = “TIME WEIGHT CONSTANT FOR SMALL STEPS” ;
RESM:units = “” ;
对小时间步长时的时间权重常数

float ZETATOP(Time) ;
ZETATOP:description = “ZETA AT MODEL TOP” ;
ZETATOP:units = “” ;
模式大气层顶的eta值

float CF1(Time) ;
CF1:description = “2nd order extrapolation constant” ;
CF1:units = “” ;
2阶外推常数1

float CF2(Time) ;
CF2:description = “2nd order extrapolation constant” ;
CF2:units = “” ;
2阶外推常数2

float CF3(Time) ;
CF3:description = “2nd order extrapolation constant” ;
CF3:units = “” ;
2阶外推常数3

int ITIMESTEP(Time) ;
ITIMESTEP:description = “” ;
ITIMESTEP:units = “” ;
时间步长计数

float XTIME(Time) ;
XTIME:description = “minutes since simulation start” ;
XTIME:units = “” ;
已经模拟时间的长度

float QVAPOR(Time, bottom_top, south_north, west_east) ;
QVAPOR:description = “Water vapor mixing ratio” ;
QVAPOR:units = “kg kg-1” ;
水汽混合比

float QCLOUD(Time, bottom_top, south_north, west_east) ;
QCLOUD:description = “Cloud water mixing ratio” ;
QCLOUD:units = “kg kg-1” ;
云水混合比

float QRAIN(Time, bottom_top, south_north, west_east) ;
QRAIN:description = “Rain water mixing ratio” ;
QRAIN:units = “kg kg-1” ;
雨水混合比

float LANDMASK(Time, south_north, west_east) ;
LANDMASK:description = “LAND MASK (1 FOR LAND, 0 FOR WATER)” ;
LANDMASK:units = “” ;
陆面指标（1是陆地，0是水体）

float TSLB(Time, soil_layers_stag, south_north, west_east) ;
TSLB:description = “SOIL TEMPERATURE” ;
TSLB:units = “K” ;
各层土壤温度

float SMOIS(Time, soil_layers_stag, south_north, west_east) ;
SMOIS:description = “SOIL MOISTURE” ;
SMOIS:units = “m3 m-3” ;
各层土壤湿度

float SH2O(Time, soil_layers_stag, south_north, west_east) ;
SH2O:description = “SOIL LIQUID WATER” ;
SH2O:units = “m3 m-3” ;
各层土壤液态水含量？

float SEAICE(Time, south_north, west_east) ;
SEAICE:description = “SEA ICE FLAG” ;
SEAICE:units = “” ;
海冰标志

float XICEM(Time, south_north, west_east) ;
XICEM:description = “SEA ICE FLAG (PREVIOUS STEP)” ;
XICEM:units = “” ;
上一步时的海冰标志

float SFROFF(Time, south_north, west_east) ;
SFROFF:description = “SURFACE RUNOFF” ;
SFROFF:units = “mm” ;
地表径流

float UDROFF(Time, south_north, west_east) ;
UDROFF:description = “UNDERGROUND RUNOFF” ;
UDROFF:units = “mm” ;
地下径流

int IVGTYP(Time, south_north, west_east) ;
IVGTYP:description = “DOMINANT VEGETATION CATEGORY” ;
IVGTYP:units = “” ;
占主导的植被种类

int ISLTYP(Time, south_north, west_east) ;
ISLTYP:description = “DOMINANT SOIL CATEGORY” ;
ISLTYP:units = “” ;
占主导的土壤种类

float VEGFRA(Time, south_north, west_east) ;
VEGFRA:description = “VEGETATION FRACTION” ;
VEGFRA:units = “” ;
植被比例

float GRDFLX(Time, south_north, west_east) ;
GRDFLX:description = “GROUND HEAT FLUX” ;
GRDFLX:units = “W m-2” ;
地面热通量

float SNOW(Time, south_north, west_east) ;
SNOW:description = “SNOW WATER EQUIVALENT” ;
SNOW:units = “kg m-2” ;
雪水当量

float SNOWH(Time, south_north, west_east) ;
SNOWH:description = “PHYSICAL SNOW DEPTH” ;
SNOWH:units = “m” ;
实际雪厚

float RHOSN(Time, south_north, west_east) ;
RHOSN:description = “ SNOW DENSITY” ;
RHOSN:units = “kg m-3” ;
雪密度

float CANWAT(Time, south_north, west_east) ;
CANWAT:description = “CANOPY WATER” ;
CANWAT:units = “kg m-2” ;
冠层中的水

float SST(Time, south_north, west_east) ;
SST:description = “SEA SURFACE TEMPERATURE” ;
SST:units = “K” ;
海表面温度

float QNDROPSOURCE(Time, bottom_top, south_north, west_east) ;
QNDROPSOURCE:description = “Droplet number source” ;
QNDROPSOURCE:units = “ /kg/s” ;
水滴数源

float MAPFAC_M(Time, south_north, west_east) ;
MAPFAC_M:description = “Map scale factor on mass grid” ;
MAPFAC_M:units = “” ;
质量格点处的地图比例系数

float MAPFAC_U(Time, south_north, west_east_stag) ;
MAPFAC_U:description = “Map scale factor on u-grid” ;
MAPFAC_U:units = “” ;
u-格点处的地图比例系数

float MAPFAC_V(Time, south_north_stag, west_east) ;
MAPFAC_V:description = “Map scale factor on v-grid” ;
MAPFAC_V:units = “” ;
v-格点处的地图比例系数

float MAPFAC_MX(Time, south_north, west_east) ;
MAPFAC_MX:description = “Map scale factor on mass grid, x direction” ;
MAPFAC_MX:units = “” ;
x方向上质量格点处的地图比例系数

float MAPFAC_MY(Time, south_north, west_east) ;
MAPFAC_MY:description = “Map scale factor on mass grid, y direction” ;
MAPFAC_MY:units = “” ;
y方向上质量格点处的地图比例系数

float MAPFAC_UX(Time, south_north, west_east_stag) ;
MAPFAC_UX:description = “Map scale factor on u-grid, x direction” ;
MAPFAC_UX:units = “” ;
x方向上u-格点处的地图比例系数

float MAPFAC_UY(Time, south_north, west_east_stag) ;
MAPFAC_UY:description = “Map scale factor on u-grid, y direction” ;
MAPFAC_UY:units = “” ;
y方向上u-格点处的地图比例系数

float MAPFAC_VX(Time, south_north_stag, west_east) ;
MAPFAC_VX:description = “Map scale factor on v-grid, x direction” ;
MAPFAC_VX:units = “” ;
x方向上v-格点处的地图比例系数

float MF_VX_INV(Time, south_north_stag, west_east) ;
MF_VX_INV:description = “Inverse map scale factor on v-grid, x direction”
MF_VX_INV:units = “” ;
1除以x方向上v-格点处的地图比例系数

float MAPFAC_VY(Time, south_north_stag, west_east) ;
MAPFAC_VY:description = “Map scale factor on v-grid, y direction” ;
MAPFAC_VY:units = “” ;
y方向上v-格点处的地图比例系数

float F(Time, south_north, west_east) ;
F:description = “Coriolis sine latitude term” ;
F:units = “s-1” ;
科氏力中sin(Ω)的部分，Ω为纬度

float E(Time, south_north, west_east) ;
E:description = “Coriolis cosine latitude term” ;
E:units = “s-1” ;
科氏力中cos(Ω)的部分，Ω为纬度

float SINALPHA(Time, south_north, west_east) ;
SINALPHA:description = “Local sine of map rotation” ;
SINALPHA:units = “” ;
局地sin（地图旋转角）

float COSALPHA(Time, south_north, west_east) ;
COSALPHA:description = “Local cosine of map rotation” ;
COSALPHA:units = “” ;
局地cos（地图旋转角）

float HGT(Time, south_north, west_east) ;
HGT:description = “Terrain Height” ;
HGT:units = “m” ;
地形高度

float HGT_SHAD(Time, south_north, west_east) ;
HGT_SHAD:description = “Height of orographic shadow” ;
HGT_SHAD:units = “m” ;
山岳背光坡的高度

float TSK(Time, south_north, west_east) ;
TSK:description = “SURFACE SKIN TEMPERATURE” ;
TSK:units = “K” ;
地表温度

float P_TOP(Time) ;
P_TOP:description = “PRESSURE TOP OF THE MODEL” ;
P_TOP:units = “Pa” ;
模式顶的气压

float MAX_MSTFX(Time) ;
MAX_MSTFX:description = “Max map factor in domain” ;
MAX_MSTFX:units = “” ;
区域内的最大地图比例系数

float RAINC(Time, south_north, west_east) ;
RAINC:description = “ACCUMULATED TOTAL CUMULUS PRECIPITATION” ;
RAINC:units = “mm” ;
累积的积云对流降水

float RAINNC(Time, south_north, west_east) ;
RAINNC:description = “ACCUMULATED TOTAL GRID SCALE PRECIPITATION” ;
RAINNC:units = “mm” ;
累积的格点降水

float PRATEC(Time, south_north, west_east) ;
PRATEC:description = “PRECIP RATE FROM CUMULUS SCHEME” ;
PRATEC:units = “mm s-1” ;
由积云方案算的对流降水率

float RAINCV(Time, south_north, west_east) ;
RAINCV:description = “TIME-STEP CUMULUS PRECIPITATION” ;
RAINCV:units = “mm” ;
计算对流降水的时间步长

float SNOWNC(Time, south_north, west_east) ;
SNOWNC:description = “ACCUMULATED TOTAL GRID SCALE SNOW AND ICE” ;
SNOWNC:units = “mm” ;
累积的格点降雪和冰量

float GRAUPELNC(Time, south_north, west_east) ;
GRAUPELNC:description = “ACCUMULATED TOTAL GRID SCALE GRAUPEL” ;
GRAUPELNC:units = “mm” ;
累积的格点降雪丸量

float EDT_OUT(Time, south_north, west_east) ;
EDT_OUT:description = “EDT FROM GD SCHEME” ;
EDT_OUT:units = “” ;
Grell-Devenyi方案计算的edt场

float SWDOWN(Time, south_north, west_east) ;
SWDOWN:description = “DOWNWARD SHORT WAVE FLUX AT GROUND SURFACE” ;
SWDOWN:units = “W m-2” ;
地表处的向下短波辐射通量

float GLW(Time, south_north, west_east) ;
GLW:description = “DOWNWARD LONG WAVE FLUX AT GROUND SURFACE” ;
GLW:units = “W m-2” ;
地表处的向下长波辐射通量

float OLR(Time, south_north, west_east) ;
OLR:description = “TOA OUTGOING LONG WAVE” ;
OLR:units = “W m-2” ;
TOA处的向上长波辐射通量

float XLAT(Time, south_north, west_east) ;
XLAT:description = “LATITUDE, SOUTH IS NEGATIVE” ;
XLAT:units = “degree_north” ;
质量点的纬度（南半球为负值）

float XLONG(Time, south_north, west_east) ;
XLONG:description = “LONGITUDE, WEST IS NEGATIVE” ;
XLONG:units = “degree_east” ;
质量点的经度（西半球为负值）

float XLAT_U(Time, south_north, west_east_stag) ;
XLAT_U:description = “LATITUDE, SOUTH IS NEGATIVE” ;
XLAT_U:units = “degree_north” ;
U-格点的纬度（南半球为负值）

float XLONG_U(Time, south_north, west_east_stag) ;
XLONG_U:description = “LONGITUDE, WEST IS NEGATIVE” ;
XLONG_U:units = “degree_east” ;
U-格点的经度（西半球为负值）

float XLAT_V(Time, south_north_stag, west_east) ;
XLAT_V:description = “LATITUDE, SOUTH IS NEGATIVE” ;
XLAT_V:units = “degree_north” ;
V-格点的纬度（南半球为负值）

float XLONG_V(Time, south_north_stag, west_east) ;
XLONG_V:description = “LONGITUDE, WEST IS NEGATIVE” ;
XLONG_V:units = “degree_east” ;
V-格点的经度（西半球为负值）

float ALBEDO(Time, south_north, west_east) ;
ALBEDO:description = “ALBEDO” ;
ALBEDO:units = “-“ ;
反照率

float ALBBCK(Time, south_north, west_east) ;
ALBBCK:description = “BACKGROUND ALBEDO” ;
ALBBCK:units = “” ;
背景反照率

float EMISS(Time, south_north, west_east) ;
EMISS:description = “SURFACE EMISSIVITY” ;
EMISS:units = “” ;
地面发射率

float TMN(Time, south_north, west_east) ;
TMN:description = “SOIL TEMPERATURE AT LOWER BOUNDARY” ;
TMN:units = “K” ;
更低边界处的土壤温度

float XLAND(Time, south_north, west_east) ;
XLAND:description = “LAND MASK (1 FOR LAND, 2 FOR WATER)” ;
XLAND:units = “” ;
陆面指标（1是陆地，2是水体）

float UST(Time, south_north, west_east) ;
UST:description = “U* IN SIMILARITY THEORY” ;
UST:units = “m s-1” ;
相似性理论中的摩擦速度

float PBLH(Time, south_north, west_east) ;
PBLH:description = “PBL HEIGHT” ;
PBLH:units = “m” ;
行星边界层高度

float HFX(Time, south_north, west_east) ;
HFX:description = “UPWARD HEAT FLUX AT THE SURFACE” ;
HFX:units = “W m-2” ;
地表面处向上的热量通量（感热通量）

float QFX(Time, south_north, west_east) ;
QFX:description = “UPWARD MOISTURE FLUX AT THE SURFACE” ;
QFX:units = “kg m-2 s-1” ;
地表面处向上的水汽通量

float LH(Time, south_north, west_east) ;
LH:description = “LATENT HEAT FLUX AT THE SURFACE” ;
LH:units = “W m-2” ;
地表面处的潜热通量

float SNOWC(Time, south_north, west_east) ;
SNOWC:description = “FLAG INDICATING SNOW COVERAGE (1 FOR SNOW COVER)” ;
SNOWC:units = “” ;
雪盖标志（1是有雪）

降水处理

WRF后处理：降雨量的说明以及降雨的绘制_wrf模拟降水量偏小-CSDN博客

1
2
3
RAINC: ACCUMULATED TOTAL CUMULUS PRECIPITATION
RAINNC: ACCUMULATED TOTAL GRID SCALE PRECIPITATION
RAINSH: ACCUMULATED SHALLOW CUMULUS PRECIPITATION

RAINC: 积云深对流过程产生的累积降水量，也就是模式中的积云对流参数化方案导致的降雨( cu_physics)。对于高分辨率的模拟，比如dx<5km，通常会将积云对流参数化方案关闭，此时RAINC为0。

RANNC: 此类降雨来源于云微物理参数化方案(mp_physics)，如大尺度抬升过程产生的凝结等微物理过程降水，也就是非对流产生的降水。

RAINSH: 积云对流参数化方案主要是反映深对流的降水过程，但是一些积云对流参数化方案，能够支持浅对流导致的降水，此时总降水还需要加上RAINSH。WRF中支持浅对流的参数化方案(cu_physics)有以下几种：KF，SAS，G3，BMJ，Tiedtke。WRF中也有独立于深对流过程的浅对流方案，通过namelist中设置shcu_physics。一般情况下，浅对流产生的降水量较小。

云量

Total cloud fraction (definition & calculation method) | WRF & MPAS-A Support Forum (ucar.edu)

wrf.cloudfrac — wrf-python 1.3.4.1 documentation

The UPP (Unified Post-Processing) column max method is a technique used to calculate the total cloud fraction from atmospheric model output files, such as WRF (Weather Research and Forecasting) model output. Let’s break down how it works:

  1. Objective:

    • The goal is to determine the total cloud fraction within each grid cell of the model domain.

  2. Vertical Levels:

    • Atmospheric models represent the atmosphere in vertical layers (levels). These levels can be categorized as low, mid, and high.
    • For cloud fraction calculations, we consider these vertical levels.

  3. Cloud Fraction Variable:

    • In WRF output files (e.g., wrfout_d0*), there is a variable called CLDFRA (cloud fraction).
    • This variable provides the fraction of each grid cell covered by clouds at different vertical levels.

  4. Methodology:

    • For each grid cell, we examine the cloud fraction values at all vertical levels.
    • The UPP column max method selects the maximum cloud fraction value across all vertical levels within that grid cell.
    • This maximum value represents the total cloud fraction for that specific grid cell.

  5. Calculation:

    • Suppose we have cloud fraction values for low, mid, and high levels (denoted as CLDFRA_low, CLDFRA_mid, and CLDFRA_high).

    • The total cloud fraction (CLDFRA_total) is calculated as:

      CLDFRAtotal​=max(CLDFRAlow​,CLDFRAmid​,CLDFRAhigh​)

  6. Application:

    • Researchers and meteorologists use this method to analyze cloud cover patterns, study cloud processes, and assess model performance.
    • It’s essential to consider the limitations and assumptions associated with this approach when interpreting results.

单位为mm的变量

205

1
GRAUPELNC

206

1
HAILNC

201 

1
RAINC

203

1
RAINNC

202

1
RAINSH

150

1
SFROFF

204

1
SNOWNC

151

1
UDROFF

instantaneous & accumulated

instantaneous output
accumulation over hours 0-3 = 10
output at hr 03 = 10
accumulation over hours 3-6 = 5
output at hr 06 = 5
accumulation over hours 6-9 = 7
output at hr 09 = 7

accumulated output
accumulation over hours 0-3 = 10
output at hr 03 = 10
accumulation over hours 3-6 = 5
output at hr 06 = 15
accumulation over hours 6-9 = 7
output at hr 09 = 22

地表能量平衡的诊断

Confusing names for variables in wrfout | WRF & MPAS-A Support Forum (ucar.edu)

WRF surface net radiation | WRF & MPAS-A Support Forum (ucar.edu)

How to calculate Net radiation in WRF | WRF & MPAS-A Support Forum (ucar.edu)

Sensible heat flux | WRF & MPAS-A Support Forum (ucar.edu)

WRF energy and moisture budget? | WRF & MPAS-A Support Forum (ucar.edu)

n general, the following calculations can be used for radiative flux calculations:

$GSW=(1-ALBEDO)SWDOWN$: net shortwave radiative flux
GLW: downward longwave radiative flux
$LWUPFLUX = STBOLTTSK*4$: where STBOLT is the Stefan-Boltzmann constant
$netLW = EMISS(GLW-LWUPFLUX)$: net long wave radiative flux

net radiative flux: $RNET=GSW+netLW$

Surface energy budget can be computed as $RNET-GRDFLX-HFX-LH = 0$

世纪巨坑：

如果从陆面的角度拆解潜热和感热的话：

$HFX/SH=F_{veg}SHG+(1-F_{veg})SHB+SHC$

$LH=F_{veg}EVG+(1-F_{veg})EVB+EVC+TR$

其中$F_{veg}$为植被分数，读取VEGFRA变量并除以100.

HFX就是感热通量。