模式结果和观测的对比

首先对模式的海表温度进行对比

(base) [chengxl@login02 10yearmean]$ ls sst_*.nc
sst_2000.nc sst_2002.nc sst_2004.nc sst_2006.nc sst_2008.nc

sst_2001.nc sst_2003.nc sst_2005.nc sst_2007.nc sst_2009.nc

首先将观测资料转化为年平均

 ls sst.200* |xargs -I{} cdo -b f32 yearmean {} sst.2000-2009ym.nc

(base) [chengxl@login02 10yearmean]$ ls sst.2000-2009ym.nc

sst.2000-2009ym.nc

接着对分辨率进行修改

(base) [chengxl@login02 10yearmean]$ ncdump -h sst.2000-2009ym.nc

netcdf sst.2000-2009ym {
dimensions:
time = UNLIMITED ; // (1 currently)
bnds = 2 ;
longitude = 480 ;
latitude = 241 ;

可以看到是241x480的

我们将其转化为和模式一样的分辨率，查看一下模式的分辨率

(base) [chengxl@login02 run]$ ncdump -h MMEAN2000-2009ym.nc 
netcdf MMEAN2000-2009ym {

dimensions:

time = UNLIMITED ; // (1 currently)

bnds = 2 ;

lon = 362 ;

lat = 196 ;

lev1 = 31 ;

lev = 30 ;

可以看见模式的是196x362的

现在用cdo改下分辨率

cdo remapbil,r362x196 sst.2000-2009ym.nc  sst.2000-2009ym_remap.nc
estep [0.11s 32MB].
cdo remapbil: Bilinear weights from lonlat (480x241) to lonlat (362x196) grid, with source mask (76542)

cdo remapbil: Processed 115680 values from 1 variable over 1 timestep [0.11s 32MB].

现在我们明确一下文件的位置

obs

/data/chengxl/obs_duibiyanjiu/10yearmean/sst.2000-2009ym_remap.nc

ctrl

/data/chengxl/CAS-ESM2.0/cas-esm/run/atmocn_hist_test/run/MMEAN2000-2009ym.nc

exp1

/data/chengxl/CAS-ESM2.0-test2/run/HIST_model_test_finial/run/MMEAN2000-2009ym.nc

现在编写一下ncl 的脚本来绘图分析

begin
f_obs = addfile("/data/chengxl/obs_duibiyanjiu/10yearmean/sst.2000-2009ym_remap.nc","r")f_ctrl  = addfile("/data/chengxl/CAS-ESM2.0/cas-esm/run/atmocn_hist_test/run/MMEAN2000-2009ym.nc","r" ) f_exp1  = addfile("/data/chengxl/CAS-ESM2.0-test2/run/HIST_model_test_finial/run/MMEAN2000-2009ym.nc","r") lat = f_ctrl->lat 
lon = f_ctrl->lon var_ctrl = f_ctrl->ts(0,0,:,:)var_ctrl@lat   = lat var_ctrl@lon   = lonvar_exp1 = f_exp1->ts(0,0,:,:)var_exp1_ctrl = var_exp1 - var_ctrl 
;copy_VarMeta(var_ctrl,var_exp1_ctrl)var_obs =short2flt( f_obs->sst(0,:,:) )- 273.15
var_obs!0 = "lat"
var_obs!1 = "lon"
var_obs@lat = f_obs->lat
var_obs@lon = f_obs->lon
copy_VarCoords(var_ctrl,var_obs)
printVarSummary(var_ctrl)
printVarSummary(var_obs)
var_ctrl_obs = var_ctrl - var_obs 
copy_VarCoords(var_ctrl,var_ctrl_obs)var_exp1_obs = var_exp1 - var_obs 
;copy_VarCoords(var_ctrl , var_exp1_obs)print(f_obs->lat)
print(f_ctrl->lat)wks = gsn_open_wks( "png","sst_ctrl_obs")plot = new(2,graphic)res                      = Trueres@cnFillOn             = True                 ; turn on colorres@cnFillMode           = "RasterFill"         ; turn on raster moderes@cnLinesOn            = False                ; turn off contour linesres@mpLandFillColor      = "white"              ; draw landmasses in whiteres@cnLevelSelectionMode = "ExplicitLevels"res@cnLevels             = (/-1.8,0,5,10,15,20,22.5,25,27.5,30/)res@mpCenterLonF         = 180res@gsnDraw             = Falseres@gsnFrame            = Falseres@mpFillOn            = False            ; no need   res@cnLevelSelectionMode= "ManualLevels"   ; manual set levelsres@cnMinLevelValF      = -3.0res@cnMaxLevelValF      = 30.0res@cnLevelSpacingF     = 1.5              ; 20 contour levels        res@cnFillOn            = True             ; color fill plotres@cnFillPalette       = "BlAqGrYeOrRe"res@cnLinesOn           = Falseres@cnLineLabelsOn      = Falseres@cnInfoLabelOn       = Falseres@lbLabelBarOn        = True            ; turn off individual label bars
;
; Formatting the labelbar strings helps make the two sets of labelbars
; match better. Even though the labelbar is turned off, it is internally
;  still generated.
;res@lbLabelStrings      = sprintf("%4.1f",ispan(-30,370,15)*0.1)res@gsnLeftString       = "SST"res@gsnRightString      = "~S~o~N~C"res@gsnCenterString     = "(a)OBS"plot(0) = gsn_csm_contour_map(wks,var_obs(:,:359),res)res@cnMinLevelValF      = -5.res@cnMaxLevelValF      =  5.res@cnLevelSpacingF     =  1.delete(res@cnFillPalette) res@cnFillPalette       ="BlueWhiteOrangeRed";"nrl_sirkes"; "BlueDarkRed18";"BlueWhiteOrangeRed"    ; select a color map with white in the middle;---Formatting the labelbar strings helps make the two sets of labelbars match betterres@lbLabelStrings      := sprintf("%4.1f",ispan(-5,5,1))res@gsnCenterString     = "(b)CTRL - OBS"res@lbLabelBarOn        = True            ; turn off individual label barsplot(1) = gsn_csm_contour_map(wks,var_ctrl_obs(:,:359),res);---Panel the two sets of plots. Note no special resources need to be set.pres = Truepres@gsnMaximize    = True                ; maximize plotsgsn_panel(wks,plot,(/2,1/),pres)end

然后我发现一个棘手的难题，观测资料的坐标系和模式的坐标系不一样，

观察发现，模式是从北纬90度到南纬78度

观测再分析师从南纬90度到北纬90度，这对于作差来说根本没那么简单。

我也找不到什么NCL 的内置函数，说实话NCL真的很不灵活，然后我决定自己写一个把纬度转换下

不出意外，一下就成功了，但是南纬还是存在一些问题

其实可以只画北半球

这样稍微看起来正常一点

再看看插值方法能不能做。

翻出之前的一个插值的NCL代码

begin
dirfile = systemfunc("ls wrf*");
f = addfiles(dirfile,"r")
ListSetType(f,"cat")
t=f[:]->T2
lon_wrfout = f[:]->XLONG(0,0,:)
lat_wrfout = f[:]->XLAT(0,:,0)
t_mean = dim_avg_n_Wrap(t,0)
printVarSummary(t_mean)
f_demHigh = addfile("DemHiRes.nc","r")
lon_high = f_demHigh->lon
lat_high = f_demHigh->lat
qsort(lon_wrfout)
qsort(lat_wrfout)
qsort(lon_high)
qsort(lat_high)
t_high_regrid = linint2_Wrap(lon_wrfout,lat_wrfout,t_mean,False,lon_high,lat_high,0)
system("rm -f MaYubin_HiRes.nc")
fout = addfile("MaYubin_HiRes.nc","c")
fout->t_mean_wrfout_regrid = t_high_regrid
end

我现在想要把再分析的向模式的坐标转化

lon_obs = f_obs->lon
lat_obs = f_ctrl->lat
lon_ctrl = f_ctrl->lon
lat_ctrl = f_ctrl->lat
qsort(lon_obs)
qsort(lat_obs)
qsort(lon_ctrl)
qsort(lat_ctrl)var_regrid = linint2_Wrap(lon_obs,lat_obs,var_obs,False,lon_ctrl,lat_ctrl,0)fout = addfile("sst_reanalysis_regrid_20002009.nc","c")
fout->sst = var_regrid