diff --git a/.readthedocs.yml b/.readthedocs.yml index 4fe8d6300d..c48ee502d8 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -8,6 +8,10 @@ build: os: ubuntu-22.04 tools: python: "3.11" + jobs: + build: + html: + - sphinx-build -M html docs $READTHEDOCS_OUTPUT -j auto # Extra formats # PDF generation is failing for now; disabled on 2020-12-02 diff --git a/ac/makedep b/ac/makedep index 4903d88274..65a79044c0 100755 --- a/ac/makedep +++ b/ac/makedep @@ -389,7 +389,9 @@ def scan_fortran_file(src_file, defines=None): cpp_defines = defines if defines is not None else [] - cpp_macros = dict([t.split('=') for t in cpp_defines]) + cpp_macros = dict( + [t.split('=') if '=' in t else (t, None) for t in cpp_defines] + ) cpp_group_stack = [] with io.open(src_file, 'r', errors='replace') as file: @@ -454,13 +456,14 @@ def scan_fortran_file(src_file, defines=None): # Activate a new macro (ignoring the value) match = re_cpp_define.match(line) if match: + # TODO: Tokenize this, don't hunt for `(` in `macro`. tokens = line.strip()[1:].split(maxsplit=2) macro = tokens[1] value = tokens[2] if tokens[2:] else None if '(' in macro: # TODO: Actual handling of function macros macro, arg = macro.split('(', maxsplit=1) - value = '(' + arg + value + value = '(' + arg + value if value else '(' + arg cpp_macros[macro] = value # Deactivate a macro diff --git a/docs/.gitignore b/docs/.gitignore index e8b6a0513b..a9246f1e36 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -16,3 +16,8 @@ xml # Citation output bib*.aux citelist.doc* + + +# Python bytecode cache from the vendored _ext/autodoc_doxygen extension +__pycache__ +*.pyc diff --git a/docs/Doxyfile_rtd b/docs/Doxyfile_rtd index 479c03e1b4..b7c217769d 100644 --- a/docs/Doxyfile_rtd +++ b/docs/Doxyfile_rtd @@ -872,7 +872,8 @@ RECURSIVE = YES # Note that relative paths are relative to the directory from which doxygen is # run. -EXCLUDE = ../src/equation_of_state/TEOS10 +EXCLUDE = ../src/equation_of_state/TEOS10 \ + ../src/parameterizations/CVmix # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded @@ -1954,7 +1955,7 @@ XML_OUTPUT = xml # The default value is: YES. # This tag requires that the tag GENERATE_XML is set to YES. -XML_PROGRAMLISTING = NO +XML_PROGRAMLISTING = YES #--------------------------------------------------------------------------- # Configuration options related to the DOCBOOK output diff --git a/docs/Makefile b/docs/Makefile index 5729f73d90..69f4fddf1c 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -25,7 +25,7 @@ endif # Internal variables. PAPEROPT_a4 = -D latex_paper_size=a4 PAPEROPT_letter = -D latex_paper_size=letter -ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . +ALLSPHINXOPTS = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) # the i18n builder cannot share the environment and doctrees with the others I18NSPHINXOPTS = $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) . @@ -65,38 +65,45 @@ $(BUILDDIR): mkdir -p $@ html: $(DOXYGENBIN) $(BUILDDIR) - $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html + $(SPHINXBUILD) -M html . $(BUILDDIR) $(ALLSPHINXOPTS) -j auto +# Optional equation post-processing (was a patch in jr3cermak/sphinx@v3.2.1mom6.4 +# applied to sphinx/cmd/build.py); now invoked from the Makefile so we can run +# against stock upstream Sphinx. The nortd target below has the same hook. +ifeq ($(UPDATEHTMLEQS), Y) + @echo "Post processing equations." + @python3 ./postProcessEquations.py -d $(BUILDDIR) -p html -b sphinx -s index.html $(UPDATEHTMLEQSVERBOSE) +endif @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)/html." dirhtml: - $(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml + $(SPHINXBUILD) -M dirhtml . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml." singlehtml: - $(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml + $(SPHINXBUILD) -M singlehtml . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml." pickle: - $(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle + $(SPHINXBUILD) -M pickle . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished; now you can process the pickle files." json: - $(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json + $(SPHINXBUILD) -M json . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished; now you can process the JSON files." htmlhelp: - $(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp + $(SPHINXBUILD) -M htmlhelp . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished; now you can run HTML Help Workshop with the" \ ".hhp project file in $(BUILDDIR)/htmlhelp." qthelp: - $(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp + $(SPHINXBUILD) -M qthelp . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished; now you can run "qcollectiongenerator" with the" \ ".qhcp project file in $(BUILDDIR)/qthelp, like this:" @@ -105,7 +112,7 @@ qthelp: @echo "# assistant -collectionFile $(BUILDDIR)/qthelp/MOM6.qhc" devhelp: - $(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp + $(SPHINXBUILD) -M devhelp . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished." @echo "To view the help file:" @@ -114,48 +121,48 @@ devhelp: @echo "# devhelp" epub: - $(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub + $(SPHINXBUILD) -M epub . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The epub file is in $(BUILDDIR)/epub." latex: $(DOXYGENBIN) $(BUILDDIR) - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) -M latex . $(BUILDDIR) $(ALLSPHINXOPTS) @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)." latexpdf: $(DOXYGENBIN) $(BUILDDIR) - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) -M latex . $(BUILDDIR) $(ALLSPHINXOPTS) @echo "Running LaTeX files through pdflatex..." $(MAKE) -C $(BUILDDIR)/latex LATEXMKOPTS="-f -silent" all-pdf @echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex." latexpdfja: - $(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex + $(SPHINXBUILD) -M latex . $(BUILDDIR) $(ALLSPHINXOPTS) @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." text: - $(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text + $(SPHINXBUILD) -M text . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The text files are in $(BUILDDIR)/text." man: - $(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man + $(SPHINXBUILD) -M man . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The manual pages are in $(BUILDDIR)/man." texinfo: - $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + $(SPHINXBUILD) -M texinfo . $(BUILDDIR) $(ALLSPHINXOPTS) @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)." info: - $(SPHINXBUILD) -b texinfo $(ALLSPHINXOPTS) $(BUILDDIR)/texinfo + $(SPHINXBUILD) -M texinfo . $(BUILDDIR) $(ALLSPHINXOPTS) @echo "Running Texinfo files through makeinfo..." make -C $(BUILDDIR)/texinfo info @echo "makeinfo finished; the Info files are in $(BUILDDIR)/texinfo." @@ -166,28 +173,28 @@ gettext: @echo "Build finished. The message catalogs are in $(BUILDDIR)/locale." changes: - $(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes + $(SPHINXBUILD) -M changes . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "The overview file is in $(BUILDDIR)/changes." linkcheck: - $(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck + $(SPHINXBUILD) -M linkcheck . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Link check complete; look for any errors in the above output " \ "or in $(BUILDDIR)/linkcheck/output.txt." doctest: - $(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest + $(SPHINXBUILD) -M doctest . $(BUILDDIR) $(ALLSPHINXOPTS) @echo "Testing of doctests in the sources finished, look at the " \ "results in $(BUILDDIR)/doctest/output.txt." xml: - $(SPHINXBUILD) -b xml $(ALLSPHINXOPTS) $(BUILDDIR)/xml + $(SPHINXBUILD) -M xml . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The XML files are in $(BUILDDIR)/xml." pseudoxml: - $(SPHINXBUILD) -b pseudoxml $(ALLSPHINXOPTS) $(BUILDDIR)/pseudoxml + $(SPHINXBUILD) -M pseudoxml . $(BUILDDIR) $(ALLSPHINXOPTS) @echo @echo "Build finished. The pseudo-XML files are in $(BUILDDIR)/pseudoxml." diff --git a/docs/README.md b/docs/README.md index 1bec91f288..8ccf773aad 100644 --- a/docs/README.md +++ b/docs/README.md @@ -249,9 +249,11 @@ added to the RTD administrative interface if you want to trigger post processing If `UPDATEHTMLEQSVERBOSE` is set to `-v` this will turn on verbose printing for the post processor. -NOTE: These options affect solo doxygen html processing only for the `make nortd` option. For sphinx, they are -utilized in the sphinx python module to handle post processing and are not part of the Makefile. This was done -as RTD runs sphinx processing directly using `sphinx-build` and not the Makefile. +The post-processing hook is invoked from the `Makefile` for both the +`make nortd` (doxygen html) and `make html` (sphinx html) targets. RTD, +which runs `sphinx-build` directly without using this Makefile, will not +trigger the hook automatically; if equation post-processing is desired +on RTD it must be run as a separate step from the RTD build configuration. ##### PAPER @@ -316,18 +318,27 @@ pip install -r requirements.txt You may need to use `pip3` to install requirements for python3. -Requirements: -- sphinx -- sphinx-rtd-theme -- sphinx-bibtex -- sphinx-fortran -- sphinxcontrib\_autodox-doxygen -- flint -- lxml -- numpy -- future - -For machines that need to build future, numpy or lxml, these packages are required: +The full pinned set lives in `requirements.txt`. The current toolchain +uses stock upstream Sphinx 8.x from PyPI: + +- `sphinx>=8,<9` +- `sphinx-rtd-theme` +- `sphinxcontrib-bibtex` +- `lxml` (used by the vendored `_ext/autodoc_doxygen` extension to + parse the Doxygen XML output) +- `numpy` +- `sphinx-fortran` from the upstream `VACUMM/sphinx-fortran` repository, + pinned to a specific commit (upstream has not cut a PyPI release past + 1.1.1 but the master branch has continued fixes) +- `six`, still required at module load time by the pinned `sphinx-fortran` + commit + +The `sphinxcontrib-autodoc_doxygen` Sphinx extension that was previously +pulled in as a separate fork is now vendored in-tree under +`docs/_ext/autodoc_doxygen/`; nothing extra needs to be installed for it. + +For machines that need to build numpy or lxml from source, these packages +are required: - Cython - wheel @@ -342,8 +353,8 @@ PDF generation requires the following packages ### doxygen You may choose to download the [source](https://www.doxygen.nl/download.html). - -Latest is `doxygen-1.8.20.src.tar.gz`. +The example below uses `1.8.20` but you can substitute any compatible +release tarball. ```bash tar xzf doxygen-1.8.20.src.tar.gz @@ -356,15 +367,17 @@ sudo make install ``` The makefile for doxygen attempts to install the compiled version into /usr/local/bin. -You can link to a specific executable within the virtual environment. At this point we -also recommend renaming `doxygen` to `doxygen-1.8.20` within `/usr/local/bin`. +You can link to a specific executable within the virtual environment. -NOTE: The makefile for the documentation framework will attempt to compile a local doxygen -binary of version 1.8.13 if a binary cannot be found in the `$PATH`. +NOTE: If a doxygen binary is not found in `$PATH`, the documentation Makefile +will attempt to clone and compile its own copy of doxygen from source. The +default release pulled in this fallback path is set by `DOXYGEN_RELEASE` in +the Makefile (currently `Release_1_8_13`); pass a different value to override. #### Testing -A lot of manual testing has been completed using the following versions: +The toolchain is known to work with doxygen versions in the 1.8.x and 1.9.x +series. Older manual testing was performed against: * 1.8.13 * 1.8.14 * 1.8.19 @@ -376,7 +389,9 @@ The [Read the Docs](https://readthedocs.org/) (RTD) site uses a virtual machine (VM) for processing documentation. The VM architecture is type x86\_64. The default version for doxygen is 1.8.13 on the RTD VM. -NOTE: Using modified python modules on RTD is possible through careful crafting of the requirements.txt file. It is impossible to replace system binaries or compile code on RTD. It is possible to ship replacement binaries that can be run from the repo. For security reasons, a binary cannot be included in the MOM6 repository. +NOTE: It is impossible to replace system binaries or compile code on RTD. +It is possible to ship replacement binaries that can be run from the repo. +For security reasons, a binary cannot be included in the MOM6 repository. #### Logfiles @@ -387,19 +402,37 @@ Most websites force download of `*.log` files. # Credits +## 2026 +The documentation toolchain was modernized to run against stock upstream +Sphinx 8.x. The four-fork chain that the build had been carrying since +2020 was reduced to: + +- A vendored copy of the Doxygen-to-Sphinx bridge under + `docs/_ext/autodoc_doxygen/`, originally derived from the `0.7.13` + release of `jr3cermak/sphinxcontrib-autodoc_doxygen`. The vendored + version is ported to the Sphinx 8 API and lives in-tree where it can + be debugged and edited like any other project source. +- A pinned commit of upstream `VACUMM/sphinx-fortran` master. +- A small monkey-patch in `conf.py` for `sphinx.util.math.wrap_displaymath` + that suppresses Sphinx's default outer wrapping when the source already + supplies its own LaTeX environment, replacing the only functional change + the previous Sphinx fork carried. +- A small monkey-patch in `conf.py` for + `sphinxfortran.fortran_domain.FortranDomain.merge_domaindata` to fix a + parallel-build bug in upstream sphinx-fortran. To be removed when + upstream merges a fix. + +The `flint` dependency (formerly used to patch Doxygen's incomplete +parsing of Fortran functions with `result()` clauses) was dropped after +verifying empirically that nothing in the build pipeline references it. + ## 2020 -The documentation pipeline was upgraded by [Rob Cermak](https://github.com/jr3cermak) and [Marshall Ward](https://github.com/marshallward). Four modified python modules are required -to process the MOM6 documentation. The versions are tagged and placed into the production version of `requirements.txt`. Development versions may be found in the respective `dev` branches. - -| Source | Modified | Version | Development | -| ------ | -------- | ------- | ----------- | -| [sphinx](https://github.com/sphinx-doc/sphinx) | [sphinx-3.2.1mom6.4](https://github.com/jr3cermak/sphinx) | B:3.2.1mom6.4 | B:dev | -| [sphinxcontrib-autodoc-doxygen](https://github.com/rmcgibbo/sphinxcontrib-autodoc_doxygen) | [sphinxcontrib-autodoc-doxygen](https://github.com/jr3cermak/sphinxcontrib-autodoc_doxygen) | T:0.7.13 | B:dev | -| [sphinx-fortran](https://github.com/VACUMM/sphinx-fortran) | [sphinx-fortran](https://github.com/jr3cermak/sphinx-fortran) | T:1.2.2 | B:dev | -| [flint](https://github.com/marshallward/flint) | [flint](https://github.com/jr3cermak/flint) | T:0.0.1 | B:dev | -| [MOM6](https://github.com/NOAA-GFDL/MOM6) | [esmg-docs](https://github.com/ESMG/MOM6/tree/esmg-docs) | [esmg-docs](https://github.com/ESMG/MOM6/tree/esmg-docs) | B:[esmg-test](https://github.com/jr3cermak/MOM6/tree/esmg-test) | - -T: tag B: branch +The documentation pipeline was upgraded by [Rob Cermak](https://github.com/jr3cermak) +and [Marshall Ward](https://github.com/marshallward). The pipeline at that time +required four modified Python modules (`sphinx`, `sphinxcontrib-autodoc_doxygen`, +`sphinx-fortran`, `flint`), all forked under the `jr3cermak` GitHub account. +That toolchain was retired during the 2026 modernization above; this entry +is retained as historical context. ## 2017 The sphinx documentation of MOM6 is made possible by modifications by [Angus Gibson](https://github.com/angus-g) to two packages, [sphinx-fortran](https://github.com/angus-g/sphinx-fortran) and [autodoc\_doxygen](https://github.com/angus-g/sphinxcontrib-autodoc_doxygen). diff --git a/docs/_ext/autodoc_doxygen/__init__.py b/docs/_ext/autodoc_doxygen/__init__.py new file mode 100644 index 0000000000..7c4a3afc80 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/__init__.py @@ -0,0 +1,101 @@ +import os.path +from lxml import etree as ET +from sphinx.errors import ExtensionError + + +def set_doxygen_xml(app): + """Load all doxygen XML files from the app config variable + `app.config.doxygen_xml` which should be a path to a directory + containing doxygen xml output. If the configured path is relative, + it is resolved against `app.confdir` rather than the current working + directory -- Sphinx may have any cwd by the time builder-inited fires, + and in particular RTD runs sphinx-build from the repo root. + """ + doxygen_xml = app.config.doxygen_xml + if not os.path.isabs(doxygen_xml): + doxygen_xml = os.path.join(app.confdir, doxygen_xml) + + err = ExtensionError( + '[autodoc_doxygen] No doxygen ' + 'xml output found in doxygen_xml="%s"' % doxygen_xml) + + if not os.path.isdir(doxygen_xml): + raise err + + files = [os.path.join(doxygen_xml, f) + for f in os.listdir(doxygen_xml) + if f.lower().endswith('.xml') and not f.startswith('._')] + if len(files) == 0: + raise err + + setup.DOXYGEN_ROOT = ET.ElementTree(ET.Element('root')).getroot() + for file in files: + root = ET.parse(file).getroot() + for node in root: + setup.DOXYGEN_ROOT.append(node) + + +def get_doxygen_root(): + """Get the root element of the doxygen XML document. + """ + if not hasattr(setup, 'DOXYGEN_ROOT'): + setup.DOXYGEN_ROOT = ET.Element("root") # dummy + return setup.DOXYGEN_ROOT + + +def get_doxygen_id_index(): + """Return a dict mapping every ``@id`` in the merged doxygen tree to + the element that owns it. Built lazily on first use and memoized + on the :func:`setup` function object. + + Profiling a serial build at full MOM6 input (XML_PROGRAMLISTING=YES, + 109 MB merged tree) showed ``xmlutils.visit_ref`` burning 250 s of + self time -- 27% of total wall clock -- in a single ``findall('.//* + [@id=X]')`` call that linearly scanned the entire merged tree once + per ```` in prose. This index turns that scan into an O(1) + dict lookup. Same shape of fix as the scanNode `//` -> `.//` patch + in commit 8a217135e. + """ + if not hasattr(setup, 'DOXYGEN_ID_INDEX'): + root = get_doxygen_root() + index = {} + for el in root.iter(): + eid = el.get('id') + if eid is not None: + index[eid] = el + setup.DOXYGEN_ID_INDEX = index + return setup.DOXYGEN_ID_INDEX + + +def setup(app): + import sphinx + from .autodoc import ( + DoxygenMethodDocumenter, + DoxygenTypeDocumenter, + DoxygenModuleDocumenter, + ) + from .autosummary import DoxygenAutosummary, DoxygenAutoEnum + from .autosummary.generate import process_generate_options + from .autodoxysource import AutoDoxySourceDirective + + app.connect("builder-inited", set_doxygen_xml) + app.connect("builder-inited", process_generate_options) + + app.setup_extension('sphinx.ext.autodoc') + app.setup_extension('sphinx.ext.autosummary') + + app.add_autodocumenter(DoxygenModuleDocumenter) + app.add_autodocumenter(DoxygenMethodDocumenter) + app.add_autodocumenter(DoxygenTypeDocumenter) + + app.add_config_value("doxygen_xml", "", 'env') + # Used in autodoc_doxygen/autosummary/generate.py + app.add_config_value('autosummary_toctree', '', 'html') + + app.add_directive('autodoxysummary', DoxygenAutosummary) + app.add_directive('autodoxyenum', DoxygenAutoEnum) + app.add_directive('autodoxysource', AutoDoxySourceDirective) + + app.add_css_file('autodoxysource.css') + + return {'version': sphinx.__display_version__, 'parallel_read_safe': True} diff --git a/docs/_ext/autodoc_doxygen/autodoc.py b/docs/_ext/autodoc_doxygen/autodoc.py new file mode 100644 index 0000000000..1b62eeca47 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autodoc.py @@ -0,0 +1,610 @@ +import re + +from docutils.parsers.rst import directives +from lxml import etree as ET +from sphinx.ext.autodoc import Documenter, members_option, ALL +from sphinx.errors import ExtensionError +from sphinx.util import logging + +from . import get_doxygen_root +from .autodoxysource import get_source_link +from .xmlutils import format_xml_paragraph, flatten + +logger = logging.getLogger(__name__) + + +class DoxygenDocumenter(Documenter): + # Variables to store the names of the object being documented. modname and fullname are redundant, + # and objpath is always the empty list. This is inelegant, but we need to work with the superclass. + + fullname = None # example: "OpenMM::NonbondedForce" or "OpenMM::NonbondedForce::methodName"" + modname = None # example: "OpenMM::NonbondedForce" or "OpenMM::NonbondedForce::methodName"" + objname = None # example: "NonbondedForce" or "methodName" + objpath = [] # always the empty list + object = None # the xml node for the object + # This allows section headers in autogenerated content + titles_allowed = True + + option_spec = { + 'members': members_option, + } + + # original + #def __init__(self, directive, name, indent=u'', id=None): + # super(DoxygenDocumenter, self).__init__(directive, name, indent) + # if id is not None: + # self.parse_id(id) + # new + def __init__(self, directive, name, indent=u'', id=None, brief=False, parent=None): + super().__init__(directive, name, indent) + + self.parent = parent + if id is not None: + self.parse_id(id) + self.brief = brief + + def parse_id(self, id): + return False + + def parse_name(self): + """Determine what module to import and what attribute to document. + Returns True and sets *self.modname*, *self.objname*, *self.fullname*, + if parsing and resolving was successful. + """ + # To view the context and order in which all of these methods get called, + # See, Documenter.generate(). That's the main "entry point" that first + # calls parse_name(), follwed by import_object(), format_signature(), + # add_directive_header(), and then add_content() (which calls get_doc()) + + # methods in the superclass sometimes use '.' to join namespace/class + # names with method names, and we don't want that. + self.name = self.name.replace('.', '::') + self.fullname = self.name + self.modname = self.fullname + self.objpath = [] + + if '::' in self.name: + parts = self.name.split('::') + self.objname = parts[-1] + else: + self.objname = self.name + + return True + + def document_members(self, all_members=False): + """Generate reST for member documentation. + If *all_members* is True, do all members, else those given by + *self.options.members*. + """ + want_all = all_members or self.options.inherited_members or \ + self.options.members is ALL + # new + members = all_members + # find out which members are documentable + # removed + #members_check_module, members = self.get_object_members(want_all) + + # remove members given by exclude-members + if self.options.exclude_members: + members = [(membername, member) for (membername, member) in members + if membername not in self.options.exclude_members] + + # document non-skipped members + memberdocumenters = [] + for (mname, member, isattr) in self.filter_members(members, want_all): + classes = [cls for cls in self.env.app.registry.documenters.values() + if cls.can_document_member(member, mname, isattr, self)] + if not classes: + # don't know how to document this member + continue + + # prefer the documenter with the highest priority + classes.sort(key=lambda cls: cls.priority) + + # change + #documenter = classes[-1](self.directive, mname, indent=self.indent, id=member.get('id')) + documenter = classes[-1](self.directive, mname, indent=self.indent, + id=member.get('id'), brief=self.brief, + parent=self.object) + memberdocumenters.append((documenter, isattr)) + + for documenter, isattr in memberdocumenters: + documenter.generate( + all_members=True, real_modname=self.real_modname, + #check_module=members_check_module and not isattr) + # modified + check_module=False and not isattr) + + # reset current objects + self.env.temp_data['autodoc:module'] = None + self.env.temp_data['autodoc:class'] = None + +# Copy of DoxygenClassDocumenter -> DoxygenModuleDocumenter +class DoxygenModuleDocumenter(DoxygenDocumenter): + # change + #objtype = 'doxyclass' + #directivetype = 'class' + #domain = 'cpp' + objtype = 'doxymodule' + directivetype = 'module' + domain = 'f' + priority = 100 + + option_spec = { + 'members': members_option, + 'methods': directives.flag, + 'types': directives.flag, + } + + @classmethod + def can_document_member(cls, member, membername, isattr, parent): + # this method is only called from Documenter.document_members + # when a higher level documenter (module or namespace) is trying + # to choose the appropriate documenter for each of its lower-level + # members. Currently not implemented since we don't have a higher-level + # doumenter like a DoxygenNamespaceDocumenter. + return False + + def import_object(self): + """Import the object and set it as *self.object*. In the call sequence, this + is executed right after parse_name(), so it can use *self.fullname*, *self.objname*, + and *self.modname*. + + Returns True if successful, False if an error occurred. + """ + # change + #xpath_query = './/compoundname[text()="%s"]/..' % self.fullname + xpath_query = './compounddef/compoundname[text()="%s"]/..' % self.fullname + match = get_doxygen_root().xpath(xpath_query) + if len(match) != 1: + # change + #raise ExtensionError('[autodoc_doxygen] could not find class (fullname="%s"). I tried' + raise ExtensionError('[autodoc_doxygen] could not find module (fullname="%s"). I tried' + 'the following xpath: "%s"' % (self.fullname, xpath_query)) + + self.object = match[0] + if self.env.app.verbosity > 0: print("[debug] xpath(%s) match(%s)" % (xpath_query,match[0].items())) + return True + + # todo: typo: report upstream + #def format_signaure(self): + def format_signature(self): + if self.env.app.verbosity > 0: print("[debug] DoxygenModuleDocumenter format_signature called") + #return '' + + def format_name(self): + return self.fullname + + # change + #def get_doc(self): + # detaileddescription = self.object.find('detaileddescription') + # doc = [format_xml_paragraph(detaileddescription)] + # encoding depricated? + # called via add_content + #def get_doc(self, encoding): + def get_doc(self): + if self.brief: + description = self.object.find('briefdescription') + else: + description = self.object.find('detaileddescription') + # use the brief description if there's no content in the + # detailed description + if not len(description) and not description.text.strip(): + description = self.object.find('briefdescription') + + #if self.env.app.verbosity > 0: print("[debug] get_doc(%s)(%s)" % (self.brief, description.items())) + doc = [format_xml_paragraph(description, self.env.config.sphinx_build_mode, + verbosity=self.env.app.verbosity)] + + #if self.env.app.verbosity > 0: + # if self.name == 'mom_ice_shelf': + + if not any(len(d.strip()) for d in doc[0]): + doc.append(['', '']) + + if self.brief: + # More references need to be unique across all pages + doc.append(['`More... `_' % (self.name), '']) + + return doc + + def get_object_members(self, want_all): + # change + pass + #all_members = self.object.xpath('.//sectiondef[@kind="public-func" ' + # 'or @kind="public-static-func"]/memberdef[@kind="function"]') + # + #if want_all: + # return False, ((m.find('name').text, m) for m in all_members) + #else: + # if not self.options.members: + # return False, [] + # else: + # return False, ((m.find('name').text, m) for m in all_members + # if m.find('name').text in self.options.members) + + def filter_members(self, members, want_all): + ret = [] + for (membername, member) in members: + ret.append((membername, member, False)) + return ret + + # change function arguments + #def document_members(self, all_members=False): + def document_members(self, member_type, all_members=False): + if member_type == 'func': + all_members = self.object.xpath('./sectiondef[@kind="func" ' + 'or @kind="public-static-func"]/memberdef[@kind="function"]') + + members = [(m.find('name').text, m) for m in all_members] + + elif member_type == 'type': + classes = self.object.findall('./innerclass') + members = [] + for c in classes: + + class_obj = get_doxygen_root().find('./compounddef[@id="%s"]' % c.get('refid')) + if class_obj.get('kind') == 'type': + members.append((class_obj.find('compoundname').text, class_obj)) + + if self.env.app.verbosity > 0: print("[debug] members(%s)" % (members)) + # Calls the plain DoxygenDocumenter to generate the rst + super().document_members(all_members=members) + # change + # super(DoxygenClassDocumenter, self).document_members(all_members=all_members) + # Uncomment to view the generated rst for the class. + # print('\n'.join(self.directive.result)) + + # added functions + + def add_title(self, title, char='='): + sourcename = self.get_sourcename() + + self.add_line(u'', sourcename) + self.add_line(char * len(title), sourcename) + self.add_line(title, sourcename) + self.add_line(char * len(title), sourcename) + self.add_line(u'', sourcename) + + # This generates the autogenerated content for all the module + # functions + def generate(self, more_content=None, real_modname=None, + check_module=False, all_members=False): + if not self.parse_name(): + logger.warning("don't know which module to import for autodocumenting %r" % self.name) + return + + if not self.import_object(): + return + + self.real_modname = real_modname or self.get_real_modname() + + # we can't import anything, since we're not Python + self.analyzer = None + + if check_module and not self.check_module(): + return + + sourcename = self.get_sourcename() + + # add title + title = '%s module reference' % self.format_name() + if self.env.app.verbosity > 0: + print("[debug] add_title:%s" % (title)) + self.add_title(title, char='=') + + # module directive + self.add_line(u'.. f:module:: %s' % self.format_name(), sourcename) + self.add_line(u'', sourcename) + + # brief description + self.brief = True + self.add_content(more_content) + + # we want a brief description of types/functions here + + # detailed description + #self.add_line(u'.. _`More...`:', sourcename) + self.add_line('','') + self.add_line(u'.. _DETA%s:' % (self.name), '') + self.add_title('Detailed Description', char='-') + self.brief = False + self.add_content(None) + + if 'types' in self.options: + self.add_title('Type Documentation', char='-') + self.document_members('type', all_members) + + # member doc + if 'methods' in self.options: + self.add_title('Function/Subroutine Documentation', char='-') + self.document_members('func', all_members) + + # [source] link at the bottom of the module page + src = get_source_link(self.object) + if src: + docname, line = src + file_id = docname.rsplit('/', 1)[-1] + self.add_line(u'', sourcename) + self.add_line( + u'`[source] <../source/%s.html#L%s>`__' % (file_id, line), + sourcename) + + if self.env.app.verbosity > 0: + if self.real_modname == 'mom_eos': + # Result: self.directive.result.data + pass + +class DoxygenMethodDocumenter(DoxygenDocumenter): + objtype = 'doxymethod' + directivetype = 'function' + # See if this can be modified to work for cpp and f? + # Put a breakpoint here + #domain = 'cpp' + domain = 'f' + priority = 100 + + @classmethod + def can_document_member(cls, member, membername, isattr, parent): + if ET.iselement(member) and member.tag == 'memberdef' and member.get('kind') == 'function': + return True + return False + + # add function + def add_directive_header(self, sig): + """Add the directive header and options to the generated content.""" + domain = self.domain + # use the field, without other information (e.g. public) + # functions are badly formed out of this function + if self.env.app.verbosity > 0: + print("[debug] add_directive_header sig(%s)" % (sig)) + if sig is not None and sig.find('adv_dyn') >= 0: + pass + name = self.format_name() + sourcename = self.get_sourcename() + + # determine which directive to use from the typefield + typefield = self.get_typefield() + if 'subroutine' in typefield: + directive = 'subroutine' + else: + directive = 'function' + + if self.env.app.verbosity > 0: + print("[debug] DoxygenMethodDocumenter directive(%s) name(%s)" % (directive,name)) + if name.find('eos_domain') >= 0 or name.find('mom_state_is') >= 0: + pass + + self.add_line(u'.. %s:%s:: %s%s' % (domain, directive, name, sig), + sourcename) + + def parse_id(self, id): + # added + # try to search our parent node instead of the entire tree + parent = self.parent + if parent is None: + parent = get_doxygen_root() + + xp = './/*[@id="%s"]' % id + # original + #match = get_doxygen_root().xpath(xp) + match = parent.xpath(xp) + if len(match) > 0: + match = match[0] + self.fullname = match.find('./definition').text.split()[-1] + self.modname = self.fullname + self.objname = match.find('./name').text + self.object = match + return False + + def import_object(self): + if ET.iselement(self.object): + # self.object already set from DoxygenDocumenter.parse_name(), + # caused by passing in the `id` of the node instead of just a + # classname or method name + return True + + return False + # original + xpath_query = ('.//compoundname[text()="%s"]/../sectiondef[@kind="public-func"]' + '/memberdef[@kind="function"]/name[text()="%s"]/..') % tuple(self.fullname.rsplit('::', 1)) + match = get_doxygen_root().xpath(xpath_query) + if len(match) == 0: + raise ExtensionError('[autodoc_doxygen] could not find method (modname="%s", objname="%s"). I tried ' + 'the following xpath: "%s"' % (tuple(self.fullname.rsplit('::', 1)) + (xpath_query,))) + self.object = match[0] + + def get_doc(self): + doc = [format_xml_paragraph(self.object.find('briefdescription'), self.env.config.sphinx_build_mode, + verbosity=self.env.app.verbosity)] + # debug + if self.object.find('name').text == 'eos_domain': + pass + #detaileddescription = self.object.find('detaileddescription') + #doc = [format_xml_paragraph(detaileddescription,self.env.config.sphinx_build_mode)] + + # add parameter documentation (in detaileddescription) for main function documentation + if not self.brief: + doc += [format_xml_paragraph(self.object.find('detaileddescription'), self.env.config.sphinx_build_mode, + verbosity=self.env.app.verbosity)] + + if self.object.find('name').text == 'eos_domain': + pass + # File location + # ffile = self.object.find('location').get('file') + # add references/referencedby + references = self.object.findall('references') + for ref in references: + name = ref.text + doc.append([':callto: :f:func:`%s <%s>`' % (name, name.split('::')[-1])]) + referencedby = self.object.findall('referencedby') + for ref in referencedby: + name = ref.text + doc.append([':calledfrom: :f:func:`%s <%s>`' % (name, name.split('::')[-1])]) + + # If a document returns a string with :return: that is our signal to use flint to + # repair the documentation a bit. Doxygen cannot parse functions of this syntax: + #!> This subroutine returns a two point integer array indicating the domain of i-indices + #!! to work on in EOS calls based on information from a hor_index type + #function EOS_domain(HI, halo) result(EOSdom) + # type(hor_index_type), intent(in) :: HI !< The horizontal index structure + # integer, optional, intent(in) :: halo !< The halo size to work on; missing is equivalent to 0. + # integer, dimension(2) :: EOSdom !< The index domain that the EOS will work on, taking into account + # !! that the arrays inside the EOS routines will start at 1. + # + # ! Local variables + # integer :: halo_sz + # + # halo_sz = 0 ; if (present(halo)) halo_sz = halo + # + # EOSdom(1) = HI%isc - (HI%isd-1) - halo_sz + # EOSdom(2) = HI%iec - (HI%isd-1) + halo_sz + # + #end function EOS_domain + #if self.object.find('name').text == 'eos_domain': + + if not self.brief: + src = get_source_link(self.object) + if src: + docname, line = src + file_id = docname.rsplit('/', 1)[-1] + doc.append(['`[source] <../source/%s.html#L%s>`__' % (file_id, line)]) + + return doc + + # added function + def get_typefield(self): + return ' '.join(self.object.find('definition').text.split()[:-1]) + + def format_name(self): + #def text(el): + # if el.text is not None: + # return el.text + # return '' + + #def tail(el): + # if el.tail is not None: + # return el.tail + # return '' + + #rtype_el = self.object.find('type') + #rtype_el_ref = rtype_el.find('ref') + #if rtype_el_ref is not None: + # rtype = text(rtype_el) + text(rtype_el_ref) + tail(rtype_el_ref) + #else: + # rtype = rtype_el.text + + # replaced above with + + # we just want to get the bare part of the "type" field + # i.e. subroutine or function + typefield = self.get_typefield() + + if typefield is None: + rtype = None + elif 'function' in typefield: + # Don't include the return type in the directive name. + # sphinx-fortran's f_sig_re regex only matches + # "function name(args)" or bare "name(args)"; a return + # type prefix like "real name(args)" causes the regex + # to fail, leaving the function unregistered and + # unclickable. The f:function directive already labels + # the entry as "function". + rtype = None + else: + rtype = 'subroutine' if 'subroutine' in typefield else 'unknown' + + signame = (rtype and (rtype + ' ') or '') + self.objname + return self.format_template_name() + signame + + def format_template_name(self): + types = [e.text for e in self.object.findall('templateparamlist/param/type')] + if len(types) == 0: + return '' + return 'template <%s>\n' % ','.join(types) + + def format_signature(self): + args = self.object.find('argsstring').text + if self.env.app.verbosity > 0: + print("[debug] DoxygenMethodDocumenter format_signature called (%s)" % (args)) + if args is not None and args.find('adv_dyn')>=0: + pass + return args + + def document_members(self, all_members=False): + pass + +# add class + +class DoxygenTypeDocumenter(DoxygenDocumenter): + objtype = 'doxytype' + directivetype = 'type' + domain = 'f' + priority = 100 + + @classmethod + def can_document_member(cls, member, membername, isattr, parent): + if ET.iselement(member) and member.tag == 'compounddef' and member.get('kind') == 'type': + return True + return False + + def import_object(self): + if ET.iselement(self.object): + return True + return False + + def parse_id(self, id): + self.object = get_doxygen_root().find('./compounddef[@id="%s"]' % id) + self.fullname = self.object.find('compoundname').text + self.modname, self.objname = self.fullname.rsplit('::') + + return False + + def format_name(self): + return self.objname + + def add_directive_header(self, sig): + """Add the directive header and options to the generated content.""" + domain = self.domain + directive = 'type' + name = self.format_name() + sourcename = self.get_sourcename() + + self.add_line(u'.. %s:%s:: %s' % (domain, directive, name), + sourcename) + + def get_doc(self): + desc = [format_xml_paragraph(self.object.find('briefdescription'), + self.env.config.sphinx_build_mode, verbosity=self.env.app.verbosity)] + + for member in self.object.findall('./sectiondef/memberdef'): + name = member.find('name').text + type_el = member.find('type') + full_type = ''.join(type_el.itertext()).strip() if type_el is not None else '' + + if member.get('prot') == 'private': + if full_type: + full_type += ', private' + else: + full_type = 'private' + + field = ':typefield %s:' % name + if full_type: + field += ' ``%s``' % full_type + + brief = member.find('briefdescription/para') + if brief is not None: + field += ' ' + brief.text + + desc.append([field]) + + src = get_source_link(self.object) + if src: + docname, line = src + file_id = docname.rsplit('/', 1)[-1] + desc.append(['`[source] <../source/%s.html#L%s>`__' % (file_id, line)]) + + return desc + + def document_members(self, all_members=False): + pass diff --git a/docs/_ext/autodoc_doxygen/autodoxysource.py b/docs/_ext/autodoc_doxygen/autodoxysource.py new file mode 100644 index 0000000000..a12af14b11 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autodoxysource.py @@ -0,0 +1,231 @@ +"""autodoxysource directive: render doxygen as +syntax-highlighted source with clickable cross-references.""" + +from docutils import nodes +from docutils.parsers.rst import Directive +from sphinx import addnodes +from sphinx.util import logging + +from . import get_doxygen_root + +logger = logging.getLogger(__name__) + +# Map doxygen highlight classes to CSS classes for the source listing. +_HL_CSS = { + 'comment': 'f-hl-comment', + 'normal': 'f-hl-normal', + 'keyword': 'f-hl-keyword', + 'keywordtype': 'f-hl-keywordtype', + 'keywordflow': 'f-hl-keywordflow', + 'stringliteral': 'f-hl-stringliteral', + 'preprocessor': 'f-hl-preprocessor', +} + +# Lazily-built index: doxygen refid -> (refdomain, reftype, reftarget) +_REFID_INDEX = None + +# Index: absolute file path -> doxygen file-ID (for [source] links) +_FILE_PATH_INDEX = None + + +def _build_refid_index(): + """Walk the merged doxygen tree once, building refid -> xref info.""" + global _REFID_INDEX + _REFID_INDEX = {} + root = get_doxygen_root() + + for cd in root.findall('./compounddef'): + cid = cd.get('id') + kind = cd.get('kind') + cn = cd.find('compoundname') + if cn is None or cn.text is None: + continue + name = cn.text + + if kind == 'namespace': + _REFID_INDEX[cid] = ('f', 'mod', name) + elif kind in ('type', 'struct'): + # compoundname is "parent::typename" -> "parent/typename" + _REFID_INDEX[cid] = ('f', 'type', name.replace('::', '/')) + elif kind == 'interface': + _REFID_INDEX[cid] = ('f', 'func', name.replace('::', '/')) + # kind='file' / 'page' / 'dir' -> no xref + + # Index memberdefs inside this compound + for md in cd.iter('memberdef'): + mid = md.get('id') + mk = md.get('kind') + mn_el = md.find('name') + if mn_el is None or mn_el.text is None: + continue + qualified = name + '/' + mn_el.text + if mk == 'function': + _REFID_INDEX[mid] = ('f', 'func', qualified) + # variables / enums aren't useful xref targets here + + +def _resolve_ref(refid): + """O(1) lookup of a doxygen refid to (refdomain, reftype, reftarget).""" + global _REFID_INDEX + if _REFID_INDEX is None: + _build_refid_index() + return _REFID_INDEX.get(refid) + + +def _build_file_path_index(): + """Build filepath -> file_id index from compounddef[@kind='file'].""" + global _FILE_PATH_INDEX + _FILE_PATH_INDEX = {} + root = get_doxygen_root() + for cd in root.findall('./compounddef[@kind="file"]'): + fid = cd.get('id') + loc = cd.find('location') + if loc is not None and loc.get('file'): + _FILE_PATH_INDEX[loc.get('file')] = fid + + +def get_source_link(xml_node): + """Given a doxygen XML node (compounddef or memberdef), return + (docname, line) for a [source] link, or None if unavailable. + + *docname* is the Sphinx document name (e.g. 'api/generated/source/MOM_8F90'). + *line* is the source line number string. + """ + global _FILE_PATH_INDEX + if _FILE_PATH_INDEX is None: + _build_file_path_index() + + loc = xml_node.find('location') + if loc is None: + return None + + filepath = loc.get('file') + line = loc.get('line', '1') + if not filepath: + return None + + file_id = _FILE_PATH_INDEX.get(filepath) + if file_id is None: + return None + + docname = 'api/generated/source/' + file_id + return (docname, line) + + +class AutoDoxySourceDirective(Directive): + """.. autodoxysource:: + + Render the source listing for a doxygen file compound, with + per-line anchors, syntax highlighting, and clickable identifiers + that link to the Sphinx API documentation. + """ + required_arguments = 1 + optional_arguments = 0 + has_content = False + + def run(self): + file_id = self.arguments[0] + root = get_doxygen_root() + + compounddef = root.find('./compounddef[@id="%s"]' % file_id) + if compounddef is None: + logger.warning('autodoxysource: compounddef not found for %s', + file_id) + return [nodes.paragraph('', 'Source listing not available.')] + + programlisting = compounddef.find('.//programlisting') + if programlisting is None: + logger.warning('autodoxysource: no programlisting in %s', file_id) + return [nodes.paragraph('', 'Source listing not available.')] + + # Outer wrapper. support_smartquotes=False propagates to all + # descendants via sphinx.util.nodes.is_smartquotable, so every + # text node inside the listing is skipped by the smartquotes + # transform -- a meaningful speedup on 329 large source pages. + table = nodes.container(classes=['autodoxysource']) + table['support_smartquotes'] = False + + for codeline in programlisting.findall('codeline'): + lineno = codeline.get('lineno', '') + + # Per-line container carries the #L anchor via ``ids``; + # the previous standalone ``target`` node was dropped to + # cut a node per line. The inner ``source-code`` inline + # wrapper is kept because Sphinx's HTML writer asserts + # that ``reference`` nodes (produced when pending_xrefs + # resolve) have a ``TextElement`` parent. + line_node = nodes.container(classes=['source-line']) + if lineno: + line_node['ids'] = ['L' + lineno] + + line_node += nodes.inline(lineno, lineno, + classes=['source-lineno']) + + code_node = nodes.inline(classes=['source-code']) + for hl in codeline: + if hl.tag != 'highlight': + continue + css = _HL_CSS.get(hl.get('class', 'normal'), 'f-hl-normal') + _walk_highlight(hl, css, code_node) + line_node += code_node + + table += line_node + + return [table] + + +def _walk_highlight(hl, css_class, parent): + """Walk mixed content of a element, appending nodes to + *parent*. + + Text fragments are coalesced: consecutive characters with the same + CSS class (including ```` expansions and tail text) are + flushed as a single ``inline`` node rather than one per fragment. + Before coalescing, a typical line emitted 10-30 nodes; after, it + emits closer to 3-5. Node count drives pickling, transform walks, + and HTML writing cost linearly for the source-listing pages. + """ + buf = [] + + def flush(): + if buf: + text = ''.join(buf) + parent.append(nodes.inline(text, text, classes=[css_class])) + buf.clear() + + if hl.text: + buf.append(hl.text) + + for child in hl: + if child.tag == 'sp': + buf.append(' ') + elif child.tag == 'ref': + flush() + _emit_ref(child, css_class, parent) + elif child.text: + buf.append(child.text) + + if child.tail: + buf.append(child.tail) + + flush() + + +def _emit_ref(ref_el, css_class, parent): + """Emit a pending_xref (or plain text fallback) for a element.""" + ref_text = ref_el.text or '' + refid = ref_el.get('refid', '') + + resolved = _resolve_ref(refid) + if resolved: + refdomain, reftype, reftarget = resolved + inner = nodes.inline(ref_text, ref_text, classes=[css_class]) + xref = addnodes.pending_xref( + '', inner, + refdomain=refdomain, + reftype=reftype, + reftarget=reftarget, + ) + parent += xref + else: + parent += nodes.inline(ref_text, ref_text, classes=[css_class]) diff --git a/docs/_ext/autodoc_doxygen/autosummary/__init__.py b/docs/_ext/autodoc_doxygen/autosummary/__init__.py new file mode 100644 index 0000000000..fec98be964 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autosummary/__init__.py @@ -0,0 +1,260 @@ +import re +import operator +from functools import reduce +from itertools import count, groupby + +from docutils import nodes +from docutils.parsers.rst import directives +from docutils.statemachine import StringList, ViewList +from sphinx import addnodes +from sphinx.ext.autosummary import Autosummary, autosummary_table +from sphinx.util import logging + +from .. import get_doxygen_root +from ..autodoc import DoxygenMethodDocumenter, DoxygenModuleDocumenter +from ..xmlutils import format_xml_paragraph + +logger = logging.getLogger(__name__) + + +def import_by_name(name, env=None, prefixes=None, i=0): + """Get xml documentation for a class/method with a given name. + If there are multiple classes or methods with that name, you + can use the `i` kwarg to pick which one. + """ + if prefixes is None: + prefixes = [None] + + # angus + if env is not None: + parent = env.ref_context.get('cpp:parent_symbol') + parent_symbols = [] + while parent is not None and parent.identifier is not None: + parent_symbols.insert(0, str(parent.identifier)) + parent = parent.parent + prefixes.append('::'.join(parent_symbols)) + + # orig + if env is not None: + #if env.ref_context != None and env.ref_context != {}: + # print("[debug] parents: %s" % env.ref_context) + parents = env.ref_context.get('cpp:parent_key') + if parents is not None: + parent_symbols = [p[0].get_display_string() for p in parents] + prefixes.append('::'.join(parent_symbols)) + + # unmodified + tried = [] + for prefix in prefixes: + try: + if prefix: + prefixed_name = '::'.join([prefix, name]) + else: + prefixed_name = name + return _import_by_name(prefixed_name, i=i) + except ImportError: + tried.append(prefixed_name) + raise ImportError('no module named %s' % ' or '.join(tried)) + +def _import_by_name(name, i=0): + root = get_doxygen_root() + name = name.replace('.', '::') + + if '::' in name: + xpath_query = ( + './compounddef/compoundname[text()="%s"]/../' + 'sectiondef[@kind="func"]/memberdef[@kind="function"]/' + 'name[text()="%s"]/..') % tuple(name.rsplit('::', 1)) + m = root.xpath(xpath_query) + if len(m) > 0: + obj = m[i] + full_name = '.'.join(name.rsplit('::', 1)) + return full_name, obj, full_name, '' + + xpath_query = ('./compounddef/compoundname[text()="%s"]/..' % name) + m = root.xpath(xpath_query) + if len(m) > 0: + obj = m[i] + return (name, obj, name, '') + + raise ImportError() + +def get_documenter(obj, full_name): + if obj.tag == 'memberdef' and obj.get('kind') == 'function': + return DoxygenMethodDocumenter + elif obj.tag == 'compounddef': + return DoxygenModuleDocumenter + + raise NotImplementedError(obj.tag) + + +class DoxygenAutosummary(Autosummary): + # add + option_spec = { + 'toctree': directives.unchanged, + 'nosignatures': directives.flag, + 'template': directives.unchanged, + 'kind': directives.unchanged, + 'generate': directives.flag, + } + + def get_items(self, names): + """Try to import the given names, and return a list of + ``[(name, signature, summary_string, real_name), ...]``. + """ + env = self.state.document.settings.env + items = [] + + # Add "generate" directive + if 'generate' in self.options: + # don't generate a summary for pages + if self.options['kind'] == 'page': + return [] + + modules = get_doxygen_root().xpath('./compound[@kind="namespace"]') + names = [m.find('name').text for m in modules] + + # TODO: silently fail when there are no fortran files provided? + try: + names_and_counts = reduce(operator.add, + [tuple(zip(g, count())) for _, g in groupby(names)]) # type: List[(Str, Int)] + except: + return items + + for name, i in names_and_counts: + display_name = name + if name.startswith('~'): + name = name[1:] + display_name = name.split('::')[-1] + + try: + real_name, obj, parent, modname = import_by_name(name, env=env, i=i) + except ImportError: + logger.warning('failed to import %s' % name) + items.append((name, '', '', name)) + continue + + self.bridge.result = StringList() # initialize for each documenter + documenter = get_documenter(obj, parent)(self.bridge, real_name, + id=obj.get('id'), + brief=True, + parent=obj.find('..')) + if not documenter.parse_name(): + logger.warning('failed to parse name %s' % real_name) + items.append((display_name, '', '', real_name)) + continue + if not documenter.import_object(): + logger.warning('failed to import object %s' % real_name) + items.append((display_name, '', '', real_name)) + continue + if documenter.options.members and not documenter.check_module(): + continue + # -- Grab the signature + sig = documenter.format_signature() + + # -- Grab the summary + documenter.add_content(None) + doc = list(documenter.process_doc([self.bridge.result.data])) + + while doc and not doc[0].strip(): + doc.pop(0) + + # If there's a blank line, then we can assume the first sentence / + # paragraph has ended, so anything after shouldn't be part of the + # summary + for i, piece in enumerate(doc): + if not piece.strip(): + doc = doc[:i] + break + + # Try to find the "first sentence", which may span multiple lines + m = re.search(r"^([A-Z].*?\.)(?:\s|$)", " ".join(doc).strip()) + if m: + summary = m.group(1).strip() + elif doc: + summary = doc[0].strip() + else: + summary = '' + + items.append((display_name, sig, summary, real_name)) + + return items + + def get_tablespec(self): + table_spec = addnodes.tabular_col_spec() + table_spec['spec'] = 'll' + + table = autosummary_table('') + real_table = nodes.table('', classes=['longtable']) + table.append(real_table) + group = nodes.tgroup('', cols=2) + real_table.append(group) + group.append(nodes.colspec('', colwidth=10)) + group.append(nodes.colspec('', colwidth=90)) + body = nodes.tbody('') + group.append(body) + + def append_row(*column_texts): + row = nodes.row('') + for text in column_texts: + node = nodes.paragraph('') + vl = ViewList() + vl.append(text, '') + self.state.nested_parse(vl, 0, node) + try: + if isinstance(node[0], nodes.paragraph): + node = node[0] + except IndexError: + pass + row.append(nodes.entry('', node)) + body.append(row) + return table, table_spec, append_row + + def get_table(self, items): + """Generate a proper list of table nodes for autosummary:: directive. + + *items* is a list produced by :meth:`get_items`. + """ + table, table_spec, append_row = self.get_tablespec() + for name, sig, summary, real_name in items: + # required for cpp autolink + # original + #qualifier = 'cpp:any' + #full_name = real_name.replace('.', '::') + kind = self.options['kind'] + qualifier = 'f:' + kind + # modified + #col1 = ':%s:`%s <%s>`' % (qualifier, name, full_name) + col1 = ':%s:`%s`' % (qualifier, name) + col2 = summary + append_row(col1, col2) + + # RemovedInSphinx40Warning: Autosummary.result is deprecated + # debug: find alternative + #self.result.append(' .. rubric: sdsf', 0) + # removed + #self.bridge.result.append(' .. rubric: sdsf', 0) + # debug + return [table_spec, table] + + +class DoxygenAutoEnum(DoxygenAutosummary): + + def get_items(self, names): + env = self.state.document.settings.env + self.name = names[0] + + real_name, obj, parent, modname = import_by_name(self.name, env=env) + names = [n.text for n in obj.findall('./enumvalue/name')] + descriptions = [format_xml_paragraph(d) for d in obj.findall('./enumvalue/detaileddescription')] + return zip(names, descriptions) + + def get_table(self, items): + table, table_spec, append_row = self.get_tablespec() + for name, description in items: + col1 = ':strong:`' + name + '`' + while description and not description[0].strip(): + description.pop(0) + col2 = ' '.join(description) + append_row(col1, col2) + return [nodes.rubric('', 'Enum: %s' % self.name), table] diff --git a/docs/_ext/autodoc_doxygen/autosummary/generate.py b/docs/_ext/autodoc_doxygen/autosummary/generate.py new file mode 100644 index 0000000000..41d0453c03 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autosummary/generate.py @@ -0,0 +1,391 @@ +import codecs +import os +import re +import sys + +from jinja2 import FileSystemLoader +from jinja2.sandbox import SandboxedEnvironment +from sphinx.jinja2glue import BuiltinTemplateLoader +from sphinx.util.osutil import ensuredir + +from . import import_by_name, get_doxygen_root +from ..xmlutils import format_xml_paragraph + + +def is_type(node): + def_node = get_doxygen_root().find('./compounddef[@id="%s"]' % node.get('refid')) + return def_node.get('kind') == 'type' + +def generate_autosummary_docs(sources, output_dir=None, suffix='.rst', + #base_path=None, builder=None, template_dir=None): + # add toctree argument + base_path=None, builder=None, template_dir=None, toctree=None, + build_mode=None): + + showed_sources = list(sorted(sources)) + if len(showed_sources) > 20: + showed_sources = showed_sources[:10] + ['...'] + showed_sources[-10:] + print('[autosummary] generating autosummary for: %s' % + ', '.join(showed_sources)) + + if output_dir: + print('[autosummary] writing to %s' % output_dir) + + if base_path is not None: + sources = [os.path.join(base_path, filename) for filename in sources] + + # create our own templating environment + template_dirs = [os.path.join(os.path.dirname(__file__), 'templates')] + + if builder is not None: + # allow the user to override the templates + template_loader = BuiltinTemplateLoader() + template_loader.init(builder, dirs=template_dirs) + else: + if template_dir: + template_dirs.insert(0, template_dir) + template_loader = FileSystemLoader(template_dirs) + #template_env = SandboxedEnvironment(loader=template_loader) + # modified + template_env = SandboxedEnvironment(loader=template_loader, + trim_blocks=True, lstrip_blocks=True) + + # read + items = find_autosummary_in_files(sources) + + # keep track of new files + new_files = [] + + for name, path, template_name in sorted(set(items), key=str): + # replace + path = path or output_dir or os.path.abspath(toctree) + # debug + + # this is extra? + #if path is None: + # # The corresponding autosummary:: directive did not have + # # a :toctree: option + # print("[debug] directive did not have a :toctree: option") + # continue + + #path = output_dir or os.path.abspath(path) + if builder.app.verbosity > 0: + print("[debug] checking path: %s" % (path)) + ensuredir(path) + + try: + name, obj, parent, mod_name = import_by_name(name) + except ImportError as e: + print('WARNING [autosummary] failed to import %r: %s' % (name, e), file=sys.stderr) + continue + + fn = os.path.join(path, name + suffix).replace('::', '.') + + # skip it if it exists + if os.path.isfile(fn): + continue + + # removed? + #new_files.append(fn) + + if template_name is None: + if obj.tag == 'compounddef' and obj.get('kind') in ['namespace', 'module']: + template_name = 'doxymodule.rst' + elif obj.tag == 'compounddef' and obj.get('kind') == 'page': + template_name = 'doxypage.rst' + else: + raise NotImplementedError('No template for %s (%s %s)' % (obj.items(), obj.tag, obj.get('kind'))) + + if builder.app.verbosity > 0: + print("[debug] template:%s kind: %s obj.items():%s" % (template_name, obj.get('kind'), obj.items())) + with open(fn, 'w') as f: + template = template_env.get_template(template_name) + # The ns keys feed into the template + ns = {} + if obj.tag == 'compounddef' and obj.get('kind') == 'namespace': + ns['methods'] = [e.text for e in obj.findall('./sectiondef[@kind="func"]/memberdef[@kind="function"]/name')] + ns['types'] = [e.text for e in obj.findall('./innerclass') if is_type(e)] + ns['objtype'] = 'namespace' + elif obj.tag == 'compounddef' and obj.get('kind') == 'page': + if builder.app.verbosity > 0: + print("[debug] xml parsing for %s" % (obj.get('id'))) + ns['title'] = obj.find('title').text + ns['underline'] = len(ns['title']) * '=' + #ns['text'] = format_xml_paragraph(obj.find('detaileddescription'),build_mode) + ns = format_xml_paragraph(obj.find('detaileddescription'), build_mode, nsOrig=ns, verbosity=builder.app.verbosity) + #if obj.get('id') == 'Specifics': + else: + raise NotImplementedError(obj) + + parts = name.split('::') + mod_name, obj_name = '::'.join(parts[:-1]), parts[-1] + + ns['fullname'] = name + ns['module'] = mod_name + ns['objname'] = obj_name + ns['name'] = parts[-1] + if not('underline' in ns): + ns['underline'] = len(name) * '=' + + rendered = template.render(**ns) + f.write(rendered) + # debug: date/time caching hack + # f.write('\n..\n {}'.format(datetime.datetime.now())) + + # descend recursively to new files + if new_files: + generate_autosummary_docs(new_files, output_dir=output_dir, + suffix=suffix, base_path=base_path, builder=builder, + #template_dir=template_dir) + # add toctree argument + template_dir=template_dir, toctree=toctree) + + +def find_autosummary_in_files(filenames): + """Find out what items are documented in source/*.rst. + + See `find_autosummary_in_lines`. + """ + # todo: break when this doesn't exist + # look for modules and standalone documentation pages, but *not* the index page + # itself (which it links to from itself for some reason...) + documented = [] + for filename in filenames: + with codecs.open(filename, 'r', encoding='utf-8', errors='ignore') as f: + lines = f.read().splitlines() + documented.extend(find_autosummary_in_lines(lines, filename=filename)) + + return documented + + +def find_autosummary_in_lines(lines, module=None, filename=None): + """Find out what items appear in autosummary:: directives in the + given lines. + + Returns a list of (name, toctree, template) where *name* is a name + of an object and *toctree* the :toctree: path of the corresponding + autosummary directive (relative to the root of the file name), and + *template* the value of the :template: option. *toctree* and + *template* ``None`` if the directive does not have the + corresponding options set. + """ + + # add generate_arg_re + autosummary_re = re.compile(r'^(\s*)\.\.\s+autodoxysummary::\s*') + toctree_arg_re = re.compile(r'^\s+:toctree:\s*(.*?)\s*$') + template_arg_re = re.compile(r'^\s+:template:\s*(.*?)\s*$') + kind_arg_re = re.compile(r'^\s+:kind:\s*(.*?)\s*$') + generate_arg_re = re.compile(r'^\s+:generate:\s*$') + autosummary_item_re = re.compile(r'^\s+(~?[_a-zA-Z][a-zA-Z0-9_.:]*)\s*.*?') + + documented = [] + + toctree = None + template = None + in_autosummary = False + generate = False + base_indent = "" + + for line in lines: + if in_autosummary: + m = toctree_arg_re.match(line) + if m: + toctree = m.group(1) + if filename: + toctree = os.path.join(os.path.dirname(filename), + toctree) + continue + + m = template_arg_re.match(line) + if m: + template = m.group(1).strip() + continue + + # add + + m = generate_arg_re.match(line) + if m: + generate = True + continue + + m = kind_arg_re.match(line) + if m and generate: + kind = m.group(1).strip() + xpath = None + if kind == 'mod': + xpath = './compound[@kind="namespace"]' + elif kind == 'page': + xpath = './compound[@kind="page" and not(@refid="indexpage")]' + + if xpath is not None: + results = get_doxygen_root().xpath(xpath) + for result in results: + documented.append((result.find('name').text, toctree, template)) + + continue + + # end add + + if line.strip().startswith(':'): + continue # skip options + + m = autosummary_item_re.match(line) + if m: + name = m.group(1).strip() + if name.startswith('~'): + name = name[1:] + documented.append((name, toctree, template)) + continue + + if not line.strip() or line.startswith(base_indent + " "): + continue + + in_autosummary = False + + m = autosummary_re.match(line) + if m: + in_autosummary = True + base_indent = m.group(1) + toctree = None + template = None + # add + generate = False + continue + + return documented + + +def _generate_source_stubs(app): + """Generate one :orphan: stub per doxygen file compound under + api/generated/source/, each invoking ``.. autodoxysource::``.""" + root = get_doxygen_root() + source_dir = os.path.join(app.srcdir, 'api', 'generated', 'source') + ensuredir(source_dir) + + template_dirs = [os.path.join(os.path.dirname(__file__), 'templates')] + template_loader = FileSystemLoader(template_dirs) + template_env = SandboxedEnvironment(loader=template_loader, + trim_blocks=True, lstrip_blocks=True) + template = template_env.get_template('doxysource.rst') + + files = root.findall('./compounddef[@kind="file"]') + count = 0 + for cd in files: + file_id = cd.get('id') + if file_id is None: + continue + # Only generate if there is a programlisting + if cd.find('.//programlisting') is None: + continue + + fn = os.path.join(source_dir, file_id + '.rst') + if os.path.isfile(fn): + continue + + # Title from the location filename, or fall back to file_id + loc = cd.find('location') + if loc is not None and loc.get('file'): + title = os.path.basename(loc.get('file')) + else: + title = file_id + + rendered = template.render( + title=title, + underline='=' * len(title), + file_id=file_id, + ) + with open(fn, 'w') as f: + f.write(rendered) + count += 1 + + if count: + print('[autodoxysource] generated %d source stubs in %s' % + (count, source_dir)) + + +def _generate_function_index(app): + """Generate api/functions.rst listing every function/subroutine + across all namespace compounds, with cross-reference links to + the function's entry on its module page.""" + root = get_doxygen_root() + fn = os.path.join(app.srcdir, 'api', 'functions.rst') + if os.path.isfile(fn): + return + + entries = [] + for cd in root.findall('./compounddef[@kind="namespace"]'): + modname = cd.find('compoundname') + if modname is None or modname.text is None: + continue + mod = modname.text + for md in cd.findall('.//sectiondef[@kind="func"]/memberdef[@kind="function"]'): + name_el = md.find('name') + if name_el is None or name_el.text is None: + continue + name = name_el.text + brief_el = md.find('briefdescription/para') + brief = '' + if brief_el is not None and brief_el.text: + brief = brief_el.text.strip().replace('|', r'\|') + qualified = '%s/%s' % (mod, name) + entries.append((name, mod, qualified, brief)) + + entries.sort(key=lambda e: e[0].lower()) + + lines = [ + '.. _Functions:', + '', + '=========', + 'Functions', + '=========', + '', + '.. list-table::', + ' :widths: 30 30 40', + ' :header-rows: 1', + '', + ' * - Name', + ' - Module', + ' - Description', + ] + for name, mod, qualified, brief in entries: + lines.append(' * - :f:func:`%s <%s>`' % (name, qualified)) + lines.append(' - :f:mod:`%s`' % mod) + lines.append(' - %s' % brief) + + lines.append('') + + with open(fn, 'w') as f: + f.write('\n'.join(lines)) + print('[autodoxysource] generated function index with %d entries at %s' % + (len(entries), fn)) + + +def process_generate_options(app): + genfiles = app.config.autosummary_generate + # add + toctree = app.config.autosummary_toctree + # This is important to handle \htmlonly and \latexonly directives + sphinx_build_mode = app.config.sphinx_build_mode + + if genfiles and not hasattr(genfiles, '__len__'): + env = app.builder.env + genfiles = [os.fspath(env.doc2path(x, base=None)) for x in env.found_docs + if os.path.isfile(env.doc2path(x))] + + if not genfiles: + return + + ext = list(app.config.source_suffix)[0] + genfiles = [genfile + (not genfile.endswith(ext) and ext or '') + for genfile in genfiles] + + generate_autosummary_docs(genfiles, builder=app.builder, + # add toctree argument + # suffix=ext, base_path=app.srcdir) + suffix=ext, base_path=app.srcdir, toctree=toctree, build_mode=sphinx_build_mode) + + # Generate source browser stubs + _generate_source_stubs(app) + + # Generate function index + _generate_function_index(app) diff --git a/docs/_ext/autodoc_doxygen/autosummary/templates/doxymodule.rst b/docs/_ext/autodoc_doxygen/autosummary/templates/doxymodule.rst new file mode 100644 index 0000000000..2c89736453 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autosummary/templates/doxymodule.rst @@ -0,0 +1,34 @@ +.. autodoxymodule:: {{ fullname }} + :members: + {% if methods %} + :methods: + {% endif %} + {% if types %} + :types: + {% endif %} + + {% if types %} + ---------- + Data Types + ---------- + + .. autodoxysummary:: + :kind: type + + {% for item in types %} + ~{{ item }} + {% endfor %} + {% endif %} + + {% if methods %} + --------------------- + Functions/Subroutines + --------------------- + + .. autodoxysummary:: + :kind: func + + {% for item in methods %} + ~{{ fullname }}::{{ item }} + {% endfor %} + {% endif %} diff --git a/docs/_ext/autodoc_doxygen/autosummary/templates/doxypage.rst b/docs/_ext/autodoc_doxygen/autosummary/templates/doxypage.rst new file mode 100644 index 0000000000..96320dc499 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autosummary/templates/doxypage.rst @@ -0,0 +1,21 @@ +.. _{{ name }}: +{# comment +When the name is provided, we get "(INFO/1) Duplicate implicit target name:" +without this, we get undefined reference. This needs to be fixed later. +#} + +{{ underline }} +{{ title }} +{{ underline }} + +{% for line in text %} +{{ line }} +{% endfor %} +{% if footnotes %} + +.. rubric:: Footnotes + +{% for line in footnotes %} +.. [#] {{ line }} +{% endfor %} +{% endif %} diff --git a/docs/_ext/autodoc_doxygen/autosummary/templates/doxysource.rst b/docs/_ext/autodoc_doxygen/autosummary/templates/doxysource.rst new file mode 100644 index 0000000000..cef1c4cab7 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/autosummary/templates/doxysource.rst @@ -0,0 +1,6 @@ +:orphan: + +{{ title }} +{{ underline }} + +.. autodoxysource:: {{ file_id }} diff --git a/docs/_ext/autodoc_doxygen/xmlutils.py b/docs/_ext/autodoc_doxygen/xmlutils.py new file mode 100644 index 0000000000..98459f5da9 --- /dev/null +++ b/docs/_ext/autodoc_doxygen/xmlutils.py @@ -0,0 +1,997 @@ +import re + +from . import get_doxygen_root, get_doxygen_id_index + + +def flatten(xmlnode): + # this.textchild0.textchild0.tail... + + t = '' + + # text of this node + if xmlnode.text is not None: + t += xmlnode.text + + # process all children recursively + for n in xmlnode: + t += ' ' + t += flatten(n) + if n.tail is not None: + t += ' ' + t += n.tail + + return t + +def format_xml_paragraph(xmlnode,build_mode,nsOrig=None,verbosity=0): + """Format an Doxygen XML segment (principally a detaileddescription) + as a paragraph for inclusion in the rst document + + Parameters + ---------- + xmlnode + + Returns + ------- + lines + A list of lines. + """ + # Here we are operating on the entire document for the template + # This helps support \footnotes{} + if nsOrig is not None: + xmlParagraphFormatter = _DoxygenXmlParagraphFormatter() + xmlParagraphFormatter.setNS(nsOrig) + xmlParagraphFormatter.setVerbosity(verbosity) + xmlParagraphFormatter.generic_visit(xmlnode,build_mode=build_mode) + xmlParagraphFormatter.ns['text'] = [l.rstrip() for l in xmlParagraphFormatter.lines] + return xmlParagraphFormatter.ns + else: + # Return processing for typically ns['text'] only + # Expand to allow setting of options + xmlParagraphFormatter = _DoxygenXmlParagraphFormatter() + xmlParagraphFormatter.setVerbosity(verbosity) + xmlParagraphFormatter.generic_visit(xmlnode,build_mode=build_mode) + return [l.rstrip() for l in xmlParagraphFormatter.lines] + +class _DoxygenXmlParagraphFormatter(object): + # This class follows the model of the stdlib's ast.NodeVisitor for tree traversal + # where you dispatch on the element type to a different method for each node + # during the traverse. + + # It's supposed to handle paragraphs, references, preformatted text (code blocks), and lists. + + def __init__(self): + self.ns = {} + self.lines = [''] + self.continue_line = False + # We need to track specified math labels and place them prior to the ".. math::" blocks + self.math_labels = [] + self.build_mode = None + self.verbosity = 0 + self.indent = -1 + self.options = [] + + # new + def setNS(self, ns): + self.ns = ns + + def setVerbosity(self, verbosity): + self.verbosity = verbosity + if self.verbosity > 0: print("[debug] verbosity = %s" % (self.verbosity)) + + def visit_latexonly(self, node): + if not(self.build_mode in ('latexpdf','latex')): + return + + text = node.text + if text == None: + return + + # Convert \\ref{tag} to :ref:` ` and the sphinx latex processor + # converts it to a proper label reference. + ref_match = re.search('\\\\ref{(.*?)}', text) + if ref_match is not None: + tag_string = ref_match.groups()[0] + #val = [' :ref:`%s`' % tag_string] + val = [':latex:`\\ref{%s}`' % tag_string] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + # If we have then skip DoxyImage provided material + if 'skipDoxyImage' in self.options: + if text.find('DoxyImage') >= 0: + return + + # At this point, just pass everything through to latex + self.concat_text(':latex:`%s`' % (text)) + + return + + # new + # Newer versions of doxygen utilize tag in XML + # Doxygen 1.8.13 leaves all this in see: para_eqref + def visit_htmlonly(self, node): + if self.build_mode != 'html': + return + + text = node.text + if text == None: + return + + # Check for \eqref2{tag,txt} and convert to :ref:`tag`_ + eqref_match = re.search('\\\\eqref2{(.*?)}', text) + if eqref_match is not None: + tag_string = eqref_match.groups()[0] + if tag_string.find(',') >= 0: + fc = tag_string.find(',') + val = [':math:numref:`%s` - %s' % (tag_string[0:fc],tag_string[fc+1:])] + else: + val = [':math:numref:`%s`' % tag_string] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + # This supports \footnotes{} + if text.find('title=') >= 0: + text = text.replace('\n',' ') + title_match = re.search('title="(.*)"', text) + if title_match: + title_string = title_match.groups()[0] + # Recover \cite that have been converted to @cite to :cite:`%s` + if title_string.find('@cite') >= 0: + citeCommand = '@cite ([\w\-\_]+)' + m = re.search(citeCommand, title_string) + while m: + replStr = title_string[m.start():m.end()] + newStr = ':cite:`%s`' % (m.groups()[0]) + title_string = title_string.replace(replStr, newStr) + m = re.search(citeCommand, title_string) + if 'footnotes' in self.ns: + self.ns['footnotes'].append(title_string) + else: + self.ns['footnotes'] = [title_string] + + val = ["[#]_"] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + # Check for \eqref{ replace with :ref:`tag`_ + # Post processing of equations will place a link into the HTML + eqref_match = re.search('\\\\eqref{(.*?)}', text) + if eqref_match is not None: + tag_string = eqref_match.groups()[0] + val = [':math:numref:`%s`' % tag_string] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + # undefined + if self.verbosity > 0: + print("[debug] WARNING: Uncaptured htmlonly string (%s)" % text) + + # new + # reStructured text only permits one label per math:: block + def emit_math_labels(self): + if len(self.math_labels) == 0: + return + + if self.verbosity > 0: print("[debug] inserting math labels") + + math_block_idx = -1 + for idx in range(len(self.lines)-1,0,-1): + if self.lines[idx].startswith('.. math::'): + math_block_idx = idx + break + + # Add new label right after the math:: block + if math_block_idx >=0: + new_lines = self.lines[0:math_block_idx+1] + new_label = " :label: %s" % (self.math_labels[0]) + new_lines.append(new_label) + #new_lines.append('') + self.blank_line() + new_lines = new_lines + self.lines[math_block_idx+1:] + self.lines = new_lines + + self.math_labels = [] + + # Add appropriate implicit labels from anchors + def visit_anchor(self, node): + if self.verbosity > 0: + print("[debug] anchor id(%s)" % (node.get('id'))) + citeID = node.get('id') + if citeID.find('_1CITE') == 0: + citeID = "citeref_%s" % (citeID) + implicitLink = '.. _%s:' % (citeID) + self.lines.append(implicitLink) + #self.lines.append('') + self.blank_line() + + # Original + def visit(self, node): + method = 'visit_' + node.tag + if self.verbosity > 0: print("[debug] method=%s" % (method)) + if len(self.math_labels) > 0 and node.tag != 'formula': + self.emit_math_labels() + visitor = getattr(self, method, self.generic_visit) + return visitor(node) + + def generic_visit(self, node, build_mode=None): + if build_mode: + if self.verbosity > 2: print("[debug] Setting build mode: %s" % (build_mode)) + self.build_mode = build_mode + # Perform a scan for htmlonly or latexonly to prevent double processing of + # references + if not('scanned' in self.options): + self.options.append('scanned') + self.scanNode(node) + for child in node.getchildren(): + self.visit(child) + return self + + # Scan the node and set appropriate options + def scanNode(self, node): + # NOTE: these XPath expressions must use './/' rather than '//'. + # In XPath, '//foo' is an abbreviation for /descendant-or-self::node()/foo + # starting from the *document root*, not from `node`. Because our + # autodoc_doxygen extension concatenates every doxygen XML file into a + # single merged tree (see set_doxygen_xml), every `node` passed here is + # a small subtree (e.g. a single ) whose owner + # document is the *entire* MOM6 doxygen output. Using '//' here + # therefore scans the whole merged tree on every call, which made + # scanNode the dominant cost of `make html` -- 75% of single-threaded + # build time at full MOM6 input scale, quadratic in the tree size. + # Using './/' scans only descendants of the actual node, which is what + # was intended and makes each call O(local subtree size). + xp = node.xpath('.//latexonly') + if len(xp) > 0: + self.options.append('latexonly') + xp = node.xpath('.//htmlonly') + if len(xp) > 0: + self.options.append('htmlonly') + + if 'latexonly' in self.options: + xp = node.xpath('.//image[@type="latex"]') + if len(xp) > 0: + self.options.append('skipDoxyImage') + + def visit_ref(self, node): + refid = node.get('refid') + name_node = None + ream_name = None + kind = None + + # O(1) lookup via the lazily-built id -> element index. + # The previous findall('.//*[@id=X]') on the merged tree was the + # single largest cost in `make html` under XML_PROGRAMLISTING=YES + # -- see docs/REMAINING_TASKS.md / profile notes. + hit = get_doxygen_id_index().get(refid) + if self.verbosity > 0: print("[debug] refid(%s) kindref(%s) ref(%s)" % + (refid, node.get('kindref'), hit)) + if hit is not None: + ref = hit + kind = ref.get('kind') + if self.verbosity > 0: print("[debug] ref(%s)" % ref.items()) + if ref.tag == 'memberdef': + parent = ref.xpath('./ancestor::compounddef/compoundname')[0].text + name = ref.find('./name').text + real_name = parent + '::' + name + elif ref.tag in ('compounddef', 'enumvalue'): + if kind == 'page': + # :ref: works, but requires an explicit tag placed at the top of pages + # that generates an INFO message. FIX LATER. + val = [':ref:`%s`' % ref.get('id')] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + name_node = ref.find('./name') + real_name = name_node.text if name_node is not None else '' + elif ref.tag in ('anchor','sect1','sect2','sect3','sect4'): + # If _1CITEREF_ this is a doxygen processed citation + if refid.find('_1CITEREF_') >= 0: + citation = refid[18:] + val = [':cite:`%s`' % (citation)] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + # Capture sectional links + + # Treat the rest of these as general links + if refid.find('_1') >= 0: + reftext = node.text + reftext = reftext.strip() + refid2 = refid[refid.find('_1')+2:] + if reftext != '' and reftext != refid2: + if self.verbosity > 0: print("[debug] refid2(%s) reftext(%s)" % (refid2,reftext)) + val = [':ref:`%s<%s>`' % (reftext,refid)] + else: + if self.verbosity > 0: print("[debug] refid(%s)" % (refid)) + val = [':ref:`%s`' % refid] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + else: + print('[error] Unimplemented anchor tag: %s' % (ref.tag)) + raise NotImplementedError(ref.tag) + else: + print('[error] Unimplemented tag: %s' % (ref.tag)) + raise NotImplementedError(ref.tag) + else: + real_name = None + + + # Older doxygen support 1.8.13 for citation references + if node.get('kindref') == 'member' and refid.find('_1CITEREF_') >= 0: + citation = refid[18:] + val = [':cite:`%s`' % (citation)] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + # if kind='file' treat as file references + if kind == 'file': + # for now treat these as text + # TODO: references to code + val = ['``%s``' % node.text] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + return + + #debug + code_type = 'f' + if code_type == 'f': + val = [':%s:func:`%s' % (code_type, node.text)] + else: + val = [':%s:any:`' % code_type, node.text] + if real_name: + val.extend((' <%s>`' % (real_name))) + else: + val.append('`') + if node.tail is not None: + val.append(node.tail) + + if self.verbosity > 0: print("[debug] kind(%s) real_name(%s) node_name(%s)" % + (kind, real_name, name_node)) + #self.lines[-1] += ''.join(val) + self.concat_text(''.join(val)) + + # add visit_ulink + def visit_ulink(self, node): + self.para_text('`%s <%s>`_' % (node.text, node.get('url'))) + + # add visit_emphasis + def visit_emphasis(self, node): + self.para_text('*%s*' % node.text) + + # add role_text + def role_text(self, node, role): + # Is this even used? + if self.verbosity > 0: + print("[debug] role_text") + # XXX we should probably escape preceeding whitespace... + # but there's no backward equivalent of `tail` + text = ' :%s:`%s`' % (role, node.text) + + if node.tail is not None and not node.tail.startswith(' '): + # escape following whitespace + text += '\\' + + text += ' ' # interpretered text needs surrounding whitespace + self.para_text(text) + + # add visit_image + def visit_image(self, node): + + # Filter activity based on build type and type of image + image_type = node.get('type') + if image_type == 'html' and self.build_mode != 'html': + return + if image_type == 'latex' and not(self.build_mode in ('latexpdf','latex')): + return + + if self.verbosity > 0: + print("[debug] image type(%s) mode(%s)" % (image_type, self.build_mode)) + + # node.text is None for an empty element (no caption text); + # treat that the same as an empty caption and emit `.. image::` rather + # than `.. figure::`. The fork's original code crashed with + # AttributeError on these. Doxygen produces empty elements for + # cases like an image referenced from a `\image` command with no caption. + if node.text and node.text.strip(): + type = 'figure' + else: + type = 'image' + + self.lines.append('.. %s:: /images/%s' % (type, node.get('name'))) + + if type in 'figure': + # NOTE: Escaped math equations do not play nicely with "literal strings" in python! + + caption = node.text + + # Detect simple math commands and replace them with sphinx :math: directives + mathCommand = '\\\\f\$(.*?)\\\\f\$' + m = re.search(mathCommand, caption) + while m: + replStr = caption[m.start():m.end()] + newStr = ':math:`%s`' % (m.groups()[0]) + caption = caption.replace(replStr, newStr) + m = re.search(mathCommand, caption) + + # Only html needs to be double escaped + if image_type == 'html': + caption = node.text.replace('\\','\\\\') + + #if caption.find('Phi') >= 0: + + # Search for $[math]$ and convert to \([math]\) for html + # Use :math: for latex + mathCommand = '\$(.*?)\$' + m = re.search(mathCommand, caption) + while m: + replStr = caption[m.start():m.end()] + if image_type == 'html': + newStr = '\\\\(%s\\\\)' % (m.groups()[0]) + else: + newStr = ':math:`%s`' % (m.groups()[0]) + caption = caption.replace(replStr, newStr) + m = re.search(mathCommand, caption) + + # For html, scan for \\f and remove that too + if image_type == 'html': + caption = caption.replace('\\f','') + + if self.verbosity > 0: + # For math we have to double the number of escapes so we pass an + # escape from RST to HTML. + print("[debug] caption text(%s)" % (caption)) + self.lines.extend(['', " %s" % (caption), '']) + + # add visit_superscript + def visit_superscript(self, node): + self.role_text(node, 'superscript') + + # add visit_subscript + def visit_subscript(self, node): + self.role_text(node, 'subscript') + + # add visit_sup + # Support for doxygen 1.8.13 as it passes everything to + # Support for doxygen \footnote{} + def visit_sup(self, node): + + # Skip if we detect htmlonly or latexonly + if self.para_ignore(): + return + + title_string = node.get('title') + if title_string: + citeCommand = '@cite ([\w\-\_]+)' + m = re.search(citeCommand, title_string) + while m: + replStr = title_string[m.start():m.end()] + newStr = ':cite:`%s`' % (m.groups()[0]) + title_string = title_string.replace(replStr, newStr) + m = re.search(citeCommand, title_string) + + if 'footnotes' in self.ns: + self.ns['footnotes'].append(title_string) + else: + self.ns['footnotes'] = [title_string] + + val = ["[#]_"] + #self.lines[-1] += ''.join(val) + self.concat_text(val[0]) + + # Ignore duplicates provided by xmlonly if we detect latexonly or htmlonly + def para_ignore(self): + + if 'latexonly' in self.options or 'htmlonly' in self.options: + return True + return False + + # add replace any references of \eqref, \eqref2, \eqref4 + # Doxygen 1.8.13 + # html: use eqref2; remove eqref4 + # latex: use eqref4; remove eqref2 + # with appropriate replacements + # Remove duplicates here if latexonly or htmlonly is detected + def para_eqref(self, text): + + chg = True + while text.find('\\\\eqref2') >= 0 and chg: + chg = False + m = re.search('\\\\eqref2{(.*?)}', text) + if m: + ref = m.groups()[0] + fullRef = '\\\\eqref2{%s}' % (ref) + if ref.find(',') >= 0: + i = ref.find(',') + sphinxRef = ':math:numref:`%s` - %s' % (ref[0:i],ref[i+1:]) + if self.build_mode in ('latexpdf','latex'): + sphinxRef = '' + if self.para_ignore(): + sphinxRef = '' + text = text.replace(fullRef, sphinxRef) + chg = True + + chg = True + while text.find('\\\\eqref') >= 0 and chg: + chg = False + m = re.search('\\\\eqref{(.*?)}', text) + if m: + ref = m.groups()[0] + fullRef = '\\\\eqref{%s}' % (ref) + if self.build_mode in ('latexpdf','latex'): + sphinxRef = ':latex:`\\ref{%s}`' % ref + else: + sphinxRef = ':math:numref:`%s`' % (ref) + if self.para_ignore(): + sphinxRef = '' + text = text.replace(fullRef, sphinxRef) + chg = True + + chg = True + while text.find('\\\\eqref4') >= 0 and chg: + chg = False + m = re.search('\\\\eqref4{(.*?)}', text) + if m: + ref = m.groups()[0] + fullRef = '\\\\eqref4{%s}' % (ref) + sphinxRef = ':latex:`\\ref{%s}`' % (ref) + if self.build_mode in ('html'): + sphinxRef = '' + if self.para_ignore(): + sphinxRef = '' + text = text.replace(fullRef, sphinxRef) + chg = True + + return text + + # Assistant for ensuring there is blank lines between directives + # It makes sure we do not overly add blank lines + def blank_line(self): + if len(self.lines) == 0: + return + + if self.lines[-1] == '': + return + + self.lines.append('') + + # Assistant for putting sentences together + def concat_text(self, text): + if len(self.lines) == 0: + self.lines.append(text) + return + + lastLine = self.lines[-1] + + if len(lastLine) == 0: + self.lines[-1] = text + return + + lastChar = lastLine[-1] + newText = text + if len(newText) == 0: + return + + firstChar = newText[0] + + # Emphasis + if lastChar == "*" or firstChar == "*": + newText = " %s" % (newText) + firstChar = " " + + # whitespace after :cite:`tag` + if lastChar == '`': + if (firstChar >= 'a' and firstChar <= 'z') or (firstChar >= 'A' and firstChar <= 'Z') or firstChar in ['(','[','{']: + newText = " %s" % (newText) + firstChar = " " + + # whitespace before :commands: + if firstChar == ':': + if (lastChar >= 'a' and lastChar <= 'z') or (lastChar >= 'A' and lastChar <= 'Z') or lastChar in [',','.','=']: + newText = " %s" % (newText) + firstChar = " " + + # Footnotes and any items that end with _ + if newText == '[#]_': + newText = " %s" % (newText) + if lastChar == '_': + if len(lastLine) > 3: + if lastLine[-4:] == "[#]_" and firstChar != '.': + newText = " %s" % (newText) + firstChar = " " + else: + newText = " %s" % (newText) + + # Inline text check for space before (``) + if len(newText) >= 2: + if newText[0:2] == "``" and lastChar != ' ': + newText = " %s" % (newText) + firstChar = " " + + self.lines[-1] += newText + return + + # add para_text parser + # Doxygen 1.8.13 support for \eqref \eqref2 + def para_text(self, text): + + if text is not None: + if text.find('Some time later') >= 0: + a = 0 + + if text.find('eqref') >= 0: + text = self.para_eqref(text) + if self.continue_line: + if len(self.lines) >= 1: + # If we are in a continue_line situation but already + # have a linefeed, do an append instead + if self.lines[-1] == '': + self.lines.append(text) + return + self.concat_text(text) + else: + self.lines.append(text) + + def visit_para(self, node): + + self.para_text(node.text) + + # visit children and append tail + for child in node.getchildren(): + self.visit(child) + self.continue_line = True + + if child.tail is not None: + self.para_text(child.tail.lstrip()) + + # replaced + #if node.text is not None: + # if self.continue_line: + # self.lines[-1] += node.text + # else: + # self.lines.append(node.text) + #self.generic_visit(node) + + self.continue_line = False + #self.lines.append('') + self.blank_line() + + # add visit_formula + def visit_formula(self, node): + text = node.text + + # Remove the faked link for pdf version + if self.build_mode in ('latexpdf','latex'): + label_match = re.search(' \\\\label{(html:.*?)}.*?\\\\\\\\', text) + if label_match: + replace_string = label_match.group() + text = text.replace(replace_string,'') + + # detect inline or block math + if text.startswith('\\[') or not text.startswith('$'): + if text.startswith('\\['): + text = text[2:-2] + + # if we are emitting a math block and we have + # pending math labels, go back and emit those + # first. + if len(self.math_labels) > 0: + self.emit_math_labels() + + self.blank_line() + if '\n' in text: + self.lines.append('.. math::') + self.lines.append('') + for mathline in text.split('\n'): + self.lines.append(' ' + mathline) + else: + self.lines.append('.. math:: ' + text) + self.blank_line() + # Math blocks require an explicit blank line as well? + #self.lines.append('') + self.continue_line = False + else: + inline = ':math:`' + node.text.strip()[1:-1].strip() + '`' + if self.continue_line: + #self.lines[-1] += inline + self.concat_text(inline) + else: + self.lines.append(inline) + + self.continue_line = True + + # detect \label{html:tag} blocks + if text.find('\\label') >= 0: + # If we have a big block of equations, supply one label + label_matches = re.findall('\\\label{html:(.*?)?}',text) + if len(label_matches) > 0: + [self.math_labels.append(i) for i in label_matches] + else: + label_matches = re.findall('\\\label{(.*?)?}',text) + if len(label_matches) > 0: + [self.math_labels.append(i) for i in label_matches] + if self.verbosity > 0: + # For math we have to double the number of escapes so we pass an + # escape from RST to HTML. + print("[debug] math_labels(%s)" % (label_matches)) + + def visit_parametername(self, node): + if 'direction' in node.attrib: + direction = '[%s] ' % node.get('direction') + else: + direction = '' + + param_name = node.text or '' + + # Look up the parameter's Fortran type from the parent + # memberdef's elements. The inside + # only carries names and descriptions; + # the types live on the siblings of . + param_type = '' + memberdefs = node.xpath('./ancestor::memberdef') + if memberdefs: + for p in memberdefs[0].findall('param'): + defname = p.find('defname') + if defname is not None and defname.text == param_name: + type_el = p.find('type') + if type_el is not None: + param_type = ''.join(type_el.itertext()).strip() + break + + # Prepend the type as literal text in the description rather + # than using :param type name: (sphinx-fortran's regex can't + # handle commas in the type) or :type name: (sphinx-fortran's + # xref resolver crashes on % in dimension expressions). + if param_type: + self.lines.append(':param %s: ``%s`` %s' % (param_name, param_type, direction)) + else: + self.lines.append(':param %s: %s' % (param_name, direction)) + self.continue_line = True + + def visit_parameterlist(self, node): + lines = [l for l in type(self)().generic_visit(node).lines if l != ''] + # replaced + #self.lines.extend([':parameters:', ''] + ['* %s' % l for l in lines] + ['']) + self.lines.extend([''] + lines + ['']) + + # TODO: Doxygen generates a simplesect for functions with + # a specified return argument. For now, we leave as + # :returns undefined: + # marker so we can fix up the document using flint. + # Supports doxygen /sa or /see command + def visit_simplesect(self, node): + if self.verbosity > 0: + print("[debug] simplesect kind(%s)" % (node.get('kind'))) + + # Do nothing for \note for now + + # fortran function handling + if node.get('kind') == 'return': + self.lines.append(':returns undefined: ') + self.continue_line = True + self.generic_visit(node) + + # Add bold text psudo section for \see, \sa roughly acts like doxygen + if node.get('kind') in ('see', 'sa'): + see_also_label = "See also" + #self.lines.append('') + self.blank_line() + self.lines.append('**%s**' % (see_also_label)) + #self.lines.append('') + self.blank_line() + #self.lines.append('') + self.generic_visit(node) + # add + + def visit_sect(self, node, char): + """Generic visit section""" + title_node = node.find('title') + if title_node is not None: + title = title_node.text + # Filter html data (possibly if we see a <, / and >) + if self.verbosity > 0: + print("[debug] visit_sect id(%s) title(%s)" % (node.get('id'),title)) + if title.find('<') >=0 and title.find('>') >=0 and title.find('/') >=0: + html_match = False + # Filter => `` + if title.find("") >= 0: + title = title.replace('','``') + title = title.replace('','`` ') + html_match = True + if not(html_match) and self.verbosity > 0: + print("[debug] unmatched html (%s)" % (title)) + # Add a implicit lable for the sections + implicitLink = '.. _%s:' % (node.get('id')) + self.lines.append(implicitLink) + #self.lines.append('') + self.blank_line() + self.lines.append(title) + self.lines.append(len(title) * char) + #self.lines.append('') + self.blank_line() + + self.generic_visit(node) + + def visit_sect1(self, node): + self.visit_sect(node, '=') + + def visit_sect2(self, node): + self.visit_sect(node, '-') + + def visit_sect3(self, node): + self.visit_sect(node, '^') + + def visit_sect4(self, node): + self.visit_sect(node, '"') + + # add end + + # allows us to handle nested ordered lists + def visit_orderedlist(self, node): + self.indent = self.indent + 1 + self.generic_visit(node) + #self.lines.append('') + self.blank_line() + self.indent = self.indent - 1 + + # allows us to handle nested itemized lists + def visit_itemizedlist(self, node): + self.indent = self.indent + 1 + self.generic_visit(node) + #self.lines.append('') + self.blank_line() + self.indent = self.indent - 1 + + # Source of citation and numbering + def visit_listitem(self, node): + #char = '*' if node.getparent().tag == 'itemizedlist' else '#.' + if node.getparent().tag == 'itemizedlist': + #self.lines.append('') + char = '*' + else: + char = '#.' + if self.verbosity > 1: print("[debug] listitem indent = %s" % (self.indent)) + self.lines.append(' '*(self.indent*2) + char + ' ') + # replaced + #self.lines.append(' - ') + self.continue_line = True + # recursion + self.generic_visit(node) + + # add + def preformat_text(self, lines): + self.lines.extend(('::', '')) + self.lines.extend([' ' + l for l in lines]) + #self.lines.append('') + self.blank_line() + + def visit_preformatted(self, node): + segment = [node.text if node.text is not None else ''] + for n in node.getchildren(): + segment.append(n.text) + if n.tail is not None: + segment.append(n.tail) + + lines = ''.join(segment).split('\n') + # add line + self.preformat_text(lines) + # extra? no effect + #self.lines.extend(('.. code-block:: C++', '')) + #self.lines.extend([' ' + l for l in lines]) + + # add + def visit_programlisting(self, node): + lines = [] + for n in node.getchildren(): + lines.append(flatten(n)) + self.preformat_text(lines) + + #add + def visit_verbatim(self, node): + self.visit_preformatted(node) + + def visit_computeroutput(self, node): + c = node.find('preformatted') + if c is not None: + return self.visit_preformatted(c) + # add + # I don't think we can put links inside + # computeroutput text... + #self.lines[-1] += '``' + flatten(node) + '`` ' + self.concat_text('``' + flatten(node) + '``') + # omitted + #return self.visit_preformatted(node) + + def visit_xrefsect(self, node): + if node.find('xreftitle').text == 'Deprecated': + sublines = type(self)().generic_visit(node).lines + self.lines.extend(['.. admonition:: Deprecated'] + [' ' + s for s in sublines]) + return + # add - if not depricated + title = node.find('xreftitle').text + sublines = type(self)().generic_visit(node).lines + self.lines.extend(['.. admonition:: %s' % title] + [' ' + s for s in sublines]) + #else: + # raise ValueError(node) + + def visit_subscript(self, node): + #self.lines[-1] += '\ :sub:`%s` %s' % (node.text, node.tail) + self.concat_text(':sub:`%s` %s' % (node.text, node.tail)) + + def visit_table(self, node): + # save the number of columns + cols = int(node.get('cols')) + table = [] + # save the current output + lines = self.lines + + # get width of each column + widths = [0] * cols + + # build up the table contents + for row_node in node.findall('row'): + row = [] + for i, entry in enumerate(row_node.getchildren()): + self.lines = [''] + self.generic_visit(entry) + row.append(self.lines) + + # find width of this entry (including leading and trailing space) + widths[i] = max(widths[i], max([len(line) for line in self.lines]) + 2) + + table.append(row) + + def append_row(row): + # find number of lines in row + num_lines = max([len(e) for e in row]) + lines = [] + + for k in range(num_lines): + line = '|' + for i, e in enumerate(row): + if k < len(e): + # this is a valid line + line += ' ' + e[k] + # pad rest of line + line += ' ' * (widths[i] - len(e[k]) - 1) + else: + # invalid line, just fill with spaces + line += ' ' * widths[i] + + line += '|' + + lines.append(line) + + return lines + + self.lines = lines + # start with a blank + #self.lines.append('') + self.blank_line() + + # usual separator line + sep = '+' + for width in widths: + sep += '-' * width + sep += '+' + + self.lines.append(sep) + + # header row + self.lines.extend(append_row(table[0])) + # header separator uses '=' instead of '-' + self.lines.append(sep.replace('-', '=')) + + # loop over body rows + for row in table[1:]: + self.lines.extend(append_row(row)) + self.lines.append(sep) + + # end with a blank + #self.lines.append('') + self.blank_line() diff --git a/docs/_static/autodoxysource.css b/docs/_static/autodoxysource.css new file mode 100644 index 0000000000..81f42425a9 --- /dev/null +++ b/docs/_static/autodoxysource.css @@ -0,0 +1,94 @@ +/* Source browser styling for autodoxysource directive. + Font stack and sizes match the sphinx_rtd_theme's code-block rules + so source listings look consistent with fenced code blocks. */ + +.autodoxysource { + font-family: SFMono-Regular, Menlo, Monaco, Consolas, "Liberation Mono", + "Courier New", Courier, monospace; + font-size: 12px; + line-height: 1.2; + background: #f8f8f8; + border: 1px solid #e1e4e5; + border-radius: 4px; + padding: 12px 0; + overflow-x: auto; + min-width: fit-content; + max-width: calc(6.5em + 120ch + 24px); +} + +.autodoxysource .source-line { + display: block; + white-space: nowrap; + margin: 0; + padding: 0 12px 0 0; + line-height: 1.2; +} + +.autodoxysource .source-code { + white-space: pre; +} + +.autodoxysource .source-line:target { + background-color: #ffffcc; +} + +.autodoxysource .source-lineno { + display: inline-block; + width: 4.5em; + text-align: right; + padding-right: 0.8em; + margin-right: 0.8em; + color: #858585; + border-right: 1px solid #e6e9ea; + text-decoration: none; + -webkit-user-select: none; + user-select: none; +} + +.autodoxysource .source-lineno:hover { + color: #333; +} + +/* Doxygen highlight classes */ +.f-hl-comment { + color: #408080; + font-style: italic; +} + +.f-hl-keyword { + color: #008000; + font-weight: bold; +} + +.f-hl-keywordtype { + color: #b00040; +} + +.f-hl-keywordflow { + color: #008000; + font-weight: bold; +} + +.f-hl-stringliteral { + color: #ba2121; +} + +.f-hl-preprocessor { + color: #bc7a00; +} + +.f-hl-normal { + color: #333; +} + +/* Make xref links within source subtle */ +.autodoxysource a.reference { + color: inherit; + text-decoration: none; + border-bottom: 1px dotted #999; +} + +.autodoxysource a.reference:hover { + border-bottom: 1px solid #333; + color: #2980b9; +} diff --git a/docs/api/files.rst b/docs/api/files.rst new file mode 100644 index 0000000000..b90fd69b86 --- /dev/null +++ b/docs/api/files.rst @@ -0,0 +1,11 @@ +.. _Files: + +============ +Source Files +============ + +.. toctree:: + :maxdepth: 1 + :glob: + + generated/source/* diff --git a/docs/apiref.rst b/docs/apiref.rst index 7f94e1c4b9..2e884ab1a6 100644 --- a/docs/apiref.rst +++ b/docs/apiref.rst @@ -11,3 +11,5 @@ The complete API documentation is generated with doxygen and can be found at htt :maxdepth: 1 api/modules + api/functions + api/files diff --git a/docs/conf.py b/docs/conf.py index 4407d88356..50bb5a9e52 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -20,11 +20,144 @@ # 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('.')) +sys.path.insert(0, os.path.abspath('_ext')) # -- Custom configuration values and roles ----------------------------------- from docutils import nodes +# -- Monkey-patch: fix sphinx-fortran's broken parallel-build merge ---------- +# +# Upstream VACUMM/sphinx-fortran (as of commit reachable from master, 2025-10) +# ships a FortranDomain.merge_domaindata() that has two bugs which together +# cause every f-domain object to be lost when sphinx-build is run with -j > 1: +# +# 1. The function references a name `outNames` that does not exist +# (typo for `ourNames`); accessing it raises NameError, which Sphinx's +# parallel worker error path swallows silently. +# +# 2. Even with the typo fixed, the unpack `for name, docname in +# otherdata['modules'].items()` is wrong, because `modules` values are +# 4-tuples `(docname, synopsis, platform, deprecated)` and `objects` +# values are 2-tuples `(docname, type)`, not bare docnames. +# +# Symptom: with `make html` (which passes -j 4), env.domaindata['f'] ends up +# with 0 modules and 0 objects, every :f:func:/:f:type:/:f:mod: cross- +# reference fails to resolve, and pages like f-modindex.html disappear. +# +# We patch this in-process so the existing -j 4 build still works. The fix +# is small and obvious; once it lands upstream we should pin sphinx-fortran +# to a post-fix commit and remove this patch. +# +# TODO(piece-2): submit upstream PR to VACUMM/sphinx-fortran with this fix, +# then drop this monkey-patch and pin requirements.txt to a +# post-fix commit. +def _patch_sphinx_fortran_merge_domaindata(): + try: + from sphinxfortran.fortran_domain import FortranDomain + except ImportError: + return + def merge_domaindata(self, docnames, otherdata): + ourNames = self.data['modules'] + for name, data in otherdata['modules'].items(): + # data == (docname, synopsis, platform, deprecated) + if data[0] in docnames and name not in ourNames: + ourNames[name] = data + ourNames = self.data['objects'] + for name, data in otherdata['objects'].items(): + # data == (docname, type) + if data[0] in docnames and name not in ourNames: + ourNames[name] = data + FortranDomain.merge_domaindata = merge_domaindata +_patch_sphinx_fortran_merge_domaindata() + +# -- Monkey-patch: stop sphinx.util.math.wrap_displaymath from double-wrapping +# LaTeX environments that the source already provides -------------------- +# +# MOM6's documentation contains a lot of math written directly with explicit +# LaTeX environments (`\begin{equation}`, `\begin{eqnarray}`, `\begin{align}`) +# inside `.. math::` directives. By default, Sphinx's wrap_displaymath() takes +# the body of a single-part `.. math::` directive and wraps it in +# `\begin{split}...\end{split}` and then again in `\begin{equation}` (or the +# starred unnumbered form). When the body already contains its own +# `\begin{equation}` etc., that produces nested LaTeX environments and +# pdflatex chokes. +# +# Sphinx's official workaround is the `:nowrap:` option on every affected +# `.. math::` directive, which would mean editing every math directive +# upstream in the MOM6 source tree. The jr3cermak/sphinx fork that the docs +# build previously depended on instead patched wrap_displaymath() so that any +# part containing one of these begin-environments is emitted verbatim and the +# outer wrapping is suppressed. We replicate that behavior here as a +# function-level monkey-patch on stock upstream Sphinx, so we can build +# against unmodified Sphinx 8.x. +# +# The detection is intentionally a plain substring search, matching the +# fork's behavior — `begin{equation`, `begin{eqnarray`, and `begin{align` all +# also match their starred and `aligned`/`equation*` variants. This is the +# same heuristic the fork shipped and what the existing MOM6 sources expect. +# +# Sphinx upstream issue tracking the same problem has been open since 2017 +# (sphinx-doc/sphinx#3785). A faithful upstream PR would be the right +# long-term fix, but that conversation is much older than this MOM6 upgrade +# work, so we are not blocking on it. +# +# TODO(piece-3): consider submitting a cleaner version upstream and dropping +# this patch when/if it lands. +def _patch_sphinx_wrap_displaymath(): + import sphinx.util.math as _sm + + def wrap_displaymath(text, label, numbering): + def is_equation(part): + return part.strip() + + if label is None: + labeldef = '' + else: + labeldef = r'\label{%s}' % label + numbering = True + + parts = list(filter(is_equation, text.split('\n\n'))) + + # Detect parts that already supply their own LaTeX environment. + nowrap = any( + ('begin{equation' in p) or + ('begin{eqnarray' in p) or + ('begin{align' in p) + for p in parts + ) + + equations = [] + if len(parts) == 0: + return '' + elif len(parts) == 1: + if numbering: + begin = r'\begin{equation}' + labeldef + end = r'\end{equation}' + else: + begin = r'\begin{equation*}' + labeldef + end = r'\end{equation*}' + if nowrap: + equations.append('%s\n' % parts[0]) + else: + equations.append('\\begin{split}%s\\end{split}\n' % parts[0]) + else: + if numbering: + begin = r'\begin{align}%s\!\begin{aligned}' % labeldef + end = r'\end{aligned}\end{align}' + else: + begin = r'\begin{align*}%s\!\begin{aligned}' % labeldef + end = r'\end{aligned}\end{align*}' + equations.extend('%s\\\\\n' % part.strip() for part in parts) + + if nowrap: + begin = '' + end = '' + + return '%s\n%s%s' % (begin, ''.join(equations), end) + + _sm.wrap_displaymath = wrap_displaymath +_patch_sphinx_wrap_displaymath() + def setup(app): app.add_config_value('sphinx_build_mode', '', 'env') app.add_role('latex', latexPassthru) @@ -137,13 +270,18 @@ def latexPassthru(name, rawtext, text, lineno, inliner, options={}, content=[]): extensions = [ 'sphinxcontrib.bibtex', 'sphinx.ext.ifconfig', - 'sphinxcontrib.autodoc_doxygen', + 'autodoc_doxygen', 'sphinxfortran.fortran_domain', ] bibtex_bibfiles = ['ocean.bib', 'references.bib', 'zotero.bib'] autosummary_generate = ['api/modules.rst', 'api/pages.rst'] -doxygen_xml = 'xml' +# Absolute path so the autodoc_doxygen extension can find the doxygen XML +# output regardless of what cwd Sphinx has at builder-inited time. This +# previously broke on RTD, where `sphinx-build -M html docs ...` runs from +# the repo root rather than from `docs/`, and the extension's os.path.isdir +# check resolved "xml" against the wrong directory. +doxygen_xml = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'xml') # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -182,7 +320,17 @@ def latexPassthru(name, rawtext, text, lineno, inliner, options={}, content=[]): # List of patterns, relative to source directory, that match files and # directories to ignore when looking for source files. -exclude_patterns = ['_build', 'details', 'src', 'Thumbs.db', '.DS_Store'] +exclude_patterns = [ + '_build', '_build.*', + 'details', 'src', 'Thumbs.db', '.DS_Store', + # Local virtualenvs that may sit alongside the docs source. Sphinx walks + # the entire source tree by default and otherwise picks up LICENSE.rst, + # README.rst, autosummary template files, etc. from inside site-packages + # and reports them as "isn't included in any toctree". + 'venv', 'venv.*', 'venv-*', + # Vendored extension's Jinja2 templates are not real .rst documents. + '_ext/*/templates', '_ext/*/*/templates', +] # The reST default role (used for this markup: `text`) to use for all # documents. @@ -241,7 +389,7 @@ def latexPassthru(name, rawtext, text, lineno, inliner, options={}, content=[]): # 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'] +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 diff --git a/docs/requirements.txt b/docs/requirements.txt index 539de3fc0e..ef24c6989b 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,20 +1,16 @@ -git+https://github.com/jr3cermak/sphinx.git@v3.2.1mom6.4 -git+https://github.com/jr3cermak/sphinxcontrib-autodoc_doxygen.git@0.7.13#egg=sphinxcontrib-autodoc_doxygen -git+https://github.com/jr3cermak/sphinx-fortran.git@1.2.2#egg=sphinx-fortran -git+https://github.com/jr3cermak/flint.git@0.0.1#egg=flint +sphinx>=8,<9 sphinx-rtd-theme sphinxcontrib-bibtex -# requirements.txt not working from sphinx-fortran +# lxml is required by the vendored docs/_ext/autodoc_doxygen extension +# (used to parse Doxygen XML output). +lxml numpy + +# Upstream VACUMM/sphinx-fortran has not cut a PyPI release past 1.1.1 but +# its master branch has had continued fixes through 2025. Pinned to a +# specific commit for reproducibility. See Notes-sphinx-upgrade.md (piece 2). +# NB: this commit contains a broken FortranDomain.merge_domaindata that we +# monkey-patch in conf.py to keep parallel builds (-j > 1) working. sphinx- +# fortran also still imports `six` at module load time on this commit. +git+https://github.com/VACUMM/sphinx-fortran.git@b14f438c1cc74d1dbcd5acd9a330c3b509caab56#egg=sphinx-fortran six -future -# Old Sphinx requires an old Jinja2 -jinja2<3.1 -sphinxcontrib_applehelp<1.0.8 -sphinxcontrib_devhelp<1.0.6 -sphinxcontrib_htmlhelp<2.0.5 -sphinxcontrib_qthelp<1.0.7 -sphinxcontrib_serializinghtml<1.0.7 -alabaster<0.7.14 -# setuptools 82.0.0 removed pkg_resources -setuptools<82.0.0 \ No newline at end of file diff --git a/src/ALE/MOM_remapping.F90 b/src/ALE/MOM_remapping.F90 index 78d31a48f1..d32b957717 100644 --- a/src/ALE/MOM_remapping.F90 +++ b/src/ALE/MOM_remapping.F90 @@ -2102,23 +2102,23 @@ logical function remapping_unit_tests(verbose, num_comp_samp) integer :: om4 ! Loop parameter, 0 or 1 integer :: ntests ! Number of iterations when brute force testing character(len=4) :: om4_tag ! Generated label - type(PCM) :: PCM - type(PLM_CW) :: PLM_CW - type(PLM_hybgen) :: PLM_hybgen - type(MPLM_WA) :: MPLM_WA - type(EMPLM_WA) :: EMPLM_WA - type(MPLM_WA_poly) :: MPLM_WA_poly - type(EMPLM_WA_poly) :: EMPLM_WA_poly - type(PLM_CWK) :: PLM_CWK - type(MPLM_CWK) :: MPLM_CWK - type(EMPLM_CWK) :: EMPLM_CWK - type(PPM_H4_2019) :: PPM_H4_2019 - type(PPM_H4_2018) :: PPM_H4_2018 - type(PPM_CW) :: PPM_CW - type(PPM_hybgen) :: PPM_hybgen - type(PPM_CWK) :: PPM_CWK - type(EPPM_CWK) :: EPPM_CWK - type(PLM_WLS) :: PLM_WLS + type(PCM) :: PCM_instance + type(PLM_CW) :: PLM_CW_instance + type(PLM_hybgen) :: PLM_hybgen_instance + type(MPLM_WA) :: MPLM_WA_instance + type(EMPLM_WA) :: EMPLM_WA_instance + type(MPLM_WA_poly) :: MPLM_WA_poly_instance + type(EMPLM_WA_poly) :: EMPLM_WA_poly_instance + type(PLM_CWK) :: PLM_CWK_instance + type(MPLM_CWK) :: MPLM_CWK_instance + type(EMPLM_CWK) :: EMPLM_CWK_instance + type(PPM_H4_2019) :: PPM_H4_2019_instance + type(PPM_H4_2018) :: PPM_H4_2018_instance + type(PPM_CW) :: PPM_CW_instance + type(PPM_hybgen) :: PPM_hybgen_instance + type(PPM_CWK) :: PPM_CWK_instance + type(EPPM_CWK) :: EPPM_CWK_instance + type(PLM_WLS) :: PLM_WLS_instance call test%set( verbose=verbose ) ! Sets the verbosity flag in test ! call test%set( stop_instantly=.true. ) ! While debugging @@ -2732,23 +2732,23 @@ logical function remapping_unit_tests(verbose, num_comp_samp) 3, (/0.,0.,0./), (/0.,0.,0./) ) if (verbose) write(test%stdout,*) '- - - - - - - - - - Recon1d PCM tests - - - - - - - - -' - call test%test( PCM%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test') - call test%test( MPLM_WA%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA unit test') - call test%test( EMPLM_WA%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA unit test') - call test%test( MPLM_WA_poly%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA_poly unit test') - call test%test( EMPLM_WA_poly%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA_poly unit test') - call test%test( PLM_hybgen%unit_tests(verbose, test%stdout, test%stderr), 'PLM_hybgen unit test') - call test%test( PLM_CW%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CW unit test') - call test%test( PLM_CWK%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CWK unit test') - call test%test( MPLM_CWK%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_CWK unit test') - call test%test( EMPLM_CWK%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_CWK unit test') - call test%test( PPM_H4_2019%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test') - call test%test( PPM_H4_2018%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2018 unit test') - call test%test( PPM_hybgen%unit_tests(verbose, test%stdout, test%stderr), 'PPM_hybgen unit test') - call test%test( PPM_CW%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CW unit test') - call test%test( PPM_CWK%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CWK unit test') - call test%test( EPPM_CWK%unit_tests(verbose, test%stdout, test%stderr), 'EPPM_CWK unit test') - call test%test( PLM_WLS%unit_tests(verbose, test%stdout, test%stderr), 'PLM_WLS unit test') + call test%test( PCM_instance%unit_tests(verbose, test%stdout, test%stderr), 'PCM unit test') + call test%test( MPLM_WA_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA unit test') + call test%test( EMPLM_WA_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA unit test') + call test%test( MPLM_WA_poly_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_WA_poly unit test') + call test%test( EMPLM_WA_poly_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_WA_poly unit test') + call test%test( PLM_hybgen_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_hybgen unit test') + call test%test( PLM_CW_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CW unit test') + call test%test( PLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_CWK unit test') + call test%test( MPLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'MPLM_CWK unit test') + call test%test( EMPLM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'EMPLM_CWK unit test') + call test%test( PPM_H4_2019_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2019 unit test') + call test%test( PPM_H4_2018_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_H4_2018 unit test') + call test%test( PPM_hybgen_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_hybgen unit test') + call test%test( PPM_CW_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CW unit test') + call test%test( PPM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'PPM_CWK unit test') + call test%test( EPPM_CWK_instance%unit_tests(verbose, test%stdout, test%stderr), 'EPPM_CWK unit test') + call test%test( PLM_WLS_instance%unit_tests(verbose, test%stdout, test%stderr), 'PLM_WLS unit test') ! Randomized, brute force tests ntests = 3000 diff --git a/src/ALE/Recon1d_PLM_WLS.F90 b/src/ALE/Recon1d_PLM_WLS.F90 index 24d6988f24..f32ff59b73 100644 --- a/src/ALE/Recon1d_PLM_WLS.F90 +++ b/src/ALE/Recon1d_PLM_WLS.F90 @@ -231,8 +231,7 @@ logical function check_reconstruction(this, h, u) enddo ! Create a perturbable reconstruction - call perturbed%init( this%n, h_neglect=this%h_neglect ) - call perturbed%reconstruct( h, u ) ! Should reproduce "this" + perturbed = this ! Complete copy of this ! Check the copy is identical do k = 1, this%n if ( abs( perturbed%u_mean(k) - this%u_mean(k) ) > 0. ) check_reconstruction = .true. @@ -240,6 +239,7 @@ logical function check_reconstruction(this, h, u) if ( abs( perturbed%ur(k) - this%ur(k) ) > 0. ) check_reconstruction = .true. if ( abs( perturbed%slp(k) - this%slp(k) ) > 0. ) check_reconstruction = .true. enddo + ! The !DIR$ NOINLINE directive would be needed here to avoid ifort -O2 changing answers ! Now perturb the slope. The local error should not decrease. do k = 1, this%n slp = this%slp(k) * ( 1.0 + 1. * epsilon(slp) ) @@ -272,6 +272,7 @@ real function LS_error(this, k, h, u) real :: h_l0, h_r0, hc0 ! Thickness of left, right, center cells with h_neglect added [H] real :: hx2l, hx2r ! Contributions to denominator, [H3] real :: hxyl, hxyr ! Contributions to numerator, [H2 A] + real :: slp ! The PLM slopes (difference across cell) [A] integer :: km1, kp1 km1 = max(1, k-1) @@ -287,11 +288,12 @@ real function LS_error(this, k, h, u) h_l0 = h_l + this%h_neglect h_r0 = h_r + this%h_neglect - hxyl = ( h_l * 0.5 * ( h_c + h_l ) ) * ( u_c - u_l ) - hxyr = ( h_r * 0.5 * ( h_c + h_r ) ) * ( u_r - u_c ) - hx2l = h_l0 * 0.25 * ( h_c + h_l0 )**2 - hx2r = h_r0 * 0.25 * ( h_c + h_r0 )**2 - LS_error = h_c * ( ( hx2l + hx2r ) * this%slp(k) - h(k) * ( hxyl + hxyr ) )**2 + hxyl = ( h_l * ( h_c + h_l ) ) * ( u_c - u_l ) + hxyr = ( h_r * ( h_c + h_r ) ) * ( u_r - u_c ) + hx2l = h_l0 * ( h_c + h_l0 )**2 + hx2r = h_r0 * ( h_c + h_r0 )**2 + slp = 2. * h_c * ( hxyr + hxyl ) / ( hx2l + hx2r ) + LS_error = h_c * ( ( hx2l + hx2r ) * ( this%slp(k) - slp ) )**2 LS_error = LS_error / ( hc0 * ( hx2l + hx2r ) ) end function LS_error @@ -446,22 +448,22 @@ end function unit_tests !! then !! \f{align}{ !! G(b) &= 2 < h e^2 > \\\\ -!! &= 2 b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 4 b \frac{ < h x' y' > }{ h_{k} } + 2 < h' {y'}^2 > +!! &= 2 b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 4 b \frac{ < h x' y' > }{ h_{k} } + 2 < h {y'}^2 > !! \;\; . !! \f} !! !! If we denote the value of b that yields the minimum value as \f$ b^* \f$ then !! \f[ -!! G(b^*) = < h {y'}^2 > - \frac{ < h x' y' >^2 }{ < h {x'}^2 > } +!! G(b^*) = 2 < h {y'}^2 > - \frac{ 2 < h x' y' >^2 }{ < h {x'}^2 > } !! \;\; . !! \f] !! !! Let !! \f{align}{ !! G''(b) &= G(b) - G(b^*) \\\\ -!! &= b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 2 b \frac{ < h x' y' > }{ h_{k} } -!! + \frac{ < h x' y' > }{ < h {x'}^2 > } \\\\ -!! &= \frac{ \left( b < h {x'}^2 > - h_{k} < h x' y' > \right)^2 }{ h_{k} < h {x'}^2 > } +!! &= 2 b^2 \frac{ < h {x'}^2 > }{ h_{k}^2 } - 4 b \frac{ < h x' y' > }{ h_{k} } +!! + 2 \frac{ < h x' y' >^2 }{ < h {x'}^2 > } \\\\ +!! &= 2 \frac{ \left( b < h {x'}^2 > - h_{k} < h x' y' > \right)^2 }{ h_{k} < h {x'}^2 > } !! \;\; . !! \f} !! Minimizing \f$ G''(b) \f$ is equivalent to minimizing \f$ G(b) \f$ for the same data. diff --git a/src/ALE/Recon1d_PPM_H4_2019.F90 b/src/ALE/Recon1d_PPM_H4_2019.F90 index 26985be644..2dc3315eaf 100644 --- a/src/ALE/Recon1d_PPM_H4_2019.F90 +++ b/src/ALE/Recon1d_PPM_H4_2019.F90 @@ -365,10 +365,10 @@ subroutine end_value_h4(dz, u, Csys) Wt(2,4) = -4.0 * I_h1234 * (I_h23 * (I_h123 + I_h234)) ! Wt*h1^3 > -4* (h1/h23)*(1+h1/h234) Wt(3,4) = 4.0 * I_denom ! = 4.0*I_h1234 * I_h234 * I_h34 ! Wt*h1^3 < 4 * (h1/h234)*(h1/h34) - Csys(1) = ((u(1) + Wt(1,1) * (u(2)-u(1))) + Wt(2,1) * (u(3)-u(2))) + Wt(3,1) * (u(4)-u(3)) - Csys(2) = (Wt(1,2) * (u(2)-u(1)) + Wt(2,2) * (u(3)-u(2))) + Wt(3,2) * (u(4)-u(3)) - Csys(3) = (Wt(1,3) * (u(2)-u(1)) + Wt(2,3) * (u(3)-u(2))) + Wt(3,3) * (u(4)-u(3)) - Csys(4) = (Wt(1,4) * (u(2)-u(1)) + Wt(2,4) * (u(3)-u(2))) + Wt(3,4) * (u(4)-u(3)) + Csys(1) = ((u(1) + (Wt(1,1) * (u(2)-u(1)))) + (Wt(2,1) * (u(3)-u(2)))) + (Wt(3,1) * (u(4)-u(3))) + Csys(2) = ((Wt(1,2) * (u(2)-u(1))) + (Wt(2,2) * (u(3)-u(2)))) + (Wt(3,2) * (u(4)-u(3))) + Csys(3) = ((Wt(1,3) * (u(2)-u(1))) + (Wt(2,3) * (u(3)-u(2)))) + (Wt(3,3) * (u(4)-u(3))) + Csys(4) = ((Wt(1,4) * (u(2)-u(1))) + (Wt(2,4) * (u(3)-u(2)))) + (Wt(3,4) * (u(4)-u(3))) ! endif ! End of non-uniform layer thickness branch. diff --git a/src/ALE/Recon1d_PPM_hybgen.F90 b/src/ALE/Recon1d_PPM_hybgen.F90 index 9d26e27eed..343da5d8e8 100644 --- a/src/ALE/Recon1d_PPM_hybgen.F90 +++ b/src/ALE/Recon1d_PPM_hybgen.F90 @@ -28,7 +28,7 @@ module Recon1d_PPM_hybgen !! The source for the methods ultimately used by this class are: !! - init() -> recon1d_ppm_cw.init() !! - reconstruct() *locally defined -!! - average() -> recon1d_ppm_cw.average() +!! - average() *locally defined but calls recon1d_ppm_cw.average() !! - f() -> recon1d_ppm_cw.f() !! - dfdx() -> recon1d_ppm_cw.dfdx() !! - check_reconstruction() *locally defined @@ -42,6 +42,8 @@ module Recon1d_PPM_hybgen contains !> Implementation of the PPM_hybgen reconstruction procedure :: reconstruct => reconstruct + !> Implementation of the PPM_hybgen average over an interval [A] + procedure :: average => average !> Implementation of check reconstruction for the PPM_hybgen reconstruction procedure :: check_reconstruction => check_reconstruction !> Implementation of unit tests for the PPM_hybgen reconstruction @@ -51,212 +53,85 @@ module Recon1d_PPM_hybgen contains -!> Calculate a 1D PPM_hybgen reconstructions based on h(:) and u(:) +!> Calculate a 1D PPM_hybgen reconstruction based on h(:) and u(:) +!! +!! First pass: hybgen_ppm_coefs() computes initial edge estimates with CW monotonicity. +!! Second pass: applies OM4-era bound_edge_values() and check_discontinuous_edge_values(), +!! then the standard CW PPM limiter (post-2018 expressions, answer_date=99991231). +!! This reproduces bit-for-bit the behavior of the old-style PPM_HYBGEN scheme. subroutine reconstruct(this, h, u) class(PPM_hybgen), intent(inout) :: this !< This reconstruction real, intent(in) :: h(*) !< Grid spacing (thickness) [typically H] real, intent(in) :: u(*) !< Cell mean values [A] ! Local variables - real :: h0, h1, h2, h3 ! Cell thickness h(k-2), h(k-1), h(k), h(k+1) in K loop [H] - real :: h01_h112, h23_h122 ! Approximately 2/3 [nondim] - real :: h112, h122 ! Approximately 3 h [H] - real :: ddh ! Approximately 0 [nondim] - real :: I_h12, I_h01, I_h0123 ! Reciprocals of d12 and sum(h) [H-1] - real :: dul, dur ! Left and right cell PLM slopes [A] - real :: u0, u1, u2 ! Far left, left, and right cell values [A] - real :: edge ! Edge value between cell k-1 and k [A] - real :: u_min, u_max ! Minimum and maximum value across edge [A] - real :: a6 ! Colella and Woodward curvature [A] - real :: du, duc ! Difference between edges across cell [A] - real :: slp(this%n) ! PLM slope [A] - real :: sigma_l, sigma_c, sigma_r ! Left, central and right slope estimates as - ! differences across the cell [A] - real :: slope_x_h ! retained PLM slope times half grid step [A] - real :: edge_l, edge_r ! Edge values (left and right) [A] - real :: expr1, expr2 ! Temporary expressions [A2] - real :: u0_avg ! avg value at given edge [A] - integer :: k, n, km1, kp1 + integer :: k, n + real :: ppoly_e(this%n, 2) ! PPM edge values [A] + real :: u_l, u_c, u_r ! Left, center, right cell averages [A] + real :: edge_l, edge_r ! Left and right edge values [A] + real :: expr1, expr2 ! Temporary expressions [A2] n = this%n - ! First populate the PLM reconstructions - slp(1) = 0. - do k = 2, n-1 - h0 = max( this%h_neglect, h(k-1) ) - h1 = max( this%h_neglect, h(k) ) - h2 = max( this%h_neglect, h(k+1) ) - dul = u(k) - u(k-1) - dur = u(k+1) - u(k) - h112 = ( 2.0 * h0 + h1 ) - h122 = ( h1 + 2.0 * h2 ) - I_h01 = 1. / ( h0 + h1 ) - I_h12 = 1. / ( h1 + h2 ) - h01_h112 = ( 2.0 * h0 + h1 ) / ( h0 + h1 ) ! When uniform -> 3/2 - h23_h122 = ( 2.0 * h2 + h1 ) / ( h2 + h1 ) ! When uniform -> 3/2 - if ( dul * dur > 0.) then - du = ( h1 / ( h1 + ( h0 + h2 ) ) ) * ( h112 * dur * I_h12 + h122 * dul * I_h01 ) - slp(k) = sign( min( abs(2.0 * dul), abs(du), abs(2.0 * dur) ), du) - else - slp(k) = 0. - endif - enddo - slp(n) = 0. - - this%ul(1) = u(1) ! PCM - this%ur(1) = u(1) ! PCM - this%ul(2) = u(1) ! PCM - do K = 3, n-1 ! K=3 is interface between cells 2 and 3 - h0 = max( this%h_neglect, h(k-2) ) - h1 = max( this%h_neglect, h(k-1) ) - h2 = max( this%h_neglect, h(k) ) - h3 = max( this%h_neglect, h(k+1) ) - h01_h112 = ( h0 + h1 ) / ( 2. * h1 + h2 ) ! When uniform -> 2/3 - h23_h122 = ( h2 + h3 ) / ( h1 + 2. * h2 ) ! When uniform -> 2/3 - ddh = h01_h112 - h23_h122 ! When uniform -> 0 - I_h12 = 1.0 / ( h1 + h2 ) ! When uniform -> 1/(2h) - I_h0123 = 1.0 / ( ( h0 + h1 ) + ( h2 + h3 ) ) ! When uniform -> 1/(4h) - dul = slp(k-1) - dur = slp(k) - u1 = u(k-1) - u2 = u(k) - edge = I_h12 * ( h2 * u1 + h1 * u2 ) & ! 1/2 u1 + 1/2 u2 - + I_h0123 * ( 2.0 * h1 * h2 * I_h12 * ( u2 - u1 ) * ddh & ! 0 - + ( h2 * dul * h23_h122 - h1 * dur * h01_h112 ) ) ! 1/6 dul - 1/6 dur - this%ur(k-1) = edge - this%ul(k) = edge - enddo - this%ur(n-1) = u(n) ! PCM - this%ur(n) = u(n) ! PCM - this%ul(n) = u(n) ! PCM - - do K = 2, n-1 ! K=2 is interface between cells 1 and 2 - u0 = u(k-1) - u1 = u(k) - u2 = u(k+1) - a6 = 3.0 * ( ( u1 - this%ul(k) ) + ( u1 - this%ur(k) ) ) - a6 = 6.0 * u1 - 3.0 * ( this%ul(k) + this%ur(k) ) - du = this%ur(k) - this%ul(k) - if ( ( u2 - u1 ) * ( u1 - u0 ) <= 0.0 ) then ! Large scale extrema - this%ul(k) = u1 - this%ur(k) = u1 - elseif ( du * a6 > du * du ) then ! Extrema on right - edge = 3.0 * u1 - 2.0 * this%ur(k) ! Subject to round off - ! u_min = min( u0, u1 ) - ! u_max = max( u0, u1 ) - ! edge = max( min( edge, u_max), u_min ) - this%ul(k) = edge - elseif ( du * a6 < - du * du ) then ! Extrema on left - edge = 3.0 * u1 - 2.0 * this%ul(k) ! Subject to round off - ! u_min = min( u1, u2 ) - ! u_max = max( u1, u2 ) - ! edge = max( min( edge, u_max), u_min ) - this%ur(k) = edge - endif - enddo - - ! ### Note that the PPM_HYBGEM option calculated the CW PPM coefficients and then - ! invoked the OM4-era limiters afterwards, effectively doing the limiters twice. - ! This second pass does change answers! + ! First pass: compute initial edge estimates using the hybgen algorithm with CW monotonicity + call hybgen_ppm_coefs(u, h, ppoly_e, n, this%h_neglect) - ! Loop on cells to bound edge value - do k = 1, n + ! Second pass: apply OM4-era PPM limiters (post-2018 answers via answer_date=99991231) + call bound_edge_values(n, h, u, ppoly_e, this%h_neglect, answer_date=99991231) + call check_discontinuous_edge_values(n, u, ppoly_e) - ! For the sake of bounding boundary edge values, the left neighbor of the left boundary cell - ! is assumed to be the same as the left boundary cell and the right neighbor of the right - ! boundary cell is assumed to be the same as the right boundary cell. This effectively makes - ! boundary cells look like extrema. - km1 = max(1,k-1) ; kp1 = min(k+1,N) - - slope_x_h = 0.0 - sigma_l = ( u(k) - u(km1) ) - if ( (h(km1) + h(kp1)) + 2.0*h(k) > 0. ) then - sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) ) - else - sigma_c = 0. - endif - sigma_r = ( u(kp1) - u(k) ) - - ! The limiter is used in the local coordinate system to each cell, so for convenience store - ! the slope times a half grid spacing. (See White and Adcroft JCP 2008 Eqs 19 and 20) - if ( (sigma_l * sigma_r) > 0.0 ) & - slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) - - ! Limit the edge values - if ( (u(km1)-this%ul(k)) * (this%ul(k)-u(k)) < 0.0 ) then - this%ul(k) = u(k) - sign( min( abs(slope_x_h), abs(this%ul(k)-u(k)) ), slope_x_h ) - endif - - if ( (u(kp1)-this%ur(k)) * (this%ur(k)-u(k)) < 0.0 ) then - this%ur(k) = u(k) + sign( min( abs(slope_x_h), abs(this%ur(k)-u(k)) ), slope_x_h ) - endif - - ! Finally bound by neighboring cell means in case of roundoff - this%ul(k) = max( min( this%ul(k), max(u(km1), u(k)) ), min(u(km1), u(k)) ) - this%ur(k) = max( min( this%ur(k), max(u(kp1), u(k)) ), min(u(kp1), u(k)) ) - - enddo ! loop on interior edges - - do k = 1, n-1 - if ( (this%ul(k+1) - this%ur(k)) * (u(k+1) - u(k)) < 0.0 ) then - u0_avg = 0.5 * ( this%ur(k) + this%ul(k+1) ) - u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) ) - this%ur(k) = u0_avg - this%ul(k+1) = u0_avg - endif - enddo ! end loop on interior edges - - ! Loop on interior cells to apply the standard - ! PPM limiter (Colella & Woodward, JCP 84) + ! Apply the standard CW PPM limiter (Colella & Woodward, JCP 84) on interior cells do k = 2, n-1 - - ! Get cell averages - u0 = u(k-1) - u1 = u(k) - u2 = u(k+1) - - edge_l = this%ul(k) - edge_r = this%ur(k) - - if ( (u2 - u1)*(u1 - u0) <= 0.0) then - ! Flatten extremum - edge_l = u1 - edge_r = u1 + u_l = u(k-1) ; u_c = u(k) ; u_r = u(k+1) + edge_l = ppoly_e(k,1) ; edge_r = ppoly_e(k,2) + if ( (u_r - u_c)*(u_c - u_l) <= 0.0 ) then + edge_l = u_c ; edge_r = u_c else - expr1 = 3.0 * (edge_r - edge_l) * ( (u1 - edge_l) + (u1 - edge_r)) + expr1 = 3.0 * (edge_r - edge_l) * ( (u_c - edge_l) + (u_c - edge_r) ) expr2 = (edge_r - edge_l) * (edge_r - edge_l) if ( expr1 > expr2 ) then - ! Place extremum at right edge of cell by adjusting left edge value - edge_l = u1 + 2.0 * ( u1 - edge_r ) - edge_l = max( min( edge_l, max(u0, u1) ), min(u0, u1) ) ! In case of round off + edge_l = u_c + 2.0 * ( u_c - edge_r ) + edge_l = max( min( edge_l, max(u_l, u_c) ), min(u_l, u_c) ) elseif ( expr1 < -expr2 ) then - ! Place extremum at left edge of cell by adjusting right edge value - edge_r = u1 + 2.0 * ( u1 - edge_l ) - edge_r = max( min( edge_r, max(u2, u1) ), min(u2, u1) ) ! In case of round off + edge_r = u_c + 2.0 * ( u_c - edge_l ) + edge_r = max( min( edge_r, max(u_r, u_c) ), min(u_r, u_c) ) endif endif - ! This checks that the difference in edge values is representable - ! and avoids overshoot problems due to round off. - !### The 1.e-60 needs to have units of [A], so this dimensionally inconsistent. - if ( abs( edge_r - edge_l ) Average between xa and xb for cell k of a PPM_hybgen reconstruction [A] +!! +!! Calls the parent PPM_CW average and then clamps the result to [min(ul,ur), max(ul,ur)]. +!! This replicates the force_bounds_in_subcell behavior of the equivalent old-style PPM_HYBGEN +!! scheme. +real function average(this, k, xa, xb) + class(PPM_hybgen), intent(in) :: this !< This reconstruction + integer, intent(in) :: k !< Cell number + real, intent(in) :: xa !< Start of averaging interval on element (0 to 1) + real, intent(in) :: xb !< End of averaging interval on element (0 to 1) + real :: u_lo, u_hi ! Bounds on the sub-cell average given by the edge values [A] + + average = this%PPM_CW%average(k, xa, xb) + u_lo = min(this%ul(k), this%ur(k)) + u_hi = max(this%ul(k), this%ur(k)) + average = max(u_lo, min(u_hi, average)) + +end function average + !> Checks the PPM_hybgen reconstruction for consistency logical function check_reconstruction(this, h, u) class(PPM_hybgen), intent(in) :: this !< This reconstruction @@ -402,4 +277,177 @@ end function unit_tests !> \namespace recon1d_ppm_hybgen !! +! ============================================================================ +! Private subroutines copied from phased-out modules to avoid dependencies. +! These reproduce bit-for-bit the results of the original functions they replace. +! ============================================================================ + +!> Set up edge values for PPM reconstruction using the hybgen (HYCOM) algorithm. +!! +!! Copied from MOM_hybgen_remap.hybgen_ppm_coefs(). +!! Original code by Tim Campbell (MSU, 2002) and Alan Wallcraft (NRL, 2007). +subroutine hybgen_ppm_coefs(s, h_src, edges, nk, thin, PCM_lay) + integer, intent(in) :: nk !< The number of input layers + real, intent(in) :: s(nk) !< The input scalar fields [A] + real, intent(in) :: h_src(nk) !< The input grid layer thicknesses [H ~> m or kg m-2] + real, intent(out) :: edges(nk,2) !< The PPM interpolation edge values [A] + real, intent(in) :: thin !< A negligible layer thickness [H ~> m or kg m-2] + logical, optional, intent(in) :: PCM_lay(nk) !< If true for a layer, use PCM remapping + + real :: dp(nk) ! Input grid layer thicknesses, but with a minimum thickness given by thin [H ~> m or kg m-2] + logical :: PCM_layer(nk) ! True for layers that should use PCM remapping + real :: da ! Difference between the unlimited scalar edge value estimates [A] + real :: a6 ! Scalar field differences that are proportional to the curvature [A] + real :: slk, srk ! Differences between adjacent cell averages of scalars [A] + real :: sck ! Scalar differences across a cell [A] + real :: as(nk) ! Scalar field difference across each cell [A] + real :: al(nk), ar(nk) ! Scalar field at the left and right edges of a cell [A] + real :: h112(nk+1), h122(nk+1) ! Combinations of thicknesses [H ~> m or kg m-2] + real :: I_h12(nk+1) ! Inverses of combinations of thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: h2_h123(nk) ! A ratio of a layer thickness to the sum of 3 adjacent thicknesses [nondim] + real :: I_h0123(nk) ! Inverse of the sum of 4 adjacent thicknesses [H-1 ~> m-1 or m2 kg-1] + real :: h01_h112(nk+1) ! A ratio of sums of adjacent thicknesses [nondim] + real :: h23_h122(nk+1) ! A ratio of sums of adjacent thicknesses [nondim] + integer :: k + + do k=1,nk ; dp(k) = max(h_src(k), thin) ; enddo + + if (present(PCM_lay)) then + do k=1,nk ; PCM_layer(k) = (PCM_lay(k) .or. dp(k) <= thin) ; enddo + else + do k=1,nk ; PCM_layer(k) = (dp(k) <= thin) ; enddo + endif + + do k=2,nk + h112(K) = 2.*dp(k-1) + dp(k) + h122(K) = dp(k-1) + 2.*dp(k) + I_h12(K) = 1.0 / (dp(k-1) + dp(k)) + enddo + do k=2,nk-1 + h2_h123(k) = dp(k) / (dp(k) + (dp(k-1)+dp(k+1))) + enddo + do K=3,nk-1 + I_h0123(K) = 1.0 / ((dp(k-2) + dp(k-1)) + (dp(k) + dp(k+1))) + h01_h112(K) = (dp(k-2) + dp(k-1)) / (2.0*dp(k-1) + dp(k)) + h23_h122(K) = (dp(k) + dp(k+1)) / (dp(k-1) + 2.0*dp(k)) + enddo + + as(1) = 0. + do k=2,nk-1 + if (PCM_layer(k)) then + as(k) = 0.0 + else + slk = s(k)-s(k-1) + srk = s(k+1)-s(k) + if (slk*srk > 0.) then + sck = h2_h123(k)*( h112(K)*srk*I_h12(K+1) + h122(K+1)*slk*I_h12(K) ) + as(k) = sign(min(abs(2.0*slk), abs(sck), abs(2.0*srk)), sck) + else + as(k) = 0. + endif + endif + enddo + as(nk) = 0. + al(1) = s(1) + ar(1) = s(1) + al(2) = s(1) + do K=3,nk-1 + al(k) = (dp(k)*s(k-1) + dp(k-1)*s(k)) * I_h12(K) & + + I_h0123(K)*( 2.*dp(k)*dp(k-1)*I_h12(K)*(s(k)-s(k-1)) * & + ( h01_h112(K) - h23_h122(K) ) & + + (dp(k)*as(k-1)*h23_h122(K) - dp(k-1)*as(k)*h01_h112(K)) ) + ar(k-1) = al(k) + enddo + ar(nk-1) = s(nk) + al(nk) = s(nk) + ar(nk) = s(nk) + do k=2,nk-1 + if ((PCM_layer(k)) .or. ((s(k+1)-s(k))*(s(k)-s(k-1)) <= 0.)) then + al(k) = s(k) + ar(k) = s(k) + else + da = ar(k)-al(k) + a6 = 6.0*s(k) - 3.0*(al(k)+ar(k)) + if (da*a6 > da*da) then + al(k) = 3.0*s(k) - 2.0*ar(k) + elseif (da*a6 < -da*da) then + ar(k) = 3.0*s(k) - 2.0*al(k) + endif + endif + enddo + do k=1,nk + edges(k,1) = al(k) + edges(k,2) = ar(k) + enddo + +end subroutine hybgen_ppm_coefs + +!> Bound edge values by the averages of the neighboring cells. +!! +!! Copied from regrid_edge_values.bound_edge_values(). +subroutine bound_edge_values(N, h, u, edge_val, h_neglect, answer_date) + integer, intent(in) :: N !< Number of cells + real, dimension(N), intent(in) :: h !< Cell widths [H] + real, dimension(N), intent(in) :: u !< Cell averages [A] + real, dimension(N,2), intent(inout) :: edge_val !< Edge values [A] + real, intent(in) :: h_neglect !< A negligibly small width [H] + integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use + + real :: sigma_l, sigma_c, sigma_r + real :: slope_x_h + logical :: use_2018_answers + integer :: k, km1, kp1 + + use_2018_answers = .true. ; if (present(answer_date)) use_2018_answers = (answer_date < 20190101) + + do k = 1,N + km1 = max(1,k-1) ; kp1 = min(k+1,N) + slope_x_h = 0.0 + if (use_2018_answers) then + sigma_l = 2.0 * ( u(k) - u(km1) ) / ( h(k) + h_neglect ) + sigma_c = 2.0 * ( u(kp1) - u(km1) ) / ( h(km1) + 2.0*h(k) + h(kp1) + h_neglect ) + sigma_r = 2.0 * ( u(kp1) - u(k) ) / ( h(k) + h_neglect ) + if ( (sigma_l * sigma_r) > 0.0 ) & + slope_x_h = 0.5 * h(k) * sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) + elseif ( ((h(km1) + h(kp1)) + 2.0*h(k)) > 0.0 ) then + sigma_l = ( u(k) - u(km1) ) + sigma_c = ( u(kp1) - u(km1) ) * ( h(k) / ((h(km1) + h(kp1)) + 2.0*h(k)) ) + sigma_r = ( u(kp1) - u(k) ) + if ( (sigma_l * sigma_r) > 0.0 ) & + slope_x_h = sign( min(abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c ) + endif + if ( (u(km1)-edge_val(k,1)) * (edge_val(k,1)-u(k)) < 0.0 ) then + edge_val(k,1) = u(k) - sign( min( abs(slope_x_h), abs(edge_val(k,1)-u(k)) ), slope_x_h ) + endif + if ( (u(kp1)-edge_val(k,2)) * (edge_val(k,2)-u(k)) < 0.0 ) then + edge_val(k,2) = u(k) + sign( min( abs(slope_x_h), abs(edge_val(k,2)-u(k)) ), slope_x_h ) + endif + edge_val(k,1) = max( min( edge_val(k,1), max(u(km1), u(k)) ), min(u(km1), u(k)) ) + edge_val(k,2) = max( min( edge_val(k,2), max(u(kp1), u(k)) ), min(u(kp1), u(k)) ) + enddo + +end subroutine bound_edge_values + +!> Replace discontinuous edge values with their average when not monotonic. +!! +!! Copied from regrid_edge_values.check_discontinuous_edge_values(). +subroutine check_discontinuous_edge_values(N, u, edge_val) + integer, intent(in) :: N !< Number of cells + real, dimension(N), intent(in) :: u !< Cell averages [A] + real, dimension(N,2), intent(inout) :: edge_val !< Edge values [A] + + integer :: k + real :: u0_avg + + do k = 1,N-1 + if ( (edge_val(k+1,1) - edge_val(k,2)) * (u(k+1) - u(k)) < 0.0 ) then + u0_avg = 0.5 * ( edge_val(k,2) + edge_val(k+1,1) ) + u0_avg = max( min( u0_avg, max(u(k), u(k+1)) ), min(u(k), u(k+1)) ) + edge_val(k,2) = u0_avg + edge_val(k+1,1) = u0_avg + endif + enddo + +end subroutine check_discontinuous_edge_values + end module Recon1d_PPM_hybgen diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index a9ae24464c..f26013c1a5 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1862,12 +1862,14 @@ subroutine ALE_regridding_and_remapping(CS, G, GV, US, u, v, h, tv, dtdia, Time_ call remap_dyn_split_RK2_aux_vars(G, GV, CS%dyn_split_RK2_CSp, h_old_u, h_old_v, h_new_u, h_new_v, CS%ALE_CSp) endif - if (associated(CS%OBC)) then + if (associated(CS%OBC) .or. associated(CS%visc%Kv_shear_Bu)) then call pass_var(h, G%Domain, complete=.false.) call pass_var(h_new, G%Domain, complete=.true.) - call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell) endif + if (associated(CS%OBC)) & + call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell) + call remap_vertvisc_aux_vars(G, GV, CS%visc, h, h_new, CS%ALE_CSp, CS%OBC) if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, clock=id_clock_pass, halo=1) diff --git a/src/core/MOM_boundary_update.F90 b/src/core/MOM_boundary_update.F90 index bc06bc5844..70e35f7274 100644 --- a/src/core/MOM_boundary_update.F90 +++ b/src/core/MOM_boundary_update.F90 @@ -44,6 +44,8 @@ module MOM_boundary_update logical :: use_shelfwave = .false. !< If true, use the shelfwave open boundary. logical :: use_dyed_channel = .false. !< If true, use the dyed channel open boundary. logical :: debug_OBCs = .false. !< If true, write verbose OBC values for debugging purposes. + logical :: value_update_bug = .true. !< If true, recover a bug that OBC segment data does not + !! update if all segments use 'value' and none uses 'file'. integer :: nk_OBC_debug = 0 !< The number of layers of OBC segment data to write out !! in full when DEBUG_OBCS is true. !>@{ Pointers to the control structures for named OBC specifications @@ -87,6 +89,9 @@ subroutine call_OBC_register(G, GV, US, param_file, CS, OBC, tr_Reg) call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "OBC_VALUE_UPDATE_BUG", CS%value_update_bug, & + "If true, recover a bug that OBC segment data does not update if all segments "//& + "use 'value' and none uses 'file'.", default=.true.) call get_param(param_file, mdl, "USE_FILE_OBC", CS%use_files, & "If true, use external files for the open boundary.", & default=.false.) @@ -169,10 +174,13 @@ subroutine update_OBC_data(OBC, G, GV, US, tv, h, CS, Time) call shelfwave_set_OBC_data(OBC, CS%shelfwave_OBC_CSp, G, GV, US, h, Time) if (CS%use_dyed_channel) & call dyed_channel_update_flow(OBC, CS%dyed_channel_OBC_CSp, G, GV, US, h, Time) - if (OBC%any_needs_IO_for_data .or. OBC%add_tide_constituents) then - call read_OBC_segment_data(G, GV, US, OBC, tv, h, Time) - call update_OBC_segment_data(G, GV, US, OBC, h, Time) + + if (.not. OBC%user_BCs_set_globally) then + if (OBC%any_needs_IO_for_data) call read_OBC_segment_data(G, GV, US, OBC, tv, h, Time) + if ((.not.CS%value_update_bug) .or. (OBC%any_needs_IO_for_data .or. OBC%add_tide_constituents)) & + call update_OBC_segment_data(G, GV, US, OBC, h, Time) endif + if (CS%debug_OBCs) call chksum_OBC_segments(OBC, G, GV, US, CS%nk_OBC_debug) end subroutine update_OBC_data diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index 9b1a9686bb..ee58249dd1 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -1338,10 +1338,11 @@ subroutine remap_dyn_split_RK2_aux_vars(G, GV, CS, h_old_u, h_old_v, h_new_u, h_ call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%u_av, CS%v_av) call pass_vector(CS%u_av, CS%v_av, G%Domain, complete=.false.) call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%CAu_pred, CS%CAv_pred) - call pass_vector(CS%CAu_pred, CS%CAv_pred, G%Domain, complete=.true.) + call pass_vector(CS%CAu_pred, CS%CAv_pred, G%Domain, complete=.false.) endif call ALE_remap_velocities(ALE_CSp, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, CS%diffu, CS%diffv) + call pass_vector(CS%diffu, CS%diffv, G%Domain, complete=.true.) end subroutine remap_dyn_split_RK2_aux_vars diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 0af4d64609..ee8d59cef3 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -6,7 +6,8 @@ module MOM_open_boundary use MOM_array_transform, only : rotate_array, rotate_array_pair -use MOM_coms, only : sum_across_PEs, Set_PElist, Get_PElist, PE_here, num_PEs +use MOM_coms, only : sum_across_PEs, any_across_PEs +use MOM_coms, only : Set_PElist, Get_PElist, PE_here, num_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_debugging, only : hchksum, uvchksum, chksum use MOM_diag_mediator, only : diag_ctrl, time_type @@ -229,11 +230,8 @@ module MOM_open_boundary logical :: on_pe !< true if any portion of the segment is located in this PE's data domain logical :: temp_segment_data_exists !< true if temperature data arrays are present logical :: salt_segment_data_exists !< true if salinity data arrays are present - real, allocatable :: Cg(:,:) !< The external gravity wave speed [L T-1 ~> m s-1] - !! at OBC-points. real, allocatable :: Htot(:,:) !< The total column thickness [H ~> m or kg m-2] at OBC-points. real, allocatable :: dZtot(:,:) !< The total column vertical extent [Z ~> m] at OBC segment faces. - real, allocatable :: h(:,:,:) !< The cell thickness [H ~> m or kg m-2] at OBC segment faces real, allocatable :: normal_vel(:,:,:) !< The layer velocity normal to the OB !! segment [L T-1 ~> m s-1]. real, allocatable :: tangential_vel(:,:,:) !< The layer velocity tangential to the OB segment @@ -336,7 +334,6 @@ module MOM_open_boundary logical :: update_OBC = .false. !< Is OBC data time-dependent logical :: update_OBC_seg_data = .false. !< Is it the time for OBC segment data update for fields that !! require less frequent update - logical :: needs_IO_for_data = .false. !< Is any i/o needed for OBCs on the current PE logical :: any_needs_IO_for_data = .false. !< Is any i/o needed for OBCs globally integer :: vorticity_config !< An integer indicating OBC relative vorticity configuration integer :: strain_config !< An integer indicating OBC strain configuration @@ -971,7 +968,6 @@ end subroutine open_boundary_setup_vert !> Determine which physical fields are required for this segment based on boundary-condition type !! and segment orientation. Also enable groups of physical fields required by tides or thermodynamics. !! Note the tidal group could be further narrowed based on modes. -!! This subroutine could turn into a TBP for OBC_segment_type. subroutine segment_determine_required_fields(segment, tides, temp_salt) type(OBC_segment_type), intent(inout) :: segment !< OBC segment logical, optional, intent(in) :: tides !< Switch for tidal variables @@ -1034,6 +1030,28 @@ integer function find_phys_field_index(name) endif ; enddo end function find_phys_field_index +!> Set global flag OBC%any_needs_IO_for_data. +subroutine OBC_any_IO(OBC) + type(ocean_OBC_type), intent(inout) :: OBC !< Open boundary control structure + + ! Local variables + integer :: m, n + logical :: use_IO + + use_IO = .false. + do n=1,OBC%number_of_segments + do m=1,OBC%segment(n)%num_fields + if (OBC%segment(n)%field(m)%use_IO) then + use_IO = .true. + exit + endif + enddo + if (use_IO) exit + enddo + + OBC%any_needs_IO_for_data = any_across_PEs(use_IO) +end subroutine OBC_any_IO + !> Allocate data (buffer_src, buffer_dst and dz_src) for a field at an OBC segment. subroutine allocate_segment_field_data(field, OBC, segment, US, inputdir, filename, varname, & suffix, value, turns, nz) @@ -1060,6 +1078,8 @@ subroutine allocate_segment_field_data(field, OBC, segment, US, inputdir, filena integer :: dim ! Loop index for siz/siz_check integer :: nk_dst ! k-axis size of buffer_dst + if (.not. segment%on_pe) return + isd = segment%HI%isd ; ied = segment%HI%ied ; IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB jsd = segment%HI%jsd ; jed = segment%HI%jed ; JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB nk_dst = nz @@ -1070,10 +1090,9 @@ subroutine allocate_segment_field_data(field, OBC, segment, US, inputdir, filena ! The scale factor for tracers may also be set in register_segment_tracer, and a constant input ! value is rescaled there. field%scale = scale_factor_from_name(field%name, US, segment%tr_Reg) + field%use_IO = (trim(filename) /= 'none') - if (trim(filename) /= 'none') then - field%use_IO = .true. - + if (field%use_IO) then full_filename = trim(inputdir) // trim(filename) full_varname = trim(varname) // trim(suffix) @@ -1133,8 +1152,6 @@ subroutine allocate_segment_field_data(field, OBC, segment, US, inputdir, filena init_value_dst = 0.0 else ! This data is not being read from a file. - field%use_IO = .false. - field%value = field%scale * value ! Change the sign of the specified velocities, depending on the number of quarter turns of the grid. if ( ( ((field%name == 'U') .or. (field%name == 'Uamp')) .and. & @@ -1200,7 +1217,6 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) integer :: current_pe integer, dimension(1) :: single_pelist type(external_tracers_segments_props), pointer :: obgc_segments_props_list =>NULL() - integer :: IO_needs(2) ! Sums to determine global OBC data use and update patterns. logical :: check_ts_needed ! Check if temperature and salinity are explicitly specified. integer :: idx character(len=256) :: routine_name ! Name of this subroutine @@ -1209,6 +1225,8 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) routine_name = trim(mdl) // ', initialize_segment_data' + OBC%update_OBC = .true. ! Data is time-dependent if not using user BC. + check_ts_needed = use_temperature .and. (.not. OBC%ts_needed_bug) call get_param(PF, mdl, "INPUTDIR", inputdir, default=".") @@ -1274,7 +1292,8 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) write(mesg,'("OBC segment ",I0," has an unknown input field: ",a)') n, trim(phys_inputs(m)) call MOM_error(FATAL, trim(routine_name) // ", " // trim(mesg)) endif - if (.not. segment%field(idx)%required) then + if ((.not. segment%field(idx)%required) .and. & + ((.not. (idx == F_T .or. idx == F_S)) .or. check_ts_needed)) then write(mesg,'("OBC segment ",I0," has an unnecessary field: ",a)') & n, trim(phys_inputs(m)) call MOM_error(WARNING, trim(mesg)) @@ -1304,10 +1323,10 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) enddo ! Allocate BGC tracer fields - obgc_segments_props_list => OBC%obgc_segments_props !pointer to the head node + obgc_segments_props_list => OBC%obgc_segments_props ! pointer to the head node do m = NUM_PHYS_FIELDS+1, segment%num_fields segment%field(m)%bgc_tracer = .true. - ! Query the obgc segment properties by traversing the linkedlist + ! Query the obgc segment properties by traversing the linked list call get_obgc_segments_props(obgc_segments_props_list, bgc_input, filename, varname, & segment%field(m)%resrv_lfac_in, segment%field(m)%resrv_lfac_out) ! Make sure the obgc tracer is not specified in the MOM6 param file too. @@ -1322,16 +1341,6 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) inputdir, filename, varname, suffix, 0.0, turns, GV%ke) enddo - !! - ! CODE HERE FOR OTHER OPTIONS (CLAMPED, NUDGED,..) - !! - do m=1,segment%num_fields - if (segment%field(m)%use_IO) then - OBC%update_OBC = .true. ! Data is assumed to be time-dependent if we are reading from file - OBC%needs_IO_for_data = .true. ! At least one segment is using I/O for OBC data - endif - enddo - ! write(stderr, '(A)') trim(suffix)//" segment checksum" if (OBC%debug) call chksum_OBC_segment_data(OBC%segment(n_seg), GV, US, OBC%nk_OBC_debug, n) @@ -1340,12 +1349,7 @@ subroutine initialize_segment_data(GV, US, OBC, PF, turns, use_temperature) call Set_PElist(saved_pelist) ! Determine global IO data requirement patterns. - IO_needs(1) = 0 ; if (OBC%needs_IO_for_data) IO_needs(1) = 1 - IO_needs(2) = 0 ; if (OBC%update_OBC) IO_needs(2) = 1 - call sum_across_PES(IO_needs, 2) - OBC%any_needs_IO_for_data = (IO_needs(1) > 0) - OBC%update_OBC = (IO_needs(2) > 0) - + call OBC_any_IO(OBC) end subroutine initialize_segment_data !> Determine whether a particular field is descretized at the normal-velocity faces of an open @@ -4142,7 +4146,6 @@ subroutine allocate_OBC_segment_data(OBC, segment) ! Allocate dZtot with extra values at the end to avoid segmentation faults in cases where ! it is interpolated to OBC vorticity points. allocate(segment%dZtot(IsdB:IedB,jsd-1:jed+1), source=0.0) - allocate(segment%h(IsdB:IedB,jsd:jed,OBC%ke), source=0.0) allocate(segment%SSH(IsdB:IedB,jsd:jed), source=0.0) allocate(segment%tidal_elev(IsdB:IedB,jsd:jed), source=0.0) if (segment%radiation) & @@ -4187,7 +4190,6 @@ subroutine allocate_OBC_segment_data(OBC, segment) ! Allocate dZtot with extra values at the end to avoid segmentation faults in cases where ! it is interpolated to OBC vorticity points. allocate(segment%dZtot(isd-1:ied+1,JsdB:JedB), source=0.0) - allocate(segment%h(isd:ied,JsdB:JedB,OBC%ke), source=0.0) allocate(segment%SSH(isd:ied,JsdB:JedB), source=0.0) allocate(segment%tidal_elev(isd:ied,JsdB:JedB), source=0.0) if (segment%radiation) & @@ -4236,7 +4238,6 @@ subroutine deallocate_OBC_segment_data(segment) if (allocated(segment%Htot)) deallocate(segment%Htot) if (allocated(segment%dZtot)) deallocate(segment%dZtot) - if (allocated(segment%h)) deallocate(segment%h) if (allocated(segment%SSH)) deallocate(segment%SSH) if (allocated(segment%tidal_elev)) deallocate(segment%tidal_elev) if (allocated(segment%rx_norm_rad)) deallocate(segment%rx_norm_rad) @@ -4677,7 +4678,7 @@ subroutine read_OBC_segment_data(G, GV, US, OBC, tv, h, Time) enddo ! endd segment loop end subroutine read_OBC_segment_data -!> Update the OBC values on the segments. +!> Update OBC segment velocities, gradient, SSH and the external fields %t of thickness/tracer reservoirs. subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure @@ -4692,10 +4693,10 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB integer :: is_seg, ie_seg, js_seg, je_seg ! Orientation-agnostic loop ranges integer :: i_offset_in, j_offset_in ! Indexing offset for interior cells + integer :: F_G, F_VN, F_VNAMP, F_VNPHASE, F_VT, F_VTAMP, F_VTPHASE ! Field indices real :: ramp_value ! If OBC%ramp is True, where we are on the ramp from 0 to 1, or 1 otherwise [nondim]. real :: time_delta ! Time since tidal reference date [T ~> s] real :: tidal_amp, tidal_phase ! Tidal amplitude [Z ~> m] and phase [rad] - integer :: F_G, F_VT, F_VTAMP, F_VTPHASE if (.not. associated(OBC)) return if (OBC%user_BCs_set_globally) return @@ -4710,6 +4711,12 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) if (.not. segment%on_pe) cycle ! continue to next segment if not in data domain + ! Segment indices are on q points: + ! | x | x | x | x | jsd/jed (if southern boundary) + ! |-----------|-----------|-----------|-----------| JsdB/JedB + ! IsdB isd ied IedB + ! | x | x | x | x | jsd/jed (if northern boundary) + isd = segment%HI%isd ; ied = segment%HI%ied ; IsdB = segment%HI%IsdB ; IedB = segment%HI%IedB jsd = segment%HI%jsd ; jed = segment%HI%jed ; JsdB = segment%HI%JsdB ; JedB = segment%HI%JedB i_offset_in = ied - IedB ! = 0 if East, South, North; = 1 if West @@ -4718,100 +4725,67 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) if (segment%is_E_or_W) then is_seg = IsdB ; ie_seg = is_seg js_seg = jsd ; je_seg = jed + F_VN = F_U ; F_VNAMP = F_UAMP ; F_VNPHASE = F_UPHASE + F_VT = F_V ; F_VTAMP = F_VAMP ; F_VTPHASE = F_VPHASE ; F_G = F_VX else is_seg = isd ; ie_seg = ied js_seg = JsdB ; je_seg = js_seg + F_VN = F_V ; F_VNAMP = F_VAMP ; F_VNPHASE = F_VPHASE + F_VT = F_U ; F_VTAMP = F_UAMP ; F_VTPHASE = F_UPHASE ; F_G = F_UY endif - ! Calculate auxiliary fields at staggered locations. - ! Segment indices are on q points: - ! - ! |-----------|------------|-----------|-----------| J_obc - ! Is_obc Ie_obc - ! - ! i2 has to start at Is_obc+1 and end at Ie_obc. - ! j2 is J_obc and jshift has to be +1 at both the north and south. - - ! Calculate auxiliary fields at staggered locations - segment%Htot(:,:) = 0.0 - do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg - segment%h(i,j,k) = h(i+i_offset_in,j+j_offset_in,k) - segment%Htot(i,j) = segment%Htot(i,j) + segment%h(i,j,k) - enddo ; enddo ; enddo - - ! Update segment velocities, gradient, SSH and thickness/tracer reserviors - if (segment%is_E_or_W .and. allocated(segment%field(F_U)%buffer_dst)) then + ! Update normal velocity, transport. Split by orientation for now because of G%dyCu and G%dxCv. + if (allocated(segment%field(F_VN)%buffer_dst)) then ! Update tidal normal velocity segment%tidal_vn(:,:) = 0.0 if (OBC%add_tide_constituents) then - do c=1,OBC%n_tide_constituents ; do j=jsd,jed ; do I=IsdB,IedB - tidal_amp = OBC%tide_fn(c) * segment%field(F_UAMP)%buffer_dst(I,j,c) - tidal_phase = (time_delta * OBC%tide_frequencies(c) - segment%field(F_UPHASE)%buffer_dst(I,j,c)) & + do c=1,OBC%n_tide_constituents ; do j=js_seg,je_seg ; do i=is_seg,ie_seg + tidal_amp = OBC%tide_fn(c) * segment%field(F_VNAMP)%buffer_dst(i,j,c) + tidal_phase = (time_delta * OBC%tide_frequencies(c) - segment%field(F_VNPHASE)%buffer_dst(i,j,c)) & + (OBC%tide_eq_phases(c) + OBC%tide_un(c)) - segment%tidal_vn(I,j) = segment%tidal_vn(I,j) + tidal_amp * cos(tidal_phase) + segment%tidal_vn(i,j) = segment%tidal_vn(i,j) + tidal_amp * cos(tidal_phase) enddo ; enddo ; enddo endif + segment%Htot(:,:) = 0.0 segment%normal_trans_bt(:,:) = 0.0 - do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - segment%normal_vel(I,j,k) = segment%field(F_U)%buffer_dst(I,j,k) + segment%tidal_vn(I,j) - segment%normal_trans(I,j,k) = segment%normal_vel(I,j,k) * segment%h(I,j,k) * G%dyCu(I,j) - segment%normal_trans_bt(I,j) = segment%normal_trans_bt(I,j) + segment%normal_trans(I,j,k) - enddo ; enddo ; enddo - - do j=jsd,jed ; do I=IsdB,IedB - segment%normal_vel_bt(I,j) = segment%normal_trans_bt(I,j) & - / (max(segment%Htot(I,j), 1.e-12 * GV%m_to_H) * G%dyCu(I,j)) - enddo ; enddo - - if (allocated(segment%nudged_normal_vel)) then - do k=1,nz ; do j=jsd,jed ; do I=IsdB,IedB - segment%nudged_normal_vel(I,j,k) = segment%normal_vel(I,j,k) + if (segment%is_E_or_W) then + do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%Htot(i,j) = segment%Htot(i,j) + h(i+i_offset_in,j+j_offset_in,k) + segment%normal_vel(i,j,k) = segment%field(F_VN)%buffer_dst(i,j,k) + segment%tidal_vn(i,j) + segment%normal_trans(i,j,k) = & + segment%normal_vel(i,j,k) * h(i+i_offset_in,j+j_offset_in,k) * G%dyCu(i,j) + segment%normal_trans_bt(i,j) = segment%normal_trans_bt(i,j) + segment%normal_trans(i,j,k) enddo ; enddo ; enddo - endif - endif - - if (segment%is_N_or_S .and. allocated(segment%field(F_V)%buffer_dst)) then - ! Update tidal normal velocity - segment%tidal_vn(:,:) = 0.0 - if (OBC%add_tide_constituents) then - do c=1,OBC%n_tide_constituents ; do J=JsdB,JedB ; do i=isd,ied - tidal_amp = OBC%tide_fn(c) * segment%field(F_VAMP)%buffer_dst(i,J,c) - tidal_phase = (time_delta * OBC%tide_frequencies(c) - segment%field(F_VPHASE)%buffer_dst(i,J,c)) & - + (OBC%tide_eq_phases(c) + OBC%tide_un(c)) - segment%tidal_vn(i,J) = segment%tidal_vn(i,J) + tidal_amp * cos(tidal_phase) + do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%normal_vel_bt(i,j) = segment%normal_trans_bt(i,j) & + / (max(segment%Htot(i,j), 1.e-12 * GV%m_to_H) * G%dyCu(i,j)) + enddo ; enddo + else + do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%Htot(i,j) = segment%Htot(i,j) + h(i+i_offset_in,j+j_offset_in,k) + segment%normal_vel(i,j,k) = segment%field(F_VN)%buffer_dst(i,j,k) + segment%tidal_vn(i,j) + segment%normal_trans(i,j,k) = & + segment%normal_vel(i,j,k) * h(i+i_offset_in,j+j_offset_in,k) * G%dxCv(i,j) + segment%normal_trans_bt(i,j) = segment%normal_trans_bt(i,j) + segment%normal_trans(i,j,k) enddo ; enddo ; enddo + do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%normal_vel_bt(i,j) = segment%normal_trans_bt(i,j) & + / (max(segment%Htot(i,j), 1.e-12 * GV%m_to_H) * G%dxCv(i,j)) + enddo ; enddo endif - segment%normal_trans_bt(:,:) = 0.0 - do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied - segment%normal_vel(i,J,k) = segment%field(F_V)%buffer_dst(i,J,k) + segment%tidal_vn(i,J) - segment%normal_trans(i,J,k) = segment%normal_vel(i,J,k) * segment%h(i,J,k) * G%dxCv(i,J) - segment%normal_trans_bt(i,J) = segment%normal_trans_bt(i,J) + segment%normal_trans(i,J,k) - enddo ; enddo ; enddo - - do J=JsdB,JedB ; do i=isd,ied - segment%normal_vel_bt(i,J) = segment%normal_trans_bt(i,J) & - / (max(segment%Htot(i,J), 1.e-12 * GV%m_to_H) * G%dxCv(i,J)) - enddo ; enddo - if (allocated(segment%nudged_normal_vel)) then - do k=1,nz ; do J=JsdB,JedB ; do i=isd,ied - segment%nudged_normal_vel(i,J,k) = segment%normal_vel(i,J,k) + do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg + segment%nudged_normal_vel(i,j,k) = segment%normal_vel(i,j,k) enddo ; enddo ; enddo endif endif - if (((segment%is_E_or_W .and. allocated(segment%field(F_V)%buffer_dst)) .or. & - (segment%is_N_or_S .and. allocated(segment%field(F_U)%buffer_dst))) .and. & - allocated(segment%tangential_vel)) then - if (segment%is_E_or_W) then - F_VT = F_V ; F_VTAMP = F_VAMP ; F_VTPHASE = F_VPHASE - else - F_VT = F_U ; F_VTAMP = F_UAMP ; F_VTPHASE = F_UPHASE - endif - segment%tidal_vt(:,:) = 0.0 + ! Update tangential velocity + if (allocated(segment%tangential_vel) .and. allocated(segment%field(F_VT)%buffer_dst)) then ! Update tidal tangential velocity + segment%tidal_vt(:,:) = 0.0 if (OBC%add_tide_constituents) then do c=1,OBC%n_tide_constituents ; do J=JsdB,JedB ; do I=IsdB,IedB tidal_amp = OBC%tide_fn(c) * segment%field(F_VTAMP)%buffer_dst(I,J,c) @@ -4832,12 +4806,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) endif endif - if (segment%is_E_or_W) then - F_G = F_VX - else - F_G = F_UY - endif - + ! Update tangential gradient dvdx and dudy if (allocated(segment%tangential_grad) .and. allocated(segment%field(F_G)%buffer_dst)) then do k=1,nz ; do J=JsdB,JedB ; do I=IsdB,IedB segment%tangential_grad(I,J,k) = segment%field(F_G)%buffer_dst(I,J,k) @@ -4850,6 +4819,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) endif endif + ! Update SSH if (allocated(segment%field(F_Z)%buffer_dst)) then ! Update tidal SSH segment%tidal_elev(:,:) = 0.0 @@ -4868,13 +4838,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, h, Time) enddo ; enddo endif - ! Set the thickness reservoir data. + ! Update thickness registry if (OBC%thickness_x_reservoirs_used .or. OBC%thickness_y_reservoirs_used) then do k=1,nz ; do j=js_seg,je_seg ; do i=is_seg,ie_seg segment%h_Reg%h(i,j,k) = h(i+i_offset_in,j+j_offset_in,k) enddo ; enddo ; enddo endif + ! Update tracer registry do m = NUM_PHYS_FIELDS-1, segment%num_fields ! F_T = NUM_PHYS_FIELDS-1 and F_S = NUM_PHYS_FIELDS if (.not. allocated(segment%field(m)%buffer_dst) .or. & (segment%field(m)%bgc_tracer .and. (.not. OBC%update_OBC_seg_data))) then @@ -6650,7 +6621,6 @@ subroutine rotate_OBC_config(OBC_in, G_in, OBC, G, turns) ! These are set by initialize_segment_data OBC%brushcutter_mode = OBC_in%brushcutter_mode OBC%update_OBC = OBC_in%update_OBC - OBC%needs_IO_for_data = OBC_in%needs_IO_for_data OBC%any_needs_IO_for_data = OBC_in%any_needs_IO_for_data OBC%update_OBC_seg_data = OBC_in%update_OBC_seg_data @@ -6988,7 +6958,6 @@ subroutine write_OBC_info(OBC, G, GV, US) if (OBC%user_BCs_set_globally) call MOM_mesg("user_BCs_set_globally", verb=1) if (OBC%update_OBC) call MOM_mesg("update_OBC", verb=1) if (OBC%update_OBC_seg_data) call MOM_mesg("update_OBC_seg_data", verb=1) - if (OBC%needs_IO_for_data) call MOM_mesg("needs_IO_for_data", verb=1) if (OBC%any_needs_IO_for_data) call MOM_mesg("any_needs_IO_for_data", verb=1) if (OBC%zero_biharmonic) call MOM_mesg("zero_biharmonic", verb=1) if (OBC%brushcutter_mode) call MOM_mesg("brushcutter_mode", verb=1) @@ -7140,7 +7109,6 @@ subroutine chksum_OBC_segment_data(segment, GV, US, nk, nseg_out) if (allocated(segment%Htot)) call write_2d_array_vals("Htot"//trim(sn), segment%Htot, dir, nk, unscale=GV%H_to_mks) if (allocated(segment%dZtot)) call write_2d_array_vals("dZtot"//trim(sn), segment%dZtot, dir, nk, unscale=US%Z_to_m) if (allocated(segment%SSH)) call write_2d_array_vals("SSH"//trim(sn), segment%SSH, dir, nk, unscale=US%Z_to_m) - if (allocated(segment%h)) call write_3d_array_vals("h"//trim(sn), segment%h, dir, nk, unscale=GV%H_to_mks) if (allocated(segment%normal_vel)) & call write_3d_array_vals("normal_vel"//trim(sn), segment%normal_vel, dir, nk, unscale=norm*US%L_T_to_m_s) if (allocated(segment%normal_vel_bt)) & diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 86c15681cd..86d6015598 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -16,6 +16,7 @@ module MOM_diag_mediator use MOM_diag_manager_infra, only : send_data_infra, MOM_diag_field_add_attribute, EAST, NORTH use MOM_diag_manager_infra, only : register_diag_field_infra, register_static_field_infra use MOM_diag_manager_infra, only : get_MOM_diag_field_id, DIAG_FIELD_NOT_FOUND +use MOM_diag_manager_infra, only : diag_send_complete_infra use MOM_diag_remap, only : diag_remap_ctrl, diag_remap_update, diag_remap_calc_hmask use MOM_diag_remap, only : diag_remap_init, diag_remap_end, diag_remap_do_remap use MOM_diag_remap, only : vertically_reintegrate_diag_field, vertically_interpolate_diag_field @@ -46,7 +47,7 @@ module MOM_diag_mediator public set_axes_info, post_data, register_diag_field, time_type public post_data_3d_by_column, post_data_3d_final public post_product_u, post_product_sum_u, post_product_v, post_product_sum_v -public set_masks_for_axes +public set_masks_for_axes, MOM_diag_send_complete ! post_data_1d_k is a deprecated interface that can be replaced by a call to post_data, but ! it is being retained for backward compatibility to older versions of the ocean_BGC code. public post_data_1d_k @@ -73,6 +74,11 @@ module MOM_diag_mediator module procedure post_data_3d, post_data_2d, post_data_1d_k, post_data_0d end interface post_data +!> Registers a non-array scalar diagnostic, returning an integer handle +interface register_scalar_field + module procedure register_scalar_field_CS, register_scalar_field_axes +end interface register_scalar_field + !> Down sample a field interface downsample_field module procedure downsample_field_2d, downsample_field_3d @@ -260,7 +266,7 @@ module MOM_diag_mediator integer :: ie !< The end i-index of cell centers within the computational domain integer :: js !< The start j-index of cell centers within the computational domain integer :: je !< The end j-index of cell centers within the computational domain - ! These give the memory-domain sizes, and can be start at any value on each PE. + ! These give the memory-domain sizes, and can start at any value on each PE. integer :: isd !< The start i-index of cell centers within the data domain integer :: ied !< The end i-index of cell centers within the data domain integer :: jsd !< The start j-index of cell centers within the data domain @@ -329,7 +335,7 @@ module MOM_diag_mediator real, dimension(:,:,:), pointer :: S => null() !< The salinities needed for remapping [S ~> ppt] type(EOS_type), pointer :: eqn_of_state => null() !< The equation of state type type(thermo_var_ptrs), pointer :: tv => null() !< A structure with thermodynamic variables that are - !! are used to convert thicknesses to vertical extents + !! used to convert thicknesses to vertical extents type(ocean_grid_type), pointer :: G => null() !< The ocean grid type type(verticalGrid_type), pointer :: GV => null() !< The model's vertical ocean grid type(unit_scale_type), pointer :: US => null() !< A dimensional unit scaling type @@ -380,62 +386,56 @@ subroutine set_axes_info(G, GV, US, param_file, diag_cs, set_vertical) real, allocatable, dimension(:) :: IaxB, iax ! Index-based integer and half-integer i-axis labels [nondim] real, allocatable, dimension(:) :: JaxB, jax ! Index-based integer and half-integer j-axis labels [nondim] - set_vert = .true. ; if (present(set_vertical)) set_vert = set_vertical - if (diag_cs%index_space_axes) then allocate(IaxB(G%IsgB:G%IegB)) - do i=G%IsgB, G%IegB - Iaxb(i)=real(i) + do I=G%IsgB,G%IegB + Iaxb(I) = real(I) enddo allocate(iax(G%isg:G%ieg)) - do i=G%isg, G%ieg - iax(i)=real(i)-0.5 + do i=G%isg,G%ieg + iax(i) = real(i)-0.5 enddo allocate(JaxB(G%JsgB:G%JegB)) - do j=G%JsgB, G%JegB - JaxB(j)=real(j) + do J=G%JsgB,G%JegB + JaxB(J) = real(J) enddo allocate(jax(G%jsg:G%jeg)) - do j=G%jsg, G%jeg - jax(j)=real(j)-0.5 + do j=G%jsg,G%jeg + jax(j) = real(j)-0.5 enddo endif ! Horizontal axes for the native grids - if (G%symmetric) then - if (diag_cs%index_space_axes) then - id_xq = diag_axis_init('iq', IaxB(G%isgB:G%iegB), 'none', 'x', & - 'q point grid-space longitude', G%Domain, position=EAST) - id_yq = diag_axis_init('jq', JaxB(G%jsgB:G%jegB), 'none', 'y', & - 'q point grid space latitude', G%Domain, position=NORTH) + if (diag_cs%index_space_axes) then + if (G%symmetric) then + id_xq = diag_axis_init('Iq', IaxB(G%IsgB:G%IegB), 'none', 'x', & + 'Boundary (q) point grid-space longitude', G%Domain, position=EAST) + id_yq = diag_axis_init('Jq', JaxB(G%JsgB:G%JegB), 'none', 'y', & + 'Boundary (q) point grid-space latitude', G%Domain, position=NORTH) else - id_xq = diag_axis_init('xq', G%gridLonB(G%isgB:G%iegB), G%x_axis_units, 'x', & - 'q point nominal longitude', G%Domain, position=EAST) - id_yq = diag_axis_init('yq', G%gridLatB(G%jsgB:G%jegB), G%y_axis_units, 'y', & - 'q point nominal latitude', G%Domain, position=NORTH) - endif - else - if (diag_cs%index_space_axes) then id_xq = diag_axis_init('Iq', IaxB(G%isg:G%ieg), 'none', 'x', & - 'q point grid-space longitude', G%Domain, position=EAST) + 'Boundary (q) point grid-space longitude', G%Domain, position=EAST) id_yq = diag_axis_init('Jq', JaxB(G%jsg:G%jeg), 'none', 'y', & - 'q point grid space latitude', G%Domain, position=NORTH) + 'Boundary (q) point grid-space latitude', G%Domain, position=NORTH) + endif + id_xh = diag_axis_init('ih', iax(G%isg:G%ieg), 'none', 'x', & + 'Tracer (h) point grid-space longitude', G%Domain) + id_yh = diag_axis_init('jh', jax(G%jsg:G%jeg), 'none', 'y', & + 'Tracer (h) point grid-space latitude', G%Domain) + else + if (G%symmetric) then + id_xq = diag_axis_init('xq', G%gridLonB(G%IsgB:G%IegB), G%x_axis_units, 'x', & + 'q point nominal longitude', G%Domain, position=EAST) + id_yq = diag_axis_init('yq', G%gridLatB(G%JsgB:G%JegB), G%y_axis_units, 'y', & + 'q point nominal latitude', G%Domain, position=NORTH) else id_xq = diag_axis_init('xq', G%gridLonB(G%isg:G%ieg), G%x_axis_units, 'x', & 'q point nominal longitude', G%Domain, position=EAST) id_yq = diag_axis_init('yq', G%gridLatB(G%jsg:G%jeg), G%y_axis_units, 'y', & 'q point nominal latitude', G%Domain, position=NORTH) endif - endif - - if (diag_cs%index_space_axes) then - id_xh = diag_axis_init('ih', iax(G%isg:G%ieg), 'none', 'x', & - 'h point grid-space longitude', G%Domain) - id_yh = diag_axis_init('jh', jax(G%jsg:G%jeg), 'none', 'y', & - 'h point grid space latitude', G%Domain) - else id_xh = diag_axis_init('xh', G%gridLonT(G%isg:G%ieg), G%x_axis_units, 'x', & 'h point nominal longitude', G%Domain) id_yh = diag_axis_init('yh', G%gridLatT(G%jsg:G%jeg), G%y_axis_units, 'y', & @@ -1439,7 +1439,7 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) real, dimension(:,:), pointer :: locfield ! The field being offered in arbitrary unscaled units [a] real, dimension(:,:), pointer :: locmask ! A pointer to the data mask to use [nondim] logical :: used ! The return value of send_data is not used for anything. - logical :: is_stat + logical :: is_stat, not_static integer :: cszi, cszj, dszi, dszj integer :: isv, iev, jsv, jev, i, j, isv_o, jsv_o real, dimension(:,:), allocatable, target :: locfield_dsamp ! A downsampled version of locfield [a] @@ -1453,6 +1453,7 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) locfield => NULL() locmask => NULL() is_stat = .false. ; if (present(is_static)) is_stat = is_static + not_static = .not. is_stat ! Determine the proper array indices, noting that because of the (:,:) ! declaration of field, symmetric arrays are using a SW-grid indexing, @@ -1502,12 +1503,14 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) if (present(mask)) then locmask => mask - elseif (.NOT. is_stat .and. associated(diag%axes)) then + elseif (not_static .and. associated(diag%axes)) then + ! If we were to decide to allow masking of static diagnostics, we could do so by changing the line above to + ! elseif (associated(diag%axes) .and. (diag_CS%mask_static_diags .or. not_static)) then if (associated(diag%axes%mask2d)) locmask => diag%axes%mask2d endif dl = 1 - if (.NOT. is_stat .and. associated(diag%axes)) & + if (not_static .and. associated(diag%axes)) & dl = diag%axes%downsample_level ! Static field downsampling is not supported yet. ! Downsample the diag field and mask as appropriate. if (dl > 1) then @@ -1523,6 +1526,8 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) locmask => diag%axes%dsamp(dl)%mask2d endif endif + if (associated(locmask)) call assert(size(locfield) == size(locmask), & + 'post_data_2d_low: mask size mismatch: '//trim(diag%debug_str)) if (diag_cs%diag_as_chksum) then ! Append timestep to mesg @@ -1547,22 +1552,15 @@ subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) endif else if (is_stat) then - if (present(mask)) then - call assert(size(locfield) == size(locmask), & - 'post_data_2d_low is_stat: mask size mismatch: '//trim(diag%debug_str)) + if (associated(locmask)) then used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=locmask) - !elseif (associated(diag%axes%mask2d)) then - ! used = send_data(diag%fms_diag_id, locfield, & - ! is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=diag%axes%mask2d) else used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev) endif elseif (diag_cs%ave_enabled) then if (associated(locmask)) then - call assert(size(locfield) == size(locmask), & - 'post_data_2d_low: mask size mismatch: '//trim(diag%debug_str)) used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, & time=diag_cs%time_end, weight=diag_cs%time_int, rmask=locmask) @@ -1757,7 +1755,7 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) character(len=300) :: mesg logical :: used ! The return value of send_data is not used for anything. logical :: staggered_in_x, staggered_in_y - logical :: is_stat + logical :: is_stat, not_static integer :: cszi, cszj, dszi, dszj integer :: isv, iev, jsv, jev, ks, ke, i, j, k, isv_c, jsv_c, isv_o, jsv_o real, dimension(:,:,:), allocatable, target :: locfield_dsamp ! A downsampled version of locfield [a] @@ -1771,6 +1769,7 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) locfield => NULL() locmask => NULL() is_stat = .false. ; if (present(is_static)) is_stat = is_static + not_static = .not. is_stat ! Determine the proper array indices, noting that because of the (:,:) ! declaration of field, symmetric arrays are using a SW-grid indexing, @@ -1838,12 +1837,14 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) if (present(mask)) then locmask => mask - elseif (associated(diag%axes%mask3d)) then - locmask => diag%axes%mask3d + elseif (associated(diag%axes) .and. (not_static)) then + ! If we were to decide to allow masking of static diagnostics, we could do so by changing the line above to + ! elseif (associated(diag%axes) .and. (diag_CS%mask_static_diags .or. not_static)) then + if (associated(diag%axes%mask3d)) locmask => diag%axes%mask3d endif dl = 1 - if (.NOT. is_stat) & + if (not_static .and. associated(diag%axes)) & dl = diag%axes%downsample_level ! Static field downsampling is not supported yet. ! Downsample the diag field and mask as appropriate. if (dl > 1) then @@ -1859,6 +1860,8 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) locmask => diag%axes%dsamp(dl)%mask3d endif endif + if (associated(locmask)) call assert(size(locfield) == size(locmask), & + 'post_data_3d_low: mask size mismatch: '//trim(diag%debug_str)) if (diag%fms_diag_id>0) then if (diag_cs%diag_as_chksum) then @@ -1884,22 +1887,15 @@ subroutine post_data_3d_low(diag, field, diag_cs, is_static, mask) endif else if (is_stat) then - if (present(mask)) then - call assert(size(locfield) == size(locmask), & - 'post_data_3d_low is_stat: mask size mismatch: '//diag%debug_str) + if (associated(locmask)) then used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=locmask) - !elseif (associated(diag%axes%mask2d)) then - ! used = send_data(diag%fms_diag_id, locfield, & - ! is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, rmask=diag%axes%mask2d) else used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev) endif elseif (diag_cs%ave_enabled) then if (associated(locmask)) then - call assert(size(locfield) == size(locmask), & - 'post_data_3d_low: mask size mismatch: '//diag%debug_str) used = send_data_infra(diag%fms_diag_id, locfield, & is_in=isv, ie_in=iev, js_in=jsv, je_in=jev, & time=diag_cs%time_end, weight=diag_cs%time_int, rmask=locmask) @@ -2263,7 +2259,7 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time character(len=256) :: msg, cm_string character(len=256) :: new_module_name character(len=480) :: module_list, var_list - character(len=16) :: dimensions + character(len=24) :: dimensions integer :: num_modnm, num_varnm logical :: active @@ -2287,6 +2283,14 @@ integer function register_diag_field(module_name, field_name, axes_in, init_time axes => diag_cs%axesCui elseif (axes_in%id == diag_cs%axesCvi%id) then axes => diag_cs%axesCvi + elseif (axes_in%id == diag_cs%axesT1%id) then + axes => diag_cs%axesT1 + elseif (axes_in%id == diag_cs%axesB1%id) then + axes => diag_cs%axesB1 + elseif (axes_in%id == diag_cs%axesCu1%id) then + axes => diag_cs%axesCu1 + elseif (axes_in%id == diag_cs%axesCv1%id) then + axes => diag_cs%axesCv1 else allocate(axes) axes = axes_in @@ -2614,7 +2618,7 @@ logical function register_diag_field_expand_cmor(dm_id, module_name, field_name, posted_cmor_long_name = "not provided" ! ! If attributes are present for MOM variable names, use them first for the register_diag_field - ! call for CMOR verison of the variable + ! call for CMOR version of the variable if (present(units)) posted_cmor_units = units if (present(standard_name)) posted_cmor_standard_name = standard_name if (present(long_name)) posted_cmor_long_name = long_name @@ -2971,11 +2975,51 @@ subroutine attach_cell_methods(id, axes, ostring, cell_methods, & end subroutine attach_cell_methods -!> Registers a scalar diagnostic, returning an integer handle -function register_scalar_field(module_name, field_name, init_time, diag_cs, & +!> Registers a non-array scalar diagnostic, returning an integer handle +function register_scalar_field_axes(module_name, field_name, axes, init_time, & + long_name, units, missing_value, range, standard_name, & + do_not_log, err_msg, interp_method, cmor_field_name, & + cmor_long_name, cmor_units, cmor_standard_name, conversion) result (register_scalar_field) + integer :: register_scalar_field !< An integer handle for a diagnostic array. + character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" + !! or "ice_shelf_model" + character(len=*), intent(in) :: field_name !< Name of the diagnostic field + type(axes_grp), target, intent(in) :: axes !< Container with up to 3 integer handles that + !! indicates axes for this field + type(time_type), intent(in) :: init_time !< Time at which a field is first available? + character(len=*), optional, intent(in) :: long_name !< Long name of a field. + character(len=*), optional, intent(in) :: units !< Units of a field. + character(len=*), optional, intent(in) :: standard_name !< Standardized name associated with a field + real, optional, intent(in) :: missing_value !< A value that indicates missing values in + !! output files, in unscaled arbitrary units [a] + real, optional, intent(in) :: range(2) !< Valid range of a variable (not used in MOM?) + !! in arbitrary units [a] + logical, optional, intent(in) :: do_not_log !< If true, do not log something (not used in MOM?) + character(len=*), optional, intent(out):: err_msg !< String into which an error message might be + !! placed (not used in MOM?) + character(len=*), optional, intent(in) :: interp_method !< If 'none' indicates the field should not + !! be interpolated as a scalar + character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name of a field + character(len=*), optional, intent(in) :: cmor_long_name !< CMOR long name of a field + character(len=*), optional, intent(in) :: cmor_units !< CMOR units of a field + character(len=*), optional, intent(in) :: cmor_standard_name !< CMOR standardized name associated with a field + real, optional, intent(in) :: conversion !< A value to multiply data by before writing to files, + !! often including factors to undo internal scaling and + !! in units of [a A-1 ~> 1] + + register_scalar_field = register_scalar_field_CS(module_name, field_name, init_time, axes%diag_cs, & long_name, units, missing_value, range, standard_name, & do_not_log, err_msg, interp_method, cmor_field_name, & cmor_long_name, cmor_units, cmor_standard_name, conversion) + +end function register_scalar_field_axes + + +!> Registers a scalar diagnostic, returning an integer handle +function register_scalar_field_CS(module_name, field_name, init_time, diag_cs, & + long_name, units, missing_value, range, standard_name, & + do_not_log, err_msg, interp_method, cmor_field_name, & + cmor_long_name, cmor_units, cmor_standard_name, conversion) result (register_scalar_field) integer :: register_scalar_field !< An integer handle for a diagnostic array. character(len=*), intent(in) :: module_name !< Name of this module, usually "ocean_model" !! or "ice_shelf_model" @@ -3042,7 +3086,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & posted_cmor_long_name = "not provided" ! If attributes are present for MOM variable names, use them as defaults for the - ! register_diag_field_infra call for CMOR verison of the variable + ! register_diag_field_infra call for CMOR version of the variable if (present(units)) posted_cmor_units = units if (present(standard_name)) posted_cmor_standard_name = standard_name if (present(long_name)) posted_cmor_long_name = long_name @@ -3084,7 +3128,7 @@ function register_scalar_field(module_name, field_name, init_time, diag_cs, & register_scalar_field = dm_id -end function register_scalar_field +end function register_scalar_field_CS !> Registers a static diagnostic, returning an integer handle function register_static_field(module_name, field_name, axes, & @@ -3130,7 +3174,7 @@ function register_static_field(module_name, field_name, axes, & integer :: dm_id, fms_id character(len=256) :: posted_cmor_units, posted_cmor_standard_name, posted_cmor_long_name character(len=9) :: axis_name - character(len=16) :: dimensions + character(len=24) :: dimensions MOM_missing_value = axes%diag_cs%missing_value if (present(missing_value)) MOM_missing_value = missing_value @@ -3186,7 +3230,7 @@ function register_static_field(module_name, field_name, axes, & posted_cmor_long_name = "not provided" ! If attributes are present for MOM variable names, use them first for the register_static_field - ! call for CMOR verison of the variable + ! call for CMOR version of the variable if (present(units)) posted_cmor_units = units if (present(standard_name)) posted_cmor_standard_name = standard_name if (present(long_name)) posted_cmor_long_name = long_name @@ -3587,7 +3631,7 @@ end subroutine diag_mediator_init subroutine diag_set_state_ptrs(h, tv, diag_cs) real, dimension(:,:,:), target, intent(in ) :: h !< the model thickness array [H ~> m or kg m-2] type(thermo_var_ptrs), target, intent(in ) :: tv !< A structure with thermodynamic variables that are - !! are used to convert thicknesses to vertical extents + !! used to convert thicknesses to vertical extents type(diag_ctrl), intent(inout) :: diag_cs !< diag mediator control structure ! Keep pointers to h, T, S needed for the diagnostic remapping @@ -4804,4 +4848,9 @@ logical function found_in_diagtable(diag, varName) end function found_in_diagtable +!> Finishes the diag manager reduction methods as needed for the time_step +subroutine MOM_diag_send_complete() + call diag_send_complete_infra() +end subroutine MOM_diag_send_complete + end module MOM_diag_mediator diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index a1fce2feed..56f01be873 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -142,6 +142,10 @@ module MOM_ice_shelf_dynamics !! 4 quadrature points surrounding the cell vertices [L-1 ~> m-1]. real, pointer, dimension(:,:,:) :: PhiC => NULL() !< The gradients of bilinear basis elements at 1 cell-centered !! quadrature point per cell [L-1 ~> m-1]. + real, pointer, dimension(:,:,:) :: Jac => NULL() !< Jacobian determinant |J_q| = a_q*d_q of the element + !! mapping at each of the 4 Gaussian quadrature points [L2 ~> m2]. + !! Equal to G%areaT only for rectangular elements; differs when + !! opposite cell edges have unequal lengths (non-rectangular quads). real, pointer, dimension(:,:,:,:,:,:) :: Phisub => NULL() !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] integer :: OD_rt_counter = 0 !< A counter of the number of contributions to OD_rt. @@ -708,9 +712,10 @@ subroutine initialize_ice_shelf_dyn(param_file, Time, ISS, CS, G, US, diag, new_ endif allocate(CS%Phi(1:8,1:4,isd:ied,jsd:jed), source=0.0) + allocate(CS%Jac(1:4,isd:ied,jsd:jed), source=0.0) do j=G%jsd,G%jed ; do i=G%isd,G%ied - call bilinear_shape_fn_grid(G, i, j, CS%Phi(:,:,i,j)) - enddo ; enddo + call bilinear_shape_fn_grid(G, i, j, CS%Phi(:,:,i,j), CS%Jac(:,i,j)) + enddo; enddo if (CS%GL_regularize) then allocate(CS%Phisub(2,2,CS%n_sub_regularize,CS%n_sub_regularize,2,2), source=0.0) @@ -2744,6 +2749,7 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, real :: ux, uy, vx, vy ! Components of velocity shears or divergence [T-1 ~> s-1] real :: uq, vq ! Interpolated velocities [L T-1 ~> m s-1] real :: strx_n, stry_n, strsh_n, dstrain_n, inner_dot_n ! Newton correction variables [T-1 ~> s-1], [T-2 ~> s-2] + real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim] integer :: iq, jq, iphi, jphi, i, j, ilq, jlq, Itgt, Jtgt, qp, qpv logical :: visc_qp4 logical :: use_newton ! Whether to apply Newton tangent stiffness corrections @@ -2824,22 +2830,26 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, inner_dot_n = CS%newton_umid(i,j)*uq + CS%newton_vmid(i,j)*vq endif + ! Ratio |J_q|/areaT corrects the uniform-area weight baked into ice_visc for + ! non-rectangular elements where opposite cell edges have unequal lengths. + jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) + do jphi=1,2 ; Jtgt = J-2+jphi ; do iphi=1,2 ; Itgt = I-2+iphi - if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & + if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = jac_wt * ice_visc(i,j,qpv) * & (((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & ((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) - if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = ice_visc(i,j,qpv) * & + if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = jac_wt * ice_visc(i,j,qpv) * & (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & ((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) ! Newton tangent stiffness correction: add (dη/dε_e^2) * (g·δε) * (g·φ_m) term if (do_newton_visc) then if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - CS%newton_visc_factor(i,j,qpv) * dstrain_n * & + jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_n * & ((2.*strx_n + stry_n) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & strsh_n * 0.5 * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - CS%newton_visc_factor(i,j,qpv) * dstrain_n * & + jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_n * & (strsh_n * 0.5 * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & (2.*stry_n + strx_n) * Phi(2*(2*(jphi-1)+iphi),qp,i,j)) endif @@ -2848,15 +2858,15 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, ilq = 1 ; if (iq == iphi) ilq = 2 jlq = 1 ; if (jq == jphi) jlq = 2 if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - ((basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))) + (jac_wt * (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq))) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - ((basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))) + (jac_wt * (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq))) ! Newton basal drag tangent stiffness: (m-1)*basal_trac/|u|^2 * u_i * (u . delta_u) contribution if (use_newton) then if (umask(Itgt,Jtgt) == 1) uret_qp(iphi,jphi,qp) = uret_qp(iphi,jphi,qp) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) + jac_wt * CS%newton_drag_coef(i,j) * CS%newton_umid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) if (vmask(Itgt,Jtgt) == 1) vret_qp(iphi,jphi,qp) = vret_qp(iphi,jphi,qp) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) + jac_wt * CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j) * inner_dot_n * (xquad(ilq) * xquad(jlq)) endif endif enddo ; enddo @@ -2883,7 +2893,8 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, Hcell(:,:) = H_node(I-1:I,J-1:J) call CG_action_subgrid_basal(Phisub, Hcell, Ucell, Vcell, & - bathyT(i,j), dens_ratio, Usub, Vsub) + bathyT(i,j), dens_ratio, Usub, Vsub, & + G%dxCv(i,j-1), G%dxCv(i,j), G%dyCu(i-1,j), G%dyCu(i,j), G%IareaT(i,j)) if (umask(I-1,J-1) == 1) uret_b(I-1,J-1,4) = uret_b(I-1,J-1,4) + (Usub(1,1) * basal_trac(i,j)) if (umask(I-1,J ) == 1) uret_b(I-1,J ,2) = uret_b(I-1,J ,2) + (Usub(1,2) * basal_trac(i,j)) @@ -2928,7 +2939,8 @@ subroutine CG_action(CS, uret, vret, u_shlf, v_shlf, Phi, Phisub, umask, vmask, end subroutine CG_action -subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, Vcontr) +subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, Vcontr, & + dxCv_S, dxCv_N, dyCu_W, dyCu_E, IareaT) real, dimension(:,:,:,:,:,:), & intent(in) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] @@ -2943,18 +2955,25 @@ subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, !! is grounded, or 0 where it is floating [L T-1 ~> m s-1]. real, dimension(2,2), intent(out) :: Vcontr !< The areal average of v-velocities where the ice shelf !! is grounded, or 0 where it is floating [L T-1 ~> m s-1]. + real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] + real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] + real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] + real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] + real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: Ucontr_sub, Vcontr_sub ! The contributions to Ucontr and Vcontr !! at each sub-cell real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: uloc_arr !The local sub-cell u-velocity [L T-1 ~> m s-1] real, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: vloc_arr !The local sub-cell v-velocity [L T-1 ~> m s-1] real, dimension(2,2) :: Ucontr_q, Vcontr_q !Contributions to a node from each quadrature point in a sub-grid cell - real :: subarea ! The fractional sub-cell area [nondim] - real :: hloc ! The local sub-cell ice thickness [Z ~> m] + real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] + real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] + real :: subarea ! The fractional sub-cell area in reference space [nondim] + real :: hloc ! The local sub-cell ice thickness [Z ~> m] integer :: nsub, i, j, qx, qy, m, n nsub = size(Phisub,3) - subarea = 1.0 / (nsub**2) + subarea = 1.0 / real(nsub)**2 uloc_arr(:,:,:,:) = 0.0 ; vloc_arr(:,:,:,:)=0.0 @@ -2971,15 +2990,25 @@ subroutine CG_action_subgrid_basal(Phisub, H, U, V, bathyT, dens_ratio, Ucontr, do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub do qy=1,2 ; do qx=1,2 + ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values + ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation + ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. + ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. + a = dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1)) + & ! (1-y) * dxCv_S + dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2)) ! + y * dxCv_N + d = dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2)) + & ! (1-x) * dyCu_W + dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2)) ! + x * dyCu_E + jac_sub_wt = 0.25 * subarea * (a * d) * IareaT + !calculate quadrature point contributions for the sub-cell, to each node - Ucontr_q(qx,qy) = Phisub(qx,qy,i,j,m,n) * uloc_arr(qx,qy,i,j) - Vcontr_q(qx,qy) = Phisub(qx,qy,i,j,m,n) * vloc_arr(qx,qy,i,j) - enddo ; enddo + Ucontr_q(qx,qy) = jac_sub_wt * Phisub(qx,qy,i,j,m,n) * uloc_arr(qx,qy,i,j) + Vcontr_q(qx,qy) = jac_sub_wt * Phisub(qx,qy,i,j,m,n) * vloc_arr(qx,qy,i,j) + enddo; enddo !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell - Ucontr_sub(i,j,m,n) = (subarea * 0.25) * ((Ucontr_q(1,1) + Ucontr_q(2,2)) + (Ucontr_q(1,2)+Ucontr_q(2,1))) - Vcontr_sub(i,j,m,n) = (subarea * 0.25) * ((Vcontr_q(1,1) + Vcontr_q(2,2)) + (Vcontr_q(1,2)+Vcontr_q(2,1))) - enddo ; enddo ; enddo ; enddo + Ucontr_sub(i,j,m,n) = (Ucontr_q(1,1) + Ucontr_q(2,2)) + (Ucontr_q(1,2)+Ucontr_q(2,1)) + Vcontr_sub(i,j,m,n) = (Vcontr_q(1,1) + Vcontr_q(2,2)) + (Vcontr_q(1,2)+Vcontr_q(2,1)) + enddo; enddo ; enddo ; enddo !sum up the sub-cell contributions to each node do n=1,2 ; do m=1,2 @@ -3070,6 +3099,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, real :: ux, uy, vx, vy ! Interpolated weight gradients [L-1 ~> m-1] real :: uq, vq + real :: jac_wt ! Per-quadrature-point metric correction |J_q|/areaT [nondim] real :: strx_n, stry_n, strsh_n ! Newton viscosity strain rates [T-1 ~> s-1] real :: dstrain_diag_u, dstrain_diag_v ! Newton viscosity diagonal correction factors [T-1 L-1 ~> s-1 m-1] real :: phi_m_sq ! Squared basis function value at quadrature point [nondim] @@ -3109,6 +3139,10 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, qp = 2*(jq-1)+iq !current quad point if (visc_qp4) qpv = qp !current quad point for viscosity + ! Ratio |J_q|/areaT corrects the uniform-area weight baked into ice_visc for + ! non-rectangular elements where opposite cell edges have unequal lengths. + jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) + ! Pre-compute Newton strain data for this QP (for viscosity diagonal correction) if (do_newton_visc) then strx_n = CS%newton_str_ux(i,j,qpv) @@ -3129,7 +3163,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, vx = 0. vy = 0. - u_diag_qp(iphi,jphi,qp) = & + u_diag_qp(iphi,jphi,qp) = jac_wt * & ice_visc(i,j,qpv) * (((4*ux+2*vy) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & ((uy+vx) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) @@ -3139,17 +3173,17 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, dstrain_diag_u = (2.*strx_n + stry_n) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & strsh_n * 0.5 * Phi(2*(2*(jphi-1)+iphi),qp,i,j) u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & - CS%newton_visc_factor(i,j,qpv) * dstrain_diag_u**2 + jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_diag_u**2 endif if (float_cond(i,j) == 0) then uq = xquad(ilq) * xquad(jlq) u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & - (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) + jac_wt * (basal_trac(i,j) * uq) * (xquad(ilq) * xquad(jlq)) ! Newton basal drag diagonal correction: newton_drag_coef * (umid_i)^2 * phi_m^2 if (CS%doing_newton) then u_diag_qp(iphi,jphi,qp) = u_diag_qp(iphi,jphi,qp) + & - CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * phi_m_sq + jac_wt * CS%newton_drag_coef(i,j) * CS%newton_umid(i,j)**2 * phi_m_sq endif endif endif @@ -3161,7 +3195,7 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, ux = 0. uy = 0. - v_diag_qp(iphi,jphi,qp) = & + v_diag_qp(iphi,jphi,qp) = jac_wt * & ice_visc(i,j,qpv) * (((uy+vx) * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j)) + & ((4*vy+2*ux) * Phi(2*(2*(jphi-1)+iphi),qp,i,j))) @@ -3170,17 +3204,17 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, dstrain_diag_v = strsh_n * 0.5 * Phi(2*(2*(jphi-1)+iphi)-1,qp,i,j) + & (2.*stry_n + strx_n) * Phi(2*(2*(jphi-1)+iphi),qp,i,j) v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & - CS%newton_visc_factor(i,j,qpv) * dstrain_diag_v**2 + jac_wt * CS%newton_visc_factor(i,j,qpv) * dstrain_diag_v**2 endif if (float_cond(i,j) == 0) then vq = xquad(ilq) * xquad(jlq) v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & - (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) + jac_wt * (basal_trac(i,j) * vq) * (xquad(ilq) * xquad(jlq)) ! Newton basal drag diagonal correction: newton_drag_coef * (vmid_i)^2 * phi_m^2 if (CS%doing_newton) then v_diag_qp(iphi,jphi,qp) = v_diag_qp(iphi,jphi,qp) + & - CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * phi_m_sq + jac_wt * CS%newton_drag_coef(i,j) * CS%newton_vmid(i,j)**2 * phi_m_sq endif endif endif @@ -3205,7 +3239,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, if (float_cond(i,j) == 1) then Hcell(:,:) = H_node(i-1:i,j-1:j) - call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground) + call CG_diagonal_subgrid_basal(Phisub, Hcell, CS%bed_elev(i,j), dens_ratio, sub_ground, & + G%dxCv(i,j-1), G%dxCv(i,j), G%dyCu(i-1,j), G%dyCu(i,j), G%IareaT(i,j)) if (CS%umask(I-1,J-1) == 1) u_diag_b(I-1,J-1,4) = u_diag_b(I-1,J-1,4) + (sub_ground(1,1) * basal_trac(i,j)) if (CS%umask(I-1,J ) == 1) u_diag_b(I-1,J ,2) = u_diag_b(I-1,J ,2) + (sub_ground(1,2) * basal_trac(i,j)) @@ -3248,7 +3283,8 @@ subroutine matrix_diagonal(CS, G, US, float_cond, H_node, ice_visc, basal_trac, end subroutine matrix_diagonal -subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd) +subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd, & + dxCv_S, dxCv_N, dyCu_W, dyCu_E, IareaT) real, dimension(:,:,:,:,:,:), & intent(in) :: Phisub !< Quadrature structure weights at subgridscale !! locations for finite element calculations [nondim] @@ -3259,17 +3295,24 @@ subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd !! of seawater [nondim] real, dimension(2,2), intent(out) :: f_grnd !< The weighted fraction of the sub-cell where the ice shelf !! is grounded [nondim] + real, intent(in) :: dxCv_S !< The cell width at the southern (v-point) edge [L ~> m] + real, intent(in) :: dxCv_N !< The cell width at the northern (v-point) edge [L ~> m] + real, intent(in) :: dyCu_W !< The cell height at the western (u-point) edge [L ~> m] + real, intent(in) :: dyCu_E !< The cell height at the eastern (u-point) edge [L ~> m] + real, intent(in) :: IareaT !< The inverse of the cell area at the tracer point [L-2 ~> m-2] real, dimension(SIZE(Phisub,3),SIZE(Phisub,3),2,2) :: f_grnd_sub ! The contributions to nodal f_grnd !! from each sub-cell integer, dimension(2,2,SIZE(Phisub,3),SIZE(Phisub,3)) :: grnd_stat !0 at floating quad points, 1 at grounded real, dimension(2,2) :: f_grnd_q !Contributions to a node from each quadrature point in a sub-grid cell - real :: subarea ! The fractional sub-cell area [nondim] - real :: hloc ! The local sub-region thickness [Z ~> m] + real :: jac_sub_wt ! Per-sub-cell-QP metric correction |J_sub|/areaT [nondim] + real :: a, d ! Interpolated cell-edge spacings at the sub-cell QP [L ~> m] + real :: subarea ! The fractional sub-cell area in reference space [nondim] + real :: hloc ! The local sub-region thickness [Z ~> m] integer :: nsub, i, j, qx, qy, m, n nsub = size(Phisub,3) - subarea = 1.0 / (nsub**2) + subarea = 1.0 / real(nsub)**2 grnd_stat(:,:,:,:)=0 @@ -3281,10 +3324,20 @@ subroutine CG_diagonal_subgrid_basal (Phisub, H_node, bathyT, dens_ratio, f_grnd do n=1,2 ; do m=1,2 ; do j=1,nsub ; do i=1,nsub do qy=1,2 ; do qx = 1,2 - f_grnd_q(qx,qy) = grnd_stat(qx,qy,i,j) * Phisub(qx,qy,i,j,m,n)**2 + ! Interpolate cell-edge metrics to the sub-cell QP using the bilinear shape function values + ! from bilinear_shape_functions_subgrid. Marginal sums of Phisub give the interpolation + ! weights: sum over k=1 nodes gives (1-y); k=2 gives y; l=1 gives (1-x); l=2 gives x. + ! This is analogous to jac_wt = CS%Jac(qp,i,j) * G%IareaT(i,j) in the regular routines. + a = dxCv_S * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,2,1)) + & ! (1-y) * dxCv_S + dxCv_N * (Phisub(qx,qy,i,j,1,2) + Phisub(qx,qy,i,j,2,2)) ! + y * dxCv_N + d = dyCu_W * (Phisub(qx,qy,i,j,1,1) + Phisub(qx,qy,i,j,1,2)) + & ! (1-x) * dyCu_W + dyCu_E * (Phisub(qx,qy,i,j,2,1) + Phisub(qx,qy,i,j,2,2)) ! + x * dyCu_E + jac_sub_wt = 0.25 * subarea * (a * d) * IareaT + + f_grnd_q(qx,qy) = jac_sub_wt * grnd_stat(qx,qy,i,j) * Phisub(qx,qy,i,j,m,n)**2 enddo ; enddo !calculate sub-cell contribution to each node by summing up quadrature point contributions from the sub-cell - f_grnd_sub(i,j,m,n) = (subarea * 0.25) * ((f_grnd_q(1,1) + f_grnd_q(2,2)) + (f_grnd_q(1,2)+f_grnd_q(2,1))) + f_grnd_sub(i,j,m,n) = (f_grnd_q(1,1) + f_grnd_q(2,2)) + (f_grnd_q(1,2)+f_grnd_q(2,1)) enddo ; enddo ; enddo ; enddo !sum up the sub-cell contributions to each node @@ -3863,12 +3916,14 @@ end subroutine bilinear_shape_functions !> This subroutine calculates the gradients of bilinear basis elements that are centered at the !! vertices of the cell using a locally orthogoal MOM6 grid. Values are calculated at !! points of gaussian quadrature. -subroutine bilinear_shape_fn_grid(G, i, j, Phi) +subroutine bilinear_shape_fn_grid(G, i, j, Phi, Jac) type(ocean_grid_type), intent(in) :: G !< The grid structure used by the ice shelf. integer, intent(in) :: i !< The i-index in the grid to work on. integer, intent(in) :: j !< The j-index in the grid to work on. real, dimension(8,4), intent(inout) :: Phi !< The gradients of bilinear basis elements at Gaussian !! quadrature points surrounding the cell vertices [L-1 ~> m-1]. + real, dimension(4), optional, intent(out) :: Jac !< Jacobian determinant |J_q| = a_q*d_q at each + !! Gaussian quadrature point [L2 ~> m2]. ! This subroutine calculates the gradients of bilinear basis elements that ! that are centered at the vertices of the cell. The values are calculated at @@ -3922,6 +3977,7 @@ subroutine bilinear_shape_fn_grid(G, i, j, Phi) Phi(2*node,qpoint) = ( (a * (2 * ynode - 3)) * xexp ) / (a*d) enddo + if (present(Jac)) Jac(qpoint) = a * d enddo end subroutine bilinear_shape_fn_grid @@ -4228,13 +4284,16 @@ subroutine ice_shelf_dyn_end(CS) deallocate(CS%u_shelf, CS%v_shelf) deallocate(CS%taudx_shelf, CS%taudy_shelf) + deallocate(CS%sx_shelf, CS%sy_shelf) deallocate(CS%t_shelf, CS%tmask) deallocate(CS%u_bdry_val, CS%v_bdry_val) deallocate(CS%u_face_mask, CS%v_face_mask) + deallocate(CS%u_flux_bdry_val, CS%v_flux_bdry_val) deallocate(CS%umask, CS%vmask) deallocate(CS%u_face_mask_bdry, CS%v_face_mask_bdry) deallocate(CS%h_bdry_val) deallocate(CS%float_cond) + if (associated(CS%calve_mask)) deallocate(CS%calve_mask) deallocate(CS%ice_visc, CS%AGlen_visc) deallocate(CS%newton_visc_factor, CS%newton_str_ux, CS%newton_str_vy, CS%newton_str_sh) @@ -4243,6 +4302,10 @@ subroutine ice_shelf_dyn_end(CS) deallocate(CS%OD_rt, CS%OD_av) deallocate(CS%t_bdry_val, CS%bed_elev) deallocate(CS%ground_frac, CS%ground_frac_rt) + if (associated(CS%Jac)) deallocate(CS%Jac) + if (associated(CS%Phi)) deallocate(CS%Phi) + if (associated(CS%Phisub)) deallocate(CS%Phisub) + if (associated(CS%PhiC)) deallocate(CS%PhiC) deallocate(CS) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 6ed2b8242d..671d265a0e 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -133,7 +133,7 @@ module MOM_set_visc ! Allocatable data arrays real, allocatable, dimension(:,:) :: cdrag_u !< The spatially varying quadratic drag coefficient [nondim] real, allocatable, dimension(:,:) :: cdrag_v !< The spatially varying quadratic drag coefficient [nondim] - real, allocatable, dimension(:,:) :: tideamp !< RMS tidal amplitude at h points [Z T-1 ~> m s-1] + real, allocatable, dimension(:,:) :: tideamp !< RMS tidal amplitude at h points [L T-1 ~> m s-1] ! Diagnostic arrays real, allocatable, dimension(:,:) :: bbl_u !< BBL mean U current [L T-1 ~> m s-1] real, allocatable, dimension(:,:) :: bbl_v !< BBL mean V current [L T-1 ~> m s-1] @@ -3298,7 +3298,7 @@ subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS allocate(CS%tideamp(isd:ied,jsd:jed), source=0.0) filename = trim(CS%inputdir) // trim(tideamp_file) call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename) - call MOM_read_data(filename, tideamp_var, CS%tideamp, G%domain, scale=US%m_to_Z*US%T_to_s) + call MOM_read_data(filename, tideamp_var, CS%tideamp, G%domain, scale=US%m_s_to_L_T) call pass_var(CS%tideamp,G%domain) endif endif diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90 index c7d0d4a255..dda8bf7df1 100644 --- a/src/tracer/MOM_tracer_registry.F90 +++ b/src/tracer/MOM_tracer_registry.F90 @@ -305,7 +305,7 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u character(len=48) :: flux_units ! The units for fluxes, either ! [units] m3 s-1 or [units] kg s-1. character(len=48) :: conv_units ! The units for flux convergences, either - ! [units] m2 s-1 or [units] kg s-1. + ! [units] m s-1 or [units] kg m-2 s-1. character(len=48) :: unit2 ! The dimensions of the tracer squared character(len=72) :: cmorname ! The CMOR name of this tracer. character(len=120) :: cmor_longname ! The CMOR long name of that variable. @@ -379,9 +379,9 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u y_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) Tr%id_hbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_hbd_diffy", & diag%axesCvL, Time, trim(flux_longname)//" diffusive meridional " //& - "flux from the horizontal boundary diffusion scheme", trim(flux_units), & - v_extensive=.true., & - x_cell_method='sum', conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) + "flux from the horizontal boundary diffusion scheme", & + trim(flux_units), v_extensive=.true., x_cell_method='sum', & + conversion=(US%L_to_m**2)*Tr%flux_scale*US%s_to_T) else Tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", & diag%axesCuL, Time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), & @@ -487,37 +487,32 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u if (Tr%diag_form == 1) then Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + conv_units, conversion=Tr%conv_scale*US%s_to_T, v_extensive=.true.) Tr%id_dfxy_cont_2d = register_diag_field("ocean_model", & trim(shortnm)//'_dfxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated neutral diffusion tracer content "//& - "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum') + "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T) Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & diag%axesTL, Time, "Horizontal boundary diffusion tracer content tendency for "//& trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + conv_units, conversion=Tr%conv_scale*US%s_to_T, v_extensive=.true.) Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", & trim(shortnm)//'_hbdxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated horizontal boundary diffusion tracer content "//& - "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum') + "tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T) else cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& trim(lowercase(flux_longname))//& ' content due to parameterized mesoscale neutral diffusion' Tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', & diag%axesTL, Time, "Neutral diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, & + conv_units, conversion=Tr%conv_scale*US%s_to_T, v_extensive=.true., & cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff', & cmor_long_name=trim(cmor_var_lname), & - cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & - x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + cmor_standard_name=trim(cmor_long_std(cmor_var_lname))) cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//& trim(lowercase(flux_longname))//& @@ -528,20 +523,17 @@ subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE, u "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & cmor_field_name=trim(Tr%cmor_tendprefix)//'pmdiff_2d', & cmor_long_name=trim(cmor_var_lname), & - cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), & - x_cell_method='sum', y_cell_method='sum') + cmor_standard_name=trim(cmor_long_std(cmor_var_lname))) Tr%id_hbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_hbdxy_cont_tendency', & diag%axesTL, Time, & "Horizontal boundary diffusion tracer content tendency for "//trim(shortnm), & - conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum', v_extensive=.true.) + conv_units, conversion=Tr%conv_scale*US%s_to_T, v_extensive=.true.) Tr%id_hbdxy_cont_2d = register_diag_field("ocean_model", & trim(shortnm)//'_hbdxy_cont_tendency_2d', & diag%axesT1, Time, "Depth integrated horizontal boundary diffusion of tracer "//& - "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T, & - x_cell_method='sum', y_cell_method='sum') + "content tendency for "//trim(shortnm), conv_units, conversion=Tr%conv_scale*US%s_to_T) endif Tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', & diag%axesTL, Time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), &