NCL计算Heat budget，出的图很奇怪，本人已死，请求大神帮助

燕灵飞雪

计算heat budget，用的资料是SODA2.2.4，时间为1960-2008年，经度为25.75°S-25.75°N，纬度为119,75°E-299.75°W。
所用公式是图片上的公式，程序如下：

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"

begin
  f1=addfile("E:/work/data/SODA2.2.4/temp.nc","r")
  f2=addfile("E:/work/data/SODA2.2.4/uwnd.nc","r")
  f3=addfile("E:/work/data/SODA2.2.4/vwnd.nc","r")
  f4=addfile("E:/work/data/SODA2.2.4/w.nc","r")
  temp=f1->temp(:,:,11:92,60:341)                                      ;经纬度取了20°S-20°N，150°E-70°W范围。经向上共82个格点，纬向上共282个格点
  uwnd=f2->u(:,{5:45},11:92,60:341)            
  vwnd=f3->v(:,{5:45},11:92,60:341)                                  ;u、v的层次取5-45m，即混合层的上层
  w=f4->w(:,{45},11:92,60:341)                                           ;w的层次l取45m
  t1=dim_avg_n_Wrap(temp(:,{5:45},:,:),1)                         ;计算5-45m上t的平均
  t2=dim_avg_n_Wrap(temp(:,{45:85},:,:),1)                       ;计算45-85m上t的平均
  u=dim_avg_n_Wrap(uwnd,1)                                              ;计算5-45m上u的平均
  v=dim_avg_n_Wrap(vwnd,1)                                              ;计算5-45m上v的平均

  t1ave=clmMonTLL(t1)
  t1a=calcMonAnomTLL(t1,t1ave)
  t2ave=clmMonTLL(t2)
  t2a=calcMonAnomTLL(t2,t2ave)
  uave=clmMonTLL(u)
  ua=calcMonAnomTLL(u,uave)
  vave=clmMonTLL(v)
  va=calcMonAnomTLL(v,vave)
  wave=clmMonTLL(w)
  wa=calcMonAnomTLL(w,wave)

;-----------------------------------------------------计算差分----------------------------------------------------------
  lon=t1ave(0,0,:)
  lat=t1ave(0,:,0)
  
  dTavedX = new ((/12,82,282/), typeof(t1ave), t1ave@_FillValue)
  dTadX = new ((/588,82,282/), typeof(t1a), t1a@_FillValue)
  dlon = 0.5*0.0174533
  do j=0,81
  do i=1,280        
    dX = 6378388.*cos(0.0174533*lat(j))*dlon
    dTavedX(:,j,0)=(t1ave(:,j,1)-t1ave(:,j,0))/dX
    dTavedX(:,j,i)=(t1ave(:,j,i+1)-t1ave(:,j,i-1))/(2*dX)
    dTavedX(:,j,281)=(t1ave(:,j,281)-t1ave(:,j,280))/dX
    dTadX(:,j,0)=(t1a(:,j,1)-t1a(:,j,0))/dX
    dTadX(:,j,i)=(t1a(:,j,i+1)-t1a(:,j,i-1))/(2*dX)
    dTadX(:,j,281)=(t1a(:,j,281)-t1a(:,j,280))/dX
  end do
  end do

  dTavedY = new ((/12,82,282/), typeof(t1ave), t1ave@_FillValue)
  dTadY = new ((/588,82,282/), typeof(t1a), t1a@_FillValue)
  dY = 6378388.*(51.5/103*0.0174533)
  do i=0,281
  do j=1,80
    dTavedY(:,0,i)=(t1ave(:,1,i)-t1ave(:,0,i))/dY
    dTavedY(:,j,i)=(t1ave(:,j+1,i)-t1ave(:,j-1,i))/(2*dY)
    dTavedY(:,81,i)=(t1ave(:,81,i)-t1ave(:,80,i))/dY
    dTadY(:,0,i)=(t1a(:,1,i)-t1a(:,0,i))/dY
    dTadY(:,j,i)=(t1a(:,j+1,i)-t1a(:,j-1,i))/(2*dY)
    dTadY(:,81,i)=(t1a(:,81,i)-t1a(:,80,i))/dY
  end do
  end do
  copy_VarMeta(t1ave, dTavedX)
  copy_VarMeta(t1ave, dTavedY)
  copy_VarMeta(t1a, dTadX)
  copy_VarMeta(t1a, dTadY)
