2
$\begingroup$

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)
Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ esrl.noaa.gov/psd/data/composites/day/calculation.html You can make an attempt using the information present on that page. Show us where you are stuck and then we can try to interject at that point. $\endgroup$
    – gansub
    Dec 22 '18 at 3:38
  • $\begingroup$ how will you define a layer ? A layer is not a pressure surface. Usually it is some average of two surfaces numerically. That link I provided gives an idea on how they calculate layers $\endgroup$
    – gansub
    Dec 22 '18 at 4:16
  • $\begingroup$ Yes, you are right. Actually, I have searched the mailing list, and have not found any useful info to me. $\endgroup$
    – Love_qq_xq
    Dec 22 '18 at 4:17
  • $\begingroup$ I used the vibeta function in NCL to get the moisture flux of the entire layer from 1100 to 100. $\endgroup$
    – Love_qq_xq
    Dec 22 '18 at 4:19
  • $\begingroup$ no that is not a definition of a numerical atmospherical layer. A layer can be defined in many ways. Typically the surface in consideration and then above and below. You take an average of the values of those three surfaces. Then you integrate across the layers or do a Simpson's summation across layers to integrate or some other numerical integration scheme $\endgroup$
    – gansub
    Dec 22 '18 at 4:43
5
$\begingroup$

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 $$

Improve this answer
$\endgroup$
2
  • $\begingroup$ Dear @gansub, your answer is not clear for me. Would you be more clear please? For instance, is this correct?: dp = 1000 - 925, 925 - 850, etc qu = 100/g* ( q_1000*u_1000*(1000 - 925) + q_925*u_925*(925 - 850) + q_850*u_850*(850 - 700) + q_700*u_700*(700 - 500) ) $\endgroup$
    – Sergio
    Feb 24 '19 at 14:45
  • $\begingroup$ Thanks for the update. However, I calculated the VIMFT in this way: qu = (((u925*q925 + u1000*q1000)/2)*(1000 - 925)) + (((u850*q850 + u925*q925)/2)*(925 - 850)) +... And then the same for other layers and qv. I generated vector plots in agreement with the literature. However, the link that you mention for calculating the divergence is not working for me. For the diverence I should use the resulting qu and qv, right? Many thanks!!! $\endgroup$
    – Sergio
    Feb 27 '19 at 15:25

Not the answer you're looking for? Browse other questions tagged or ask your own question.