From 4404e55ab87105da2c75c4c1a7703a4596988c1f Mon Sep 17 00:00:00 2001 From: Ben Steffensmeier Date: Thu, 6 Jun 2013 15:19:27 -0500 Subject: [PATCH] Issue #2043 Port Dcape to java. Change-Id: Ibcdd38fc47bb57eb6b2c6d0921d9790c2f1893ba Former-commit-id: 6be12509e5b214f702f9981b7f6f3d035e33061f [formerly 6be12509e5b214f702f9981b7f6f3d035e33061f [formerly eddf7c71bb23057fe4e68a7462a634c6d81799d9]] Former-commit-id: 660b04ba32ceca229a292fd9b98fca1816f2bda8 Former-commit-id: 98dfa3f8b344a8cbe09e90ed49611062c7f132bd --- .../derivedParameters/functions/Dcape.py | 96 +---- .../python/function/AdiabeticTemperature.java | 5 +- .../derivparam/python/function/CalcTw.java | 136 +++++++ .../derivparam/python/function/CapeFunc.java | 21 +- .../derivparam/python/function/Constants.java | 64 ++++ .../derivparam/python/function/DCapeFunc.java | 334 ++++++++++++++++++ 6 files changed, 551 insertions(+), 105 deletions(-) create mode 100644 cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CalcTw.java create mode 100644 cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/Constants.java create mode 100644 cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/DCapeFunc.java diff --git a/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Dcape.py b/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Dcape.py index 69413dd048..bb996fe7ab 100644 --- a/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Dcape.py +++ b/cave/com.raytheon.uf.viz.derivparam.python/localization/derivedParameters/functions/Dcape.py @@ -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 diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java index 27e88a0457..6b715a07d7 100644 --- a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java +++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/AdiabeticTemperature.java @@ -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); } diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CalcTw.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CalcTw.java new file mode 100644 index 0000000000..8d308f51dc --- /dev/null +++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CalcTw.java @@ -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. + * + * 
+ * 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.
+ *
+ * + * 
+ * 
+ * SOFTWARE HISTORY
+ * 
+ * Date         Ticket#    Engineer    Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 06, 2013  2043       bsteffen    Ported from meteolib fortran
+ * 
+ *
+ * + * @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; + + } +} diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java index dfe0b72abc..6d15e7af4b 100644 --- a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java +++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/CapeFunc.java @@ -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) { diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/Constants.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/Constants.java new file mode 100644 index 0000000000..3dbe94953d --- /dev/null +++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/Constants.java @@ -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. + * + * 
+ * 
+ * SOFTWARE HISTORY
+ * 
+ * Date         Ticket#    Engineer    Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 06, 2013  2043       bsteffen    Ported from meteolib C
+ * 
+ *
+ * + * @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; +} diff --git a/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/DCapeFunc.java b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/DCapeFunc.java new file mode 100644 index 0000000000..e562efbe2d --- /dev/null +++ b/cave/com.raytheon.uf.viz.derivparam.python/src/com/raytheon/uf/viz/derivparam/python/function/DCapeFunc.java @@ -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. + * + * 
+ * 
+ * SOFTWARE HISTORY
+ * 
+ * Date         Ticket#    Engineer    Description
+ * ------------ ---------- ----------- --------------------------
+ * Jun 05, 2013  2043       bsteffen    Ported from meteolib C
+ * 
+ *
+ * + * @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); + } + + +}