问题简介
我们都知道，在运行WRF模式之前，需要使用其WPS前处理程序，制作对应的模拟域与气象场输入，而WPS主要由3个主程序组成，分别为：
geogrid.exe定义模式的模拟域，并将静态地理数据插值到模式网格中。ungrib.exe从Grib的格式的文件读取气象数据。metgrid.exe将读取的气象数据，水平插值到geogrid.exe定义的模拟域网格中。
通过以上介绍不难看出，整个WPS的处理流程并不复杂，而最后我们的目的就是要将模拟域与气象场结合，输入到WRF中。
在运行WPS过程中，metgrid.exe需要最后运行，因为它是将其他两个主程序处理后的数据进行插值结合的程序，即：从geogrid.exe中获得模拟域与地理信息，从ungrib.exe中获得气象数据，而metgrid.exe从两个程序中获取信息的方式则为：读取geogrid.exe和ungrib.exe运行生成的中间文件（intermediate file），即geo_dem和 FILE打头的文件。
那么，就有了问题，如果我们想添加额外的气象数据，而这个气象数据本身并不支持grib格式，无法用ungrib.exe处理，我们应当怎么做呢？
这时候，我们就需要自己根据我们的数据，撰写可以被metgrid.exe所读取的中间文件，并修改metgrid.TBL文件和namelist.wps，使得我们添加的数据与变量可以被metgrid.exe处理。
中间文件可分为三种格式：WPS、WRFSI、MM5。而WRF-ARW-online网站提供了三种格式的信息与描述： intermediate file这里主要介绍最简单与最常用的WPS FROMAT。
WPS intermediate file
在将数据写入WPS中间格式时，二维字段被写入实值的矩形数组。三维数组必须跨垂直维度拆分为独立编写的二维数组。还应注意的是，对于全局数据集，必须使用高斯投影或柱面等距投影，对于区域数据集，可以使用墨卡托投影、兰伯特共形投影、极坐标投影或柱面等距投影。
WRF的官网手册第三章，给出了如何将单个二维数组写成中间格式的Fortran代码：
integer :: version ! Format version (must =5 for WPS format)
integer :: nx, ny ! x- and y-dimensions of 2-d array
integer :: iproj ! Code for projection of data in array:
! 0 = cylindrical equidistant
! 1 = Mercator
! 3 = Lambert conformal conic
! 4 = Gaussian (global only!)
! 5 = Polar stereographic
real :: nlats ! Number of latitudes north of equator
! (for Gaussian grids)
real :: xfcst ! Forecast hour of data
real :: xlvl ! Vertical level of data in 2-d array
real :: startlat, startlon ! Lat/lon of point in array indicated by
! startloc string
real :: deltalat, deltalon ! Grid spacing, degrees
real :: dx, dy ! Grid spacing, km
real :: xlonc ! Standard longitude of projection
real :: truelat1, truelat2 ! True latitudes of projection
real :: earth_radius ! Earth radius, km
real, dimension(nx,ny) :: slab ! The 2-d array holding the data
logical :: is_wind_grid_rel ! Flag indicating whether winds are
! relative to source grid (TRUE) or
! relative to earth (FALSE)
character (len=8) :: startloc ! Which point in array is given by
! startlat/startlon; set either
! to 'SWCORNER' or 'CENTER '
character (len=9) :: field ! Name of the field
character (len=24) :: hdate ! Valid date for data YYYY:MM:DD_HH:00:00
character (len=25) :: units ! Units of data
character (len=32) :: map_source ! Source model / originating center
character (len=46) :: desc ! Short description of data
! 1) WRITE FORMAT VERSION
write(unit=ounit) version
! 2) WRITE METADATA
! Cylindrical equidistant
if (iproj == 0) then
write(unit=ounit) hdate, xfcst, map_source, field, &
units, desc, xlvl, nx, ny, iproj
write(unit=ounit) startloc, startlat, startlon, &
deltalat, deltalon, earth_radius
! Mercator
else if (iproj == 1) then
write(unit=ounit) hdate, xfcst, map_source, field, &
units, desc, xlvl, nx, ny, iproj
write(unit=ounit) startloc, startlat, startlon, dx, dy, &
truelat1, earth_radius
! Lambert conformal
else if (iproj == 3) then
write(unit=ounit) hdate, xfcst, map_source, field, &
units, desc, xlvl, nx, ny, iproj
write(unit=ounit) startloc, startlat, startlon, dx, dy, &
xlonc, truelat1, truelat2, earth_radius
! Gaussian
else if (iproj == 4) then
write(unit=ounit) hdate, xfcst, map_source, field, &
units, desc, xlvl, nx, ny, iproj
write(unit=ounit) startloc, startlat, startlon, &
nlats, deltalon, earth_radius
! Polar stereographic
else if (iproj == 5) then
write(unit=ounit) hdate, xfcst, map_source, field, &
units, desc, xlvl, nx, ny, iproj
write(unit=ounit) startloc, startlat, startlon, dx, dy, &
xlonc, truelat1, earth_radius
end if
! 3) WRITE WIND ROTATION FLAG
write(unit=ounit) is_wind_grid_rel
! 4) WRITE 2-D ARRAY OF DATA
write(unit=ounit) slab
可以看出在撰写WPS intermediate file时，我们主要需要知道和选择以下参数：dx,dy(维度)、投影、地球半径、预报时间、中间经纬度、垂直层以及一些对于变量的描述等等。而对于三维数组，我们只能将其拆成二维，一层一层的去写。
Fortran撰写WPS intermediate file固然方便，但需要注意的是，Fortran读取二进制的气象数据十分方便，但当读取nc、hdf5这种格式的数据时，就不太友好了，这是，我们可以使用NCL中的wrf_wps_write_int 函数。
ncl官网给出了函数的用法与例子，我们可以通过例子进行学习wrf-wps-write-int
确认中间格式是否正确：rd_intermediate.exe
当我们写出l了我们需要的中间文件时，如何确认我们格式是否正确呢？WPS也提供了rd_intermediate.exe程序帮助我们检验读取。
在/WPS/util文件夹下，输入命令：
./rd_intermediate.exe your intermediate file path >log
查看log文件，如果出现： SUCCESSFUL COMPLETION OF PROGRAM RD_INTERMEDIATE即表示格式正确。
下面是一个正确格式的log示例：
================================================
FIELD = GHT
UNITS = m DESCRIPTION = Height
DATE = 2018-03-01_12:00:00 FCST = 0.000000
SOURCE = ECMWF
LEVEL = 100000.000000
I,J DIMS = 1440, 181
IPROJ = 0 PROJECTION = LAT LON
REF_X, REF_Y = 1.000000, 1.000000
REF_LAT, REF_LON = 90.000008, -180.000015
DLAT, DLON = -0.250000, 0.250000
EARTH_RADIUS = 6367.470215
DATA(1,1)=254.153992
================================================
FIELD = TT
UNITS = K DESCRIPTION = Temperature
DATE = 2018-03-01_12:00:00 FCST = 0.000000
SOURCE = ECMWF
LEVEL = 100000.000000
I,J DIMS = 1440, 181
IPROJ = 0 PROJECTION = LAT LON
REF_X, REF_Y = 1.000000, 1.000000
REF_LAT, REF_LON = 90.000008, -180.000015
DLAT, DLON = -0.250000, 0.250000
EARTH_RADIUS = 6367.470215
DATA(1,1)=250.097992
================================================
SUCCESSFUL COMPLETION OF PROGRAM RD_INTERMEDIATE
修改METGRID.TBL
METGRID.TBL文件是一个文本文件，它定义了要用网格插值的每个气象场的参数。换句话说，只有参数在TBL文件中被定义，它才能被metgrid.exe文件读取插值。
每个字段的参数在单独的节中定义，节由一行等号分隔(例如，’ ============== ')。在每个部分中都有规范，每个规范的形式都是关键字=value。一些关键字在一个节中是必需的，而另一些则是可选的;有些关键字与其他关键字互斥。
关于METGRID.TBL的详细描述可参见官方手册；METGRID.TBL，下面是METGRID.TBL的一个示例：
========================================
name=SOIL_LEVELS
derived=yes #是否从其他变量处生成
z_dim_name=num_soilt_levels #垂直层名称
flag_in_output=FLAG_SOIL_LEVELS #全局属性输出名称
fill_lev=all:vertical_index; level_template=SOILT #指定插值层与方式
========================================
name=PRES
z_dim_name=num_metgrid_levels
derived=yes
mandatory=yes # MUST HAVE THIS FIELD
fill_lev=all:PRESSURE
fill_lev=200100:PSFC(200100)
fill_lev=all:vertical_index; level_template=TT
========================================
修改namelist.wps
写好了中间文件后，我们还必须修改namelist,wps的&metgrid部分，使得metgrid.exe读取我们的中间文件：
/
&metgrid
fg_name = 'FILE','SEAICE'#读取以FILE和SEAICE开头的文件
/
实例-海冰场添加
数据信息
下面简单给出一个实例，我们想要给WRF添加海冰场，数据为nc格式，其基本描述为：
dimensions:
x = 304 ;
y = 448 ;
variables:
float sea_ice_thickness(y, x) ;
sea_ice_thickness:units = "Meters" ;
sea_ice_thickness:long_name = "Sea ice thickness" ;
float snow_depth(y, x) ;
snow_depth:units = "Meters" ;
snow_depth:long_name = "Snow depth" ;
float snow_density(y, x) ;
snow_density:units = "Kg / Meters^3" ;
snow_density:long_name = "Snow density" ;
float lat(y, x) ;
lat:units = "Degrees" ;
lat:long_name = "Latitude" ;
float lon(y, x) ;
lon:units = "Degrees" ;
lon:long_name = "Longitude" ;
float freeboard(y, x) ;
freeboard:units = "Meters" ;
freeboard:long_name = "Ice freeboard" ;
float roughness(y, x) ;
roughness:units = "Meters" ;
roughness:long_name = "Ice surface roughness" ;
float ice_con(y, x) ;
ice_con:units = "Percent x 100" ;
ice_con:long_name = "Sea ice concentration" ;
可以看到有7个变量，我们想要提取其中的sea_ice、snow_depeth、ice_coe三个变量写出。
NCL脚本
首先读取数据，并进行简单处理：
begin
nrows = 448
ncols = 304
ymdd1=getenv("20180310")
ymdd2=getenv("2018-03-10")
;print(ymdd1)
data_filename="/public/home/zhangzilu/RDEFT4/RDEFT4_"+ymdd1+".nc"
print(data_filename)
f=addfile(data_filename,"r")
sea_ice_thickness=f->sea_ice_thickness(:,:)
snow_depth=f->snow_depth(:,:)
snow_density=f->snow_density(:,:)
freeboard=f->freeboard(:,:)
roughness=f->roughness(:,:)
ice_con=f->ice_con(:,:)
printVarSummary(sea_ice_thickness)
at=f->lat
lon=f->lon
lat@units="degrees_N"
lon@units="degrees_E"
lat=(/ lat(::-1,:) /)
sea_ice_thickness=(/sea_ice_thickness(::-1,:)/)
snow_depth=(/snow_depth(::-1,:)/)
snow_density=(/snow_density(::-1,:)/)
freeboard=(/snow_depth(::-1,:)/)
roughness=(/freeboard(::-1,:)/)
ice_con=(/ice_con(::-1,:)/)
print(max(lon)+" " +min(lon))
print(max(lat)+" " +min(lat))
sea_ice_thickness@lat2d=lat
sea_ice_thickness@lon2d=lon
snow_depth@lat2d=lat
snow_depth@lon2d=lon
snow_density@lat2d=lat
snow_density@lon2d=lon
freeboard@lat2d=lat
freeboard@lon2d=lon
roughness@lat2d=lat
roughness@lon2d=lon
ice_con@lat2d=lat
ice_con@lon2d=lon
ce_con=where(ice_con.eq.-9999.,-1.E30,ice_con)
sea_ice_thickness=where(sea_ice_thickness.eq.-9999.,-1.E30,sea_ice_thickness)
snow_depth=where(snow_depth.eq.-9999.,-1.E30,snow_depth)
snow_density=where(snow_density.eq.-9999.,-1.E30,snow_density)
freeboard=where(freeboard.eq.-9999.,-1.E30,freeboard)
roughness=where(roughness.eq.-9999.,-1.E30,ice_con)
snow_depth=where(snow_depth.eq.-999.0,-1.E30,snow_depth)
snow_depth=where(snow_depth.eq.-4.995,-1.E30,snow_depth)
roughness=where(roughness.eq.-999.0,-1.E30,roughness)
freeboard=where(freeboard.lt.-250.0,-1.E30,freeboard)
sea_ice_thickness@_FillValue=-1.E30
snow_depth@_FillValue=-1.E30
snow_density@_FillValue=-1.E30
freeboard@_FillValue=-1.E30
roughness@_FillValue=-1.E30
ice_con@_FillValue=-1.E30
end
将读取的数据输出为WPS intermediate file:
undef ("write_to_intermediate")
function write_to_intermediate(vname:string,ddate:string,vvv:numeric,uunits:string,ddesc:string,mmap:string
;定义函数，输入变量：varname(字符串)、ddate（字符串）、变量值、单位、变量描述、数据来源
local DATE,pnew,DATE1,output_file_name,hh,j ;定义函数内变量
begin
hh=(/"00","12"/)
; 定义输入变量与输出文件名
do j=0,1
DATE = ddate+"_"+hh(j)
DATE1=DATE+":00:00"
WPS_IM_root_name = vname
output_file_name = WPS_IM_root_name + ":" + DATE
earth_radius = 6371.23
system("/rm " + output_file_name)
FIELD_ICE = vname
UNITS_ICE = uunits
DESC_ICE = ddesc
pnew =200100
opt = True
opt@projection = 0 ;Cylindrical Equidistant (Lat/lon) projection
opt@date = DATE1
opt@level = 200100
opt@map_source = mmap
opt@startloc = "SWCORNER" ; 8 chars exact
opt@startlon = -180.
opt@startlat = 90.
opt@deltalat=-0.25
opt@deltalon=0.25
; opt@deltalat=-1.0
; opt@deltalon=1.0
opt@is_wind_earth_relative = False
wrf_wps_write_int(WPS_IM_root_name,FIELD_ICE,UNITS_ICE,DESC_ICE,vvv(:,:),opt); 输出
delete(opt)
end do
return(output_file_name)
end
bb=write_to_intermediate("SEAICE",ymdd2,ice_con2,"fraction","sea ice concentration", "RDEFT4 SEA ICE concentration")
cc=write_to_intermediate("ICEDEPTH",ymdd2,sea_ice_thickness1,"M","ice depth", "RDEFT4 25 km SEA ICE THICKNESS")
dd=write_to_intermediate("SNOWSI",ymdd2,snow_depth1,"M","snow depth", "RDEFT4 25 km SNOW DEPTH")
输入metgrid
如下添加：
导出文件如下：
修改namelist.wps，运行metgrid.exe，成功读取后，met_em文件会多出FLAG_SNOWSI=1等字样。
以上，完成。
2124
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