;--------------------------------------------------------------------------------------------------------------------------

;-------------------------------------------------------dT/dt------------------------------------------------------------
  dt=t1a
  do i=0,48
  do j=1,10
  k=12*i+j
  dt(12*i,:,:)=t1a(12*i+1,:,:)-t1a(12*i,:,:)
  dt(k,:,:)=(t1a(k+1,:,:)-t1a(k-1,:,:))/2
  dt(12*i+11,:,:)=t1a(12*i+11,:,:)-t1a(12*i+10,:,:)
  end do
  end do
  copy_VarMeta(t1a, dt)
;--------------------------------------------------------------------------------------------------------------------------

;-----------------------------------------------------线性部分----------------------------------------------------------
  uatc=ua
  vatc=va
  ucta=ua
  vcta=va
  watc=wa
  wcta=wa  
  do i=0,11
  do j=0,48
  k=12*j+i
    uatc(k,:,:) = -ua(k,:,:)*dTavedX(i,:,:)
    vatc(k,:,:) = -va(k,:,:)*dTavedY(i,:,:)
    ucta(k,:,:) = -uave(i,:,:)*dTadX(k,:,:)
    vcta(k,:,:) = -vave(i,:,:)*dTadY(k,:,:)
    watc(k,:,:) = -wa(k,:,:)*((t1ave(i,:,:)-t2ave(i,:,:))/40)
    wcta(k,:,:) = -wave(i,:,:)*((t1a(k,:,:)-t2a(k,:,:))/40)
  end do
  end do
  copy_VarMeta(t1a, uatc)
  copy_VarMeta(t1a, vatc)
  copy_VarMeta(t1a, ucta)
  copy_VarMeta(t1a, vcta)
  copy_VarMeta(t1a, watc)
  copy_VarMeta(t1a, wcta)
;------------------------------------------------------------------------------------------------------------------------
  
;----------------------------------------------非线性部分------------------------------------------------------------
  uata=ua
  vata=va
  wata=wa
  do i=0,587
  uata(i,:,:) = -ua(i,:,:)*dTadX(i,:,:)
  vata(i,:,:) = -va(i,:,:)*dTadY(i,:,:)  
  wata(i,:,:) = -wa(i,:,:)*((t1a(i,:,:)-t2a(i,:,:))/40)
  end do
  copy_VarMeta(t1a, uata)
  copy_VarMeta(t1a, vata)
  copy_VarMeta(t1a, wata)
;-------------------------------------------------------------------------------------------------------------------------

;-----------------------------------------w<0的垂直项部分值设成0------------------------------------------------
  do i=0,587
  do j=0,81
  do k=0,281
    if( .not. ismissing(w(i,j,k))) then               ;先排除w中的缺测值
      if(w(i,j,k) .lt. 0) then
        wcta(i,j,k)=0
        watc(i,j,k)=0
        wata(i,j,k)=0
      end if
    end if
  end do
  end do
  end do
;---------------------------------------------------------------------------------------------------------------------------

  fout=addfile("E:/work/data/SODA2.2.4/dt.nc","c")
  fout->T=dt

  fout=addfile("E:/work/data/SODA2.2.4/uatc.nc","c")
  fout->UATC=uatc
  fout=addfile("E:/work/data/SODA2.2.4/vatc.nc","c")
  fout->VATC=vatc
  fout=addfile("E:/work/data/SODA2.2.4/ucta.nc","c")
  fout->UCTA=ucta
  fout=addfile("E:/work/data/SODA2.2.4/vcta.nc","c")
  fout->VCTA=vcta
  fout=addfile("E:/work/data/SODA2.2.4/watc.nc","c")
  fout->WATC=watc
  fout=addfile("E:/work/data/SODA2.2.4/wcta.nc","c")
  fout->WCTA=wcta  
  
  fout=addfile("E:/work/data/SODA2.2.4/uata.nc","c")
  fout->UATA=uata
  fout=addfile("E:/work/data/SODA2.2.4/vata.nc","c")
  fout->VATA=vata
  fout=addfile("E:/work/data/SODA2.2.4/wata.nc","c")
  fout->WATA=wata



  
