问题简介
我们都知道,在运行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等字样。
以上,完成。