forked from 3rdparty/wrf-python
65 changed files with 31413 additions and 1039 deletions
Binary file not shown.
@ -0,0 +1,230 @@
@@ -0,0 +1,230 @@
|
||||
# Makefile for Sphinx documentation
|
||||
#
|
||||
|
||||
# You can set these variables from the command line.
|
||||
SPHINXOPTS =
|
||||
SPHINXBUILD = sphinx-build
|
||||
PAPER =
|
||||
BUILDDIR = build
|
||||
|
||||
# User-friendly check for sphinx-build
|
||||
ifeq ($(shell which $(SPHINXBUILD) >/dev/null 2>&1; echo $$?), 1) |
||||
$(error The '$(SPHINXBUILD)' command was not found. Make sure you have Sphinx installed, then set the SPHINXBUILD environment variable to point to the full path of the '$(SPHINXBUILD)' executable. Alternatively you can add the directory with the executable to your PATH. If you don\'t have Sphinx installed, grab it from http://sphinx-doc.org/)
|
||||
endif |
||||
|
||||
# Internal variables.
|
||||
PAPEROPT_a4 = -D latex_paper_size=a4
|
||||
PAPEROPT_letter = -D latex_paper_size=letter
|
||||
ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
|
||||
# the i18n builder cannot share the environment and doctrees with the others
|
||||
I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
|
||||
|
||||
.PHONY: help |
||||
help: |
||||
@echo "Please use \`make <target>' where <target> is one of"
|
||||
@echo " html to make standalone HTML files"
|
||||
@echo " dirhtml to make HTML files named index.html in directories"
|
||||
@echo " singlehtml to make a single large HTML file"
|
||||
@echo " pickle to make pickle files"
|
||||
@echo " json to make JSON files"
|
||||
@echo " htmlhelp to make HTML files and a HTML help project"
|
||||
@echo " qthelp to make HTML files and a qthelp project"
|
||||
@echo " applehelp to make an Apple Help Book"
|
||||
@echo " devhelp to make HTML files and a Devhelp project"
|
||||
@echo " epub to make an epub"
|
||||
@echo " epub3 to make an epub3"
|
||||
@echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
|
||||
@echo " latexpdf to make LaTeX files and run them through pdflatex"
|
||||
@echo " latexpdfja to make LaTeX files and run them through platex/dvipdfmx"
|
||||
@echo " text to make text files"
|
||||
@echo " man to make manual pages"
|
||||
@echo " texinfo to make Texinfo files"
|
||||
@echo " info to make Texinfo files and run them through makeinfo"
|
||||
@echo " gettext to make PO message catalogs"
|
||||
@echo " changes to make an overview of all changed/added/deprecated items"
|
||||
@echo " xml to make Docutils-native XML files"
|
||||
@echo " pseudoxml to make pseudoxml-XML files for display purposes"
|
||||
@echo " linkcheck to check all external links for integrity"
|
||||
@echo " doctest to run all doctests embedded in the documentation (if enabled)"
|
||||
@echo " coverage to run coverage check of the documentation (if enabled)"
|
||||
@echo " dummy to check syntax errors of document sources"
|
||||
|
||||
.PHONY: clean |
||||
clean: |
||||
rm -rf $(BUILDDIR)/*
|
||||
|
||||
.PHONY: html |
||||
html: |
||||
$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
|
||||
@echo
|
||||
@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
|
||||
|
||||
.PHONY: dirhtml |
||||
dirhtml: |
||||
$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
|
||||
@echo
|
||||
@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
|
||||
|
||||
.PHONY: singlehtml |
||||
singlehtml: |
||||
$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
|
||||
@echo
|
||||
@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
|
||||
|
||||
.PHONY: pickle |
||||
pickle: |
||||
$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
|
||||
@echo
|
||||
@echo "Build finished; now you can process the pickle files."
|
||||
|
||||
.PHONY: json |
||||
json: |
||||
$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
|
||||
@echo
|
||||
@echo "Build finished; now you can process the JSON files."
|
||||
|
||||
.PHONY: htmlhelp |
||||
htmlhelp: |
||||
$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
|
||||
@echo
|
||||
@echo "Build finished; now you can run HTML Help Workshop with the" \
|
||||
".hhp project file in $(BUILDDIR)/htmlhelp."
|
||||
|
||||
.PHONY: qthelp |
||||
qthelp: |
||||
$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
|
||||
@echo
|
||||
@echo "Build finished; now you can run "qcollectiongenerator" with the" \
|
||||
".qhcp project file in $(BUILDDIR)/qthelp, like this:"
|
||||
@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/wrf-python.qhcp"
|
||||
@echo "To view the help file:"
|
||||
@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/wrf-python.qhc"
|
||||
|
||||
.PHONY: applehelp |
||||
applehelp: |
||||
$(SPHINXBUILD) -b applehelp $(ALLSPHINXOPTS) $(BUILDDIR)/applehelp
|
||||
@echo
|
||||
@echo "Build finished. The help book is in $(BUILDDIR)/applehelp."
|
||||
@echo "N.B. You won't be able to view it unless you put it in" \
|
||||
"~/Library/Documentation/Help or install it in your application" \
|
||||
"bundle."
|
||||
|
||||
.PHONY: devhelp |
||||
devhelp: |
||||
$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
|
||||
@echo
|
||||
@echo "Build finished."
|
||||
@echo "To view the help file:"
|
||||
@echo "# mkdir -p $$HOME/.local/share/devhelp/wrf-python"
|
||||
@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/wrf-python"
|
||||
@echo "# devhelp"
|
||||
|
||||
.PHONY: epub |
||||
epub: |
||||
$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
|
||||
@echo
|
||||
@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
|
||||
|
||||
.PHONY: epub3 |
||||
epub3: |
||||
$(SPHINXBUILD) -b epub3 $(ALLSPHINXOPTS) $(BUILDDIR)/epub3
|
||||
@echo
|
||||
@echo "Build finished. The epub3 file is in $(BUILDDIR)/epub3."
|
||||
|
||||
.PHONY: latex |
||||
latex: |
||||
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
|
||||
@echo
|
||||
@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
|
||||
@echo "Run \`make' in that directory to run these through (pdf)latex" \
|
||||
"(use \`make latexpdf' here to do that automatically)."
|
||||
|
||||
.PHONY: latexpdf |
||||
latexpdf: |
||||
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
|
||||
@echo "Running LaTeX files through pdflatex..."
|
||||
$(MAKE) -C $(BUILDDIR)/latex all-pdf
|
||||
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
|
||||
|
||||
.PHONY: latexpdfja |
||||
latexpdfja: |
||||
$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
|
||||
@echo "Running LaTeX files through platex and dvipdfmx..."
|
||||
$(MAKE) -C $(BUILDDIR)/latex all-pdf-ja
|
||||
@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
|
||||
|
||||
.PHONY: text |
||||
text: |
||||
$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
|
||||
@echo
|
||||
@echo "Build finished. The text files are in $(BUILDDIR)/text."
|
||||
|
||||
.PHONY: man |
||||
man: |
||||
$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
|
||||
@echo
|
||||
@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
|
||||
|
||||
.PHONY: texinfo |
||||
texinfo: |
||||
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
|
||||
@echo
|
||||
@echo "Build finished. The Texinfo files are in $(BUILDDIR)/texinfo."
|
||||
@echo "Run \`make' in that directory to run these through makeinfo" \
|
||||
"(use \`make info' here to do that automatically)."
|
||||
|
||||
.PHONY: info |
||||
info: |
||||
$(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo
|
||||
@echo "Running Texinfo files through makeinfo..."
|
||||
make -C $(BUILDDIR)/texinfo info
|
||||
@echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo."
|
||||
|
||||
.PHONY: gettext |
||||
gettext: |
||||
$(SPHINXBUILD) -b gettext $(I18NSPHINXOPTS) $(BUILDDIR)/locale
|
||||
@echo
|
||||
@echo "Build finished. The message catalogs are in $(BUILDDIR)/locale."
|
||||
|
||||
.PHONY: changes |
||||
changes: |
||||
$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
|
||||
@echo
|
||||
@echo "The overview file is in $(BUILDDIR)/changes."
|
||||
|
||||
.PHONY: linkcheck |
||||
linkcheck: |
||||
$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
|
||||
@echo
|
||||
@echo "Link check complete; look for any errors in the above output " \
|
||||
"or in $(BUILDDIR)/linkcheck/output.txt."
|
||||
|
||||
.PHONY: doctest |
||||
doctest: |
||||
$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
|
||||
@echo "Testing of doctests in the sources finished, look at the " \
|
||||
"results in $(BUILDDIR)/doctest/output.txt."
|
||||
|
||||
.PHONY: coverage |
||||
coverage: |
||||
$(SPHINXBUILD) -b coverage $(ALLSPHINXOPTS) $(BUILDDIR)/coverage
|
||||
@echo "Testing of coverage in the sources finished, look at the " \
|
||||
"results in $(BUILDDIR)/coverage/python.txt."
|
||||
|
||||
.PHONY: xml |
||||
xml: |
||||
$(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml
|
||||
@echo
|
||||
@echo "Build finished. The XML files are in $(BUILDDIR)/xml."
|
||||
|
||||
.PHONY: pseudoxml |
||||
pseudoxml: |
||||
$(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml
|
||||
@echo
|
||||
@echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml."
|
||||
|
||||
.PHONY: dummy |
||||
dummy: |
||||
$(SPHINXBUILD) -b dummy $(ALLSPHINXOPTS) $(BUILDDIR)/dummy
|
||||
@echo
|
||||
@echo "Build finished. Dummy builder generates no files."
|
@ -0,0 +1,281 @@
@@ -0,0 +1,281 @@
|
||||
@ECHO OFF |
||||
|
||||
REM Command file for Sphinx documentation |
||||
|
||||
if "%SPHINXBUILD%" == "" ( |
||||
set SPHINXBUILD=sphinx-build |
||||
) |
||||
set BUILDDIR=build |
||||
set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% source |
||||
set I18NSPHINXOPTS=%SPHINXOPTS% source |
||||
if NOT "%PAPER%" == "" ( |
||||
set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS% |
||||
set I18NSPHINXOPTS=-D latex_paper_size=%PAPER% %I18NSPHINXOPTS% |
||||
) |
||||
|
||||
if "%1" == "" goto help |
||||
|
||||
if "%1" == "help" ( |
||||
:help |
||||
echo.Please use `make ^<target^>` where ^<target^> is one of |
||||
echo. html to make standalone HTML files |
||||
echo. dirhtml to make HTML files named index.html in directories |
||||
echo. singlehtml to make a single large HTML file |
||||
echo. pickle to make pickle files |
||||
echo. json to make JSON files |
||||
echo. htmlhelp to make HTML files and a HTML help project |
||||
echo. qthelp to make HTML files and a qthelp project |
||||
echo. devhelp to make HTML files and a Devhelp project |
||||
echo. epub to make an epub |
||||
echo. epub3 to make an epub3 |
||||
echo. latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter |
||||
echo. text to make text files |
||||
echo. man to make manual pages |
||||
echo. texinfo to make Texinfo files |
||||
echo. gettext to make PO message catalogs |
||||
echo. changes to make an overview over all changed/added/deprecated items |
||||
echo. xml to make Docutils-native XML files |
||||
echo. pseudoxml to make pseudoxml-XML files for display purposes |
||||
echo. linkcheck to check all external links for integrity |
||||
echo. doctest to run all doctests embedded in the documentation if enabled |
||||
echo. coverage to run coverage check of the documentation if enabled |
||||
echo. dummy to check syntax errors of document sources |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "clean" ( |
||||
for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i |
||||
del /q /s %BUILDDIR%\* |
||||
goto end |
||||
) |
||||
|
||||
|
||||
REM Check if sphinx-build is available and fallback to Python version if any |
||||
%SPHINXBUILD% 1>NUL 2>NUL |
||||
if errorlevel 9009 goto sphinx_python |
||||
goto sphinx_ok |
||||
|
||||
:sphinx_python |
||||
|
||||
set SPHINXBUILD=python -m sphinx.__init__ |
||||
%SPHINXBUILD% 2> nul |
||||
if errorlevel 9009 ( |
||||
echo. |
||||
echo.The 'sphinx-build' command was not found. Make sure you have Sphinx |
||||
echo.installed, then set the SPHINXBUILD environment variable to point |
||||
echo.to the full path of the 'sphinx-build' executable. Alternatively you |
||||
echo.may add the Sphinx directory to PATH. |
||||
echo. |
||||
echo.If you don't have Sphinx installed, grab it from |
||||
echo.http://sphinx-doc.org/ |
||||
exit /b 1 |
||||
) |
||||
|
||||
:sphinx_ok |
||||
|
||||
|
||||
if "%1" == "html" ( |
||||
%SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The HTML pages are in %BUILDDIR%/html. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "dirhtml" ( |
||||
%SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "singlehtml" ( |
||||
%SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "pickle" ( |
||||
%SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished; now you can process the pickle files. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "json" ( |
||||
%SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished; now you can process the JSON files. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "htmlhelp" ( |
||||
%SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished; now you can run HTML Help Workshop with the ^ |
||||
.hhp project file in %BUILDDIR%/htmlhelp. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "qthelp" ( |
||||
%SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished; now you can run "qcollectiongenerator" with the ^ |
||||
.qhcp project file in %BUILDDIR%/qthelp, like this: |
||||
echo.^> qcollectiongenerator %BUILDDIR%\qthelp\wrf-python.qhcp |
||||
echo.To view the help file: |
||||
echo.^> assistant -collectionFile %BUILDDIR%\qthelp\wrf-python.ghc |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "devhelp" ( |
||||
%SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "epub" ( |
||||
%SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The epub file is in %BUILDDIR%/epub. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "epub3" ( |
||||
%SPHINXBUILD% -b epub3 %ALLSPHINXOPTS% %BUILDDIR%/epub3 |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The epub3 file is in %BUILDDIR%/epub3. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "latex" ( |
||||
%SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished; the LaTeX files are in %BUILDDIR%/latex. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "latexpdf" ( |
||||
%SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex |
||||
cd %BUILDDIR%/latex |
||||
make all-pdf |
||||
cd %~dp0 |
||||
echo. |
||||
echo.Build finished; the PDF files are in %BUILDDIR%/latex. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "latexpdfja" ( |
||||
%SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex |
||||
cd %BUILDDIR%/latex |
||||
make all-pdf-ja |
||||
cd %~dp0 |
||||
echo. |
||||
echo.Build finished; the PDF files are in %BUILDDIR%/latex. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "text" ( |
||||
%SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The text files are in %BUILDDIR%/text. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "man" ( |
||||
%SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The manual pages are in %BUILDDIR%/man. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "texinfo" ( |
||||
%SPHINXBUILD% -b texinfo %ALLSPHINXOPTS% %BUILDDIR%/texinfo |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The Texinfo files are in %BUILDDIR%/texinfo. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "gettext" ( |
||||
%SPHINXBUILD% -b gettext %I18NSPHINXOPTS% %BUILDDIR%/locale |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The message catalogs are in %BUILDDIR%/locale. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "changes" ( |
||||
%SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.The overview file is in %BUILDDIR%/changes. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "linkcheck" ( |
||||
%SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Link check complete; look for any errors in the above output ^ |
||||
or in %BUILDDIR%/linkcheck/output.txt. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "doctest" ( |
||||
%SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Testing of doctests in the sources finished, look at the ^ |
||||
results in %BUILDDIR%/doctest/output.txt. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "coverage" ( |
||||
%SPHINXBUILD% -b coverage %ALLSPHINXOPTS% %BUILDDIR%/coverage |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Testing of coverage in the sources finished, look at the ^ |
||||
results in %BUILDDIR%/coverage/python.txt. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "xml" ( |
||||
%SPHINXBUILD% -b xml %ALLSPHINXOPTS% %BUILDDIR%/xml |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The XML files are in %BUILDDIR%/xml. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "pseudoxml" ( |
||||
%SPHINXBUILD% -b pseudoxml %ALLSPHINXOPTS% %BUILDDIR%/pseudoxml |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. The pseudo-XML files are in %BUILDDIR%/pseudoxml. |
||||
goto end |
||||
) |
||||
|
||||
if "%1" == "dummy" ( |
||||
%SPHINXBUILD% -b dummy %ALLSPHINXOPTS% %BUILDDIR%/dummy |
||||
if errorlevel 1 exit /b 1 |
||||
echo. |
||||
echo.Build finished. Dummy builder generates no files. |
||||
goto end |
||||
) |
||||
|
||||
:end |
@ -0,0 +1,288 @@
@@ -0,0 +1,288 @@
|
||||
# -*- coding: utf-8 -*- |
||||
# |
||||
# wrf-python documentation build configuration file, created by |
||||
# sphinx-quickstart on Wed Jun 29 14:57:16 2016. |
||||
# |
||||
# This file is execfile()d with the current directory set to its |
||||
# containing dir. |
||||
# |
||||
# Note that not all possible configuration values are present in this |
||||
# autogenerated file. |
||||
# |
||||
# All configuration values have a default; values that are commented out |
||||
# serve to show the default. |
||||
|
||||
import sys |
||||
import os |
||||
|
||||
# If extensions (or modules to document with autodoc) are in another directory, |
||||
# add these directories to sys.path here. If the directory is relative to the |
||||
# documentation root, use os.path.abspath to make it absolute, like shown here. |
||||
#sys.path.insert(0, os.path.abspath('.')) |
||||
|
||||
# -- General configuration ------------------------------------------------ |
||||
|
||||
# If your documentation needs a minimal Sphinx version, state it here. |
||||
#needs_sphinx = '1.0' |
||||
|
||||
# Add any Sphinx extension module names here, as strings. They can be |
||||
# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom |
||||
# ones. |
||||
extensions = [ |
||||
'sphinx.ext.autodoc', 'sphinx.ext.napoleon' |
||||
] |
||||
|
||||
# Add any paths that contain templates here, relative to this directory. |
||||
templates_path = ['_templates'] |
||||
|
||||
# The suffix(es) of source filenames. |
||||
# You can specify multiple suffix as a list of string: |
||||
# source_suffix = ['.rst', '.md'] |
||||
source_suffix = '.rst' |
||||
|
||||
# The encoding of source files. |
||||
#source_encoding = 'utf-8-sig' |
||||
|
||||
# The master toctree document. |
||||
master_doc = 'index' |
||||
|
||||
# General information about the project. |
||||
project = u'wrf-python' |
||||
copyright = u'2016, Bill Ladwig' |
||||
author = u'Bill Ladwig' |
||||
|
||||
# The version info for the project you're documenting, acts as replacement for |
||||
# |version| and |release|, also used in various other places throughout the |
||||
# built documents. |
||||
# |
||||
# The short X.Y version. |
||||
version = u'0.0.1' |
||||
# The full version, including alpha/beta/rc tags. |
||||
release = u'0.0.1' |
||||
|
||||
# The language for content autogenerated by Sphinx. Refer to documentation |
||||
# for a list of supported languages. |
||||
# |
||||
# This is also used if you do content translation via gettext catalogs. |
||||
# Usually you set "language" from the command line for these cases. |
||||
language = None |
||||
|
||||
# There are two options for replacing |today|: either, you set today to some |
||||
# non-false value, then it is used: |
||||
#today = '' |
||||
# Else, today_fmt is used as the format for a strftime call. |
||||
#today_fmt = '%B %d, %Y' |
||||
|
||||
# List of patterns, relative to source directory, that match files and |
||||
# directories to ignore when looking for source files. |
||||
# This patterns also effect to html_static_path and html_extra_path |
||||
exclude_patterns = [] |
||||
|
||||
# The reST default role (used for this markup: `text`) to use for all |
||||
# documents. |
||||
#default_role = None |
||||
|
||||
# If true, '()' will be appended to :func: etc. cross-reference text. |
||||
#add_function_parentheses = True |
||||
|
||||
# If true, the current module name will be prepended to all description |
||||
# unit titles (such as .. function::). |
||||
#add_module_names = True |
||||
|
||||
# If true, sectionauthor and moduleauthor directives will be shown in the |
||||
# output. They are ignored by default. |
||||
#show_authors = False |
||||
|
||||
# The name of the Pygments (syntax highlighting) style to use. |
||||
pygments_style = 'sphinx' |
||||
|
||||
# A list of ignored prefixes for module index sorting. |
||||
#modindex_common_prefix = [] |
||||
|
||||
# If true, keep warnings as "system message" paragraphs in the built documents. |
||||
#keep_warnings = False |
||||
|
||||
# If true, `todo` and `todoList` produce output, else they produce nothing. |
||||
todo_include_todos = False |
||||
|
||||
|
||||
# -- Options for HTML output ---------------------------------------------- |
||||
|
||||
# The theme to use for HTML and HTML Help pages. See the documentation for |
||||
# a list of builtin themes. |
||||
html_theme = 'alabaster' |
||||
|
||||
# Theme options are theme-specific and customize the look and feel of a theme |
||||
# further. For a list of options available for each theme, see the |
||||
# documentation. |
||||
#html_theme_options = {} |
||||
|
||||
# Add any paths that contain custom themes here, relative to this directory. |
||||
#html_theme_path = [] |
||||
|
||||
# The name for this set of Sphinx documents. |
||||
# "<project> v<release> documentation" by default. |
||||
#html_title = u'wrf-python v0.0.1' |
||||
|
||||
# A shorter title for the navigation bar. Default is the same as html_title. |
||||
#html_short_title = None |
||||
|
||||
# The name of an image file (relative to this directory) to place at the top |
||||
# of the sidebar. |
||||
#html_logo = None |
||||
|
||||
# The name of an image file (relative to this directory) to use as a favicon of |
||||
# the docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32 |
||||
# pixels large. |
||||
#html_favicon = None |
||||
|
||||
# Add any paths that contain custom static files (such as style sheets) here, |
||||
# relative to this directory. They are copied after the builtin static files, |
||||
# so a file named "default.css" will overwrite the builtin "default.css". |
||||
html_static_path = ['_static'] |
||||
|
||||
# Add any extra paths that contain custom files (such as robots.txt or |
||||
# .htaccess) here, relative to this directory. These files are copied |
||||
# directly to the root of the documentation. |
||||
#html_extra_path = [] |
||||
|
||||
# If not None, a 'Last updated on:' timestamp is inserted at every page |
||||
# bottom, using the given strftime format. |
||||
# The empty string is equivalent to '%b %d, %Y'. |
||||
#html_last_updated_fmt = None |
||||
|
||||
# If true, SmartyPants will be used to convert quotes and dashes to |
||||
# typographically correct entities. |
||||
#html_use_smartypants = True |
||||
|
||||
# Custom sidebar templates, maps document names to template names. |
||||
#html_sidebars = {} |
||||
|
||||
# Additional templates that should be rendered to pages, maps page names to |
||||
# template names. |
||||
#html_additional_pages = {} |
||||
|
||||
# If false, no module index is generated. |
||||
#html_domain_indices = True |
||||
|
||||
# If false, no index is generated. |
||||
#html_use_index = True |
||||
|
||||
# If true, the index is split into individual pages for each letter. |
||||
#html_split_index = False |
||||
|
||||
# If true, links to the reST sources are added to the pages. |
||||
#html_show_sourcelink = True |
||||
|
||||
# If true, "Created using Sphinx" is shown in the HTML footer. Default is True. |
||||
#html_show_sphinx = True |
||||
|
||||
# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True. |
||||
#html_show_copyright = True |
||||
|
||||
# If true, an OpenSearch description file will be output, and all pages will |
||||
# contain a <link> tag referring to it. The value of this option must be the |
||||
# base URL from which the finished HTML is served. |
||||
#html_use_opensearch = '' |
||||
|
||||
# This is the file name suffix for HTML files (e.g. ".xhtml"). |
||||
#html_file_suffix = None |
||||
|
||||
# Language to be used for generating the HTML full-text search index. |
||||
# Sphinx supports the following languages: |
||||
# 'da', 'de', 'en', 'es', 'fi', 'fr', 'hu', 'it', 'ja' |
||||
# 'nl', 'no', 'pt', 'ro', 'ru', 'sv', 'tr', 'zh' |
||||
#html_search_language = 'en' |
||||
|
||||
# A dictionary with options for the search language support, empty by default. |
||||
# 'ja' uses this config value. |
||||
# 'zh' user can custom change `jieba` dictionary path. |
||||
#html_search_options = {'type': 'default'} |
||||
|
||||
# The name of a javascript file (relative to the configuration directory) that |
||||
# implements a search results scorer. If empty, the default will be used. |
||||
#html_search_scorer = 'scorer.js' |
||||
|
||||
# Output file base name for HTML help builder. |
||||
htmlhelp_basename = 'wrf-pythondoc' |
||||
|
||||
# -- Options for LaTeX output --------------------------------------------- |
||||
|
||||
latex_elements = { |
||||
# The paper size ('letterpaper' or 'a4paper'). |
||||
#'papersize': 'letterpaper', |
||||
|
||||
# The font size ('10pt', '11pt' or '12pt'). |
||||
#'pointsize': '10pt', |
||||
|
||||
# Additional stuff for the LaTeX preamble. |
||||
#'preamble': '', |
||||
|
||||
# Latex figure (float) alignment |
||||
#'figure_align': 'htbp', |
||||
} |
||||
|
||||
# Grouping the document tree into LaTeX files. List of tuples |
||||
# (source start file, target name, title, |
||||
# author, documentclass [howto, manual, or own class]). |
||||
latex_documents = [ |
||||
(master_doc, 'wrf-python.tex', u'wrf-python Documentation', |
||||
u'Bill Ladwig', 'manual'), |
||||
] |
||||
|
||||
# The name of an image file (relative to this directory) to place at the top of |
||||
# the title page. |
||||
#latex_logo = None |
||||
|
||||
# For "manual" documents, if this is true, then toplevel headings are parts, |
||||
# not chapters. |
||||
#latex_use_parts = False |
||||
|
||||
# If true, show page references after internal links. |
||||
#latex_show_pagerefs = False |
||||
|
||||
# If true, show URL addresses after external links. |
||||
#latex_show_urls = False |
||||
|
||||
# Documents to append as an appendix to all manuals. |
||||
#latex_appendices = [] |
||||
|
||||
# If false, no module index is generated. |
||||
#latex_domain_indices = True |
||||
|
||||
|
||||
# -- Options for manual page output --------------------------------------- |
||||
|
||||
# One entry per manual page. List of tuples |
||||
# (source start file, name, description, authors, manual section). |
||||
man_pages = [ |
||||
(master_doc, 'wrf-python', u'wrf-python Documentation', |
||||
[author], 1) |
||||
] |
||||
|
||||
# If true, show URL addresses after external links. |
||||
#man_show_urls = False |
||||
|
||||
|
||||
# -- Options for Texinfo output ------------------------------------------- |
||||
|
||||
# Grouping the document tree into Texinfo files. List of tuples |
||||
# (source start file, target name, title, author, |
||||
# dir menu entry, description, category) |
||||
texinfo_documents = [ |
||||
(master_doc, 'wrf-python', u'wrf-python Documentation', |
||||
author, 'wrf-python', 'One line description of project.', |
||||
'Miscellaneous'), |
||||
] |
||||
|
||||
# Documents to append as an appendix to all manuals. |
||||
#texinfo_appendices = [] |
||||
|
||||
# If false, no module index is generated. |
||||
#texinfo_domain_indices = True |
||||
|
||||
# How to display URL addresses: 'footnote', 'no', or 'inline'. |
||||
#texinfo_show_urls = 'footnote' |
||||
|
||||
# If true, do not generate a @detailmenu in the "Top" node's menu. |
||||
#texinfo_no_detailmenu = False |
@ -0,0 +1,36 @@
@@ -0,0 +1,36 @@
|
||||
.. wrf-python documentation master file, created by |
||||
sphinx-quickstart on Wed Jun 29 14:57:16 2016. |
||||
You can adapt this file completely to your liking, but it should at least |
||||
contain the root `toctree` directive. |
||||
|
||||
Welcome to wrf-python's documentation! |
||||
====================================== |
||||
|
||||
Contents: |
||||
|
||||
.. toctree:: |
||||
:maxdepth: 2 |
||||
|
||||
Extraction Routine |
||||
------------------ |
||||
.. autofunction:: wrf.getvar |
||||
|
||||
Interpolation Routines |
||||
---------------------- |
||||
|
||||
.. autofunction:: wrf.interplevel |
||||
|
||||
.. autofunction:: wrf.vertcross |
||||
|
||||
.. autofunction:: wrf.interpline |
||||
|
||||
.. autofunction:: wrf.vinterp |
||||
|
||||
|
||||
Indices and tables |
||||
================== |
||||
|
||||
* :ref:`genindex` |
||||
* :ref:`modindex` |
||||
* :ref:`search` |
||||
|
@ -0,0 +1,126 @@
@@ -0,0 +1,126 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & |
||||
vs, w, uh, tem1, tem2) |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: uh |
||||
|
||||
INTEGER, INTENT(IN) :: nx, ny, nz, nzp1 |
||||
REAL(KIND=8), DIMENSION(nx,ny,nzp1), INTENT(IN) :: zp |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: mapfct |
||||
REAL(KIND=8), INTENT(IN) :: dx, dy |
||||
REAL(KIND=8), INTENT(IN) :: uhmnhgt, uhmxhgt |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: us |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: vs |
||||
REAL(KIND=8), DIMENSION(nx,ny,nzp1), INTENT(IN) :: w |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: uh |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: tem1 |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: tem2 |
||||
|
||||
!NCLEND |
||||
|
||||
! Misc local variables |
||||
INTEGER :: i, j, k, kbot, ktop |
||||
REAL(KIND=8) :: twodx, twody, wgtlw, sum, wmean, wsum, wavg |
||||
REAL(KIND=8) :: helbot, heltop, wbot, wtop |
||||
REAL(KIND=8) :: zbot, ztop |
||||
|
||||
! Initialize arrays |
||||
uh = 0.0 |
||||
tem1 = 0.0 |
||||
|
||||
! Calculate vertical component of helicity at scalar points |
||||
! us: u at scalar points |
||||
! vs: v at scalar points |
||||
|
||||
twodx = 2.0*dx |
||||
twody = 2.0*dy |
||||
DO k=2,nz-2 |
||||
DO j=2,ny-1 |
||||
DO i=2,nx-1 |
||||
wavg = 0.5*(w(i,j,k)+w(i,j,k+1)) |
||||
tem1(i,j,k) = wavg * ((vs(i+1,j,k) - vs(i-1,j,k))/(twodx * mapfct(i,j)) - & |
||||
(us(i,j+1,k) - us(i,j-1,k))/(twody * mapfct(i,j))) |
||||
tem2(i,j,k) = 0.5*(zp(i,j,k) + zp(i,j,k+1)) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
! Integrate over depth uhminhgt to uhmxhgt AGL |
||||
! |
||||
! WRITE(6,'(a,f12.1,a,f12.1,a)') & |
||||
! 'Calculating UH from ',uhmnhgt,' to ',uhmxhgt,' m AGL' |
||||
DO j=2,ny-2 |
||||
DO i=2,nx-2 |
||||
zbot = zp(i,j,2) + uhmnhgt |
||||
ztop = zp(i,j,2) + uhmxhgt |
||||
! |
||||
! Find wbar, weighted-mean vertical velocity in column |
||||
! Find w at uhmnhgt AGL (bottom) |
||||
! |
||||
DO k=2,nz-3 |
||||
IF(zp(i,j,k) > zbot) EXIT |
||||
END DO |
||||
kbot = k |
||||
wgtlw = (zp(i,j,kbot) - zbot)/(zp(i,j,kbot) - zp(i,j,kbot-1)) |
||||
wbot = (wgtlw*w(i,j,kbot-1)) + ((1.-wgtlw)*w(i,j,kbot)) |
||||
|
||||
! Find w at uhmxhgt AGL (top) |
||||
DO k=2,nz-3 |
||||
IF(zp(i,j,k) > ztop) EXIT |
||||
END DO |
||||
ktop = k |
||||
wgtlw = (zp(i,j,ktop) - ztop)/(zp(i,j,ktop) - zp(i,j,ktop-1)) |
||||
wtop = (wgtlw*w(i,j,ktop-1)) + ((1.-wgtlw)*w(i,j,ktop)) |
||||
|
||||
! First part, uhmnhgt to kbot |
||||
wsum = 0.5*(w(i,j,kbot) + wbot) * (zp(i,j,kbot) - zbot) |
||||
|
||||
! Integrate up through column |
||||
DO k=(kbot+1),(ktop-1) |
||||
wsum = wsum + 0.5*(w(i,j,k) + w(i,j,k-1)) * (zp(i,j,k) - zp(i,j,k-1)) |
||||
END DO |
||||
|
||||
! Last part, ktop-1 to uhmxhgt |
||||
wsum = wsum + 0.5*(wtop + w(i,j,ktop-1)) * (ztop - zp(i,j,ktop-1)) |
||||
wmean = wsum/(uhmxhgt - uhmnhgt) |
||||
|
||||
IF (wmean > 0.) THEN ! column updraft, not downdraft |
||||
|
||||
! Find helicity at uhmnhgt AGL (bottom) |
||||
DO k=2,nz-3 |
||||
IF (tem2(i,j,k) > zbot) EXIT |
||||
END DO |
||||
kbot = k |
||||
wgtlw = (tem2(i,j,kbot) - zbot)/(tem2(i,j,kbot) - tem2(i,j,kbot-1)) |
||||
helbot = (wgtlw*tem1(i,j,kbot-1)) + ((1.-wgtlw)*tem1(i,j,kbot)) |
||||
|
||||
! Find helicity at uhmxhgt AGL (top) |
||||
DO k=2,nz-3 |
||||
IF (tem2(i,j,k) > ztop) EXIT |
||||
END DO |
||||
ktop = k |
||||
wgtlw = (tem2(i,j,ktop) - ztop)/(tem2(i,j,ktop) - tem2(i,j,ktop-1)) |
||||
heltop = (wgtlw*tem1(i,j,ktop-1)) + ((1.-wgtlw)*tem1(i,j,ktop)) |
||||
|
||||
! First part, uhmnhgt to kbot |
||||
sum = 0.5*(tem1(i,j,kbot) + helbot) * (tem2(i,j,kbot) - zbot) |
||||
|
||||
! Integrate up through column |
||||
DO k=(kbot+1),(ktop-1) |
||||
sum = sum + 0.5*(tem1(i,j,k) + tem1(i,j,k-1)) * (tem2(i,j,k) - tem2(i,j,k-1)) |
||||
END DO |
||||
|
||||
! Last part, ktop-1 to uhmxhgt |
||||
uh(i,j) = sum + 0.5*(heltop + tem1(i,j,ktop-1)) * (ztop - tem2(i,j,ktop-1)) |
||||
END IF |
||||
END DO |
||||
END DO |
||||
|
||||
uh = uh * 1000. ! Scale according to Kain et al. (2008) |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCALCUH |
@ -0,0 +1,56 @@
@@ -0,0 +1,56 @@
|
||||
! NCLFORTSTART |
||||
SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew) |
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: lowc, midc, highc |
||||
|
||||
INTEGER nz, ns, ew |
||||
REAL(KIND=8), DIMENSION(ew, ns, nz), INTENT(IN) :: pres, rh |
||||
REAL(KIND=8), DIMENSION(ew, ns), INTENT(OUT) :: lowc, midc, highc |
||||
|
||||
! NCLEND |
||||
|
||||
INTEGER i, j, k |
||||
INTEGER kchi, kcmi, kclo |
||||
|
||||
! Remove compiler warnings |
||||
kchi = 0 |
||||
kcmi = 0 |
||||
kclo = 0 |
||||
|
||||
DO j = 1,ns |
||||
DO i = 1,ew |
||||
DO k = 1,nz-1 |
||||
IF ( pres(i,j,k) .GT. 97000. ) kclo=k |
||||
IF ( pres(i,j,k) .GT. 80000. ) kcmi=k |
||||
IF ( pres(i,j,k) .GT. 45000. ) kchi=k |
||||
END DO |
||||
|
||||
DO k = 1,nz-1 |
||||
IF (k .GE. kclo .AND. k .LT. kcmi) THEN |
||||
lowc(i,j) = AMAX1(rh(i,j,k), lowc(i,j)) |
||||
ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud |
||||
midc(i,j) = AMAX1(rh(i,j,k), midc(i,j)) |
||||
ELSE if (k .GE. kchi) THEN ! high cloud |
||||
highc(i,j) = AMAX1(rh(i,j,k), highc(i,j)) |
||||
END IF |
||||
END DO |
||||
|
||||
lowc(i,j) = 4.0 * lowc(i,j)/100. - 3.0 |
||||
midc(i,j) = 4.0 * midc(i,j)/100. - 3.0 |
||||
highc(i,j) = 2.5 * highc(i,j)/100. - 1.5 |
||||
|
||||
lowc(i,j) = amin1(lowc(i,j), 1.0) |
||||
lowc(i,j) = amax1(lowc(i,j), 0.0) |
||||
midc(i,j) = amin1(midc(i,j), 1.0) |
||||
midc(i,j) = amax1(midc(i,j), 0.0) |
||||
highc(i,j) = amin1(highc(i,j), 1.0) |
||||
highc(i,j) = amax1(highc(i,j), 0.0) |
||||
|
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCLOUDFRAC |
@ -1,10 +1,66 @@
@@ -1,10 +1,66 @@
|
||||
! These are chosen to match the wrf module_model_constants.F where |
||||
! applicable |
||||
MODULE constants |
||||
INTEGER :: ERRLEN=512 |
||||
REAL(KIND=8), PARAMETER :: P1000MB=100000.D0 |
||||
REAL(KIND=8), PARAMETER :: R_D=287.D0 |
||||
REAL(KIND=8), PARAMETER :: CP=7.D0*R_D/2.D0 |
||||
REAL(KIND=8), PARAMETER :: R=287.04D0 |
||||
REAL(KIND=8), PARAMETER :: G=9.81D0 |
||||
REAL(KIND=8), PARAMETER :: GAMMA=0.0065D0 |
||||
INTEGER :: ALGERR=64 |
||||
|
||||
REAL(KIND=8), PARAMETER :: WRF_EARTH_RADIUS = 6370000.D0 |
||||
REAL(KIND=8), PARAMETER :: T_BASE = 300.0 |
||||
REAL(KIND=8), PARAMETER :: PI = 3.1415926535897932384626433D0 |
||||
REAL(KIND=8), PARAMETER :: RAD_PER_DEG = PI/180.D0 |
||||
REAL(KIND=8), PARAMETER :: DEG_PER_RAD = 180.D0/PI |
||||
REAL(KIND=8), PARAMETER :: DEFAULT_FILL = 9.9692099683868690D36 |
||||
|
||||
REAL(KIND=8), PARAMETER :: P1000MB = 100000.D0 |
||||
! j/k/kg |
||||
REAL(KIND=8), PARAMETER :: RD = 287.D0 |
||||
REAL(KIND=8), PARAMETER :: RV = 461.6D0 |
||||
!REAL(KIND=8), PARAMETER :: RV = 461.5D0 |
||||
! j/k/kg note: not using bolton's value of 1005.7 |
||||
REAL(KIND=8), PARAMETER :: CP = 1004.5D0 |
||||
!REAL(KIND=8), PARAMETER :: CP = 1004D0 |
||||
|
||||
REAL(KIND=8), PARAMETER :: G = 9.81D0 |
||||
REAL(KIND=8), PARAMETER :: USSALR = 0.0065D0 ! deg C per m |
||||
|
||||
REAL(KIND=8), PARAMETER :: CELKEL = 273.15D0 |
||||
REAL(KIND=8), PARAMETER :: CELKEL_TRIPLE = 273.16D0 |
||||
!REAL(KIND=8), PARAMETER :: GRAV = 9.81D0 |
||||
! hpa |
||||
REAL(KIND=8), PARAMETER :: EZERO = 6.112D0 |
||||
REAL(KIND=8), PARAMETER :: ESLCON1 = 17.67D0 |
||||
REAL(KIND=8), PARAMETER :: ESLCON2 = 29.65D0 |
||||
REAL(KIND=8), PARAMETER :: EPS = 0.622D0 |
||||
REAL(KIND=8), PARAMETER :: GAMMA = RD/CP |
||||
! cp_moist=cp*(1.+cpmd*qvp) |
||||
REAL(KIND=8), PARAMETER :: CPMD = .887D0 |
||||
! rgas_moist=rgas*(1.+rgasmd*qvp) |
||||
REAL(KIND=8), PARAMETER :: RGASMD = .608D0 |
||||
! gamma_moist=gamma*(1.+gammamd*qvp) |
||||
REAL(KIND=8), PARAMETER :: GAMMAMD = RGASMD - CPMD |
||||
REAL(KIND=8), PARAMETER :: TLCLC1 = 2840.D0 |
||||
REAL(KIND=8), PARAMETER :: TLCLC2 = 3.5D0 |
||||
REAL(KIND=8), PARAMETER :: TLCLC3 = 4.805D0 |
||||
REAL(KIND=8), PARAMETER :: TLCLC4 = 55.D0 |
||||
! k |
||||
REAL(KIND=8), PARAMETER :: THTECON1 = 3376.D0 |
||||
REAL(KIND=8), PARAMETER :: THTECON2 = 2.54D0 |
||||
REAL(KIND=8), PARAMETER :: THTECON3 = .81D0 |
||||
|
||||
REAL(KIND=8), PARAMETER :: ABSCOEFI = .272D0 ! cloud ice absorption coefficient in m^2/g |
||||
REAL(KIND=8), PARAMETER :: ABSCOEF = .145D0 ! cloud water absorption coefficient in m^2/g |
||||
|
||||
REAL(KIND=8), PARAMETER :: GAMMA_SEVEN = 720.D0 |
||||
REAL(KIND=8), PARAMETER :: RHOWAT = 1000.D0 |
||||
REAL(KIND=8), PARAMETER :: RHO_R = RHOWAT |
||||
REAL(KIND=8), PARAMETER :: RHO_S = 100.D0 |
||||
REAL(KIND=8), PARAMETER :: RHO_G = 400.D0 |
||||
REAL(KIND=8), PARAMETER :: ALPHA = 0.224D0 |
||||
|
||||
REAL(KIND=8), PARAMETER :: SCLHT = RD*256.D0/G |
||||
REAL(KIND=8), PARAMETER :: EXPON = RD*USSALR/G |
||||
REAL(KIND=8), PARAMETER :: EXPONI = 1./EXPON |
||||
|
||||
|
||||
END MODULE constants |
||||
|
||||
|
@ -0,0 +1,51 @@
@@ -0,0 +1,51 @@
|
||||
! Theta-e |
||||
!NCLFORTSTART |
||||
SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh) |
||||
USE constants, ONLY : EPS, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, TLCLC4, & |
||||
THTECON1, THTECON2, THTECON3 |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: eth |
||||
|
||||
! Input variables |
||||
! Sizes |
||||
INTEGER,INTENT(IN) :: miy, mjx, mkzh |
||||
! Qvapor [g/kg] |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: qvp |
||||
! Temperature [K] |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk |
||||
! full pressure (=P+PB) [hPa] |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs |
||||
! Output variable |
||||
! equivalent potential temperature [K] |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: eth |
||||
|
||||
!NCLEND |
||||
|
||||
! local variables |
||||
REAL(KIND=8) :: q |
||||
REAL(KIND=8) :: t |
||||
REAL(KIND=8) :: p |
||||
REAL(KIND=8) :: e |
||||
REAL(KIND=8) :: tlcl |
||||
INTEGER :: i, j, k |
||||
|
||||
DO k = 1,mkzh |
||||
DO j = 1,mjx |
||||
DO i = 1,miy |
||||
q = MAX(qvp(i,j,k), 1.D-15) |
||||
t = tmk(i,j,k) |
||||
p = prs(i,j,k)/100. |
||||
e = q*p/(EPS + q) |
||||
tlcl = TLCLC1/(LOG(t**TLCLC2/e) - TLCLC3) + TLCLC4 |
||||
eth(i,j,k) = t*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q))* & |
||||
EXP((THTECON1/tlcl - THTECON2)*q*(1.D0 + THTECON3*q)) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DEQTHECALC |
@ -0,0 +1,617 @@
@@ -0,0 +1,617 @@
|
||||
!====================================================================== |
||||
! |
||||
! !IROUTINE: VIRTUAL -- Calculate virtual temperature (K) |
||||
! |
||||
! !DESCRIPTION: |
||||
! |
||||
! This function returns a single value of virtual temperature in |
||||
! K, given temperature in K and mixing ratio in kg/kg. For an |
||||
! array of virtual temperatures, use subroutine VIRTUAL_TEMP. |
||||
! |
||||
! !INPUT: |
||||
! RATMIX - water vapor mixing ratio (kg/kg) |
||||
! TEMP - temperature (K) |
||||
! |
||||
! !OUTPUT: |
||||
! TV - Virtual temperature (K) |
||||
! |
||||
|
||||
! NCLFORTSTART |
||||
REAL(KIND=8) FUNCTION TVIRTUAL(temp, ratmix) |
||||
USE constants, ONLY : EPS |
||||
|
||||
!f2py threadsafe |
||||
|
||||
IMPLICIT NONE |
||||
REAL(KIND=8), INTENT(IN) :: temp, ratmix |
||||
|
||||
! NCLEND |
||||
|
||||
TVIRTUAL = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix)) |
||||
|
||||
RETURN |
||||
|
||||
END FUNCTION TVIRTUAL |
||||
|
||||
! NCLFORTSTART |
||||
REAL(KIND=8) FUNCTION TONPSADIABAT(thte,prs,psadithte,psadiprs,psaditmk,gamma,& |
||||
errstat, errmsg) |
||||
USE constants, ONLY : ALGERR |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: cape, cin |
||||
|
||||
IMPLICIT NONE |
||||
REAL(KIND=8), INTENT(IN) :: thte |
||||
REAL(KIND=8), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(150), INTENT(IN) :: psadithte |
||||
REAL(KIND=8), DIMENSION(150), INTENT(IN) :: psadiprs |
||||
REAL(KIND=8), DIMENSION(150,150), INTENT(IN) :: psaditmk |
||||
REAL(KIND=8), INTENT(IN) :: gamma |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
! NCLEND |
||||
|
||||
REAL(KIND=8) :: fracjt |
||||
REAL(KIND=8) :: fracjt2 |
||||
REAL(KIND=8) :: fracip |
||||
REAL(KIND=8) :: fracip2 |
||||
|
||||
INTEGER :: ip, ipch, jt, jtch |
||||
|
||||
! This function gives the temperature (in K) on a moist adiabat |
||||
! (specified by thte in K) given pressure in hPa. It uses a |
||||
! lookup table, with data that was generated by the Bolton (1980) |
||||
! formula for theta_e. |
||||
|
||||
! First check if pressure is less than min pressure in lookup table. |
||||
! If it is, assume parcel is so dry that the given theta-e value can |
||||
! be interpretted as theta, and get temperature from the simple dry |
||||
! theta formula. |
||||
|
||||
IF (prs .LE. psadiprs(150)) THEN |
||||
TONPSADIABAT = thte * (prs/1000.D0)**gamma |
||||
RETURN |
||||
END IF |
||||
|
||||
! Otherwise, look for the given thte/prs point in the lookup table. |
||||
|
||||
jt = -1 |
||||
DO jtch = 1, 150-1 |
||||
IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN |
||||
jt = jtch |
||||
EXIT |
||||
!GO TO 213 |
||||
END IF |
||||
END DO |
||||
|
||||
ip = -1 |
||||
DO ipch = 1, 150-1 |
||||
IF (prs .LE. psadiprs(ipch) .AND. prs .GT. psadiprs(ipch+1)) THEN |
||||
ip = ipch |
||||
EXIT |
||||
!GO TO 215 |
||||
END IF |
||||
END DO |
||||
|
||||
IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN |
||||
! Set the error and return |
||||
TONPSADIABAT = -1 |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "capecalc3d: Outside of lookup table bounds. prs,thte=", prs, thte |
||||
RETURN |
||||
END IF |
||||
|
||||
fracjt = (thte-psadithte(jt)) / (psadithte(jt+1)-psadithte(jt)) |
||||
fracjt2 = 1.D0 - fracjt |
||||
fracip = (psadiprs(ip)-prs) / (psadiprs(ip)-psadiprs(ip+1)) |
||||
fracip2 = 1.D0 - fracip |
||||
|
||||
IF (psaditmk(ip,jt) .GT. 1D9 .OR. psaditmk(ip+1,jt) .GT. 1D9 .OR. & |
||||
psaditmk(ip,jt+1) .GT. 1D9 .OR. psaditmk(ip+1,jt+1) .GT. 1D9) THEN |
||||
! Set the error and return |
||||
TONPSADIABAT = -1 |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "capecalc3d: Tried to access missing temperature in lookup table. ", & |
||||
"Prs and Thte probably unreasonable. prs,thte=", prs, thte |
||||
RETURN |
||||
END IF |
||||
|
||||
TONPSADIABAT = fracip2*fracjt2*psaditmk(ip,jt) + fracip*fracjt2*psaditmk(ip+1,jt) + & |
||||
fracip2*fracjt*psaditmk(ip,jt+1) + fracip*fracjt*psaditmk(ip+1,jt+1) |
||||
|
||||
RETURN |
||||
|
||||
END FUNCTION TONPSADIABAT |
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, fname, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR |
||||
|
||||
!f2py threadsafe |
||||
|
||||
REAL(KIND=8), DIMENSION(150), INTENT(INOUT) :: psadithte, psadiprs |
||||
REAL(KIND=8), DIMENSION(150,150), INTENT(INOUT) :: psaditmk |
||||
CHARACTER(LEN=*), INTENT(IN) :: fname |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
!NCLEND |
||||
|
||||
! Locals |
||||
INTEGER :: iustnlist, i, nthte, nprs, ip, jt |
||||
|
||||
! FNAME = 'psadilookup.dat' |
||||
iustnlist = 33 |
||||
OPEN (UNIT=iustnlist, FILE=fname, FORM='formatted', STATUS='old') |
||||
|
||||
DO i = 1,14 |
||||
READ (iustnlist, FMT=*) |
||||
END DO |
||||
|
||||
READ (iustnlist, FMT=*) nthte, nprs |
||||
|
||||
IF (nthte .NE. 150 .OR. nprs .NE. 150) THEN |
||||
errstat = ALGERR |
||||
errmsg = "Number of pressure or theta_e levels in lookup table file not 150" |
||||
RETURN |
||||
END IF |
||||
|
||||
READ (iustnlist, FMT="(5D15.7)") (psadithte(jt),jt=1,nthte) |
||||
READ (iustnlist, FMT="(5D15.7)") (psadiprs(ip),ip=1,nprs) |
||||
READ (iustnlist, FMT="(5D15.7)") ((psaditmk(ip,jt),ip=1,nprs),jt=1,nthte) |
||||
|
||||
CLOSE (iustnlist) |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DLOOKUP_TABLE |
||||
|
||||
|
||||
! Historically, this routine calculated the pressure at full sigma |
||||
! levels when RIP was specifically designed for MM4/MM5 output. |
||||
! With the new generalized RIP (Feb '02), this routine is still |
||||
! intended to calculate a set of pressure levels that bound the |
||||
! layers represented by the vertical grid points, although no such |
||||
! layer boundaries are assumed to be defined. The routine simply |
||||
! uses the midpoint between the pressures of the vertical grid |
||||
! points as the bounding levels. The array only contains mkzh |
||||
! levels, so the pressure of the top of the uppermost layer is |
||||
! actually excluded. The kth value of pf is the lower bounding |
||||
! pressure for the layer represented by kth data level. At the |
||||
! lower bounding level of the lowest model layer, it uses the |
||||
! surface pressure, unless the data set is pressure-level data, in |
||||
! which case it assumes the lower bounding pressure level is as far |
||||
! below the lowest vertical level as the upper bounding pressure |
||||
! level is above. |
||||
SUBROUTINE DPFCALC(prs, sfp, pf, miy, mjx, mkzh, ter_follow) |
||||
|
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: sfp |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: pf |
||||
INTEGER, INTENT(IN) :: ter_follow,miy,mjx,mkzh |
||||
|
||||
INTEGER :: i,j,k |
||||
|
||||
! do j=1,mjx-1 Artifact of MM5 |
||||
DO j = 1,mjx |
||||
! do i=1,miy-1 staggered grid |
||||
DO i = 1,miy |
||||
DO k = 1,mkzh |
||||
IF (k .EQ. mkzh) THEN |
||||
! terrain-following data |
||||
IF (ter_follow .EQ. 1) THEN |
||||
pf(i,j,k) = sfp(i,j) |
||||
! pressure-level data |
||||
ELSE |
||||
pf(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1)) |
||||
END IF |
||||
ELSE |
||||
pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k)) |
||||
END IF |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DPFCALC |
||||
|
||||
!====================================================================== |
||||
! |
||||
! !IROUTINE: capecalc3d -- Calculate CAPE and CIN |
||||
! |
||||
! !DESCRIPTION: |
||||
! |
||||
! If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2, |
||||
! or J/kg) for every grid point in the entire 3D domain (treating |
||||
! each grid point as a parcel). If i3dflag=0, then it |
||||
! calculates CAPE and CIN only for the parcel with max theta-e in |
||||
! the column, (i.e. something akin to Colman's MCAPE). By "parcel", |
||||
! we mean a 500-m deep parcel, with actual temperature and moisture |
||||
! averaged over that depth. |
||||
! |
||||
! In the case of i3dflag=0, |
||||
! CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of |
||||
! the cape and cin arrays. Also, if i3dflag=0, LCL and LFC heights |
||||
! are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array. |
||||
! |
||||
|
||||
|
||||
! Important! The z-indexes must be arranged so that mkzh (max z-index) is the |
||||
! surface pressure. So, pressure must be ordered in ascending order before |
||||
! calling this routine. Other variables must be ordered the same (p,tk,q,z). |
||||
|
||||
! Also, be advised that missing data values are not checked during the computation. |
||||
! Also also, Pressure must be hPa |
||||
|
||||
! NCLFORTSTART |
||||
SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& |
||||
cmsg,miy,mjx,mkzh,i3dflag,ter_follow,& |
||||
psafile, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & |
||||
EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & |
||||
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: cape, cin |
||||
|
||||
INTEGER, INTENT(IN) :: miy,mjx,mkzh,i3dflag,ter_follow |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: qvp |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: ght |
||||
REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: ter |
||||
REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) ::sfp |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: cape |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(OUT) :: cin |
||||
REAL(KIND=8), INTENT(IN) :: cmsg |
||||
CHARACTER(LEN=*), INTENT(IN) :: psafile |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
! NCLFORTEND |
||||
|
||||
! local variables |
||||
INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar, kpar1, kpar2 |
||||
REAL(KIND=8) :: davg, ethmax, q, t, p, e, eth, tlcl, zlcl |
||||
REAL(KIND=8) :: pavg, tvirtual, p1, p2, pp1, pp2, th, totthe, totqvp, totprs |
||||
REAL(KIND=8) :: cpm, deltap, ethpari, gammam, ghtpari, qvppari, prspari, tmkpari |
||||
REAL(KIND=8) :: facden, fac1, fac2, qvplift, tmklift, tvenv, tvlift, ghtlift |
||||
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat |
||||
REAL(KIND=8) :: benamin, dz, pup, pdn |
||||
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh) :: prsf |
||||
REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs |
||||
REAL(KIND=8), DIMENSION(150,150) :: psaditmk |
||||
LOGICAL :: elfound |
||||
|
||||
! To remove compiler warnings |
||||
tmkpari = 0 |
||||
qvppari = 0 |
||||
klev = 0 |
||||
klcl = 0 |
||||
kel = 0 |
||||
|
||||
|
||||
! the comments were taken from a mark stoelinga email, 23 apr 2007, |
||||
! in response to a user getting the "outside of lookup table bounds" |
||||
! error message. |
||||
|
||||
! tmkpari - initial temperature of parcel, k |
||||
! values of 300 okay. (not sure how much from this you can stray.) |
||||
|
||||
! prspari - initial pressure of parcel, hpa |
||||
! values of 980 okay. (not sure how much from this you can stray.) |
||||
|
||||
! thtecon1, thtecon2, thtecon3 |
||||
! these are all constants, the first in k and the other two have |
||||
! no units. values of 3376, 2.54, and 0.81 were stated as being |
||||
! okay. |
||||
|
||||
! tlcl - the temperature at the parcel's lifted condensation level, k |
||||
! should be a reasonable atmospheric temperature around 250-300 k |
||||
! (398 is "way too high") |
||||
|
||||
! qvppari - the initial water vapor mixing ratio of the parcel, |
||||
! kg/kg (should range from 0.000 to 0.025) |
||||
! |
||||
|
||||
! calculated the pressure at full sigma levels (a set of pressure |
||||
! levels that bound the layers represented by the vertical grid points) |
||||
|
||||
CALL DPFCALC(prs, sfp, prsf, miy, mjx, mkzh, ter_follow) |
||||
|
||||
! before looping, set lookup table for getting temperature on |
||||
! a pseudoadiabat. |
||||
|
||||
CALL DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, psafile, errstat, errmsg) |
||||
|
||||
IF (errstat .NE. 0) THEN |
||||
RETURN |
||||
END IF |
||||
|
||||
! do j=1,mjx-1 |
||||
DO j = 1,mjx |
||||
! do i=1,miy-1 |
||||
DO i = 1,miy |
||||
cape(i,j,1) = 0.d0 |
||||
cin(i,j,1) = 0.d0 |
||||
|
||||
IF (i3dflag .EQ. 1) THEN |
||||
kpar1 = 2 |
||||
kpar2 = mkzh |
||||
ELSE |
||||
! find parcel with max theta-e in lowest 3 km agl. |
||||
ethmax = -1.d0 |
||||
DO k = mkzh,1,-1 |
||||
IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN |
||||
q = MAX(qvp(i,j,k), 1.d-15) |
||||
t = tmk(i,j,k) |
||||
p = prs(i,j,k) |
||||
e = q*p/(EPS + q) |
||||
tlcl = TLCLC1 / (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4 |
||||
eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*& |
||||
EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q)) |
||||
IF (eth .GT. ethmax) THEN |
||||
klev = k |
||||
ethmax = eth |
||||
END IF |
||||
END IF |
||||
END DO |
||||
kpar1 = klev |
||||
kpar2 = klev |
||||
|
||||
! Establish average properties of that parcel |
||||
! (over depth of approximately davg meters) |
||||
|
||||
! davg=.1 |
||||
davg = 500.d0 |
||||
pavg = davg*prs(i,j,kpar1)*& |
||||
G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1))) |
||||
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh)) |
||||
p1 = p2 - pavg |
||||
totthe = 0.d0 |
||||
totqvp = 0.d0 |
||||
totprs = 0.d0 |
||||
DO k = mkzh,2,-1 |
||||
IF (prsf(i,j,k) .LE. p1) EXIT !GOTO 35 |
||||
IF (prsf(i,j,k-1) .GE. p2) CYCLE !GOTO 34 |
||||
p = prs(i,j,k) |
||||
pup = prsf(i,j,k) |
||||
pdn = prsf(i,j,k-1) |
||||
q = MAX(qvp(i,j,k),1.d-15) |
||||
th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q)) |
||||
pp1 = MAX(p1,pdn) |
||||
pp2 = MIN(p2,pup) |
||||
IF (pp2 .GT. pp1) THEN |
||||
deltap = pp2 - pp1 |
||||
totqvp = totqvp + q*deltap |
||||
totthe = totthe + th*deltap |
||||
totprs = totprs + deltap |
||||
END IF |
||||
! 34 CONTINUE |
||||
END DO |
||||
! 35 CONTINUE |
||||
qvppari = totqvp/totprs |
||||
tmkpari = (totthe/totprs)*& |
||||
(prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1))) |
||||
END IF |
||||
|
||||
DO kpar = kpar1, kpar2 |
||||
|
||||
! Calculate temperature and moisture properties of parcel |
||||
! (note, qvppari and tmkpari already calculated above for 2d case.) |
||||
|
||||
IF (i3dflag .EQ. 1) THEN |
||||
qvppari = qvp(i,j,kpar) |
||||
tmkpari = tmk(i,j,kpar) |
||||
END IF |
||||
prspari = prs(i,j,kpar) |
||||
ghtpari = ght(i,j,kpar) |
||||
gammam = GAMMA * (1.d0 + GAMMAMD*qvppari) |
||||
cpm = CP * (1.d0 + CPMD*qvppari) |
||||
|
||||
e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari)) |
||||
tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 |
||||
ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*& |
||||
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari)) |
||||
zlcl = ghtpari + (tmkpari - tlcl)/(G/cpm) |
||||
|
||||
! Calculate buoyancy and relative height of lifted parcel at |
||||
! all levels, and store in bottom up arrays. add a level at the lcl, |
||||
! and at all points where buoyancy is zero. |
||||
! |
||||
! For arrays that go bottom to top |
||||
kk = 0 |
||||
ilcl = 0 |
||||
|
||||
IF (ghtpari .GE. zlcl) THEN |
||||
! Initial parcel already saturated or supersaturated. |
||||
ilcl = 2 |
||||
klcl = 1 |
||||
END IF |
||||
|
||||
k = kpar |
||||
DO WHILE (k .GE. 1)!k = kpar, 1, -1 |
||||
!DO k = kpar, 1, -1 |
||||
! For arrays that go bottom to top |
||||
! 33 kk = kk + 1 |
||||
kk = kk + 1 |
||||
|
||||
! Model level is below lcl |
||||
IF (ght(i,j,k) .LT. zlcl) THEN |
||||
qvplift = qvppari |
||||
tmklift = tmkpari - G/cpm*(ght(i,j,k) - ghtpari) |
||||
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) |
||||
tvlift = tvirtual(tmklift, qvplift) |
||||
ghtlift = ght(i,j,k) |
||||
ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN |
||||
! This model level and previous model level straddle the lcl, |
||||
! so first create a new level in the bottom-up array, at the lcl. |
||||
tmklift = tlcl |
||||
qvplift = qvppari |
||||
facden = ght(i,j,k) - ght(i,j,k+1) |
||||
fac1 = (zlcl-ght(i,j,k+1))/facden |
||||
fac2 = (ght(i,j,k)-zlcl)/facden |
||||
tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1 |
||||
qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1 |
||||
tvenv = tvirtual(tmkenv, qvpenv) |
||||
tvlift = tvirtual(tmklift, qvplift) |
||||
ghtlift = zlcl |
||||
ilcl = 1 |
||||
ELSE |
||||
tmklift = TONPSADIABAT(ethpari, prs(i,j,k), psadithte, psadiprs,& |
||||
psaditmk, GAMMA, errstat, errmsg) |
||||
eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) |
||||
qvplift = EPS*eslift/(prs(i,j,k) - eslift) |
||||
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) |
||||
tvlift = tvirtual(tmklift, qvplift) |
||||
ghtlift = ght(i,j,k) |
||||
END IF |
||||
! Buoyancy |
||||
buoy(kk) = G*(tvlift - tvenv)/tvenv |
||||
zrel(kk) = ghtlift - ghtpari |
||||
|
||||
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0d0)) THEN |
||||
! Parcel ascent curve crosses sounding curve, so create a new level |
||||
! in the bottom-up array at the crossing. |
||||
kk = kk + 1 |
||||
buoy(kk) = buoy(kk-1) |
||||
zrel(kk) = zrel(kk-1) |
||||
buoy(kk-1) = 0.d0 |
||||
zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& |
||||
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) |
||||
END IF |
||||
|
||||
IF (ilcl .EQ. 1) THEN |
||||
klcl = kk |
||||
ilcl = 2 |
||||
!GOTO 33 |
||||
CYCLE |
||||
END IF |
||||
|
||||
k = k - 1 |
||||
END DO |
||||
|
||||
kmax = kk |
||||
IF (kmax .GT. 150) THEN |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax |
||||
RETURN |
||||
END IF |
||||
|
||||
! If no lcl was found, set klcl to kmax. it is probably not really |
||||
! at kmax, but this will make the rest of the routine behave |
||||
! properly. |
||||
IF (ilcl .EQ. 0) klcl=kmax |
||||
|
||||
! Get the accumulated buoyant energy from the parcel's starting |
||||
! point, at all levels up to the top level. |
||||
benaccum(1) = 0.0d0 |
||||
benamin = 9d9 |
||||
DO k = 2,kmax |
||||
dz = zrel(k) - zrel(k-1) |
||||
benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k)) |
||||
IF (benaccum(k) .LT. benamin) THEN |
||||
benamin = benaccum(k) |
||||
END IF |
||||
END DO |
||||
|
||||
! Determine equilibrium level (el), which we define as the highest |
||||
! level of non-negative buoyancy above the lcl. note, this may be |
||||
! the top level if the parcel is still buoyant there. |
||||
|
||||
elfound = .FALSE. |
||||
DO k = kmax,klcl,-1 |
||||
IF (buoy(k) .GE. 0.d0) THEN |
||||
! k of equilibrium level |
||||
kel = k |
||||
elfound = .TRUE. |
||||
EXIT |
||||
!GOTO 50 |
||||
END IF |
||||
END DO |
||||
|
||||
! If we got through that loop, then there is no non-negative |
||||
! buoyancy above the lcl in the sounding. in these situations, |
||||
! both cape and cin will be set to -0.1 j/kg. (see below about |
||||
! missing values in v6.1.0). also, where cape is |
||||
! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so |
||||
! that the zero contour in either the cin or cape fields will |
||||
! circumscribe regions of non-zero cape. |
||||
|
||||
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
||||
! value of this function. at that time we decided to change -0.1 |
||||
! to a more appropriate missing value, which is passed into this |
||||
! routine as cmsg. |
||||
|
||||
! cape(i,j,kpar) = -0.1d0 |
||||
! cin(i,j,kpar) = -0.1d0 |
||||
IF (.NOT. elfound) THEN |
||||
cape(i,j,kpar) = cmsg |
||||
cin(i,j,kpar) = cmsg |
||||
klfc = kmax |
||||
CYCLE |
||||
END IF |
||||
|
||||
! GOTO 102 |
||||
|
||||
! 50 CONTINUE |
||||
|
||||
! If there is an equilibrium level, then cape is positive. we'll |
||||
! define the level of free convection (lfc) as the point below the |
||||
! el, but at or above the lcl, where accumulated buoyant energy is a |
||||
! minimum. the net positive area (accumulated buoyant energy) from |
||||
! the lfc up to the el will be defined as the cape, and the net |
||||
! negative area (negative of accumulated buoyant energy) from the |
||||
! parcel starting point to the lfc will be defined as the convective |
||||
! inhibition (cin). |
||||
|
||||
! First get the lfc according to the above definition. |
||||
benamin = 9d9 |
||||
klfc = kmax |
||||
DO k = klcl,kel |
||||
IF (benaccum(k) .LT. benamin) THEN |
||||
benamin = benaccum(k) |
||||
klfc = k |
||||
END IF |
||||
END DO |
||||
|
||||
! Now we can assign values to cape and cin |
||||
|
||||
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0) |
||||
cin(i,j,kpar) = MAX(-benamin, 0.1d0) |
||||
|
||||
! cin is uninteresting when cape is small (< 100 j/kg), so set |
||||
! cin to -0.1 (see note about missing values in v6.1.0) in |
||||
! that case. |
||||
|
||||
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
||||
! value of this function. at that time we decided to change -0.1 |
||||
! to a more appropriate missing value, which is passed into this |
||||
! routine as cmsg. |
||||
|
||||
! IF (cape(i,j,kpar).lt.100.d0) cin(i,j,kpar) = -0.1d0 |
||||
IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg |
||||
! 102 CONTINUE |
||||
|
||||
END DO |
||||
|
||||
IF (i3dflag .EQ. 0) THEN |
||||
cape(i,j,mkzh) = cape(i,j,kpar1) |
||||
cin(i,j,mkzh) = cin(i,j,kpar1) |
||||
! meters agl |
||||
cin(i,j,mkzh-1) = zrel(klcl) + ghtpari - ter(i,j) |
||||
! meters agl |
||||
cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j) |
||||
|
||||
ENDIF |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
END SUBROUTINE DCAPECALC3D |
@ -0,0 +1,125 @@
@@ -0,0 +1,125 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz) |
||||
USE constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: ctt |
||||
|
||||
INTEGER, INTENT(IN) :: nz, ns, ew, haveqci |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp |
||||
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter |
||||
REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt |
||||
|
||||
!NCLEND |
||||
|
||||
! REAL(KIND=8) :: znfac(nz) |
||||
|
||||
! LOCAL VARIABLES |
||||
INTEGER i,j,k,ripk |
||||
!INTEGER :: mjx,miy,mkzh |
||||
REAL(KIND=8) :: vt,opdepthu,opdepthd,dp |
||||
REAL(KIND=8) :: ratmix,arg1,arg2,agl_hgt |
||||
REAL(KIND=8) :: fac,prsctt |
||||
!REAL(KIND=8) :: eps,ussalr,rgas,grav,abscoefi,abscoef,celkel,wrfout |
||||
!REAL(KIND=8) :: ght(ew,ns,nz),stuff(ew,ns) |
||||
!REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf(ns,ew,nz),p1,p2 |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf |
||||
REAL(KIND=8) :: p1, p2 |
||||
|
||||
!mjx = ew |
||||
!miy = ns |
||||
!mkzh = nz |
||||
|
||||
prsctt = 0 ! removes the warning |
||||
|
||||
! Calculate the surface pressure |
||||
DO j=1,ns |
||||
DO i=1,ew |
||||
ratmix = .001d0*qvp(i,j,1) |
||||
arg1 = EPS + ratmix |
||||
arg2 = EPS * (1. + ratmix) |
||||
vt = tk(i,j,1) * arg1/arg2 !Virtual temperature |
||||
agl_hgt = ght(i,j,nz) - ter(i,j) |
||||
arg1 = -G / (RD * USSALR) |
||||
pf(i,j,nz) = prs(i,j,1) * (vt / (vt + USSALR*(agl_hgt)))**(arg1) |
||||
END DO |
||||
END DO |
||||
|
||||
DO k=1,nz-1 |
||||
DO j=1,ns |
||||
DO i=1,ew |
||||
ripk = nz-k+1 |
||||
pf(i,j,k) = .5d0 * (prs(i,j,ripk) + prs(i,j,ripk-1)) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
DO j=1,ns |
||||
DO i=1,ew |
||||
opdepthd = 0.d0 |
||||
k = 0 |
||||
|
||||
! Integrate downward from model top, calculating path at full |
||||
! model vertical levels. |
||||
|
||||
!20 opdepthu=opdepthd |
||||
|
||||
DO k=1, nz |
||||
opdepthu = opdepthd |
||||
!k=k+1 |
||||
ripk = nz - k + 1 |
||||
|
||||
IF (k .EQ. 1) THEN |
||||
dp = 200.d0 * (pf(i,j,1) - prs(i,j,nz)) ! should be in Pa |
||||
ELSE |
||||
dp = 100.d0 * (pf(i,j,k) - pf(i,j,k-1)) ! should be in Pa |
||||
END IF |
||||
|
||||
IF (haveqci .EQ. 0) then |
||||
IF (tk(i,j,k) .LT. CELKEL) then |
||||
! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed |
||||
opdepthd = opdepthu + ABSCOEFI*qcw(i,j,k) * dp/G |
||||
ELSE |
||||
opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G |
||||
END IF |
||||
ELSE |
||||
opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk)) * dp/G |
||||
END IF |
||||
|
||||
IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN |
||||
!GOTO 20 |
||||
CYCLE |
||||
|
||||
ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN |
||||
prsctt = prs(i,j,1) |
||||
EXIT |
||||
ELSE |
||||
fac = (1. - opdepthu) / (opdepthd - opdepthu) |
||||
prsctt = pf(i,j,k-1) + fac*(pf(i,j,k) - pf(i,j,k-1)) |
||||
prsctt = MIN(prs(i,j,1), MAX(prs(i,j,nz), prsctt)) |
||||
EXIT |
||||
END IF |
||||
END DO |
||||
|
||||
DO k=2,nz |
||||
ripk = nz-k+1 |
||||
p1 = prs(i,j,ripk+1) |
||||
p2 = prs(i,j,ripk) |
||||
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN |
||||
fac = (prsctt - p1) / (p2 - p1) |
||||
arg1 = fac * (tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL |
||||
ctt(i,j) = tk(i,j,ripk+1) + arg1 |
||||
!GOTO 40 |
||||
EXIT |
||||
END IF |
||||
END DO |
||||
END DO |
||||
END DO |
||||
! 30 CONTINUE |
||||
! 40 CONTINUE |
||||
! 190 CONTINUE |
||||
RETURN |
||||
|
||||
END SUBROUTINE wrfcttcalc |
@ -0,0 +1,122 @@
@@ -0,0 +1,122 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DCOMPUTEABSVORT(av, u, v, msfu, msfv, msft, cor, dx, dy, nx, ny, nz,& |
||||
nxp1, nyp1) |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: av |
||||
|
||||
INTEGER, INTENT(IN) :: nx, ny, nz, nxp1, nyp1 |
||||
REAL(KIND=8), DIMENSION(nxp1,ny,nz), INTENT(IN) :: u |
||||
REAL(KIND=8), DIMENSION(nx,nyp1,nz), INTENT(IN) :: v |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: av |
||||
REAL(KIND=8), DIMENSION(nxp1,ny), INTENT(IN):: msfu |
||||
REAL(KIND=8), DIMENSION(nx,nyp1), INTENT(IN) :: msfv |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: msft |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: cor |
||||
REAL(KIND=8) :: dx, dy |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: jp1, jm1, ip1, im1, i, j, k |
||||
REAL(KIND=8) :: dsy, dsx, dudy, dvdx, avort |
||||
REAL(KIND=8) :: mm |
||||
|
||||
! PRINT*,'nx,ny,nz,nxp1,nyp1' |
||||
! PRINT*,nx,ny,nz,nxp1,nyp1 |
||||
DO k = 1,nz |
||||
DO j = 1,ny |
||||
jp1 = MIN(j+1, ny) |
||||
jm1 = MAX(j-1, 1) |
||||
DO i = 1,nx |
||||
ip1 = MIN(i+1, nx) |
||||
im1 = MAX(i-1, 1) |
||||
! PRINT *,jp1,jm1,ip1,im1 |
||||
dsx = (ip1 - im1) * dx |
||||
dsy = (jp1 - jm1) * dy |
||||
mm = msft(i,j)*msft(i,j) |
||||
! PRINT *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1) |
||||
dudy = 0.5D0 * (u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & |
||||
u(i,jm1,k)/msfu(i,jm1) - u(i+1,jm1,k)/msfu(i+1,jm1))/dsy*mm |
||||
dvdx = 0.5D0 * (v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & |
||||
v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm |
||||
avort = dvdx - dudy + cor(i,j) |
||||
av(i,j,k) = avort*1.D5 |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCOMPUTEABSVORT |
||||
|
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DCOMPUTEPV(pv, u, v, theta, prs, msfu, msfv, msft, cor, dx, dy, nx, & |
||||
ny, nz, nxp1, nyp1) |
||||
USE constants, ONLY : G |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: pv |
||||
|
||||
INTEGER,INTENT(IN) :: nx, ny, nz, nxp1, nyp1 |
||||
REAL(KIND=8), DIMENSION(nxp1,ny,nz), INTENT(IN) :: u |
||||
REAL(KIND=8), DIMENSION(nx,nyp1,nz), INTENT(IN) :: v |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: theta |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: pv |
||||
REAL(KIND=8), DIMENSION(nxp1,ny), INTENT(IN) :: msfu |
||||
REAL(KIND=8), DIMENSION(nx,nyp1), INTENT(IN) :: msfv |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: msft |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: cor |
||||
REAL(KIND=8) :: dx,dy |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: kp1, km1, jp1, jm1, ip1, im1, i, j, k |
||||
REAL(KIND=8) :: dsy, dsx, dp, dudy, dvdx, dudp, dvdp, dthdp, avort |
||||
REAL(KIND=8) :: dthdx, dthdy, mm |
||||
|
||||
! PRINT*,'nx,ny,nz,nxp1,nyp1' |
||||
! PRINT*,nx,ny,nz,nxp1,nyp1 |
||||
DO k = 1,nz |
||||
kp1 = MIN(k+1, nz) |
||||
km1 = MAX(k-1, 1) |
||||
DO J = 1,ny |
||||
jp1 = MIN(j+1, ny) |
||||
jm1 = MAX(j-1, 1) |
||||
DO i = 1,nx |
||||
ip1 = MIN(i+1, nx) |
||||
im1 = MAX(i-1, 1) |
||||
! PRINT *,jp1,jm1,ip1,im1 |
||||
dsx = (ip1 - im1) * dx |
||||
dsy = (jp1 - jm1) * dy |
||||
mm = msft(i,j)*msft(i,j) |
||||
! PRINT *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1) |
||||
dudy = 0.5D0 * (u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & |
||||
u(i,jm1,k)/msfu(i,jm1) - u(i+1,jm1,k)/msfu(i+1,jm1))/dsy*mm |
||||
dvdx = 0.5D0 * (v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & |
||||
v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm |
||||
avort = dvdx - dudy + cor(i,j) |
||||
dp = prs(i,j,kp1) - prs(i,j,km1) |
||||
dudp = 0.5D0 * (u(i,j,kp1) + u(i+1,j,kp1) - u(i,j,km1) - u(i+1,j,km1))/dp |
||||
dvdp = 0.5D0 * (v(i,j,kp1) + v(i,j+1,kp1) - v(i,j,km1) - v(i,J+1,km1))/dp |
||||
dthdp = (theta(i,j,kp1) - theta(i,j,km1))/dp |
||||
dthdx = (theta(ip1,j,k) - theta(im1,j,k))/dsx * msft(i,j) |
||||
dthdy = (theta(i,jp1,k) - theta(i,jm1,k))/dsy * msft(i,j) |
||||
pv(i,j,k) = -G * (dthdp*avort - dvdp*dthdx + dudp*dthdy)*10000.D0 |
||||
! if(i.eq.300 .and. j.eq.300) then |
||||
! PRINT*,'avort,dudp,dvdp,dthdp,dthdx,dthdy,pv' |
||||
! PRINT*,avort,dudp,dvdp,dthdp,dthdx,dthdy,pv(i,j,k) |
||||
! endif |
||||
pv(i,j,k) = pv(i,j,k)*1.D2 |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCOMPUTEPV |
@ -0,0 +1,31 @@
@@ -0,0 +1,31 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) |
||||
USE constants, ONLY : RD |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: pw |
||||
|
||||
INTEGER, INTENT(IN) :: nx, ny, nz, nzh |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv |
||||
REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht |
||||
REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pw |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: i, j, k |
||||
!REAL(KIND=8),PARAMETER :: R=287.06 |
||||
|
||||
pw = 0 |
||||
DO k=1,nz |
||||
DO j=1,ny |
||||
DO i=1,nx |
||||
pw(i,j) = pw(i,j) + ((p(i,j,k)/(RD*tv(i,j,k))) * qv(i,j,k) * (ht(i,j,k+1) - ht(i,j,k))) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCOMPUTEPW |
@ -0,0 +1,84 @@
@@ -0,0 +1,84 @@
|
||||
! NCLFORTSTART |
||||
SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) |
||||
USE constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: sreh |
||||
|
||||
INTEGER, INTENT(IN) :: miy, mjx, mkzh |
||||
REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: u, v, ght |
||||
REAL(KIND=8), INTENT(IN) :: top |
||||
REAL(KIND=8), DIMENSION(miy,mjx), INTENT(IN) :: ter |
||||
REAL(KIND=8), DIMENSION(miy,mjx), INTENT(OUT) :: sreh |
||||
|
||||
! NCLEND |
||||
|
||||
! This helicity code was provided by Dr. Craig Mattocks, and |
||||
! verified by Cindy Bruyere to produce results equivalent to |
||||
! those generated by RIP4. (The code came from RIP4?) |
||||
|
||||
REAL(KIND=8) :: dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr |
||||
REAL(KIND=8) :: cu, cv, x, sum |
||||
INTEGER :: i, j, k, k10, k3, ktop |
||||
!REAL(KIND=8), PARAMETER :: DTR=PI/180.d0, DPR=180.d0/PI |
||||
|
||||
DO j = 1, mjx-1 |
||||
DO i = 1, miy-1 |
||||
sdh = 0.d0 |
||||
su = 0.d0 |
||||
sv = 0.d0 |
||||
k3 = 0 |
||||
k10 = 0 |
||||
ktop = 0 |
||||
DO k = mkzh, 2, -1 |
||||
IF (((ght(i,j,k) - ter(i,j)) .GT. 10000.D0) .AND. (k10 .EQ. 0)) THEN |
||||
k10 = k |
||||
EXIT |
||||
ENDIF |
||||
IF (((ght(i,j,k) - ter(i,j)) .GT. top) .AND. (ktop .EQ. 0)) THEN |
||||
ktop = k |
||||
ENDIF |
||||
IF (((ght(i,j,k) - ter(i,j)) .GT. 3000.D0) .AND. (k3 .EQ. 0)) THEN |
||||
k3 = k |
||||
ENDIF |
||||
END DO |
||||
|
||||
IF (k10 .EQ. 0) THEN |
||||
k10 = 2 |
||||
ENDIF |
||||
DO k = k3, k10, -1 |
||||
dh = ght(i,j,k-1) - ght(i,j,k) |
||||
sdh = sdh + dh |
||||
su = su + 0.5D0*dh*(u(i,j,k-1) + u(i,j,k)) |
||||
sv = sv + 0.5D0*dh*(v(i,j,k-1) + v(i,j,k)) |
||||
END DO |
||||
ua = su / sdh |
||||
va = sv / sdh |
||||
asp = SQRT(ua*ua + va*va) |
||||
IF (ua .EQ. 0.D0 .AND. va .EQ. 0.D0) THEN |
||||
adr = 0.D0 |
||||
ELSE |
||||
adr = DEG_PER_RAD * (PI + ATAN2(ua,va)) |
||||
ENDIF |
||||
bsp = 0.75D0 * asp |
||||
bdr = adr + 30.D0 |
||||
IF (bdr .GT. 360.D0) THEN |
||||
bdr = bdr - 360.D0 |
||||
ENDIF |
||||
cu = -bsp * SIN(bdr * RAD_PER_DEG) |
||||
cv = -bsp * COS(bdr * RAD_PER_DEG) |
||||
sum = 0.D0 |
||||
DO k = mkzh-1, ktop, -1 |
||||
x = ((u(i,j,k) - cu) * (v(i,j,k) - v(i,j,k+1))) - & |
||||
((v(i,j,k) - cv) * (u(i,j,k) - u(i,j,k+1))) |
||||
sum = sum + x |
||||
END DO |
||||
sreh(i,j) = -sum |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DCALRELHL |
@ -0,0 +1,269 @@
@@ -0,0 +1,269 @@
|
||||
!====================================================================== |
||||
! |
||||
! !IROUTINE: WETBULBCALC -- Calculate wet bulb temperature (C) |
||||
! |
||||
! !DESCRIPTION: |
||||
! |
||||
! Calculates wet bulb temperature in C, given pressure in |
||||
! temperature in K and mixing ratio in kg/kg. |
||||
! |
||||
! !INPUT: |
||||
! nx - index for x dimension |
||||
! ny - index for y dimension |
||||
! nz - index for z dimension |
||||
! prs - pressure (mb) |
||||
! tmk - temperature (K) |
||||
! qvp - water vapor mixing ratio (kg/kg) |
||||
! |
||||
! !OUTPUT: |
||||
! twb - Wet bulb temperature (C) |
||||
! |
||||
! !ASSUMPTIONS: |
||||
! |
||||
! !REVISION HISTORY: |
||||
! 2009-March - Mark T. Stoelinga - from RIP4.5 |
||||
! 2010-August - J. Schramm |
||||
! 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
||||
! |
||||
! !INTERFACE: |
||||
! ------------------------------------------------------------------ |
||||
|
||||
! NCLFORTSTART |
||||
SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, & |
||||
EPS, TLCLC4, THTECON1, THTECON2, THTECON3 |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: twb |
||||
|
||||
INTEGER, INTENT(IN) :: nx, ny, nz |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: tmk |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: qvp |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: twb |
||||
CHARACTER(LEN=*), INTENT(IN) :: psafile |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: i, j, k |
||||
INTEGER :: jtch, jt, ipch, ip |
||||
REAL(KIND=8) :: q, t, p, e, tlcl, eth |
||||
REAL(KIND=8) :: fracip, fracip2, fracjt, fracjt2 |
||||
REAL(KIND=8), DIMENSION(150) :: PSADITHTE, PSADIPRS |
||||
REAL(KIND=8), DIMENSION(150,150) :: PSADITMK |
||||
REAL(KIND=8) :: tonpsadiabat |
||||
|
||||
|
||||
! Before looping, set lookup table for getting temperature on |
||||
! a pseudoadiabat. |
||||
|
||||
CALL DLOOKUP_TABLE(PSADITHTE, PSADIPRS, PSADITMK, psafile, errstat, errmsg) |
||||
|
||||
IF (errstat .NE. 0) THEN |
||||
RETURN |
||||
END IF |
||||
|
||||
DO k=1,nz |
||||
DO j=1,ny |
||||
DO i=1,nx |
||||
q = DMAX1(qvp(i,j,k), 1.D-15) |
||||
t = tmk(i,j,k) |
||||
p = prs(i,j,k)/100. |
||||
e = q*p/(EPS + q) |
||||
tlcl = TLCLC1/(DLOG(t**TLCLC2/e) - TLCLC3) + TLCLC4 |
||||
eth = t*(1000./p)**(GAMMA*(1. + GAMMAMD*q))*& |
||||
EXP((THTECON1/tlcl - THTECON2)*q*(1. + THTECON3*q)) |
||||
|
||||
! Now we need to find the temperature (in K) on a moist adiabat |
||||
! (specified by eth in K) given pressure in hPa. It uses a |
||||
! lookup table, with data that was generated by the Bolton (1980) |
||||
! formula for theta_e. |
||||
|
||||
! First check if pressure is less than min pressure in lookup table. |
||||
! If it is, assume parcel is so dry that the given theta-e value can |
||||
! be interpretted as theta, and get temperature from the simple dry |
||||
! theta formula. |
||||
|
||||
IF (p .LE. PSADIPRS(150)) THEN |
||||
tonpsadiabat = eth*(p/1000.)**GAMMA |
||||
ELSE |
||||
! Otherwise, look for the given thte/prs point in the lookup table. |
||||
jt=-1 |
||||
DO jtch=1,150-1 |
||||
IF (eth .GE. PSADITHTE(jtch) .AND. eth .LT. PSADITHTE(jtch+1)) THEN |
||||
jt = jtch |
||||
EXIT |
||||
ENDIF |
||||
END DO |
||||
|
||||
ip=-1 |
||||
DO ipch=1,150-1 |
||||
IF (p .LE. PSADIPRS(ipch) .AND. p .GT. PSADIPRS(ipch+1)) THEN |
||||
ip = ipch |
||||
EXIT |
||||
ENDIF |
||||
END DO |
||||
|
||||
IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "Outside of lookup table bounds. prs,thte=", p, eth |
||||
RETURN |
||||
ENDIF |
||||
|
||||
fracjt = (eth - PSADITHTE(jt))/(PSADITHTE(jt+1) - PSADITHTE(jt)) |
||||
fracjt2 = 1. - fracjt |
||||
fracip = (PSADIPRS(ip) - p)/(PSADIPRS(ip) - PSADIPRS(ip+1)) |
||||
fracip2 = 1. - fracip |
||||
|
||||
IF (PSADITMK(ip,jt) .GT. 1e9 .OR. PSADITMK(ip+1,jt) .GT. 1e9 .OR. & |
||||
PSADITMK(ip,jt+1) .GT. 1e9 .OR. PSADITMK(ip+1,jt+1) .GT. 1e9) THEN |
||||
!PRINT*,'Tried to access missing tmperature in lookup table.' |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "Prs and Thte probably unreasonable. prs, thte=", p, eth |
||||
RETURN |
||||
ENDIF |
||||
|
||||
tonpsadiabat = fracip2*fracjt2*PSADITMK(ip,jt) + & |
||||
fracip*fracjt2*PSADITMK(ip+1,jt) + & |
||||
fracip2*fracjt*PSADITMK(ip,jt+1) + & |
||||
fracip*fracjt*PSADITMK(ip+1,jt+1) |
||||
ENDIF |
||||
|
||||
twb(i,j,k) = tonpsadiabat |
||||
|
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE WETBULBCALC |
||||
|
||||
|
||||
! !IROUTINE: omgcalc -- Calculate omega (dp/dt) |
||||
! |
||||
! !DESCRIPTION: |
||||
! |
||||
! Calculate approximate omega, based on vertical velocity w (dz/dt). |
||||
! It is approximate because it cannot take into account the vertical |
||||
! motion of pressure surfaces. |
||||
! |
||||
! !INPUT: |
||||
! mx - index for x dimension |
||||
! my - index for y dimension |
||||
! mx - index for vertical dimension |
||||
! qvp - water vapor mixing ratio (kg/kg) |
||||
! tmk - temperature (K) |
||||
! www - vertical velocity (m/s) |
||||
! prs - pressure (Pa) |
||||
! |
||||
! !OUTPUT: |
||||
! omg - omega (Pa/sec) |
||||
! |
||||
! !ASSUMPTIONS: |
||||
! |
||||
! !REVISION HISTORY: |
||||
! 2009-March - Mark T. Stoelinga - from RIP4.5 |
||||
! 2010-August - J. Schramm |
||||
! 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
||||
! |
||||
! ------------------------------------------------------------------ |
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE OMGCALC(qvp, tmk, www, prs, omg, mx, my, mz) |
||||
USE constants, ONLY : G, RD, EPS |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: omg |
||||
|
||||
INTEGER,INTENT(IN) :: mx, my, mz |
||||
REAL(KIND=8), INTENT(IN), DIMENSION(mx,my,mz) :: qvp |
||||
REAL(KIND=8), INTENT(IN), DIMENSION(mx,my,mz) :: tmk |
||||
REAL(KIND=8), INTENT(IN), DIMENSION(mx,my,mz) :: www |
||||
REAL(KIND=8), INTENT(IN), DIMENSION(mx,my,mz) :: prs |
||||
REAL(KIND=8), INTENT(OUT), DIMENSION(mx,my,mz) :: omg |
||||
|
||||
!NCLEND |
||||
|
||||
! Local variables |
||||
INTEGER :: i, j, k |
||||
!REAL(KIND=8), PARAMETER :: GRAV=9.81, RGAS=287.04, EPS=0.622 |
||||
|
||||
DO k=1,mz |
||||
DO j=1,my |
||||
DO i=1,mx |
||||
omg(i,j,k) = -G*prs(i,j,k)/& |
||||
(RD*((tmk(i,j,k)*(EPS + qvp(i,j,k)))/& |
||||
(EPS*(1. + qvp(i,j,k)))))*www(i,j,k) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE OMGCALC |
||||
|
||||
!====================================================================== |
||||
! |
||||
! !IROUTINE: VIRTUAL_TEMP -- Calculate virtual temperature (K) |
||||
! |
||||
! !DESCRIPTION: |
||||
! |
||||
! Calculates virtual temperature in K, given temperature |
||||
! in K and mixing ratio in kg/kg. |
||||
! |
||||
! !INPUT: |
||||
! NX - index for x dimension |
||||
! NY - index for y dimension |
||||
! NZ - index for z dimension |
||||
! RATMIX - water vapor mixing ratio (kg/kg) |
||||
! TEMP - temperature (K) |
||||
! |
||||
! !OUTPUT: |
||||
! TV - Virtual temperature (K) |
||||
! |
||||
! !ASSUMPTIONS: |
||||
! |
||||
! !REVISION HISTORY: |
||||
! 2009-March - Mark T. Stoelinga - from RIP4.5 |
||||
! 2010-August - J. Schramm |
||||
! 2014-March - A. Jaye - modified to run with NCL and ARW wrf output |
||||
! |
||||
! ------------------------------------------------------------------ |
||||
!NCLFORTSTART |
||||
SUBROUTINE VIRTUAL_TEMP(temp, ratmix, tv, nx, ny, nz) |
||||
USE constants, ONLY : EPS |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: tv |
||||
|
||||
INTEGER, INTENT(IN) :: nx, ny, nz |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: temp |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: ratmix |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: tv |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: i,j,k |
||||
!REAL(KIND=8),PARAMETER :: EPS = 0.622D0 |
||||
|
||||
DO k=1,nz |
||||
DO j=1,ny |
||||
DO i=1,nx |
||||
tv(i,j,k) = temp(i,j,k)*(EPS + ratmix(i,j,k))/(EPS*(1.D0 + ratmix(i,j,k))) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE VIRTUAL_TEMP |
||||
|
@ -0,0 +1,147 @@
@@ -0,0 +1,147 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx, ny, nz) |
||||
USE constants, ONLY : GAMMA_SEVEN, RHOWAT, RHO_R, RHO_S, RHO_G, ALPHA, & |
||||
CELKEL, PI, RD |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: dbz |
||||
|
||||
! Arguments |
||||
INTEGER, INTENT(IN) :: nx, ny, nz |
||||
INTEGER, INTENT(IN) :: sn0, ivarint, iliqskin |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: dbz |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: prs |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: tmk |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qvp |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qra |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qsn |
||||
REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(INOUT) :: qgr |
||||
|
||||
!NCLEND |
||||
|
||||
! Local Variables |
||||
INTEGER :: i, j, k |
||||
REAL(KIND=8) :: temp_c, virtual_t |
||||
REAL(KIND=8) :: gonv, ronv, sonv |
||||
REAL(KIND=8) :: factor_g, factor_r, factor_s |
||||
REAL(KIND=8) :: factorb_g, factorb_s |
||||
REAL(KIND=8) :: rhoair, z_e |
||||
|
||||
! Constants used to calculate variable intercepts |
||||
REAL(KIND=8), PARAMETER :: R1 = 1.D-15 |
||||
REAL(KIND=8), PARAMETER :: RON = 8.D6 |
||||
REAL(KIND=8), PARAMETER :: RON2 = 1.D10 |
||||
REAL(KIND=8), PARAMETER :: SON = 2.D7 |
||||
REAL(KIND=8), PARAMETER :: GON = 5.D7 |
||||
REAL(KIND=8), PARAMETER :: RON_MIN = 8.D6 |
||||
REAL(KIND=8), PARAMETER :: RON_QR0 = 0.00010D0 |
||||
REAL(KIND=8), PARAMETER :: RON_DELQR0 = 0.25D0*RON_QR0 |
||||
REAL(KIND=8), PARAMETER :: RON_CONST1R = (RON2-RON_MIN)*0.5D0 |
||||
REAL(KIND=8), PARAMETER :: RON_CONST2R = (RON2+RON_MIN)*0.5D0 |
||||
|
||||
! Constant intercepts |
||||
REAL(KIND=8), PARAMETER :: RN0_R = 8.D6 |
||||
REAL(KIND=8), PARAMETER :: RN0_S = 2.D7 |
||||
REAL(KIND=8), PARAMETER :: RN0_G = 4.D6 |
||||
|
||||
! Force all Q arrays to be 0.0 or greater. |
||||
DO k = 1,nz |
||||
DO j = 1,ny |
||||
DO i = 1,nx |
||||
IF (qvp(i,j,k).LT.0.0) THEN |
||||
qvp(i,j,k) = 0.0 |
||||
END IF |
||||
IF (qra(i,j,k).LT.0.0) THEN |
||||
qra(i,j,k) = 0.0 |
||||
END IF |
||||
IF (qsn(i,j,k).LT.0.0) THEN |
||||
qsn(i,j,k) = 0.0 |
||||
END IF |
||||
IF (qgr(i,j,k).LT.0.0) THEN |
||||
qgr(i,j,k) = 0.0 |
||||
END IF |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
! Input pressure is Pa, but we need hPa in calculations |
||||
|
||||
IF (sn0 .EQ. 0) THEN |
||||
DO k = 1,nz |
||||
DO j = 1,ny |
||||
DO i = 1,nx |
||||
IF (tmk(i,j,k) .LT. CELKEL) THEN |
||||
qsn(i,j,k) = qra(i,j,k) |
||||
qra(i,j,k) = 0.D0 |
||||
END IF |
||||
END DO |
||||
END DO |
||||
END DO |
||||
END IF |
||||
|
||||
factor_r = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_R))**1.75D0 |
||||
factor_s = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_S))**1.75D0*(RHO_S/RHOWAT)**2*ALPHA |
||||
factor_g = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_G))**1.75D0*(RHO_G/RHOWAT)**2*ALPHA |
||||
|
||||
DO k = 1,nz |
||||
DO j = 1,ny |
||||
DO i = 1,nx |
||||
virtual_t = tmk(i,j,k)*(0.622D0 + qvp(i,j,k))/(0.622D0*(1.D0 + qvp(i,j,k))) |
||||
rhoair = prs(i,j,k) / (RD*virtual_t) |
||||
|
||||
! Adjust factor for brightband, where snow or graupel particle |
||||
! scatters like liquid water (alpha=1.0) because it is assumed to |
||||
! have a liquid skin. |
||||
IF (iliqskin .EQ. 1 .AND. tmk(i,j,k) .GT. CELKEL) THEN |
||||
factorb_s = factor_s/ALPHA |
||||
factorb_g = factor_g/ALPHA |
||||
ELSE |
||||
factorb_s = factor_s |
||||
factorb_g = factor_g |
||||
END IF |
||||
|
||||
! Calculate variable intercept parameters |
||||
IF (ivarint .EQ. 1) THEN |
||||
|
||||
temp_c = DMIN1(-0.001D0, tmk(i,j,k)-CELKEL) |
||||
sonv = DMIN1(2.0D8, 2.0D6*EXP(-0.12D0*temp_c)) |
||||
|
||||
gonv = gon |
||||
IF (qgr(i,j,k) .GT. R1) THEN |
||||
gonv = 2.38D0 * (PI*RHO_G/(rhoair*qgr(i,j,k)))**0.92D0 |
||||
gonv = MAX(1.D4, MIN(gonv,GON)) |
||||
END IF |
||||
|
||||
ronv = RON2 |
||||
IF (qra(i,j,k) .GT. R1) THEN |
||||
ronv = RON_CONST1R*TANH((RON_QR0-qra(i,j,k))/RON_DELQR0) + RON_CONST2R |
||||
END IF |
||||
|
||||
ELSE |
||||
ronv = RN0_R |
||||
sonv = RN0_S |
||||
gonv = RN0_G |
||||
END IF |
||||
|
||||
! Total equivalent reflectivity factor (z_e, in mm^6 m^-3) is |
||||
! the sum of z_e for each hydrometeor species: |
||||
|
||||
z_e = factor_r*(rhoair*qra(i,j,k))**1.75D0/ronv**.75D0 + & |
||||
factorb_s*(rhoair*qsn(i,j,k))**1.75D0/sonv**.75D0 + & |
||||
factorb_g* (rhoair*qgr(i,j,k))**1.75D0/gonv**.75D0 |
||||
|
||||
! Adjust small values of Z_e so that dBZ is no lower than -30 |
||||
z_e = MAX(z_e, .001D0) |
||||
|
||||
! Convert to dBZ |
||||
dbz(i,j,k) = 10.D0*LOG10(z_e) |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE CALCDBZ |
||||
|
@ -0,0 +1,524 @@
@@ -0,0 +1,524 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE ROTATECOORDS(ilat, ilon, olat, olon, lat_np, lon_np, lon_0, direction) |
||||
USE constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: olat, olon |
||||
|
||||
REAL(KIND=8), INTENT(IN) :: ilat, ilon |
||||
REAL(KIND=8), INTENT(OUT) :: olat, olon |
||||
REAL(KIND=8), INTENT(IN) :: lat_np, lon_np, lon_0 |
||||
INTEGER, INTENT(IN) :: direction |
||||
|
||||
! NCLFORTEND |
||||
|
||||
! >=0, default : computational -> geographical |
||||
! < 0 : geographical -> computational |
||||
|
||||
REAL(KIND=8) :: rlat, rlon |
||||
REAL(KIND=8) :: phi_np, lam_np, lam_0, dlam |
||||
REAL(KIND=8) :: sinphi, cosphi, coslam, sinlam |
||||
!REAL(KIND=8), PARAMETER :: PI=3.141592653589793D0 |
||||
!REAL(KIND=8), PARAMETER :: RAD_PER_DEG=PI/180.D0 |
||||
!REAL(KIND=8), PARAMETER :: DEG_PER_RAD=180.D0/PI |
||||
|
||||
!convert all angles to radians |
||||
phi_np = lat_np*RAD_PER_DEG |
||||
lam_np = lon_np*RAD_PER_DEG |
||||
lam_0 = lon_0*RAD_PER_DEG |
||||
rlat = ilat*RAD_PER_DEG |
||||
rlon = ilon*RAD_PER_DEG |
||||
|
||||
IF (direction .LT. 0) THEN |
||||
! the equations are exactly the same except for one |
||||
! small difference with respect to longitude ... |
||||
dlam = pi - lam_0 |
||||
ELSE |
||||
dlam = lam_np |
||||
END IF |
||||
|
||||
sinphi = COS(phi_np)*COS(rlat)*COS(rlon - dlam) + SIN(phi_np)*SIN(rlat) |
||||
cosphi = SQRT(1.D0 - sinphi*sinphi) |
||||
coslam = SIN(phi_np)*COS(rlat)*COS(rlon - dlam) - COS(phi_np)*SIN(rlat) |
||||
sinlam = COS(rlat)*SIN(rlon - dlam) |
||||
|
||||
IF (cosphi.NE.0.D0) THEN |
||||
coslam = coslam/cosphi |
||||
sinlam = sinlam/cosphi |
||||
END IF |
||||
|
||||
olat = DEG_PER_RAD*ASIN(sinphi) |
||||
olon = DEG_PER_RAD*(ATAN2(sinlam,coslam) - dlam - lam_0 + lam_np) |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE ROTATECOORDS |
||||
|
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& |
||||
pole_lat, pole_lon, knowni, knownj, dx, dy, latinc,& |
||||
loninc, lat, lon, loc, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS |
||||
|
||||
! Converts input lat/lon values to the cartesian (i,j) value |
||||
! for the given projection. |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: loc |
||||
|
||||
|
||||
INTEGER, INTENT(IN) :: map_proj |
||||
REAL(KIND=8), INTENT(IN) :: stdlon |
||||
REAL(KIND=8), INTENT(IN) ::lat1, lon1, pole_lat, pole_lon, knowni, knownj |
||||
REAL(KIND=8), INTENT(IN) ::dx, dy, latinc, loninc |
||||
REAL(KIND=8), INTENT(INOUT) :: lat, lon, truelat1, truelat2 ! these might get modified |
||||
REAL(KIND=8), DIMENSION(2), INTENT(OUT) :: loc |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
!NCLEND |
||||
|
||||
REAL(KIND=8) :: deltalon1 |
||||
REAL(KIND=8) :: tl1r |
||||
REAL(KIND=8) :: clain, dlon, rsw, deltalon, deltalat |
||||
REAL(KIND=8) :: reflon, scale_top, ala1, alo1, ala, alo, rm, polei, polej |
||||
! earth radius divided by dx |
||||
REAL(KIND=8) :: rebydx |
||||
REAL(KIND=8) :: ctl1r, arg, cone, hemi |
||||
REAL(KIND=8) :: i, j |
||||
REAL(KIND=8) :: lat1n, lon1n, olat, olon |
||||
|
||||
! Contants |
||||
!REAL(KIND=8),PARAMETER :: PI=3.141592653589793D0 |
||||
!REAL(KIND=8),PARAMETER :: RAD_PER_DEG=PI/180.D0 |
||||
!REAL(KIND=8),PARAMETER :: DEG_PER_RAD=180.D0/PI |
||||
!REAL(KIND=8),PARAMETER :: RE_M=6370000.D0 ! radius of earth in meters |
||||
|
||||
! lat1 ! sw latitude (1,1) in degrees (-90->90n) |
||||
! lon1 ! sw longitude (1,1) in degrees (-180->180e) |
||||
! dx ! grid spacing in meters at truelats |
||||
! dlat ! lat increment for lat/lon grids |
||||
! dlon ! lon increment for lat/lon grids |
||||
! stdlon ! longitude parallel to y-axis (-180->180e) |
||||
! truelat1 ! first true latitude (all projections) |
||||
! truelat2 ! second true lat (lc only) |
||||
! hemi ! 1 for nh, -1 for sh |
||||
! cone ! cone factor for lc projections |
||||
! polei ! computed i-location of pole point |
||||
! polej ! computed j-location of pole point |
||||
! rsw ! computed radius to sw corner |
||||
! knowni ! x-location of known lat/lon |
||||
! knownj ! y-location of known lat/lon |
||||
! re_m ! radius of spherical earth, meters |
||||
! rebydx ! earth radius divided by dx |
||||
|
||||
errstat = INT(dy) ! remove compiler warning since dy not used |
||||
errstat = 0 |
||||
|
||||
rebydx = WRF_EARTH_RADIUS/dx |
||||
|
||||
! Get rid of compiler warnings |
||||
i=0 |
||||
j=0 |
||||
|
||||
hemi = 1.0D0 |
||||
IF (truelat1 .LT. 0.0D0) THEN |
||||
hemi = -1.0D0 |
||||
END IF |
||||
|
||||
! mercator |
||||
IF (map_proj.EQ.3) THEN |
||||
|
||||
! preliminary variables |
||||
clain = COS(RAD_PER_DEG*truelat1) |
||||
dlon = dx/(WRF_EARTH_RADIUS*clain) |
||||
|
||||
! compute distance from equator to origin, and store in |
||||
! the rsw tag. |
||||
rsw = 0.D0 |
||||
IF (lat1 .NE. 0.D0) THEN |
||||
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon |
||||
END IF |
||||
|
||||
deltalon = lon - lon1 |
||||
IF (deltalon .LT. -180.D0) deltalon = deltalon + 360.D0 |
||||
IF (deltalon .GT. 180.D0) deltalon = deltalon - 360.D0 |
||||
i = knowni + (deltalon/(dlon*DEG_PER_RAD)) |
||||
j = knownj + (DLOG(TAN(0.5D0*((lat + 90.D0)*RAD_PER_DEG))))/dlon - rsw |
||||
|
||||
! ps |
||||
ELSE IF (map_proj .EQ. 2) THEN |
||||
|
||||
reflon = stdlon + 90.D0 |
||||
|
||||
! compute numerator term of map scale factor |
||||
scale_top = 1.D0 + hemi*SIN(truelat1*RAD_PER_DEG) |
||||
|
||||
! compute radius to lower-left (sw) corner |
||||
ala1 = lat1*RAD_PER_DEG |
||||
rsw = rebydx*COS(ala1)*scale_top/(1.D0 + hemi*SIN(ala1)) |
||||
|
||||
! find the pole point |
||||
alo1 = (lon1 - reflon)*RAD_PER_DEG |
||||
polei = knowni - rsw*COS(alo1) |
||||
polej = knownj - hemi*rsw*SIN(alo1) |
||||
|
||||
! find radius to desired point |
||||
ala = lat*RAD_PER_DEG |
||||
rm = rebydx*COS(ala)*scale_top/(1.D0 + hemi*SIN(ala)) |
||||
alo = (lon - reflon)*RAD_PER_DEG |
||||
i = polei + rm*COS(alo) |
||||
j = polej + hemi*rm*SIN(alo) |
||||
|
||||
! lambert |
||||
ELSE IF (map_proj .EQ. 1) THEN |
||||
|
||||
IF (ABS(truelat2) .GT. 90.D0) THEN |
||||
truelat2 = truelat1 |
||||
END IF |
||||
|
||||
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN |
||||
cone = (DLOG(COS(truelat1*RAD_PER_DEG))-DLOG(COS(truelat2*RAD_PER_DEG)))/& |
||||
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0))-& |
||||
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0))) |
||||
ELSE |
||||
cone = SIN(ABS(truelat1)*RAD_PER_DEG) |
||||
END IF |
||||
|
||||
! compute longitude differences and ensure we stay |
||||
! out of the forbidden "cut zone" |
||||
deltalon1 = lon1 - stdlon |
||||
IF (deltalon1 .GT. +180.D0) deltalon1 = deltalon1 - 360.D0 |
||||
IF (deltalon1 .LT. -180.D0) deltalon1 = deltalon1 + 360.D0 |
||||
|
||||
! convert truelat1 to radian and compute cos for later use |
||||
tl1r = truelat1*RAD_PER_DEG |
||||
ctl1r = COS(tl1r) |
||||
|
||||
! compute the radius to our known lower-left (sw) corner |
||||
rsw = rebydx*ctl1r/cone*(TAN((90.D0*hemi - lat1)*RAD_PER_DEG/2.D0)/& |
||||
TAN((90.D0*hemi - truelat1)*RAD_PER_DEG/2.D0))**cone |
||||
|
||||
! find pole point |
||||
arg = cone*(deltalon1*RAD_PER_DEG) |
||||
polei = hemi*knowni - hemi*rsw*SIN(arg) |
||||
polej = hemi*knownj + rsw*COS(arg) |
||||
|
||||
! compute deltalon between known longitude and standard |
||||
! lon and ensure it is not in the cut zone |
||||
deltalon = lon - stdlon |
||||
IF (deltalon .GT. +180.D0) deltalon = deltalon - 360.D0 |
||||
IF (deltalon .LT. -180.D0) deltalon = deltalon + 360.D0 |
||||
|
||||
! radius to desired point |
||||
rm = rebydx*ctl1r/cone*(TAN((90.D0*hemi - lat)*RAD_PER_DEG/2.D0)/& |
||||
TAN((90.D0*hemi - truelat1)*RAD_PER_DEG/2.D0))**cone |
||||
|
||||
arg = cone*(deltalon*RAD_PER_DEG) |
||||
i = polei + hemi*rm*SIN(arg) |
||||
j = polej - rm*COS(arg) |
||||
|
||||
! finally, if we are in the southern hemisphere, flip the |
||||
! i/j values to a coordinate system where (1,1) is the sw |
||||
! corner (what we assume) which is different than the |
||||
! original ncep algorithms which used the ne corner as |
||||
! the origin in the southern hemisphere (left-hand vs. |
||||
! right-hand coordinate?) |
||||
i = hemi*i |
||||
j = hemi*j |
||||
|
||||
! lat-lon |
||||
ELSE IF (map_proj .EQ. 6) THEN |
||||
|
||||
IF (pole_lat .NE. 90.D0) THEN |
||||
CALL ROTATECOORDS(lat, lon, olat, olon, pole_lat, pole_lon, stdlon, -1) |
||||
lat = olat |
||||
lon = olon + stdlon |
||||
END IF |
||||
|
||||
! make sure center lat/lon is good |
||||
IF (pole_lat .NE. 90.D0) THEN |
||||
CALL ROTATECOORDS(lat1, lon1, olat, olon, pole_lat, pole_lon, stdlon, -1) |
||||
lat1n = olat |
||||
lon1n = olon + stdlon |
||||
deltalat = lat - lat1n |
||||
deltalon = lon - lon1n |
||||
ELSE |
||||
deltalat = lat - lat1 |
||||
deltalon = lon - lon1 |
||||
END IF |
||||
|
||||
! compute i/j |
||||
i = deltalon/loninc |
||||
j = deltalat/latinc |
||||
|
||||
i = i + knowni |
||||
j = j + knownj |
||||
|
||||
ELSE |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "Do not know map projection ", map_proj |
||||
RETURN |
||||
END IF |
||||
|
||||
loc(1) = j |
||||
loc(2) = i |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DLLTOIJ |
||||
|
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& |
||||
pole_lat, pole_lon, knowni, knownj, dx, dy, latinc,& |
||||
loninc, ai, aj, loc, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS |
||||
|
||||
! converts input lat/lon values to the cartesian (i,j) value |
||||
! for the given projection. |
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: loc |
||||
|
||||
INTEGER, INTENT(IN) :: map_proj |
||||
REAL(KIND=8), INTENT(IN) :: stdlon |
||||
REAL(KIND=8), INTENT(IN) :: lat1, lon1, pole_lat, pole_lon, knowni, knownj |
||||
REAL(KIND=8), INTENT(IN) :: dx, dy, latinc, loninc, ai, aj |
||||
REAL(KIND=8), INTENT(INOUT) :: truelat1, truelat2 |
||||
REAL(KIND=8), DIMENSION(2), INTENT(OUT) :: loc |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
!NCLEND |
||||
|
||||
REAL(KIND=8) :: gi2 |
||||
REAL(KIND=8) :: arccos |
||||
REAL(KIND=8) :: deltalon1 |
||||
REAL(KIND=8) :: tl1r |
||||
REAL(KIND=8) :: clain, dlon, rsw, deltalon, deltalat |
||||
REAL(KIND=8) :: reflon, scale_top, ala1, alo1, polei, polej |
||||
! earth radius divided by dx |
||||
REAL(KIND=8) :: rebydx |
||||
REAL(KIND=8) :: ctl1r, cone, hemi |
||||
|
||||
!REAL(KIND=8),PARAMETER :: PI = 3.141592653589793D0 |
||||
!REAL(KIND=8),PARAMETER :: RAD_PER_DEG = PI/180.D0 |
||||
!REAL(KIND=8),PARAMETER :: DEG_PER_RAD = 180.D0/PI |
||||
!REAL(KIND=8),PARAMETER :: RE_M = 6370000.D0 ! radius of sperical earth |
||||
|
||||
REAL(KIND=8) :: inew, jnew, r, r2 |
||||
REAL(KIND=8) :: chi, chi1, chi2 |
||||
REAL(KIND=8) :: xx, yy, lat, lon |
||||
|
||||
REAL(KIND=8) :: olat, olon, lat1n, lon1n |
||||
|
||||
|
||||
! lat1 ! sw latitude (1,1) in degrees (-90->90n) |
||||
! lon1 ! sw longitude (1,1) in degrees (-180->180e) |
||||
! dx ! grid spacing in meters at truelats |
||||
! dlat ! lat increment for lat/lon grids |
||||
! dlon ! lon increment for lat/lon grids |
||||
! stdlon ! longitude parallel to y-axis (-180->180e) |
||||
! truelat1 ! first true latitude (all projections) |
||||
! truelat2 ! second true lat (lc only) |
||||
! hemi ! 1 for nh, -1 for sh |
||||
! cone ! cone factor for lc projections |
||||
! polei ! computed i-location of pole point |
||||
! polej ! computed j-location of pole point |
||||
! rsw ! computed radius to sw corner |
||||
! knowni ! x-location of known lat/lon |
||||
! knownj ! y-location of known lat/lon |
||||
! re_m ! radius of spherical earth, meters |
||||
! rebydx ! earth radius divided by dx |
||||
|
||||
errstat = INT(dy) ! Remove compiler warning since dy not used |
||||
errstat = 0 |
||||
|
||||
rebydx = WRF_EARTH_RADIUS/dx |
||||
|
||||
hemi = 1.0D0 |
||||
IF (truelat1 .LT. 0.0D0) THEN |
||||
hemi = -1.0D0 |
||||
END IF |
||||
|
||||
! mercator |
||||
IF (map_proj .EQ. 3) THEN |
||||
! preliminary variables |
||||
clain = COS(RAD_PER_DEG*truelat1) |
||||
dlon = dx/(WRF_EARTH_RADIUS*clain) |
||||
|
||||
! compute distance from equator to origin, and store in |
||||
! the rsw tag. |
||||
rsw = 0.D0 |
||||
IF (lat1 .NE. 0.D0) THEN |
||||
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon |
||||
END IF |
||||
|
||||
lat = 2.0D0*ATAN(EXP(dlon*(rsw + aj - knownj)))*DEG_PER_RAD - 90.D0 |
||||
lon = (ai - knowni)*dlon*DEG_PER_RAD + lon1 |
||||
IF (lon .GT. 180.D0) lon = lon - 360.D0 |
||||
IF (lon .LT. -180.D0) lon = lon + 360.D0 |
||||
|
||||
! ps |
||||
ELSE IF (map_proj .EQ. 2) THEN |
||||
|
||||
! compute the reference longitude by rotating 90 degrees to |
||||
! the east to find the longitude line parallel to the |
||||
! positive x-axis. |
||||
reflon = stdlon + 90.D0 |
||||
|
||||
! compute numerator term of map scale factor |
||||
scale_top = 1.D0 + hemi*SIN(truelat1*RAD_PER_DEG) |
||||
|
||||
! compute radius to known point |
||||
ala1 = lat1*RAD_PER_DEG |
||||
rsw = rebydx*COS(ala1)*scale_top/(1.D0 + hemi*SIN(ala1)) |
||||
|
||||
! find the pole point |
||||
alo1 = (lon1 - reflon)*RAD_PER_DEG |
||||
polei = knowni - rsw*COS(alo1) |
||||
polej = knownj - hemi*rsw*SIN(alo1) |
||||
|
||||
! compute radius to point of interest |
||||
xx = ai - polei |
||||
yy = (aj - polej)*hemi |
||||
r2 = xx**2 + yy**2 |
||||
|
||||
! now the magic code |
||||
IF (r2 .EQ. 0.D0) THEN |
||||
lat = hemi*90.D0 |
||||
lon = reflon |
||||
ELSE |
||||
gi2 = (rebydx*scale_top)**2.D0 |
||||
lat = DEG_PER_RAD*hemi*ASIN((gi2 - r2)/(gi2 + r2)) |
||||
arccos = ACOS(xx/SQRT(r2)) |
||||
IF (yy .GT. 0) THEN |
||||
lon = reflon + DEG_PER_RAD*arccos |
||||
ELSE |
||||
lon = reflon - DEG_PER_RAD*arccos |
||||
END IF |
||||
END IF |
||||
|
||||
! convert to a -180 -> 180 east convention |
||||
IF (lon .GT. 180.D0) lon = lon - 360.D0 |
||||
IF (lon .LT. -180.D0) lon = lon + 360.D0 |
||||
|
||||
! !lambert |
||||
ELSE IF (map_proj .EQ. 1) THEN |
||||
|
||||
IF (ABS(truelat2) .GT. 90.D0) THEN |
||||
truelat2 = truelat1 |
||||
END IF |
||||
|
||||
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN |
||||
cone = (DLOG(COS(truelat1*RAD_PER_DEG)) - DLOG(COS(truelat2*RAD_PER_DEG)))/& |
||||
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0)) - & |
||||
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0))) |
||||
ELSE |
||||
cone = SIN(ABS(truelat1)*RAD_PER_DEG) |
||||
END IF |
||||
|
||||
! compute longitude differences and ensure we stay out of the |
||||
! forbidden "cut zone" |
||||
deltalon1 = lon1 - stdlon |
||||
IF (deltalon1 .GT. +180.D0) deltalon1 = deltalon1 - 360.D0 |
||||
IF (deltalon1 .LT. -180.D0) deltalon1 = deltalon1 + 360.D0 |
||||
|
||||
! convert truelat1 to radian and compute cos for later use |
||||
tl1r = truelat1*RAD_PER_DEG |
||||
ctl1r = COS(tl1r) |
||||
|
||||
! compute the radius to our known point |
||||
rsw = rebydx*ctl1r/cone*(TAN((90.D0*hemi - lat1)*RAD_PER_DEG/2.D0)/& |
||||
TAN((90.D0*hemi - truelat1)*RAD_PER_DEG/2.D0))**cone |
||||
|
||||
! find pole point |
||||
alo1 = cone*(deltalon1*RAD_PER_DEG) |
||||
polei = hemi*knowni - hemi*rsw*SIN(alo1) |
||||
polej = hemi*knownj + rsw*COS(alo1) |
||||
|
||||
chi1 = (90.D0 - hemi*truelat1)*RAD_PER_DEG |
||||
chi2 = (90.D0 - hemi*truelat2)*RAD_PER_DEG |
||||
|
||||
! see if we are in the southern hemispere and flip the |
||||
! indices if we are. |
||||
inew = hemi*ai |
||||
jnew = hemi*aj |
||||
|
||||
! compute radius**2 to i/j location |
||||
reflon = stdlon + 90.D0 |
||||
xx = inew - polei |
||||
yy = polej - jnew |
||||
r2 = (xx*xx + yy*yy) |
||||
r = sqrt(r2)/rebydx |
||||
|
||||
! convert to lat/lon |
||||
IF (r2 .EQ. 0.D0) THEN |
||||
lat = hemi*90.D0 |
||||
lon = stdlon |
||||
ELSE |
||||
lon = stdlon + DEG_PER_RAD*ATAN2(hemi*xx,yy)/cone |
||||
lon = dmod(lon + 360.D0, 360.D0) |
||||
IF (chi1 .EQ. chi2) THEN |
||||
chi = 2.0D0*ATAN((r/TAN(chi1))**(1.D0/cone)*TAN(chi1*0.5D0)) |
||||
ELSE |
||||
chi = 2.0D0*ATAN((r*cone/SIN(chi1))**(1.D0/cone)*TAN(chi1*0.5D0)) |
||||
END IF |
||||
lat = (90.0D0 - chi*DEG_PER_RAD)*hemi |
||||
END IF |
||||
|
||||
IF (lon .GT. +180.D0) lon = lon - 360.D0 |
||||
IF (lon .LT. -180.D0) lon = lon + 360.D0 |
||||
|
||||
! !lat-lon |
||||
ELSE IF (map_proj .EQ. 6) THEN |
||||
inew = ai - knowni |
||||
jnew = aj - knownj |
||||
|
||||
IF (inew .LT. 0.D0) inew = inew + 360.D0/loninc |
||||
IF (inew .GE. 360.D0/dx) inew = inew - 360.D0/loninc |
||||
|
||||
! compute deltalat and deltalon |
||||
deltalat = jnew*latinc |
||||
deltalon = inew*loninc |
||||
|
||||
IF (pole_lat .NE. 90.D0) THEN |
||||
CALL ROTATECOORDS(lat1, lon1, olat, olon, pole_lat, pole_lon, stdlon, -1) |
||||
lat1n = olat |
||||
lon1n = olon + stdlon |
||||
lat = deltalat + lat1n |
||||
lon = deltalon + lon1n |
||||
ELSE |
||||
lat = deltalat + lat1 |
||||
lon = deltalon + lon1 |
||||
END IF |
||||
|
||||
IF (pole_lat .NE. 90.D0) THEN |
||||
lon = lon - stdlon |
||||
CALL ROTATECOORDS(lat, lon, olat, olon, pole_lat, pole_lon, stdlon, 1) |
||||
lat = olat |
||||
lon = olon |
||||
END IF |
||||
|
||||
IF (lon .LT. -180.D0) lon = lon + 360.D0 |
||||
IF (lon .GT. 180.D0) lon = lon - 360.D0 |
||||
|
||||
ELSE |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "Do not know map projection ", map_proj |
||||
RETURN |
||||
END IF |
||||
|
||||
loc(1) = lat |
||||
loc(2) = lon |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE DIJTOLL |
@ -0,0 +1,445 @@
@@ -0,0 +1,445 @@
|
||||
!NCLFORTSTART |
||||
SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw) |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: out |
||||
|
||||
INTEGER, INTENT(IN) :: idir, ew, ns, nz, icorsw |
||||
REAL(KIND=8), INTENT(IN) :: delta |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: in |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(OUT) :: out |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz) :: lvprs |
||||
REAL(KIND=8), DIMENSION(ew,ns) :: cor |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: i, j, k, ripk, k300 |
||||
|
||||
k300 = 0 ! removes the warning |
||||
|
||||
DO j=1,ns |
||||
DO i=1,ew |
||||
IF (icorsw .EQ. 1 .AND. cor(i,j) .LT. 0.) THEN |
||||
DO k=1,nz |
||||
in(i,j,k) = -in(i,j,k) |
||||
END DO |
||||
END IF |
||||
|
||||
! First find k index that is at or below (height-wise) |
||||
! the 300 hPa level. |
||||
DO k = 1,nz |
||||
ripk = nz-k+1 |
||||
IF (lvprs(i,j,k) .LE. 300.D0) THEN |
||||
k300 = k |
||||
EXIT |
||||
END IF |
||||
END DO |
||||
|
||||
DO k = k300,1,-1 |
||||
IF (idir .EQ. 1) THEN |
||||
out(i,j,k) = MIN(in(i,j,k), in(i,j,k+1) + delta) |
||||
ELSE IF (idir .EQ. -1) THEN |
||||
out(i,j,k) = MAX(in(i,j,k), in(i,j,k+1) - delta) |
||||
END IF |
||||
END DO |
||||
|
||||
DO k = k300+1, nz |
||||
IF (idir .EQ. 1) THEN |
||||
out(i,j,k) = MAX(in(i,j,k), in(i,j,k-1) - delta) |
||||
ELSE IF (idir .EQ. -1) THEN |
||||
out(i,j,k) = MIN(in(i,j,k), in(i,j,k-1) + delta) |
||||
END IF |
||||
END DO |
||||
END DO |
||||
END DO |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE wrf_monotonic |
||||
|
||||
|
||||
!NCLFORTSTART |
||||
FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, SCLHT |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
|
||||
INTEGER, INTENT(IN) :: icase |
||||
REAL(KIND=8), INTENT(IN) :: wvalp0, wvalp1, vlev, vcp0, vcp1 |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
REAL(KIND=8) :: wrf_intrp_value |
||||
|
||||
!NCLEND |
||||
|
||||
REAL(KIND=8) :: valp0, valp1, rvalue |
||||
REAL(KIND=8) :: chkdiff |
||||
|
||||
!REAL(KIND=8), PARAMETER :: RGAS=287.04d0 |
||||
!REAL(KIND=8), PARAMETER :: USSALR=0.0065d0 |
||||
!REAL(KIND=8), PARAMETER :: SCLHT=RGAS*256.d0/9.81d0 |
||||
|
||||
errstat = 0 |
||||
|
||||
valp0 = wvalp0 |
||||
valp1 = wvalp1 |
||||
IF ( icase .EQ. 2) THEN !GHT |
||||
valp0=EXP(-wvalp0/SCLHT) |
||||
valp1=EXP(-wvalp1/SCLHT) |
||||
END IF |
||||
|
||||
chkdiff = vcp1 - vcp0 |
||||
IF(chkdiff .EQ. 0) THEN |
||||
errstat = ALGERR |
||||
errmsg = "bad difference in vcp's" |
||||
wrf_intrp_value = 0 |
||||
RETURN |
||||
!PRINT *,"bad difference in vcp's" |
||||
!STOP |
||||
END IF |
||||
|
||||
rvalue = (vlev - vcp0)*(valp1 - valp0)/(vcp1 - vcp0) + valp0 |
||||
IF (icase .EQ. 2) THEN !GHT |
||||
wrf_intrp_value = -SCLHT*LOG(rvalue) |
||||
ELSE |
||||
wrf_intrp_value = rvalue |
||||
END IF |
||||
|
||||
RETURN |
||||
|
||||
END FUNCTION wrf_intrp_value |
||||
|
||||
! NOTES: |
||||
! vcarray is the array holding the values for the vertical coordinate. |
||||
! It will always come in with the dimensions of the staggered U and V grid. |
||||
|
||||
!NCLFORTSTART |
||||
SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& |
||||
sfp, smsfp, vcarray, interp_levels, numlevels,& |
||||
icase, ew, ns, nz, extrap, vcor, logp, rmsg,& |
||||
errstat, errmsg) |
||||
USE constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & |
||||
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & |
||||
CELKEL, EPS, USSALR |
||||
|
||||
IMPLICIT NONE |
||||
|
||||
!f2py threadsafe |
||||
!f2py intent(in,out) :: dataout |
||||
|
||||
INTEGER, INTENT(IN) :: ew, ns, nz, icase, extrap |
||||
INTEGER, INTENT(IN) :: vcor, numlevels, logp |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: datain, pres, tk, qvp |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght |
||||
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: terrain, sfp, smsfp |
||||
REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout |
||||
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray |
||||
REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels |
||||
REAL(KIND=8), INTENT(IN) :: rmsg |
||||
INTEGER, INTENT(INOUT) :: errstat |
||||
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
||||
|
||||
!NCLEND |
||||
|
||||
INTEGER :: nreqlvs, ripk !njx,niy |
||||
INTEGER :: i, j, k, kupper !itriv |
||||
INTEGER :: ifound, isign !miy,mjx |
||||
REAL(KIND=8), DIMENSION(ew,ns) :: tempout |
||||
REAL(KIND=8) :: rlevel, vlev, diff |
||||
REAL(KIND=8) :: tmpvlev |
||||
REAL(KIND=8) :: vcp1, vcp0, valp0, valp1 |
||||
! REAL(KIND=8) :: cvc |
||||
REAL(KIND=8) :: vclhsl, vctophsl !qvlhsl,ttlhsl |
||||
REAL(KIND=8) :: wrf_intrp_value |
||||
REAL(KIND=8) :: plhsl, zlhsl, ezlhsl, tlhsl, psurf, pratio, tlev |
||||
REAL(KIND=8) :: ezsurf, psurfsm, zsurf, qvapor, vt |
||||
REAL(KIND=8) :: ezlev, plev, zlev, ptarget, dpmin, dp |
||||
REAL(KIND=8) :: pbot, zbot, tbotextrap, e |
||||
REAL(KIND=8) :: tlcl, gammam |
||||
CHARACTER(LEN=1) :: cvcord |
||||
|
||||
!REAL(KIND=8), PARAMETER :: RGAS = 287.04d0 !J/K/kg |
||||
!REAL(KIND=8), PARAMETER :: RGASMD = .608d0 |
||||
!REAL(KIND=8), PARAMETER :: USSALR = .0065d0 ! deg C per m |
||||
!REAL(KIND=8), PARAMETER :: SCLHT = RGAS*256.d0/9.81d0 |
||||
!REAL(KIND=8), PARAMETER :: EPS = 0.622d0 |
||||
!REAL(KIND=8), PARAMETER :: RCONST = -9.81d0/(RGAS * USSALR) |
||||
!REAL(KIND=8), PARAMETER :: EXPON = RGAS*USSALR/9.81d0 |
||||
!REAL(KIND=8), PARAMETER :: EXPONI = 1./EXPON |
||||
!REAL(KIND=8), PARAMETER :: TLCLC1 = 2840.d0 |
||||
!REAL(KIND=8), PARAMETER :: TLCLC2 = 3.5d0 |
||||
!REAL(KIND=8), PARAMETER :: TLCLC3 = 4.805d0 |
||||
!REAL(KIND=8), PARAMETER :: TLCLC4 = 55.d0 |
||||
!REAL(KIND=8), PARAMETER :: THTECON1 = 3376.d0 ! K |
||||
!REAL(KIND=8), PARAMETER :: THTECON2 = 2.54d0 |
||||
!REAL(KIND=8), PARAMETER :: THTECON3 = 0.81d0 |
||||
!REAL(KIND=8), PARAMETER :: CP = 1004.d0 |
||||
!REAL(KIND=8), PARAMETER :: CPMD = 0.887d0 |
||||
!REAL(KIND=8), PARAMETER :: GAMMA = RGAS/CP |
||||
!REAL(KIND=8), PARAMETER :: GAMMAMD = RGASMD-CPMD |
||||
!REAL(KIND=8), PARAMETER :: CELKEL = 273.16d0 |
||||
|
||||
! Removes the warnings for uninitialized variables |
||||
cvcord = '' |
||||
plev = 0 |
||||
zlev = 0 |
||||
vlev = 0 |
||||
errstat = 0 |
||||
|
||||
IF (vcor .EQ. 1) THEN |
||||
cvcord = 'p' |
||||
ELSE IF ((vcor .EQ. 2) .OR. (vcor .EQ. 3)) THEN |
||||
cvcord = 'z' |
||||
ELSE IF ((vcor .EQ. 4) .OR. (vcor .EQ. 5)) THEN |
||||
cvcord = 't' |
||||
END IF |
||||
|
||||
!miy = ns |
||||
!mjx = ew |
||||
!njx = ew |
||||
!niy = ns |
||||
|
||||
DO j = 1,ns |
||||
DO i = 1,ew |
||||
tempout(i,j) = rmsg |
||||
END DO |
||||
END DO |
||||
|
||||
DO nreqlvs = 1,numlevels |
||||
IF (cvcord .EQ. 'z') THEN |
||||
! Convert rlevel to meters from km |
||||
rlevel = interp_levels(nreqlvs) * 1000.D0 |
||||
vlev = EXP(-rlevel/SCLHT) |
||||
ELSE IF (cvcord .EQ. 'p') THEN |
||||
vlev = interp_levels(nreqlvs) |
||||
ELSE IF (cvcord .EQ. 't') THEN |
||||
vlev = interp_levels(nreqlvs) |
||||
END IF |
||||
|
||||
DO j=1,ns |
||||
DO i=1,ew |
||||
! Get the interpolated value that is within the model domain |
||||
ifound = 0 |
||||
DO k = 1,nz-1 |
||||
ripk = nz-k+1 |
||||
vcp1 = vcarray(i,j,ripk-1) |
||||
vcp0 = vcarray(i,j,ripk) |
||||
valp0 = datain(i,j,ripk) |
||||
valp1 = datain(i,j,ripk-1) |
||||
|
||||
IF ((vlev .GE. vcp0 .AND. vlev .LE. vcp1) .OR. & |
||||
(vlev .LE. vcp0 .AND. vlev .GE. vcp1)) THEN |
||||
! print *,i,j,valp0,valp1 |
||||
IF ((valp0 .EQ. rmsg) .OR. (valp1 .EQ. rmsg)) THEN |
||||
tempout(i,j) = rmsg |
||||
ifound = 1 |
||||
ELSE |
||||
IF (logp .EQ. 1) THEN |
||||
vcp1 = LOG(vcp1) |
||||
vcp0 = LOG(vcp0) |
||||
IF (vlev .EQ. 0.0D0) THEN |
||||
errstat = ALGERR |
||||
WRITE(errmsg, *) "Pres=0. Unable to take log of 0." |
||||
RETURN |
||||
!PRINT *,"Pressure value = 0" |
||||
!PRINT *,"Unable to take log of 0" |
||||
!STOP |
||||
END IF |
||||
tmpvlev = LOG(vlev) |
||||
ELSE |
||||
tmpvlev = vlev |
||||
END IF |
||||
tempout(i,j) = wrf_intrp_value(valp0, valp1, tmpvlev, vcp0, & |
||||
vcp1, icase, errstat, errmsg) |
||||
IF (errstat .NE. 0) THEN |
||||
RETURN |
||||
END IF |
||||
|
||||
! print *,"one ",i,j,tempout(i,j) |
||||
ifound = 1 |
||||
END IF |
||||
!GOTO 115 ! EXIT |
||||
EXIT |
||||
END IF |
||||
END DO !end for the k loop |
||||
!115 CONTINUE |
||||
|
||||
IF (ifound .EQ. 1) THEN !Grid point is in the model domain |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
|
||||
!If the user has requested no extrapolatin then just assign |
||||
!all values above or below the model level to rmsg. |
||||
IF (extrap .EQ. 0) THEN |
||||
tempout(i,j) = rmsg |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
|
||||
! The grid point is either above or below the model domain |
||||
! First we will check to see if the grid point is above the |
||||
! model domain. |
||||
vclhsl = vcarray(i,j,1) !lowest model level |
||||
vctophsl = vcarray(i,j,nz) !highest model level |
||||
diff = vctophsl - vclhsl |
||||
isign = NINT(diff/ABS(diff)) |
||||
|
||||
IF (isign*vlev .GE. isign*vctophsl) THEN |
||||
! Assign the highest model level to the out array |
||||
tempout(i,j) = datain(i,j,nz) |
||||
! print *,"at warn",i,j,tempout(i,j) |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
|
||||
! Only remaining possibility is that the specified level is below |
||||
! lowest model level. If lowest model level value is missing, |
||||
! set interpolated value to missing. |
||||
|
||||
IF (datain(i,j,1) .EQ. rmsg) THEN |
||||
tempout(i,j) = rmsg |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
|
||||
! If the field comming in is not a pressure,temperature or height |
||||
! field we can set the output array to the value at the lowest |
||||
! model level. |
||||
|
||||
tempout(i,j) = datain(i,j,1) |
||||
|
||||
! For the special cases of pressure on height levels or height on |
||||
! pressure levels, or temperature-related variables on pressure or |
||||
! height levels, perform a special extrapolation based on |
||||
! US Standard Atmosphere. Here we calcualate the surface pressure |
||||
! with the altimeter equation. This is how RIP calculates the |
||||
! surface pressure. |
||||
IF (icase .GT. 0) THEN |
||||
plhsl = pres(i,j,1) * 0.01D0 !pressure at lowest model level |
||||
zlhsl = ght(i,j,1) !grid point height a lowest model level |
||||
ezlhsl = EXP(-zlhsl/SCLHT) |
||||
tlhsl = tk(i,j,1) !temperature in K at lowest model level |
||||
zsurf = terrain(i,j) |
||||
qvapor = MAX((qvp(i,j,1)*.001D0),1.e-15) |
||||
! virtual temperature |
||||
! vt = tlhsl * (eps + qvapor)/(eps*(1.0 + qvapor)) |
||||
! psurf = plhsl * (vt/(vt+USSALR * (zlhsl-zsurf)))**rconst |
||||
psurf = sfp(i,j) |
||||
psurfsm = smsfp(i,j) |
||||
ezsurf = EXP(-zsurf/SCLHT) |
||||
|
||||
! The if for checking above ground |
||||
IF ((cvcord .EQ. 'z' .AND. vlev .LT. ezsurf) .OR. & |
||||
(cvcord .EQ. 'p' .AND. vlev .LT. psurf)) THEN |
||||
|
||||
! We are below the lowest data level but above the ground. |
||||
! Use linear interpolation (linear in prs and exp-height). |
||||
|
||||
IF (cvcord .EQ. 'p') THEN |
||||
plev = vlev |
||||
ezlev = ((plev - plhsl)*& |
||||
ezsurf + (psurf - plev)*ezlhsl)/(psurf - plhsl) |
||||
zlev = -SCLHT*LOG(ezlev) |
||||
IF (icase .EQ. 2) THEN |
||||
tempout(i,j) = zlev |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
|
||||
ELSE IF (cvcord .EQ. 'z') THEN |
||||
ezlev = vlev |
||||
zlev = -SCLHT*LOG(ezlev) |
||||
plev = ((ezlev - ezlhsl)*& |
||||
psurf + (ezsurf - ezlev)*plhsl)/(ezsurf - ezlhsl) |
||||
IF (icase .EQ. 1) THEN |
||||
tempout(i,j) = plev |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
END IF |
||||
|
||||
ELSE !else for checking above ground |
||||
ptarget = psurfsm - 150.D0 |
||||
dpmin = 1.e4 |
||||
DO k=1,nz |
||||
ripk = nz-k+1 |
||||
dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget) |
||||
IF (dp .GT. dpmin) THEN |
||||
!GOTO 334 ! EXIT |
||||
EXIT |
||||
END IF |
||||
dpmin = MIN(dpmin, dp) |
||||
END DO |
||||
!334 |
||||
kupper = k-1 |
||||
|
||||
ripk = nz - kupper + 1 |
||||
pbot = MAX(plhsl,psurf) |
||||
zbot = MIN(zlhsl,zsurf) |
||||
pratio = pbot/(pres(i,j,ripk) * 0.01D0) |
||||
tbotextrap = tk(i,j,ripk)*(pratio)**EXPON |
||||
! virtual temperature |
||||
vt = tbotextrap * (EPS + qvapor)/(EPS*(1.0D0 + qvapor)) |
||||
IF (cvcord .EQ. 'p') THEN |
||||
plev = vlev |
||||
zlev = zbot + vt/USSALR*(1. - (vlev/pbot)**EXPON) |
||||
IF (icase .EQ. 2) THEN |
||||
tempout(i,j) = zlev |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
ELSE IF (cvcord .EQ. 'z') THEN |
||||
zlev = -sclht*LOG(vlev) |
||||
plev = pbot*(1. + USSALR/vt*(zbot - zlev))**EXPONI |
||||
IF (icase .EQ. 1) THEN |
||||
tempout(i,j) = plev |
||||
!GOTO 333 ! CYCLE |
||||
CYCLE |
||||
END IF |
||||
END IF |
||||
END IF !end if for checking above ground |
||||
END IF !for icase gt 0 |
||||
|
||||
IF (icase .GT. 2) THEN !extrapolation for temperature |
||||
tlev = tlhsl + (zlhsl - zlev)*USSALR |
||||
qvapor = MAX(qvp(i,j,1), 1.e-15) |
||||
gammam = GAMMA*(1. + GAMMAMD*qvapor) |
||||
IF (icase .EQ. 3) THEN |
||||
tempout(i,j) = tlev - CELKEL |
||||
ELSE IF (icase .EQ. 4) THEN |
||||
tempout(i,j) = tlev |
||||
! Potential temperature - theta |
||||
ELSE IF (icase .EQ. 5) THEN |
||||
tempout(i,j) = tlev*(1000.D0/plev)**gammam |
||||
! extraolation for equivalent potential temperature |
||||
ELSE IF (icase .EQ. 6) THEN |
||||
e = qvapor*plev/(EPS + qvapor) |
||||
tlcl = TLCLC1/(LOG(tlev**TLCLC2/e) - TLCLC3) + TLCLC4 |
||||
tempout(i,j)=tlev*(1000.D0/plev)**(gammam)*& |
||||
EXP((THTECON1/tlcl - THTECON2)*& |
||||
qvapor*(1. + THTECON3*qvapor)) |
||||
END IF |
||||
END IF |
||||
|
||||
!333 CONTINUE |
||||
|
||||
END DO |
||||
END DO |
||||
|
||||
! print *,"----done----",interp_levels(nreqlvs) |
||||
DO j = 1,ns |
||||
DO i = 1,ew |
||||
dataout(i,j,nreqlvs) = tempout(i,j) |
||||
END DO |
||||
END DO |
||||
|
||||
END DO !end for the nreqlvs |
||||
|
||||
RETURN |
||||
|
||||
END SUBROUTINE wrf_vintrp |
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,60 @@
@@ -0,0 +1,60 @@
|
||||
C NCLFORTSTART |
||||
subroutine cloud_frac(pres,rh,lowc,midc,highc,nz,ns,ew) |
||||
|
||||
implicit none |
||||
integer nz,ns,ew |
||||
real pres(ew,ns,nz),rh(ew,ns,nz) |
||||
real lowc(ew,ns),midc(ew,ns),highc(ew,ns) |
||||
C NCLEND |
||||
|
||||
integer i,j,k |
||||
integer kchi,kcmi,kclo |
||||
|
||||
|
||||
DO j = 1,ns |
||||
DO i = 1,ew |
||||
DO k = 1,nz-1 |
||||
|
||||
c if((pres(i,j,k) .ge. 45000. ) .and. |
||||
c & (pres(i,j,k) .lt. 80000.)) then |
||||
c kchi = k |
||||
|
||||
c else if((pres(i,j,k) .ge. 80000.) .and. |
||||
c & (pres(i,j,k) .lt. 97000.)) then |
||||
c kcmi = k |
||||
|
||||
c else if (pres(i,j,k) .ge. 97000.) then |
||||
c kclo = k |
||||
c end if |
||||
IF ( pres(i,j,k) .gt. 97000. ) kclo=k |
||||
IF ( pres(i,j,k) .gt. 80000. ) kcmi=k |
||||
IF ( pres(i,j,k) .gt. 45000. ) kchi=k |
||||
|
||||
end do |
||||
|
||||
DO k = 1,nz-1 |
||||
IF ( k .ge. kclo .AND. k .lt. kcmi ) then |
||||
lowc(i,j) = AMAX1(rh(i,j,k),lowc(i,j)) |
||||
else IF ( k .ge. kcmi .AND. k .lt. kchi ) then !! mid cloud |
||||
midc(i,j) = AMAX1(rh(i,j,k),midc(i,j)) |
||||
else if ( k .ge. kchi ) then !! high cloud |
||||
highc(i,j) = AMAX1(rh(i,j,k),highc(i,j)) |
||||
end if |
||||
END DO |
||||
|
||||
|
||||
lowc(i,j) = 4.0 * lowc(i,j)/100.-3.0 |
||||
midc(i,j) = 4.0 * midc(i,j)/100.-3.0 |
||||
highc(i,j) = 2.5 * highc(i,j)/100.-1.5 |
||||
|
||||
lowc(i,j) = amin1(lowc(i,j),1.0) |
||||
lowc(i,j) = amax1(lowc(i,j),0.0) |
||||
midc(i,j) = amin1(midc(i,j),1.0) |
||||
midc(i,j) = amax1(midc(i,j),0.0) |
||||
highc(i,j) = amin1(highc(i,j),1.0) |
||||
highc(i,j) = amax1(highc(i,j),0.0) |
||||
|
||||
END DO |
||||
END DO |
||||
return |
||||
end |
@ -0,0 +1,170 @@
@@ -0,0 +1,170 @@
|
||||
!WRF:MODEL_LAYER:CONSTANTS |
||||
! |
||||
|
||||
MODULE module_model_constants |
||||
|
||||
! 2. Following are constants for use in defining real number bounds. |
||||
|
||||
! A really small number. |
||||
|
||||
REAL , PARAMETER :: epsilon = 1.E-15 |
||||
|
||||
! 4. Following is information related to the physical constants. |
||||
|
||||
! These are the physical constants used within the model. |
||||
|
||||
! JM NOTE -- can we name this grav instead? |
||||
REAL , PARAMETER :: g = 9.81 ! acceleration due to gravity (m {s}^-2) |
||||
|
||||
#if ( NMM_CORE == 1 ) |
||||
REAL , PARAMETER :: r_d = 287.04 |
||||
REAL , PARAMETER :: cp = 1004.6 |
||||
#else |
||||
REAL , PARAMETER :: r_d = 287. |
||||
REAL , PARAMETER :: cp = 7.*r_d/2. |
||||
#endif |
||||
|
||||
REAL , PARAMETER :: r_v = 461.6 |
||||
REAL , PARAMETER :: cv = cp-r_d |
||||
REAL , PARAMETER :: cpv = 4.*r_v |
||||
REAL , PARAMETER :: cvv = cpv-r_v |
||||
REAL , PARAMETER :: cvpm = -cv/cp |
||||
REAL , PARAMETER :: cliq = 4190. |
||||
REAL , PARAMETER :: cice = 2106. |
||||
REAL , PARAMETER :: psat = 610.78 |
||||
REAL , PARAMETER :: rcv = r_d/cv |
||||
REAL , PARAMETER :: rcp = r_d/cp |
||||
REAL , PARAMETER :: rovg = r_d/g |
||||
REAL , PARAMETER :: c2 = cp * rcv |
||||
real , parameter :: mwdry = 28.966 ! molecular weight of dry air (g/mole) |
||||
|
||||
REAL , PARAMETER :: p1000mb = 100000. |
||||
REAL , PARAMETER :: t0 = 300. |
||||
REAL , PARAMETER :: p0 = p1000mb |
||||
REAL , PARAMETER :: cpovcv = cp/(cp-r_d) |
||||
REAL , PARAMETER :: cvovcp = 1./cpovcv |
||||
REAL , PARAMETER :: rvovrd = r_v/r_d |
||||
|
||||
REAL , PARAMETER :: reradius = 1./6370.0e03 |
||||
|
||||
REAL , PARAMETER :: asselin = .025 |
||||
! REAL , PARAMETER :: asselin = .0 |
||||
REAL , PARAMETER :: cb = 25. |
||||
|
||||
REAL , PARAMETER :: XLV0 = 3.15E6 |
||||
REAL , PARAMETER :: XLV1 = 2370. |
||||
REAL , PARAMETER :: XLS0 = 2.905E6 |
||||
REAL , PARAMETER :: XLS1 = 259.532 |
||||
|
||||
REAL , PARAMETER :: XLS = 2.85E6 |
||||
REAL , PARAMETER :: XLV = 2.5E6 |
||||
REAL , PARAMETER :: XLF = 3.50E5 |
||||
|
||||
REAL , PARAMETER :: rhowater = 1000. |
||||
REAL , PARAMETER :: rhosnow = 100. |
||||
REAL , PARAMETER :: rhoair0 = 1.28 |
||||
! |
||||
! Now namelist-specified parameter: ccn_conc - RAS |
||||
! REAL , PARAMETER :: n_ccn0 = 1.0E8 |
||||
! |
||||
REAL , PARAMETER :: piconst = 3.1415926535897932384626433 |
||||
REAL , PARAMETER :: DEGRAD = piconst/180. |
||||
REAL , PARAMETER :: DPD = 360./365. |
||||
|
||||
REAL , PARAMETER :: SVP1=0.6112 |
||||
REAL , PARAMETER :: SVP2=17.67 |
||||
REAL , PARAMETER :: SVP3=29.65 |
||||
REAL , PARAMETER :: SVPT0=273.15 |
||||
REAL , PARAMETER :: EP_1=R_v/R_d-1. |
||||
REAL , PARAMETER :: EP_2=R_d/R_v |
||||
REAL , PARAMETER :: KARMAN=0.4 |
||||
REAL , PARAMETER :: EOMEG=7.2921E-5 |
||||
REAL , PARAMETER :: STBOLT=5.67051E-8 |
||||
|
||||
REAL , PARAMETER :: prandtl = 1./3.0 |
||||
! constants for w-damping option |
||||
REAL , PARAMETER :: w_alpha = 0.3 ! strength m/s/s |
||||
REAL , PARAMETER :: w_beta = 1.0 ! activation cfl number |
||||
|
||||
REAL , PARAMETER :: pq0=379.90516 |
||||
REAL , PARAMETER :: epsq2=0.2 |
||||
REAL , PARAMETER :: a2=17.2693882 |
||||
REAL , PARAMETER :: a3=273.16 |
||||
REAL , PARAMETER :: a4=35.86 |
||||
REAL , PARAMETER :: epsq=1.e-12 |
||||
REAL , PARAMETER :: p608=rvovrd-1. |
||||
!#if ( NMM_CORE == 1 ) |
||||
REAL , PARAMETER :: climit=1.e-20 |
||||
REAL , PARAMETER :: cm1=2937.4 |
||||
REAL , PARAMETER :: cm2=4.9283 |
||||
REAL , PARAMETER :: cm3=23.5518 |
||||
! REAL , PARAMETER :: defc=8.0 |
||||
! REAL , PARAMETER :: defm=32.0 |
||||
REAL , PARAMETER :: defc=0.0 |
||||
REAL , PARAMETER :: defm=99999.0 |
||||
REAL , PARAMETER :: epsfc=1./1.05 |
||||
REAL , PARAMETER :: epswet=0.0 |
||||
REAL , PARAMETER :: fcdif=1./3. |
||||
#if ( HWRF == 1 ) |
||||
REAL , PARAMETER :: fcm=0.0 |
||||
#else |
||||
REAL , PARAMETER :: fcm=0.00003 |
||||
#endif |
||||
REAL , PARAMETER :: gma=-r_d*(1.-rcp)*0.5 |
||||
REAL , PARAMETER :: p400=40000.0 |
||||
REAL , PARAMETER :: phitp=15000.0 |
||||
REAL , PARAMETER :: pi2=2.*3.1415926, pi1=3.1415926 |
||||
REAL , PARAMETER :: plbtm=105000.0 |
||||
REAL , PARAMETER :: plomd=64200.0 |
||||
REAL , PARAMETER :: pmdhi=35000.0 |
||||
REAL , PARAMETER :: q2ini=0.50 |
||||
REAL , PARAMETER :: rfcp=0.25/cp |
||||
REAL , PARAMETER :: rhcrit_land=0.75 |
||||
REAL , PARAMETER :: rhcrit_sea=0.80 |
||||
REAL , PARAMETER :: rlag=14.8125 |
||||
REAL , PARAMETER :: rlx=0.90 |
||||
REAL , PARAMETER :: scq2=50.0 |
||||
REAL , PARAMETER :: slopht=0.001 |
||||
REAL , PARAMETER :: tlc=2.*0.703972477 |
||||
REAL , PARAMETER :: wa=0.15 |
||||
REAL , PARAMETER :: wght=0.35 |
||||
REAL , PARAMETER :: wpc=0.075 |
||||
REAL , PARAMETER :: z0land=0.10 |
||||
#if ( HWRF == 1 ) |
||||
REAL , PARAMETER :: z0max=0.01 |
||||
#else |
||||
REAL , PARAMETER :: z0max=0.008 |
||||
#endif |
||||
REAL , PARAMETER :: z0sea=0.001 |
||||
!#endif |
||||
|
||||
|
||||
! Earth |
||||
|
||||
! The value for P2SI *must* be set to 1.0 for Earth |
||||
! Although, now we may not need this declaration here (see above) |
||||
!REAL , PARAMETER :: P2SI = 1.0 |
||||
|
||||
! Orbital constants: |
||||
|
||||
INTEGER , PARAMETER :: PLANET_YEAR = 365 |
||||
REAL , PARAMETER :: OBLIQUITY = 23.5 |
||||
REAL , PARAMETER :: ECCENTRICITY = 0.014 |
||||
REAL , PARAMETER :: SEMIMAJORAXIS = 1.0 ! In AU |
||||
! Don't know the following values, so we'll fake them for now |
||||
REAL , PARAMETER :: zero_date = 0.0 ! Time of perihelion passage |
||||
! Fraction into the year (from perhelion) of the |
||||
! occurrence of the Northern Spring Equinox |
||||
REAL , PARAMETER :: EQUINOX_FRACTION= 0.0 |
||||
|
||||
! 2012103 |
||||
#if (EM_CORE == 1) |
||||
! for calls to set_tiles |
||||
INTEGER, PARAMETER :: ZONE_SOLVE_EM = 1 |
||||
INTEGER, PARAMETER :: ZONE_SFS = 2 |
||||
#endif |
||||
|
||||
CONTAINS |
||||
SUBROUTINE init_module_model_constants |
||||
END SUBROUTINE init_module_model_constants |
||||
END MODULE module_model_constants |
@ -0,0 +1,380 @@
@@ -0,0 +1,380 @@
|
||||
undef("set_mp_wrf_map_resources") |
||||
function set_mp_wrf_map_resources(in_file[1]:file,opt_args[1]:logical) |
||||
|
||||
begin |
||||
; |
||||
opts = opt_args ; Make a copy of the resource list |
||||
|
||||
; Set some resources depending on what kind of map projection is |
||||
; chosen. |
||||
; |
||||
; MAP_PROJ = 0 : "CylindricalEquidistant" |
||||
; MAP_PROJ = 1 : "LambertConformal" |
||||
; MAP_PROJ = 2 : "Stereographic" |
||||
; MAP_PROJ = 3 : "Mercator" |
||||
; MAP_PROJ = 6 : "Lat/Lon" |
||||
|
||||
if(isatt(in_file,"MAP_PROJ")) |
||||
|
||||
; CylindricalEquidistant |
||||
if(in_file@MAP_PROJ .eq. 0) |
||||
projection = "CylindricalEquidistant" |
||||
opts@mpProjection = projection |
||||
opts@mpGridSpacingF = 45 |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; LambertConformal projection |
||||
if(in_file@MAP_PROJ .eq. 1) |
||||
projection = "LambertConformal" |
||||
opts@mpProjection = projection |
||||
opts@mpLambertParallel1F = get_res_value_keep(opts, "mpLambertParallel1F",in_file@TRUELAT1) |
||||
opts@mpLambertParallel2F = get_res_value_keep(opts, "mpLambertParallel2F",in_file@TRUELAT2) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
opts@mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
opts@mpLambertMeridianF = get_res_value_keep(opts, "mpLambertMeridianF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; Stereographic projection |
||||
if(in_file@MAP_PROJ .eq. 2) |
||||
projection = "Stereographic" |
||||
opts@mpProjection = projection |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@CEN_LAT) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; Mercator projection |
||||
if(in_file@MAP_PROJ .eq. 3) |
||||
projection = "Mercator" |
||||
opts@mpProjection = projection |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; global WRF CylindricalEquidistant |
||||
if(in_file@MAP_PROJ .eq. 6) |
||||
projection = "CylindricalEquidistant" |
||||
opts@mpProjection = projection |
||||
opts@mpGridSpacingF = 45 |
||||
|
||||
if (isatt(in_file,"POLE_LON") .and. isatt(in_file,"POLE_LAT") .and. isatt(in_file,"STAND_LON")) then |
||||
|
||||
if (in_file@POLE_LON .eq. 0 .and. in_file@POLE_LAT .eq. 90) then |
||||
; not rotated |
||||
|
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",180 - in_file@STAND_LON) |
||||
|
||||
else |
||||
; rotated |
||||
|
||||
southern = False ; default to northern hemisphere |
||||
if (in_file@POLE_LON .eq. 0.0) then |
||||
southern = True |
||||
else if (in_file@POLE_LON .ne. 180) then |
||||
if (isatt(in_file,"CEN_LAT") .and. in_file@CEN_LAT .lt. 0.0) then |
||||
southern = True ; probably but not necessarily true -- no way to tell for sure |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
if (.not. southern) then |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 90.0 - in_file@POLE_LAT) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", -in_file@STAND_LON) |
||||
else |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@POLE_LAT - 90) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", 180 - in_file@STAND_LON) |
||||
end if |
||||
|
||||
end if |
||||
|
||||
else if (isatt(in_file,"ref_lon") .and. isatt(in_file,"ref_lat")) then |
||||
;; this is definitely true for NMM grids but unlikely for ARW grids especially if ref_x and ref_y are set |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", in_file@REF_LAT) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", in_file@REF_LON) |
||||
|
||||
else if (isatt(in_file,"cen_lat") .and. isatt(in_file,"cen_lon")) then |
||||
;; these usually specifiy the center of the coarse domain --- not necessarily the center of the projection |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF",in_file@CEN_LAT) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF",in_file@CEN_LON) |
||||
|
||||
else |
||||
;; default values for global grid |
||||
opts@mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF", 0.0) |
||||
opts@mpCenterLonF = get_res_value_keep(opts, "mpCenterLonF", 180.0) |
||||
|
||||
end if |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
end if |
||||
|
||||
return(opts) ; Return. |
||||
|
||||
end |
||||
|
||||
|
||||
undef("wrf_map_resources") |
||||
function wrf_map_resources(in_file[1]:file,map_args[1]:logical) |
||||
local lat, lon, x1, x2, y1, y2, dims, ii, jj, southern |
||||
begin |
||||
; |
||||
; This function sets resources for a WRF map plot, basing the projection on |
||||
; the MAP_PROJ attribute in the given file. It's intended to be callable |
||||
; by users who need to set mpXXXX resources for other plotting scripts. |
||||
; |
||||
|
||||
; Set some resources depending on what kind of map projection is |
||||
; chosen. |
||||
; |
||||
; MAP_PROJ = 0 : "CylindricalEquidistant" |
||||
; MAP_PROJ = 1 : "LambertConformal" |
||||
; MAP_PROJ = 2 : "Stereographic" |
||||
; MAP_PROJ = 3 : "Mercator" |
||||
; MAP_PROJ = 6 : "Lat/Lon" |
||||
|
||||
if(isatt(in_file,"MAP_PROJ")) |
||||
|
||||
; CylindricalEquidistant |
||||
if(in_file@MAP_PROJ .eq. 0) |
||||
map_args@mpProjection = "CylindricalEquidistant" |
||||
map_args@mpGridSpacingF = 45 |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; LambertConformal projection |
||||
if(in_file@MAP_PROJ .eq. 1) |
||||
map_args@mpProjection = "LambertConformal" |
||||
map_args@mpLambertParallel1F = get_res_value_keep(map_args, "mpLambertParallel1F",in_file@TRUELAT1) |
||||
map_args@mpLambertParallel2F = get_res_value_keep(map_args, "mpLambertParallel2F",in_file@TRUELAT2) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
map_args@mpLambertMeridianF = get_res_value_keep(map_args, "mpLambertMeridianF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
map_args@mpLambertMeridianF = get_res_value_keep(map_args, "mpLambertMeridianF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; Stereographic projection |
||||
if(in_file@MAP_PROJ .eq. 2) |
||||
map_args@mpProjection = "Stereographic" |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@CEN_LAT) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; Mercator projection |
||||
if(in_file@MAP_PROJ .eq. 3) |
||||
map_args@mpProjection = "Mercator" |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0) |
||||
if(isatt(in_file,"STAND_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@STAND_LON) |
||||
else |
||||
if(isatt(in_file,"CEN_LON")) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON) |
||||
else |
||||
print("ERROR: Found neither STAND_LON or CEN_LON in file") |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
; global WRF CylindricalEquidistant |
||||
if(in_file@MAP_PROJ .eq. 6) |
||||
print ("YES, THIS WORKED") |
||||
projection = "CylindricalEquidistant" |
||||
map_args@mpProjection = projection |
||||
map_args@mpGridSpacingF = 45 |
||||
|
||||
;; according to the docs if POLE_LON is 0 then the projection center is in the southern hemisphere |
||||
;; if POLE_LON is 180 the projection center is in the northern hemisphere |
||||
;; otherwise you can't tell for sure -- CEN_LAT does not have to be the projection center but hopefully |
||||
;; it is in the same hemisphere. The same is true for REF_LAT except that if REF_Y is specified REF_LAT might |
||||
;; be in a corner or somewhere else and therefore it is even less reliable |
||||
;; |
||||
|
||||
if (isatt(in_file,"POLE_LON") .and. isatt(in_file,"POLE_LAT") .and. isatt(in_file,"STAND_LON")) then |
||||
|
||||
if (in_file@POLE_LON .eq. 0 .and. in_file@POLE_LAT .eq. 90) then |
||||
; not rotated |
||||
|
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",180 - in_file@STAND_LON) |
||||
|
||||
else |
||||
; rotated |
||||
|
||||
southern = False ; default to northern hemisphere |
||||
if (in_file@POLE_LON .eq. 0.0) then |
||||
southern = True |
||||
else if (in_file@POLE_LON .ne. 180) then |
||||
if (isatt(in_file,"CEN_LAT") .and. in_file@CEN_LAT .lt. 0.0) then |
||||
southern = True ; probably but not necessarily true -- no way to tell for sure |
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
if (.not. southern) then |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 90.0 - in_file@POLE_LAT) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", -in_file@STAND_LON) |
||||
else |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@POLE_LAT - 90) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", 180 - in_file@STAND_LON) |
||||
end if |
||||
|
||||
end if |
||||
|
||||
else if (isatt(in_file,"ref_lon") .and. isatt(in_file,"ref_lat")) then |
||||
;; this is definitely true for NMM grids but unlikely for ARW grids especially if ref_x and ref_y are set |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", in_file@REF_LAT) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", in_file@REF_LON) |
||||
|
||||
else if (isatt(in_file,"cen_lat") .and. isatt(in_file,"cen_lon")) then |
||||
;; these usually specifiy the center of the coarse domain --- not necessarily the center of the projection |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF",in_file@CEN_LAT) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF",in_file@CEN_LON) |
||||
|
||||
else |
||||
;; default values for global grid |
||||
map_args@mpCenterLatF = get_res_value_keep(map_args, "mpCenterLatF", 0.0) |
||||
map_args@mpCenterLonF = get_res_value_keep(map_args, "mpCenterLonF", 180.0) |
||||
|
||||
end if |
||||
end if |
||||
end if |
||||
|
||||
end if |
||||
|
||||
else |
||||
|
||||
return(map_args) |
||||
|
||||
end if |
||||
|
||||
map_args@mpNestTime = get_res_value_keep(map_args, "mpNestTime",0) |
||||
|
||||
if(isfilevar(in_file,"XLAT")) |
||||
lat = in_file->XLAT(map_args@mpNestTime,:,:) |
||||
lon = in_file->XLONG(map_args@mpNestTime,:,:) |
||||
else |
||||
lat = in_file->XLAT_M(map_args@mpNestTime,:,:) |
||||
lon = in_file->XLONG_M(map_args@mpNestTime,:,:) |
||||
end if |
||||
dims = dimsizes(lat) |
||||
|
||||
do ii = 0, dims(0)-1 |
||||
do jj = 0, dims(1)-1 |
||||
if ( lon(ii,jj) .lt. 0.0) then |
||||
lon(ii,jj) = lon(ii,jj) + 360. |
||||
end if |
||||
end do |
||||
end do |
||||
|
||||
map_args@start_lat = lat(0,0) |
||||
map_args@start_lon = lon(0,0) |
||||
map_args@end_lat = lat(dims(0)-1,dims(1)-1) |
||||
map_args@end_lon = lon(dims(0)-1,dims(1)-1) |
||||
|
||||
; end_lon must be greater than start_lon, or errors are thrown |
||||
if (map_args@end_lon .le. map_args@start_lon) then |
||||
map_args@end_lon = map_args@end_lon + 360.0 |
||||
end if |
||||
|
||||
|
||||
; Set some resources common to all map projections. |
||||
map_args = set_mp_resources(map_args) |
||||
|
||||
if ( isatt(map_args,"ZoomIn") .and. map_args@ZoomIn ) then |
||||
y1 = 0 |
||||
x1 = 0 |
||||
y2 = dims(0)-1 |
||||
x2 = dims(1)-1 |
||||
if ( isatt(map_args,"Ystart") ) then |
||||
y1 = map_args@Ystart |
||||
delete(map_args@Ystart) |
||||
end if |
||||
if ( isatt(map_args,"Xstart") ) then |
||||
x1 = map_args@Xstart |
||||
delete(map_args@Xstart) |
||||
end if |
||||
if ( isatt(map_args,"Yend") ) then |
||||
if ( map_args@Yend .le. y2 ) then |
||||
y2 = map_args@Yend |
||||
end if |
||||
delete(map_args@Yend) |
||||
end if |
||||
if ( isatt(map_args,"Xend") ) then |
||||
if ( map_args@Xend .le. x2 ) then |
||||
x2 = map_args@Xend |
||||
end if |
||||
delete(map_args@Xend) |
||||
end if |
||||
|
||||
map_args@mpLeftCornerLatF = lat(y1,x1) |
||||
map_args@mpLeftCornerLonF = lon(y1,x1) |
||||
map_args@mpRightCornerLatF = lat(y2,x2) |
||||
map_args@mpRightCornerLonF = lon(y2,x2) |
||||
|
||||
if ( map_args@mpRightCornerLonF .lt. 0.0 ) then |
||||
map_args@mpRightCornerLonF = map_args@mpRightCornerLonF + 360.0 |
||||
end if |
||||
|
||||
if ( map_args@mpRightCornerLonF .le. map_args@mpRightCornerLonF ) then |
||||
map_args@mpRightCornerLonF = map_args@mpRightCornerLonF + 360.0 |
||||
end if |
||||
|
||||
delete(map_args@ZoomIn) |
||||
end if |
||||
|
||||
return(map_args) |
||||
end |
@ -0,0 +1,27 @@
@@ -0,0 +1,27 @@
|
||||
from __future__ import (absolute_import, division, print_function, |
||||
unicode_literals) |
||||
|
||||
from .constants import Constants |
||||
from .extension import _tk, _rh, _cloudfrac |
||||
from .metadecorators import set_cloudfrac_metadata |
||||
from .util import extract_vars |
||||
|
||||
@set_cloudfrac_metadata() |
||||
def get_cloudfrac(wrfnc, timeidx=0, method="cat", squeeze=True, |
||||
cache=None, meta=True): |
||||
|
||||
vars = extract_vars(wrfnc, timeidx, ("P", "PB", "QVAPOR", "T"), |
||||
method, squeeze, cache, meta=False) |
||||
|
||||
p = vars["P"] |
||||
pb = vars["PB"] |
||||
qv = vars["QVAPOR"] |
||||
t = vars["T"] |
||||
|
||||
full_p = p + pb |
||||
full_t = t + Constants.T_BASE |
||||
|
||||
tk = _tk(full_p, full_t) |
||||
rh = _rh(qv, full_p, tk) |
||||
|
||||
return _cloudfrac(full_p, rh) |
@ -0,0 +1,55 @@
@@ -0,0 +1,55 @@
|
||||
from sys import version_info |
||||
|
||||
# Dictionary python 2-3 compatibility stuff |
||||
def viewitems(d): |
||||
func = getattr(d, "viewitems", None) |
||||
if func is None: |
||||
func = d.items |
||||
return func() |
||||
|
||||
|
||||
def viewkeys(d): |
||||
func = getattr(d, "viewkeys", None) |
||||
if func is None: |
||||
func = d.keys |
||||
return func() |
||||
|
||||
|
||||
def viewvalues(d): |
||||
func = getattr(d, "viewvalues", None) |
||||
if func is None: |
||||
func = d.values |
||||
return func() |
||||
|
||||
def isstr(s): |
||||
try: |
||||
return isinstance(s, basestring) |
||||
except NameError: |
||||
return isinstance(s, str) |
||||
|
||||
|
||||
# Python 2 rounding behavior |
||||
def _round2(x, d=0): |
||||
p = 10 ** d |
||||
return float(floor((x * p) + copysign(0.5, x)))/p |
||||
|
||||
|
||||
def py2round(x, d=0): |
||||
if version_info >= (3,): |
||||
return _round2(x, d) |
||||
|
||||
return round(x, d) |
||||
|
||||
|
||||
def py3range(*args): |
||||
if version_info >= (3,): |
||||
return range(*args) |
||||
|
||||
return xrange(*args) |
||||
|
||||
|
||||
def ucode(*args, **kwargs): |
||||
if version_info >= (3, ): |
||||
return str(*args, **kwargs) |
||||
|
||||
return unicode(*args, **kwargs) |
@ -0,0 +1,489 @@
@@ -0,0 +1,489 @@
|
||||
from __future__ import (absolute_import, division, print_function, |
||||
unicode_literals) |
||||
|
||||
import numpy as np |
||||
|
||||
import wrapt |
||||
|
||||
#from .destag import destagger |
||||
from .util import iter_left_indexes, npvalues |
||||
from .py3compat import py3range |
||||
from .config import xarray_enabled |
||||
from .constants import Constants |
||||
|
||||
if xarray_enabled(): |
||||
from xarray import DataArray |
||||
|
||||
|
||||
# def uvmet_left_iter(): |
||||
# """Decorator to handle iterating over leftmost dimensions when using |
||||
# multiple files and/or multiple times with the uvmet product. |
||||
# |
||||
# """ |
||||
# @wrapt.decorator |
||||
# def func_wrapper(wrapped, instance, args, kwargs): |
||||
# u = args[0] |
||||
# v = args[1] |
||||
# lat = args[2] |
||||
# lon = args[3] |
||||
# cen_long = args[4] |
||||
# cone = args[5] |
||||
# |
||||
# if u.ndim == lat.ndim: |
||||
# num_right_dims = 2 |
||||
# is_3d = False |
||||
# else: |
||||
# num_right_dims = 3 |
||||
# is_3d = True |
||||
# |
||||
# is_stag = False |
||||
# if ((u.shape[-1] != lat.shape[-1]) or |
||||
# (u.shape[-2] != lat.shape[-2])): |
||||
# is_stag = True |
||||
# |
||||
# if is_3d: |
||||
# extra_dim_num = u.ndim - 3 |
||||
# else: |
||||
# extra_dim_num = u.ndim - 2 |
||||
# |
||||
# if is_stag: |
||||
# u = destagger(u,-1) |
||||
# v = destagger(v,-2) |
||||
# |
||||
# # No special left side iteration, return the function result |
||||
# if (extra_dim_num == 0): |
||||
# return wrapped(u, v, lat, lon, cen_long, cone) |
||||
# |
||||
# # Start by getting the left-most 'extra' dims |
||||
# outdims = u.shape[0:extra_dim_num] |
||||
# extra_dims = list(outdims) # Copy the left-most dims for iteration |
||||
# |
||||
# # Append the right-most dimensions |
||||
# outdims += [2] # For u/v components |
||||
# |
||||
# #outdims += [u.shape[x] for x in py3range(-num_right_dims,0,1)] |
||||
# outdims += list(u.shape[-num_right_dims:]) |
||||
# |
||||
# output = np.empty(outdims, u.dtype) |
||||
# |
||||
# for left_idxs in iter_left_indexes(extra_dims): |
||||
# # Make the left indexes plus a single slice object |
||||
# # The single slice will handle all the dimensions to |
||||
# # the right (e.g. [1,1,:]) |
||||
# left_and_slice_idxs = tuple([x for x in left_idxs] + [slice(None)]) |
||||
# |
||||
# new_u = u[left_and_slice_idxs] |
||||
# new_v = v[left_and_slice_idxs] |
||||
# new_lat = lat[left_and_slice_idxs] |
||||
# new_lon = lon[left_and_slice_idxs] |
||||
# |
||||
# # Call the numerical routine |
||||
# result = wrapped(new_u, new_v, new_lat, new_lon, cen_long, cone) |
||||
# |
||||
# # Note: The 2D version will return a 3D array with a 1 length |
||||
# # dimension. Numpy is unable to broadcast this without |
||||
# # sqeezing first. |
||||
# result = np.squeeze(result) |
||||
# |
||||
# output[left_and_slice_idxs] = result[:] |
||||
# |
||||
# return output |
||||
# |
||||
# return func_wrapper |
||||
|
||||
def _move_dim_to_left(arr, dimidx): |
||||
|
||||
if isinstance(arr, np.ma.MaskedArray): |
||||
has_missing = True |
||||
missing = result.fill_value |
||||
|
||||
shape = list(arr.shape) |
||||
move_dim_size = shape.pop(dimidx) |
||||
|
||||
output_dims = [move_dim_size] + shape |
||||
|
||||
output = np.empty(output_dims, arr.dtype) |
||||
|
||||
if dimidx < 0: |
||||
right_ndim = -dimidx |
||||
else: |
||||
right_ndim = arr.ndim - dimidx |
||||
|
||||
rightdims = [slice(None)] * right_ndim |
||||
|
||||
for i in py3range(move_dim_size): |
||||
rightdims[0] = i |
||||
outidxs = (Ellipsis,) + tuple(rightdims) |
||||
output[i,:] = outview_array[outidxs] |
||||
|
||||
if has_missing: |
||||
output = np.ma.masked_values(output, missing) |
||||
|
||||
return output |
||||
|
||||
|
||||
def uvmet_left_iter_nocopy(alg_dtype=np.float64): |
||||
"""Decorator to handle iterating over leftmost dimensions when using |
||||
multiple files and/or multiple times with the uvmet product. |
||||
|
||||
""" |
||||
@wrapt.decorator |
||||
def func_wrapper(wrapped, instance, args, kwargs): |
||||
u = args[0] |
||||
v = args[1] |
||||
lat = args[2] |
||||
lon = args[3] |
||||
cen_long = args[4] |
||||
cone = args[5] |
||||
|
||||
orig_dtype = u.dtype |
||||
|
||||
lat_lon_fixed = False |
||||
if lat.ndim == 2: |
||||
lat_lon_fixed = True |
||||
|
||||
if lon.ndim == 2 and not lat_lon_fixed: |
||||
raise ValueError("'lat' and 'lon' shape mismatch") |
||||
|
||||
num_left_dims_u = u.ndim - 2 |
||||
num_left_dims_lat = lat.ndim - 2 |
||||
|
||||
if (num_left_dims_lat > num_left_dims_u): |
||||
raise ValueError("number of 'lat' dimensions is greater than 'u'") |
||||
|
||||
if lat_lon_fixed: |
||||
mode = 0 # fixed lat/lon |
||||
else: |
||||
if num_left_dims_u == num_left_dims_lat: |
||||
mode = 1 # lat/lon same as u |
||||
else: |
||||
mode = 2 # probably 3D with 2D lat/lon plus Time |
||||
|
||||
has_missing = False |
||||
u_arr = u |
||||
if isinstance(u, DataArray): |
||||
u_arr = npvalues(u) |
||||
|
||||
v_arr = v |
||||
if isinstance(v, DataArray): |
||||
v_arr = npvalues(v) |
||||
|
||||
umissing = Constants.DEFAULT_FILL |
||||
if isinstance(u_arr, np.ma.MaskedArray): |
||||
has_missing = True |
||||
umissing = u_arr.fill_value |
||||
|
||||
vmissing = Constants.DEFAULT_FILL |
||||
if isinstance(v_arr, np.ma.MaskedArray): |
||||
has_missing = True |
||||
vmissing = v_arr.fill_value |
||||
|
||||
uvmetmissing = umissing |
||||
|
||||
is_stag = 0 |
||||
if (u.shape[-1] != lat.shape[-1] or u.shape[-2] != lat.shape[-2]): |
||||
is_stag = 1 |
||||
# Sanity check |
||||
if (v.shape[-1] == lat.shape[-1] or v.shape[-2] == lat.shape[-2]): |
||||
raise ValueError("u is staggered but v is not") |
||||
|
||||
if (v.shape[-1] != lat.shape[-1] or v.shape[-2] != lat.shape[-2]): |
||||
is_stag = 1 |
||||
# Sanity check |
||||
if (u.shape[-1] == lat.shape[-1] or u.shape[-2] == lat.shape[-2]): |
||||
raise ValueError("v is staggered but u is not") |
||||
|
||||
|
||||
|
||||
# No special left side iteration, return the function result |
||||
if (num_left_dims_u == 0): |
||||
return wrapped(u, v, lat, lon, cen_long, cone, isstag=is_stag, |
||||
has_missing=has_missing, umissing=umissing, |
||||
vmissing=vmissing, uvmetmissing=uvmetmissing) |
||||
|
||||
# Initial output is time,nz,2,ny,nx to create contiguous views |
||||
outdims = u.shape[0:num_left_dims_u] |
||||
extra_dims = tuple(outdims) # Copy the left-most dims for iteration |
||||
|
||||
outdims += (2,) |
||||
|
||||
outdims += lat.shape[-2:] |
||||
|
||||
outview_array = np.empty(outdims, alg_dtype) |
||||
|
||||
# Final Output moves the u_v dimension to left side |
||||
output_dims = (2,) |
||||
output_dims += extra_dims |
||||
output_dims += lat.shape[-2:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
for left_idxs in iter_left_indexes(extra_dims): |
||||
left_and_slice_idxs = left_idxs + (slice(None),) |
||||
|
||||
if mode == 0: |
||||
lat_left_and_slice = (slice(None),) |
||||
elif mode == 1: |
||||
lat_left_and_slice = left_and_slice_idxs |
||||
elif mode == 2: |
||||
# Only need the left-most |
||||
lat_left_and_slice = tuple(left_idx |
||||
for left_idx in left_idxs[0:num_left_dims_lat]) |
||||
|
||||
u_output_idxs = (0,) + left_idxs + (slice(None),) |
||||
v_output_idxs = (1,) + left_idxs + (slice(None),) |
||||
u_view_idxs = left_idxs + (0, slice(None)) |
||||
v_view_idxs = left_idxs + (1, slice(None)) |
||||
|
||||
|
||||
new_u = u[left_and_slice_idxs] |
||||
new_v = v[left_and_slice_idxs] |
||||
new_lat = lat[lat_left_and_slice] |
||||
new_lon = lon[lat_left_and_slice] |
||||
outview = outview_array[left_and_slice_idxs] |
||||
|
||||
# Call the numerical routine |
||||
result = wrapped(new_u, new_v, new_lat, new_lon, cen_long, cone, |
||||
isstag=is_stag, has_missing=has_missing, |
||||
umissing=umissing, vmissing=vmissing, |
||||
uvmetmissing=uvmetmissing, outview=outview) |
||||
|
||||
# Make sure the result is the same data as what got passed in |
||||
# Can delete this once everything works |
||||
if (result.__array_interface__["data"][0] != |
||||
outview.__array_interface__["data"][0]): |
||||
raise RuntimeError("output array was copied") |
||||
|
||||
output[u_output_idxs] = ( |
||||
outview_array[u_view_idxs].astype(orig_dtype)) |
||||
output[v_output_idxs] = ( |
||||
outview_array[v_view_idxs].astype(orig_dtype)) |
||||
|
||||
if has_missing: |
||||
output = np.ma.masked_values(output, uvmetmissing) |
||||
|
||||
return output |
||||
|
||||
return func_wrapper |
||||
|
||||
|
||||
|
||||
def cape_left_iter(alg_dtype=np.float64): |
||||
"""Decorator to handle iterating over leftmost dimensions when using |
||||
multiple files and/or multiple times with the cape product. |
||||
|
||||
|
||||
""" |
||||
@wrapt.decorator |
||||
def func_wrapper(wrapped, instance, args, kwargs): |
||||
# The cape calculations use an ascending vertical pressure coordinate |
||||
|
||||
new_args = list(args) |
||||
new_kwargs = dict(kwargs) |
||||
|
||||
p_hpa = args[0] |
||||
tk = args[1] |
||||
qv = args[2] |
||||
ht = args[3] |
||||
ter = args[4] |
||||
sfp = args[5] |
||||
missing = args[6] |
||||
i3dflag = args[7] |
||||
ter_follow = args[8] |
||||
|
||||
is2d = False if i3dflag != 0 else True |
||||
|
||||
# Need to order in ascending pressure order |
||||
flip = False |
||||
bot_idxs = (0,) * p_hpa.ndim |
||||
top_idxs = list(bot_idxs) |
||||
top_idxs[-3] = -1 |
||||
top_idxs = tuple(top_idxs) |
||||
|
||||
if p_hpa[bot_idxs] > p_hpa[top_idxs]: |
||||
flip = True |
||||
p_hpa = np.ascontiguousarray(p_hpa[...,::-1,:,:]) |
||||
tk = np.ascontiguousarray(tk[...,::-1,:,:]) |
||||
qv = np.ascontiguousarray(qv[...,::-1,:,:]) |
||||
ht = np.ascontiguousarray(ht[...,::-1,:,:]) |
||||
new_args[0] = p_hpa |
||||
new_args[1] = tk |
||||
new_args[2] = qv |
||||
new_args[3] = ht |
||||
|
||||
num_left_dims = p_hpa.ndim - 3 |
||||
orig_dtype = p_hpa.dtype |
||||
|
||||
# No special left side iteration, build the output from the cape,cin |
||||
# result |
||||
if (num_left_dims == 0): |
||||
cape, cin = wrapped(*new_args, **new_kwargs) |
||||
|
||||
output_dims = (2,) |
||||
output_dims += p_hpa.shape[-3:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
if flip and not is2d: |
||||
output[0,:] = cape[::-1,:,:] |
||||
output[1,:] = cin[::-1,:,:] |
||||
else: |
||||
output[0,:] = cape[:] |
||||
output[1,:] = cin[:] |
||||
|
||||
return output |
||||
|
||||
|
||||
# Initial output is ...,cape_cin,nz,ny,nx to create contiguous views |
||||
outdims = p_hpa.shape[0:num_left_dims] |
||||
extra_dims = tuple(outdims) # Copy the left-most dims for iteration |
||||
|
||||
outdims += (2,) # cape_cin |
||||
|
||||
outdims += p_hpa.shape[-3:] |
||||
|
||||
outview_array = np.empty(outdims, alg_dtype) |
||||
|
||||
# Create the output array where the leftmost dim is the product type |
||||
output_dims = (2,) |
||||
output_dims += extra_dims |
||||
output_dims += p_hpa.shape[-3:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
for left_idxs in iter_left_indexes(extra_dims): |
||||
left_and_slice_idxs = left_idxs + (slice(None),) |
||||
cape_idxs = left_idxs + (0, slice(None)) |
||||
cin_idxs = left_idxs + (1, slice(None)) |
||||
|
||||
cape_output_idxs = (0,) + left_idxs + (slice(None),) |
||||
cin_output_idxs = (1,) + left_idxs + (slice(None),) |
||||
view_cape_reverse_idxs = left_idxs + (0, slice(None,None,-1), |
||||
slice(None)) |
||||
view_cin_reverse_idxs = left_idxs + (1, slice(None,None,-1), |
||||
slice(None)) |
||||
|
||||
new_args[0] = p_hpa[left_and_slice_idxs] |
||||
new_args[1] = tk[left_and_slice_idxs] |
||||
new_args[2] = qv[left_and_slice_idxs] |
||||
new_args[3] = ht[left_and_slice_idxs] |
||||
new_args[4] = ter[left_and_slice_idxs] |
||||
new_args[5] = sfp[left_and_slice_idxs] |
||||
capeview = outview_array[cape_idxs] |
||||
cinview = outview_array[cin_idxs] |
||||
|
||||
# Call the numerical routine |
||||
new_kwargs["capeview"] = capeview |
||||
new_kwargs["cinview"] = cinview |
||||
|
||||
cape, cin = wrapped(*new_args, **new_kwargs) |
||||
|
||||
# Make sure the result is the same data as what got passed in |
||||
# Can delete this once everything works |
||||
if (cape.__array_interface__["data"][0] != |
||||
capeview.__array_interface__["data"][0]): |
||||
raise RuntimeError("output array was copied") |
||||
|
||||
|
||||
if flip and not is2d: |
||||
output[cape_output_idxs] = ( |
||||
outview_array[view_cape_reverse_idxs].astype(orig_dtype)) |
||||
output[cin_output_idxs] = ( |
||||
outview_array[view_cin_reverse_idxs].astype(orig_dtype)) |
||||
else: |
||||
output[cape_output_idxs] = ( |
||||
outview_array[cape_idxs].astype(orig_dtype)) |
||||
output[cin_output_idxs] = ( |
||||
outview_array[cin_idxs].astype(orig_dtype)) |
||||
|
||||
return output |
||||
|
||||
return func_wrapper |
||||
|
||||
def cloudfrac_left_iter(alg_dtype=np.float64): |
||||
"""Decorator to handle iterating over leftmost dimensions when using |
||||
multiple files and/or multiple times with the cloudfrac product. |
||||
|
||||
|
||||
""" |
||||
@wrapt.decorator |
||||
def func_wrapper(wrapped, instance, args, kwargs): |
||||
# The cape calculations use an ascending vertical pressure coordinate |
||||
new_args = list(args) |
||||
new_kwargs = dict(kwargs) |
||||
|
||||
p = args[0] |
||||
rh = args[1] |
||||
|
||||
num_left_dims = p.ndim - 3 |
||||
orig_dtype = p.dtype |
||||
|
||||
# No special left side iteration, build the output from the cape,cin |
||||
# result |
||||
if (num_left_dims == 0): |
||||
low, med, high = wrapped(*new_args, **new_kwargs) |
||||
|
||||
output_dims = (3,) |
||||
output_dims += p.shape[-2:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
output[0,:] = low[:] |
||||
output[1,:] = med[:] |
||||
output[2,:] = high[:] |
||||
|
||||
return output |
||||
|
||||
|
||||
# Initial output is ...,cape_cin,nz,ny,nx to create contiguous views |
||||
outdims = p.shape[0:num_left_dims] |
||||
extra_dims = tuple(outdims) # Copy the left-most dims for iteration |
||||
|
||||
outdims += (3,) # low_mid_high |
||||
|
||||
outdims += p.shape[-2:] |
||||
|
||||
outview_array = np.empty(outdims, alg_dtype) |
||||
|
||||
# Create the output array where the leftmost dim is the cloud type |
||||
output_dims = (3,) |
||||
output_dims += extra_dims |
||||
output_dims += p.shape[-2:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
for left_idxs in iter_left_indexes(extra_dims): |
||||
left_and_slice_idxs = left_idxs + (slice(None),) |
||||
low_idxs = left_idxs + (0, slice(None)) |
||||
med_idxs = left_idxs + (1, slice(None)) |
||||
high_idxs = left_idxs + (2, slice(None)) |
||||
|
||||
low_output_idxs = (0,) + left_idxs + (slice(None),) |
||||
med_output_idxs = (1,) + left_idxs + (slice(None),) |
||||
high_output_idxs = (2,) + left_idxs + (slice(None),) |
||||
|
||||
new_args[0] = p[left_and_slice_idxs] |
||||
new_args[1] = rh[left_and_slice_idxs] |
||||
|
||||
lowview = outview_array[low_idxs] |
||||
medview = outview_array[med_idxs] |
||||
highview = outview_array[high_idxs] |
||||
|
||||
new_kwargs["lowview"] = lowview |
||||
new_kwargs["medview"] = medview |
||||
new_kwargs["highview"] = highview |
||||
|
||||
low, med, high = wrapped(*new_args, **new_kwargs) |
||||
|
||||
# Make sure the result is the same data as what got passed in |
||||
# Can delete this once everything works |
||||
if (low.__array_interface__["data"][0] != |
||||
lowview.__array_interface__["data"][0]): |
||||
raise RuntimeError("output array was copied") |
||||
|
||||
output[low_output_idxs] = ( |
||||
outview_array[low_idxs].astype(orig_dtype)) |
||||
output[med_output_idxs] = ( |
||||
outview_array[med_idxs].astype(orig_dtype)) |
||||
output[high_output_idxs] = ( |
||||
outview_array[high_idxs].astype(orig_dtype)) |
||||
|
||||
return output |
||||
|
||||
return func_wrapper |
||||
|
@ -1,229 +0,0 @@
@@ -1,229 +0,0 @@
|
||||
from __future__ import (absolute_import, division, print_function, |
||||
unicode_literals) |
||||
|
||||
import numpy as np |
||||
|
||||
import wrapt |
||||
|
||||
#from .destag import destagger |
||||
from .util import iter_left_indexes, py3range, npvalues |
||||
from .config import xarray_enabled |
||||
from .constants import Constants |
||||
|
||||
if xarray_enabled(): |
||||
from xarray import DataArray |
||||
|
||||
|
||||
# def uvmet_left_iter(): |
||||
# """Decorator to handle iterating over leftmost dimensions when using |
||||
# multiple files and/or multiple times with the uvmet product. |
||||
# |
||||
# """ |
||||
# @wrapt.decorator |
||||
# def func_wrapper(wrapped, instance, args, kwargs): |
||||
# u = args[0] |
||||
# v = args[1] |
||||
# lat = args[2] |
||||
# lon = args[3] |
||||
# cen_long = args[4] |
||||
# cone = args[5] |
||||
# |
||||
# if u.ndim == lat.ndim: |
||||
# num_right_dims = 2 |
||||
# is_3d = False |
||||
# else: |
||||
# num_right_dims = 3 |
||||
# is_3d = True |
||||
# |
||||
# is_stag = False |
||||
# if ((u.shape[-1] != lat.shape[-1]) or |
||||
# (u.shape[-2] != lat.shape[-2])): |
||||
# is_stag = True |
||||
# |
||||
# if is_3d: |
||||
# extra_dim_num = u.ndim - 3 |
||||
# else: |
||||
# extra_dim_num = u.ndim - 2 |
||||
# |
||||
# if is_stag: |
||||
# u = destagger(u,-1) |
||||
# v = destagger(v,-2) |
||||
# |
||||
# # No special left side iteration, return the function result |
||||
# if (extra_dim_num == 0): |
||||
# return wrapped(u, v, lat, lon, cen_long, cone) |
||||
# |
||||
# # Start by getting the left-most 'extra' dims |
||||
# outdims = u.shape[0:extra_dim_num] |
||||
# extra_dims = list(outdims) # Copy the left-most dims for iteration |
||||
# |
||||
# # Append the right-most dimensions |
||||
# outdims += [2] # For u/v components |
||||
# |
||||
# #outdims += [u.shape[x] for x in py3range(-num_right_dims,0,1)] |
||||
# outdims += list(u.shape[-num_right_dims:]) |
||||
# |
||||
# output = np.empty(outdims, u.dtype) |
||||
# |
||||
# for left_idxs in iter_left_indexes(extra_dims): |
||||
# # Make the left indexes plus a single slice object |
||||
# # The single slice will handle all the dimensions to |
||||
# # the right (e.g. [1,1,:]) |
||||
# left_and_slice_idxs = tuple([x for x in left_idxs] + [slice(None)]) |
||||
# |
||||
# new_u = u[left_and_slice_idxs] |
||||
# new_v = v[left_and_slice_idxs] |
||||
# new_lat = lat[left_and_slice_idxs] |
||||
# new_lon = lon[left_and_slice_idxs] |
||||
# |
||||
# # Call the numerical routine |
||||
# result = wrapped(new_u, new_v, new_lat, new_lon, cen_long, cone) |
||||
# |
||||
# # Note: The 2D version will return a 3D array with a 1 length |
||||
# # dimension. Numpy is unable to broadcast this without |
||||
# # sqeezing first. |
||||
# result = np.squeeze(result) |
||||
# |
||||
# output[left_and_slice_idxs] = result[:] |
||||
# |
||||
# return output |
||||
# |
||||
# return func_wrapper |
||||
|
||||
def uvmet_left_iter_nocopy(alg_dtype=np.float64): |
||||
"""Decorator to handle iterating over leftmost dimensions when using |
||||
multiple files and/or multiple times with the uvmet product. |
||||
|
||||
""" |
||||
@wrapt.decorator |
||||
def func_wrapper(wrapped, instance, args, kwargs): |
||||
u = args[0] |
||||
v = args[1] |
||||
lat = args[2] |
||||
lon = args[3] |
||||
cen_long = args[4] |
||||
cone = args[5] |
||||
|
||||
orig_dtype = u.dtype |
||||
|
||||
lat_lon_fixed = False |
||||
if lat.ndim == 2: |
||||
lat_lon_fixed = True |
||||
|
||||
if lon.ndim == 2 and not lat_lon_fixed: |
||||
raise ValueError("'lat' and 'lon' shape mismatch") |
||||
|
||||
num_left_dims_u = u.ndim - 2 |
||||
num_left_dims_lat = lat.ndim - 2 |
||||
|
||||
if (num_left_dims_lat > num_left_dims_u): |
||||
raise ValueError("number of 'lat' dimensions is greater than 'u'") |
||||
|
||||
if lat_lon_fixed: |
||||
mode = 0 # fixed lat/lon |
||||
else: |
||||
if num_left_dims_u == num_left_dims_lat: |
||||
mode = 1 # lat/lon same as u |
||||
else: |
||||
mode = 2 # probably 3D with 2D lat/lon plus Time |
||||
|
||||
has_missing = False |
||||
u_arr = u |
||||
if isinstance(u, DataArray): |
||||
u_arr = npvalues(u) |
||||
|
||||
v_arr = v |
||||
if isinstance(v, DataArray): |
||||
v_arr = npvalues(v) |
||||
|
||||
umissing = Constants.DEFAULT_FILL |
||||
if isinstance(u_arr, np.ma.MaskedArray): |
||||
has_missing = True |
||||
umissing = u_arr.fill_value |
||||
|
||||
vmissing = Constants.DEFAULT_FILL |
||||
if isinstance(v_arr, np.ma.MaskedArray): |
||||
has_missing = True |
||||
vmissing = v_arr.fill_value |
||||
|
||||
uvmetmissing = umissing |
||||
|
||||
is_stag = 0 |
||||
if (u.shape[-1] != lat.shape[-1] or u.shape[-2] != lat.shape[-2]): |
||||
is_stag = 1 |
||||
# Sanity check |
||||
if (v.shape[-1] == lat.shape[-1] or v.shape[-2] == lat.shape[-2]): |
||||
raise ValueError("u is staggered but v is not") |
||||
|
||||
if (v.shape[-1] != lat.shape[-1] or v.shape[-2] != lat.shape[-2]): |
||||
is_stag = 1 |
||||
# Sanity check |
||||
if (u.shape[-1] == lat.shape[-1] or u.shape[-2] == lat.shape[-2]): |
||||
raise ValueError("v is staggered but u is not") |
||||
|
||||
|
||||
|
||||
# No special left side iteration, return the function result |
||||
if (num_left_dims_u == 0): |
||||
return wrapped(u, v, lat, lon, cen_long, cone, isstag=is_stag, |
||||
has_missing=has_missing, umissing=umissing, |
||||
vmissing=vmissing, uvmetmissing=uvmetmissing) |
||||
|
||||
# Initial output is time,nz,2,ny,nx to create contiguous views |
||||
outdims = u.shape[0:num_left_dims_u] |
||||
extra_dims = tuple(outdims) # Copy the left-most dims for iteration |
||||
|
||||
outdims += (2,) |
||||
|
||||
outdims += lat.shape[-2:] |
||||
|
||||
outview_array = np.empty(outdims, alg_dtype) |
||||
|
||||
for left_idxs in iter_left_indexes(extra_dims): |
||||
left_and_slice_idxs = left_idxs + (slice(None),) |
||||
|
||||
if mode == 0: |
||||
lat_left_and_slice = (slice(None),) |
||||
elif mode == 1: |
||||
lat_left_and_slice = left_and_slice_idxs |
||||
elif mode == 2: |
||||
# Only need the left-most |
||||
lat_left_and_slice = tuple(left_idx |
||||
for left_idx in left_idxs[0:num_left_dims_lat]) |
||||
|
||||
|
||||
new_u = u[left_and_slice_idxs] |
||||
new_v = v[left_and_slice_idxs] |
||||
new_lat = lat[lat_left_and_slice] |
||||
new_lon = lon[lat_left_and_slice] |
||||
outview = outview_array[left_and_slice_idxs] |
||||
|
||||
# Call the numerical routine |
||||
result = wrapped(new_u, new_v, new_lat, new_lon, cen_long, cone, |
||||
isstag=is_stag, has_missing=has_missing, |
||||
umissing=umissing, vmissing=vmissing, |
||||
uvmetmissing=uvmetmissing, outview=outview) |
||||
|
||||
# Make sure the result is the same data as what got passed in |
||||
# Can delete this once everything works |
||||
if (result.__array_interface__["data"][0] != |
||||
outview.__array_interface__["data"][0]): |
||||
raise RuntimeError("output array was copied") |
||||
|
||||
# Need to reshape this so that u_v is left dim, then time (or others), |
||||
# then nz, ny, nz |
||||
output_dims = (2,) |
||||
output_dims += extra_dims |
||||
output_dims += lat.shape[-2:] |
||||
output = np.empty(output_dims, orig_dtype) |
||||
|
||||
output[0,:] = outview_array[...,0,:,:].astype(orig_dtype) |
||||
output[1,:] = outview_array[...,1,:,:].astype(orig_dtype) |
||||
|
||||
if has_missing: |
||||
output = np.ma.masked_values(output, uvmetmissing) |
||||
|
||||
return output |
||||
|
||||
return func_wrapper |
||||
|
@ -1,21 +1,26 @@
@@ -1,21 +1,26 @@
|
||||
import numpy as np |
||||
|
||||
import wrf._wrffortran |
||||
errlen = int(wrf._wrffortran.constants.errlen) |
||||
|
||||
|
||||
a = np.ones((3,3,3)) |
||||
b = np.zeros((3,3,3,3)) |
||||
c = np.zeros(errlen, "c") |
||||
errstat = np.array(0) |
||||
errmsg = np.zeros(512, "c") |
||||
errmsg = np.zeros(errlen, "c") |
||||
|
||||
c[:] = "Test String" |
||||
|
||||
|
||||
for i in xrange(2): |
||||
outview = b[i,:] |
||||
outview = outview.T |
||||
q = wrf._wrffortran.testfunc(a,outview,errstat=errstat,errstr=errmsg) |
||||
q[1,1,1] = 100 |
||||
q = wrf._wrffortran.testfunc(a,outview,c,errstat=errstat,errmsg=errmsg) |
||||
print errstat |
||||
|
||||
|
||||
|
||||
print errstat |
||||
print b |
||||
str_bytes = (bytes(c).decode("utf-8") for c in errmsg[:]) |
||||
str_bytes = (bytes(char).decode("utf-8") for char in errmsg[:]) |
||||
print repr(errmsg) |
||||
print "".join(str_bytes).strip() |
Loading…
Reference in new issue