GRAPES MESO模式学习笔记07——脚本运行之产品生成

meteorology
目录

产品包括但不限于各种不同格式的数据文件（GRIB2格式、MICAPS格式）、图形、推送产品等。

当前GRAPES模式积分结果为GRADS格式，需要转换成GRIB2格式、MICAPS等格式。我的同事们开发了一整套格式转换和插值的工具软件，用在我们的实时业务运行中。
GRAPES模式中的绘图脚本使用GRADS，但我们还有一套独立的作业系统，使用NCL来绘图。
推送产品目前使用共享目录（IBM的GPFS文件系统）和ftp传输方式。
下面先介绍下将GRADS数据转换成GRIB2的脚本。

GRADS to GRIB2

数据转码工具包使用datapost，我用的是GRAPESV3.3的业务系统中用的2.5版本，从大型机的某个目录中拷贝而来。最新版本为2.7，在原有基础上增加集合预报产品的压码功能。我目前使用2.5版本，该软件包提供三个命令
grads2grib，将grads转为grib1或grib2格式
grads2micaps，将grads转为micaps格式
gribpost，将grib1、grib2或micaps格式转换为grads格式
前两个为压码，第三个为解码。
将GRADS转换为GRIB2需要使用grads2grib.exe命令。

输入

namelist.input    参数文件
grib_para_file    grib参数文件
cma.grib2    grib2的模板文件
grads数据文件以及对应的ctl文件

输入文件说明

datapost中有详细的说明文档，可以参阅。我简要说明下

参数文件namelist.input

软件包的示例：
[code]
&input_file
grads_ctl = “./h500.ctl”
/
&model_parameter
model_name = “TEST”
init_time = “2013032200”
start_time = “2013032200”
end_time = “2013032200”
/
&grads_parameter
grads_map_type = 0 ! 0:lon-lat map(GRAPES,AREM) 1:lambert map(MM5,WRF)
nx_grads = 502 ! east-west direction
ny_grads = 330 ! south-north direction
nz_grads = 1
dx_grads = 0.15
dy_grads = 0.15
slon_grads = 70.0
elon_grads = 145.15
slat_grads = 15.0
elat_grads = 64.35
!
corner_lon_grads = 120.0
corner_lat_grads = 42.0
corner_i_grads = 40
corner_j_grads = 35
STD_LON_grads = 120.0 !the same as corner_lon
TRUE_LAT1_grads = 30.0,
TRUE_LAT2_grads = 60.0
MAP_DX_grads = 15000.
MAP_DY_grads = 15000.
/
&grib_parameter1
grib_map_type = 0
nx_grib = 502
ny_grib = 330
dx_grib = 0.15
dy_grib = -0.15
slon_grib = 70.0
elon_grib = 145.15
slat_grib = 64.35
elat_grib = 15.0
!
corner_lon_grib = 120.0
corner_lat_grib = 42.0
corner_i_grib = 40
corner_j_grib = 35
STD_LON_grib = 120.0 !the same as corner_lon
TRUE_LAT1_grib = 30.0,
TRUE_LAT2_grib = 60.0
MAP_DX_grib = 15000.
MAP_DY_grib = 15000.
/
&grib_parameter2
grib_type = 2 ! 1: grib1, 2: grib2, 12: grib1 and grib2
grib_para_file = “./grib_para_file”
dt_grib = 9999 ! 9999: only 1 time level, others: normal value
file_status = 0 ! 0: create a new file, 1: append an old file
quality_control = 0 ! 0: NO, 1: YES
/
[/code]
参数文件中主要设置时间、网格参数等。
目前单个GRIB2文件只有一个时间的输出，所以起始时间和结束时间一致，为起报时间加上预报时效。CTL文件内容只有起始时间，业务中ctl文件名会有起报时间和预报时效。
网格参数中corner系列没看懂。经纬度网格下，其余参数与CTL中xdef、ydef对应。

GRIB参数文件

