纬圈环流（Walker）结果异常

顺利毕业！

更新：
已解决
本质上是理解公式有误，具体可参考一个文献，提到了Hadley和Walker环流的公式，比较详细（Variability of the extent of the Hadley circulation in the southern hemisphere: a regional perspective）
计算环流前需要先对矢量风即U和V进行处理，主要用到两个函数：uv2dvF_Wrap和dv2uvF_Wrap，具体用法可参考NCL官网
如果变量存在缺测值，则需要先利用NCL的poisson_grid_fill函数把缺测值补上（上述两个函数不允许缺测）
主要注意点就这些，希望能帮到同样没有大气动力学基础的小白，如果有不足之处还请各位前辈指出~
;;;;===================

最近在尝试分析Walker环流，打算用流函数来表示。
看到一个博主分享了两种方案来计算经圈环流（Hadley）：一是NCL自带的 zonal_mpsi函数，二是自己做垂直积分计算。这样还可以用于计算纬圈流函数。（博主链接如下：辐散风和旋转风 (renqlsysu.github.io)）ps：博主好多可以学习的帖子，赞！！！
于是我就按照方案二去写NCL脚本，但是气候态的Walker环流就和正常情况下（西太上升流，东太下降流）不一致，自查了好久都没发现问题，还请各位前辈帮忙看看代码是哪里出了问题！！！

begin

;;;;矢量风数据读取及处理
f1             = addfile("./***_U_V_chazhi.nc", "r")                          ;气候态数据
f2             = addfile("./***_PS.nc", "r")         
U_ctrl         = f1->U_ctrl_vinth2p(:,:,{100:1000},{-5:5},:)   
Uc_ave         = dim_avg_n_Wrap(U_ctrl,0)
Uc_JJA         = month_to_season(Uc_ave, "JJA")
Uc_JJAave      = dim_avg_n_Wrap(Uc_JJA,0)                   ;JJA平均,（lev_p,lat,lon）; lev_p is top-to-bottom (hpa)
Uc_Za          = dim_avg_n_Wrap(Uc_JJAave,1)                ;对lat纬进行平均
Uc_Za          = Uc_Za(lev_p|::-1,lon|:)                            ;lev_p is bottom-to-top (hpa)               

f3             = addfile("./*****_U_V_chazhi.nc", "r")                             ;强迫数据，分析的是JJA1
f4             = addfile("./*****_PS.nc", "r")  

U_foc         = f3->U_forcing_vinth2p(:,12:23,{100:1000},{-5:5},:)        ;取JJA1
Uf_ave         = dim_avg_n_Wrap(U_foc,0)
Uf_JJA         = month_to_season(Uf_ave, "JJA")
Uf_JJAave      = Uf_JJA(0,:,:,:)                                                          ;JJA1 lev_p,lat,lon                              
Uf_Za          = dim_avg_n_Wrap(Uf_JJAave,1)                                   ;对lat纬进行平均
Uf_Za          = Uf_Za(lev_p|::-1,lon|:)

;;;;垂直积分
pi             = atan(1.0)*4                                ;圆周率
a              = 6378388.                                   ;地球半径
g              = 9.80665
lat            = f1->lat({-5:5})
lon            = f1->lon
plev           = (/1000,975,950,925,900,875,850,825,800,775,750,700,650,600,550,500,450 \
,400,350,300,250,225,200,175,150,125,100/)*100             ;单位为Pa
ptop           = min(plev)
ps_ctrl        = f2->PS_ctrl(:,:,{-5:5},:)                             ;单位为Pa
pc_ave         = dim_avg_n_Wrap(ps_ctrl,0)
pc_JJA         = month_to_season(pc_ave, "JJA")
pc             = dim_avg_n_Wrap(pc_JJA,0)
dp_ctrl        = dpres_plevel_Wrap(plev, pc, ptop, 0)
dpc            = dim_avg_n_Wrap(dp_ctrl,1)                  ;对lat纬进行平均

nlev           = dimsizes(plev)
nlon           = 144
msf_c          = new((/nlev,nlon/), float)
do i=0,nlev-1,1
    msf_c(i,:) = dim_sum_n(Uc_Za(i:nlev-1,:)*dpc(i:nlev-1,:),0)   
    msf_c(i,:) = msf_c(i,:)*2*pi*a/g
end do  

Uc             = msf_c /10000000000.                        ; *10的-10次 方便画图
copy_VarCoords(Uc_Za, Uc)

ps_foc         = f4->PS_foc(:,12:23,{-5:5},:)               ;JJA1
pf_ave         = dim_avg_n_Wrap(ps_foc,0)
pf_JJA         = month_to_season(pf_ave, "JJA")
pf               = pf_JJA(0,:,:)
dp_foc         = dpres_plevel_Wrap(plev, pf, ptop, 0)
dpf            = dim_avg_n_Wrap(dp_foc,1)                  ;对lat纬进行平均

msf_f          = new((/nlev,nlon/), float)
do i=0,nlev-1,1
    msf_f(i,:) = dim_sum_n(Uf_Za(i:nlev-1,:)*dpf(i:nlev-1,:),0)
    msf_f(i,:) = msf_f(i,:)*2*pi*a/g
