<p>So i don't know for how many timesteps you've integrated. And I don't really understand what you have done to obtain the temperature <span class="math-container">$\rm T[t,s,y,x]$</span> in your second box, but this seems a bit like a redundant operation: the end result is <span class="math-container">$T:=T$</span>, so of course your right side looks like the left side, in terms of the temperature colourmap. I don't understand where the different scaling comes from though.</p>

<p>Your first numerical recipe is a whole different beast:
Approximating a differential in a first order symmetric fashion as <span class="math-container">$dT/dx \approx (T[t,y,x+1]-T[t,y,x-1])/(2*dx) $</span> is known to be numerically unstable. This is why when you want to stick to finite differences, one uses <a href="https://en.wikipedia.org/wiki/Upwind_scheme" rel="nofollow noreferrer">upwinding</a>.</p>

<p>Also your finite differences are cartesian ones, while you project them on a sphere, so to get physically correct results, you need to take the correct differentials in spherical coordinates into account.</p>

CC BY-SA 4.0
2019-03-20T19:32:10.680Z
I edited the question to remove the $s$ dimention in $T$. It is only a function of time and space $T$[t,y,x].
ouranos
CC BY-SA 4.0
2019-03-20T19:39:00.260Z
Three fields are given. T, u and v: temperature, zonal and meridional components of the horizontal wind. That is where I obtained T. I don't understand your comment _"...but this seems a bit like a redundant operation: the end result is :="_ Could you elaborate more? I am not calculating T. T is given. I am scaling u and v, multiplying each by T, making uT and vT. T is in module, because it contains negative values from anomalies. That is what I called "scaling wind components by temperature". The magnitude of the scaled wind vector is then given by $\sqrt{uT^{2}+vT^{2}}$.
ouranos
CC BY-SA 4.0
2019-03-20T19:43:55.983Z
Centered finite differences are very commonly used in this field GrADS has been successfuly using it for decades... http://cola.gmu.edu/grads/gadoc/gradfunchcurl.html You are probably right about the "upwinding". But that should not be the main source of the striking difference between the middle and rightmost figures. I am asking about my assumption that the classical temperature advection should be similar to the magnitude of the temperature-scaled wind.
ouranos
CC BY-SA 4.0
2019-03-20T19:44:44.250Z
The integration time should not matter, since this is not a prognostic quantity. The question is still valid for one time step.
ouranos
CC BY-SA 4.0
2019-03-20T19:46:38.200Z
The grid is rectangular, but $dy$ is changed as a function of latitude, as shown in the code. In addition, the given T, u and v fields have been weighted by the square root of the cosine of the latitude, prior to the calculations shown here.
ouranos
CC BY-SA 4.0
2019-03-22T16:03:30.517Z
@ouranos: I guess we speak a slightly different language. 'prognostic' means time-evolving for me, which you are doing when advecting the temperature. And I was pointing out, that the scheme you're using is unstable, and produces exponential growth of perturbations, which could manifest in like what you see in your center picture, if enough timesteps have been done. Furthermore I was pointing out that the right-hand method doesn't make sense to me, but that might not mean much.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-03-22T16:05:05.087Z
@ouranos: Other than that, if you're still stuck with that problem, I suggest to reduce the size of the problem: See what happens if you would advect only one temperature bump with your scheme along one direction. This would give you a cleaner setup and an easier way to identify what's wrong.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-03-25T17:20:12.367Z
I meant that there aren't any derivatives with respect to time. Thus, the number of time steps cannot influence. The calculation at a time step, does not depend on the previous or on the next time steps.
ouranos
CC BY-SA 4.0
2019-03-25T18:18:44.883Z
@ouranos: Well, you're solving the r.h.s of an advection equation in several steps, so yeah, you're solving a time-dependent problem. The flow field may not evolve, but you still can violate a CFL condition somewhere in your domain, when the advected flow overshoots a cell size.
AtmosphericPrisonEscape
CC BY-SA 4.0
2019-03-20T18:36:04.330Z
2019-03-20T18:36:04.330Z
AtmosphericPrisonEscape
<p>I've come to a solution. I have followed <a href="https://earthscience.stackexchange.com/users/1066/gansub">@gansub</a>'s suggestion to check what there is already there, before reinventing the wheel. I found MetPy's code very (too) well structured, so the code for advection was sparse among several different functions; not so immediate to trace back and put a function together. However, GrADS brings a very nice example on how to easily build such a function in their documentation page for the <code>cdiff</code> function <a href="http://cola.gmu.edu/grads/gadoc/gradfunccdiff.html" rel="nofollow noreferrer">here</a>. So I rewrote their code in Python, and it now works well. I wrote a small package available <a href="https://github.com/davidmnielsen/mygrads" rel="nofollow noreferrer">here</a>, if anyone is interested in the working solution. </p>

