模式结果和观测的对比

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

(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   = lon

var_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                      = True
  res@cnFillOn             = True                 ; turn on color
  res@cnFillMode           = "RasterFill"         ; turn on raster mode
  res@cnLinesOn            = False                ; turn off contour lines
  res@mpLandFillColor      = "white"              ; draw landmasses in white
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels             = (/-1.8,0,5,10,15,20,22.5,25,27.5,30/)

  res@mpCenterLonF         = 180

   
    res@gsnDraw             = False
    res@gsnFrame            = False
    res@mpFillOn            = False            ; no need   
    res@cnLevelSelectionMode= "ManualLevels"   ; manual set levels
    res@cnMinLevelValF      = -3.0
    res@cnMaxLevelValF      = 30.0
    res@cnLevelSpacingF     = 1.5              ; 20 contour levels        
    res@cnFillOn            = True             ; color fill plot
    res@cnFillPalette       = "BlAqGrYeOrRe"
    res@cnLinesOn           = False
    res@cnLineLabelsOn      = False
    res@cnInfoLabelOn       = False
    res@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 better
    res@lbLabelStrings      := sprintf("%4.1f",ispan(-5,5,1))
     res@gsnCenterString     = "(b)CTRL - OBS"

    res@lbLabelBarOn        = True            ; turn off individual label bars
    plot(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 = True
    pres@gsnMaximize    = True                ; maximize plots
    gsn_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

成功了一部分，先这样吧。由于资料差异有点大，我还是只给出观测，就不给差别好了。