GRIB参数文件由多行文本组成，每行为一个单元。形如
[code]
33 100 1000 -1000 1000 0 1 0 U|10hpa东西风
33 100 2000 -1000 1000 0 1 0 U|20hpa东西风
33 100 3000 -1000 1000 0 1 0 U|30hpa东西风
33 100 5000 -1000 1000 0 1 0 U|50hpa东西风
33 100 7000 -1000 1000 0 1 0 U|70hpa东西风
33 100 10000 -1000 1000 0 1 0 U|100hpa东西风
…
[/code]
每行对应输出的GRIB2文件中的一个消息，使用grib_ls查看输出的grib2文件：
[code]
rmf.gra.2014061012000.grb2
edition centre date dataType gridType step levelType lev short_name packingType valuesCount
2 babj 20140610 an regular_ll 0 pl 1000 u grid_jpeg 165660
2 babj 20140610 an regular_ll 0 pl 2000 u grid_jpeg 165660
2 babj 20140610 an regular_ll 0 pl 3000 u grid_jpeg 165660
2 babj 20140610 an regular_ll 0 pl 5000 u grid_jpeg 165660
2 babj 20140610 an regular_ll 0 pl 7000 u grid_jpeg 165660
2 babj 20140610 an regular_ll 0 pl 10000 u grid_jpeg 165660
…
[/code]
可以看到每个grib2消息对应GRIB参数文件中的一行，且次序相同。
引用datapost的文档说明GRIB参数文件的每行记录含义：

每行记录： elecode,levType,lev,min,max,bas1,exp,bas2,var_name|description
elecode: 该要素在grib文件中的要素代码,见codeTable.txt
levType: 该要素在grib文件中的类型代码,等压面:100 地面:1 海平面:2 10米风或2米温度:105
lev: 层次,若levType值为100，则lev单位为‘Pa’，若levType的值为1或者2，则lev的值为0，若levType的值为105，则lev为10或者2
min: 最小值，用于质量控制
max: 最大值，用于质量控制
bas1: 变换系数之一： x_new=(x_org+bas1)*exp+bas2
exp: 变换系数之一： x_new=(x_org+bas1)*exp+bas2
bas2: 变换系数之一： x_new=(x_org+bas1)*exp+bas2
var_name: 变量名，大写，要与CTL文件中的变量名一致
describe: 对变量的描述

要素代码说明表格示例如下：

<td width="72">
  discipline
</td>

<td width="72">
  Category
</td>

<td width="72">
  number
</td>

<td width="104">
  物理意义
</td>

<td width="135">
  单位
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
</td>

<td width="104">
  地面气压
</td>

<td>
  pa
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
  1
</td>

<td width="104">
  海平面气压
</td>

<td>
  pa
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
  5
</td>

<td width="104">
  位势高度（包括地形高度）
</td>

<td>
  gpm
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
  6
</td>

<td width="104">
  几何高度
</td>

<td>
  M
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="104">
  温度（包括地面温度和2米温度）
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  2
</td>

<td width="104">
  位温
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="104">
  假相当位温
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  4
</td>

<td width="104">
  最高温度（包括2米最高温度）
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  5
</td>

<td width="104">
  最低温度（包括2米最低温度）
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  6
</td>

<td width="104">
  露点温度
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  7
</td>

<td width="104">
  温度露点差
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
  19
</td>

<td width="72">
</td>

<td width="104">
  能见度
</td>

<td>
  m
</td>

<td width="72">
</td>

<td width="72">
  7
</td>

<td width="72">
</td>

<td width="104">
  气块抬升指数
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
</td>

<td width="72">
  9
</td>

<td width="104">
  24<span>小时变温</span>
</td>

<td>
  K
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
  8
</td>

<td width="104">
  24<span>小时变压</span>
</td>

<td>
  pa
</td>

<td width="72">
</td>

<td width="72">
  3
</td>

<td width="72">
  9