CC BY-SA 4.0
2019-12-17T19:02:07.923Z
2019-12-17T19:02:07.923Z
ouranos
<p>Advection of a scalar quantity, such as temperature (<em>T</em>), by the horizontal wind, is defined as follows:</p>

<p><span class="math-container">$-\textbf{U}\cdot\nabla T$</span></p>

<p>where <span class="math-container">$\textbf{U}$</span> is the horizontal wind vector, <span class="math-container">$\textbf{U}=(u,v)$</span>, being <span class="math-container">$u$</span> and <span class="math-container">$v$</span> its zonal and meridional components, respectively (<a href="http://glossary.ametsoc.org/wiki/Advection" rel="nofollow noreferrer">source</a>).</p>

<p>To compute it with gridded data, it is convenient to use finite differences (or similar numerical approaches), as pointed <a href="https://earthscience.stackexchange.com/questions/9634/how-to-compute-du-dx-and-dv-dy-in-moisture-flux-convergence/10667#10667">here</a>.</p>

<p>Using this method, the following lines written in Python calculate the advection of temperature (TAdv):</p>



<pre><code>dy=111000 # [m]
lonres=lon[1]-lon[0] # constant
tadv=np.zeros(np.shape(T)); tadv.fill(np.nan)
for t in range(np.shape(T)[0]):
    for x in np.arange(1,len(lon)-1):
        for y in np.arange(1,len(lat)-1):
            dx = abs(111000*np.cos(lat[y]*(2*np.pi/360))*lonres)
            tadv[t,y,x] = -(u[t,y,x]*(T[t,y,x+1]-T[t,y,x-1])/(2*dx) +\
                              v[t,y,x]*(T[t,y+1,x]-T[t,y-1,x])/(2*dy))
</code></pre>

<p>I would expect a similar result from scaling the wind components by the temperature, at each grid cell, which is done as follows:</p>



<pre><code>ut  = np.zeros(np.shape(T))
vt  = np.zeros(np.shape(T))
UT  = np.zeros(np.shape(T))
for t in range(np.shape(T)[0]):
    for x in range(len(lon)):
        for y in range(len(lat)):
            ut[t,y,x] = u[t,y,x]*abs(T[t,y,x])
            vt[t,y,x] = v[t,y,x]*abs(T[t,y,x])
            UT[t,y,x] = ((ut[t,y,x]**2 + vt[t,y,x]**2)**0.5)*np.sign(T[t,y,x])
</code></pre>

<p>The magnitude of the scaled wind vectors (UT) would be then similar to this "amount of transport of T". It is multiplied by the sign of T (-1 or 1), because T are anomalies (there is T&lt;0), in order to keep its the sign.</p>

<p>However, the results are totally different. What I am missing here? Is it a mistake in this assumption (that both results should look fairly similar)? Or is it a coding problem?</p>

<p>In the figure below, I show the original T, u and v fields (left), and the results of the two above calculations: Tadv with u and v (middle), and the magnitude of the scaled vectors, UT, with the scaled vectors themselves uT and vT (right). TAdv is multiplied by 24*3600, so that the units are K/day. </p>

<p><a href="https://i.sstatic.net/VxJpt.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/VxJpt.png" alt="figure"></a></p>

<p>Here is a link to a Dropbox folder, where the full script and sample data are available for download: <a href="https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0" rel="nofollow noreferrer">https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0</a></p>