end

guizi

要是 时间长度，可以变成1980年-2014年，建议改用 godas 数据，SODA数据 我做的时候，确实 有问题，余项太大，结果错的离谱，并且 godas 自己的热通量有问题，建议ncep 的热通量，会好很多

燕灵飞雪

为了检验算的结果对不对，我就想挑97年典型的厄尔尼诺事件画一下它的非线性分布，结果是这样的

燕灵飞雪

正确的应该是这样的，从人家文献里截下来的：

燕灵飞雪

下面是画97年nonlinear heating的程序：

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"

begin
  f1=addfile("E:/work/data/SODA2.2.4/uata.nc","r")
  f2=addfile("E:/work/data/SODA2.2.4/vata.nc","r")
  f3=addfile("E:/work/data/SODA2.2.4/wata.nc","r")
  uata=f1->UATA
  vata=f2->VATA
  wata=f3->WATA
  nonlin=uata*1e6+vata*1e6+wata*10      ;这边的量级是一个很大的问题T_T，我总是不知道该乘以多少。据我输出数据看，uata和vata是1*10-7~10-9量级，wata是1*10-2~10-5量级。所以我把uata\vata乘了个1e6，wata乘了个10
  copy_VarMeta(uata,nonlin)
  nonlinseason=month_to_seasonN(nonlin, (/"DJF","MAM","JJA","SON"/))  
  
  wks = gsn_open_wks("ps","E:/work/asymmetry/97elnino")
  res = True
  res@gsnDraw = False ; don't draw
  res@gsnFrame = False ; don't advance frame
  res@gsnAddCyclic = False
  res@gsnRightString = ""; 关闭右边字符
  res@gsnLeftString = "" ; 关闭左边字符
  res@mpLimitMode = "LatLon"
  res@mpMinLatF = -20         
  res@mpMaxLatF = 20
  res@mpMinLonF =160
  res@mpMaxLonF =270
  res@mpCenterLonF=179.5
  res@mpFillOn = False   ;关闭地图填色
  cnres=res
  cnres@vpWidthF  = 0.9         
  cnres@vpHeightF = 0.5
  cnres@tiMainFont = "helvetica"
  cnres@tiMainOffsetYF = 0.02  
  cnres@tiMainFontHeightF = 0.02
  cnres@cnFillOn       = True
  cnres@cnLinesOn      = True
  cnres@pmLabelBarOrthogonalPosF = 0.05 ; 移动
  cnres@lbLabelFontHeightF = 0.01     
  cnres@cnLineLabelsOn = True
  cnres@cnLineLabelDensityF = 5     
  cnres@cnLineLabelFontHeightF = 0.01     
  cnres@cnFillDrawOrder = "PreDraw"
  cnres@cnLevelSelectionMode = "ExplicitLevels"   
  cnres@cnLevels    = (/-0.8,-0.6,-0.4,-0.2,0.,0.2,0.4,0.6,0.8,1.0,1.2,1.4/)
  cnres@cnFillColors = (/48,64,79,100,127,135,150,175,191,207,223,230,240/)
  cnres@lbLabelBarOn = True  
  cnres@gsnCenterStringFontHeightF=0.03
  cnres@cnInfoLabelOn  = False    ;去掉countour from ... to ...
  ;cnres@lbOrientation = "vertical"  
  cnres@gsnCenterString = "nonlinear heating 1997"
  
   map=gsn_csm_contour_map(wks,nonlinseason(0,37,:,:),cnres)  
   draw(map)  
   
  
end

燕灵飞雪

求求大家帮我看看哪里有问题下周就要汇报可是老是这种图我已经要崩溃了

weiguang1233

加油{:5_213:}{:5_213:}{:5_213:}

燕灵飞雪

weiguang1233 发表于 2015-12-13 21:22
加油

光哥啊。。我眼巴巴地看好不容易回复多了一条。。

风往北吹

{:eb502:}

二爷名声在外

@兰溪之水 黄董有时间麻烦帮忙看一下，么么哒

燕灵飞雪

二爷名声在外 发表于 2015-12-13 22:13
@兰溪之水 黄董有时间麻烦帮忙看一下，么么哒

谢二爷！感激啊