NCL如何画wrfout水汽通量及水汽通量散度，分享一下脚本

lsk123

请教一下，可以分享一下NCL画wrfout水汽通量及水汽通量散度的脚本吗？谢谢大家

manzx

begin

     fshum          = addfile("specifichumidity.nc","r")
     fu             = addfile("u.nc","r")
     fv             = addfile("v.nc","r")
     fsp            = addfile("surfacepressure.nc", "r")
     
     pp             = fu->level({1000:300})     ;hPa
     p              = pp*100                    ;Pa
     u              = short2flt(fu->u)          ;m/s
     v              = short2flt(fv->v)          ;m/s
     q              = short2flt(fshum->q)       ;kg/kg
     ps             = short2flt(fsp->sp)        ;Pa=(kg/(m*s2))

     time           = fu->time
     timeARR        = cd_calendar(time, -5)
     yr_arr         = timeARR(:,0)
     mon_arr        = timeARR(:,1)

; for units g/kg, q should be multiplied 1000
; however, the water flux is supposed to be down to 10^-5
     q              = q*0.01
     q@units        = "g/kg"                    ;g/kg

;================================================================================================
;================================================================================================
;===================================================================================================
;--------------计算水汽通量 moisture flux=(dp(Pa)*u(m/s)*q(kg/kg))/g(m/s^2)  ----------------------



;---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)
;---Reference: http://www.cgd.ucar.edu/cas/catalog/newbudgets/
     ptop           = 30000                      ; Pa
     g              = 9.80665                    ; m/s2
     dp             = dpres_plevel_Wrap(p, ps, ptop, 0)   ;layer pressure thickness
                                                          ;这里变量单位必须全部统一成pa或者hpa
     ;[time|225]x[plev|8]x[latitude|161]x[longitude|281]
     dpg            = dp/g   
     dpg@long_name  = "Layer Mass Weighting"     ; 层质量权重
     dpg@units      = "kg/m2"     
     copy_VarCoords(u,dpg)
     ;[time|225]x[level|8]x[latitude|161]x[longitude|281]


;---Moisture flux components at each pressure level  每层水汽通量分量
     uq             = u*q                                  ; (:,:,:,:)
     uq@long_name   = "Zonal Moisture Flux [uq]"           ; 纬向水分通量
     uq@units       = "["+u@units+"]["+q@units+"]"         
     copy_VarCoords(u,uq)
     ;[time|225]x[level|8]x[latitude|161]x[longitude|281]

     vq             = v*q                                  ; (:,:,:,:)
     vq@long_name   = "Meridional Moisture Flux [vq]"      ; 纬向水分通量
     vq@units       = "["+v@units+"]["+q@units+"]"
     copy_VarCoords(v,vq)      
     ;[time|225]x[level|8]x[latitude|161]x[longitude|281]


;---Integrated mass weighted moisture flux components
;iuq for: integrated u*q   [time|225]x[latitude|161]x[longitude|281]
;ivq for: integrated v:q   [time|225]x[latitude|161]x[longitude|281]
;    按层质量权重做积分，在level维度做水汽分量的累加积分
     uq_dpg         = uq*dpg                                
                    ; mass weighted 'uq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
     iuq            = dim_sum_n(uq_dpg, 1)                 ;kg/m*s
     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)
     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)                 ;kg/m*s
     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)
     delete(vq_dpg)


     printMinMax(iuq,0)
     printMinMax(ivq,0)

end

lsk123

manzx 发表于 2023-11-21 19:21
begin

     fshum          = addfile("specifichumidity.nc","r")

好的，谢谢大佬，这个.nc全部换成wrfout.nc就好了吗？wrfout数据文件格式好像有点不太一样

manzx

lsk123 发表于 2023-11-22 10:23
好的，谢谢大佬，这个.nc全部换成wrfout.nc就好了吗？wrfout数据文件格式好像有点不太一样

嘶，我没理解这个问题的意思？
我这个脚本用的是era5的数据不是wrfout，变量要找到wrfout里对应的比湿、UV风和气压

lsk123

manzx 发表于 2023-11-22 17:10
嘶，我没理解这个问题的意思？
我这个脚本用的是era5的数据不是wrfout，变量要找到wrfout里对应的比湿、 ...

嗯嗯，好的，感谢感谢

95niande77

manzx 发表于 2023-11-21 19:21
begin

     fshum          = addfile("specifichumidity.nc","r")

感恩！{:eb511:}{:eb511:}{:eb511:}{:eb511:}{:eb511:}