ncl计算整层水汽通量和一个指数的相关系数时青藏高原地区的资料缺失

周扬澄

感觉可能是积分的部分有问题？求助
fil = "6789.txt"
score = asciiread(fil,65,"float")

nlat = 73
nlon = 144
nyear = 65  
nlev = 8         
p = (/ 1000.,925.,850.,700.,600.,500.,400.,300./)
linlog = 1
pbot   = 1100.                                             
ptop   = 300.
u_file = addfile("uwnd.mon.mean.nc","r")
v_file = addfile("vwnd.mon.mean.nc","r")
p_file = addfile("pres.mon.mean.nc","r")
shum_file = addfile("shum.mon.mean.nc","r")

lat = u_file->lat
lon = u_file->lon
lev = u_file->level
time = u_file->time
ftime = cd_calendar(time, -1)

ts=195001  ;start
te=201412  ;end

recs = ind(ts.eq.ftime)   ;record start
rece = ind(te.eq.ftime)

u = u_file->uwnd(recs:rece,0:7,:,:)
v  = v_file->vwnd(recs:rece,0:7,:,:)
pres  = p_file->pres(recs:rece,:,:)
shum  = shum_file->shum(recs:rece,:,:,:)
printVarSummary(u)

uJJA = month_to_season(u, "JJA")
vJJA = month_to_season(v, "JJA")
shumJJA = month_to_season(shum, "JJA")
presJJA = month_to_season(pres, "JJA")
copy_VarCoords(u,uJJA)
copy_VarCoords(v,vJJA)
copy_VarCoords(shum,shumJJA)
copy_VarCoords(pres,presJJA)
printVarSummary(uJJA)
printVarSummary(shumJJA)

qu = new((/nyear,nlev,nlat,nlon/),float)
qv = new((/nyear,nlev,nlat,nlon/),float)

;do i=0,7
    qu = shumJJA*uJJA
    qv = shumJJA*vJJA
;end do
         
;do j=8,nlev+9-1
;    qu(:,j,:,:)= 0.0
;    qv(:,j,:,:)= 0.0
;end do

copy_VarCoords(uJJA,qu)
copy_VarCoords(vJJA,qv)
printVarSummary(qu)
printVarSummary(presJJA)   

uu = vibeta (p,qu(time|:,lat|:,lon|:,level|:),linlog,presJJA,pbot,ptop)/9.8  ; returns u(time,lat,lon)
vv = vibeta (p,qv(time|:,lat|:,lon|:,level|:),linlog,presJJA,pbot,ptop)/9.8  
uv = sqrt(uu^2+vv^2)
copy_VarCoords(uJJA(:,0,:,:),uu)
copy_VarCoords(vJJA(:,0,:,:),vv)
copy_VarCoords(vJJA(:,0,:,:),uv)
;u@_FillValue=-9.99E+33  
;v@_FillValue=-9.99E+33  
         
printVarSummary(uv)
; asciiwrite ("allqu.txt" , u)
;delete(uu)
;delete(vv)
;delete(pres)
;delete(shum)
;************************************************
;求相关系数
;************************************************
r = escorc(score,uv(lat|:,lon|:,time|:))
copy_VarMeta(uv(0,:,:), r)
printVarSummary(r)
;print(r)

n = dimsizes(score)
df = (n-2)
t = r*sqrt((n-2)/(1-r^2))

;************************************************
;t-test
;************************************************
;通过0.05检验的打点
pp = student_t(t, df)
psig = 0.01
;print("t="+t+"p="+p)

tn = 73*144
polylat_d_s = new(tn,double)
polylon_d_s = new(tn,double)

kds = 0
do i = 0,73-1
    do j = 0,144-1
       if (.not.ismissing(pp(i,j)).and.pp(i,j).le.psig)then
         polylat_d_s(kds) = lat(i)
         polylon_d_s(kds) = lon(j)
         kds = kds+1
       end if
    end do
