Issue #2043 Port Dcape to java.
Change-Id: Ibcdd38fc47bb57eb6b2c6d0921d9790c2f1893ba Former-commit-id:6be12509e5[formerly6be12509e5[formerly eddf7c71bb23057fe4e68a7462a634c6d81799d9]] Former-commit-id:660b04ba32Former-commit-id:98dfa3f8b3
This commit is contained in:
Ben Steffensmeier
parent
abfd0b9ead
commit
4404e55ab8
6 changed files with 551 additions and 105 deletions
|
|
@ -22,13 +22,12 @@
|
|||
#
|
||||
# Date Ticket# Engineer Description
|
||||
# ------------ ---------- ----------- --------------------------
|
||||
# 17 Feb 2010 #4502 jelkins Initial Creation.
|
||||
# Feb 17, 2010 4502 jelkins Initial Creation.
|
||||
# Jun 05, 2013 2043 bsteffen Ported from meteolib C
|
||||
|
||||
# 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 DCapeFunc
|
||||
|
||||
|
||||
def execute(threeDtemperature, threeDdewpoint, pressure, potentialTemperature, specificHumidity,maxEvaporation,maxRelativeHumidity,useVirtualTemp):
|
||||
""" Calculate Downdraft Convective Available Potential Energy
|
||||
|
|
@ -62,89 +61,6 @@ def execute(threeDtemperature, threeDdewpoint, pressure, potentialTemperature, s
|
|||
pressureValue = pressure
|
||||
pressure = zeros(potentialTemperature.shape, potentialTemperature.dtype)
|
||||
pressure[:] = pressureValue
|
||||
threeDshape = threeDpressure.shape
|
||||
return DCapeFunc.dcapeFunc(useVirtualTemp, threeDpressure, threeDtemperature[0], threeDdewpoint[0], pressure, potentialTemperature, specificHumidity, threeDshape[1], threeDshape[2], threeDshape[0], maxEvaporation, maxRelativeHumidity).__numpy__[0]
|
||||
|
||||
return dcapeFunc(useVirtualTemp, threeDpressure, threeDtemperature[0], threeDdewpoint[0], pressure, potentialTemperature, specificHumidity, maxEvaporation, maxRelativeHumidity)
|
||||
|
||||
from com.raytheon.edex.meteoLib import MeteoLibUtil
|
||||
import numpy as np
|
||||
import ctypes as ct
|
||||
|
||||
def dcapeFunc(usetv, p_dat, t_dat, td_dat, p0, th0, sh0, max_evap, max_rh):
|
||||
""" Use the native dcapeFunc function
|
||||
"""
|
||||
|
||||
# define the c_float_ptr type
|
||||
c_float_ptr = ct.POINTER(ct.c_float)
|
||||
|
||||
# determine the input dimensions
|
||||
dataShape = t_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)
|
||||
t_dat = np.copy(t_dat)
|
||||
td_dat = np.copy(td_dat)
|
||||
p0 = np.copy(p0)
|
||||
th0 = np.copy(th0)
|
||||
sh0 = np.copy(sh0)
|
||||
p_dat.resize((nz, nx * ny,))
|
||||
t_dat.resize((nz, nx * ny,))
|
||||
td_dat.resize((nz, nx * ny,))
|
||||
p0.resize((p0.size,))
|
||||
th0.resize((th0.size,))
|
||||
sh0.resize((sh0.size,))
|
||||
|
||||
# load the library
|
||||
meteoLibPath = MeteoLibUtil.getSoPath()
|
||||
meteoLib = np.ctypeslib.load_library(meteoLibPath,"")
|
||||
dcapeFunc = meteoLib.dcapeFunc
|
||||
|
||||
# "define" the capeFunc signature
|
||||
dcapeFunc.restype = None # return type
|
||||
dcapeFunc.argtypes = [ct.c_float,
|
||||
ct.POINTER(c_float_ptr),
|
||||
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,
|
||||
ct.c_float,
|
||||
ct.c_float,
|
||||
c_float_ptr]
|
||||
|
||||
# result arrays
|
||||
dcape_dat = np.zeros(gridArea,p_dat.dtype)
|
||||
|
||||
dcapeFuncArgs = [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(t_dat))(*[row.ctypes.data_as(c_float_ptr) for row in t_dat]),
|
||||
(c_float_ptr*len(td_dat))(*[row.ctypes.data_as(c_float_ptr) for row in td_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),
|
||||
ct.c_float(max_evap),
|
||||
ct.c_float(max_rh),
|
||||
dcape_dat.ctypes.data_as(c_float_ptr)]
|
||||
|
||||
dcapeFunc(*dcapeFuncArgs)
|
||||
|
||||
# resize the cape data to the appropriate grid size
|
||||
dcape_dat.resize((ny,nx))
|
||||
|
||||
return dcape_dat
|
||||
|
|
|
|||
|
|
@ -19,6 +19,9 @@
|
|||
**/
|
||||
package com.raytheon.uf.viz.derivparam.python.function;
|
||||
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c0;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c2;
|
||||
import static java.lang.Math.exp;
|
||||
|
||||
/**
|
||||
|
|
@ -44,7 +47,7 @@ import static java.lang.Math.exp;
|
|||
public class AdiabeticTemperature {
|
||||
|
||||
public static double adiabatic_te(double temp, double press) {
|
||||
double e = exp(26.660820 - 0.0091379024 * temp - 6106.396 / temp);
|
||||
double e = exp(c0 - c1 * temp - c2 / temp);
|
||||
e = 0.622 * e / (press - e);
|
||||
return temp * exp(2740.0 * e / temp);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -0,0 +1,136 @@
|
|||
/**
|
||||
* 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.Constants.c0;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.f;
|
||||
import static java.lang.Math.abs;
|
||||
import static java.lang.Math.exp;
|
||||
import static java.lang.Math.log;
|
||||
import static java.lang.Math.max;
|
||||
import static java.lang.Math.min;
|
||||
import static java.lang.Math.sqrt;
|
||||
|
||||
/**
|
||||
* Routine to calculate wetbulb from temperature, and relative humidity.
|
||||
*
|
||||
* <pre>
|
||||
* Inputs/Outputs:
|
||||
*
|
||||
* Variable Var Type I/O Description
|
||||
* ---------- ---------- ----- -------------
|
||||
* p RA I Pressure (mb)
|
||||
* t RA I Temperature (K)
|
||||
* rh RA I Relative humidity [range: 0. - 100.]
|
||||
* tw RA O wet-bulb temp (K)
|
||||
*
|
||||
*
|
||||
* User Notes:
|
||||
*
|
||||
* 1. No quality control is performed in this routine.
|
||||
* </pre>
|
||||
*
|
||||
* <pre>
|
||||
*
|
||||
* SOFTWARE HISTORY
|
||||
*
|
||||
* Date Ticket# Engineer Description
|
||||
* ------------ ---------- ----------- --------------------------
|
||||
* Jun 06, 2013 2043 bsteffen Ported from meteolib fortran
|
||||
*
|
||||
* </pre>
|
||||
*
|
||||
* @author bsteffen
|
||||
* @version 1.0
|
||||
*/
|
||||
|
||||
public class CalcTw {
|
||||
|
||||
public static double calctw(double p, double t, double rh) {
|
||||
double rhqc = min(100.0, max(1.0, rh));
|
||||
double b = c1 * t + c2 / t - log(rhqc / 100.0);
|
||||
double td = (b - sqrt(b * b - c_1)) / c_2;
|
||||
return mytw(t, td, p);
|
||||
}
|
||||
|
||||
/*
|
||||
* This function takes temperature in degrees K, dewpoint in degrees K and
|
||||
* pressure in millibars and returns the isobaric wet-bulb temperature in
|
||||
* degrees K using an iterative technique. For a given guess for the wet
|
||||
* bulb temp, one tries to do an energy balance, matching cp*(T-Tw) to
|
||||
* (esat(Tw)-esat(Td))*eps*L/p*.
|
||||
*/
|
||||
public static double mytw(double k, double kd, double p) {
|
||||
// Special cases of Td >= T or a ridiculously low T.
|
||||
if (kd >= k) {
|
||||
return (k + kd) / 2;
|
||||
} else if (k < 100) {
|
||||
return k;
|
||||
}
|
||||
|
||||
// Special case of a ridiculously high saturation vapor pressure.
|
||||
double ew = c0 - c1 * k - c2 / k;
|
||||
if (ew > 10.0) {
|
||||
return (k + kd) / 2;
|
||||
}
|
||||
ew = exp(ew);
|
||||
|
||||
// Kw is our current guess for wet-bulb, ed the vapor pressure
|
||||
// corresponding to the depoint. Deal with case of a ridiculously small
|
||||
// dewpoint vapor pressure.
|
||||
double kdx = kd;
|
||||
double ed = c0 - c1 * kdx - c2 / kdx;
|
||||
while (ed < -50.0) {
|
||||
kdx = kdx + 10;
|
||||
ed = c0 - c1 * kdx - c2 / kdx;
|
||||
}
|
||||
ed = exp(ed);
|
||||
double fp = p * f;
|
||||
double s = (ew - ed) / (k - kdx);
|
||||
double kw = (k * fp + kdx * s) / (fp + s);
|
||||
|
||||
// At each step of the iteration, esat(Tw)-esat(Td) is compared to
|
||||
// (T-Tw)*p/(eps*L). When that difference is less than one part in
|
||||
// 10000 of esat(Tw), or ten iterations have been done, the iteration
|
||||
// stops.
|
||||
// This is basically trying to find the value of Kw where de is 0. The
|
||||
// value s is the derivative of de with respect to Kw, a fairly standard
|
||||
// numerical technique for finding the zero value of a function.
|
||||
for (int l = 1; l <= 10; l += 1) {
|
||||
ew = c0 - c1 * kw - c2 / kw;
|
||||
if (ew < -50.0 || ew > 10.0) {
|
||||
break;
|
||||
}
|
||||
ew = exp(ew);
|
||||
double de = fp * (k - kw) + ed - ew;
|
||||
if (abs(de / ew) < 1e-5) {
|
||||
continue;
|
||||
}
|
||||
s = ew * (c1 - c2 / (kw * kw)) - fp;
|
||||
kw = kw - de / s;
|
||||
}
|
||||
return kw;
|
||||
|
||||
}
|
||||
}
|
||||
|
|
@ -20,6 +20,13 @@
|
|||
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.Constants.c0;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.kapa;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.kapa_1;
|
||||
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;
|
||||
|
|
@ -58,20 +65,6 @@ import jep.INumpyable;
|
|||
|
||||
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) {
|
||||
|
|
|
|||
|
|
@ -0,0 +1,64 @@
|
|||
/**
|
||||
* 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;
|
||||
|
||||
/**
|
||||
* Consolidated constants from various meteolib functions.
|
||||
*
|
||||
* <pre>
|
||||
*
|
||||
* SOFTWARE HISTORY
|
||||
*
|
||||
* Date Ticket# Engineer Description
|
||||
* ------------ ---------- ----------- --------------------------
|
||||
* Jun 06, 2013 2043 bsteffen Ported from meteolib C
|
||||
*
|
||||
* </pre>
|
||||
*
|
||||
* @author bsteffen
|
||||
* @version 1.0
|
||||
*/
|
||||
|
||||
public class Constants {
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double c0 = 26.66082;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double c1 = 0.0091379024;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double c2 = 6106.396;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double c_1 = 223.1986;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double c_2 = 0.0182758048;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double kapa = 0.286;
|
||||
|
||||
// from meteoLib capeFunc.c
|
||||
public static final double kapa_1 = 3.498257;
|
||||
|
||||
// from meteoLib calctw.f
|
||||
public static final double f = 0.0006355;
|
||||
}
|
||||
|
|
@ -0,0 +1,334 @@
|
|||
/**
|
||||
* 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.CalcTw.mytw;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c0;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_1;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.c_2;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.kapa;
|
||||
import static com.raytheon.uf.viz.derivparam.python.function.Constants.kapa_1;
|
||||
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 com.raytheon.uf.common.python.PythonNumpyFloatArray;
|
||||
|
||||
/**
|
||||
* We input theta and specific humidity for surface conditions because these are
|
||||
* things that can be arithemitically averaged for a mixed layer. In order for a
|
||||
* dcape to be calculated, there must be positive bouyancy somewhere in the
|
||||
* column based on the input surface conditions, neglecting any cap. Sinking
|
||||
* parcel starts from the minimum thetaE level that is above the condensation
|
||||
* pressure and below 400mb and the highest positive rising parcel bouyancy.
|
||||
* max_evap is the limit to how much water can be evaporated into the parcel as
|
||||
* it decends, in Kg/Kg. max_rh is the desired RH (%) as the parcel reaches the
|
||||
* surface. Will initially evaporate up to one-third of max_evap into the
|
||||
* sinking parcel at the start, afterward attempting to trend the RH toward
|
||||
* max_rh at the ground. If usetv=1, buoyancy is done with virtual temp, if
|
||||
* usetv=0 with temp.
|
||||
*
|
||||
* 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 05, 2013 2043 bsteffen Ported from meteolib C
|
||||
*
|
||||
* </pre>
|
||||
*
|
||||
* @author bsteffen
|
||||
* @version 1.0
|
||||
*/
|
||||
|
||||
public class DCapeFunc {
|
||||
|
||||
public static PythonNumpyFloatArray dcapeFunc(float usetv, float[] p_dat,
|
||||
float[] t_dat, float[] td_dat, float[] p0, float[] th0,
|
||||
float[] sh0, int nx, int ny, int nz, float max_evap, float max_rh) {
|
||||
int n2 = nx * ny;
|
||||
int nzm = nz - 1;
|
||||
|
||||
float[] dcape = new float[n2];
|
||||
|
||||
for (int i = 0; i < n2; i += 1) {
|
||||
|
||||
// Calculate the ascending parcel start equivalent temp, virtual
|
||||
// temp, and press at LCL, and the initial virtual temperature.
|
||||
// Initialize pm and wm, which now will be pressure at and min
|
||||
// environmetal virtual theta E.
|
||||
double tec = 0;
|
||||
double tvc = 0;
|
||||
double pc = 0;
|
||||
double tv = 0;
|
||||
double pp0 = p0[i];
|
||||
if (Double.isNaN(pp0) || Double.isNaN(th0[i])
|
||||
|| Double.isNaN(sh0[i])) {
|
||||
tec = tvc = pc = Double.NaN;
|
||||
} else {
|
||||
double t0 = th0[i] * pow(pp0 / 1000, kapa);
|
||||
tv = t0 * (1 + usetv * 0.000608 * sh0[i]);
|
||||
double b = c0 - log(pp0 / (622.0 / sh0[i] + 0.378));
|
||||
double td = (b - sqrt(b * b - c_1)) / c_2;
|
||||
td -= (t0 - td)
|
||||
* (-0.37329638 + 41.178204 / t0 + 0.0015945203 * td);
|
||||
pc = pp0 * pow(td / t0, kapa_1);
|
||||
tec = adiabatic_te(td, pc);
|
||||
tvc = td * (1 + usetv * 0.000608 * sh0[i]);
|
||||
}
|
||||
|
||||
// Now calculate the virtual temperature of the accending parcel at
|
||||
// the pressures in the input data.
|
||||
double[] tvp = new double[nz];
|
||||
|
||||
for (int k = 0; k < nz; k++) {
|
||||
float pp = p_dat[k * n2 + i];
|
||||
if (Double.isNaN(pc) || Double.isNaN(tec) || Double.isNaN(tvc)
|
||||
|| Double.isNaN(pp)) {
|
||||
tvp[k] = Double.NaN;
|
||||
} else if (pp > pc) {
|
||||
tvp[k] = tvc * pow(pp / pc, kapa);
|
||||
} else {
|
||||
double t0 = tec * pow(pp / pc, kapa);
|
||||
t0 = temp_of_te(t0, pp);
|
||||
tvp[k] = t0
|
||||
* pp
|
||||
/ (pp - usetv
|
||||
* exp(25.687958917 - c1 * t0 - c2 / t0));
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate environment virtual temp, where we force the
|
||||
// environment to be no cooler than dry adiabatic from the ascending
|
||||
// parcel start. Find pressure of min environmetal virtual theta E
|
||||
// above condensation pressure...record temperature and dewpoint
|
||||
// there. Since we do not need the accending parcel temps to
|
||||
// complete the dcape calc, we will put the environmental virtual
|
||||
// temp into the that storage.
|
||||
double wm = Double.NaN;
|
||||
double pm = Double.NaN;
|
||||
double tm = 0;
|
||||
double tdm = 0;
|
||||
for (int k = nzm; k >= 0; k -= 1) {
|
||||
float pp = p_dat[k * n2 + i];
|
||||
float tt = t_dat[k * n2 + i];
|
||||
float td3 = td_dat[k * n2 + i];
|
||||
|
||||
if (Double.isNaN(tvc) || Double.isNaN(pc) || Double.isNaN(pp)
|
||||
|| Double.isNaN(tvp[k]) || Double.isNaN(tt)
|
||||
|| Double.isNaN(td3)) {
|
||||
tvp[k] = Double.NaN;
|
||||
continue;
|
||||
}
|
||||
double t0 = tt;
|
||||
double eee = exp(26.186004814 - c1 * td3 - c2 / td3);
|
||||
double qd = eee / (pp - 0.60771703 * eee);
|
||||
eee = (1 + usetv * 0.608 * qd);
|
||||
double thve = t0 * eee;
|
||||
double pr = pow(pp / pc, kapa);
|
||||
if (thve < tvc * pr) {
|
||||
thve = tvc * pr;
|
||||
t0 = thve / eee;
|
||||
}
|
||||
if (tvp[k] <= thve && Double.isNaN(wm) || pp > pc
|
||||
&& (pm >= 400 || Double.isNaN(pm))) {
|
||||
if (Double.isNaN(pm) && pp < pc) {
|
||||
pm = pc;
|
||||
}
|
||||
tvp[k] = thve;
|
||||
continue;
|
||||
}
|
||||
tvp[k] = thve;
|
||||
thve = (thve + 2529 * qd) * pow(1000 / (pp), kapa);
|
||||
if (thve > wm && pm >= 400) {
|
||||
continue;
|
||||
}
|
||||
wm = thve;
|
||||
pm = pp;
|
||||
tm = t0;
|
||||
tdm = td3;
|
||||
}
|
||||
|
||||
// Here we will reuse our condensation level storage for
|
||||
// the level above current. This loop performs the actual dcape
|
||||
// calculation.
|
||||
double rhm = 0;
|
||||
double qq = 0;
|
||||
double tve1 = tvc;
|
||||
double pp1 = pc;
|
||||
double tvp1 = tec;
|
||||
for (int k = nzm; k >= 0; k -= 1) {
|
||||
double tve = tvp[k];
|
||||
float pp = p_dat[k * n2 + i];
|
||||
|
||||
if (k == nzm) {
|
||||
dcape[i] = Float.NaN;
|
||||
pp1 = tvp1 = tve1 = Double.NaN;
|
||||
}
|
||||
if (k == nzm) {
|
||||
dcape[i] = Float.NaN;
|
||||
pp1 = tvp1 = tve1 = Double.NaN;
|
||||
}
|
||||
if (Double.isNaN(pm) || Double.isNaN(pp0) || Double.isNaN(tv)) {
|
||||
continue;
|
||||
} else if (Double.isNaN(pp0)) {
|
||||
;
|
||||
} else if (pp1 >= pp0) {
|
||||
continue;
|
||||
} else if (Double.isNaN(tve1)) {
|
||||
;
|
||||
} else if (Double.isNaN(pp) || Double.isNaN(tve)) {
|
||||
continue;
|
||||
} else if (pp <= pm) {
|
||||
;
|
||||
} else if (Double.isNaN(wm)) {
|
||||
dcape[i] = 0;
|
||||
} else {
|
||||
|
||||
// Now we finally have the data we need for calculating
|
||||
// the dcape contribution for this layer. If we have not
|
||||
// made any dcape calculations to this point, initialize
|
||||
// the decent parcel.
|
||||
if (Double.isNaN(dcape[i])) {
|
||||
dcape[i] = 0;
|
||||
double eee = exp(26.186004814 - c1 * tdm - c2 / tdm);
|
||||
double qd = eee / (pm - 0.60771703 * eee);
|
||||
double qw = qd + max_evap / 3;
|
||||
double t0 = tm - 2529 * max_evap / 3;
|
||||
eee = exp(26.186004814 - c1 * t0 - c2 / t0);
|
||||
double qs = eee / (pm - 0.60771703 * eee);
|
||||
if (qs >= qw) {
|
||||
wm = max_evap - max_evap / 3;
|
||||
tm = t0;
|
||||
rhm = qw / qs;
|
||||
double b = c0 - log(qw * pm / (0.622 - 0.378 * qw));
|
||||
tdm = (b - sqrt(b * b - c_1)) / c_2;
|
||||
} else {
|
||||
tm = tdm = mytw(tm, tdm, pm);
|
||||
rhm = 1.0;
|
||||
eee = exp(26.186004814 - c1 * tm - c2 / tm);
|
||||
qw = eee / (pm - 0.60771703 * eee);
|
||||
wm = max_evap - (qw - qd);
|
||||
}
|
||||
qq = qw;
|
||||
tvp1 = tm * (1 + usetv * 0.608 * qw);
|
||||
pp1 = pm;
|
||||
}
|
||||
|
||||
// Deal with reaching the surface, add in top of layer part.
|
||||
double pb, dlnp, thve;
|
||||
if (pp > pp0) {
|
||||
pb = pp0;
|
||||
dlnp = log(pb / pp1);
|
||||
thve = tv;
|
||||
} else {
|
||||
pb = pp;
|
||||
dlnp = log(pb / pp1);
|
||||
thve = tve;
|
||||
}
|
||||
double up = -dlnp * 287 * 0.5 * (tvp1 - tve1);
|
||||
if (up < -dcape[i]) {
|
||||
dcape[i] = 0;
|
||||
} else {
|
||||
dcape[i] += up;
|
||||
}
|
||||
|
||||
// Deal with letting parcel fall to pb
|
||||
double pr = pow(pb / pp1, kapa);
|
||||
if (wm <= 0)
|
||||
tvp1 *= pr;
|
||||
else {
|
||||
double rhmx = rhm + (pb - pp1) * (max_rh - rhm)
|
||||
/ (pp0 - pp1);
|
||||
double t0 = tm * pr;
|
||||
double eee = exp(26.186004814 - c1 * t0 - c2 / t0);
|
||||
double qs = eee / (pb - 0.60771703 * eee);
|
||||
if (qq / qs > rhmx) {
|
||||
tm = t0;
|
||||
double b = c0 - log(qq * pb / (0.622 - 0.378 * qq));
|
||||
tdm = (b - sqrt(b * b - c_1)) / c_2;
|
||||
tvp1 *= pr;
|
||||
rhm = qq / qs;
|
||||
} else {
|
||||
double qd = (rhmx * qs - qq)
|
||||
/ sqrt(1000 * (rhmx * qs + qq));
|
||||
if (qd > wm) {
|
||||
qd = wm;
|
||||
}
|
||||
double qw = qq + qd;
|
||||
double td = t0 - 2529 * wm;
|
||||
eee = exp(26.186004814 - c1 * td - c2 / td);
|
||||
qs = eee / (pb - 0.60771703 * eee);
|
||||
if (qs >= qw) {
|
||||
tm = td;
|
||||
rhm = qw / qs;
|
||||
double b = c0
|
||||
- log(qw * pb / (0.622 - 0.378 * qw));
|
||||
tdm = (b - sqrt(b * b - c_1)) / c_2;
|
||||
} else {
|
||||
double b = c0
|
||||
- log(qq * pb / (0.622 - 0.378 * qq));
|
||||
tdm = (b - sqrt(b * b - c_1)) / c_2;
|
||||
tm = tdm = mytw(t0, tdm, pb);
|
||||
rhm = 1.0;
|
||||
eee = exp(26.186004814 - c1 * tm - c2 / tm);
|
||||
qw = eee / (pb - 0.60771703 * eee);
|
||||
qd = qw - qq;
|
||||
}
|
||||
wm -= qd;
|
||||
qq = qw;
|
||||
tvp1 = tm * (1 + usetv * 0.608 * qw);
|
||||
}
|
||||
}
|
||||
|
||||
// Add contribution of bottom of layer.
|
||||
double dn = -dlnp * 287 * 0.5 * (tvp1 - thve);
|
||||
if (dn < -dcape[i])
|
||||
dcape[i] = 0;
|
||||
else
|
||||
dcape[i] += dn;
|
||||
|
||||
}
|
||||
|
||||
// Make current layer top next layer bottom.
|
||||
tve1 = tve;
|
||||
pp1 = pp;
|
||||
}
|
||||
}
|
||||
return new PythonNumpyFloatArray(dcape, ny, nx);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue