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[@]}doecho ${var}for hour in ${hlist[@]}doecho ${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"donedone
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/bashvar3d=("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->longitudeOpt                = 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   = Truelatlon_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                      = TrueOpt@InterpMethod         = "bilinear"     ; defaultOpt@ForceOverwrite       = TrueOpt@PrintTimings         = TrueESMF_regrid_gen_weights(srcGridName,dstGridName,wgtFileName,Opt)delete(Opt);---------------------------------
;----------------------------------------------------------------------
; Apply the weights to a given variable
;----------------------------------------------------------------------Opt                = TrueOpt@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 variableFIELD_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             = 200100wrf_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-1opt@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