end do
;print(kds)


;************************************************
;plot
;************************************************
  wks = gsn_open_wks("pdf","YRI&Moisture Transport")                ; send graphics to PNG file
  res = True                                    ; plot mods desired
  res@gsnDraw              = False              ; don't draw
  res@gsnFrame             = False              ; don't advance frame


;  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
;  res@cnMinLevelValF       = -10.               ; set min contour level
;  res@cnMaxLevelValF       =  35.               ; set max contour level
;  res@cnLevelSpacingF      =   3.               ; set contour spacing

  res@gsnLeftString ="" ;clear up the long_name informational
  res@gsnCenterString ="Correlation between YRI and Moisture Transport" ;setup plot title

;************************************************
  res@gsnMaximize           = True     ; Make plot fit the frame
  res@cnFillOn              = True     ; turn on color fill
  res@cnLinesOn             = False    ; turn of contour lines
;  res@cnLevelSpacingF       =   0.15              ; set contour spacing
  res@cnFillPalette         = "BlueRed"
;  res@cnFillPalette         = "cmp_b2r"
  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels             = (/-0.45,-0.3,-0.15,0,0.15,0.3,0.45/)
  res@gsnAddCyclic          = True    ; data already has cyclic point
;
; Note that the gsn_csm_*map* templates automatically set
; res@mpLimitMode="LatLon" for you. If you are plotting a
; different projection, you may have to set this resource.
;
  res@mpMinLatF            = -10     
  res@mpMaxLatF            = 50.
  res@mpMinLonF            = 70.
  res@mpMaxLonF            = 170.
;  res@mpCenterLonF         =  180
  res@gsnLeftString        = ""       ; clear the left string
  res@gsnRightString       = ""       ; clear the right string

   res@cnLineLabelsOn = False
   res@cnInfoLabelOn = False
   res@lbLabelBarOn = True
   res@pmTickMarkDisplayMode = "Always"
   res@cnLevelSelectionMode = "ExplicitLevels"
   res@cnFillOn             = True  

   ; res@cnFillDrawOrder="postDraw"
   ; res@cnLevels    = (/-tflag1_hl_90,tflag1_hl_90/)
   ; res@cnMonoFillPattern    = False            ; want multiple patterns
   ; res@cnFillPatterns       = (/17,-1,17/) ; the patterns
   ; res@cnMonoFillScale      = False            ; want different densities
   ; res@cnFillScales         = (/1,1,1/) ; change densities         
   ; res@cnMonoFillColor      = True
   ; res@cnFillDotSizeF       = 0.0045

   ; plot2(0) = gsn_csm_contour(wks,tvalue1_hl,res)
   plot = gsn_csm_contour_map(wks,r,res)  ; create plot
   ; overlay(plot, plot2(0))
  

   polyres                   = True
   polyres@gsMarkerOpacityF   = 1
   polyres@gsMarkerColor     = "black"
   polyres@gsMarkerIndex     = 16          ; polymarker style
   polyres@gsMarkerSizeF     = 0.005         ; polymarker size


   dum_s     = gsn_add_polymarker(wks,plot,polylon_d_s,polylat_d_s,polyres)


  opt           = True
  opt@color     = "blue"
  opt@thickness = 4
  ;opt@dashP    = 1
  add_Changjiang_Huanghe(wks,plot,opt)

  draw(plot)
  frame(wks)



end
D:\chengxu\YRI&Moisture Transport

愛、星晴

我说一下我猜测的可能原因，我看到你的pbot设置的是1100hPa，其实青藏高原的平均海拔最低都得600hPa,这样在高原地区如果从1100hPa开始积分会加入很多缺测值，可能会对你积分函数的结果造成影响。

shuaguaxyy

楼主你可以参考一下这个帖子
http://bbs.06climate.com/forum.p ... E%C6%FB%CD%A8%C1%BF

帅帅f

学习学习··