立即注册 登录
气象家园 返回首页

籽悠的个人空间 http://bbs.06climate.com/?30355 [收藏] [复制] [分享] [RSS]

日志

ncl 请教

已有 102 次阅读2017-3-18 10:46

一开始都能运行成功,但改个语句就不行了。不知道哪里出问题了。只能换重新复制个文本再运行,但总不是办法,请教!文本如下

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
GPP     = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/GPP-obs-2000-2004-avg.nc","r")
GPPrat  = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/vegfrac.nc","r")
oz  = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/ozone.nc","r")

noveg  = GPPrat->noveg
nt     = GPPrat->nt
bt     = GPPrat->bt
c3c4   = GPPrat->c3c4
lon    = GPPrat->lon
lat    = GPPrat->lat

gppj        = GPP->gppj
land        = GPP->landmask
area        = GPP->area

time       = oz->time
date       = oz->date
olat        = oz->lat
olon        = oz->lon
ted        = oz->time_bnds
O3         = oz->O3_SRF

gppj    = where(gppj.ne.0,gppj,gppj@_FillValue)
bt!0 = "lat"
bt!1 = "lon"
bt&lat = lat
bt&lon = lon
bt&lat@units   = "degrees_north"
bt&lon@units   = "degrees_east"

bt1      = dble2flt(bt)
gppj1    = gppj
gppj1    = gppj*bt1

c3c4!0 = "lat"
c3c4!1 = "lon"
c3c4&lat = lat
c3c4&lon = lon
c3c4&lat@units   = "degrees_north"
c3c4&lon@units   = "degrees_east"

c3c41      = dble2flt(c3c4)
gppj2    = gppj
gppj2    = gppj*c3c41


nt!0 = "lat"
nt!1 = "lon"
nt&lat = lat
nt&lon = lon
nt&lat@units   = "degrees_north"
nt&lon@units   = "degrees_east"
nt1       = dble2flt(nt)
gppj3    = gppj
gppj3    = gppj*nt1

noveg!0 = "lat"
noveg!1 = "lon"
noveg&lat = lat
noveg&lon = lon
noveg&lat@units   = "degrees_north"
noveg&lon@units   = "degrees_east"
noveg1      = dble2flt(noveg)
gppj4    = gppj
gppj4    = gppj*noveg1

zong0    = gppj
zong0    = gppj1+gppj2+gppj3+gppj4

O3             = O3*1000000000
O3@units   = "ppb"
oa             = O3

;-----不同PFT年损伤率-----
tnu    = 12
bf0    = 1
bf1    = 0.9006
bf2    = 0.7681
bf3    = 0.7915
do i  = 0,tnu-1
oa(i,:,:)   = where(O3(i,:,:).lt.30,bf0,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.30 .and. O3(i,:,:).lt.60,bf1,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.60 .and. O3(i,:,:).lt.90,bf2,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.90,bf3,oa(i,:,:))
end do


o1         = dim_avg_n_Wrap(oa,0)
LAT    = fspan(-90,90,96)
LON    = fspan(0,357.5,144)
oa1     = area_hi2lores_Wrap(o1&lon,o1&lat,o1,True,1,LON,LAT,False)

tnu    = 12
cf0    = 1
cf1    = 0.8723
cf2    = 0.8459
cf3    = 0.8200
do i  = 0,tnu-1
oa(i,:,:)   = where(O3(i,:,:).lt.30,cf0,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.30 .and. O3(i,:,:).lt.60,cf1,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.60 .and. O3(i,:,:).lt.90,cf2,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.90,cf3,oa(i,:,:))
end do
o2         = dim_avg_n_Wrap(oa,0)
LAT       = fspan(-90,90,96)
LON      = fspan(0,357.5,144)
oa2       = area_hi2lores_Wrap(o2&lon,o2&lat,o2,True,1,LON,LAT,False)

tnu    = 12
nf0    = 1
nf1    = 1.0479
nf2    = 1.0398
nf3    = 0.8740
do i  = 0,tnu-1
oa(i,:,:)   = where(O3(i,:,:).lt.30,nf0,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.30 .and. O3(i,:,:).lt.60,nf1,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.60 .and. O3(i,:,:).lt.90,nf2,oa(i,:,:))
oa(i,:,:)   = where(O3(i,:,:).gt.90,nf3,oa(i,:,:))
end do
o3         = dim_avg_n_Wrap(oa,0)
LAT     = fspan(-90,90,96)
LON     = fspan(0,357.5,144)
oa3      = area_hi2lores_Wrap(o3&lon,o3&lat,o3,True,1,LON,LAT,False)

gppj1   = gppj1/oa1
gppj2   = gppj2/oa2
gppj3   = gppj3/oa3
gppj4   = gppj4/oa1

zong1      = gppj
land_only = gppj
cha        = gppj
zong1    = gppj1+gppj2+gppj3+gppj4
cha      = zong0-zong1
land_only   = mask(cha,land,1)
wks         = gsn_open_wks("png","cha")
cmap       = read_colormap_file("BlAqGrYeOrRe")         ; read colormap file
res          = True
   res@cnFillPalette           = cmap(20:76,:)
   res@cnFillOn                 = True     ; turn on color fill
  res@cnLinesOn               = False    ; turn off contour lines
    res@gsnLeftString         = "diff(MTE_GPP&GPP~B~pot~N~)"
  res@gsnRightString         = "g C m~S~-2"
  plot       = gsn_csm_contour_map(wks,land_only,res)

  draw(plot)
  frame(wks)
end

评论 (0 个评论)

facelist doodle 涂鸦板

您需要登录后才可以评论 登录 | 立即注册

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

返回顶部