[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[IDV #NVU-263991]: Error when calculating temperature tendency

> Hi,

The sigma is always position for sure, but the lap(-adv(TA)) does have the 
negative
values on the GFS dataset I tested over. I don't fully understand the physics 
of 
this formula, my suggestion is one of the two options:

1) use the abs() function over the results of lap().

2) use the substitute(data, low, high, newValue) to replace those negative 
values to newvalues.

Yuan
> 
> This request is a little urgent. I am about to finish my project and I am
> stuck on calculating one term of my equations. I am trying to calculate
> (R/p)/sigma * lap(-adv(T)), where sigma = (-R/p)*(T/theta)*(partial theta/
> partial p), which is the temperature advection term in the traditional form
> of the QG omega equation. I can compute the sigma factor separately and
> obtain positive values for each lat/lon. I can also compute the
> Laplacian of the temperature advection without any problem and multiply it
> by the R/p factor. However, when I try to put it together, i.e., calculate
> (R/p) * lap(-adv(TA)) over sigma, I get the sign of the Laplacian inverted
> in the central regions of the resulting field. If I change the computed
> sigma by a positive constant value with the same units, then I get
> apparently correct results for the quotient (no sign changes). I tried to
> save the fields in cache and then calculate the quotient from the cached
> fields but the problem still remains. I don't know how can a quotient over
> a positive term can change the sign of the dividend. I am attaching the
> script for the temperature advection term and the sigma factor, and the
> bundle with which I am working. I would like to know if there is any way to
> fix that.
> 
> Thank you very much.
> 
> Paula
> 
> # SIGMA FACTOR
> def sigma_V2(t, thta, thtaT, thtaB):
> from ucar.visad.quantities import DryAirGasConstant
> Rd = DryAirGasConstant.newReal()
> p = extractPressureFromNWPGrid(thta)
> dthtadp = GridMath.partial(thta,2)
> minus_Rd_over_p = - Rd / p
> T_over_thta = quo(t, thta)
> return mul(mul(minus_Rd_over_p, T_over_thta),dthtadp)
> 
> # TEMPERATURE ADVECTION TERM
> def TA_V2(z, t, thta, thtaT, thtaB):
> from ucar.visad.quantities import DryAirGasConstant
> Rd = DryAirGasConstant.newReal()
> p = extractPressureFromNWPGrid(t)
> Rd_over_p = Rd / p
> factor = quo(Rd_over_p, sigma_V2(t, thta, thtaT, thtaB))
> minus_TA = dot(geo(z), grad(t))
> return mul(factor, lap(minus_TA))
> 
> # - lap(TA)
> def minus_lapTA(z, t):
> return lap(dot(geo(z), grad(t)))
> 
> El mar, 31 may 2022 a las 23:30, Paula Caballero (<address@hidden>)
> escribió:
> 
> > Hi,
> >
> > I sent you another email explaining that I was trying to fix that very
> > same problem using that function (newUnit()) and that I was getting an
> > error message. This answer has already helpt me to fix that error; I have
> > realized I was not writing the " " in the last two arguments of the
> > function. Now it works; hence there is no need for you to reply to the
> > other email.
> >
> > Thank you very much!
> >
> > Paula
> >
> > El mar, 31 may 2022 a las 22:27, Unidata IDV Support (<
> > address@hidden>) escribió:
> >
> >> > Hi,
> >> >
> >> > I'm trying to calculate this:
> >> >
> >> > def TEST_LHomega(t, V, omega, theta, theta1, theta2):
> >> > T_over_theta = quo(t, theta)
> >> > dtheta = sub(theta1, theta2)
> >> > p1 = extractPressureFromNWPGrid(theta1)
> >> > p2 = extractPressureFromNWPGrid(theta2)
> >> > dp = sub(p1, p2)
> >> > dthetadp = quo(dtheta,dp)
> >> >
> >> > return mul(mul(T_over_theta, dthetadp),omega)
> >> >
> >> > where: t=temperature (K), theta = potential temperature (K), V=true wind
> >> > vector (m/s), omega=vertical pressure velocity (Pa/s).
> >> >
> >> > The final units should be K/s but I'm getting K*m/(kg*s^2), it seems
> >> that
> >> > the units from the omega field are not reported correctly. How can I
> >> fix it?
> >> >
> >> > Thank you,
> >> >
> >> > Paula
> >> >
> >>
> >> If only the unit of omega field incorrected, you can use the newUnit(..)
> >> to set the
> >> correct unit. Example: newUnit(grid, "topo", "m")
> >>
> >> Yuan
> >>
> >> Ticket Details
> >> ===================
> >> Ticket ID: NVU-263991
> >> Department: Support IDV
> >> Priority: High
> >> Status: Open
> >> ===================
> >> NOTE: All email exchanges with Unidata User Support are recorded in the
> >> Unidata inquiry tracking system and then made publicly available through
> >> the web.  If you do not want to have your interactions made available in
> >> this way, you must let us know in each email you send to us.
> >>
> >>
> >>
> 
> 


Ticket Details
===================
Ticket ID: NVU-263991
Department: Support IDV
Priority: Urgent
Status: Open
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.