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

news/2024/7/8 1:20:15 标签: 学习, 驱动开发

动力降尺度

国际耦合模式比较计划（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_namewrf_nameunitsdimdescnotes
psPSFCPa2dsurface pressure
pslPMSLPa2dMean sea-level pressure
zgGHTm3dgeopotential height
taTTK3dair temperature
tasTTK2d2-m temerature
uaUUm s-13du-component wind;
uasUUm s-12d`10-m u-component wind
vaVVm s-13dv-component wind
vasVVm s-12d10-m v-component wind
husSPECHUMDkg kg-13dspecific humidity
hussSPECHUMDkg kg-12d10-m specific humidity
tsSKINTEMPK2dSkin temperature
tslST000010K2d0-10cm soil temperature
tosSSTK2dSea temperatureoptional
mrsosSM000010m3 m-32d0-10cm soil moisture
snwSNOWkg m-22dsnow massoptional
sicSEAICE12dseaiceoptional

下载数据命名一般为以下格式：
在这里插入图片描述
变量名称_时间分辨率_层级_模式名称_情景名称_年份时间，在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时记得修改与其保持一致。

要点（坑）

主要的坑在于数据本身。

  1. 3D数据垂直层不一致，ua与va数据仅有7层，而ta hus zg等数据却又有8层，因此在撰写文件时必须注意对应层数保持一致。
  2. 经纬度格点不一致，注意经纬度各点的类别，在插值时注意区分。
  3. CMIP6数据本身并不是再分析资料，而是预测气候变化的模型输出，常常会出现数据量与WRF所需不对应的问题，请注意各个变量描述。
  4. 除了2D与2D以外，也可以使用CMIP6的2D静态数据，如LANDSEA，步骤类似，主要是通过namelist.wps中的constant_name来设置读取。

相关代码与数据已经发布在Github上，请查询：Write_CMIP_to_wps_int


http://www.niftyadmin.cn/n/5536022.html

相关文章

c# 的 goto

c# 的 goto

搞循环感觉没什么必要 int number 0; Console.WriteLine("请输入一个数字（输入-1结束）:"); start: // 标签 number int.Parse(Console.ReadLine()); if (number -1) { Console.WriteLine("程序结束。"); } else { Cons…
阅读更多...
数据结构（期末）

数据结构（期末）

目录 逻辑结构 存储结构 算法有以下五个特性 算法数据结构 程序 时间复杂度 空间复杂度 数据元素是数据的基本单位 数据项是数据的最小单位 数据结构是带有结构的各数据元素的集合 时间复杂度例题 线性表 关于带头节点的单链表及不带头节点的单链表 栈和队列 栈…
阅读更多...
阿里模型调用体验

阿里模型调用体验

引言 随着人工智能技术的飞速发展，大型模型已成为推动技术进步的关键因素之一。阿里大模型作为国内领先的人工智能技术之一，其在多个领域的应用展示了强大的潜力。本文将通过调用案例，简单解析阿里大模型在特定场景中的应用及其效果。 1.导…
阅读更多...
如何更改 Python pip 源为国内源

如何更改 Python pip 源为国内源

在使用 Python 安装包工具 pip 时，经常会遇到下载速度慢的问题。这通常是因为默认使用的官方源 https://pypi.org/simple 在国内访问速度较慢。为了提高下载速度，我们可以将 pip 源更改为国内的镜像源。本文将介绍如何临时和永久地更改 pip 源为国内源。…
阅读更多...
【Unix】SunOS/Oracle Solaris系统介绍

【Unix】SunOS/Oracle Solaris系统介绍

一.SunOS系统介绍 SunOS 是由 Sun Microsystems 开发的 Unix 操作系统。它最初是为 Sun 的 SPARC 架构计算机设计的，后来也支持了 Intel x86 架构。SunOS 是基于 UNIX System V 4.1 版本，并且随着时间的发展，SunOS 经历了多个版本迭代&#…
阅读更多...
Linux安装ftp、Java的FTP上传下载文件工具类

Linux安装ftp、Java的FTP上传下载文件工具类

Linux安装ftp、Java的FTP上传下载文件工具类 文章说明Linux安装vsftpdJava的工具类 文章说明 网上找到说linux安装ftp，采用vsftpd，在后续的配置中少了一些说明，给我折磨了许久，写下这篇文章来记录 Linux安装vsftpd 命令非常简单&a…
阅读更多...
数据同步软件有哪些

数据同步软件有哪些

数据同步软件有哪些呢？随着企业规模的扩大，企业数据也积累得越来越多，万一发生宕机风险，那么这个损失将不可估量。所以为了容灾备用，我们往往需要将数据同步到另一台备胎服务器上，进行冗余。 那么需要同步的…
阅读更多...
2024-07-04 base SAS programming学习笔记8(HTML）

2024-07-04 base SAS programming学习笔记8(HTML）

当使用ODS来进行结果或数据集输出的时候，可以同时设置多个ODS 命令，同时输出到多个不同的文件。使用_ALL_ 表示关闭所有的ODS输出窗口，比如： ods html file(body)"html-file-pathname"; ods html file"pdf-file-pa…
阅读更多...
最新文章

Copyright @ 2022~2023