Merge "Issue #2043 switch capeFunc to java. Change-Id: I35db8657838d1e11180248ebfb5011181ea09659" into development
Former-commit-id:03b2db312e[formerly82341a9fc1] [formerly03b2db312e[formerly82341a9fc1] [formerly3be654e028[formerly aa178f96c90b1c13ade6eb7db297583e5241601c]]] Former-commit-id:3be654e028Former-commit-id:9ef5384cd3[formerlya765656f1d] Former-commit-id:93010b6cb2
This commit is contained in:
Richard Peter
Gerrit Code Review
commit
0751d58e02
4 changed files with 846 additions and 106 deletions
|
|
@ -18,11 +18,7 @@
|
|||
# further licensing information.
|
||||
##
|
||||
|
||||
# determine the location of the python ported meteolib
|
||||
# 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
|
||||
from com.raytheon.uf.viz.derivparam.python.function import CapeFunc
|
||||
|
||||
def execute(*args):
|
||||
return __execute(*args)[0]
|
||||
|
|
@ -93,104 +89,9 @@ def __execute(*args):
|
|||
pressureValue = pressure
|
||||
pressure = zeros(temperatureValues.shape, temperatureValues.dtype)
|
||||
pressure[:] = pressureValue
|
||||
|
||||
return capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure)
|
||||
|
||||
import numpy as np
|
||||
import ctypes as ct
|
||||
from com.raytheon.edex.meteoLib import MeteoLibUtil
|
||||
|
||||
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
|
||||
if upperTerminationPressure is None:
|
||||
threeDshape = pressureValues.shape
|
||||
return CapeFunc.capeFunc(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
|
||||
else:
|
||||
threeDshape = pressureValues.shape
|
||||
return CapeFunc.capeFuncTop(useVirtualTemp, pressureValues, temperatureValues, pressure, potentialTemperature, specificHumidity, upperTerminationPressure, threeDshape[1], threeDshape[2], threeDshape[0]).__numpy__
|
||||
|
|
|
|||
|
|
@ -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