WRF学习——使用CMIP6数据驱动WRF/基于ncl与vdo的CMIP6数据处理

动力降尺度

国际耦合模式比较计划（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脚本实现。

CMIP6数据准备

本次使用的模式为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

#suffix="_6hrPlevPt_MPI-ESM1-2-HR_ssp126_r1i1p1f1_gn_201501010600-202001010000.nc"
startdate=$1
enddate=$2
#date1=$1
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")

#echo ${date1}

do
date1=`date -d "${startdate}" +%Y-%m-%d`
echo ${date1}
  for var in ${var2d[@]}
  do
  echo ${var}
   for hour in ${hlist[@]}
   do
   echo ${hour}
   out_filename=${var}_${date1}_${hour}.nc
   #input_filename=${var}${suffix}
   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"
   done
 done
startdate=`date -d "+1 day ${startdate}" +%Y%m%d`
done

在运行时，输入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")

#date1=$1
startdate=$1
enddate=$2
hlist=("00" "06" "12" "18")
while [[ ${startdate} -lt ${enddate} ]]
do
for hour in ${hlist[@]}
do
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}
#cdo merge tas_${suffix} uas_${suffix} vas_${suffix} ts_${suffix} huss_${suffix} tsl_${suffix} psl_${suffix} snw_${suffix} 2D_${outputfile}
done
startdate=`date -d "+1 day ${startdate}" +%Y%m%d`
echo ${startdate}
done

运行后可得到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)
end

在这里定义了一个regird_mpi函数，在使用是输入文件名以及要插值的变量即可，根据经纬度网格点不同，可将该函数中的curvilinear_to_SCRIP和rectilinear_to_SCRIP相互替换。

WRF中间文件撰写

最后，要将插值后的变量数据，转写为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