气候态平均水汽通量和水汽通量散度（NCEP2)

徵言片羽

    大神们已经写过不少关于水汽通量和水汽通量散度的帖子，我也是借鉴ncl官网的例子和他人的代码，写出了下面的代码。与文献中贴出的图基本一致，可以认为是正确的。
    需要注意的是，数据的单位与结果密切相关，必须统一，否则画图时容易出错。其他问题大家一起讨论。废话不多说，直接贴代码：

load "$NCARG_ROOT/lib/ncarg/nclscripts/cnmap/cnmap.ncl"
begin  
  ptop = 300              ; 'shum' upper level                  
  ptop@units = "hPa"
  g    = 9.80665          ; m/s2

  ymStrt   = 198101     
  ymLast   = 201012

  ;---ESRL: CDC data
  diri = "D:/data/NCEP2/"

;NCEP2 monthly data
  filq = "shum.mon.mean.nc"            
  filu = "uwnd.mon.mean.nc"              
  filv = "vwnd.mon.mean.nc"
  filps= "pres.sfc.mon.mean.nc"

  pthu = diri+filu      
  pthv = diri+filv
  pthq = diri+filq     
  pthps= diri+filps   

  fu   = addfile(pthu ,"r")
  fv   = addfile(pthv ,"r")
  fq   = addfile(pthq ,"r")
  fps  = addfile(pthps,"r")

;time choice
  YYYYMM  = cd_calendar(fu->time, -1)    ; ymd: human readable
  iStrt  = ind(YYYYMM.eq.ymStrt)
  iLast  = ind(YYYYMM.eq.ymLast)

  u0    = fu->uwnd(iStrt:iLast,{1000:ptop},:,:); m/s, (time,level,lat,lon)
  v0    = fv->vwnd(iStrt:iLast,{1000:ptop},:,:)
  q0    = fq->shum(iStrt:iLast,{1000:ptop},:,:)        ; [kg/kg], 1000-300 levels only
  ps0   = fps->pres(iStrt:iLast,:,:)           ; Pa=>[kg/(m-s2)], (time,lat,lon)

;---Vertical levels
  ptop = ptop*100
  ptop@units = "Pa"

  plev = q0&level                        ; hPa  
  ;plev = plev(::-1)
  plev = plev*100                     ; [100000,...,30000] Pa [kg/(m-s2)]
  plev@units = "Pa"
  print(plev)
;---Change [kg/kg] to [g/kg]; not necessary: but common units for q
  q0    = q0*1000            
  q0@units = "g/kg"
;---Divergence function [used later] requires S->N grid order
  u1    = u0(:,:,::-1,:)  
  v1    = v0(:,:,::-1,:)
  q1    = q0(:,:,::-1,:)     
  ps1   =ps0(:,  ::-1,:)
;Calculate climatological data
  uclm  = clmMonTLLL(u1)
  vclm  = clmMonTLLL(v1)
  qclm  = clmMonTLLL(q1)
  psclm = clmMonTLL(ps1)   
;autumn average
  u = month_to_season(uclm, "SON")
  v = month_to_season(vclm, "SON")  
  q = month_to_season(qclm, "SON")  
  ps = month_to_season(psclm, "SON")     
;---Layer thickness: ; Pa=>[kg/(m-s2)], (time,level,lat,lon)
;---Mass weighting:  (dp/g) => [Pa/(m/s2)] => (Pa-s2)/m => [kg/(m-s2)][s2/m] =>  (kg/m2)
  dp   = dpres_plevel_Wrap(plev, ps, ptop, 0) ; Pa; layar thickness

  dpg  = dp/g   
  dpg@long_name = "Layer Mass Weighting"
  dpg@units     = "kg/m2"                     ; dp/g, Pa/(m s-2), reduce to kg m-2

;---Moisture flux components at each pressure level
  printVarSummary(u)
  printMinMax(u, False)
  ;print(u&lev)
  printVarSummary(q)
  printMinMax(q, False)
  ;print(q&lev)

  uq   = u*q                                  ; (:,:,:,:)
  uq@long_name = "Zonal Moisture Flux [uq]"
  uq@units = "["+u@units+"]["+q@units+"]"     ; [m/s][g/kg]     
  copy_VarCoords(u,uq)                        ; (time,level,lat,lon)

  vq   = v*q                                  ; (:,:,:,:)
  vq@long_name = "Meridional Moisture Flux [vq]"
  vq@units = "["+v@units+"]["+q@units+"]"
  copy_VarCoords(v,vq)                        ; (time,level,lat,lon)

