本文相关视频讲解:
CMIP数据驱动WRF
国际耦合模式比较计划(CMIP)为研究不同情景下的气候变化提供了大量的模拟数据,而在实际研究中,全球气候模式输出的数据空间分辨率往往较低(>100Km,缺乏区域气候特征,为了更好地研究不同情景下,某一区域的气候变化特征,我们往往需要更高分辨率的模拟数据。此时便需要对全球模式输出的数据进行降尺度研究,将大尺度信息变量与小尺度信息变量建立联系,获得更多的小尺度的变量。通常,降尺度可分为(1)动力降尺度 (2)统计降尺度 (3)二者结合降尺度
动力降尺度通常是将全球模式输出的数据作为驱动场,输入至区域气候模式中,从而获取描述区域气候特征的更高分辨率数据。WRF作为最常用的中尺度天气预测模式,将其与最新的CMIP6全球模式数据结合,是动力降尺度最常用的方法。
下面以CMIP6数据中的MPI-ESM2-HR数据为例,展示使用CMIP6数据驱动WRF的基本步骤。
对于CMIP6驱动WRF,已有老师在github上上传了基于python的工具包,习惯在LINUX下使用python的用户可以尝试:
cmip6-to-wrfnterm
基本思路都是类似的,只不过本文主要是基于服务器现有的NCL、CDO与shell脚本实现。
本次使用的模式为MPI-ESM2-HR数据,空间分辨率为100km,选择原因:数据全面、分辨率好、应用较多,可自己根据需求去官网下载数据。
本次所需下载并驱动的变量有:
v_name
|
wrf_name
|
units
|
dim
|
desc
|
notes
|
ps
|
PSFC
|
Pa
|
2d
|
surface pressure
|
|
psl
|
PMSL
|
Pa
|
2d
|
Mean sea-level pressure
|
|
zg
|
GHT
|
m
|
3d
|
geopotential height
|
|
ta
|
TT
|
K
|
3d
|
air temperature
|
|
tas
|
TT
|
K
|
2d
|
2-m temerature
|
|
ua
|
UU
|
m s-1
|
3d
|
u-component wind;
|
|
uas
|
UU
|
m s-1
|
2d
|
`10-m u-component wind
|
|
va
|
VV
|
m s-1
|
3d
|
v-component wind
|
|
vas
|
VV
|
m s-1
|
2d
|
10-m v-component wind
|
|
hus
|
SPECHUMD
|
kg kg-1
|
3d
|
specific humidity
|
|
huss
|
SPECHUMD
|
kg kg-1
|
2d
|
10-m specific humidity
|
|
ts
|
SKINTEMP
|
K
|
2d
|
Skin temperature
|
|
tsl
|
ST000010
|
K
|
2d
|
0-10cm soil temperature
|
|
tos
|
SST
|
K
|
2d
|
Sea temperature
|
optional
|
mrsos
|
SM000010
|
m3 m-3
|
2d
|
0-10cm soil moisture
|
|
snw
|
SNOW
|
kg m-2
|
2d
|
snow mass
|
optional
|
sic
|
SEAICE
|
1
|
2d
|
seaice
|
optional
|
下载数据命名一般为以下格式:
变量名称_时间分辨率_层级_模式名称_情景名称_年份时间,在MPI-ECSM-HR变量的对应需要去自行查询变量表格、
首先由于下载的数据通常是10年一个,致使在实际使用时,我们首先要将需要对应时段的变量提取出来,同样的,为了后续使用shell脚本更好地进行批量处理,这些变量应当以一种特定的格式命名。
使用cdo 的selvar seldate功能,并结合shell脚本,便能很好的完成:
#!/bin/bash
startdate=$1
enddate=$2
hlist=("00" "06" "12" "18")
varlist=("ta" "ua" "va" "zg" "tas" "uas" "vas" "ts" "tsl" "snw" "hus" "huss" "psl")
var2d=("tas" "uas" "vas" "ts" "tsl" "snw" "huss" "psl")
var3d=("ta" "ua" "va" "zg" "hus")
date1=`date -d "${startdate}" +%Y-%m-%d`
echo ${date1}
for var in ${var2d[@]}
echo ${var}
for hour in ${hlist[@]}
echo ${hour}
out_filename=${var}_${date1}_${hour}.nc
input_filename=`ls ${var}_*ssp*`
echo ${out_filename}
echo ${input_filename}
cdo -seldate,${date1} -selhour,${hour} -selname,${var} ${input_filename} ${out_filename}
echo "cdo done"
startdate=`date -d "+1 day ${startdate}" +%Y%m%d`
在运行时,输入bash run.sh startdate enddate便可将相应时间段提取并输出为变量_时间的形式,如:
注意:下载的CMIP6数据变量存在2D与3D区别,在处理3D变量,如ua va时,cdo还应当加上sellevel-提取对应的层数,这是由于该数据中ua va的垂直层仅有7层,而ta hus zg则有28层,在后续运行WRF时,3D数据的层次应当保持一致!!!
处理好后的数据,为了方便,根据2d与3d的不同将其合并:
#!/bin/bash
var3d=("ta" "ua" "va" "zg" "hus")
var2d=("tas" "uas" "vas" "ts" "tsl" "snw" "huss" "psl")
startdate=$1
enddate=$2
hlist=("00" "06" "12" "18")
while [[ ${startdate} -lt ${enddate} ]]
for hour in ${hlist[@]}
date1=`date -d "${startdate}" +%Y-%m-%d`
echo ${date1}
echo ${hour}
suffix=${date1}_${hour}.nc
echo ${suffix}
outputfile=MPI_HR_${date1}_${hour}_00:00.nc
echo ${outputfile}
cdo merge ta_${suffix} ua_${suffix} va_${suffix} zg_${suffix} hus_${suffix} 3D_${outputfile}
startdate=`date -d "+1 day ${startdate}" +%Y%m%d`
echo ${startdate}
运行后可得到2d_MPI_HR和3dMPI_HR文件。
我们应当注意的是,CMIP6的经纬度很多时候并不是等经纬度间距的,比如我下载的数据就是100km,在海洋上分辨率有时达到50km。这就使得我们在运行之前,首先要将其插值到均一的lat/lon坐标下,否则WRF将很难处理。
值得注意的是,CMIP6数据可分为大气与海洋两部分,而大气的经纬度与海洋的经纬度则存在差异,比如,大气的经纬度数据为一维数据,海洋则以二维数据给出,因此插值时需要分开处理。请在插值前弄清楚变量的经纬度网格。
tas经纬度,以一维表征
tos经纬度网格为曲线网格,经纬度为二维形式
对于两种在ncl中使用不同的插值函数即可,对一维使用rectilinear_to_SCRIP
将经纬度转为映射文件,在使用ESMF_regrid_with_weights
插值,对二维曲线网格,使用curvilinear_to_SCRIP
函数,再使用ESMF_regrid_with_weights
插值。
以下为插值的代码示例:
undef ("regrid_MPI")
function regrid_MPI(fname:string,inputv:numeric)
local regrid_var,MPI_var,lat,lon,inputf
begin
inputf=addfile(fname,"r")
lat=inputf->latitude
lon=inputf->longitude
Opt = True
Opt@SrcRegional = True
Opt@ForceOverwrite = True
Opt@PrintTimings = True
Opt@Title = "MPI-ESM1"
Opt@CopyVarAtts = True
;Opt@GridMask = where(.not.ismissing(zg),1,0)
Opt@CopyVarCoords = False
srcGridName = "SCRIP_MPI-ESM1_grid"+".nc"
curvilinear_to_SCRIP(srcGridName, lat,lon, Opt)
delete(Opt)
;----------------------------------------------------------------------
; Convert destination grid to a SCRIP convention file.
;----------------------------------------------------------------------
dstGridName = "dst_SCRIP.nc"
Opt = True
Opt@LLCorner = (/ -90.d, 0.d/)
Opt@URCorner = (/ 90.d,360.d/)
Opt@ForceOverwrite = True
Opt@PrintTimings = True
latlon_to_SCRIP(dstGridName,"1x1",Opt)
;---Clean up
delete(Opt)
;----------------------------------------------------------------------
; Generate the weights that take you from the NCEP grid to a
; 1x1 degree grid.
;----------------------------------------------------------------------
wgtFileName = "MPI_2_Rect.nc"
Opt = True
Opt@InterpMethod = "bilinear" ; default
Opt@ForceOverwrite = True
Opt@PrintTimings = True
ESMF_regrid_gen_weights(srcGridName,dstGridName,wgtFileName,Opt)
delete(Opt)
;---------------------------------
;----------------------------------------------------------------------
; Apply the weights to a given variable
;----------------------------------------------------------------------
Opt = True
Opt@PrintTimings = True
;---In V6.1.0, coordinates and attributes are copied automatically
regrid_var = ESMF_regrid_with_weights(inputv,wgtFileName,Opt)
;printVarSummary(regrid_var)
return(regrid_var)
在这里定义了一个regird_mpi函数,在使用是输入文件名以及要插值的变量即可,根据经纬度网格点不同,可将该函数中的curvilinear_to_SCRIP
和rectilinear_to_SCRIP
相互替换。
最后,要将插值后的变量数据,转写为WPS的中间文件,该中间文件可直接被metgrid.exe读取,并生成met_em文件。
ncl中就有现成的函数,需要注意的是,撰写时变量的FIELD应当包括在METGRID.TBL中相同,否则metgrid无法识别。
如果数据是2d,则直接撰写,如果数据为3d,则根据层数,循环一层层写:
; for 2d variable
FIELD_ICE = "SEAICE"
UNITS_ICE = "1"
DESC_ICE = "ocean seaice"
FIELD_ST = "SST"
UNITS_ST = "K"
DESC_ST = "sea surface temperature"
re_sic=regrid_MPI(data_filename,sic)
; re_tos=regrid_MPI(data_filename,tos)
opt = True
opt@projection = 0 ; "Equidistant_Lat_Lon"
opt@date = DATE1
opt@map_source = "1×1"
opt@startloc = "SWCORNER" ; 8 chars exact
opt@startlon = 0
opt@startlat = -90
opt@deltalon = 1
opt@deltalat = 1
;opt@is_wind_earth_rel = False
opt@is_wind_earth_relative = False
opt@level = 200100
wrf_wps_write_int(WPS_IM_root_name,FIELD_ICE,UNITS_ICE,\
DESC_ICE,re_sic,opt)
pnew2=(/925,850,700,600,500,250,50/)*100
; For 3D variables
do jlev=0,NLEV2-1
opt@level = pnew2(jlev)
wrf_wps_write_int(WPS_IM_root_name,FIELD_U,UNITS_U,\
DESC_U,UonP(jlev,:,:),opt)
wrf_wps_write_int(WPS_IM_root_name,FIELD_V,UNITS_V,\
DESC_V,VonP(jlev,:,:),opt)
end do
最后使用WPS文件夹下的util/./rd_intermediate.exe 确认是否读取成功,关于这一部分,可参考我以前的博客:撰写WPS intermediate file添加海冰场
将输出的中间文件链接至WPS文件夹,修改namelist.wps文件夹&metgrid部分的fg_name,使其读取我们撰写的中间文件。
之后的步骤就和普通运行WRF一样了,请注意由于数据本身的性质,土壤层数与垂直层数量较少,在设置namelist时记得修改与其保持一致。
主要的坑在于数据本身。
- 3D数据垂直层不一致,ua与va数据仅有7层,而ta hus zg等数据却又有8层,因此在撰写文件时必须注意对应层数保持一致。
- 经纬度格点不一致,注意经纬度各点的类别,在插值时注意区分。
- CMIP6数据本身并不是再分析资料,而是预测气候变化的模型输出,常常会出现数据量与WRF所需不对应的问题,请注意各个变量描述。
- 除了2D与2D以外,也可以使用CMIP6的2D静态数据,如LANDSEA,步骤类似,主要是通过namelist.wps中的constant_name来设置读取。
相关代码与数据已经发布在Github上,请查询:Write_CMIP_to_wps_int
偏差校正数据具有基于 ERA5 的平均气候和年际方差,但具有来自 18 个 CMIP6 模型的集合平均值的非线性趋势。该数据集涵盖 1979-2014 年的历史时间段以及 2015-2100 年的未来情景(SSP245 和 SSP585),水平网格间距为(1.25° × 1.25°),每隔六小时一次。我们的评估表明,在气候平均值、年际方差和极端事件方面,偏差校正数据的质量优于单独的 CMIP6 模型。该数据集将有助于对地球未来气候、大气环境、水文、农业、风力发电等进行动态降尺度预测。
全球气候模型(Global Climate Model, GCM),亦称全球环流模型或全球大气模型,是一种数值模型,被广泛用于模拟地球的气候系统。GCM利用一系列的数学公式来描绘气候系统的各个主要组成部分,包括大气、海洋、冻土以及地表和海洋表面的生物地理过程。GCM的空间和时间精度可以根据需要进行调整。这些模型为我们提供了理解气候系统运行机制的途径,为预测气候变化趋势、评估气候变化对人类社会和生态系统的影响以及制定应对气候变化的策略提供了关键工具。
《第六次国际耦合模式比较计划(CMIP6)评述》周天军
链接: http://www.climatechange.cn/CN/10.12006/j.issn.1673-1719.2019.193
这篇文章介绍了CMIP的发展历程,介绍了CMIP6的组织思路、其核心实验和23个比较子计划,读这篇文章,可以对CMIP6有个整体概念。文章部分截图如下:
CMIP6数据被广泛应用于全球和地区的气候变化研究、极端天气和气候事件研究、气候变化影响和风险评估、气候变化的不确定性研究、气候反馈和敏感性研究以及气候政策和决策支持等多个领域。这些数据为我们理解和预测气候变化,评估气候变化的影响和风险,以及制定有效的气候政策和决策提供了关键的信息和工具。
注意:CMIP6数据网站属国外网站,比较卡,需耐心等待及其下载。
这里对一般情况下,我们用到的检索条件做个简要的介绍,Variable下选择我们需要的变量,如降雨、降雪、均温等等,Frequency下选择我们需要的尺度,如日尺度、月尺度,小时尺度等等,Experiment ID 下选择我们需要的实验,如历史模拟(historical),各种ssp情景等等,一般来说,选了这三个检索条件后,就先点击search,在结果中再点击Source ID,查看有哪些模型符合这个检索
因此,对于某一具体的未来时段,可以通过计算过去和现在气候的差值(即 delta),并将其应用到未来的气候预测上,来预估未来的气候状态。,也被称为全球环流模型或全球大气模型,是一种用于模拟地球的气候系统的数值模型。温度指数:这些指数主要用于度量温度的极端情况,例如热日数(TX90p,年中最高气温超过90百分位数的天数)、冷日数(TN10p,年中最低气温低于10百分位数的天数)、热夜数(TN90p,年中最低气温超过90百分位数的天数)、冷夜数(TN10p,年中最低气温低于10百分位数的天数)等。