Merge "Issue #2043 switch capeFunc to java. Change-Id: I35db8657838d1e11180248ebfb5011181ea09659" into development
Former-commit-id:03b2db312e[formerly82341a9fc1] [formerly3be654e028[formerly aa178f96c90b1c13ade6eb7db297583e5241601c]] Former-commit-id:3be654e028Former-commit-id:a765656f1d
This commit is contained in:
Richard Peter
Gerrit Code Review
commit
9ef5384cd3
4 changed files with 846 additions and 106 deletions
|
|
@ -18,11 +18,7 @@
|
||||||
# further licensing information.
|
# further licensing information.
|
||||||
##
|
##
|
||||||
|
|
||||||
# determine the location of the python ported meteolib
|
from com.raytheon.uf.viz.derivparam.python.function import CapeFunc
|
||||||
# TODO put the meteolib location in the java MasterDerivScript instead of having
|
|
||||||
# to manually determine its location like this
|
|
||||||
import sys
|
|
||||||
from numpy import zeros
|
|
||||||
|
|
||||||
def execute(*args):
|
def execute(*args):
|
||||||
return __execute(*args)[0]
|
return __execute(*args)[0]
|
||||||
|
|
@ -93,104 +89,9 @@ def __execute(*args):
|
||||||
pressureValue = pressure
|
pressureValue = pressure
|
||||||
pressure = zeros(temperatureValues.shape, temperatureValues.dtype)
|
pressure = zeros(temperatureValues.shape, temperatureValues.dtype)
|
||||||
pressure[:] = pressureValue
|
pressure[:] = pressureValue
|
||||||
|
if upperTerminationPressure is None:
|
||||||
return capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure)
|
threeDshape = pressureValues.shape
|
||||||
|
return CapeFunc.capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
|
||||||
import numpy as np
|
else:
|
||||||
import ctypes as ct
|
threeDshape = pressureValues.shape
|
||||||
from com.raytheon.edex.meteoLib import MeteoLibUtil
|
return CapeFunc.capeFuncTop(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
|
||||||
|
|
||||||
def capeFunc(usetv, p_dat, tve_dat, p0, th0, sh0, ptop = None):
|
|
||||||
""" Use the native capeFunc function
|
|
||||||
"""
|
|
||||||
|
|
||||||
# define the c_float_ptr type
|
|
||||||
c_float_ptr = ct.POINTER(ct.c_float)
|
|
||||||
|
|
||||||
# determine the input dimensions
|
|
||||||
dataShape = tve_dat.shape
|
|
||||||
|
|
||||||
nx = dataShape[2] if len(dataShape) > 2 else None
|
|
||||||
ny = dataShape[1] if len(dataShape) > 1 else None
|
|
||||||
nz = dataShape[0]
|
|
||||||
|
|
||||||
gridArea = nx * ny
|
|
||||||
|
|
||||||
# flatten all input arrays
|
|
||||||
p_dat = np.copy(p_dat)
|
|
||||||
tve_dat = np.copy(tve_dat)
|
|
||||||
p0 = np.copy(p0)
|
|
||||||
th0 = np.copy(th0)
|
|
||||||
sh0 = np.copy(sh0)
|
|
||||||
|
|
||||||
p_dat[np.isnan(p_dat)] = 1e37
|
|
||||||
tve_dat[np.isnan(tve_dat)] = 1e37
|
|
||||||
p0[np.isnan(p0)] = 1e37
|
|
||||||
th0[np.isnan(th0)] = 1e37
|
|
||||||
sh0[np.isnan(sh0)] = 1e37
|
|
||||||
|
|
||||||
p_dat.resize((nz, nx * ny,))
|
|
||||||
tve_dat.resize((nz, nx * ny,))
|
|
||||||
p0.resize((p0.size,))
|
|
||||||
th0.resize((th0.size,))
|
|
||||||
sh0.resize((sh0.size,))
|
|
||||||
|
|
||||||
if ptop != None:
|
|
||||||
ptop.resize((ptop.size,))
|
|
||||||
|
|
||||||
|
|
||||||
# load the library
|
|
||||||
meteoLibPath = MeteoLibUtil.getSoPath()
|
|
||||||
meteoLib = np.ctypeslib.load_library(meteoLibPath,"")
|
|
||||||
capeFunc = meteoLib.capeFunc if ptop == None else meteoLib.capeFuncTop
|
|
||||||
|
|
||||||
# "define" the capeFunc signature
|
|
||||||
capeFunc.restype = ct.c_int # return type
|
|
||||||
capeFunc.argtypes = [ct.c_float,
|
|
||||||
ct.POINTER(c_float_ptr),
|
|
||||||
ct.POINTER(c_float_ptr),
|
|
||||||
c_float_ptr,
|
|
||||||
c_float_ptr,
|
|
||||||
c_float_ptr,
|
|
||||||
ct.c_int,
|
|
||||||
ct.c_int,
|
|
||||||
ct.c_int,
|
|
||||||
ct.c_int,
|
|
||||||
c_float_ptr,
|
|
||||||
c_float_ptr]
|
|
||||||
|
|
||||||
if ptop != None:
|
|
||||||
capeFunc.argtypes.append(c_float_ptr)
|
|
||||||
|
|
||||||
# result arrays
|
|
||||||
cape_dat = np.zeros(gridArea,p_dat.dtype)
|
|
||||||
cin_dat = np.zeros(gridArea,p_dat.dtype)
|
|
||||||
|
|
||||||
capeFuncArgs = [ct.c_float(usetv),
|
|
||||||
|
|
||||||
# get c_style pointers to the 2D input arrays
|
|
||||||
(c_float_ptr*len(p_dat))(*[row.ctypes.data_as(c_float_ptr) for row in p_dat]),
|
|
||||||
(c_float_ptr*len(tve_dat))(*[row.ctypes.data_as(c_float_ptr) for row in tve_dat]),
|
|
||||||
|
|
||||||
p0.ctypes.data_as(c_float_ptr),
|
|
||||||
th0.ctypes.data_as(c_float_ptr),
|
|
||||||
sh0.ctypes.data_as(c_float_ptr),
|
|
||||||
ct.c_int(nx),
|
|
||||||
ct.c_int(nx),
|
|
||||||
ct.c_int(ny),
|
|
||||||
ct.c_int(nz),
|
|
||||||
cape_dat.ctypes.data_as(c_float_ptr),
|
|
||||||
cin_dat.ctypes.data_as(c_float_ptr)]
|
|
||||||
|
|
||||||
if ptop != None:
|
|
||||||
capeFuncArgs.append(ptop.ctypes.data_as(c_float_ptr))
|
|
||||||
|
|
||||||
val = capeFunc(*capeFuncArgs)
|
|
||||||
if val == 1 :
|
|
||||||
raise MemoryError('Cape derived parameter ran out of memory')
|
|
||||||
exit
|
|
||||||
# resize the cape and cin data to the appropriate grid size
|
|
||||||
cape_dat.resize((ny,nx))
|
|
||||||
cin_dat.resize((ny,nx))
|
|
||||||
|
|
||||||
return cape_dat, cin_dat
|
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,51 @@
|
||||||
|
/**
|
||||||
|
* This software was developed and / or modified by Raytheon Company,
|
||||||
|
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
|
||||||
|
*
|
||||||
|
* U.S. EXPORT CONTROLLED TECHNICAL DATA
|
||||||
|
* This software product contains export-restricted data whose
|
||||||
|
* export/transfer/disclosure is restricted by U.S. law. Dissemination
|
||||||
|
* to non-U.S. persons whether in the United States or abroad requires
|
||||||
|
* an export license or other authorization.
|
||||||
|
*
|
||||||
|
* Contractor Name: Raytheon Company
|
||||||
|
* Contractor Address: 6825 Pine Street, Suite 340
|
||||||
|
* Mail Stop B8
|
||||||
|
* Omaha, NE 68106
|
||||||
|
* 402.291.0100
|
||||||
|
*
|
||||||
|
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
|
||||||
|
* further licensing information.
|
||||||
|
**/
|
||||||
|
package com.raytheon.uf.viz.derivparam.python.function;
|
||||||
|
|
||||||
|
import static java.lang.Math.exp;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* This routine calculates the equivalent tempurature of a temperature and
|
||||||
|
* pressure using the adiabatic definition, assuming saturation put a fudge
|
||||||
|
* factor into L/cp to get agreement of moist adiabats with a published
|
||||||
|
* thermodynamic diagram *
|
||||||
|
*
|
||||||
|
* <pre>
|
||||||
|
*
|
||||||
|
* SOFTWARE HISTORY
|
||||||
|
*
|
||||||
|
* Date Ticket# Engineer Description
|
||||||
|
* ------------ ---------- ----------- --------------------------
|
||||||
|
* Jun 3, 2013 2043 bsteffen Ported from meteolib C
|
||||||
|
*
|
||||||
|
* </pre>
|
||||||
|
*
|
||||||
|
* @author bsteffen
|
||||||
|
* @version 1.0
|
||||||
|
*/
|
||||||
|
|
||||||
|
public class AdiabeticTemperature {
|
||||||
|
|
||||||
|
public static double adiabatic_te(double temp, double press) {
|
||||||
|
double e = exp(26.660820 - 0.0091379024 * temp - 6106.396 / temp);
|
||||||
|
e = 0.622 * e / (press - e);
|
||||||
|
return temp * exp(2740.0 * e / temp);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,627 @@
|
||||||
|
/**
|
||||||
|
* This software was developed and / or modified by Raytheon Company,
|
||||||
|
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
|
||||||
|
*
|
||||||
|
* U.S. EXPORT CONTROLLED TECHNICAL DATA
|
||||||
|
* This software product contains export-restricted data whose
|
||||||
|
* export/transfer/disclosure is restricted by U.S. law. Dissemination
|
||||||
|
* to non-U.S. persons whether in the United States or abroad requires
|
||||||
|
* an export license or other authorization.
|
||||||
|
*
|
||||||
|
* Contractor Name: Raytheon Company
|
||||||
|
* Contractor Address: 6825 Pine Street, Suite 340
|
||||||
|
* Mail Stop B8
|
||||||
|
* Omaha, NE 68106
|
||||||
|
* 402.291.0100
|
||||||
|
*
|
||||||
|
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
|
||||||
|
* further licensing information.
|
||||||
|
**/
|
||||||
|
package com.raytheon.uf.viz.derivparam.python.function;
|
||||||
|
|
||||||
|
import static com.raytheon.uf.viz.derivparam.python.function.AdiabeticTemperature.adiabatic_te;
|
||||||
|
import static com.raytheon.uf.viz.derivparam.python.function.TempOfTe.temp_of_te;
|
||||||
|
import static java.lang.Math.exp;
|
||||||
|
import static java.lang.Math.log;
|
||||||
|
import static java.lang.Math.pow;
|
||||||
|
import static java.lang.Math.sqrt;
|
||||||
|
import jep.INumpyable;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* We input theta and specific humidity for initial parcel because these are
|
||||||
|
* things that can be arithemitically averaged for a mixed layer. If usetv=1,
|
||||||
|
* buoyancy is done with virtual temp, if usetv=0 with temp. If usetv=0, must
|
||||||
|
* supply temps in tve_dat.
|
||||||
|
*
|
||||||
|
* There are a few small changes from the original C version. First all math is
|
||||||
|
* done using doubles instead of floats since that is what the java.lang.Math
|
||||||
|
* library uses and it gets more accurate results. Second all uses of 1e37 have
|
||||||
|
* been replaces with Double.NaN to be more consistant with the rest of derived
|
||||||
|
* parameters. Finally the loops have been reordered so that each grid cell is
|
||||||
|
* iterated in the outermost loop which allows all intermediary calculations to
|
||||||
|
* be stored as doubles instead of float[nx*ny] which dramatically reduces the
|
||||||
|
* memory footprint.
|
||||||
|
*
|
||||||
|
* <pre>
|
||||||
|
*
|
||||||
|
* SOFTWARE HISTORY
|
||||||
|
*
|
||||||
|
* Date Ticket# Engineer Description
|
||||||
|
* ------------ ---------- ----------- --------------------------
|
||||||
|
* Jun 3, 2013 2043 bsteffen Ported from meteolib C
|
||||||
|
*
|
||||||
|
* </pre>
|
||||||
|
*
|
||||||
|
* @author bsteffen
|
||||||
|
* @version 1.0
|
||||||
|
*/
|
||||||
|
|
||||||
|
public class CapeFunc {
|
||||||
|
|
||||||
|
private static final float c0 = 26.66082f;
|
||||||
|
|
||||||
|
private static final float c1 = 0.0091379024f;
|
||||||
|
|
||||||
|
private static final float c2 = 6106.396f;
|
||||||
|
|
||||||
|
private static final float c_1 = 223.1986f;
|
||||||
|
|
||||||
|
private static final float c_2 = 0.0182758048f;
|
||||||
|
|
||||||
|
private static final float kapa = 0.286f;
|
||||||
|
|
||||||
|
private static final float kapa_1 = 3.498257f;
|
||||||
|
|
||||||
|
public static CapeCinPair capeFunc(float usetv, float[] p_dat,
|
||||||
|
float[] tve_dat, float[] p0, float[] th0, float[] sh0, int nx,
|
||||||
|
int ny, int nz) {
|
||||||
|
int n2 = nx * ny;
|
||||||
|
|
||||||
|
float[] cin = new float[n2];
|
||||||
|
float[] cap = new float[n2];
|
||||||
|
|
||||||
|
// Calculate the parcel equivalent temp, virtual temp, and press at LCL.
|
||||||
|
// Make working copy of sfc press, use as press below current 3d
|
||||||
|
// pressure.
|
||||||
|
for (int i = 0; i < n2; i += 1) {
|
||||||
|
double tec = 0;
|
||||||
|
double tvc = 0;
|
||||||
|
double pc = 0;
|
||||||
|
double pp0 = p0[i];
|
||||||
|
double pp1 = pp0;
|
||||||
|
if (Double.isNaN(pp0) || Double.isNaN(th0[i])
|
||||||
|
|| Double.isNaN(sh0[i])
|
||||||
|
|| sh0[i] < 0.0005) {
|
||||||
|
tec = tvc = pc = Double.NaN;
|
||||||
|
} else {
|
||||||
|
double t0 = th0[i] * pow(pp0 / 1000, kapa);
|
||||||
|
double b = c0 - log(pp0 / (622.0 / sh0[i] + 0.378));
|
||||||
|
double td = (b - sqrt(b * b - c_1)) / c_2;
|
||||||
|
double tdc = td - (t0 - td)
|
||||||
|
* (-0.37329638 + 41.178204 / t0 + 0.0015945203 * td);
|
||||||
|
pc = pp0 * pow(tdc / t0, kapa_1);
|
||||||
|
tec = adiabatic_te(tdc, pc);
|
||||||
|
tvc = td * (1 + usetv * 0.000608 * sh0[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initialize md and pmd, which will be pressure of and max Te
|
||||||
|
// delta.
|
||||||
|
double md = 0;
|
||||||
|
double pmd = 0;
|
||||||
|
|
||||||
|
// Now calculate the virtual temperature of the parcel at the
|
||||||
|
// pressures in the input data. Then difference it from the
|
||||||
|
// environmental temp, which has been tweaked to not be cooler than
|
||||||
|
// dry adiabatic from the parcel start. Record the level of max
|
||||||
|
// parcel difference.
|
||||||
|
double[] tvp = new double[nz];
|
||||||
|
|
||||||
|
for (int k = 0; k < nz; k += 1) {
|
||||||
|
float pp = p_dat[k * n2 + i];
|
||||||
|
float tve = tve_dat[k * n2 + i];
|
||||||
|
if (Double.isNaN(pc) || Double.isNaN(pp) || Double.isNaN(tve)) {
|
||||||
|
tvp[k] = Double.NaN;
|
||||||
|
} else {
|
||||||
|
double t0 = tvc * pow(pp / pc, kapa);
|
||||||
|
if (pp > pc) {
|
||||||
|
tvp[k] = t0;
|
||||||
|
} else {
|
||||||
|
double td = tec * pow(pp / pc, kapa);
|
||||||
|
tvp[k] = td = temp_of_te(td, pp);
|
||||||
|
if (usetv > 0) {
|
||||||
|
tvp[k] *= pp
|
||||||
|
/ (pp - exp(25.687958917 - c1 * td - c2
|
||||||
|
/ td));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (tve < t0) {
|
||||||
|
tvp[k] -= t0;
|
||||||
|
} else {
|
||||||
|
tvp[k] -= tve;
|
||||||
|
}
|
||||||
|
if (pp > pc || tvp[k] < md) {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
md = tvp[k];
|
||||||
|
pmd = pp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// This loop performs the actual cape and cin calculation. Here we
|
||||||
|
// will reuse storage for virt temp, equiv temp, and max delta for
|
||||||
|
// prev parcel temp, neg and pos. neg and pos are pending negative
|
||||||
|
// and positive contributions we have not yet added into the cape
|
||||||
|
// and cin yet.
|
||||||
|
double neg = 0;
|
||||||
|
double pos = 0;
|
||||||
|
cin[i] = cap[i] = Float.NaN;
|
||||||
|
|
||||||
|
double tvp1 = tvp[0];
|
||||||
|
pp1 = p_dat[0];
|
||||||
|
for (int k = 1; k < nz; k += 1) {
|
||||||
|
float pp = p_dat[k * n2 + i];
|
||||||
|
if (Double.isNaN(pp0)) {
|
||||||
|
continue;
|
||||||
|
} else if (Double.isNaN(pp1) || Double.isNaN(tvp1)) {
|
||||||
|
;
|
||||||
|
} else if (pp >= pp1 || Double.isNaN(tvp[k])) {
|
||||||
|
continue;
|
||||||
|
} else if (pp >= pp0) {
|
||||||
|
;
|
||||||
|
} else {
|
||||||
|
// Now we finally have the data we need for calculating
|
||||||
|
// the cape/cin contribution for this layer.
|
||||||
|
if (Double.isNaN(cap[i])) {
|
||||||
|
cap[i] = cin[i] = 0;
|
||||||
|
}
|
||||||
|
if (pmd == 0) {
|
||||||
|
continue; // No parcel delta>0, we're done.
|
||||||
|
}
|
||||||
|
|
||||||
|
// First deal with possibility of bottom lvl being below the
|
||||||
|
// initial parcel.
|
||||||
|
double dlnp;
|
||||||
|
double dn;
|
||||||
|
if (pp1 > pp0) {
|
||||||
|
dlnp = log(pp0 / pp);
|
||||||
|
dn = 0;
|
||||||
|
} else {
|
||||||
|
dlnp = log(pp1 / pp);
|
||||||
|
dn = dlnp * 287 * tvp1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Now deal with the fact that not allowing superadiabatic
|
||||||
|
// layers means no cape below condensation pressure.
|
||||||
|
double up;
|
||||||
|
if (pp1 >= pc) {
|
||||||
|
if (dn > 0) {
|
||||||
|
dn = 0;
|
||||||
|
}
|
||||||
|
if (tvp[k] <= 0) {
|
||||||
|
up = dlnp * 287 * tvp[k];
|
||||||
|
} else if (pp >= pc) {
|
||||||
|
up = 0;
|
||||||
|
} else {
|
||||||
|
up = log(pc / pp) * 287 * tvp[k];
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
up = dlnp * 287 * tvp[k];
|
||||||
|
}
|
||||||
|
|
||||||
|
// Deal with where the break point is.
|
||||||
|
double b = up * dn >= 0 ? 0.5 : up / (up - dn);
|
||||||
|
up *= b;
|
||||||
|
dn *= (1 - b);
|
||||||
|
|
||||||
|
// Now consider this layer's contribution, taking into
|
||||||
|
// account transitions between positive and negative
|
||||||
|
// acceleration.
|
||||||
|
if (up == 0 && dn == 0) {
|
||||||
|
;
|
||||||
|
// Continuing deceleration.
|
||||||
|
} else if (up <= 0
|
||||||
|
&& (dn < 0 || dn == 0 && (pp < pmd || pos == 0))) {
|
||||||
|
neg -= up + dn;
|
||||||
|
|
||||||
|
// Continuing upward acceleration.
|
||||||
|
} else if (up >= 0
|
||||||
|
&& (dn > 0 || dn == 0 && (pp < pmd || neg == 0))) {
|
||||||
|
pos += up + dn;
|
||||||
|
if (pp > pmd && cap[i] + pos <= cin[i] + neg) {
|
||||||
|
; // no net cape and below max delta
|
||||||
|
} else if (pp > pmd || cap[i] == 0) {
|
||||||
|
// below max delta or cape uninitialized
|
||||||
|
cap[i] += pos;
|
||||||
|
cin[i] += neg;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else if (pos >= neg) {
|
||||||
|
// cape initialized and net positive contribution
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
} else if (up > 0 && dn <= 0) {
|
||||||
|
// Transition to upward acceleration.
|
||||||
|
neg += -dn;
|
||||||
|
if (pp1 <= pmd) {
|
||||||
|
// above max delta, only use net pos contribution
|
||||||
|
pos += up;
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
} else if (pp <= pmd) {
|
||||||
|
// straddle max delta, force cape initialization
|
||||||
|
if (cap[i] == 0) {
|
||||||
|
cin[i] += neg;
|
||||||
|
cap[i] += pos;
|
||||||
|
} else if (neg > pos) {
|
||||||
|
cin[i] += neg - pos;
|
||||||
|
} else {
|
||||||
|
cap[i] += pos - neg;
|
||||||
|
}
|
||||||
|
cap[i] += up;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else if (cap[i] + pos + up <= cin[i] + neg) {
|
||||||
|
// no net cape to this point
|
||||||
|
if (cap[i] + pos > 0) {
|
||||||
|
// reinitialize if there was cape before
|
||||||
|
cin[i] -= cap[i] + pos;
|
||||||
|
pos = cap[i] = 0;
|
||||||
|
}
|
||||||
|
cin[i] += neg;
|
||||||
|
pos += up;
|
||||||
|
neg = 0;
|
||||||
|
} else if (cap[i] == 0) { // initialize cape
|
||||||
|
cap[i] += pos + up;
|
||||||
|
cin[i] += neg;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else { // what remains, only use net pos contribution
|
||||||
|
pos += up;
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
// Transition to decceleration.
|
||||||
|
pos += dn;
|
||||||
|
if (pp1 <= pmd) {
|
||||||
|
// above max delta, only use net pos contribution
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
neg += -up;
|
||||||
|
} else if (cap[i] + pos <= cin[i] + neg - up) {
|
||||||
|
// no net cape to this point
|
||||||
|
if (cap[i] > 0) {
|
||||||
|
// reinitialize if there was cape before
|
||||||
|
cin[i] -= cap[i] + pos;
|
||||||
|
pos = cap[i] = 0;
|
||||||
|
}
|
||||||
|
cin[i] += neg - up;
|
||||||
|
pos = neg = 0;
|
||||||
|
} else if (cap[i] == 0) { // initialize cape
|
||||||
|
cap[i] += pos;
|
||||||
|
cin[i] += neg - up;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else { // what remains, only use net pos contribution
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
neg += -up;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// Make current layer top next layer bottom.
|
||||||
|
tvp1 = tvp[k];
|
||||||
|
pp1 = pp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return new CapeCinPair(ny, nx, cap, cin);
|
||||||
|
}
|
||||||
|
|
||||||
|
// In this version we stop the computation at some arbitrary upper level,
|
||||||
|
// ptop.
|
||||||
|
public static CapeCinPair capeFuncTop(float usetv, float[] p_dat,
|
||||||
|
float[] tve_dat, float[] p0, float[] th0, float[] sh0,
|
||||||
|
float[] ptop, int nx, int ny, int nz) {
|
||||||
|
int n2 = nx * ny;
|
||||||
|
|
||||||
|
float[] cin = new float[n2];
|
||||||
|
float[] cap = new float[n2];
|
||||||
|
|
||||||
|
// Calculate the parcel equivalent temp, virtual temp, and press at LCL.
|
||||||
|
// Make working copy of sfc press, use as press below current 3d
|
||||||
|
// pressure.
|
||||||
|
for (int i = 0; i < n2; i += 1) {
|
||||||
|
double tec = 0;
|
||||||
|
double tvc = 0;
|
||||||
|
double pc = 0;
|
||||||
|
double pp0 = p0[i];
|
||||||
|
double pp1 = pp0;
|
||||||
|
double pfin = ptop[i];
|
||||||
|
if (Double.isNaN(pp0) || Double.isNaN(th0[i])
|
||||||
|
|| Double.isNaN(sh0[i]) || sh0[i] < 0.0005
|
||||||
|
|| pp0 < pfin) {
|
||||||
|
tec = tvc = pc = Double.NaN;
|
||||||
|
} else {
|
||||||
|
double t0 = th0[i] * pow(pp0 / 1000, kapa);
|
||||||
|
double b = c0 - log(pp0 / (622.0 / sh0[i] + 0.378));
|
||||||
|
double td = (b - sqrt(b * b - c_1)) / c_2;
|
||||||
|
double tdc = td - (t0 - td)
|
||||||
|
* (-0.37329638 + 41.178204 / t0 + 0.0015945203 * td);
|
||||||
|
pc = pp0 * pow(tdc / t0, kapa_1);
|
||||||
|
tec = adiabatic_te(tdc, pc);
|
||||||
|
tvc = td * (1 + usetv * 0.000608 * sh0[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Initialize md and pmd, which will be pressure of and max Te
|
||||||
|
// delta.
|
||||||
|
double md = 0;
|
||||||
|
double pmd = 0;
|
||||||
|
|
||||||
|
// Now calculate the virtual temperature of the parcel at the
|
||||||
|
// pressures in the input data. Then difference it from the
|
||||||
|
// environmental temp, which has been tweaked to not be cooler than
|
||||||
|
// dry adiabatic from the parcel start. Record the level of max
|
||||||
|
// parcel difference.
|
||||||
|
double[] tvp = new double[nz];
|
||||||
|
|
||||||
|
for (int k = 0; k < nz; k += 1) {
|
||||||
|
float pp = p_dat[k * n2 + i];
|
||||||
|
float tve = tve_dat[k * n2 + i];
|
||||||
|
if (Double.isNaN(pc) || Double.isNaN(pp) || Double.isNaN(tve)
|
||||||
|
|| pp1 <= pfin) {
|
||||||
|
tvp[k] = Double.NaN;
|
||||||
|
} else {
|
||||||
|
pp1 = pp;
|
||||||
|
double t0 = tvc * pow(pp / pc, kapa);
|
||||||
|
if (pp > pc) {
|
||||||
|
tvp[k] = t0;
|
||||||
|
} else {
|
||||||
|
double td = tec * pow(pp / pc, kapa);
|
||||||
|
tvp[k] = td = temp_of_te(td, pp);
|
||||||
|
if (usetv > 0) {
|
||||||
|
tvp[k] *= pp
|
||||||
|
/ (pp - exp(25.687958917 - c1 * td - c2
|
||||||
|
/ td));
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (tve < t0) {
|
||||||
|
tvp[k] -= t0;
|
||||||
|
} else {
|
||||||
|
tvp[k] -= tve;
|
||||||
|
}
|
||||||
|
if (pp > pc || pp < pfin || tvp[k] < md) {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
md = tvp[k];
|
||||||
|
pmd = pp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// This loop performs the actual cape and cin calculation. Here we
|
||||||
|
// will reuse storage for virt temp, equiv temp, and max delta for
|
||||||
|
// prev parcel temp, neg and pos. neg and pos are pending negative
|
||||||
|
// and positive contributions we have not yet added into the cape
|
||||||
|
// and cin yet.
|
||||||
|
double neg = 0;
|
||||||
|
double pos = 0;
|
||||||
|
cin[i] = cap[i] = Float.NaN;
|
||||||
|
|
||||||
|
double tvp1 = tvp[0];
|
||||||
|
pp1 = p_dat[i];
|
||||||
|
for (int k = 1; k < nz; k += 1) {
|
||||||
|
float pp = p_dat[k * n2 + i];
|
||||||
|
if (Double.isNaN(pp0)) {
|
||||||
|
continue;
|
||||||
|
} else if (Double.isNaN(pp1) || Double.isNaN(tvp1)) {
|
||||||
|
;
|
||||||
|
} else if (pp >= pp1 || Double.isNaN(tvp[k])) {
|
||||||
|
continue;
|
||||||
|
} else if (pp >= pp0) {
|
||||||
|
;
|
||||||
|
} else {
|
||||||
|
// Now we finally have the data we need for calculating
|
||||||
|
// the cape/cin contribution for this layer.
|
||||||
|
if (Double.isNaN(cap[i])) {
|
||||||
|
cap[i] = cin[i] = 0;
|
||||||
|
}
|
||||||
|
if (pmd == 0) {
|
||||||
|
continue; // No parcel delta>0, we're done.
|
||||||
|
}
|
||||||
|
|
||||||
|
// First deal with possibility of bottom lvl being below the
|
||||||
|
// initial parcel and/or hitting the top of the computation.
|
||||||
|
double dlnp = 0;
|
||||||
|
double dn;
|
||||||
|
if (pp < pfin) {
|
||||||
|
double b = log(pp1 / pp);
|
||||||
|
dlnp = log(pp1 / pfin);
|
||||||
|
tvp[k] = tvp1 + (dlnp / b) * (tvp[k] - tvp1);
|
||||||
|
}
|
||||||
|
if (pp1 > pp0) {
|
||||||
|
if (pp < pfin) {
|
||||||
|
dlnp = log(pp0 / pfin);
|
||||||
|
} else {
|
||||||
|
dlnp = log(pp0 / pp);
|
||||||
|
}
|
||||||
|
dn = 0;
|
||||||
|
} else {
|
||||||
|
if (pp >= pfin) {
|
||||||
|
dlnp = log(pp1 / pp);
|
||||||
|
}
|
||||||
|
dn = dlnp * 287 * tvp1;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Now deal with the fact that not allowing superadiabatic
|
||||||
|
// layers means no cape below condensation pressure.
|
||||||
|
double up;
|
||||||
|
if (pp1 >= pc) {
|
||||||
|
if (dn > 0) {
|
||||||
|
dn = 0;
|
||||||
|
}
|
||||||
|
if (tvp[k] <= 0) {
|
||||||
|
up = dlnp * 287 * tvp[k];
|
||||||
|
} else if (pp >= pc) {
|
||||||
|
up = 0;
|
||||||
|
} else if (pp < pfin) {
|
||||||
|
up = log(pc / pfin) * 287 * tvp[k];
|
||||||
|
} else {
|
||||||
|
up = log(pc / pp) * 287 * tvp[k];
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
up = dlnp * 287 * tvp[k];
|
||||||
|
}
|
||||||
|
|
||||||
|
// Deal with where the break point is.
|
||||||
|
double b = up * dn >= 0 ? 0.5 : up / (up - dn);
|
||||||
|
up *= b;
|
||||||
|
dn *= (1 - b);
|
||||||
|
|
||||||
|
// Now consider this layer's contribution, taking into
|
||||||
|
// account transitions between positive and negative
|
||||||
|
// acceleration.
|
||||||
|
if (up == 0 && dn == 0) {
|
||||||
|
;
|
||||||
|
// Continuing deceleration.
|
||||||
|
} else if (up <= 0
|
||||||
|
&& (dn < 0 || dn == 0 && (pp < pmd || pos == 0))) {
|
||||||
|
neg -= up + dn;
|
||||||
|
|
||||||
|
// Continuing upward acceleration.
|
||||||
|
} else if (up >= 0
|
||||||
|
&& (dn > 0 || dn == 0 && (pp < pmd || neg == 0))) {
|
||||||
|
pos += up + dn;
|
||||||
|
if (pp > pmd && cap[i] + pos <= cin[i] + neg) {
|
||||||
|
; // no net cape and below max delta
|
||||||
|
} else if (pp > pmd || cap[i] == 0) {
|
||||||
|
// below max delta or cape uninitialized
|
||||||
|
cap[i] += pos;
|
||||||
|
cin[i] += neg;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else if (pos >= neg) {
|
||||||
|
// cape initialized and net positive contribution
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
} else if (up > 0 && dn <= 0) {
|
||||||
|
// Transition to upward acceleration.
|
||||||
|
neg += -dn;
|
||||||
|
if (pp1 <= pmd) {
|
||||||
|
// above max delta, only use net pos contribution
|
||||||
|
pos += up;
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
} else if (pp <= pmd) {
|
||||||
|
// straddle max delta, force cape initialization
|
||||||
|
if (cap[i] == 0) {
|
||||||
|
cin[i] += neg;
|
||||||
|
cap[i] += pos;
|
||||||
|
} else if (neg > pos) {
|
||||||
|
cin[i] += neg - pos;
|
||||||
|
} else {
|
||||||
|
cap[i] += pos - neg;
|
||||||
|
}
|
||||||
|
cap[i] += up;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else if (cap[i] + pos + up <= cin[i] + neg) {
|
||||||
|
// no net cape to this point
|
||||||
|
if (cap[i] + pos > 0) {
|
||||||
|
// reinitialize if there was cape before
|
||||||
|
cin[i] -= cap[i] + pos;
|
||||||
|
pos = cap[i] = 0;
|
||||||
|
}
|
||||||
|
cin[i] += neg;
|
||||||
|
pos += up;
|
||||||
|
neg = 0;
|
||||||
|
} else if (cap[i] == 0) { // initialize cape
|
||||||
|
cap[i] += pos + up;
|
||||||
|
cin[i] += neg;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else { // what remains, only use net pos contribution
|
||||||
|
pos += up;
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
// Transition to decceleration.
|
||||||
|
pos += dn;
|
||||||
|
if (pp1 <= pmd) {
|
||||||
|
// above max delta, only use net pos contribution
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
neg += -up;
|
||||||
|
} else if (cap[i] + pos <= cin[i] + neg - up) {
|
||||||
|
// no net cape to this point
|
||||||
|
if (cap[i] > 0) {
|
||||||
|
// reinitialize if there was cape before
|
||||||
|
cin[i] -= cap[i] + pos;
|
||||||
|
pos = cap[i] = 0;
|
||||||
|
}
|
||||||
|
cin[i] += neg - up;
|
||||||
|
pos = neg = 0;
|
||||||
|
} else if (cap[i] == 0) { // initialize cape
|
||||||
|
cap[i] += pos;
|
||||||
|
cin[i] += neg - up;
|
||||||
|
neg = pos = 0;
|
||||||
|
} else { // what remains, only use net pos contribution
|
||||||
|
if (pos >= neg) {
|
||||||
|
cap[i] += (pos - neg);
|
||||||
|
neg = pos = 0;
|
||||||
|
}
|
||||||
|
neg += -up;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// Make current layer top next layer bottom.
|
||||||
|
tvp1 = tvp[k];
|
||||||
|
pp1 = pp;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return new CapeCinPair(ny, nx, cap, cin);
|
||||||
|
}
|
||||||
|
|
||||||
|
public static class CapeCinPair implements INumpyable {
|
||||||
|
|
||||||
|
private final int nx;
|
||||||
|
|
||||||
|
private final int ny;
|
||||||
|
|
||||||
|
private final float[] cape;
|
||||||
|
|
||||||
|
private final float[] cin;
|
||||||
|
|
||||||
|
public CapeCinPair(int nx, int ny, float[] cape, float[] cin) {
|
||||||
|
this.nx = nx;
|
||||||
|
this.ny = ny;
|
||||||
|
this.cape = cape;
|
||||||
|
this.cin = cin;
|
||||||
|
}
|
||||||
|
|
||||||
|
@Override
|
||||||
|
public Object[] getNumPy() {
|
||||||
|
return new Object[] { cape, cin };
|
||||||
|
}
|
||||||
|
|
||||||
|
@Override
|
||||||
|
public int getNumpyX() {
|
||||||
|
return nx;
|
||||||
|
}
|
||||||
|
|
||||||
|
@Override
|
||||||
|
public int getNumpyY() {
|
||||||
|
return ny;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,161 @@
|
||||||
|
/**
|
||||||
|
* This software was developed and / or modified by Raytheon Company,
|
||||||
|
* pursuant to Contract DG133W-05-CQ-1067 with the US Government.
|
||||||
|
*
|
||||||
|
* U.S. EXPORT CONTROLLED TECHNICAL DATA
|
||||||
|
* This software product contains export-restricted data whose
|
||||||
|
* export/transfer/disclosure is restricted by U.S. law. Dissemination
|
||||||
|
* to non-U.S. persons whether in the United States or abroad requires
|
||||||
|
* an export license or other authorization.
|
||||||
|
*
|
||||||
|
* Contractor Name: Raytheon Company
|
||||||
|
* Contractor Address: 6825 Pine Street, Suite 340
|
||||||
|
* Mail Stop B8
|
||||||
|
* Omaha, NE 68106
|
||||||
|
* 402.291.0100
|
||||||
|
*
|
||||||
|
* See the AWIPS II Master Rights File ("Master Rights File.pdf") for
|
||||||
|
* further licensing information.
|
||||||
|
**/
|
||||||
|
package com.raytheon.uf.viz.derivparam.python.function;
|
||||||
|
|
||||||
|
import static com.raytheon.uf.viz.derivparam.python.function.AdiabeticTemperature.adiabatic_te;
|
||||||
|
import static java.lang.Math.sqrt;
|
||||||
|
|
||||||
|
/**
|
||||||
|
* This routine calculates the saturation tempurature of an equivalent
|
||||||
|
* temperature at given pressure using the adiabatic definition
|
||||||
|
*
|
||||||
|
* <pre>
|
||||||
|
*
|
||||||
|
* SOFTWARE HISTORY
|
||||||
|
*
|
||||||
|
* Date Ticket# Engineer Description
|
||||||
|
* ------------ ---------- ----------- --------------------------
|
||||||
|
* Jun 3, 2013 2043 bsteffen Ported from meteolib C
|
||||||
|
*
|
||||||
|
* </pre>
|
||||||
|
*
|
||||||
|
* @author bsteffen
|
||||||
|
* @version 1.0
|
||||||
|
*/
|
||||||
|
|
||||||
|
public class TempOfTe {
|
||||||
|
|
||||||
|
private static final int tmin = 193;
|
||||||
|
|
||||||
|
private static final int tmax = 333;
|
||||||
|
|
||||||
|
private static final int nval = 1 + tmax - tmin;
|
||||||
|
|
||||||
|
private static final double[] Te1000 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te850 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te700 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te600 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te500 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te350 = new double[nval];
|
||||||
|
|
||||||
|
private static final double[] Te200 = new double[nval];
|
||||||
|
static {
|
||||||
|
int i = 0;
|
||||||
|
for (int t = tmin; t <= tmax; t += 1) {
|
||||||
|
Te1000[i] = adiabatic_te(t, 1000);
|
||||||
|
Te850[i] = adiabatic_te(t, 850);
|
||||||
|
Te700[i] = adiabatic_te(t, 700);
|
||||||
|
Te600[i] = adiabatic_te(t, 600);
|
||||||
|
Te500[i] = adiabatic_te(t, 500);
|
||||||
|
Te350[i] = adiabatic_te(t, 350);
|
||||||
|
Te200[i] = adiabatic_te(t, 200);
|
||||||
|
i += 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
public static double temp_of_te(double te, double press) {
|
||||||
|
double[] TeLookup = null;
|
||||||
|
double base;
|
||||||
|
/* find correct table, check for beyond bounds of table */
|
||||||
|
if (press <= 250) {
|
||||||
|
TeLookup = Te200;
|
||||||
|
base = 200;
|
||||||
|
} else if (press <= 400) {
|
||||||
|
TeLookup = Te350;
|
||||||
|
base = 350;
|
||||||
|
} else if (press <= 550) {
|
||||||
|
TeLookup = Te500;
|
||||||
|
base = 500;
|
||||||
|
} else if (press <= 650) {
|
||||||
|
TeLookup = Te600;
|
||||||
|
base = 600;
|
||||||
|
} else if (press <= 750) {
|
||||||
|
TeLookup = Te700;
|
||||||
|
base = 700;
|
||||||
|
} else if (press <= 900) {
|
||||||
|
TeLookup = Te850;
|
||||||
|
base = 850;
|
||||||
|
} else {
|
||||||
|
TeLookup = Te1000;
|
||||||
|
base = 1000;
|
||||||
|
}
|
||||||
|
if (te < TeLookup[1]) {
|
||||||
|
return te;
|
||||||
|
}
|
||||||
|
if (te >= TeLookup[nval - 1]) {
|
||||||
|
return Double.NaN;
|
||||||
|
}
|
||||||
|
/* use table to get first guesses for value of temp */
|
||||||
|
double t1 = tmin;
|
||||||
|
double t2 = (int) te;
|
||||||
|
if (t2 > tmax) {
|
||||||
|
t2 = tmax;
|
||||||
|
}
|
||||||
|
double t;
|
||||||
|
while (t2 - t1 >= 3) {
|
||||||
|
t = (int) ((t1 + t2) / 2);
|
||||||
|
if (TeLookup[(int) t - tmin] > te) {
|
||||||
|
t2 = t;
|
||||||
|
} else if (TeLookup[(int) t - tmin] < te) {
|
||||||
|
t1 = t;
|
||||||
|
} else {
|
||||||
|
if (t1 < t - 1) {
|
||||||
|
t1 = t - 1;
|
||||||
|
}
|
||||||
|
if (t2 > t + 1) {
|
||||||
|
t2 = t + 1;
|
||||||
|
}
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
double w = sqrt(base / press);
|
||||||
|
t1 = (1 - w) * TeLookup[(int) t1 - tmin] + w * t1;
|
||||||
|
t2 = (1 - w) * TeLookup[(int) t2 - tmin] + w * t2;
|
||||||
|
|
||||||
|
/* Iterate to find the exact solution */
|
||||||
|
double d1 = te - adiabatic_te(t1, press);
|
||||||
|
double d2 = adiabatic_te(t2, press) - te;
|
||||||
|
w = d2 / (d1 + d2);
|
||||||
|
t = w * t1 + (1 - w) * t2;
|
||||||
|
double d = adiabatic_te(t, press) - te;
|
||||||
|
int i = 0;
|
||||||
|
while (i++ < 10) {
|
||||||
|
if (d > 0.01) {
|
||||||
|
d2 = d;
|
||||||
|
t2 = t;
|
||||||
|
} else if (d < -0.01) {
|
||||||
|
d1 = -d;
|
||||||
|
t1 = t;
|
||||||
|
} else {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
w = d2 / (d1 + d2);
|
||||||
|
t = w * t1 + (1 - w) * t2;
|
||||||
|
d = adiabatic_te(t, press) - te;
|
||||||
|
}
|
||||||
|
return t;
|
||||||
|
}
|
||||||
|
}
|
||||||
Loading…
Add table
Reference in a new issue