WRF学习笔记之四:撰写WPS intermediate file添加海冰场/NCL学习/WRF进阶:如何向WRF添加额外气象场?

问题简介

我们都知道,在运行WRF模式之前,需要使用其WPS前处理程序,制作对应的模拟域与气象场输入,而WPS主要由3个主程序组成,分别为:

  1. geogrid.exe定义模式的模拟域,并将静态地理数据插值到模式网格中。
  2. ungrib.exe 从Grib的格式的文件读取气象数据。
  3. metgrid.exe将读取的气象数据,水平插值到geogrid.exe定义的模拟域网格中。

通过以上介绍不难看出,整个WPS的处理流程并不复杂,而最后我们的目的就是要将模拟域与气象场结合,输入到WRF中。
在运行WPS过程中,metgrid.exe需要最后运行,因为它是将其他两个主程序处理后的数据进行插值结合的程序,即:从geogrid.exe中获得模拟域与地理信息,从ungrib.exe中获得气象数据,而metgrid.exe从两个程序中获取信息的方式则为:读取geogrid.exeungrib.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等字样。
在这里插入图片描述
以上,完成。

  • 5
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 4
    评论
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值