end do

Uf             = msf_f /10000000000.
copy_VarCoords(Uf_Za,Uf)

Ua             = Uf - Uc
copy_VarMeta(Uc, Ua)

olev_p          = Uc_Za&lev_p
olev_p!0        = "lev_p"
olev_p@longname = "lev_p"
olev_p@units    = "hPa"
olev_p&lev_p    = olev_p
Ua!0            = "lev_p"
Ua&lev_p        = olev_p

;;;;draw
wks = gsn_open_wks("eps","./Walker_5S-5N_JJA1")
gsn_define_colormap(wks,"MPL_RdBu")

res                              = True
res@pmTickMarkDisplayMode        = "Always"         ;显示经纬度°
res@trYReverse                   = True
res@gsnDraw                      = False
res@gsnFrame                     = False
res@tmYROn                       = False
res@tmXTOn                       = False
res@trYMinF                      = 100

res@vpHeightF                    = 0.3
res@vpWidthF                     = 0.7

res@tmXBMode                     = "Explicit"
res@tmXBValues                   = (/"60","120","180","240","300"/)
res@tmXBLabels                   = (/"60~S~o~N~ E","120~S~o~N~ E","180~S~o~N","120~S~o~N~ W","60~S~o~N~ W"/)
res@tmXBLabelDeltaF              = -0.5
res@tmYLLabelDeltaF              = -0.4
res@tmYLMode                     = "Explicit"
res@tmYLValues                   = (/"1000","700","500","400","300","200","100"/)
res@tmYLLabels                   = (/"1000","700","500","400","300","200","100"/)

res@tiMainString                 = " "
res@tiYAxisString                = "Level (hPa)"
res@tiYAxisFontHeightF           = 0.022
res@tiYAxisOffsetXF              = 0.01                             ;平行于x轴移动标签
res@gsnLeftString                = "(a) Walker circulation"
res@gsnLeftStringOrthogonalPosF = -0.01       ;调节字体上下位置
res@gsnLeftStringFontHeightF    = 0.02     ;调节字体大小
res@gsnRightString               = "Tam1850_JJA1"
res@gsnRightStringOrthogonalPosF  = -0.005
res@gsnRightStringFontHeightF     = 0.02

res@tmXBMajorLengthF             = 0.008             ;设置主刻度线的长短
res@tmYLMajorLengthF             = 0.008
res@tmXBMajorOutwardLengthF      = 0.008
res@tmYLMajorOutwardLengthF      = 0.008
res@tmYROn                       = False
res@tmXTOn                       = False
res@tmXBLabelFontHeightF         = 0.02
res@tmYLLabelFontHeightF         = 0.02

res1                              = res
res1@cnFillOn                     = True
res1@cnLinesOn                    = False
res1@cnLineLabelsOn               = False
res1@lbLabelBarOn                 = True
res1@cnLevelSelectionMode         = "ManualLevels"
res1@cnMinLevelValF               = -40.                 ; -4~4 0.5 -14~14 1.
res1@cnMaxLevelValF               = 40.
res1@cnLevelSpacingF              = 4
res1@tmYRMode                     = "Automatic"          ; turn off special labels on right axis
;;lb
res1@gsnSpreadColorStart         = 129  ;选取色表颜色
res1@gsnSpreadColorEnd           = 2
res1@pmLabelBarHeightF            = 0.1      ;色标长宽
res1@pmLabelBarWidthF             = 0.85
;res1@lbBoxLinesOn                 = False
res1@lbBoxLineColor               = "black"    ;色标边框颜色
;res1@lbBoxLineThicknessF          = 0.5        ;色标边框粗细
res1@pmLabelBarOrthogonalPosF     = 0.05      ;色标垂直方向的位置
res1@pmLabelBarParallelPosF       = 0.57      ;色标水平方向的位置
res1@lbLabelFontHeightF           = 0.02       ;设置色标数字的字体大小
res1@lbLabelStride                = 2           ;数字间隔
res1@lbTitleOn                    = True
res1@lbTitleString                = "10~S~10~N~ kg/s"
res1@lbTitlePosition              = "Right"
res1@lbTitleDirection             = "Across"
res1@lbTitleJust                  = "CenterLeft"
res1@lbTitleFontHeightF           = 0.02

plot  = gsn_csm_pres_hgt(wks, Ua, res1)   ; plaace holder

res2                              = res
res2@cnFillOn                     = False
res2@cnLinesOn                    = True
res2@cnLineLabelsOn               = False
res2@lbLabelBarOn                 = False
res2@gsnContourNegLineDashPattern = 11
res2@gsnContourPosLineDashPattern = 0
res2@cnInfoLabelOn                = False
res2@tmYLOn                       = False
res2@tmXBOn                       = False
res2@cnLineDrawOrder              = "PostDraw"
plot1                             = gsn_csm_contour(wks,Uc,res2)

overlay(plot, plot1)
draw(plot)
frame(wks)



end

kongdechao

友友，求分享更正后的代码