CC BY-SA 4.0
2019-03-22T10:43:58.553Z
@gansub, the full code and sample data can be downloaded from this Dropbox folder: https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0 
Thank you for your help.
ouranos
CC BY-SA 4.0
2019-03-22T13:03:30.123Z
Sure. Done, link is edited into question.
ouranos
CC BY-SA 4.0
2019-03-25T17:15:25.957Z
The (0,0) position is at 0° lon and 20°N lat (0,20). This corner is not shown in this plot. Actually, all the boudaries are out of the figure domain.
ouranos
<meteorology><temperature><wind><numerical-modelling><atmospheric-circulation>
Temperature Advection using finite differences with gridded data
ouranos
# Temperature Advection using finite differences with gridded data

(Asked by: username_0 on 2019-03-20)

Advection of a scalar quantity, such as temperature (*T*), by the horizontal wind, is defined as follows:

$-\textbf{U}\cdot\nabla T$

where $\textbf{U}$ is the horizontal wind vector, $\textbf{U}=(u,v)$, being $u$ and $v$ its zonal and meridional components, respectively ([source](http://glossary.ametsoc.org/wiki/Advection)).

To compute it with gridded data, it is convenient to use finite differences (or similar numerical approaches), as pointed [here](https://earthscience.stackexchange.com/questions/9634/how-to-compute-du-dx-and-dv-dy-in-moisture-flux-convergence/10667#10667).

Using this method, the following lines written in Python calculate the advection of temperature (TAdv):

```
dy=111000 # [m]
lonres=lon[1]-lon[0] # constant
tadv=np.zeros(np.shape(T)); tadv.fill(np.nan)
for t in range(np.shape(T)[0]):
    for x in np.arange(1,len(lon)-1):
        for y in np.arange(1,len(lat)-1):
            dx = abs(111000*np.cos(lat[y]*(2*np.pi/360))*lonres)
            tadv[t,y,x] = -(u[t,y,x]*(T[t,y,x+1]-T[t,y,x-1])/(2*dx) +\
                              v[t,y,x]*(T[t,y+1,x]-T[t,y-1,x])/(2*dy))

```

I would expect a similar result from scaling the wind components by the temperature, at each grid cell, which is done as follows:

```
ut  = np.zeros(np.shape(T))
vt  = np.zeros(np.shape(T))
UT  = np.zeros(np.shape(T))
for t in range(np.shape(T)[0]):
    for x in range(len(lon)):
        for y in range(len(lat)):
            ut[t,y,x] = u[t,y,x]*abs(T[t,y,x])
            vt[t,y,x] = v[t,y,x]*abs(T[t,y,x])
            UT[t,y,x] = ((ut[t,y,x]**2 + vt[t,y,x]**2)**0.5)*np.sign(T[t,y,x])

```

The magnitude of the scaled wind vectors (UT) would be then similar to this "amount of transport of T". It is multiplied by the sign of T (-1 or 1), because T are anomalies (there is T<0), in order to keep its the sign.

However, the results are totally different. What I am missing here? Is it a mistake in this assumption (that both results should look fairly similar)? Or is it a coding problem?

In the figure below, I show the original T, u and v fields (left), and the results of the two above calculations: Tadv with u and v (middle), and the magnitude of the scaled vectors, UT, with the scaled vectors themselves uT and vT (right). TAdv is multiplied by 24\*3600, so that the units are K/day.

[![figure](https://i.sstatic.net/VxJpt.png)](https://i.sstatic.net/VxJpt.png)

Here is a link to a Dropbox folder, where the full script and sample data are available for download: <https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0>

**Tags:** meteorology, temperature, wind, numerical-modelling, atmospheric-circulation

> **Comments:**
> **username_0**: @gansub, the full code and sample data can be downloaded from this Dropbox folder: https://www.dropbox.com/sh/kcrb08h72jjj3rn/AAAIAUgKVrRrICAQtXhOHS2ta?dl=0  Thank you for your help.  
> **username_0**: Sure. Done, link is edited into question.  
> **username_0**: The (0,0) position is at 0° lon and 20°N lat (0,20). This corner is not shown in this plot. Actually, all the boudaries are out of the figure domain.  

---

## 2 answers

### Answer by username_1 on 2019-03-20. Score: 1

So i don't know for how many timesteps you've integrated. And I don't really understand what you have done to obtain the temperature $\rm T[t,s,y,x]$ in your second box, but this seems a bit like a redundant operation: the end result is $T:=T$, so of course your right side looks like the left side, in terms of the temperature colourmap. I don't understand where the different scaling comes from though.

Your first numerical recipe is a whole different beast:
Approximating a differential in a first order symmetric fashion as $dT/dx \approx (T[t,y,x+1]-T[t,y,x-1])/(2\*dx) $ is known to be numerically unstable. This is why when you want to stick to finite differences, one uses [upwinding](https://en.wikipedia.org/wiki/Upwind_scheme).

Also your finite differences are cartesian ones, while you project them on a sphere, so to get physically correct results, you need to take the correct differentials in spherical coordinates into account.

> **Comments:**
> **username_0**: I edited the question to remove the $s$ dimention in $T$. It is only a function of time and space $T$[t,y,x].  
> **username_0**: Three fields are given. T, u and v: temperature, zonal and meridional components of the horizontal wind. That is where I obtained T. I don't understand your comment _"...but this seems a bit like a redundant operation: the end result is :="_ Could you elaborate more? I am not calculating T. T is given. I am scaling u and v, multiplying each by T, making uT and vT. T is in module, because it contains negative values from anomalies. That is what I called "scaling wind components by temperature". The magnitude of the scaled wind vector is then given by $\sqrt{uT^{2}+vT^{2}}$.  
> **username_0**: Centered finite differences are very commonly used in this field GrADS has been successfuly using it for decades... http://cola.gmu.edu/grads/gadoc/gradfunchcurl.html You are probably right about the "upwinding". But that should not be the main source of the striking difference between the middle and rightmost figures. I am asking about my assumption that the classical temperature advection should be similar to the magnitude of the temperature-scaled wind.  
> **username_0**: The integration time should not matter, since this is not a prognostic quantity. The question is still valid for one time step.  
> **username_0**: The grid is rectangular, but $dy$ is changed as a function of latitude, as shown in the code. In addition, the given T, u and v fields have been weighted by the square root of the cosine of the latitude, prior to the calculations shown here.  
> **username_1**: @username_0: I guess we speak a slightly different language. 'prognostic' means time-evolving for me, which you are doing when advecting the temperature. And I was pointing out, that the scheme you're using is unstable, and produces exponential growth of perturbations, which could manifest in like what you see in your center picture, if enough timesteps have been done. Furthermore I was pointing out that the right-hand method doesn't make sense to me, but that might not mean much.  
> **username_1**: @username_0: Other than that, if you're still stuck with that problem, I suggest to reduce the size of the problem: See what happens if you would advect only one temperature bump with your scheme along one direction. This would give you a cleaner setup and an easier way to identify what's wrong.  
> **username_0**: I meant that there aren't any derivatives with respect to time. Thus, the number of time steps cannot influence. The calculation at a time step, does not depend on the previous or on the next time steps.  
> **username_1**: @username_0: Well, you're solving the r.h.s of an advection equation in several steps, so yeah, you're solving a time-dependent problem. The flow field may not evolve, but you still can violate a CFL condition somewhere in your domain, when the advected flow overshoots a cell size.  

---

### Answer by username_0 on 2019-12-17. Score: 1

I've come to a solution. I have followed [@gansub](https://earthscience.stackexchange.com/users/1066/gansub)'s suggestion to check what there is already there, before reinventing the wheel. I found MetPy's code very (too) well structured, so the code for advection was sparse among several different functions; not so immediate to trace back and put a function together. However, GrADS brings a very nice example on how to easily build such a function in their documentation page for the `cdiff` function [here](http://cola.gmu.edu/grads/gadoc/gradfunccdiff.html). So I rewrote their code in Python, and it now works well. I wrote a small package available [here](https://github.com/davidmnielsen/mygrads), if anyone is interested in the working solution.


---