;---Integrated mass weighted moisture flux components

  uq_dpg = uq*dpg                ; mass weighted 'uq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
  iuq    = dim_sum_n(uq_dpg, 1)
  iuq@long_name = "Integrated Zonal UQ [uq*dpg]"
  iuq@LONG_NAME = "Sum: Mass Weighted Integrated Zonal Moisture Flux [uq*dpg]"
  iuq@units     = "[m/s][g/kg]"
  copy_VarCoords(u(:,0,:,:), iuq); (time,lat,lon)
  delete(uq_dpg)

  vq_dpg = vq*dpg                ; mass weighted 'vq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
  ivq    = dim_sum_n(vq_dpg, 1)
  ivq@long_name = "Integrated Meridional VQ [vq*dpg]"
  ivq@LONG_NAME = "Sum: Mass Weighted Integrated Meridional Moisture Flux [vq*dpg]"
  ivq@units     = "[m/s][g/kg]"
  copy_VarCoords(v(:,0,:,:), ivq); (time,lat,lon)
  delete(vq_dpg)

  ;---Divergence of moisture flux: uv2dvF => global 'fixed' rectilinear grid

  duvq  = uv2dvF_Wrap(uq, vq)    ; (time,level,lat,lon)
  duvq@long_name = "Divergence of Moisture Flux"
  duvq@units     = "g/(kg-s)"    ; (1/m)*[(m/s)(g/kg)] => [g/(kg-s)]

;---Mass weighted integration [sum] of the divergence of moisture flux

  duvq_dpg = duvq*dpg            ;  [g/(kg-s)][kg/m2] => [g/(m2-s)]
  iduvq    = dim_sum_n(duvq_dpg, 1)
  iduvq@long_name = "Integrated Mass Wgt MFC"
  iduvq@LONG_NAME = "Integrated Mass Weighted Moisture Flux Convergence"
  iduvq@units     = "g/(m2-s)"
  copy_VarCoords(u(:,0,:,:), iduvq)      ; (time,lat,lon)
  delete(duvq_dpg)

  VIMFC =  iduvq           ; keep meta data                        
  ;VIMFC = -VIMFC           ; Note the preceding -1 [negative precedes integration]
  VIMFC@long_name = "VIMFC"

  PRINT_RESULT = True
if (PRINT_RESULT) then
    printVarSummary(iuq)                 ; (time,lat,lon)
    printMinMax(iuq,0)
    print("-----")
    printVarSummary(ivq)                 ; (time,lat,lon)
    printMinMax(ivq,0)
    print("-----")
    printVarSummary(duvq)                ; (time,lev,lat,lon)
    printMinMax(duvq,0)
    print("-----")
    printVarSummary(iduvq)               ; (time,lat,lon)
    printMinMax(iduvq,0)
    print("-----")
end if

  scl5  = 1e5                                  ; arbitrary: used for nicer plot values
  sclab5= "(10~S~-5~N~)"                       ; used later   
  SCLAB5= "(10~S~5~N~)"           

  scl6  = 1e6  
  sclab6= "(10~S~-6~N~)"         
  SCLAB6= "(10~S~6~N~)"         

  plot := new(2,graphic)

  wks   = gsn_open_wks("png","ncep2-climate-integrated-vapor-flux")        ; send graphics to PNG file

  resd                 = True
  resd@cnFillOn        = True                  ; color
  resd@cnLinesOn       = False                 ; turn off contour lines

  resd@cnLevelSelectionMode = "ManualLevels"   ; set manual contour levels
  resd@cnMinLevelValF  = -15.                  ; set min contour level
  resd@cnMaxLevelValF  =  15.                  ; set max contour level
  resd@cnLevelSpacingF =   1.                  ; set contour spacing
;resd@cnFillPalette   = "cmocean_balance"     ; NCL 6.5.0
  resd@cnFillPalette   = "ViBlGrWhYeOrRe"

  resd@mpFillOn        = False                 ; turn off map fill
  resd@vcRefMagnitudeF = 3.                    ; make vectors larger
  resd@vcRefLengthF    = 0.025                 ; reference vector length
  resd@vcGlyphStyle    = "CurlyVector"         ; turn on curly vectors
  resd@vcMinDistanceF  = 0.010                 ; thin the vectors
  resd@vcRefAnnoOrthogonalPosF = -1.0          ; move ref vector up
  resd@gsnLeftString   = "Divergent Wind"
  resd@gsnScalarContour= True                  ; vectors over contours

  IUQ     = iuq(0,:,:)                    ; local array: keep meta data
  IUQ     = IUQ/980                       ; scale for plot
