diff --git a/doc/source/tutorial.rst b/doc/source/tutorial.rst index 29b68b5..dc05e15 100644 --- a/doc/source/tutorial.rst +++ b/doc/source/tutorial.rst @@ -12,7 +12,7 @@ Upcoming Tutorials .. toctree:: :maxdepth: 1 - tutorials/tutorial_03_2018.rst + tutorials/wrf_workshop_2018.rst Past Tutorials @@ -22,4 +22,5 @@ Past Tutorials :maxdepth: 1 tutorials/wrf_workshop_2017.rst + tutorials/tutorial_03_2018.rst diff --git a/doc/source/tutorials/wrf_workshop_2017.rst b/doc/source/tutorials/wrf_workshop_2017.rst index 44ecbef..d01ed3b 100644 --- a/doc/source/tutorials/wrf_workshop_2017.rst +++ b/doc/source/tutorials/wrf_workshop_2017.rst @@ -1,5 +1,5 @@ -WRF Workshop 2017 -===================== +WRF Users' Workshop 2017 +========================== Welcome wrf-python tutorial attendees! diff --git a/doc/source/tutorials/wrf_workshop_2018.rst b/doc/source/tutorials/wrf_workshop_2018.rst new file mode 100644 index 0000000..9a8f8e8 --- /dev/null +++ b/doc/source/tutorials/wrf_workshop_2018.rst @@ -0,0 +1,410 @@ +WRF Users' Workshop 2018 +========================= + +Welcome WRF-Python tutorial attendees! + +The instructions below should be completed prior to arriving at the tutorial. +There will not be enough time to do this during the tutorial. + +Prerequisites +--------------- + +This tutorial assumes that you have basic knowledge of how to type commands +in to a command terminal using your preferred operating system. You +should know some basic directory commands like *cd*, *mkdir*, *cp*, *mv*. + +Regarding Python, to understand the examples in this tutorial, you +should have some experience with Python basics. Below is a list of some +Python concepts that you will see in the examples, but don't worry if you aren't +familiar with everything. + +- Opening a Python interpreter and entering commands. +- Importing packages via the import statement. +- Familiarity with some of the basic Python types: str, list, tuple, dict, bool, float, int, None. +- Creating a list, tuple, or dict with "[ ]", "( )", "{ }" syntax (e.g. my_list = [1,2,3,4,5]). +- Accessing dict/list/tuple items with the "x[ ]" syntax (e.g. my_list_item = my_list[0]). +- Slicing str/list/tuple with the ":" syntax (e.g. my_slice = my_list[1:3]). +- Using object methods and attributes with the "x.y" syntax (e.g. my_list.append(6)). +- Calling functions (e.g. result = some_function(x, y)) +- Familiarity with numpy would be helpful, as only a very brief introduction + is provided. +- Familiarity with matplotlib would be helpful, as only a very brief + introduction is provided. + +If you are completely new to Python, that shouldn't be a problem, since +most of the examples consist of basic container types and function calls. It +would be helpful to look at some introductory material before arriving at the +tutorial. If you've programmed before, picking up Python isn't too difficult. + +Here are some links: + +https://www.learnpython.org/ + +https://developers.google.com/edu/python/ + + +Step 1: Open a Command Terminal +-------------------------------- + +To begin, you will first need to know how to open a command line terminal for +your operating system. + +For Windows: + +.. code-block:: none + + WINDOWS + r + type cmd in the run window + +For Mac: + +.. code-block:: none + + Finder -> Applications -> Utilities -> Terminal + +For Linux: + +.. code-block:: none + + Try one of the following: + + CTRL + ALT + T + CTRL + ALT + F2 + + +Step 2: Download Miniconda +---------------------------- + +For this tutorial, you will need to download and install Miniconda. We are +going to use Python 2.7, but it will also work with Python 3.5+. + +Please use the appropriate link below to download Miniconda for your operating +system. + +.. note:: + + 64-bit OS recommended + +`Win64 `_ + +`Mac `_ + +`Linux `_ + +For more information, see: https://conda.io/miniconda.html + +.. note:: + + **What is Miniconda?** + + If you have used the Anaconda distribution for Python before, then you will be + familiar with Miniconda. The Anaconda Python distribution includes numerous + scientific packages out of the box, which can be difficult for users to build and + install. More importantly, Anaconda includes the conda package manager. + + The conda package manager is a utility (similar to yum or apt-get) that installs + packages from a repository of pre-compiled Python packages. These repositories + are called channels. Conda makes it easy for Python users to install and + uninstall packages, and also can be used to create isolated Python environments + (more on that later). + + Miniconda is a bare bones implementation of Anaconda and only includes the + conda package manager. Since we are going to use the conda-forge channel to + install our scientific packages, Miniconda avoids any complications between + packages provided by Anaconda and conda-forge. + + +Step 3: Install Miniconda +---------------------------- + +Windows: + + 1. Browse to the directory where you downloaded Miniconda2-latest-Windows-x86_64.exe. + + 2. Double click on Miniconda2-latest-Windows-x86_64.exe. + + 3. Follow the instructions. + + 4. Restart your command terminal. + +Mac and Linux: + + For Mac and Linux, the installer is a bash script. + + 1. Using a terminal, you need to execute the bash shell script that you downloaded by + doing:: + + bash /path/to/Miniconda2-latest-MacOSX-x86_64.sh [Mac] + + bash /path/to/Miniconda2-latest-Linux-x86_64.sh [Linux] + + 2. Follow the instructions. + + 3. At the end of the installation, it will ask if you want to add the + miniconda2 path to your bash environment. If you are unsure what to do, + you should say "yes". If you say "no", we're going to assume you know + what you are doing. + + If you said "yes", then once you restart your shell, the miniconda2 Python + will be found instead of the system Python when you type the "python" + command. If you want to undo this later, then you can edit + either ~/.bash_profile or ~/.bashrc (depending on OS used) and + comment out the line that looks similar to:: + + # added by Miniconda2 4.1.11 installer + export PATH="/path/to/miniconda2/bin:$PATH" + + 4. Restart your command terminal. + + 5. [Linux and Mac Users Only] Miniconda only works with bash. If bash is + not your default shell, then you need to activate the bash shell by typing + the following in to your command terminal:: + + bash + + 6. Verify that your system is using the correct Python interpreter by typing + the following in to your command terminal:: + + which python + + You should see the path to your miniconda installation. If not, see the + note below. + + .. note:: + + If you have already installed another Python distribution, like Enthought + Canopy, you will need to comment out any PATH entries for that distribution + in your .bashrc or .bash_profile. Otherwise, your shell environment may + pick to wrong Python installation. + + If bash is not your default shell type, and the PATH variable has been + set in .bash_profile by the miniconda installer, try executing + "bash -l" instead of the "bash" command in step 5. + + +Step 4: Set Up the Conda Environment +-------------------------------------- + +If you are new to the conda package manager, one of the nice features of conda +is that you can create isolated Python environments that prevent package +incompatibilities. This is similar to the *virtualenv* package that some +Python users may be familiar with. However, conda is not compatible with +virtualenv, so only use conda environments when working with conda. + +The name of our conda environment for this tutorial is: **tutorial_2018**. + +Follow the instructions below to create the tutorial_2018 environment. + + 1. Open a command terminal if you haven't done so. + + 2. [Linux and Mac Users Only] The conda package manager only works with bash, + so if bash is not your current shell, type:: + + bash + + 3. Add the conda-forge channel to your conda package manager. + + Type or copy the command below in to your command terminal. You should + run this command even if you have already done it in the past. + This will ensure that conda-forge is set as the highest priority channel. + + :: + + conda config --add channels conda-forge + + .. note:: + + Conda-forge is a community driven collection of packages that are + continually tested to ensure compatibility. We highly recommend using + conda-forge when working with conda. See https://conda-forge.github.io/ + for more details on this excellent project. + + 4. Create the conda environment for the tutorial. + + Type or copy this command in to your command terminal:: + + conda create -n tutorial_2018 python=2.7 matplotlib cartopy netcdf4 jupyter git ffmpeg wrf-python + + Type "y" when prompted. It will take several minutes to install everything. + + This command creates an isolated Python environment named *tutorial_2018*, and installs + the python interpreter, matplotlib, cartopy, netcdf4, jupyter, git, ffmpeg, and wrf-python + packages. + + .. note:: + + When the installation completes, your command terminal might post a message similar to: + + .. code-block:: none + + If this is your first install of dbus, automatically load on login with: + + mkdir -p ~/Library/LaunchAgents + cp /path/to/miniconda2/envs/tutorial_test/org.freedesktop.dbus-session.plist ~/Library/LaunchAgents/ + launchctl load -w ~/Library/LaunchAgents/org.freedesktop.dbus-session.plist + + This is indicating that the dbus package can be set up to automatically load on login. You + can either ignore this message or type in the commands as indicated on your command terminal. + The tutorial should work fine in either case. + + 5. Activate the conda environment. + + To activate the tutorial_2018 Python environment, type the following + in to the command terminal: + + For Linux and Mac (using bash):: + + source activate tutorial_2018 + + For Windows:: + + activate tutorial_2018 + + You should see (tutorial_2018) on your command prompt. + + To deactivate your conda environment, type the following in to the + command terminal: + + For Linux and Mac:: + + source deactivate + + For Windows:: + + deactivate tutorial_2018 + + +Step 5: Download the Student Workbook +--------------------------------------- + +The student workbook for the tutorial is available on GitHub. The tutorial_2018 +conda environment includes the git application needed to download the repository. + +These instructions download the tutorial in to your home directory. If you want +to place the tutorial in to another directory, we're going to assume you know +how to do this yourself. + +To download the student workbook, follow these instructions: + + 1. Activate the tutorial_2018 conda environment following the instructions + in the previous step (*source activate tutorial_2018* or + *activate tutorial_2018*). + + 2. Change your working directory to the home directory by typing the + following command in to the command terminal: + + For Linux and Mac:: + + cd ~ + + For Windows:: + + cd %HOMEPATH% + + 3. Download the git repository for the tutorial by typing the following + in to the command terminal:: + + git clone https://github.com/NCAR/wrf_python_tutorial.git + + 4. There may be additional changes to the tutorial after you have downloaded + it. To pull down the latest changes, type the following in to the + command terminal: + + For Linux and Mac:: + + source activate tutorial_2018 + + cd ~/wrf_python_tutorial/wrf_workshop_2018 + + git pull + + For Windows:: + + activate tutorial_2018 + + cd %HOMEPATH%\wrf_python_tutorial\wrf_workshop_2018 + + git pull + + .. note:: + + If you try the "git pull" command and it returns an error indicating + that you have made changes to the workbook, this is probably because + you ran the workbook and it contains the cell output. To fix this, + first do a checkout of the workbook, then do the pull. + + .. code-block:: none + + git checkout -- . + git pull + + +Step 6: Verify Your Environment +---------------------------------- + +Verifying that your environment is correct involves importing a few +packages and checking for errors (you may see some warnings for matplotlib +or xarray, but you can safely ignore these). + + 1. Activate the tutorial_2018 conda environment if it isn't already active + (see instructions above). + + 2. Open a python terminal by typing the following in to the command + terminal:: + + python + + 3. Now type the following in to the Python interpreter:: + + >>> import netCDF4 + >>> import matplotlib + >>> import xarray + >>> import wrf + + 4. You can exit the Python interpreter using **CTRL + D** + + +Step 7: Obtain WRF Output Files +---------------------------------- + +For this tutorial, we strongly recommend that you use your own WRF output files. +The tutorial includes an easy way to point to your own data files. The WRF +output files should all be from the same WRF run and use the same domain. +If your files are located on another system (e.g. yellowstone), then copy 2 or +3 of these files to your local computer prior to the tutorial. + +If you do not have any of your own WRF output files, then you can download the +instructor data files from a link that should have been provided to you in an +email prior to the tutorial. + +If you are using the link provided in the email for your data, you can follow +the instructions below to place your data in the default location for your +workbook. + + 1. The link in the email should take you to a location on an Amazon cloud + drive. + + 2. If you hover your mouse over the wrf_tutorial_data.zip file, you'll see + an empty check box appear next to the file name. Click this check + box. + + 3. At the bottom of the screen, you'll see a Download button next to a + cloud icon. Click this button to start the download. + + 4. The download was most likely placed in to your ~/Downloads folder + [%HOMEPATH%\\Downloads for Windows]. Using your preferred method of choice + for unzipping files, unzip this file in to your home directory. Your data + should now be in ~/wrf_tutorial_data + [%HOMEPATH%\\wrf_tutorial_data for Windows]. + + 5. Verify that you have three WRF output files in that directory. + + +Getting Help +---------------- + +If you experience problems during this installation, please send a question +to the :ref:`google-group` support mailing list. + + +We look forward to seeing you at the tutorial! diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 164d905..f181285 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -714,10 +714,10 @@ def smooth2d(field, passes, meta=True): @set_cape_alg_metadata(is2d=True, copyarg="pres_hpa") def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, missing=default_fill(np.float64), meta=True): - """Return the two-dimensional CAPE, CIN, LCL, and LFC. + """Return the two-dimensional MCAPE, MCIN, LCL, and LFC. This function calculates the maximum convective available potential - energy (CAPE), maximum convective inhibition (CIN), + energy (MCAPE), maximum convective inhibition (MCIN), lifted condensation level (LCL), and level of free convection (LFC). This function uses the RIP [Read/Interpolate/plot] code to calculate potential energy (CAPE) and convective inhibition @@ -725,18 +725,28 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, in the column (i.e. something akin to Colman's MCAPE). CAPE is defined as the accumulated buoyant energy from the level of free convection (LFC) to the equilibrium level (EL). CIN is defined as the accumulated negative - buoyant energy from the parcel starting point to the LFC. The word 'parcel' - here refers to a 500 meter deep parcel, with actual temperature and - moisture averaged over that depth. + buoyant energy from the parcel starting point to the LFC. + + The cape_2d algorithm works by first finding the maximum theta-e height + level in the lowest 3000 m. A parcel with a depth of 500 m is then + calculated and centered over this maximum theta-e height level. The + parcel's moisture and temperature characteristics are calculated by + averaging over the depth of this 500 m parcel. This 'maximum' parcel + is then used to compute MCAPE, MCIN, LCL and LFC. The leftmost dimension of the returned array represents four different quantities: - - return_val[0,...] will contain CAPE [J kg-1] - - return_val[1,...] will contain CIN [J kg-1] + - return_val[0,...] will contain MCAPE [J kg-1] + - return_val[1,...] will contain MCIN [J kg-1] - return_val[2,...] will contain LCL [m] - return_val[3,...] will contain LFC [m] + This function also supports computing MCAPE along a single vertical + column. In this mode, the *pres_hpa*, *tkel*, *qv* and *height* arguments + must be one-dimensional vertical columns, and the *terrain* and + *psfc_hpa* arguments must be scalar values + (:obj:`float`, :class:`numpy.float32` or :class:`numpy.float64`). This is the raw computational algorithm and does not extract any variables from WRF output files. Use :meth:`wrf.getvar` to both extract and compute @@ -749,6 +759,9 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, least three dimensions. The rightmost dimensions can be top_bottom x south_north x west_east or bottom_top x south_north x west_east. + When operating on only a single column of values, the vertical + column can be bottom_top or top_bottom. In this case, *terrain* + and *psfc_hpa* must be scalars. Note: @@ -774,12 +787,16 @@ def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, terrain (:class:`xarray.DataArray` or :class:`numpy.ndarray`): Terrain height in [m]. This is at least a two-dimensional array with the same dimensionality as *pres_hpa*, excluding the vertical - (bottom_top/top_bottom) dimension. + (bottom_top/top_bottom) dimension. When operating on a single + vertical column, this argument must be a scalar (:obj:`float`, + :class:`numpy.float32`, or :class:`numpy.float64`). psfc_hpa (:class:`xarray.DataArray` or :class:`numpy.ndarray`): The surface pressure in [hPa]. This is at least a two-dimensional array with the same dimensionality as *pres_hpa*, excluding the - vertical (bottom_top/top_bottom) dimension. + vertical (bottom_top/top_bottom) dimension. When operating on a + singlevertical column, this argument must be a scalar + (:obj:`float`, :class:`numpy.float32`, or :class:`numpy.float64`). Note: @@ -1111,12 +1128,16 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, True. Default is :data:`wrf.default_fill(numpy.float64)`. - opt_thresh (:obj:`float`, optional): The required amount of optical - depth (looking from top down) required to trigger a cloud top - temperature calculation. Values less than this threshold will be - treated as no cloud areas. For thin cirrus, a value of .1 might be - appropriate. For large CB clouds, 1000.0 might be appropriate. - Default is 1.0. + opt_thresh (:obj:`float`, optional): The amount of optical + depth (integrated from top down) required to trigger a cloud top + temperature calculation. The cloud top temperature is calculated at + the vertical level where this threshold is met. Vertical columns + with less than this threshold will be treated as cloud free areas. + In general, the larger the value is for this + threshold, the lower the altitude will be for the cloud top + temperature calculation, and therefore higher cloud top + temperature values. Default is 1.0, which should be sufficient for + most users. meta (:obj:`bool`): Set to False to disable metadata and return :class:`numpy.ndarray` instead of diff --git a/src/wrf/g_cape.py b/src/wrf/g_cape.py index 5f5ef10..0a43c7a 100755 --- a/src/wrf/g_cape.py +++ b/src/wrf/g_cape.py @@ -13,13 +13,13 @@ from .metadecorators import set_cape_metadata @set_cape_metadata(is2d=True) def get_2dcape(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, meta=True, _key=None, missing=default_fill(np.float64)): - """Return the 2d fields of CAPE, CIN, LCL, and LFC. + """Return the 2d fields of MCAPE, MCIN, LCL, and LFC. The leftmost dimension of the returned array represents four different quantities: - - return_val[0,...] will contain CAPE [J kg-1] - - return_val[1,...] will contain CIN [J kg-1] + - return_val[0,...] will contain MCAPE [J kg-1] + - return_val[1,...] will contain MCIN [J kg-1] - return_val[2,...] will contain LCL [m] - return_val[3,...] will contain LFC [m] diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index b35d57d..40e8e49 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -1985,7 +1985,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): p = argvals[copyarg] missing = argvals["missing"] - # Note: cape_3d supports using only a single column of data + # Note: 2D/3D cape supports using only a single column of data is1d = p.ndim == 1 # Need to squeeze the right dimensions for 1D cape @@ -1997,31 +1997,34 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outattrs = OrderedDict() + if is2d: - outname = "cape_2d" + if is1d: + outname = "cape_2d" + else: + outname = "cape_2d_column" outattrs["description"] = "mcape ; mcin ; lcl ; lfc" outattrs["units"] = "J kg-1 ; J kg-1 ; m ; m" outattrs["MemoryOrder"] = "XY" + outattrs["MemoryOrder"] = "" else: if not is1d: outname = "cape_3d" - outattrs["description"] = "cape; cin" - outattrs["units"] = "J kg-1 ; J kg-1" outattrs["MemoryOrder"] = "XYZ" else: - outname = "cape_column" - outattrs["description"] = "cape; cin" - outattrs["units"] = "J kg-1 ; J kg-1" + outname = "cape_3d_column" outattrs["MemoryOrder"] = "Z" + outattrs["description"] = "cape; cin" + outattrs["units"] = "J kg-1 ; J kg-1" if isinstance(p, DataArray): if is2d: - # Right dims - outdims[-2:] = p.dims[-2:] - # Left dims - outdims[1:-2] = p.dims[0:-3] - + if not is1d: + # Right dims + outdims[-2:] = p.dims[-2:] + # Left dims + outdims[1:-2] = p.dims[0:-3] else: if not is1d: # Right dims @@ -2030,6 +2033,7 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outdims[1:-3] = p.dims[0:-3] else: outdims[1] = p.dims[0] + outcoords = {} # Left-most is always cape_cin or cape_cin_lcl_lfc diff --git a/src/wrf/specialdec.py b/src/wrf/specialdec.py index f72899b..1795f25 100644 --- a/src/wrf/specialdec.py +++ b/src/wrf/specialdec.py @@ -229,7 +229,14 @@ def cape_left_iter(alg_dtype=np.float64): ter_follow = args[8] is2d = i3dflag == 0 - is1d = np.isscalar(sfp) + # Note: This should still work with DataArrays + is1d = np.isscalar(sfp) or np.size(sfp) == 1 + + # Make sure sfp and terrain are regular floats for 1D case + # This should also work with DataArrays + if is1d: + ter = float(ter) + sfp = float(sfp) orig_dtype = p_hpa.dtype @@ -542,39 +549,27 @@ def check_cape_args(): ter_follow = args[8] is2d = False if i3dflag != 0 else True + is1d = ((np.isscalar(sfp) or np.size(sfp) == 1) or + (np.isscalar(ter) or np.size(ter) == 1)) if not (p_hpa.shape == tk.shape == qv.shape == ht.shape): raise ValueError("arguments 0, 1, 2, 3 must be the same shape") # 2D CAPE does not allow for scalars - if is2d: - if np.isscalar(ter) or np.isscalar(sfp): - raise ValueError("arguments 4 and 5 cannot be scalars with " - "cape_2d routine") - + if not is1d: if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1: raise ValueError("arguments 4 and 5 must have " "{} dimensions".format(p_hpa.ndim-1)) - - # 3D cape is allowed to be just a vertical column with scalars - # for terrain and psfc_hpa. else: - if np.isscalar(ter) and np.isscalar(sfp): - if p_hpa.ndim != 1: - raise ValueError("arguments 0-3 " - "must be 1-dimensional when " - "arguments 4 and 5 are scalars") - if is2d: - raise ValueError("argument 7 must be 0 when using 1D data") - else: - if ((np.isscalar(ter) and not np.isscalar(sfp)) or - (not np.isscalar(ter) and np.isscalar(sfp))): - raise ValueError("arguments 4 and 5 must both be scalars") - else: - if ter.ndim != p_hpa.ndim-1 or sfp.ndim != p_hpa.ndim-1: - raise ValueError("arguments 4 and 5 " - "must have {} dimensions".format( - p_hpa.ndim-1)) + if np.size(ter) != np.size(sfp): + raise ValueError("arguments 4 and 5 must both be scalars or " + "both be arrays") + + # Only need to test p_hpa since we assured args 0-3 have same ndim + if p_hpa.ndim != 1: + raise ValueError("arguments 0-3 " + "must be 1-dimensional when " + "arguments 4 and 5 are scalars") return wrapped(*args, **kwargs) diff --git a/test/ci_tests/ci_result_file.nc b/test/ci_tests/ci_result_file.nc index 2cf9ffb..e15d100 100644 Binary files a/test/ci_tests/ci_result_file.nc and b/test/ci_tests/ci_result_file.nc differ diff --git a/test/comp_utest.py b/test/comp_utest.py index 958cd75..864bb7a 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -597,9 +597,71 @@ def test_cape3d_1d(wrfnc): nt.assert_allclose(to_np(result), to_np(ref)) return func + + +def test_cape2d_1d(wrfnc): + + def func(self): + varnames = ("T", "P", "PB", "QVAPOR", "PH", "PHB", "HGT", "PSFC") + ncvars = extract_vars(wrfnc, 0, varnames, method="cat", squeeze=True, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + ter = ncvars["HGT"] + psfc = ncvars["PSFC"] + + col_idxs = (slice(None), t.shape[-2]//2, t.shape[-1]//2) + t = t[col_idxs] + p = p[col_idxs] + pb = pb[col_idxs] + qv = qv[col_idxs] + ph = ph[col_idxs] + phb = phb[col_idxs] + ter = float(ter[col_idxs[1:]]) + psfc = float(psfc[col_idxs[1:]]) + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -1) + z = geopt_unstag/Constants.G + + # Convert pressure to hPa + p_hpa = ConversionFactors.PA_TO_HPA * full_p + psfc_hpa = ConversionFactors.PA_TO_HPA * psfc + + i3dflag = 0 + ter_follow = 1 + + result = cape_2d(p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) + + print ("RESULT", result) + + ref = getvar(wrfnc, "cape_2d") + + ref = ref[(slice(None),) + col_idxs[1:]] + + print ("REF", ref) + + nt.assert_allclose(to_np(result), to_np(ref)) + + return func if __name__ == "__main__": + from wrf import (omp_set_num_threads, omp_set_schedule, omp_get_schedule, + omp_set_dynamic, omp_get_num_procs, OMP_SCHED_STATIC) + omp_set_num_threads(omp_get_num_procs()-1) + omp_set_schedule(OMP_SCHED_STATIC, 0) + omp_set_dynamic(False) + varnames = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", @@ -639,6 +701,9 @@ if __name__ == "__main__": func = test_cape3d_1d(wrfnc) setattr(WRFVarsTest, "test_cape3d_1d", func) + func = test_cape2d_1d(wrfnc) + setattr(WRFVarsTest, "test_cape2d_1d", func) + ut.main()