Does anyone know how to calculate the vertical integrated moisture flux using NCL in WRF. It would be best if you can share a NCL script. Here is my code, I used this code to calculate the averaged moisture flux in a month. but I am confused the unit of the last output. Another question is that I am not sure if I need to divide the output by 9.8?
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
dir1 = "/Users/gerry/Desktop/WRF_EXTPCP/wrftest/"
FILES1 = systemfunc (" ls " + dir1 + "wrfout_d01_* ")
a1 = addfiles(FILES1+".nc","r")
;-------calculate total layer vapor flux--------------------------
psfc = wrf_user_getvar(a1,"PSFC",-1)*0.01 ; surface pressure, turn to hPa (; unit of PSFC in WRF is Pa)
pres = wrf_user_getvar(a1,"pressure",-1) ; air pressure hPa (time,level,lat,lon)
temp = wrf_user_getvar(a1,"tk",-1) ; K
rh = wrf_user_getvar(a1,"rh",-1) ; %
ua = wrf_user_getvar(a1,"ua",-1) ; m s-1
va = wrf_user_getvar(a1,"va",-1) ; m s-1
; printVarSummary(psfc) ; (time,lat,lon)
; printVarSummary(pres) ; (time,level,lat,lon)
; printVarSummary(temp) ; (time,level,lat,lon)
; printVarSummary(rh) ; (time,level,lat,lon)
; printVarSummary(ua) ; (time,level,lat,lon)
; printVarSummary(va) ; (time,level,lat,lon)
psfc1 = dim_avg_n_Wrap(psfc,0)
pres1 = dim_avg_n_Wrap(pres,0)
temp1 = dim_avg_n_Wrap(temp,0)
rh1 = dim_avg_n_Wrap(rh,0)
ua1 = dim_avg_n_Wrap(ua,0)
va1 = dim_avg_n_Wrap(va,0)
;https://www.ncl.ucar.edu/Document/Functions/Built-in/mixhum_ptrh.shtml
shum = mixhum_ptrh(pres1,temp1,rh1,2) ;units: (hPA,K,%,kg/kg);calculate the specific humidity(2),(1)means mixing ratio
printVarSummary(shum)
uq = ua1 * shum ;vapor flux
vq = va1 * shum
copy_VarCoords(psfc(0,:,:), psfc1)
copy_VarCoords(ua(0,:,:,:), uq)
copy_VarCoords(va(0,:,:,:), vq)
copy_VarCoords(pres(0,:,:,:), pres1)
;;--vertical integrated moisture flux
ptop = 100 ; top layer pressure
pbot = 1100
uqflux = vibeta(pres1(lat|:,lon|:,lev|:),uq(lat|:,lon|:,lev|:),2,psfc1,pbot,ptop)/9.8 ;should divide 9.8??
vqflux = vibeta(pres1(lat|:,lon|:,lev|:),vq(lat|:,lon|:,lev|:),2,psfc1,pbot,ptop)/9.8
; 1 means linear interpolation, 2 means log interpolation
qflux = sqrt(uqflux^2+vqflux^2)