;  res@gsnRightString  = "kg/(ms)"

  IVQ     = ivq(0,:,:)                    ; local array: keep meta data
  IVQ     = IVQ/980

  VIMFC=VIMFC*4


  resd@tiMainString    = "integrated and divergence of vapor flux"
  ;resd@gsnCenterString = "1000-300hPa"  
  resd@gsnLeftString   = "integrated vapor flux"     
  resd@gsnRightString  = "(10~S~-4~N~)kg/(m2-s)"   

  resd@cnMinLevelValF  =  -1.                  ; set min contour level
  resd@cnMaxLevelValF  =  1.                  ; set max contour level
  resd@cnLevelSpacingF =  0.1                  ; set contour spacing
  resd@cnFillScales         = 0.5
  resd@gsnLeftString   = "(a)"
  resd@gsnLeftStringFontHeightF     = 0.02
  resd@gsnLeftStringOrthogonalPosF  = -0.085
  resd@gsnLeftStringParallelPosF    = -0.075  ;左移
  ;resd@gsnCenterString  = "1981-2010 integrated vapor flux and moisture divergence"
  resd@tmXTOn   = False    ;turn down up scale
  resd@tmYROn   = False    ;turn down right scale

  resd@vcRefMagnitudeF         = 150.                    ; make vectors larger
  resd@mpMinLatF               = 0.                        
  resd@mpMaxLatF               = 55.
  resd@mpMinLonF               = 60.
  resd@mpMaxLonF               = 150.

  vpflux = gsn_csm_vector_scalar_map(wks,IUQ,IVQ,VIMFC(0,:,:),resd)

;>============================================================<
;                      add China map
;>------------------------------------------------------------<
cnres           = True
cnres@china     = True       ;draw china map or not
cnres@river     = False       ;draw changjiang&huanghe or not
cnres@province  = False       ;draw province boundary or not
cnres@nanhai    = False       ;draw nanhai or not
cnres@diqu      = False       ; draw diqujie or not

chinamap = add_china_map(wks,vpflux,cnres)
;>============================================================<

draw(vpflux)
frame(wks)

end
结果如下

呆妹小霸王

可以发一下你说的文献吗

伽蓝鸟

水汽通量用monthly数据来算是不准确的。。

徵言片羽

呆妹小霸王 发表于 2021-4-22 09:32
可以发一下你说的文献吗

Teleconnected influence of tropical Northwest Pacific sea surface temperature on interannual variability of autumn precipitation in Southwest China

徵言片羽

伽蓝鸟 发表于 2021-4-22 08:32
水汽通量用monthly数据来算是不准确的。。

谢谢提醒哈，不过我是在文献里看到作者是用月数据画的水汽通量和水汽通量散度的 气候态平均，就依样画葫芦画了一个。

Sharon_Liu

伽蓝鸟 发表于 2021-4-22 08:32
水汽通量用monthly数据来算是不准确的。。

请问如果计算多年平均夏季水汽通量，不能用月尺度数据嘛？应该用什么尺度的呢？

伽蓝鸟

Sharon_Liu 发表于 2021-5-30 21:32
请问如果计算多年平均夏季水汽通量，不能用月尺度数据嘛？应该用什么尺度的呢？

最好用六小时数据，勉强的话逐日数据也行。。

Sharon_Liu

伽蓝鸟 发表于 2021-5-31 08:59
最好用六小时数据，勉强的话逐日数据也行。。

如果算几十年的，多年季平均，如果按照日尺度的话，这个数据量实在太大了。。。感觉不太可行啊，这时候用月尺度的数据是不是也可呢？

Sharon_Liu

感谢大神分享，最近画水汽输送图帮助很大！
有一点疑问 IUQ = IUQ/980， VIMFC = VIMFC*4，这两处除以980和乘以4是什么意思呢？
代码前面有一个south to north 的转换，在画图时候是不是需要再转回north to south 呢？
希望大神帮忙解答一下谢谢！

毛豆薛

同样想知道楼上的问题