WRF vertical integrated moisture flux using NCL

Question

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)    

Answer 1

I am going to explain my answer based on how the authors from this paper calculate this Vertically integrated moisture flux convergence as a predictor of thunderstorms. This is the equation that is defined in their paper

$$ VIMFC = -\frac{1}{g} \int_{700\,hPa}^{1000\,hPa}\ (\frac{\partial u q}{\partial x}+\frac{\partial v q}{\partial y}) * dp $$

For a reference for the vertically integrated moisture flux please look at the equation in the reference A Hydrological Definition of Indian Monsoon Onset and Withdrawal

$$ VIMFT = -\frac{1}{g} \int_{700\,hPa}^{1000\,hPa} (q * u + q * v) * dp$$

So I will provide an algorithm for a numerical implementation of this equation and opinions can vary what exactly constitutes a "atmospheric layer" in NWP.

In this equation q is the specific humidity, u and v are the zonal and meridional velocities, p is the pressure, g is the acceleration due to gravity.

Moisture flux convergence or divergence is

$$ (\frac{\partial u q}{\partial x}+\frac{\partial v q}{\partial y})$$

Moisture flux is

$$ q * u + q * v $$

You then need to sum the horizontal moisture flux convergence(or in your case moisture flux) over the following layers 1) 1000 - 925 hPa 2) 925-850 hPa 3) 850-700 hPa So the definition of what a "layer" is subjective and you can include three pressure surfaces in it if you think a deep layer is required as opposed to a shallow layer. Numerically a layer can be thought as the average values of two pressure surfaces. As an example in your case you take the average of the velocities at the 1000 hPa and 925 hPa surfaces and the average of the specific humidity at the 1000 hPa and 925 hPa.

Once you calculate the horizontal moisture flux(or if you are calculating MFC then you need to use finite centered differences and backward and forward differences at the grid boundaries and that procedure is explained in this ESSE answer - Moisture flux convergence numerical procedure and here Finite differences for calculating MFC) over each layer you multiply that by dp and then sum up the contributions.

UPDATE

In response to the questions in the comments

$"1000 - 925"$

is not to be interpreted as a subtraction but as an average.

What needs to be done is to take an average of two pressure surfaces in the following manner

$ q_{layer} = (q_{1000} + q_{925})/2$

$ u_{layer} = (u_{1000} + u_{925})/2$

$ v_{layer} = (v_{1000} + v_{925})/2$

Once you do this for every layer multiply by the corresponding value of dp and in this case

$ dp = 1000 -925 $

This is a subtraction. Then you plug the individual values into the sum of individual contributions for each layer and then sum the layers

$ (q_{layer}* u_{layer} + q_{layer}* v_{layer}) * dp$

The units of VIMFC is typically

$$ 10^− 5 \, kg \, m^− 2 \,s^− 1 $$