How to compute du/dx and dv/dy in moisture flux convergence?

I have $u$-wind, $v$-wind and specific humidity. I would like to compute moisture flux convergence at a grid point. So I need to compute the value of $q(\frac{du}{dx} + \frac{dv}{dy})$

My question is how to compute $\frac{du}{dx}$ and $\frac{dv}{dy}$ from $u$-wind and $v$-wind at gird point $(x,y)$?

• I did not do that. Moreover, that was with NCL. I would like to do it using a program like FORTRAN code or shell script. – Kay Feb 9 '17 at 14:07
• For this calculation you probably need the u-wind and v-wind of the neighboring grid cells as well as the distances to the neighboring grid cells. – FuzzyLeapfrog Feb 9 '17 at 16:21
• See item 3 under earthscience.stackexchange.com/questions/19/… – milancurcic Feb 9 '17 at 21:11
• @milancurcic : Thank you very much. That is exactly what I was looking for. – Kay Feb 10 '17 at 0:46
• @gansub : Thank you very much for your comment. I had searched a lot in google. But my doubt got cleared here by FuzzyLeapfrog and milancurcic. – Kay Feb 10 '17 at 0:58

I imagine you have a horizontal grid (x,y) and wind components u(x,y) and v(x,y).

Normally we do 1-degree of latitude constant (y-dimension) and equal to 110000 m. So your dy depends on how many degrees per grid point you have (yresolution). However, dx will vary according to the latitude.

dy=110000*yesolution;

Now, we may use centered finite differences to compute what you want. (Below is a Matlab code, but I believe generic enough to be reproducible in other languages).

for y=2:length(lat)-1 dx=abs(110000*cos(latx(y)*(2*pi/360))*xresolution); for x=2:length(lon)-1 div(x,y,1)= (u(x+1,y,1)-u(x-1,y,1))/(2*dx) + (v(x,y+1,1)-v(x,y-1,1)/(2*dy); end end

Note that you will have an empty frame around your divergence field when x=1, x=max(x), y=1 and y=max(y), once they do not have two neighboors to compute the differences. The same is observed when you do this with hdivg function in Grads, for example.

By the way, this is based on and yields the same result as Grads (I have checked!) using cdiff to reproduce hdivg.

Also note that the cosine function here works with radians, so if you compute cosines in degrees directly (i.e. cosd function in Matlab), you should omit the term scaling by (2*pi/360), which is just a conversion.

Hope it helps!

• centered finite differences will only work for non boundary grid points. At the grid point boundaries(such as the poles) you will need to use forward and backward differences. – gansub Jan 26 '19 at 2:32
• Checking this answer back a while later, you are correct @gansub about the boudaries. Computing divergence and vorticity with GrADS, using hdivg and hcurl, also return missing values at the boundary grid cells. I guess it's fine for most applications. – ouranos Mar 22 '19 at 10:55
• it is fine if you are using most regional grids. It is not OK if you are wanting to calculate the curl of the polar vortex (right over the poles). – gansub Mar 22 '19 at 12:11

The easy answer is to use finite difference. Which one is entirely your choice. For example, I'll choose a simple centered difference scheme.

You say you have a grid, and your point ($x$,$y$) is at the ith and jth gridpoint, where i indicates the left-right index on the grid and j indicates the north-south index on the grid. Therefore $u$ and $v$ can be expressed as $u(x_{i,j},y_{i,j})$ and $v(x_{i,j},y_{i,j})$.

$$\frac{du}{dx}\approx\frac{\Delta u}{\Delta x}=\frac{u(x+\Delta x,y)-u(x-\Delta x,y)}{\Delta x}=\frac{u(x_{i+1,j},y)-u(x_{i-1,j},y)}{x_{i+1,j}-x_{i-1,j}}=\frac{u_{i+1,j}-u_{i-1,j}}{x_{i+1,j}-x_{i-1,j}}$$

A similar process can be applied to the v component, just change the index, and the variables