</td>

<td width="104">
  24<span>小时变高</span>
</td>

<td>
  gpm
</td>

<td width="72">
</td>

<td width="72">
  2
</td>

<td width="72">
  1
</td>

<td width="104">
  全风速
</td>

<td>
  m/s
</td>

<td width="72">
</td>

<td width="72">
  2
</td>

<td width="72">
  2
</td>

<td width="104">
  东西风(包括10米U）
</td>

<td>
  m/s
</td>

<td width="72">
</td>

<td width="72">
  2
</td>

<td width="72">
  3
</td>

<td width="104">
  南北风(包括10米V）
</td>

<td>
  m/s
</td>

不清楚如何计算得到的要素代码，是否由学科、分类、编号的三个数字计算得到？
层次类型，应该与GRIB2的层次类型代码相同，参见GRIB 2 CODE TABLE 4.5。
变量名是从CTL中获取要素在数据文件中位置的重要参数。不清楚如何获取变量，待以后阅读代码时分析。

GRIB2模板文件

cma.grib2的grib_dump -O输出如下（略去最后两段）：
[code collapse=”true”]
***** FILE: cma.grib2
#============== MESSAGE 1 ( length=4000174 ) ==============
1-4 identifier = GRIB
5-6 reserved = MISSING
7 discipline = 0 [Meteorological products (grib2/0.0.table) ]
8 editionNumber = 2
9-16 totalLength = 4000174
====================== SECTION_1 ( length=21, padding=0 ) ======================
1-4 section1Length = 21
5 numberOfSection = 1
6-7 identificationOfOriginatingGeneratingCentre = 98 [European Center for Medium-Range Weather Forecasts (grib1/0.table) ]
8-9 identificationOfOriginatingGeneratingSubCentre = 0
10 gribMasterTablesVersionNumber = 4 [Unknown code table entry (grib2/1.0.table) ]
11 versionNumberOfGribLocalTables = 0 [Local tables not used (grib2/1.1.table) ]
12 significanceOfReferenceTime = 1 [Start of forecast (grib2/1.2.table) ]
13-14 year = 2008
15 month = 2
16 day = 6
17 hour = 12
18 minute = 0
19 second = 0
20 productionStatusOfProcessedData = 0 [Operational products (grib2/1.3.table) ]
21 typeOfProcessedData = 2 [Analysis and forecast products (grib2/1.4.table) ]
====================== SECTION_3 ( length=72, padding=0 ) ======================
1-4 section3Length = 72
5 numberOfSection = 3
6 sourceOfGridDefinition = 0 [Specified in Code table 3.1 (grib2/3.0.table) ]
7-10 numberOfDataPoints = 496
11 numberOfOctetsForOptionalListOfNumbersDefiningNumberOfPoints = 0
12 interpretationOfListOfNumbersDefiningNumberOfPoints = 0 [There is no appended list (grib2/3.11.table) ]
13-14 gridDefinitionTemplateNumber = 0 [Latitude/longitude. Also called equidistant cylindrical, or Plate Carree (grib2/3.1.table) ]
15 shapeOfTheEarth = 0 [Earth assumed spherical with radius = 6,367,470.0 m (grib2/3.2.table) ]
16 scaleFactorOfRadiusOfSphericalEarth = MISSING
17-20 scaledValueOfRadiusOfSphericalEarth = MISSING
21 scaleFactorOfMajorAxisOfOblateSpheroidEarth = MISSING
22-25 scaledValueOfMajorAxisOfOblateSpheroidEarth = MISSING
26 scaleFactorOfMinorAxisOfOblateSpheroidEarth = MISSING
27-30 scaledValueOfMinorAxisOfOblateSpheroidEarth = MISSING
31-34 numberOfPointsAlongAParallel = 16
35-38 numberOfPointsAlongAMeridian = 31
39-42 basicAngleOfTheInitialProductionDomain = 0
43-46 subdivisionsOfBasicAngle = 0
47-50 latitudeOfFirstGridPoint = 60000000
51-54 longitudeOfFirstGridPoint = 0
55 resolutionAndComponentFlags = 48 [00110000]
56-59 latitudeOfLastGridPoint = 0
60-63 longitudeOfLastGridPoint = 30000000
64-67 iDirectionIncrement = 2000000
68-71 jDirectionIncrement = 2000000
72 scanningMode = 0 [00000000]
====================== SECTION_4 ( length=34, padding=0 ) ======================
1-4 section4Length = 34
5 numberOfSection = 4
6-7 numberOfVerticalCoordinateValues = 0
8-9 productDefinitionTemplateNumber = 0 [Analysis or forecast at a horizontal level or in a horizontal layer at a point in time (grib2/4.0.table) ]
10 parameterCategory = 0 [Temperature (grib2/4.1.0.table) ]
11 parameterNumber = 0 [Temperature (K) (grib2/4.2.0.0.table) ]
12 typeOfGeneratingProcess = 0 [Analysis (grib2/4.3.table) ]
13 backgroundGeneratingProcessIdentifier = 255
14 generatingProcessIdentifier = 130
15-16 hoursAfterReferenceTimeOfDataCutoff = 0
17 minutesAfterReferenceTimeOfDataCutoff = 0
18 indicatorOfUnitOfTimeRange = 1 [Hour (grib2/4.4.table) ]
19-22 forecastTime = 0
23 typeOfFirstFixedSurface = 103 [Specified height level above ground (m) (grib2/4.5.table) ]
24 scaleFactorOfFirstFixedSurface = MISSING
25-28 scaledValueOfFirstFixedSurface = MISSING
29 typeOfSecondFixedSurface = 255 [Missing (grib2/4.5.table) ]
30 scaleFactorOfSecondFixedSurface = MISSING
31-34 scaledValueOfSecondFixedSurface = MISSING
====================== SECTION_5 ( length=21, padding=0 ) ======================
1-4 section5Length = 21
5 numberOfSection = 5
6-9 numberOfValues = 496
10-11 dataRepresentationTemplateNumber = 0 [Grid point data – simple packing (grib2/5.0.table) ]
12-15 referenceValue = 270.467
16-17 binaryScaleFactor = -10
18-19 decimalScaleFactor = 0
20 numberOfBitsContainingEachPackedValue = 16
21 typeOfOriginalFieldValues = 0 [Floating point (grib2/5.1.table) ]
====================== SECTION_6 ( length=6, padding=0 ) ======================
1-4 section6Length = 6
5 numberOfSection = 6
6 bitMapIndicator = 255 [A bit map does not apply to this product (grib2/6.0.table) ]
[/code]
该文件作为创建GRIB2消息的样例，datapost应该使用grib api提供的从样例消息创建新grib2消息的函数，这样只需修改部分参数满足要求即可。修改的参数可以参看//www.diffchecker.com/mk36x56z，其中左边为某个实际输出的grib2文件，右边为cma.grib样例文件。

输出

grib2文件

脚本

[code]
#!/bin/ksh
# @ comment = GRAPES
# @ job_type = serial
# @ input = /dev/null
# @ output = ./out/grads2grib2_$(jobid).out
# @ error = ./out/grads2grib2_$(jobid).err
# @ initialdir =./
# @ notification = error
# @ checkpoint = no
# @ class = normal
# @ queue
#############################
# grads to grib2 orgin
#############################
set -x
set -u
set -e
date
######################
# directory setting
######################
BASE_DIR=./grapes_meso/my_grapes_meso
BASE_SOURCE_DIR=./grapes_meso/GRAPES_MESO3.3.2.4
BASE_CONDAT_DIR=./grapes_meso/my_grapes_meso/condat
BIN_DIR=${BASE_DIR}/bin
PRODS_CONDAT_DIR=${BASE_CONDAT_DIR}/prods
BASE_RUN_DIR=${BASE_DIR}/run
RUN_DIR=${BASE_RUN_DIR}/prods
POST_RUN_DIR=${BASE_RUN_DIR}/post
AN_RUN_DIR=${BASE_RUN_DIR}/an
SI_RUN_DIR=${BASE_RUN_DIR}/si
FCST_RUN_DIR=${BASE_RUN_DIR}/fcst
GRADS2GRIB_BIN_DIR=${BIN_DIR}/grads2grib.exe
############
# parameter
############
levels=26
YYYYMMDDHH=2014061012
current_forcast_hour=000
while getopts “t:b:r:g:c:” arg
do
case $arg in
t)
YYYYMMDDHH=$OPTARG
;;
b)
BCKG_DIR=$OPTARG
;;
r)
RUN_DIR=$OPTARG
;;
esac
done
model_name=GRAPES_MESO2GRIB2_ORIG
begintime=$YYYYMMDDHH
init_time=${begintime}
start_time=$( smsdate ${begintime} +${current_forcast_hour})
end_time=$start_time
cur_time_str=${begintime}${current_forcast_hour}00
#####################
# enter work space
#####################
work_space_dir=${RUN_DIR}/grib2_orig/${current_forcast_hour}
test -d ${work_space_dir} || mkdir -p ${work_space_dir}
cd ${work_space_dir}
#####################
# prepare data
#####################
ln -sf ${FCST_RUN_DIR}/postvar${cur_time_str} .
ln -sf ${FCST_RUN_DIR}/post.ctl_${cur_time_str} .
cp ${PRODS_CONDAT_DIR}/grib/cma.grib2 .
cp ${PRODS_CONDAT_DIR}/grib/grib_para_file.base.orig.grib2 grib_para_file
###########
# namelist.input
###########
ctl_file_name=post.ctl_${cur_time_str}
cat > namelist.input «eof &input_file grads_ctl = “./${ctl_file_name}” / &model_parameter model_name = “${model_name}” init_time = “${init_time}” start_time = “${start_time}” end_time = “${end_time}” / &grads_parameter grads_map_type = 0 ! 0:lon-lat map(GRAPES,AREM) 1:lambert map(MM5,WRF) nx_grads = 502 ! east-west direction ny_grads = 330 ! south-north direction nz_grads = 1 dx_grads = 0.15 dy_grads = 0.15 slon_grads = 70.0 elon_grads = 145.15 slat_grads = 15.0 elat_grads = 64.35 ! corner_lon_grads = 120.0 corner_lat_grads = 42.0 corner_i_grads = 40 corner_j_grads = 35 STD_LON_grads = 120.0 !the same as corner_lon TRUE_LAT1_grads = 30.0, TRUE_LAT2_grads = 60.0 MAP_DX_grads = 15000. MAP_DY_grads = 15000. / &grib_parameter1 grib_map_type = 0 nx_grib = 502 ny_grib = 330 dx_grib = 0.15 dy_grib = -0.15 slon_grib = 70.0 elon_grib = 145.15 slat_grib = 64.35 elat_grib = 15.0 ! corner_lon_grib = 120.0 corner_lat_grib = 42.0 corner_i_grib = 40 corner_j_grib = 35 STD_LON_grib = 120.0 !the same as corner_lon TRUE_LAT1_grib = 30.0, TRUE_LAT2_grib = 60.0 MAP_DX_grib = 15000. MAP_DY_grib = 15000. / &grib_parameter2 grib_type = 2 ! 1: grib1, 2: grib2, 12: grib1 and grib2 grib_para_file = “./grib_para_file” dt_grib = 9999 ! 9999: only 1 time level, others: normal value file_status = 0 ! 0: create a new file, 1: append an old file quality_control = 0 ! 0: NO, 1: YES / EOF ##################### # run program ##################### ${GRADS2GRIB_BIN_DIR} ###################### # END ###################### date [/code]