diff --git a/docs/ocean/design_docs/harmonic_analysis/main.tex b/docs/ocean/design_docs/harmonic_analysis/main.tex new file mode 100644 index 0000000000..b3831b8a71 --- /dev/null +++ b/docs/ocean/design_docs/harmonic_analysis/main.tex @@ -0,0 +1,439 @@ +\documentclass[11pt]{report} + +\usepackage{epsf,amsmath,amsfonts} +\usepackage{graphicx} +\usepackage{booktabs} +\usepackage{color} +\usepackage{natbib} +\usepackage{multirow} +\usepackage{hyperref} + +\setlength{\textwidth}{6.5in} +\setlength{\oddsidemargin}{0in} +\setlength{\evensidemargin}{0in} +\setlength{\textheight}{8in} +\setlength{\topmargin}{0.5in} + +\newcommand{\ds}{\displaystyle} +\setlength{\parskip}{1.2ex} +\setlength{\parindent}{0mm} +\newcommand{\bea}{\begin{eqnarray}} +\newcommand{\eea}{\end{eqnarray}} +\newcommand{\nn}{\nonumber} + +\def\Mark#1{{\textcolor{red}{\bf \textit{By Mark: } #1}}} % new text +\def\Phillip#1{{\textcolor{cyan}{\bf \textit{By PhilW: } #1}}} % new text +\def\Steven#1{{\textcolor{blue}{\bf \textit{By Steven: } #1}}} % new text +\def\Luke#1{{\textcolor{green}{\bf \textit{By Luke: } #1}}} % new text +\def\Comment#1{{\textcolor{magenta}{\bf \textit{Comment: } #1}}} % new text +\begin{document} + +\title{MPAS-O Ocean Tidal Harmonic Decomposition } +\author{Steven Brus, Brian Arbic, Kristin Barton, Nairita Pal, \\ Mark Petersen, Andrew Roberts, Joannes Westerink } + + + +\maketitle +\tableofcontents + +%----------------------------------------------------------------------- + +\chapter{Summary} + + + +%The purpose of this section is to summarize what capability is to be added to the MPAS system through this design process. It should be clear what new code will do that the current code does not. Summarizing the primary challenges with respect to software design and implementation is also appropriate for this section. Finally, this statement should contain general statement with regard to what is ``success.'' + +\section{MPAS-Ocean contributions} +Decomposing the computed sea surface height into the amplitude and phase of the known tidal constituent frequencies is a common and useful way to visualize and validate tidal model results. These types of global plots are also valuable in assessing how changes in the model (mesh resolution, bathymetry, parameterizations, etc.) effect the tides globally. + +In a harmonic decomposition, the sea surface is approximated as a harmonic series with known constituent frequencies. The amplitudes and phases associated with these frequencies must be solved for in a manner that most closely matches the time series of the actual computed sea surface height with the series approximation. Once these amplitudes and phases are known, they can be plotted globally for each tidal constituent and compared with highly accurate data-assimilated tidal models. + +This type of analysis is implemented as a new analysis member in MPAS-Ocean. This allows the time interval at which the analysis is computed and written to file to be controlled using existing code infrastructure. It also allows this capability to be easily turned off for simulations which do not include tidal forcing. + +The methodology implemented is based on the harmonic least squares method \cite{westerink_lsq} implemented in the ADCIRC model \cite{Luettich1992ADCIRC} and borrows routines from the ADCIRC source code. + +\section{Science Questions} +\begin{enumerate} + \item Assess the accuracy of MPAS-Ocean tidal simulations with observations and highly accurate, data-assimilated tidal models + \item Identify regions where the computed tidal solution is inaccurate to focus efforts on improving accuracy due to mesh resolution, bathymetry, dissipation parameterizations, etc. + \item Demonstrate how coupling with other E3SM components effects the computed tidal solution +\end{enumerate} + +%----------------------------------------------------------------------- + +\chapter{Requirements} + +\section{Requirement: The harmonic analysis should decompose the computed sea surface height into amplitude and phase values for each major tidal constituent} +The amplitudes and phases associated with the known constituent frequencies must be solved for in a manner that most closely matches the actual computed sea surface height with the series approximation. The method must be able to selectively and accurately separate closely spaced harmonics. + +\section{Requirement: Decomposition is calculated on-line during the simulations} +The harmonic analysis capability should not require writing high-frequency global time-series to file. This would result in large storage requirements and require a subsequent step to compute the decomposition. + +\section{Requirement: Harmonic analysis should be able to be restarted at any point in the simulation} +The restarted simulation should produce a result that is bit-for-bit with the uninterrupted simulation. + +\section{Requirement: The harmonic analysis should allow for general specification of the tidal constituent frequencies included in the analysis} +Each of the major constituents should have namelist options which specify whether a given constituent is used in the analysis. + +\section{Requirement: The time period over which the harmonic analysis is performed should be user-specified} +The harmonic analysis should not include the tidal-spin up period used to ramp the tidal forcing at the beginning of the simulation. Also, in some situations, additional time-series information does not improve the results of the decomposition. + +\section{Requirement: The interval at which the computed sea surface height is sampled for the harmonic analysis should be user-specified } +This sampling interval can be used to controls the balance of accuracy and efficiency in the harmonic decomposition. +%----------------------------------------------------------------------- + +\chapter{Algorithmic Formulations} +Information in sections \ref{sec:astronomical_forcing}-\ref{sec:goverining_eqns} is reproduced here from the Tides Design Document for completeness. + +\section{Astronomical forcing} +\label{sec:astronomical_forcing} +The Newtonian equilibrium tidal potential is given by $\eta_{EQ}$. For the three largest semidiurnal tidal constituents, $\eta_{EQ}$ is given by equation (5) from \citet{chassignet_primer_2018} for $i=$ M$_{2}$, S$_{2}$, N$_{2}$ (applied 3 times and summed together), viz: +\begin{equation} +\eta_{EQ,i} = A_if_i(t_{ref})(1+k_{2,i}-h_{2,i})\cos^2(\phi)\cos\left[\omega_i(t-t_{ref}) + \chi_i(t_{ref}) + \nu_i(t_{ref}) + 2\lambda\right], +\label{eq:Eq5} +\end{equation} + +where the tidal amplitude is $A$, $f(t_{ref})$ is the the nodal factor accounting for the small modulations of amplitude (computed about once per year), $f(t_{ref})$ is slow modulation of amplitude (computed about once per year), the Love numbers $k_2$ and $h_2$ respectively account for the deformation of the solid earth due to the astronomical tidal forcing and the perturbation gravitational potential due to this deformation, $\phi$ is latitude, $\omega$ is tidal frequency, $t_{ref}$ is a reference time (often taken to be the beginning of a model run), $\chi(t_{ref})$ is the astronomical argument, and $\nu(t_ref)$ is the nodal factor accounting for the small modulations of phase, and $\lambda$ is longitude. + +For the diurnal constituents, the equilibrium tide is given by equation (6) of \citet{chassignet_primer_2018}, applied twice for $j=$ K$_{1}$ and O$_{1}$, viz: + +\begin{equation} +\eta_{EQ,j} = A_j f_j(t_{ref}) (1+k_{2,j}-h_{2,j})\sin(2\phi)\cos\left[\omega_j(t-t_{ref}) + \chi_j(t_{ref}) + \nu_j(t_{ref}) + \lambda \right]. +\label{eq:Eq6} +\end{equation} + +\section{Self-attraction and loading} +\label{sec:SAL} + +Self-attraction and loading (SAL) constituent static file maps should be derived for SAL amplitude and phase from FES, for $k=$ M$_2$, S$_2$, N$_2$, K$_1$, and O$_1$. The maps can be used to compute $\eta_{SAL}$ as a sum, e.g., equation (12) from \citet{chassignet_primer_2018}, + +\begin{equation} + \eta_{k,SAL}(\phi,\lambda) = A_m(\phi,\lambda)f(t_{ref})\cos\left[\omega (t-t_{ref}) + \chi(t_{ref}) + \nu(t_{ref}) - \phi_p(\phi,\lambda)\right], +\label{eq:SAL} +\end{equation} +where $A_m(\phi,\lambda)$ is the amplitude of the SAL of the $k$ constituent as a function of latitude ($\lambda$) and longitude ($\phi$) and $\phi_p(\phi,\lambda)$ is the phase of SAL as function of lat/lon. +Equation \ref{eq:SAL} takes amplitude and phase maps that provide a prediction with amplitude and phase using nodal factors and astronomical arguments. + +The self-attraction and loading harmonic constituents $A_m$ and $\phi_p$ can be derived from a harmonic analysis of self-attraction and loading maps produced by global tidal models, e.g., TPXO (\url{http://volkov.oce.orst.edu/tides/global.html}) or TUGO-m (\url{http://sirocco.obs-mip.fr/ocean-models/tugo/}). + +\section{Incorporation into governing equations} +\label{sec:goverining_eqns} + +Once $\eta_{EQ}$ and $\eta_{SAL}$ are obtained, they can be used in an ocean model. For instance, in a shallow water model, we replace the term $\nabla \eta$ in the momentum equation with a gradient of $\eta$ referenced to the equilibrium tide and self-attraction and loading terms, viz: + +\begin{equation} + \nabla\eta \rightarrow \nabla\left( \eta - \eta_{EQ} - \eta_{SAL}\right) +\end{equation} + +In essence, the equilibrium tide and self-attraction and loading terms reset the reference against which pressure gradients are computed. + + +\begin{table} +\begin{center} +\begin{tabular}{|c|c|c|c|c|} +\hline +Constituent & $\omega\;\left(10^{-4}\,s^{-1}\right)$ & $A\;\left(\textrm{cm}\right)$ & $1+k_{2}-h_{2}$ & Period (solar days)\tabularnewline +\hline +\hline +$\textrm{M}_{m}$ & 0.026392 & 2.2191 & 0.693 & 27.5546\tabularnewline +\hline +$\textrm{M}_{f}$ & 0.053234 & 4.2041 & 0.693 & 13.6608\tabularnewline +\hline +$\textrm{Q}_{1}$ & 0.6495854 & 1.9273 & 0.695 & 1.1195\tabularnewline +\hline +$\textrm{O}_{1}$ & 0.6759774 & 10.0661 & 0.695 & 1.0758\tabularnewline +\hline +$\textrm{P}_{1}$ & 0.7252295 & 4.6848 & 0.706 & 1.0027\tabularnewline +\hline +$\textrm{K}_{1}$ & 0.7292117 & 14.1565 & 0.736 & 0.9973\tabularnewline +\hline +$\textrm{N}_{2}$ & 1.378797 & 4.6397 & 0.693 & 0.5274\tabularnewline +\hline +$\textrm{M}_{2}$ & 1.405189 & 24.2334 & 0.693 & 0.5175\tabularnewline +\hline +$\textrm{S}_{2}$ & 1.454441 & 11.2743 & 0.693 & 0.5000\tabularnewline +\hline +$\textrm{K}_{2}$ & 1.458423 & 3.0684 & 0.693 & 0.4986\tabularnewline +\hline +\end{tabular} +\label{tab:astronimcalFactors} +\caption{Constituent-dependent frequencies $\omega$, astronomical forcing amplitudes A, and Love number combinations $1 + k_2 - h_2$ used to compute equilibrium tide $\eta_{EQ}$. The periods $2\pi/\omega$ are also given. Reproduced from Table 1 of \citet{chassignet_primer_2018, arbic2004accuracy}.} +\end{center} +\end{table} + + +\section{Harmonic Least Squares Method} +The method outlined here is described in more detail in \cite{westerink_lsq}. +The computed sea-surface height $\eta(t)$ is sampled to obtain values $\eta(t_k)$ at time sampling points $t_k,~k=1,M$. The following harmonic series is used to approximate the sea surface height in the least-squares procedure: +\begin{align} + g(t) = \sum_{j=1}^N(a_j\cos\omega_j t + b_j\sin\omega_j t). +\end{align} +In this equation, $a_j$ and $b_j$ are the unknown harmonic coefficients and $\omega_j$ are the frequencies associated with the $j=1,N$ tidal constituents. The error at a fiven sampling time between the approximate harmonic series and the computed sea surface height is +\begin{align} + \varepsilon(t_k) = \sum_{j=1}^N(a_j\cos\omega_j t_k + b_j\sin\omega_j t_k) -\eta(t_k). +\end{align} +The sum of the squared error at all sampling times is +\begin{align} + E = \sum_{k=1}^M\left(\varepsilon(t_k) \right)^2. +\end{align} +The error minimization between $g(t)$ and $f(t)$ is accomplished by requiring the partial derivatives of $E$with respect to each of the $a_j$ and $b_j$ coefficients to be zero: +\begin{align} + \frac{\partial E}{\partial a_i} &= 0, \quad i=1,N \\ + \frac{\partial E}{\partial b_i} &= 0, \quad i=1,N +\end{align} +Substituting for $E$: +\begin{align} + \sum_{k=1}^M \left(2\varepsilon(t_k)\frac{\partial \varepsilon(t_k)}{\partial a_i}\right) &= 0, \quad i=1,N \\ + \sum_{k=1}^M \left(2\varepsilon(t_k)\frac{\partial \varepsilon(t_k)}{\partial b_i}\right) &= 0, \quad i=1,N. +\end{align} +Substituting for the partial derivatives of $\varepsilon$ and dividing through by 2: +\begin{align} + \sum_{k=1}^M \left(\varepsilon(t_k)\cos\omega_i t_k\right) &= 0, \quad i=1,N \\ + \sum_{k=1}^M \left(\varepsilon(t_k)\sin\omega_i t_k\right) &= 0, \quad i=1,N +\end{align} +Substituting for $\varepsilon$ and rearranging: +\begin{align} + \sum_{k=1}^M \left(\sum_{j=1}^N(a_j\cos\omega_j t_k + b_j\sin\omega_j t_k)\cos\omega_i t_k\right) &= \sum_{k=1}^M \eta(t_k)\cos\omega_i t_k, \quad i=1,N \\ + \sum_{k=1}^M \left(\sum_{j=1}^N(a_j\cos\omega_j t_k + b_j\sin\omega_j t_k)\sin\omega_i t_k\right) &= \sum_{k=1}^M \eta(t_k)\sin\omega_i t_k, \quad i=1,N +\end{align} +Rearranging further: +\begin{align} + \sum_{j=1}^N \left(a_j\sum_{k=1}^M(\cos\omega_j t_k\cos\omega_i t_k) + b_j\sum_{k=1}^M(\sin\omega_j t_k\cos\omega_i t_k)\right) &= \sum_{k=1}^M \eta(t_k)\cos\omega_i t_k, \quad i=1,N \\ + \sum_{j=1}^N \left(a_j\sum_{k=1}^M(\cos\omega_j t_k\sin\omega_i t_k) + b_j\sum_{k=1}^M(\sin\omega_j t_k\sin\omega_i t_k)\right) &= \sum_{k=1}^M \eta(t_k)\sin\omega_i t_k, \quad i=1,N +\end{align} +This leads to the system of equations (example for 2 constituents): +\begin{align} + \underbrace{\sum_{k=1}^M \begin{bmatrix} + \cos\omega_1 t_k\cos\omega_1 t_k & \sin\omega_1 t_k\cos\omega_1 t_k & \cos\omega_2 t_k\cos\omega_1 t_k & \sin\omega_2 t_k\cos\omega_1 t_k\\ + \cos\omega_1 t_k\sin\omega_1 t_k & \sin\omega_1 t_k\sin\omega_1 t_k & \cos\omega_2 t_k\sin\omega_1 t_k & \sin\omega_2 t_k\sin\omega_1 t_k\\ \cos\omega_1 t_k\cos\omega_2 t_k & \sin\omega_1 t_k\cos\omega_2 t_k & \cos\omega_2 t_k\cos\omega_2 t_k & \sin\omega_2 t_k\cos\omega_2 t_k\\ + \cos\omega_1 t_k\sin\omega_2 t_k & \sin\omega_1 t_k\sin\omega_2 t_k & \cos\omega_2 t_k\sin\omega_2 t_k & \sin\omega_2 t_k\sin\omega_2 t_k\\ + \end{bmatrix}}_{\mathbf{M}} + \underbrace{\begin{bmatrix} + a_1 \\ b_1 \\ a_2 \\ b_2 + \end{bmatrix}}_{\mathbf{a}} + = + \underbrace{\sum_{k=1}^M\begin{bmatrix} + \eta(t_k)\cos\omega_1 t_k \\ + \eta(t_k)\sin\omega_1 t_k \\ + \eta(t_k)\cos\omega_2 t_k \\ + \eta(t_k)\sin\omega_2 t_k + \end{bmatrix}}_{\mathbf{f}} +\end{align} +For each mesh cell $l=1,K$ the linear system for that cell is solved +\begin{align} + \mathbf{M}\mathbf{a}_l = \mathbf{f}_l. +\end{align} +Note that the LHS matrix $\mathbf{M}$ is the same for each cell. In addition, $\mathbf{M}$ is symmetric and in the code, only the lower triangular part is stored. The upper triangular part is filled out before the matrix is decomposed prior to solving each individual linear system for the mesh cells. +Once the $a_j$ and $b_j$ coefficients are known, the amplitude and phase the the $j^\mathrm{th}$ constituent can be calculated as follows: +\begin{align} + A_j &= \frac{\sqrt{a_j^2+b_j^2}}{f(t_{ref})}, \\ + \Phi_j &= \arctan2(a_j,b_j) + \chi(t_{ref}) + \nu(t_{ref}), +\end{align} +%----------------------------------------------------------------------- + +\chapter{Design and Implementation} + +\section{Namelist options} +\begin{verbatim} +&AM_harmonicAnalysis + config_AM_harmonicAnalysis_enable = .true. + config_AM_harmonicAnalysis_compute_interval = 00:30:00 + config_AM_harmonicAnalysis_start = 2012-11-01_00:00:00 + config_AM_harmonicAnalysis_end = 2012-12-20_00:00:00 + config_AM_harmonicAnalysis_output_stream = 'harmonicAnalysisOutput' + config_AM_harmonicAnalysis_restart_stream = 'harmonicAnalysisRestart' + config_AM_harmonicAnalysis_compute_on_startup = .false. + config_AM_harmonicAnalysis_write_on_startup = .false. + config_AM_harmonicAnalysis_use_M2 = .true. + config_AM_harmonicAnalysis_use_S2 = .true. + config_AM_harmonicAnalysis_use_N2 = .true. + config_AM_harmonicAnalysis_use_K2 = .true. + config_AM_harmonicAnalysis_use_K1 = .true. + config_AM_harmonicAnalysis_use_O1 = .true. + config_AM_harmonicAnalysis_use_Q1 = .true. + config_AM_harmonicAnalysis_use_P1 = .true. +\end{verbatim} +\section{Output stream} +\begin{verbatim} + + + + + + + + + + + + + + + + + + + + +\end{verbatim} + +\section{Restart stream} +\begin{verbatim} + + + + + +\end{verbatim} + +\section{Subroutine calls} +The subroutines involved in computing the harmonic decomposition and how they function based on throughout the harmonic analysis period is summarized in \ref{tab:subroutines} +\begin{table}[ht] +\caption{Subroutine calls for harmonic analysis AM} +\begin{center} +\begin{tabular}{ll} + \hline + \multirow{5}{*}{Initialization} & {\tt ocn\_init\_harmonic\_analysis()} \\ + &~~~~~ determine constituents selected in namelist \\ + &~~~~~ {\tt tidal\_constituent\_factors()} \\ + &~~~~~~~~~~ compute tidal constituent nodal factors \\ + &~~~~~ zero out least squares matrix and vector \\ \cline{2-2} + \multirow{2}{*}{$t < t_{HAstart}$} & {\tt ocn\_compute\_harmonic\_analysis()} \\ + &~~~~~ return \\ \cline{2-2} + \multirow{5}{*}{$t > t_{HAstart}$ \& $t < t_{HAend}$} & {\tt ocn\_compute\_harmonic\_analysis()} \\ + &~~~~~ {\tt update\_least\_squares\_LHS\_matrix()} \\ + &~~~~~~~~~~ sum $t_k$ contribution into $\mathbf{M}$ \\ + &~~~~~ {\tt update\_least\_squares\_RHS\_vector()} \\ + &~~~~~~~~~~ sum $t_k$ contribution into each $\mathbf{f}_l$\\ \cline{2-2} + \multirow{8}{*}{$t > t_{HAend}$ (first time)} & {\tt ocn\_compute\_harmonic\_analysis()} \\ + &~~~~~ {\tt harmonic\_analysis\_solve()}\\ + &~~~~~~~~~~ {\tt least\_squares\_decompose()}\\ + &~~~~~~~~~~~~~~~ decompose $\mathbf{M}$ once\\ + &~~~~~~~~~~ {\tt least\_squares\_solve()}\\ + &~~~~~~~~~~~~~~~ solve each $\mathbf{M}\mathbf{a}_l=\mathbf{f}_l$ system\\ + &~~~~~~~~~~ compute $A$ and $\Phi$ for each mesh cell and constituent\\ + &~~~~~ separate solutions into output arrays \\ \cline{2-2} + \multirow{2}{*}{$t > t_{HAend}$ (subsequent times)} & {\tt ocn\_compute\_harmonic\_analysis()} \\ + &~~~~~ return \\ + \hline +\end{tabular} +\end{center} +\label{tab:subroutines} +\end{table} + +%----------------------------------------------------------------------- + +\chapter{Testing and Validation} + +The test case used is a 3-month global tidal simulation. This includes a 15 day tidal spin-up period. The harmonic analysis begins after the completion of the first 30 days and ends after an additional period of 49 days. The horizontal mesh can be found in Figure \ref{fig:mesh}. Globally, it is a 120km mesh with 30km resolution in the northwest Atlantic Ocean and 10km resolution in the Delaware bay region. The vertical mesh is the {\tt 100layerE3SMv1} configuration. Tidal potential forcing for the 8 major constituents ($M_2$, $S_2$, $N_2$, $K_2$, $K_1$, $O_1$, $Q_1$, $P_1$) is applied. The split explicit time integrator is used with a 5 minute baroclinic time step and a 15 second barotropic time step. + +The MPAS-Ocean results are compared with the data-assimilated TPXO8 model. Given the coarse mesh resolution and the lack of tuning, the results are not very accurate. However, the comparison does show that the harmonic analysis is being computed correctly. A direct water level comparison is shown in Figure \ref{fig:water_level}. The amplitude and phase plots are shown in \ref{fig:M2Amplitude}-\ref{fig:K1Phase} + +\begin{figure} + \centering + \includegraphics[width=4in]{cellWidth.png} + \caption{USDEQU120at30cr10} + \label{fig:mesh} +\end{figure} +\begin{figure} + \centering + \includegraphics[width=3.5in]{8537121.png} + \caption{Water level comparison} + \label{fig:water_level} +\end{figure} +\begin{figure} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{M2Amplitude_comparison_compressed.png} + \caption{$M_2$ amplitude comparison} + \label{fig:M2Amplitude} + \end{minipage} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{M2Phase_comparison_compressed.png} + \caption{$M_2$ phase comparison} + \label{fig:M2Phase} + \end{minipage} +\end{figure} +\begin{figure} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{N2Amplitude_comparison_compressed.png} + \caption{$N_2$ amplitude comparison} + \label{fig:N2Amplitude} + \end{minipage} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{N2Phase_comparison_compressed.png} + \caption{$N_2$ phase comparison} + \label{fig:N2Phase} + \end{minipage} +\end{figure} +\begin{figure} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{S2Amplitude_comparison_compressed.png} + \caption{$S_2$ amplitude comparison} + \label{fig:S2Amplitude} + \end{minipage} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{S2Phase_comparison_compressed.png} + \caption{$S_2$ phase comparison} + \label{fig:S2Phase} + \end{minipage} +\end{figure} +\begin{figure} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{O1Amplitude_comparison_compressed.png} + \caption{$O_1$ amplitude comparison} + \label{fig:O1Amplitude} + \end{minipage} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{O1Phase_comparison_compressed.png} + \caption{$O_1$ phase comparison} + \label{fig:O1Phase} + \end{minipage} +\end{figure} +\begin{figure} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{K1Amplitude_comparison_compressed.png} + \caption{$K_1$ amplitude comparison} + \label{fig:K1Amplitude} + \end{minipage} + \begin{minipage}{.5\textwidth} + \centering + \includegraphics[width=2.75in]{K1Phase_comparison_compressed.png} + \caption{$K_1$ phase comparison} + \label{fig:K1Phase} + \end{minipage} +\end{figure} + +%----------------------------------------------------------------------- +\bibliographystyle{unsrtnat} +\bibliography{tides} + + +\end{document} + diff --git a/docs/ocean/design_docs/harmonic_analysis/tides.bib b/docs/ocean/design_docs/harmonic_analysis/tides.bib new file mode 100644 index 0000000000..0373ae0bff --- /dev/null +++ b/docs/ocean/design_docs/harmonic_analysis/tides.bib @@ -0,0 +1,48 @@ +@incollection{chassignet_primer_2018, + title = {A {Primer} on {Global} {Internal} {Tide} and {Internal} {Gravity} {Wave} {Continuum} {Modeling} in {HYCOM} and {MITgcm}}, + isbn = {978-1-72054-997-0}, + url = {http://purl.flvc.org/fsu/fd/FSU\_libsubv1\_scholarship\_submission\_1536242074\_55feafcc}, + language = {en}, + urldate = {2018-10-23}, + booktitle = {New {Frontiers} in {Operational} {Oceanography}}, + publisher = {GODAE OceanView}, + author = {Arbic, Brian K. and Alford, Matthew H. and Ansong, Joseph K. and Buijsman, Maarten C. and Ciotti, Robert B. and Farrar, J. Thomas and Hallberg, Robert W. and Henze, Christopher E. and Hill, Christopher N. and Luecke, Conrad A. and Menemenlis, Dimitris and Metzger, E. Joseph and Müeller, Malte and Nelson, Arin D. and Nelson, Bron C. and Ngodock, Hans E. and Ponte, Rui M. and Richman, James G. and Savage, Anna C. and Scott, Robert B. and Shriver, Jay F. and Simmons, Harper L. and Souopgui, Innocent and Timko, Patrick G. and Wallcraft, Alan J. and Zamudio, Luis and Zhao, Zhongxiang}, + editor = {Chassignet, Eric P. and Pascual, Ananda and Tintoré, Joaquin and Verron, Jacques}, + month = aug, + year = {2018}, + doi = {10.17125/gov2018.ch13}, + file = {Arbic et al. - 2018 - A Primer on Global Internal Tide and Internal Grav.pdf:/Users/pwolfram/Zotero/storage/TEDRP9CP/Arbic et al. - 2018 - A Primer on Global Internal Tide and Internal Grav.pdf:application/pdf} +} + +@article{Luettich1992ADCIRC, + author={Luettich, R.A. and Westerink, J.J. and Scheffner, N.W.}, + title = {{ADCIRC}: {An} {Advanced} {Three}-{Dimensional} {Circulation} {Model} for {Shelves}, {Coasts}, and {Estuaries}. {Report} 1. {Theory} and {Methodology} of {ADCIRC-2DDI} and {ADCIRC-3DL}.}, + language = {en}, + pages = {143}, + year={1992} +} + +@article{arbic2004accuracy, + title={The accuracy of surface elevations in forward global barotropic and baroclinic tide models}, + author={Arbic, Brian K and Garner, Stephen T and Hallberg, Robert W and Simmons, Harper L}, + journal={Deep Sea Research Part II: Topical Studies in Oceanography}, + volume={51}, + number={25-26}, + pages={3069--3101}, + year={2004}, + publisher={Elsevier} +} + +@article{westerink_lsq, +author = {Westerink, J. J. and Connor, J. J. and Stolzenbach, K. D.}, +title = {A frequency–time domain finite element model for tidal circulation based on the least-squares harmonic analysis method}, +journal = {International Journal for Numerical Methods in Fluids}, +volume = {8}, +number = {7}, +pages = {813-843}, +keywords = {Shallow Water Equations, Iterative, Harmonic Analysis, Least Squares, Finite Element, Tides}, +doi = {10.1002/fld.1650080706}, +url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650080706}, +eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.1650080706}, +year = {1988} +} diff --git a/src/core_ocean/Registry.xml b/src/core_ocean/Registry.xml index 81f8b7fbbc..be834b1c05 100644 --- a/src/core_ocean/Registry.xml +++ b/src/core_ocean/Registry.xml @@ -97,7 +97,11 @@ /> + @@ -544,6 +548,10 @@ /> + + + + + + + + + + + + @@ -1354,6 +1394,7 @@ + @@ -1760,6 +1801,7 @@ + + @@ -3008,6 +3054,27 @@ packages="tidalForcing" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/core_ocean/analysis_members/Registry_sediment_flux_index.xml b/src/core_ocean/analysis_members/Registry_sediment_flux_index.xml new file mode 100644 index 0000000000..0baed7a34a --- /dev/null +++ b/src/core_ocean/analysis_members/Registry_sediment_flux_index.xml @@ -0,0 +1,101 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/core_ocean/analysis_members/Registry_sediment_transport.xml b/src/core_ocean/analysis_members/Registry_sediment_transport.xml new file mode 100644 index 0000000000..06d0c622da --- /dev/null +++ b/src/core_ocean/analysis_members/Registry_sediment_transport.xml @@ -0,0 +1,272 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/core_ocean/analysis_members/mpas_ocn_analysis_driver.F b/src/core_ocean/analysis_members/mpas_ocn_analysis_driver.F index 86157ce16d..88781025c8 100644 --- a/src/core_ocean/analysis_members/mpas_ocn_analysis_driver.F +++ b/src/core_ocean/analysis_members/mpas_ocn_analysis_driver.F @@ -51,6 +51,9 @@ module ocn_analysis_driver use ocn_moc_streamfunction use ocn_ocean_heat_content use ocn_mixed_layer_heat_budget + use ocn_sediment_flux_index + use ocn_sediment_transport + use ocn_harmonic_analysis ! use ocn_TEM_PLATE implicit none @@ -198,6 +201,9 @@ subroutine ocn_analysis_setup_packages(configPool, packagePool, iocontext, err)! call mpas_pool_add_config(analysisMemberList, 'eddyProductVariables', 1) call mpas_pool_add_config(analysisMemberList, 'mocStreamfunction', 1) call mpas_pool_add_config(analysisMemberList, 'oceanHeatContent', 1) + call mpas_pool_add_config(analysisMemberList, 'sedimentFluxIndex', 1) + call mpas_pool_add_config(analysisMemberList, 'sedimentTransport', 1) + call mpas_pool_add_config(analysisMemberList, 'harmonicAnalysis', 1) ! call mpas_pool_add_config(analysisMemberList, 'temPlate', 1) ! DON'T EDIT BELOW HERE @@ -1057,6 +1063,8 @@ subroutine ocn_bootstrap_analysis_members(domain, analysisMemberName, iErr)!{{{ call ocn_bootstrap_pointwise_stats(domain, err_tmp) else if ( analysisMemberName(1:nameLength) == 'transectTransport' ) then call ocn_init_transect_transport(domain, err_tmp) +! else if ( analysisMemberName(1:nameLength) == 'sedimentFluxIndex' ) then +! call ocn_init_sediment_flux_index(domain, err_tmp) ! else if ( analysisMemberName(1:nameLength) == 'temPlate' ) then ! call ocn_init_TEM_PLATE(domain, err_tmp) end if @@ -1128,6 +1136,12 @@ subroutine ocn_init_analysis_members(domain, analysisMemberName, iErr)!{{{ call ocn_init_ocean_heat_content(domain, err_tmp) else if ( analysisMemberName(1:nameLength) == 'mixedLayerHeatBudget' ) then call ocn_init_mixed_layer_heat_budget(domain,err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentFluxIndex' ) then + call ocn_init_sediment_flux_index(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentTransport' ) then + call ocn_init_sediment_transport(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'harmonicAnalysis' ) then + call ocn_init_harmonic_analysis(domain, err_tmp) ! else if ( analysisMemberName(1:nameLength) == 'temPlate' ) then ! call ocn_init_TEM_PLATE(domain, err_tmp) ! rpn is third to last @@ -1224,6 +1238,12 @@ subroutine ocn_compute_analysis_members(domain, timeLevel, analysisMemberName, i call ocn_compute_ocean_heat_content(domain, timeLevel, err_tmp) else if ( analysisMemberName(1:nameLength) == 'mixedLayerHeatBudget' ) then call ocn_compute_mixed_layer_heat_budget(domain, timeLevel, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentFluxIndex' ) then + call ocn_compute_sediment_flux_index(domain, timeLevel, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentTransport' ) then + call ocn_compute_sediment_transport(domain, timeLevel, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'harmonicAnalysis' ) then + call ocn_compute_harmonic_analysis(domain, timeLevel, err_tmp) ! else if ( analysisMemberName(1:nameLength) == 'temPlate' ) then ! call ocn_compute_TEM_PLATE(domain, timeLevel, err_tmp) ! rpn is third to last @@ -1325,6 +1345,12 @@ subroutine ocn_restart_analysis_members(domain, analysisMemberName, iErr)!{{{ call ocn_restart_ocean_heat_content(domain, err_tmp) else if ( analysisMemberName(1:nameLength) == 'mixedLayerHeatBudget' ) then call ocn_restart_mixed_layer_heat_budget(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentFluxIndex' ) then + call ocn_restart_sediment_flux_index(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentTransport' ) then + call ocn_restart_sediment_transport(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'harmonicAnalysis' ) then + call ocn_restart_harmonic_analysis(domain, err_tmp) ! else if ( analysisMemberName(1:nameLength) == 'temPlate' ) then ! call ocn_restart_TEM_PLATE(domain, err_tmp) ! rpn is third to last @@ -1420,6 +1446,12 @@ subroutine ocn_finalize_analysis_members(domain, analysisMemberName, iErr)!{{{ call ocn_finalize_ocean_heat_content(domain,err_tmp) else if ( analysisMemberName(1:nameLength) == 'mixedLayerHeatBudget' ) then call ocn_finalize_mixed_layer_heat_budget(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentFluxIndex' ) then + call ocn_finalize_sediment_flux_index(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'sedimentTransport' ) then + call ocn_finalize_sediment_transport(domain, err_tmp) + else if ( analysisMemberName(1:nameLength) == 'harmonicAnalysis' ) then + call ocn_finalize_harmonic_analysis(domain, err_tmp) ! else if ( analysisMemberName(1:nameLength) == 'temPlate' ) then ! call ocn_finalize_TEM_PLATE(domain, err_tmp) ! rpn is third to last diff --git a/src/core_ocean/analysis_members/mpas_ocn_harmonic_analysis.F b/src/core_ocean/analysis_members/mpas_ocn_harmonic_analysis.F new file mode 100644 index 0000000000..4d51fd608a --- /dev/null +++ b/src/core_ocean/analysis_members/mpas_ocn_harmonic_analysis.F @@ -0,0 +1,821 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_harmonic_analysis +! +!> \brief MPAS ocean analysis mode member: harmonic_analysis +!> \author Steven Brus +!> \date July 2020 +!> \details +!> MPAS ocean analysis mode member: harmonic_analysis +!> +!> This module contains subroutines to compute the harmonic decomposition +!> of the sea-surface-height timeseries into an amplitude and phase for a +!> given number of tidal constituents. +!> +!> The subroutines that compute the harmonic decomposition are taken from the +!> ADCIRC model developed by R.A. Luettich, Jr. (University of North Carolina at Chapel Hill) +!> and J.J. Westerink (University of Notre Dame). +!> +!----------------------------------------------------------------------- + +module ocn_harmonic_analysis + + use mpas_kind_types + use mpas_derived_types + use mpas_pool_routines + use mpas_timekeeping + use mpas_stream_manager + use mpas_constants + use ocn_constants + use ocn_config + use ocn_diagnostics_routines + use ocn_vel_tidal_potential, only: char_array,tidal_constituent_factors + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_init_harmonic_analysis, & + ocn_compute_harmonic_analysis, & + ocn_restart_harmonic_analysis, & + ocn_finalize_harmonic_analysis + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + + type(char_array), dimension(37) :: constituentList + integer :: leastSquaresSolution + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_init_harmonic_analysis +! +!> \brief Initialize MPAS-Ocean analysis member +!> \author Steven Brus +!> \date July 2020 +!> \details +!> This routine conducts all initializations required for the +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_init_harmonic_analysis(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + type (dm_info) :: dminfo + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: harmonicAnalysisAMPool + + integer, pointer :: nAnalysisConstituents + real (kind=RKIND), pointer :: harmonicAnalysisStart + real (kind=RKIND), pointer :: harmonicAnalysisEnd + real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix + real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector + real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentAmplitude + real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentPhase + real (kind=RKIND), dimension(:), pointer :: analysisConstituentFrequency + real (kind=RKIND), dimension(:), pointer :: analysisConstituentNodalAmplitude + real (kind=RKIND), dimension(:), pointer :: analysisConstituentNodalPhase + real (kind=RKIND), dimension(:), allocatable :: tidalConstituentAmplitude + real (kind=RKIND), dimension(:), allocatable :: tidalConstituentLoveNumbers + real (kind=RKIND), dimension(:), allocatable :: tidalConstituentAstronomical + integer, dimension(:), allocatable :: tidalConstituentType + + type (MPAS_Time_Type) :: refTime + integer :: iCon + + err = 0 + + leastSquaresSolution = 0 + + dminfo = domain % dminfo + + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'harmonicAnalysisAM', harmonicAnalysisAMPool) + + call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentFrequency', analysisConstituentFrequency) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentNodalAmplitude', analysisConstituentNodalAmplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentNodalPhase', analysisConstituentNodalPhase) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'decomposedConstituentAmplitude', decomposedConstituentAmplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'decomposedConstituentPhase', decomposedConstituentPhase) + + call mpas_set_time(refTime, dateTimeString=config_tidal_potential_reference_time) + + nAnalysisConstituents = 0 + if (config_AM_harmonicAnalysis_use_M2) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'M2' + end if + + if (config_AM_harmonicAnalysis_use_S2) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'S2' + end if + + if (config_AM_harmonicAnalysis_use_N2) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'N2' + end if + + if (config_AM_harmonicAnalysis_use_K2) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'K2' + end if + + if (config_AM_harmonicAnalysis_use_K1) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'K1' + end if + + if (config_AM_harmonicAnalysis_use_O1) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'O1' + end if + + if (config_AM_harmonicAnalysis_use_Q1) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'Q1' + end if + + if (config_AM_harmonicAnalysis_use_P1) then + nAnalysisConstituents = nAnalysisConstituents + 1 + constituentList(nAnalysisConstituents)%constituent = 'P1' + end if + + allocate(tidalConstituentAmplitude(nAnalysisConstituents)) + allocate(tidalConstituentLoveNumbers(nAnalysisConstituents)) + allocate(tidalConstituentAstronomical(nAnalysisConstituents)) + allocate(tidalConstituentType(nAnalysisConstituents)) + + call tidal_constituent_factors(constituentList,nAnalysisConstituents,refTime, & + analysisConstituentFrequency, & + tidalConstituentAmplitude, & + tidalConstituentLoveNumbers, & + analysisConstituentNodalAmplitude, & + tidalConstituentAstronomical, & + analysisConstituentNodalPhase, & + tidalConstituentType, & + err) + + if (.not. config_do_restart) then + leastSquaresRHSVector = 0.0_RKIND + leastSquaresLHSMatrix = 0.0_RKIND + end if + + + do iCon = 1,nAnalysisConstituents + call mpas_log_write('Constituent '//constituentList(iCon)%constituent) + call mpas_log_write(' Frequency $r', realArgs=(/ analysisConstituentFrequency(iCon) /)) + call mpas_log_write(' Amplitude $r', realArgs=(/ tidalConstituentAmplitude(iCon) /)) + call mpas_log_write(' LoveNumbers $r', realArgs=(/ tidalConstituentLoveNumbers(iCon) /)) + call mpas_log_write(' NodalAmplitude $r', realArgs=(/ analysisConstituentNodalAmplitude(iCon) /)) + call mpas_log_write(' Astronomical argument $r', realArgs=(/ tidalConstituentAstronomical(iCon) /)) + call mpas_log_write(' NodalPhase $r', realArgs=(/ analysisConstituentNodalPhase(iCon) /)) + call mpas_log_write(' Type $i', intArgs=(/ tidalConstituentType(iCon) /)) + call mpas_log_write(' ') + end do + + block => block % next + end do + + end subroutine ocn_init_harmonic_analysis!}}} + +!*********************************************************************** +! +! routine ocn_compute_harmonic_analysis +! +!> \brief Compute MPAS-Ocean analysis member +!> \author Steven Brus +!> \date July 2020 +!> \details +!> This routine conducts all computation required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + integer, intent(in) :: timeLevel + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + type (dm_info) :: dminfo + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: scratchPool + type (mpas_pool_type), pointer :: diagnosticsPool + type (mpas_pool_type), pointer :: harmonicAnalysisAMPool + type (mpas_time_type) :: current_time + type (mpas_time_type) :: harmonicAnalysisStart + type (mpas_time_type) :: harmonicAnalysisEnd + + integer, pointer :: nCellsSolve + integer, pointer :: nAnalysisConstituents + real (kind=RKIND), pointer :: daysSinceStartOfSim + real (kind=RKIND), dimension(:), pointer :: ssh + real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix + real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector + real (kind=RKIND), dimension(:), pointer :: analysisConstituentFrequency + real (kind=RKIND), dimension(:), pointer :: analysisConstituentNodalAmplitude + real (kind=RKIND), dimension(:), pointer :: analysisConstituentNodalPhase + real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentAmplitude + real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentPhase + real (kind=RKIND), dimension(:), pointer :: M2Amplitude + real (kind=RKIND), dimension(:), pointer :: M2Phase + real (kind=RKIND), dimension(:), pointer :: S2Amplitude + real (kind=RKIND), dimension(:), pointer :: S2Phase + real (kind=RKIND), dimension(:), pointer :: N2Amplitude + real (kind=RKIND), dimension(:), pointer :: N2Phase + real (kind=RKIND), dimension(:), pointer :: K2Amplitude + real (kind=RKIND), dimension(:), pointer :: K2Phase + real (kind=RKIND), dimension(:), pointer :: K1Amplitude + real (kind=RKIND), dimension(:), pointer :: K1Phase + real (kind=RKIND), dimension(:), pointer :: O1Amplitude + real (kind=RKIND), dimension(:), pointer :: O1Phase + real (kind=RKIND), dimension(:), pointer :: Q1Amplitude + real (kind=RKIND), dimension(:), pointer :: Q1Phase + real (kind=RKIND), dimension(:), pointer :: P1Amplitude + real (kind=RKIND), dimension(:), pointer :: P1Phase + + integer :: iCell + integer :: iCon + real (kind=RKIND) :: time + + err = 0 + + dminfo = domain % dminfo + + block => domain % blocklist + do while (associated(block)) + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) + call mpas_pool_get_subpool(block % structs, 'harmonicAnalysisAM', harmonicAnalysisAMPool) + + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_array(diagnosticsPool, "daysSinceStartOfSim", daysSinceStartOfSim) + call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentFrequency', analysisConstituentFrequency) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentNodalAmplitude', analysisConstituentNodalAmplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'analysisConstituentNodalPhase', analysisConstituentNodalPhase) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'decomposedConstituentAmplitude', decomposedConstituentAmplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'decomposedConstituentPhase', decomposedConstituentPhase) + + ! get relevant time information + time = daysSinceStartOfSim*86400_RKIND + current_time = mpas_get_clock_time(domain % clock, MPAS_NOW, err) + call mpas_set_time(harmonicAnalysisStart, dateTimeString=config_AM_harmonicAnalysis_start) + call mpas_set_time(harmonicAnalysisEnd, dateTimeString=config_AM_harmonicAnalysis_end) + + ! exit if harmonic analysis period has not begun yet + if (current_time .lt. harmonicAnalysisStart) then + call mpas_log_write('harmonicAnalysisAM exit: before HA period') + return + end if + + ! update if within harmonic analysis period + if (current_time .le. harmonicAnalysisEnd) then + CALL update_least_squares_LHS_matrix(nAnalysisConstituents, & + time, & + analysisConstituentFrequency, & + leastSquaresLHSMatrix) + + CALL update_least_squares_RHS_vector(nAnalysisConstituents, & + time, & + nCellsSolve, & + analysisConstituentFrequency, & + ssh, & + leastSquaresRHSVector) + call mpas_log_write('harmonicAnalysisAM update') + end if + + ! solve harmonic analysis least squares system if harmonic analysis period is over + if ((current_time .ge. harmonicAnalysisEnd) .and. (leastSquaresSolution == 0)) then + call harmonic_analysis_solve(nCellsSolve, & + nAnalysisConstituents, & + leastSquaresLHSMatrix, & + leastSquaresRHSVector, & + analysisConstituentNodalAmplitude, & + analysisConstituentNodalPhase, & + decomposedConstituentAmplitude, & + decomposedConstituentPhase) + + ! copy amplitude and phase solutions into corresponding arrays + do iCon = 1,nAnalysisConstituents + if (constituentList(iCon)%constituent == 'M2') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'M2Amplitude', M2Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'M2Phase', M2Phase) + M2Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + M2Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'S2') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'S2Amplitude', S2Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'S2Phase', S2Phase) + S2Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + S2Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'N2') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'N2Amplitude', N2Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'N2Phase', N2Phase) + N2Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + N2Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'K2') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'K2Amplitude', K2Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'K2Phase', K2Phase) + K2Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + K2Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'K1') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'K1Amplitude', K1Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'K1Phase', K1Phase) + K1Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + K1Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'O1') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'O1Amplitude', O1Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'O1Phase', O1Phase) + O1Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + O1Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'Q1') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'Q1Amplitude', Q1Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'Q1Phase', Q1Phase) + Q1Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + Q1Phase(:) = decomposedConstituentPhase(iCon,:) + endif + if (constituentList(iCon)%constituent == 'P1') then + call mpas_pool_get_array(harmonicAnalysisAMPool, 'P1Amplitude', P1Amplitude) + call mpas_pool_get_array(harmonicAnalysisAMPool, 'P1Phase', P1Phase) + P1Amplitude(:) = decomposedConstituentAmplitude(iCon,:) + P1Phase(:) = decomposedConstituentPhase(iCon,:) + endif + enddo + + ! increment to indicate the system has been solved + leastSquaresSolution = leastSquaresSolution + 1 + call mpas_log_write('harmonicAnalysisAM solve') + + else if ((current_time .ge. harmonicAnalysisEnd) .and. (leastSquaresSolution > 0)) then + call mpas_log_write('harmonicAnalysisAM exit: past HA period') + end if + + block => block % next + end do + + end subroutine ocn_compute_harmonic_analysis!}}} + +!*********************************************************************** +! +! routine ocn_restart_harmonic_analysis +! +!> \brief Save restart for MPAS-Ocean analysis member +!> \author Steven Brus +!> \date July 2020 +!> \details +!> This routine conducts computation required to save a restart state +!> for the MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_restart_harmonic_analysis(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_restart_harmonic_analysis!}}} + +!*********************************************************************** +! +! routine ocn_finalize_harmonic_analysis +! +!> \brief Finalize MPAS-Ocean analysis member +!> \author Steven Brus +!> \date July 2020 +!> \details +!> This routine conducts all finalizations required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_finalize_harmonic_analysis(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + + err = 0 + + end subroutine ocn_finalize_harmonic_analysis!}}} + +!*********************************************************************** +! +! routine update_least_squares_LHS_matrix +! +!> \brief Update Left Hand Side Matrix +!> \author Steven Brus +!> \date July 2020 +!> \details +!> Adapted from the ADCIRC code by R.A. Luettich and J.J. Westerink +! +!----------------------------------------------------------------------- + + SUBROUTINE update_least_squares_LHS_matrix(nfreq,TIMELOC,hafreq,ha) !{{{ + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nfreq ! number of analysis constituents + REAL(kind=RKIND), INTENT(IN) :: TIMELOC ! model time (sec) + REAL(kind=RKIND), DIMENSION(:), INTENT(IN) :: hafreq ! analysis contituent frequencies + REAL(kind=RKIND), DIMENSION(:,:), INTENT(INOUT) :: ha ! LHS matrix + + INTEGER :: I,J,I1,I2,J1,J2 + REAL(kind=RKIND) :: TF1,TF2 + +! Update the Left Hand Side Matrix +! Note: this is a symmetric matrix and therefore only store the +! upper triangular part. The lower part will be filled out in +! subroutine least_squares_decompose() prior to decomposition + + do i=1,nfreq + do j=i,nfreq + i1=2*i-1 + i2=i1+1 + j1=2*j-1 + j2=j1+1 + tf1=hafreq(i)*TIMELOC + tf2=hafreq(j)*TIMELOC + ha(i1,j1) = ha(i1,j1) + cos(tf1)*cos(tf2) + ha(i1,j2) = ha(i1,j2) + cos(tf1)*sin(tf2) + ha(i2,j2) = ha(i2,j2) + sin(tf1)*sin(tf2) + if(i2.le.j1) ha(i2,j1) = ha(i2,j1) + sin(tf1)*cos(tf2) + end do + end do + + return + end subroutine update_least_squares_LHS_matrix!}}} + +!*********************************************************************** +! +! routine update_least_squares_RHS_vector +! +!> \brief Update Right Hand Side Vector +!> \author Steven Brus +!> \date July 2020 +!> \details +!> Adapted from the ADCIRC code by R.A. Luettich and J.J. Westerink +! +!----------------------------------------------------------------------- + + SUBROUTINE update_least_squares_RHS_vector(nfreq,TIMEUD,NP,hafreq,GLOE,GLOELV) !{{{ + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: nfreq ! number of analysis consituents + REAL(kind=RKIND), INTENT(IN) :: TIMEUD ! model time + INTEGER, INTENT(IN) :: NP ! number of mesh cells + REAL(kind=RKIND), DIMENSION(:), INTENT(IN) :: hafreq ! analysis constituent frequencies + REAL(kind=RKIND), DIMENSION(:), INTENT(IN) :: GLOE ! sea surface height + REAL(kind=RKIND), DIMENSION(:,:), INTENT(INOUT) :: GLOELV ! RHS vector + + INTEGER I,N,I1,I2 + REAL(kind=RKIND) TF1,CTF1,STF1 + +! Update the Right Hand Side Load Vectors + do i=1,nfreq + i1=2*i-1 + i2=i1+1 + tf1=hafreq(i)*TIMEUD + ctf1 = cos(tf1) + stf1 = sin(tf1) + do n=1,np + GLOELV(I1,N)=GLOELV(I1,N)+GLOE(N)*CTF1 + GLOELV(I2,N)=GLOELV(I2,N)+GLOE(N)*STF1 + end do + end do + + return + end subroutine update_least_squares_RHS_vector!}}} + +!*********************************************************************** +! +! routine harmonic_analysis_solve +! +!> \brief Solve the least squares system +!> \author Steven Brus +!> \date July 2020 +!> \details +!> Adapted from the ADCIRC code by R.A. Luettich and J.J. Westerink +! +!----------------------------------------------------------------------- + + SUBROUTINE harmonic_analysis_solve(MNP,nfreq,hmat,GLOELV,haff,haface,emagt,phaseden) !{{{ + + IMPLICIT NONE + + integer, intent(in) :: MNP ! number of mesh cells + integer, intent(in) :: nfreq ! number of analysis constituents + real(kind=RKIND), dimension(:,:), intent(in) :: hmat ! LHS matrix + real(kind=RKIND), dimension(:,:), intent(in) :: GLOELV ! RHS vector + real(kind=RKIND), dimension(:), intent(in) :: haff ! amplitude nodal factors + real(kind=RKIND), dimension(:), intent(in) :: haface ! phase nodal factors + real(kind=RKIND), dimension(:,:), intent(out) :: emagt ! constituent amplitudes + real(kind=RKIND), dimension(:,:), intent(out) :: phaseden ! constituent phases + + integer J,N,K,I,I1,I2,IT,IFR + integer mm + real(kind=RKIND) :: convrd + REAL(kind=RKIND),ALLOCATABLE :: PHASEE(:),EMAG(:) + REAL(kind=RKIND),ALLOCATABLE :: hap(:),hax(:) + REAL(kind=RKIND),ALLOCATABLE :: ha(:,:) + + convrd=180.0_RKIND/pii + + mm = 2*nfreq + ALLOCATE ( PHASEE(nfreq),EMAG(nfreq) ) + ALLOCATE ( hap(mm), hax(mm) ) + ALLOCATE ( ha(mm,mm) ) + +! Copy LHS matrix and decompose + do i = 1,mm + do j = 1,mm + ha(j,i) = hmat(j,i) + end do + end do + call least_squares_decompose(nfreq,ha) + + DO N=1,MNP + +! At each node transfer the RHS vector and solve + do k=1,mm + hap(k) = GLOELV(k,n) + end do + call least_squares_solve(nfreq,ha,hap,hax) + +! Compute amplitude and phase for each constituent making sure that the +! phase is between 0 and 360 deg. + do i=1,nfreq + i1=2*i-1 + i2=i1+1 + emag(i)=sqrt(hax(i1)*hax(i1)+hax(i2)*hax(i2)) + emagt(i,n)=emag(i)/haff(i) + if((hax(i1).eq.0.0_RKIND).and.(hax(i2).eq.0.0_RKIND)) then + phasee(i)=0.0_RKIND + else + phasee(i) = atan2(hax(i2),hax(i1)) + endif + phaseden(i,n)=phasee(i)+haface(i) + phaseden(i,n)=convrd*phaseden(i,n) + if(phaseden(i,n).lt.0.0_RKIND) phaseden(i,n)=phaseden(i,n)+360.0_RKIND + if(phaseden(i,n).ge.360.0_RKIND) phaseden(i,n)=phaseden(i,n)-360.0_RKIND + end do + + end do + + return + end subroutine harmonic_analysis_solve !}}} + + +!*********************************************************************** +! +! routine least_squares_decompose +! +!> \brief Fill out symmetric matrix and decompose +!> \author Steven Brus +!> \date July 2020 +!> \details +!> Adapted from the ADCIRC code by R.A. Luettich and J.J. Westerink +! +!----------------------------------------------------------------------- + + subroutine least_squares_decompose(nfreq,ha) !{{{ + + implicit none + + integer, intent(in) :: nfreq ! number of analysis constituents + real(kind=RKIND), dimension(:,:), intent(inout) :: ha ! LHS matrix + + integer i,j,ir,ire,k,jr,mm + + mm = 2*nfreq + +! Set up the lower triangular part of the LHS matrix + + do j=1,mm + do i=j,mm + ha(i,j)=ha(j,i) + end do + end do + +! Decomposition of matrix + + do 100 ir=1,mm + ire=ir+1 + do 20 j=ire,mm + 20 ha(ir,j)=ha(ir,j)/ha(ir,ir) + if (ire.gt.mm) goto 100 + do 40 j=ire,mm + do 40 k=ire,mm + 40 ha(k,j)=ha(k,j)-ha(k,ir)*ha(ir,j) + do 50 j=ire,mm + 50 ha(j,ir)=0.0_RKIND + 100 continue + + + end subroutine least_squares_decompose !}}} + +!*********************************************************************** +! +! routine least_squares_solve +! +!> \brief Solves system a*x=b by l*d*l(tr) decomp in full storage mode +!> \author Steven Brus +!> \date July 2020 +!> \details +!> Adapted from the ADCIRC code by R.A. Luettich and J.J. Westerink +! +!----------------------------------------------------------------------- + + subroutine least_squares_solve(nfreq,ha,hap,hax) !{{{ + + implicit none + + integer, intent(in) :: nfreq ! number of analysis constituents + real(kind=RKIND), dimension(:,:), intent(in) :: ha ! LHS matrix + real(kind=RKIND), dimension(:), intent(in) :: hap ! RHS vector + real(kind=RKIND), dimension(:), intent(out) :: hax ! solution vector + + integer idecom,i,j,ir,ire,k,jr,mm + real(kind=RKIND),allocatable :: c(:),y(:) + + mm = 2*nfreq + + allocate ( c(mm),y(mm) ) + +! solve for y by forward substitution for l*y=p + + do 120 ir=1,mm + y(ir)=hap(ir) + do 110 jr=1,ir-1 + 110 y(ir)=y(ir)-ha(jr,ir)*y(jr) + 120 continue + +! calculate c=d**(-1)*y + + do 130 ir=1,mm + 130 c(ir)=y(ir)/ha(ir,ir) + +! solve for x by back-substituting for l(tr)*x=c + + ir=mm + 140 continue + hax(ir)=c(ir) + do 150 jr=ir+1,mm + 150 hax(ir)=hax(ir)-ha(ir,jr)*hax(jr) + ir=ir-1 + if(ir.ge.1) goto 140 + + end subroutine least_squares_solve!}}} + +end module ocn_harmonic_analysis + +! vim: foldmethod=marker diff --git a/src/core_ocean/analysis_members/mpas_ocn_sediment_flux_index.F b/src/core_ocean/analysis_members/mpas_ocn_sediment_flux_index.F new file mode 100644 index 0000000000..777a6fc00d --- /dev/null +++ b/src/core_ocean/analysis_members/mpas_ocn_sediment_flux_index.F @@ -0,0 +1,318 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_sediment_flux_index +! +!> \brief MPAS ocean analysis mode member: sediment_flux_index +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/02/26 +!> \details +!> MPAS ocean analysis mode member: sediment_flux_index +!> +!----------------------------------------------------------------------- + +module ocn_sediment_flux_index + + use mpas_derived_types + use mpas_pool_routines + use mpas_dmpar + use mpas_timekeeping + use mpas_stream_manager + + use ocn_constants + use ocn_diagnostics_routines + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_init_sediment_flux_index, & + ocn_compute_sediment_flux_index, & + ocn_restart_sediment_flux_index, & + ocn_finalize_sediment_flux_index + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_init_sediment_flux_index +! +!> \brief Initialize MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/02/26 +!> \details +!> This routine conducts all initializations required for the +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_init_sediment_flux_index(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_init_sediment_flux_index!}}} + +!*********************************************************************** +! +! routine ocn_compute_sediment_flux_index +! +!> \brief Compute MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/02/26 +!> \details +!> This routine conducts all computation required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_compute_sediment_flux_index(domain, timeLevel, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + integer, intent(in) :: timeLevel + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), pointer :: sedimentFluxIndexAMPool + type (dm_info) :: dminfo + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: diagnosticsPool + + real (kind=RKIND), dimension(:,:), pointer :: velX, velY, velZ + real (kind=RKIND), dimension(:), pointer :: posX, posY + real (kind=RKIND), dimension(:), pointer :: sfiVAx, sfiVAy, sfiBx, sfiBy + logical, pointer :: on_a_sphere, use_lat_lon_coords + integer, pointer :: nCells, nVertLevels, nCellsSolve + integer k, iCell, i + err = 0 + + dminfo = domain % dminfo + + block => domain % blocklist + do while (associated(block)) + ! get dimensions + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels) + + ! get pointers to pools + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) + call mpas_pool_get_subpool(block % structs, 'sedimentFluxIndexAM', sedimentFluxIndexAMPool) + + call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentFluxIndex_use_lat_lon_coords', use_lat_lon_coords) + + if (.not. on_a_sphere) then + use_lat_lon_coords = .false. + end if + + if (use_lat_lon_coords) then + call mpas_pool_get_array(meshPool, 'lonCell', posX) + call mpas_pool_get_array(meshPool, 'latCell', posY) + call mpas_pool_get_array(diagnosticsPool, 'velocityZonal', velX) + call mpas_pool_get_array(diagnosticsPool, 'velocityMeridional', velY) + else + call mpas_pool_get_array(meshPool, 'xCell', posX) + call mpas_pool_get_array(meshPool, 'yCell', posY) + call mpas_pool_get_array(diagnosticsPool, 'velocityX', velX) + call mpas_pool_get_array(diagnosticsPool, 'velocityY', velY) + end if + + call mpas_pool_get_array(sedimentFluxIndexAMPool, 'sedimentFluxIndexVAX', sfiVAx) + call mpas_pool_get_array(sedimentFluxIndexAMPool, 'sedimentFluxIndexVAY', sfiVAy) + call mpas_pool_get_array(sedimentFluxIndexAMPool, 'sedimentFluxIndexBX', sfiBx) + call mpas_pool_get_array(sedimentFluxIndexAMPool, 'sedimentFluxIndexBY', sfiBy) + ! Computations which are functions of nCells, nEdges, or nVertices + ! must be placed within this block loop + ! Here are some example loops + do iCell = 1,nCellsSolve + sfiVAx(iCell) = (sum(velX(:,iCell))/float(nVertLevels))**3.0_RKIND + sfiVAy(iCell) = (sum(velY(:,iCell))/float(nVertLevels))**3.0_RKIND + sfiBx(iCell) = velX(1,iCell)**3.0_RKIND + sfiBy(iCell) = velY(1,iCell)**3.0_RKIND + end do + + block => block % next + end do + + end subroutine ocn_compute_sediment_flux_index!}}} + +!*********************************************************************** +! +! routine ocn_restart_sediment_flux_index +! +!> \brief Save restart for MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/02/26 +!> \details +!> This routine conducts computation required to save a restart state +!> for the MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_restart_sediment_flux_index(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_restart_sediment_flux_index!}}} + +!*********************************************************************** +! +! routine ocn_finalize_sediment_flux_index +! +!> \brief Finalize MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/02/26 +!> \details +!> This routine conducts all finalizations required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_finalize_sediment_flux_index(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_finalize_sediment_flux_index!}}} + +end module ocn_sediment_flux_index + +! vim: foldmethod=marker diff --git a/src/core_ocean/analysis_members/mpas_ocn_sediment_transport.F b/src/core_ocean/analysis_members/mpas_ocn_sediment_transport.F new file mode 100644 index 0000000000..4804684380 --- /dev/null +++ b/src/core_ocean/analysis_members/mpas_ocn_sediment_transport.F @@ -0,0 +1,545 @@ +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! and the University Corporation for Atmospheric Research (UCAR). +! +! Unless noted otherwise source code is licensed under the BSD license. +! Additional copyright and license information can be found in the LICENSE file +! distributed with this code, or at http://mpas-dev.github.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! ocn_sediment_transport +! +!> \brief MPAS ocean analysis mode member: sediment_transport +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/03/07 +!> \details +!> MPAS ocean analysis mode member: sediment_transport +!----------------------------------------------------------------------- + +module ocn_sediment_transport + + use mpas_derived_types + use mpas_pool_routines + use mpas_dmpar + use mpas_timekeeping + use mpas_stream_manager + + use ocn_constants + use ocn_diagnostics_routines + + implicit none + private + save + + !-------------------------------------------------------------------- + ! + ! Public parameters + ! + !-------------------------------------------------------------------- + + !-------------------------------------------------------------------- + ! + ! Public member functions + ! + !-------------------------------------------------------------------- + + public :: ocn_init_sediment_transport, & + ocn_compute_sediment_transport, & + ocn_restart_sediment_transport, & + ocn_finalize_sediment_transport + + !-------------------------------------------------------------------- + ! + ! Private module variables + ! + !-------------------------------------------------------------------- + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine ocn_init_sediment_transport +! +!> \brief Initialize MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/03/07 +!> \details +!> This routine conducts all initializations required for the +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_init_sediment_transport(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_init_sediment_transport!}}} + +!*********************************************************************** +! +! routine ocn_compute_sediment_transport +! +!> \brief Compute MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/03/07 +!> \details +!> This routine conducts all computation required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_compute_sediment_transport(domain, timeLevel, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + integer, intent(in) :: timeLevel + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), pointer :: sedimentTransportAMPool + type (dm_info) :: dminfo + type (block_type), pointer :: block + type (mpas_pool_type), pointer :: statePool + type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: diagnosticsPool + + real (kind=RKIND), pointer :: salpha, SD50, rho0, rhoS, Vh, Dv,Ws + real (kind=RKIND), pointer :: tau_ce, Erate, poro, tau_cd, Cd, Manning + real (kind=RKIND), dimension(:,:), pointer :: velX, velY, velZ + real (kind=RKIND), dimension(:), pointer :: posX, posY, bottomDepth, ssh + real (kind=RKIND), dimension(:), pointer :: sedFluxVAx, sedFluxVAy, sedFluxBx, sedFluxBy + real (kind=RKIND), dimension(:), pointer :: ero_flux, dep_flux + real (kind=RKIND), dimension(:), pointer :: bedld_x, bedld_y + real (kind=RKIND), dimension(:), pointer :: SSC_ref + real (kind=RKIND), dimension(:,:), pointer :: SSC + real (kind=RKIND), dimension(:,:), pointer :: zMid + + logical, pointer :: on_a_sphere, use_lat_lon_coords + logical, pointer :: bedload, suspended + integer, pointer :: nCells, nVertLevels, nCellsSolve + character (len=StrKIND), pointer :: ws_formula, bedld_formula + character (len=StrKIND), pointer :: SSC_ref_formula + + integer k, iCell, i + real ratio1, rho_R, ND50, Umag, Umagb, dstar, Chezy, Usf, RouseN + real cff1, cff2, cff3, cff4, cff5 + real phicw, w_asym !! parameters related to waves, assigned ZEROs for just now + real tau_tide, tau_wave, tau_mean, theta_mean, theta_wave, theta_ce,theta_sf,za + real velXb, velYb, velXm, velYm + real phi_x1, phi_x2, phi_x, phi_y + real, parameter :: g = 9.81_RKIND + real, parameter :: eps = 1.0E-14_RKIND + real zCell + err = 0 + + dminfo = domain % dminfo + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_alpha',salpha) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_grain_size',SD50) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_drag_coefficient',Cd) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_grain_porosity',poro) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_erate',Erate) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_tau_ce',tau_ce) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_tau_cd',tau_cd) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_Manning_coef',Manning) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_water_density',rho0) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_grain_density',rhoS) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_kinematic_viscosity',Vh) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_vertical_diffusion_coefficient',Dv) + + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_use_lat_lon_coords', use_lat_lon_coords) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_bedload', bedload) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_suspended', suspended) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_ws_formula', ws_formula) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_bedld_formula', bedld_formula) + call mpas_pool_get_config(ocnConfigs, 'config_AM_sedimentTransport_SSC_ref_formula', SSC_ref_formula) + + + block => domain % blocklist + do while (associated(block)) + ! get dimensions + call mpas_pool_get_dimension(block % dimensions, 'nCells', nCells) + call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve) + call mpas_pool_get_dimension(block % dimensions, 'nVertLevels', nVertLevels) + + ! get pointers to pools + call mpas_pool_get_subpool(block % structs, 'state', statePool) + call mpas_pool_get_subpool(block % structs, 'mesh', meshPool) + call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool) + call mpas_pool_get_subpool(block % structs, 'sedimentTransportAM', sedimentTransportAMPool) + + call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid) + call mpas_pool_get_array(statePool, 'ssh', ssh, 1) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) + call mpas_pool_get_config(meshPool, 'on_a_sphere', on_a_sphere) + if (.not. on_a_sphere) then + use_lat_lon_coords = .false. + end if + + if (use_lat_lon_coords) then + call mpas_pool_get_array(meshPool, 'lonCell', posX) + call mpas_pool_get_array(meshPool, 'latCell', posY) + call mpas_pool_get_array(diagnosticsPool, 'velocityZonal', velX) + call mpas_pool_get_array(diagnosticsPool, 'velocityMeridional', velY) + else + call mpas_pool_get_array(meshPool, 'xCell', posX) + call mpas_pool_get_array(meshPool, 'yCell', posY) + call mpas_pool_get_array(diagnosticsPool, 'velocityX', velX) + call mpas_pool_get_array(diagnosticsPool, 'velocityY', velY) + end if + + + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentFallVelocity', Ws) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentErosionFlux', ero_flux) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentDepositionFlux', dep_flux) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentFluxVAX', sedFluxVAx) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentFluxVAY', sedFluxVAy) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentFluxBX', sedFluxBx) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentFluxBY', sedFluxBy) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentBedloadX', bedld_x) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentBedloadY', bedld_y) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentBottomReferenceConcentration', SSC_ref) + call mpas_pool_get_array(sedimentTransportAMPool, 'sedimentConcentration', SSC) + +!======> COMPUTE SedimentFallVelocity [m/s] + ! sieve diameter (SD50) to nominal diameter (ND50), Raudkivi [1990] + ND50 = SD50*1.1_RKIND + ! sediment/water density comparsion -> specific graivity + rho_R = rhoS/rho0-1 + if (ws_formula == 'Goldstein-Coco2013') then + !Goldstein & Coco [2013] + Ws = (37.8_RKIND*ND50*rho_R*(1+100.0_RKIND*ND50))/ & + (0.383_RKIND+1E4*rho_R*Vh+1E2*ND50*rho_R**2) + elseif (ws_formula == 'Cheng1997') then + dstar = SD50* (rho_R*g/Vh**2.0_RKIND)**0.333_RKIND + Ws = Vh/SD50*(sqrt(25.0_RKIND+1.2_RKIND*dstar**2)-5.0_RKIND)**1.5_RKIND + elseif (ws_formula == 'VanRijn1993') then + if (SD50<1.0E-4) then + Ws = rho_R*g*SD50**2.0_RKIND/(18.0_RKIND*Vh) + elseif (SD50>1.0E-3) then + Ws = 1.1_RKIND*sqrt(rho_R*g*SD50) + else + Ws = 10.0_RKIND*Vh/SD50 *(sqrt(1.0_RKIND+0.01_RKIND*rho_R*g*SD50**3.0_RKIND/Vh**2.0_RKIND)-1) + end if + elseif (ws_formula == 'Soulsby1997') then + Ws = 10.36_RKIND*Vh/SD50 * (sqrt(1.0_RKIND+0.156_RKIND*rho_R*g*SD50**3.0_RKIND/(16.0_RKIND*Vh**2.0))-1) + else + print*,'please select the ws_formula from one of the following:' + print*, 'VanRijn1993, Soulsby1997, Cheng1997, Goldstein-Coco2013' + print*, 'Model will stop ...' + exit + end if +!<====== + +!======> COMPUTE SedimentTransportFlux [kg/m/s] + ! compute ratio1 in sedFlux = ratio1*U^3, Grawe et al. [2014] + ratio1 = salpha*Dv/Ws**2 + + do iCell = 1,nCellsSolve + velXm = sum(velX(:,iCell))/float(nVertLevels) + velYm = sum(velY(:,iCell))/float(nVertLevels) + velXb = velX(size(velX,1),iCell) + velYb = velY(size(velY,1),iCell) + + sedFluxVAx(iCell) = ratio1*velXm**3.0_RKIND + sedFluxVAy(iCell) = ratio1*velYm**3.0_RKIND + sedFluxBx(iCell) = ratio1*velXb**3.0_RKIND + sedFluxBy(iCell) = ratio1*velYb**3.0_RKIND + end do +!<====== + + cff1 = rho_R*g*SD50+eps + cff2 = sqrt(rho_R*g*SD50)*SD50*rhoS + !wave-related parameters + w_asym = 0.0_RKIND + phicw = 0.0_RKIND + + do iCell = 1,nCellsSolve + velXm = sum(velX(:,iCell))/float(nVertLevels) + velYm = sum(velY(:,iCell))/float(nVertLevels) + velXb = velX(size(velX,1),iCell) + velYb = velY(size(velY,1),iCell) + Umag = sqrt(velXm**2.0_RKIND + velYm**2.0_RKIND) + Umagb = sqrt(velXb**2.0_RKIND + velYb**2.0_RKIND) + + ! compute bottom friction Cd = [von_Karmann/(1+log(z0/h))]^2, von_karmann=0.4, z0=SD50/12, h=bottomDepth+ssh + cff3 = SD50/12.0_RKIND + cff4 = 1+log(cff3/(eps+ssh(iCell)+bottomDepth(iCell))) + Cd = (0.4_RKIND/cff4)**2.0_RKIND + ! compute shear velocity, tau= rho0*Cd*Umag**2 = rho0*Usf**2 => Usf=sqrt(Cd)*Umag + Usf = Umag * sqrt(Cd) + ! compute Rouse number, RouseN=Ws/(von_karmann*Usf) + RouseN = Ws/(0.4_RKIND*(Usf+eps)) +!! Rouse number (suspension parameter) indicates the balance between sediment suspension and settling. According to the +!! equation of it, when Ws>Usf, which also means RouseN>2.5, there should be no sediment in suspension. To avoid the +!! "Float-point exception" error when Usf is very small and RouseN is super huge, we mannually set a maximum value of +!! 25 for RouseN, which is 10 times bigger than RouseN with ideally no suspended sediment. + RouseN = min(RouseN, 25.0_RKIND) + + ! compute shear stress [m^2 s^{-2}] + tau_tide = Cd*Umagb**2.0_RKIND + tau_wave = 0.0_RKIND + tau_mean = tau_tide*(1.0_RKIND+1.2_RKIND*(tau_wave/(tau_wave+tau_tide+eps))**1.5_RKIND) + +!======> COMPUTE Suspended Transport [] <--- to be continued + if (suspended) then + ! surface erosion mass flux ero_flux in unit (kg m-2 s-1) + ero_flux(iCell) = (1.0_RKIND-poro)*Erate*MAX(tau_mean/tau_ce -1.0_RKIND, 0.0_RKIND) + end if +!<====== + +!======> COMPUTE SSC_ref: near-bottom suspended sediment concentration (reference concentration) [kg/m3] + ! Three equations: + ! 1) Lee2004,JGR; 2) Goldstein2014,Earth Surface Dynamics; 3) Zyserman-Fredsoe1994, J. Hydraul. Eng. + if (SSC_ref_formula == 'Lee2004') then + theta_sf = Usf**2.0_RKIND/cff1 + SSC_ref(iCell) = 2.58_RKIND*(theta_sf*Usf/Ws)**1.45_RKIND + elseif (SSC_ref_formula == 'Goldstein2014') then + SSC_ref(iCell) = (0.328_RKIND*Umagb/(0.0688_RKIND+1E3*SD50))**2_RKIND + elseif (SSC_ref_formula == 'Zyserman-Fredsoe1994') then + theta_sf = tau_mean/cff1 + theta_sf = max(theta_sf, 0.045_RKIND) + cff3 = 0.331_RKIND*(theta_sf-0.045_RKIND)**1.75_RKIND + cff4 = 1.0_RKIND+0.72_RKIND*(theta_sf-0.045_RKIND)**1.75_RKIND + SSC_ref(iCell) = cff3/cff4 + else + print*, 'Please pick one SSC_ref_formula from the following three:' + print*, 'Lee2004, Goldstein2014,Zyserman-Fredsoe1994' + print*, 'Model will stop ...' + exit + endif +!<====== + +!======> Compute SSC based on the Rouse Profile + ! SSC(z) = SSC_ref*[z/za * (h-za)/(h-z)]**(-b), where b is the Rouse number (or suspension parameter) + ! b = Ws/(von_karmann*usf); za is reference height (0.01m), h is bottomDepth + ! zMid in MPAS-O is the distance under free surface and is negative most of the time + ! In the above equation, z means the depth above the bottom, so z=h+zMid --> h-z = -zMid +!! Since the reference concentration is the largest SSC at height za above bed, we need make sure h+zMid>=za +!! that is to say zMid = max(zMid,za-h) + za = 0.01_RKIND + if (SSC_ref_formula == 'Zyserman-Fredsoe1994') za=2.0_RKIND*SD50 + do k=1,nVertLevels + zMid(k,iCell) = max(za-bottomDepth(iCell), zMid(k,iCell)) ! make sure -zMid>=za + zCell = bottomDepth(iCell)+zMid(k,iCell)+eps + cff3 = za/(bottomDepth(iCell)-za+eps) + if (zCell/bottomDepth(iCell) .lt. 0.5_RKIND) then + cff4 = -zMid(k,iCell)/zCell*cff3 + SSC(k,iCell) = SSC_ref(iCell)*cff4**RouseN + else + cff4 = exp(-4.0_RKIND*RouseN*(zCell/bottomDepth(iCell)-0.5_RKIND)) + SSC(k,iCell) = SSC_ref(iCell)*cff4*cff3**RouseN + endif + enddo !! k loop +!<====== + +!======> COMPUTE Bedload Transport [kg/m/s], using one of the three formulae define in bedld_formula + if (bedload) then + theta_mean = tau_mean/cff1 + theta_wave = tau_wave/cff1 + theta_ce = tau_ce/cff1 + + if (bedld_formula == 'Soulsby-Damgaard')then + cff3 = 0.5*(1.0_RKIND+SIGN(1.0_RKIND,theta_mean/theta_ce-1.0_RKIND)) + + phi_x1 = 12.0_RKIND*sqrt(theta_mean)*MAX(theta_mean-theta_ce,0.0_RKIND) + phi_x2 = 12.0_RKIND*(0.9534_RKIND+0.1907_RKIND*COS(2.0_RKIND*phicw))* & + sqrt(theta_wave)*theta_mean + & + 12.0_RKIND*(0.229_RKIND*w_asym*theta_wave**1.5_RKIND*COS(phicw)) + + if (ABS(phi_x2) .gt. phi_x1) then + phi_x = phi_x2 + else + phi_x = phi_x1 + end if + + bedld_x(iCell) = phi_x*cff2*cff3 + + phi_y = 12.0_RKIND*0.1907_RKIND*theta_wave**2* & + (theta_mean*SIN(2.0_RKIND*phicw)+1.2_RKIND*w_asym*theta_wave*SIN(phicw)) & + /(theta_wave**1.5_RKIND + 1.5_RKIND*theta_mean**1.5_RKIND+eps) + bedld_y = phi_y*cff2*cff3 + + elseif (bedld_formula == 'Meyer-Peter-Meuller') then + cff3 = 1.0_RKIND/sqrt(velXb**4.0_RKIND+velYb**4.0_RKIND+eps) + phi_x = max(8.0_RKIND*(theta_mean-0.047_RKIND), 0.0_RKIND)*velXb**2.0_RKIND + phi_y = max(8.0_RKIND*(theta_mean-0.047_RKIND), 0.0_RKIND)*velYb**2.0_RKIND + + bedld_x(iCell) = phi_x*cff2*cff3 + bedld_y(iCell) = phi_y*cff2*cff3 + + elseif (bedld_formula == 'Engelund-Hansen') then + ! Chezy = 1/Manning * water_depth^(1/6) with unit [m^0.5 s^(-1)] + ! In ocean water, Manning coeffcient varies between [0.012 0.025],see P.C. Kerr et al, 2013 JGR-Ocean + Chezy = 1.0_RKIND/Manning*abs(ssh(iCell)+bottomDepth(iCell))**(1.0_RKIND/6.0_RKIND) + Chezy = max(20.0_RKIND, Chezy) + Chezy = min(100.0_RKIND, Chezy) + cff4 = 0.05_RKIND*Umag**4.0_RKIND*rhoS + cff5 = sqrt(g)*Chezy**3.0_RKIND*rho_R**2.0_RKIND*SD50 + bedld_x(iCell) = cff4/cff5*velXb + bedld_y(iCell) = cff4/cff5*velYb + else + print*, 'pick one bedld_formula from the following three:' + print*, 'Soulsby-Damgaard, Meyer-Peter-Meuller, Engelund-Hansen' + print*, 'Model will stop ...' + exit + end if + end if + + end do !! iCell loop + + block => block % next + end do !! do while (associated(block)) loop + + end subroutine ocn_compute_sediment_transport!}}} + +!*********************************************************************** +! +! routine ocn_restart_sediment_transport +! +!> \brief Save restart for MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/03/07 +!> \details +!> This routine conducts computation required to save a restart state +!> for the MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_restart_sediment_transport(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_restart_sediment_transport!}}} + +!*********************************************************************** +! +! routine ocn_finalize_sediment_transport +! +!> \brief Finalize MPAS-Ocean analysis member +!> \author Zhendong Cao and Phillip J. Wolfram +!> \date 2019/03/07 +!> \details +!> This routine conducts all finalizations required for this +!> MPAS-Ocean analysis member. +! +!----------------------------------------------------------------------- + + subroutine ocn_finalize_sediment_transport(domain, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + type (domain_type), intent(inout) :: domain + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + end subroutine ocn_finalize_sediment_transport!}}} + +end module ocn_sediment_transport + +! vim: foldmethod=marker diff --git a/src/core_ocean/driver/mpas_ocn_core_interface.F b/src/core_ocean/driver/mpas_ocn_core_interface.F index aeae86825d..dab9fe154a 100644 --- a/src/core_ocean/driver/mpas_ocn_core_interface.F +++ b/src/core_ocean/driver/mpas_ocn_core_interface.F @@ -113,6 +113,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: thicknessFilterActive logical, pointer :: splitTimeIntegratorActive logical, pointer :: windStressBulkPKGActive + logical, pointer :: variableBottomDragPKGActive logical, pointer :: tracerBudgetActive logical, pointer :: landIcePressurePKGActive logical, pointer :: landIceFluxesPKGActive @@ -121,6 +122,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: frazilIceActive logical, pointer :: tidalForcingActive logical, pointer :: tidalPotentialForcingPKGActive + logical, pointer :: vegetationDragPKGActive logical, pointer :: inSituEOSActive logical, pointer :: variableShortwaveActive logical, pointer :: gmActive @@ -152,6 +154,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_tidal_potential_forcing logical, pointer :: config_use_GM logical, pointer :: config_use_Redi + logical, pointer :: config_use_vegetation_drag logical, pointer :: config_use_time_varying_atmospheric_forcing logical, pointer :: config_use_time_varying_land_ice_forcing @@ -160,6 +163,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ character (len=StrKIND), pointer :: config_pressure_gradient_type character (len=StrKIND), pointer :: config_sw_absorption_type + logical, pointer :: config_use_variable_drag logical, pointer :: config_use_bulk_wind_stress logical, pointer :: config_use_bulk_thickness_flux logical, pointer :: config_compute_active_tracer_budgets @@ -231,6 +235,14 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ windStressBulkPKGActive = .true. end if + ! + ! test for variable bottom drag of momentum, variableBottomDragPKG + call mpas_pool_get_package(packagePool, 'variableBottomDragPKGActive', variableBottomDragPKGActive) + call mpas_pool_get_config(configPool, 'config_use_variable_drag', config_use_variable_drag) + if ( config_use_variable_drag ) then + variableBottomDragPKGActive = .true. + end if + ! ! test for tracer budget ! @@ -297,6 +309,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ tidalPotentialForcingPKGActive = .true. end if + ! + ! test for use of vegetation drag, vegetationDragPKGActive + ! + call mpas_pool_get_package(packagePool, 'vegetationDragPKGActive', vegetationDragPKGActive) + call mpas_pool_get_config(configPool, 'config_use_vegetation_drag', config_use_vegetation_drag) + if (config_use_vegetation_drag) then + vegetationDragPKGActive = .true. + end if + ! ! test for use of gm ! diff --git a/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F b/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F index d4c49a5fd2..0daa9a3cf1 100644 --- a/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F +++ b/src/core_ocean/mode_forward/mpas_ocn_forward_mode.F @@ -56,7 +56,7 @@ module ocn_forward_mode use ocn_surface_land_ice_fluxes use ocn_frazil_forcing use ocn_tidal_forcing - use ocn_tidal_potential_forcing + use ocn_vel_tidal_potential use ocn_tracer_hmix use ocn_tracer_hmix_redi @@ -210,7 +210,7 @@ function ocn_forward_mode_init(domain, startTimeStamp) result(ierr)!{{{ ierr = ior(ierr, err_tmp) call ocn_tidal_forcing_init(err_tmp) ierr = ior(ierr, err_tmp) - call ocn_tidal_potential_forcing_init(domain,err_tmp) + call ocn_vel_tidal_potential_init(domain,err_tmp) ierr = ior(ierr, err_tmp) call ocn_tracer_hmix_init(err_tmp) diff --git a/src/core_ocean/mode_init/Registry_hurricane.xml b/src/core_ocean/mode_init/Registry_hurricane.xml index 59c4ec54ae..a024826d1b 100644 --- a/src/core_ocean/mode_init/Registry_hurricane.xml +++ b/src/core_ocean/mode_init/Registry_hurricane.xml @@ -39,4 +39,24 @@ description="Amplitude of sea level rise." possible_values="Any real number." /> + + + + + diff --git a/src/core_ocean/mode_init/Registry_tidal_boundary.xml b/src/core_ocean/mode_init/Registry_tidal_boundary.xml index 5e2134cb7b..f1aacb6613 100644 --- a/src/core_ocean/mode_init/Registry_tidal_boundary.xml +++ b/src/core_ocean/mode_init/Registry_tidal_boundary.xml @@ -51,4 +51,60 @@ description="Fraction of the domain the plug should take up initially. Only in the y direction." possible_values="Any real number between 0 and 1." /> + + + + + + + + + + + + + + diff --git a/src/core_ocean/mode_init/mpas_ocn_init_hurricane.F b/src/core_ocean/mode_init/mpas_ocn_init_hurricane.F index 29867aa54c..8b265c69f1 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_hurricane.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_hurricane.F @@ -83,6 +83,7 @@ subroutine ocn_init_setup_hurricane(domain, iErr)!{{{ type (block_type), pointer :: block_ptr type (mpas_pool_type), pointer :: meshPool + type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: statePool type (mpas_pool_type), pointer :: tracersPool type (mpas_pool_type), pointer :: verticalMeshPool @@ -104,6 +105,7 @@ subroutine ocn_init_setup_hurricane(domain, iErr)!{{{ logical, pointer :: on_a_sphere integer, dimension(:), pointer :: maxLevelCell real (kind=RKIND), dimension(:), pointer :: ssh + real (kind=RKIND), dimension(:), pointer :: bottomDrag real (kind=RKIND), dimension(:), pointer :: refBottomDepth, refZMid, & vertCoordMovementWeights, bottomDepth, bottomDepthObserved, & fCell, fEdge, fVertex, & @@ -180,6 +182,7 @@ subroutine ocn_init_setup_hurricane(domain, iErr)!{{{ do while(associated(block_ptr)) call mpas_pool_get_subpool(block_ptr % structs, 'mesh', meshPool) call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool) + call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) call mpas_pool_get_subpool(block_ptr % structs, 'verticalMesh', verticalMeshPool) call mpas_pool_get_subpool(block_ptr % structs, 'diagnostics', diagnosticsPool) call mpas_pool_get_subpool(statePool, 'tracers', tracersPool) @@ -275,16 +278,8 @@ subroutine ocn_init_setup_hurricane(domain, iErr)!{{{ end if end do - ! shift coordinates to accomodate wetting and drying - do k = 1, nVertLevels - refBottomDepth(k) = refBottomDepth(k) - globalMinBottomDepth - end do - do iCell = 1, nCellsSolve - bottomDepth(iCell) = bottomDepth(iCell) - globalMinBottomDepth - ssh(iCell) = ssh(iCell) + globalMinBottomDepth - ! make sure depth is thick enough via ssh = TOTAL_DEPTH - bottomDepth ssh(iCell) = - bottomDepth(iCell) + & max(ssh(iCell) + bottomDepth(iCell), & @@ -363,6 +358,24 @@ subroutine ocn_init_setup_hurricane(domain, iErr)!{{{ fVertex(iVertex) = 2.0_RKIND * omega * sin(latVertex(iVertex)) end do + ! Set depth-variable drag + if (config_use_variable_drag) then + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + do iCell = 1, nCellsSolve + if (bottomDepth(iCell) <= config_hurricane_land_z_limit) then + bottomDrag(iCell) = config_hurricane_land_drag + else if (config_hurricane_land_z_limit < bottomDepth(iCell) .and. bottomDepth(iCell) < config_hurricane_marsh_z_limit) then + bottomDrag(iCell) = config_hurricane_marsh_drag + else if (config_hurricane_marsh_z_limit <= bottomDepth(iCell)) then + bottomDrag(iCell) = config_hurricane_channel_drag + else + call mpas_log_write('Default value for drag is not selected' & + // ' properly for bottomDepth of $r in cell $i!', & + realArgs=(/bottomDepth(iCell)/), intArgs=(/iCell/)) + end if + end do + end if + block_ptr => block_ptr % next end do diff --git a/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F b/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F index 15f81a80f6..9cd0016b2d 100644 --- a/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F +++ b/src/core_ocean/mode_init/mpas_ocn_init_tidal_boundary.F @@ -93,7 +93,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ type (mpas_pool_type), pointer :: verticalMeshPool type (mpas_pool_type), pointer :: tracersPool - integer :: iCell, k + integer :: iCell, k, N ! Define dimensions integer, pointer :: nCellsSolve, nEdgesSolve, nVertLevels, nVertLevelsP1 @@ -101,14 +101,17 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ ! Define arrays integer, dimension(:), pointer :: maxLevelCell + integer, dimension(:), pointer :: vegetationMask + real (kind=RKIND), dimension(:), pointer :: vegetationHeight,vegetationDensity, vegetationDiameter real (kind=RKIND), dimension(:), pointer :: yCell, refBottomDepth, bottomDepth, vertCoordMovementWeights, dcEdge real (kind=RKIND), dimension(:), pointer :: tidalInputMask + real (kind=RKIND), dimension(:), pointer :: bottomDrag real (kind=RKIND), dimension(:), pointer :: ssh real (kind=RKIND), dimension(:,:), pointer :: layerThickness, restingThickness, zMid real (kind=RKIND), dimension(:,:,:), pointer :: activeTracers, debugTracers - real (kind=RKIND), dimension(:), pointer :: interfaceLocations real (kind=RKIND), parameter :: eps=1.0e-12 + real (kind=RKIND) :: cff1,cff2,cff3, dep_mark1, dep_mark2 ! intermediate variables iErr = 0 @@ -190,9 +193,19 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ end if ! Set refBottomDepth, bottomDepth, and maxLevelCell - do k = 1, nVertLevels - refBottomDepth(k) = config_tidal_boundary_right_bottom_depth * interfaceLocations(k+1) - end do + if (min(config_tidal_boundary_left_bottom_depth, config_tidal_boundary_right_bottom_depth) < 0.0_RKIND) then + ! consider the case where there is wetting / drying and vertical mesh resolution is needed "on land" + do k = 1, nVertLevels + refBottomDepth(k) = config_tidal_boundary_left_bottom_depth + & + (config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth)* interfaceLocations(k+1) + end do + else + ! assumes we just need to build vertical mesh to deepest point (e.g., no "on land' vertical mesh) + do k = 1, nVertLevels + refBottomDepth(k) = \ + max(config_tidal_boundary_left_bottom_depth, config_tidal_boundary_right_bottom_depth) * interfaceLocations(k+1) + end do + end if if (config_tidal_boundary_use_distances) then yMin = config_tidal_boundary_left_value @@ -205,6 +218,78 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ (config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth) end do + if (config_use_variable_drag) then + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + do iCell = 1, nCellsSolve + bottomDrag(iCell) = config_tidal_forcing_left_Cd_or_n & + + (yCell(iCell) - yMin) / (yMax - yMin) * & + (config_tidal_forcing_right_Cd_or_n - config_tidal_forcing_left_Cd_or_n) + end do + end if + + if (config_use_vegetation_drag) then + call mpas_pool_get_array(forcingPool, 'vegetationMask', vegetationMask) + vegetationMask = 0 + call mpas_pool_get_array(forcingPool, 'vegetationDiameter', vegetationDiameter) + vegetationDiameter = config_idealized_vegetation_diameter + call mpas_pool_get_array(forcingPool, 'vegetationHeight', vegetationHeight) + vegetationHeight = config_idealized_vegetation_height + call mpas_pool_get_array(forcingPool, 'vegetationDensity', vegetationDensity) + vegetationDensity = config_idealized_vegetation_density + endif + + if (config_use_idealized_transect) then + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + config_idealized_transect_Lmarsh = 1.0_RKIND - config_idealized_transect_Lcoast & + - config_idealized_transect_Lshore + if (config_idealized_transect_Lmarsh .lt. 0.0_RKIND) then + call mpas_log_write("Lshore+Lcoast cannot be bigger than 1.0") + iErr = 1 + endif + cff1 = yMax*config_idealized_transect_Lcoast + cff2 = yMax*config_idealized_transect_Lmarsh + cff3 = yMax*config_idealized_transect_Lshore + ! by defining the slopes and left_bottom_depth, the pre-defined right_bottom_depth won't work. + ! Need redefine it. + config_tidal_boundary_right_bottom_depth = config_tidal_boundary_left_bottom_depth + & + cff1*config_idealized_transect_Scoast + & + cff2*config_idealized_transect_Smarsh + & + cff3*config_idealized_transect_Sshore + do iCell = 1, nCellsSolve + if (yCell(iCell) .lt. cff1) then + bottomDepth(iCell) = config_tidal_boundary_left_bottom_depth & + + (yCell(iCell)-yMin)*config_idealized_transect_Scoast + dep_mark1 = bottomDepth(iCell) + bottomDrag(iCell) = config_idealized_transect_roughness + elseif (yCell(iCell) .lt. (cff1+cff2)) then + bottomDepth(iCell) = dep_mark1 + (yCell(iCell)-cff1)*config_idealized_transect_Smarsh + dep_mark2 = bottomDepth(iCell) + bottomDrag(iCell) = config_idealized_transect_roughness_marsh + if (config_use_vegetation_drag) vegetationMask(iCell) = 1 + else + if (cff2 .eq. 0.0_RKIND) dep_mark2 = dep_mark1 + bottomDepth(iCell) = dep_mark2 + (yCell(iCell)-cff1-cff2)*config_idealized_transect_Sshore + bottomDrag(iCell) = config_idealized_transect_roughness + endif + end do !! do iCell + else + ! if not config_idealized_transect, assign a constant slope + do iCell = 1, nCellsSolve + bottomDepth(iCell) = config_tidal_boundary_left_bottom_depth & + + (yCell(iCell) - yMin) / (yMax - yMin) * & + (config_tidal_boundary_right_bottom_depth - config_tidal_boundary_left_bottom_depth) + end do + end if !! if config_idealized_transect + + if (config_tidal_boundary_right_bottom_depth < config_tidal_boundary_left_bottom_depth) then + call mpas_log_write('Right boundary must be deeper than left boundary!', MPAS_LOG_CRIT) + end if + + ! Set refBottomDepth, bottomDepth, and maxLevelCell + do k = 1, nVertLevels + refBottomDepth(k) = config_tidal_boundary_right_bottom_depth * interfaceLocations(k+1) + end do + if (config_use_wetting_drying .and. config_tidal_start_dry .and. & trim(config_tidal_boundary_layer_type) == 'zstar') then do iCell = 1, nCellsSolve @@ -281,6 +366,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ if (config_use_wetting_drying .and. config_tidal_start_dry) then do iCell = 1, nCellsSolve ssh(iCell) = -bottomDepth(iCell) + config_drying_min_cell_height*maxLevelCell(iCell) + ssh(iCell) = MAX(ssh(iCell),-config_tidal_forcing_monochromatic_baseline) ! also computes zMid do k = 1, maxLevelCell(iCell) layerThickness(k,iCell) = (ssh(iCell) + bottomDepth(iCell))/maxLevelCell(iCell) @@ -307,7 +393,7 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ ! Set tidal boundary mask do iCell = 1, nCellsSolve tidalInputMask(iCell) = 0.0_RKIND - if (yCell(iCell) > (25.0e3 - dcEdgeMinGlobal/2.0_RKIND)) then + if (yCell(iCell) > (yMax - dcEdgeMinGlobal/2.0_RKIND)) then tidalInputMask(iCell) = 1.0_RKIND ! spread it over multiple cells !if (yCell(iCell) > (25.0e3 - 3*dcEdgeMinGlobal)) then @@ -315,6 +401,11 @@ subroutine ocn_init_setup_tidal_boundary(domain, iErr)!{{{ end if end do + ! check that there is some tidalInputMask + if (.not. sum(tidalInputMask) > 0) then + call mpas_log_write('Input mask for tidal case is not set!', MPAS_LOG_CRIT) + end if + ! Set salinity if ( associated(activeTracers) ) then do iCell = 1, nCellsSolve diff --git a/src/core_ocean/shared/Makefile b/src/core_ocean/shared/Makefile index 1901178222..3db7fa8c08 100644 --- a/src/core_ocean/shared/Makefile +++ b/src/core_ocean/shared/Makefile @@ -66,13 +66,13 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_framework_forcing.o \ mpas_ocn_time_varying_forcing.o \ mpas_ocn_wetting_drying.o \ - mpas_ocn_tidal_potential_forcing.o + mpas_ocn_vel_tidal_potential.o all: $(OBJS) mpas_ocn_init_routines.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_forcing.o mpas_ocn_surface_land_ice_fluxes.o -mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_tidal_potential_forcing.o +mpas_ocn_tendency.o: mpas_ocn_high_freq_thickness_hmix_del2.o mpas_ocn_tracer_surface_restoring.o mpas_ocn_thick_surface_flux.o mpas_ocn_tracer_short_wave_absorption.o mpas_ocn_tracer_advection.o mpas_ocn_tracer_hmix.o mpas_ocn_tracer_nonlocalflux.o mpas_ocn_surface_bulk_forcing.o mpas_ocn_surface_land_ice_fluxes.o mpas_ocn_tracer_surface_flux_to_tend.o mpas_ocn_tracer_interior_restoring.o mpas_ocn_tracer_exponential_decay.o mpas_ocn_tracer_ideal_age.o mpas_ocn_tracer_TTD.o mpas_ocn_vmix.o mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_frazil_forcing.o mpas_ocn_tidal_forcing.o mpas_ocn_tracer_ecosys.o mpas_ocn_tracer_DMS.o mpas_ocn_tracer_MacroMolecules.o mpas_ocn_diagnostics.o mpas_ocn_wetting_drying.o mpas_ocn_tidal_potential_forcing.o mpas_ocn_vel_tidal_potential.o mpas_ocn_diagnostics_routines.o: mpas_ocn_constants.o mpas_ocn_config.o @@ -198,6 +198,8 @@ mpas_ocn_wetting_drying.o: mpas_ocn_diagnostics.o mpas_ocn_gm.o mpas_ocn_tidal_potential_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o +mpas_ocn_vel_pressure_grad.o: mpas_ocn_constants.o + clean: $(RM) *.o *.i *.mod *.f90 diff --git a/src/core_ocean/shared/mpas_ocn_tendency.F b/src/core_ocean/shared/mpas_ocn_tendency.F index 1f129e6e0f..7b6977ec54 100644 --- a/src/core_ocean/shared/mpas_ocn_tendency.F +++ b/src/core_ocean/shared/mpas_ocn_tendency.F @@ -60,7 +60,7 @@ module ocn_tendency use ocn_vel_forcing use ocn_vmix use ocn_wetting_drying - use ocn_tidal_potential_forcing + use ocn_vel_tidal_potential implicit none private @@ -321,31 +321,6 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP ! call ocn_vel_vadv_tend(meshPool, normalVelocity, layerThicknessEdge, vertAleTransportTop, tend_normalVelocity, err) - ! - ! velocity tendency: tidal potential (if needed) - ! For RK4, subtract the tidal potential from the zMid array and store in a work array, - ! Then point zMid to work array so the tidal potential terms are included inside the grad computed in ocn_vel_pressure_grad_tend - ! zMid is pointed back to the typical value afterward, so the work array is not used in place of zMid elsewhere. - ! - if (config_use_tidal_potential_forcing) then - call ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnosticsPool, err) - end if - - if (config_use_tidal_potential_forcing .and. config_time_integrator == 'RK4') then - call mpas_pool_get_array(forcingPool, 'tidalPotentialZMid', tidalPotentialZMid) - call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', tidalPotentialEta) - call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) - - do iCell = 1, nCells - do k = 1, maxLevelCell(iCell) - tidalPotentialZMid(k,iCell) = zMid(k,iCell) - tidalPotentialEta(iCell) & - - config_self_attraction_and_loading_beta * zMid(k,iCell) - end do - end do - - call mpas_pool_get_array(forcingPool, 'tidalPotentialZMid', zMid) - end if - ! ! velocity tendency: pressure gradient ! @@ -362,10 +337,14 @@ subroutine ocn_tend_vel(tendPool, statePool, forcingPool, diagnosticsPool, meshP inSituThermalExpansionCoeff,inSituSalineContractionCoeff) endif - if (config_use_tidal_potential_forcing .and. config_time_integrator == 'RK4') then - ! point zMid back to usual array, to be safe - call mpas_pool_get_array(diagnosticsPool, 'zMid', zMid) - end if + ! + ! velocity tendency: tidal potential (if needed) + ! + call ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnosticsPool, err) + if (config_time_integrator == 'RK4') then + ! for split explicit, tidal forcing is added in barotropic subsycles + call ocn_vel_tidal_potential_tend(meshPool, forcingPool, ssh, tend_normalVelocity, err) + endif ! ! velocity tendency: del2 dissipation, \nu_2 \nabla^2 u diff --git a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F index f6fa320581..c5d80384b1 100644 --- a/src/core_ocean/shared/mpas_ocn_tidal_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_tidal_forcing.F @@ -241,7 +241,8 @@ subroutine ocn_tidal_forcing_build_array(domain, meshPool, forcingPool, diagnost ! compute the tidalHeight if (trim(config_tidal_forcing_model) == 'monochromatic') then tidalHeight = config_tidal_forcing_monochromatic_amp * & - SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim) - & + SIN(2.0_RKIND*pii/config_tidal_forcing_monochromatic_period * daysSinceStartOfSim - & + pii*config_tidal_forcing_monochromatic_phaseLag/180.0_RKIND) - & config_tidal_forcing_monochromatic_baseline !else if (trim(config_tidal_forcing_type) == 'data') then ! ! data option diff --git a/src/core_ocean/shared/mpas_ocn_time_varying_forcing.F b/src/core_ocean/shared/mpas_ocn_time_varying_forcing.F index ac0bf333a1..8009506862 100644 --- a/src/core_ocean/shared/mpas_ocn_time_varying_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_time_varying_forcing.F @@ -364,7 +364,12 @@ subroutine post_atmospheric_forcing(block)!{{{ rhoAir = 1.225_RKIND windStressCoefficientLimit = 0.0035_RKIND - ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/config_time_varying_atmospheric_forcing_ramp) + if (daysSinceStartOfSim >= config_time_varying_atmospheric_forcing_ramp_delay) then + ramp = tanh((2.0_RKIND*(daysSinceStartOfSim-config_time_varying_atmospheric_forcing_ramp_delay)) & + /config_time_varying_atmospheric_forcing_ramp) + else + ramp = 0.0_RKIND + end if !$omp parallel !$omp do schedule(runtime) private(windStressCoefficient) diff --git a/src/core_ocean/shared/mpas_ocn_tidal_potential_forcing.F b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F similarity index 80% rename from src/core_ocean/shared/mpas_ocn_tidal_potential_forcing.F rename to src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F index 7cbeabc61c..b6c5bb2115 100644 --- a/src/core_ocean/shared/mpas_ocn_tidal_potential_forcing.F +++ b/src/core_ocean/shared/mpas_ocn_vel_tidal_potential.F @@ -7,7 +7,7 @@ ! !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! -! ocn_tidal_potential_forcing +! ocn_vel_tidal_potential ! !> \brief MPAS ocean tidal potential forcing module !> \author Steven Brus @@ -19,7 +19,7 @@ ! !----------------------------------------------------------------------- -module ocn_tidal_potential_forcing +module ocn_vel_tidal_potential use mpas_kind_types use mpas_constants @@ -28,7 +28,6 @@ module ocn_tidal_potential_forcing use mpas_timekeeping use mpas_timer use ocn_constants - use ocn_config implicit none private @@ -46,8 +45,10 @@ module ocn_tidal_potential_forcing ! !-------------------------------------------------------------------- - public :: ocn_compute_tidal_potential_forcing, & - ocn_tidal_potential_forcing_init + public :: ocn_vel_tidal_potential_tend, & + ocn_compute_tidal_potential_forcing, & + ocn_vel_tidal_potential_init, & + tidal_constituent_factors !-------------------------------------------------------------------- ! @@ -60,7 +61,7 @@ module ocn_tidal_potential_forcing character(:), allocatable :: constituent end type type(char_array), dimension(37) :: constituentList - + public :: char_array !*********************************************************************** @@ -68,7 +69,101 @@ module ocn_tidal_potential_forcing !*********************************************************************** ! -! routine ocn_tidal_potential_forcing_layer_thickness +! routine ocn_vel_tidal_potential_tend +! +!> \brief Computes tendency term for tidal potential +!> \author Steven Brus +!> \date April 2020 +!> \details +!> This routine computes the tidal potential tendency for momentum +!> based on current state. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_tidal_potential_tend(meshPool, forcingPool, ssh, tend, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:), intent(in) :: & + ssh !< Input: Sea surface height + + type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information + type (mpas_pool_type), intent(in) :: forcingPool !< Input: forcinginformation + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k, cell1, cell2, nEdges + integer, dimension(:), pointer :: nEdgesArray + integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell + integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask + + real (kind=RKIND), dimension(:), pointer :: dcEdge + real (kind=RKIND), dimension(:), pointer :: tidalPotentialEta + real (kind=RKIND), pointer :: config_self_attraction_and_loading_beta + real (kind=RKIND) :: invdcEdge + real (kind=RKIND) :: potentialGrad + logical, pointer :: config_use_tidal_potential_forcing + + err = 0 + + if (.not. tidalPotentialOn) return + + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'dcEdge', dcEdge) + call mpas_pool_get_array(meshPool, 'edgeMask', edgeMask) + call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', tidalPotentialEta) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing', config_use_tidal_potential_forcing) + call mpas_pool_get_config(ocnConfigs, 'config_self_attraction_and_loading_beta', config_self_attraction_and_loading_beta) + + nEdges = nEdgesArray( 1 ) + + !$omp do schedule(runtime) private(cell1, cell2, invdcEdge, potentialGrad, k) + do iEdge=1,nEdges + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + invdcEdge = 1.0_RKIND / dcEdge(iEdge) + + potentialGrad = - gravity * invdcEdge * ( tidalPotentialEta(cell2) - tidalPotentialEta(cell1) & + + config_self_attraction_and_loading_beta*(ssh(cell2) - ssh(cell1))) + + do k=1,maxLevelEdgeTop(iEdge) + tend(k,iEdge) = tend(k,iEdge) - edgeMask(k,iEdge) * potentialGrad + end do + end do + !$omp end do + + end subroutine ocn_vel_tidal_potential_tend!}}} + +!*********************************************************************** +! +! routine ocn_compute_tidal_potential_forcing ! !> \brief Computes equilibrium tidal potential !> \author Steven Brus @@ -111,6 +206,7 @@ subroutine ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnostic !----------------------------------------------------------------- integer, pointer :: nCells + real (kind=RKIND), pointer :: tidalPotentialRamp real (kind=RKIND), dimension(:), pointer :: lonCell integer, dimension(:), pointer :: maxLevelCell @@ -132,6 +228,7 @@ subroutine ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnostic if ( .not. tidalPotentialOn ) return + call mpas_pool_get_config(ocnConfigs, 'config_tidal_potential_ramp', tidalPotentialRamp) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(meshPool, 'lonCell', lonCell) call mpas_pool_get_array(forcingPool, 'nTidalPotentialConstituents', nTidalConstituents) @@ -146,7 +243,7 @@ subroutine ocn_compute_tidal_potential_forcing(meshPool, forcingPool, diagnostic call mpas_pool_get_array(forcingPool, 'tidalPotentialEta', eta) call mpas_pool_get_array(diagnosticsPool, "daysSinceStartOfSim", daysSinceStartOfSim) - ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/config_tidal_potential_ramp) + ramp = tanh((2.0_RKIND*daysSinceStartOfSim)/tidalPotentialRamp) t = daysSinceStartOfSim*86400.0_RKIND do iCell = 1, nCells @@ -175,7 +272,7 @@ end subroutine ocn_compute_tidal_potential_forcing!}}} !*********************************************************************** ! -! routine ocn_tidal_potential_forcing_init +! routine ocn_vel_tidal_potential_init ! !> \brief Initializes ocean tidal protential forcing module. !> \author Steven Brus @@ -186,11 +283,22 @@ end subroutine ocn_compute_tidal_potential_forcing!}}} ! !----------------------------------------------------------------------- - subroutine ocn_tidal_potential_forcing_init(domain,err)!{{{ + subroutine ocn_vel_tidal_potential_init(domain,err)!{{{ type (domain_type), intent(inout) :: domain integer, intent(out) :: err !< Output: error flag + logical, pointer :: config_use_tidal_potential_forcing + logical, pointer :: config_use_tidal_potential_forcing_M2 + logical, pointer :: config_use_tidal_potential_forcing_S2 + logical, pointer :: config_use_tidal_potential_forcing_N2 + logical, pointer :: config_use_tidal_potential_forcing_K2 + logical, pointer :: config_use_tidal_potential_forcing_K1 + logical, pointer :: config_use_tidal_potential_forcing_O1 + logical, pointer :: config_use_tidal_potential_forcing_Q1 + logical, pointer :: config_use_tidal_potential_forcing_P1 + character (len=StrKIND), pointer :: config_tidal_potential_reference_time + type (block_type), pointer :: block_ptr type (mpas_pool_type), pointer :: forcingPool type (mpas_pool_type), pointer :: meshPool @@ -214,6 +322,17 @@ subroutine ocn_tidal_potential_forcing_init(domain,err)!{{{ err = 0 + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing', config_use_tidal_potential_forcing) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_M2', config_use_tidal_potential_forcing_M2) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_S2', config_use_tidal_potential_forcing_S2) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_N2', config_use_tidal_potential_forcing_N2) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_K2', config_use_tidal_potential_forcing_K2) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_K1', config_use_tidal_potential_forcing_K1) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_O1', config_use_tidal_potential_forcing_O1) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_Q1', config_use_tidal_potential_forcing_Q1) + call mpas_pool_get_config(ocnConfigs, 'config_use_tidal_potential_forcing_P1', config_use_tidal_potential_forcing_P1) + call mpas_pool_get_config(ocnConfigs, 'config_tidal_potential_reference_time', config_tidal_potential_reference_time) + block_ptr => domain % blocklist do while(associated(block_ptr)) call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) @@ -330,7 +449,7 @@ subroutine ocn_tidal_potential_forcing_init(domain,err)!{{{ end do - end subroutine ocn_tidal_potential_forcing_init!}}} + end subroutine ocn_vel_tidal_potential_init!}}} !*********************************************************************** @@ -636,7 +755,7 @@ function adjust_angle(arg) result(angle) !{{{ !*********************************************************************** -end module ocn_tidal_potential_forcing!}}} +end module ocn_vel_tidal_potential!}}} !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! vim: foldmethod=marker diff --git a/src/core_ocean/shared/mpas_ocn_vmix.F b/src/core_ocean/shared/mpas_ocn_vmix.F index 10c33f6c41..c1f25a2af4 100644 --- a/src/core_ocean/shared/mpas_ocn_vmix.F +++ b/src/core_ocean/shared/mpas_ocn_vmix.F @@ -26,6 +26,7 @@ module ocn_vmix use mpas_pool_routines use mpas_timer + use mpas_constants use ocn_constants use ocn_config use ocn_vmix_cvmix @@ -485,6 +486,388 @@ subroutine ocn_vel_vmix_tend_implicit(meshPool, dt, kineticEnergyCell, vertViscT end subroutine ocn_vel_vmix_tend_implicit!}}} +!*********************************************************************** +! +! routine ocn_vel_vmix_tend_implicit_variable +! +!> \brief Computes tendencies for implicit momentum vertical mixing +!> with spatially-variable bottom drag +!> \author Mark Petersen, Phillip J. Wolfram +!> \date September 2011, September 2019 +!> \details +!> This routine computes the tendencies for implicit vertical mixing for momentum +!> using computed coefficients from spatially-variable bottom drag. +!> Except for bottom drag coefficient, routine should be identifcal to +!> ocn_vel_vmix_tend_implicit above. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, dt, kineticEnergyCell, & !{{{ + vertViscTopOfEdge, layerThickness, & + layerThicknessEdge, normalVelocity, err) + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + kineticEnergyCell !< Input: kinetic energy at cell + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertViscTopOfEdge !< Input: vertical mixing coefficients + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: thickness at cell center + + real (kind=RKIND), dimension(:), intent(in) :: & + bottomDrag !< Input: bottomDrag at cell centeres + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + normalVelocity !< Input: velocity + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + layerThicknessEdge !< Input: thickness at edge + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k, cell1, cell2, N, nEdges + integer, pointer :: nVertLevels + real (kind=RKIND) :: implicitCd + integer, dimension(:), pointer :: nEdgesArray + + integer, dimension(:), pointer :: maxLevelEdgeTop + + integer, dimension(:,:), pointer :: cellsOnEdge + + real (kind=RKIND), dimension(:), allocatable :: A, B, C, velTemp + + err = 0 + + if(.not.velVmixOn) return + + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + nEdges = nEdgesArray( 1 ) + + allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) + A(1)=0.0_RKIND + + !$omp do schedule(runtime) + do iEdge = 1, nEdges + N = maxLevelEdgeTop(iEdge) + if (N .gt. 0) then + + ! Compute A(k), B(k), C(k) + ! layerThicknessEdge is computed in compute_solve_diag, and is not available yet, + ! so recompute layerThicknessEdge here. + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, N + layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) + end do + + ! average cell-based implicit bottom drag to edges + implicitCd = 0.5_RKIND*(bottomDrag(cell1) + bottomDrag(cell2)) + + ! A is lower diagonal term + do k = 2, N + A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & + / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! C is upper diagonal term + do k = 1, N-1 + C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & + / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! B is diagonal term + B(1) = 1.0_RKIND - C(1) + do k = 2, N-1 + B(k) = 1.0_RKIND - A(k) - C(k) + enddo + + ! Apply bottom drag boundary condition on the viscous term + ! second line uses sqrt(2.0*kineticEnergyEdge(k,iEdge)) + ! use implicitCd from spatially variable bottom drag + B(N) = 1.0_RKIND - A(N) + dt*implicitCd & + * sqrt(kineticEnergyCell(N,cell1) + kineticEnergyCell(N,cell2)) / layerThicknessEdge(N,iEdge) + + call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + + normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND + + end if + end do + !$omp end do + + deallocate(A,B,C,velTemp) + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_vmix_tend_implicit_spatially_variable!}}} + + +!*********************************************************************** +! +! routine ocn_vel_vmix_tend_implicit_variable_mannings +! +!> \brief Computes tendencies for implicit momentum vertical mixing +!> with spatially-variable bottom drag using Mannings friction +!> \author Mark Petersen, Phillip J. Wolfram +!> \date September 2011, September 2019 +!> \details +!> This routine computes the tendencies for implicit vertical mixing for momentum +!> using computed coefficients from spatially-variable bottom drag. +!> Except for bottom drag coefficient, routine should be identifcal to +!> ocn_vel_vmix_tend_implicit above. Cd uses Mannings' n values for the +!> Cd=g*n^2*h^{-1/3}. +! +!----------------------------------------------------------------------- + + subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forcingPool, bottomDrag, dt, & !{{{ + kineticEnergyCell, vertViscTopOfEdge, layerThickness, & + layerThicknessEdge, normalVelocity, ssh, bottomDepth, err) + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + type (mpas_pool_type), intent(inout) :: & + forcingPool !< Input: forcing information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + kineticEnergyCell !< Input: kinetic energy at cell + + real (kind=RKIND), dimension(:,:), intent(in) :: & + vertViscTopOfEdge !< Input: vertical mixing coefficients + + real (kind=RKIND), intent(in) :: & + dt !< Input: time step + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: thickness at cell center + + real (kind=RKIND), dimension(:), intent(in) :: & + bottomDrag !< Input: bottomDrag at cell centeres + + real (kind=RKIND), dimension(:), pointer :: ssh + + real (kind=RKIND), dimension(:), pointer :: bottomDepth + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + normalVelocity !< Input: velocity + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + layerThicknessEdge !< Input: thickness at edge + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iEdge, k, cell1, cell2, N, nEdges + integer, pointer :: nVertLevels + real (kind=RKIND) :: implicitCd + integer, dimension(:), pointer :: nEdgesArray + + integer, dimension(:), pointer :: maxLevelEdgeTop + + integer, dimension(:,:), pointer :: cellsOnEdge + + real (kind=RKIND), dimension(:), allocatable :: A, B, C, velTemp + + ! vegetation_drag + logical, pointer :: config_use_vegetation_drag + logical, pointer :: config_use_vegetation_manning_equation + real (kind=RKIND), pointer :: config_vegetation_drag_coefficient + real (kind=RKIND), dimension(:), pointer :: vegetationHeight + real (kind=RKIND), dimension(:), pointer :: vegetationDiameter + real (kind=RKIND), dimension(:), pointer ::vegetationDensity + real (kind=RKIND), dimension(:), pointer ::vegetationManning + integer, dimension(:), pointer ::vegetationMask + real (kind=RKIND) :: old_bottom_Cd, lambda, beta, alpha, total_h + real (kind=RKIND) :: inundation_depth, von_karman, cff1, cff2, cff3, cff4 + integer :: iCell, nCells + integer, pointer :: nCellsSolve + + err = 0 + von_karman = 0.4_RKIND + + if(.not.velVmixOn) return + + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve) + + call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + + call mpas_pool_get_config(ocnConfigs, 'config_use_vegetation_drag', config_use_vegetation_drag) + call mpas_pool_get_config(ocnConfigs, 'config_use_vegetation_manning_equation', config_use_vegetation_manning_equation) + call mpas_pool_get_config(ocnConfigs, 'config_vegetation_drag_coefficient', config_vegetation_drag_coefficient) + + call mpas_pool_get_array(forcingPool, 'vegetationMask', vegetationMask) + call mpas_pool_get_array(forcingPool, 'vegetationHeight', vegetationHeight) + call mpas_pool_get_array(forcingPool, 'vegetationDensity', vegetationDensity) + call mpas_pool_get_array(forcingPool, 'vegetationDiameter', vegetationDiameter) + call mpas_pool_get_array(forcingPool, 'vegetationManning', vegetationManning) + + nEdges = nEdgesArray( 1 ) + nCells = nCellsSolve + + allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),velTemp(nVertLevels)) + A(1)=0.0_RKIND + + ! Compute bottomDrag (Manning roughness) induced by vegetation + if (config_use_vegetation_drag .AND. config_use_vegetation_manning_equation) then + do iCell = 1, nCells + if (vegetationDensity(iCell) * vegetationHeight(iCell) * vegetationDiameter(iCell) .eq. 0.0_RKIND) then + vegetationMask(iCell) = 0 + endif + if (vegetationMask(iCell) .eq. 1) then + total_h = bottomDepth(iCell) + ssh(iCell) + old_bottom_Cd = gravity * bottomDrag(iCell)**2.0_RKIND * total_h**(1.0_RKIND/3.0_RKIND) + inundation_depth = MIN(vegetationHeight(iCell), total_h) + inundation_depth = MAX(inundation_depth, 1e-6) + lambda = vegetationDiameter(iCell) * vegetationDensity(iCell) + alpha = (config_vegetation_drag_coefficient*lambda/ & + (4.0_RKIND*von_karman**2.0_RKIND*inundation_depth**2.0_RKIND))**(1.0_RKIND/3.0_RKIND) + beta = 0.5_RKIND*alpha*inundation_depth*(1.0_RKIND - EXP(-2.0_RKIND*alpha*inundation_depth)) & + / (1.0_RKIND - EXP(-alpha*inundation_depth))**2.0_RKIND + cff1 = total_h**(2.0_RKIND/3.0_RKIND) & + * SQRT((0.5_RKIND*beta*lambda*config_vegetation_drag_coefficient*inundation_depth & + + old_bottom_Cd)/(gravity*total_h)) + cff2 = (alpha*inundation_depth)**2.0_RKIND/(1.0_RKIND - EXP(-alpha*inundation_depth)) + cff3 = (1.0_RKIND - EXP(-alpha*inundation_depth))/(alpha**2.0_RKIND*inundation_depth*total_h) + cff4 = LOG(total_h/inundation_depth) - (1.0_RKIND-inundation_depth/total_h) & + * (1.0_RKIND - 1.0_RKIND/(alpha*inundation_depth)) + vegetationManning(iCell) = cff1/(cff2*(cff3+cff4)) + vegetationManning(iCell) = MAX(bottomDrag(iCell), vegetationManning(iCell)) + else + vegetationManning(iCell) = bottomDrag(iCell) + endif + enddo + endif + + !$omp do schedule(runtime) + do iEdge = 1, nEdges + N = maxLevelEdgeTop(iEdge) + if (N .gt. 0) then + + ! Compute A(k), B(k), C(k) + ! layerThicknessEdge is computed in compute_solve_diag, and is not available yet, + ! so recompute layerThicknessEdge here. + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, N + layerThicknessEdge(k,iEdge) = 0.5_RKIND * (layerThickness(k,cell1) + layerThickness(k,cell2)) + end do + + ! average cell-based implicit bottom drag to edges and convert Mannings n to Cd + if (config_use_vegetation_drag .AND. config_use_vegetation_manning_equation) then + implicitCd = gravity*(0.5_RKIND*(vegetationManning(cell1) + vegetationManning(cell2)))**2.0 * & + (0.5_RKIND * (ssh(cell1) + ssh(cell2) + bottomDepth(cell1) + bottomDepth(cell2)))**(-1.0_RKIND/3.0_RKIND) + else + implicitCd = gravity*(0.5_RKIND*(bottomDrag(cell1) + bottomDrag(cell2)))**2.0 * & + (0.5_RKIND * (ssh(cell1) + ssh(cell2) + bottomDepth(cell1) + bottomDepth(cell2)))**(-1.0_RKIND/3.0_RKIND) + endif + + ! A is lower diagonal term + do k = 2, N + A(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k,iEdge) & + / (layerThicknessEdge(k-1,iEdge) + layerThicknessEdge(k,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! C is upper diagonal term + do k = 1, N-1 + C(k) = -2.0_RKIND*dt*vertViscTopOfEdge(k+1,iEdge) & + / (layerThicknessEdge(k,iEdge) + layerThicknessEdge(k+1,iEdge)) & + / layerThicknessEdge(k,iEdge) + enddo + + ! B is diagonal term + B(1) = 1.0_RKIND - C(1) + do k = 2, N-1 + B(k) = 1.0_RKIND - A(k) - C(k) + enddo + + ! Apply bottom drag boundary condition on the viscous term + ! second line uses sqrt(2.0*kineticEnergyEdge(k,iEdge)) + ! use implicitCd from spatially variable bottom drag + B(N) = 1.0_RKIND - A(N) + dt*implicitCd & + * sqrt(kineticEnergyCell(N,cell1) + kineticEnergyCell(N,cell2)) / layerThicknessEdge(N,iEdge) + + call tridiagonal_solve(A(2:N),B,C(1:N-1),normalVelocity(:,iEdge),velTemp,N) + + normalVelocity(1:N,iEdge) = velTemp(1:N) + normalVelocity(N+1:nVertLevels,iEdge) = 0.0_RKIND + + end if + end do + !$omp end do + + deallocate(A,B,C,velTemp) + + !-------------------------------------------------------------------- + + end subroutine ocn_vel_vmix_tend_implicit_spatially_variable_mannings!}}} + + !*********************************************************************** ! ! routine ocn_tracer_vmix_tend_implicit @@ -650,6 +1033,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo integer :: iCell, timeLevel, k, cell1, cell2, iEdge, nCells, nEdges integer, dimension(:), pointer :: nCellsArray, nEdgesArray + real (kind=RKIND), dimension(:), pointer :: bottomDrag, ssh, bottomDepth ! needed for depth-variable computation real (kind=RKIND), dimension(:,:), pointer :: normalVelocity, layerThickness, layerThicknessEdge, vertViscTopOfEdge, & vertDiffTopOfCell, kineticEnergyCell real (kind=RKIND), dimension(:,:), pointer :: vertViscTopOfCell, nonLocalSurfaceTracerFlux, tracerGroupSurfaceFlux @@ -678,6 +1062,7 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo call mpas_pool_get_subpool(forcingPool, 'tracersSurfaceFlux', tracersSurfaceFluxPool) call mpas_pool_get_array(statePool, 'normalVelocity', normalVelocity, timeLevel) call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, timeLevel) + call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel) call mpas_pool_get_array(diagnosticsPool, 'kineticEnergyCell', kineticEnergyCell) call mpas_pool_get_array(diagnosticsPool, 'layerThicknessEdge', layerThicknessEdge) @@ -690,11 +1075,14 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', maxLevelEdgeTop) call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge) + call mpas_pool_get_array(meshPool, 'bottomDepth', bottomDepth) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) call mpas_pool_get_dimension(meshPool, 'nEdgesArray', nEdgesArray) + call mpas_pool_get_array(forcingPool, 'bottomDrag', bottomDrag) + call mpas_timer_start('vmix coefs', .false.) call ocn_vmix_coefs(meshPool, statePool, forcingPool, diagnosticsPool, scratchPool, err, timeLevel) call mpas_timer_stop('vmix coefs') @@ -727,12 +1115,25 @@ subroutine ocn_vmix_implicit(dt, meshPool, diagnosticsPool, statePool, forcingPo call mpas_timer_start('vmix solve momentum', .false.) if (.not. (config_Rayleigh_friction .or. & config_Rayleigh_bottom_friction .or. & - config_Rayleigh_damping_depth_variable)) then + config_Rayleigh_damping_depth_variable .or. & + config_use_implicit_bottom_drag_variable .or. & + config_use_implicit_bottom_drag_variable_mannings)) then call ocn_vel_vmix_tend_implicit(meshPool, dt, kineticEnergyCell, & vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) else - call ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, & - vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + if (config_use_implicit_bottom_drag_variable) then + call ocn_vel_vmix_tend_implicit_spatially_variable(meshPool, bottomDrag, dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + else if (config_use_implicit_bottom_drag_variable_mannings) then + ! update bottomDrag via Cd=g*n^2*h^(-1/3) + call ocn_vel_vmix_tend_implicit_spatially_variable_mannings(meshPool, forcingPool, bottomDrag, & + dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, & + ssh, bottomDepth, err) + else if (config_Rayleigh_damping_depth_variable) then + call ocn_vel_vmix_tend_implicit_rayleigh(meshPool, dt, kineticEnergyCell, & + vertViscTopOfEdge, layerThickness, layerThicknessEdge, normalVelocity, err) + end if end if call mpas_timer_stop('vmix solve momentum') diff --git a/testing_and_setup/compass/ocean/drying_slope/analysis/comparison.py b/testing_and_setup/compass/ocean/drying_slope/analysis/comparison.py index 8f3805e2ae..8f09d8301d 100755 --- a/testing_and_setup/compass/ocean/drying_slope/analysis/comparison.py +++ b/testing_and_setup/compass/ocean/drying_slope/analysis/comparison.py @@ -79,9 +79,9 @@ def plot_MPASO(times, fileno, *args, **kwargs): for atime in times: # factor of 1e-16 needed to account for annoying round-off issue to get right time slices plottime = int((float(atime)/0.2 + 1e-16)*24.0) - ds = xr.open_mfdataset('output'+ fileno + '.nc') + ds = xr.open_dataset('output'+ fileno + '.nc') print('{} {} {}'.format(atime, plottime, ds.isel(Time=plottime).xtime.values)) - ds = ds.drop(np.setdiff1d([i for i in ds.variables], ['yCell','ssh'])) + ds = ds.drop_vars(np.setdiff1d([i for i in ds.variables], ['yCell','ssh'])) ymean = ds.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) x = ymean.yCell.values/1000.0 y = ymean.ssh.values @@ -93,7 +93,7 @@ def plot_MPASO(times, fileno, *args, **kwargs): def plot_tidal_forcing_comparison(): # data from MPAS-O on boundary for fileno in ['1','2']: - ds = xr.open_mfdataset('output' + fileno +'.nc') + ds = xr.open_dataset('output' + fileno +'.nc') ympas = ds.ssh.where(ds.tidalInputMask).mean('nCells').values x = np.linspace(0, 1.0, len(ds.xtime))*12.0 plt.plot(x, ympas, marker='o', label='MPAS-O forward' + fileno) diff --git a/testing_and_setup/compass/ocean/drying_slope/analysis/comparison_above_land.py b/testing_and_setup/compass/ocean/drying_slope/analysis/comparison_above_land.py new file mode 100755 index 0000000000..ec64971c02 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/analysis/comparison_above_land.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python +""" + +Drying slope comparison betewen MPAS-O, analytical, and ROMS results from + +Warner, J. C., Defne, Z., Haas, K., & Arango, H. G. (2013). A wetting and +drying scheme for ROMS. Computers & geosciences, 58, 54-61. + +Phillip J. Wolfram and Zhendong Cao +04/30/2019 + +""" + +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +import pandas as pd + +# render statically by default +plt.switch_backend('agg') + +def setup_fig(): + fig, ax = plt.subplots(nrows=1,ncols=1, sharex=True, sharey=True) + fig.text(0.04, 0.5, 'Channel depth (m)', va='center', rotation='vertical') + fig.text(0.5, 0.02, 'Along channel distance (km)', ha='center') + +def setup_subplot(): + plt.xlim(0,25) + plt.ylim(-11, 11) + ax = plt.gca() + ax.invert_yaxis() + ax.spines['top'].set_visible(False) + ax.spines['right'].set_visible(False) + + x = np.linspace(0,25,100) + y = -10 + 20.0/25.0*x + plt.plot(x, y, 'k-', lw=3) + +def plot_data(rval='0.0025', dtime='0.05', datatype='analytical', *args, **kwargs): + datafile = 'data/r' + rval + 'd' + dtime + '-' + datatype + '.csv' + data = pd.read_csv(datafile, header=None) + measured=plt.scatter(data[0], data[1], *args, **kwargs) + + +def plot_datasets(rval, times, fileno, plotdata=True): + for ii, dtime in enumerate(times): + if plotdata: + plot_data(rval=rval, dtime = dtime, datatype = 'analytical', + marker = '.', color = 'b', label='analytical') + plot_data(rval=rval, dtime = dtime, datatype = 'roms', + marker = '.', color = 'g', label='ROMS') + plot_MPASO([dtime], fileno, 'k-', lw=0.5, label='MPAS-O') + + if plotdata: + if ii == 0: + plt.legend(frameon=False, loc='lower left') + place_time_labels(times) + plt.text(0.5, 5, 'r = ' + str(rval)) + + +def place_time_labels(times): + locs = [9.3, 7.2, 4.2, 2.2, 1.2, 0.2] + for atime, ay in zip(times, locs): + plt.text(25.2, -10+20*ay, atime + ' days', size=8) + +def plot_MPASO(times, fileno, *args, **kwargs): + for atime in times: + # factor of 1e-16 needed to account for annoying round-off issue to get right time slices + plottime = int((float(atime)/0.2 + 1e-16)*24.0) + ds = xr.open_dataset('output'+ fileno + '.nc') + print('{} {} {}'.format(atime, plottime, ds.isel(Time=plottime).xtime.values)) + ds = ds.drop_vars(np.setdiff1d([i for i in ds.variables], ['yCell','ssh'])) + ymean = ds.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) + x = ymean.yCell.values/1000.0 + y = ymean.ssh.values + #print('ymin={} ymax={}\n{}\n{}'.format(y.min(), y.max(),x, y)) + plt.plot(x, -y, *args, **kwargs) + ds.close() + + +def plot_tidal_forcing_comparison(): + # data from MPAS-O on boundary + for fileno in ['']: + ds = xr.open_dataset('output' + fileno +'.nc') + ympas = ds.ssh.where(ds.tidalInputMask).mean('nCells').values + x = np.linspace(0, 1.0, len(ds.xtime))*12.0 + plt.plot(x, ympas, marker='o', label='MPAS-O forward' + fileno) + + # analytical case + x = np.linspace(0,12.0,100) + y = 20.0*np.sin(x*np.pi/12.0) - 10.0 + plt.plot(x, y, lw=3, color='black', label='analytical') + + plt.legend(frameon=False) + plt.ylabel('Tidal amplitude (m)') + plt.xlabel('Time (hrs)') + plt.suptitle('Tidal amplitude forcing (right side) for MPAS-O and analytical') + plt.savefig('tidalcomparison.png') + + +def main(): + ################ plot tidal forcing comparison ############### + plot_tidal_forcing_comparison() + ############################################################## + + ################ plot drying front comparison ############### + setup_fig() + times = ['0.50', '0.05', '0.40', '0.15', '0.30', '0.25'] + + # subplot r = 0.01 + setup_subplot() + plot_datasets(rval='0.01', times=times, fileno='', plotdata=False) + + plt.suptitle('Drying slope for land topography') + for outtype in ['.png','.pdf']: + plt.savefig('dryingslopecomparison' + outtype) + ############################################################## + + +if __name__ == "__main__": + main() diff --git a/testing_and_setup/compass/ocean/drying_slope/forward_idealized_transect.xml b/testing_and_setup/compass/ocean/drying_slope/forward_idealized_transect.xml new file mode 100644 index 0000000000..d1ea3a5578 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/forward_idealized_transect.xml @@ -0,0 +1,16 @@ + diff --git a/testing_and_setup/compass/ocean/drying_slope/forward_template.xml b/testing_and_setup/compass/ocean/drying_slope/forward_template.xml index d47fb1f1e5..99427d5c1d 100644 --- a/testing_and_setup/compass/ocean/drying_slope/forward_template.xml +++ b/testing_and_setup/compass/ocean/drying_slope/forward_template.xml @@ -23,6 +23,7 @@ + @@ -44,6 +45,7 @@ forcing.nc + diff --git a/testing_and_setup/compass/ocean/drying_slope/forward_vegetation_properties.xml b/testing_and_setup/compass/ocean/drying_slope/forward_vegetation_properties.xml new file mode 100644 index 0000000000..db1baf2b12 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/forward_vegetation_properties.xml @@ -0,0 +1,24 @@ + diff --git a/testing_and_setup/compass/ocean/drying_slope/init_idealized_transect.xml b/testing_and_setup/compass/ocean/drying_slope/init_idealized_transect.xml new file mode 100644 index 0000000000..adf984b00c --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/init_idealized_transect.xml @@ -0,0 +1,20 @@ + diff --git a/testing_and_setup/compass/ocean/drying_slope/init_template.xml b/testing_and_setup/compass/ocean/drying_slope/init_template.xml index ed4c8e56f2..41e3b30ada 100644 --- a/testing_and_setup/compass/ocean/drying_slope/init_template.xml +++ b/testing_and_setup/compass/ocean/drying_slope/init_template.xml @@ -11,6 +11,7 @@ + @@ -24,6 +25,7 @@ 0000_00:00:01 + diff --git a/testing_and_setup/compass/ocean/drying_slope/init_vegetation_properties.xml b/testing_and_setup/compass/ocean/drying_slope/init_vegetation_properties.xml new file mode 100644 index 0000000000..f072b0aefd --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/init_vegetation_properties.xml @@ -0,0 +1,16 @@ + diff --git a/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/comparison.py b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/comparison.py new file mode 100755 index 0000000000..8a6ec6ef54 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/comparison.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python +# coding: utf-8 +# Authors: Zhendong Cao, Phillip Wolfram +# Date: May 2020 + +# import libs +import numpy as np +import xarray as xr +import matplotlib.pyplot as plt +import os +#import cv2 +# plot tidal amplitude at open boundary + +fname1 = 'NoVegetation' +fname2 = 'ConstantVegManning' +fname3 = 'VegManningEquation' + +# mpas-o +ds1 = xr.open_dataset(fname1+'.nc') +ds2 = xr.open_dataset(fname2+'.nc') +ds3 = xr.open_dataset(fname3+'.nc') + +x = np.linspace(0, 1.0, len(ds1.xtime))*48.0 +ympas1 = ds1.ssh.where(ds1.tidalInputMask).mean('nCells').values +ympas2 = ds2.ssh.where(ds2.tidalInputMask).mean('nCells').values +ympas3 = ds3.ssh.where(ds3.tidalInputMask).mean('nCells').values + +plt.plot(x, ympas1, marker='o', label=fname1) +plt.plot(x, ympas2, marker='^', label=fname2) +plt.plot(x, ympas3, marker='*', label=fname3) + +# analytical +x = np.linspace(0,48.0,100) +y = 1.5*np.sin(x*np.pi/12.0) - 3.0 +plt.plot(x, y, 'k-', lw=3, label='analytical') + +plt.legend(frameon=False, loc='upper right', fontsize=8) +plt.ylabel('Tidal amplitude (m)') +plt.xlabel('Time (hrs)') +plt.suptitle('Tidal amplitude forcing for test cases and analytical') +plt.savefig('tidalcomparison.png',dpi=300) +plt.close() + +# extract bottomDepth +bottomDepth = ds1.bottomDepth.values.reshape((116, 4))[:,1] +x_along = np.linspace(0,8,116) +# plot water elevation at selected time slices +times = np.linspace(0.0,2.0,41) +ii = 0 +figname_prefix = 'Waterlevel' +for atime in times: + plt.figure(figsize=(12,6)) + plt.xlim(0, 8) + plt.ylim(-4, 0) + plottime = int((float(atime)/0.2)*24.0) + timestr = ds1.isel(Time=plottime).xtime.values + ymean = ds1.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) + x = ymean.yCell.values/1000.0 + y = ymean.ssh.values + plt.plot(x, y,'-k',label=fname1) + + ymean = ds2.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) + x = ymean.yCell.values/1000.0 + y = ymean.ssh.values + plt.plot(x, y,'-b',label=fname2) + + ymean = ds3.isel(Time=plottime).groupby('yCell').mean(dim=xr.ALL_DIMS) + x = ymean.yCell.values/1000.0 + y = ymean.ssh.values + plt.plot(x, y,'-g',label=fname3) + + plt.fill_between(x=x_along,y1=-10*np.ones(x_along.shape), y2=-bottomDepth, color='grey',zorder=9) + plt.xlabel('Cross shore distance (km)') + plt.ylabel('Bathymetry (m)') + atime='%.2f'%atime + plt.title('time='+atime+'days') + plt.legend(frameon=False, loc='upper right') + figname = '{}{:02d}.png'.format(figname_prefix,ii) + plt.savefig(figname, dpi=100) + plt.close() + print('Time '+atime+ ' done!') + ii += 1 +# make a movie +os.system("ffmpeg -r 1 -i Waterlevel%02d.png -vcodec mpeg4 -y movie.mp4") diff --git a/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_analysis.xml b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_analysis.xml new file mode 100644 index 0000000000..0e7719e952 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_analysis.xml @@ -0,0 +1,28 @@ + + + + + + + + + + + + output1.nc + vtk_output1 + allOnCells + maxEdges=0 + nVertLevels=0:10 + + + + output2.nc + vtk_output2 + allOnCells + maxEdges=0 + nVertLevels=0:10 + + + + diff --git a/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_driver.xml b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_driver.xml new file mode 100644 index 0000000000..7fa1481248 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_driver.xml @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_forward1.xml b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_forward1.xml new file mode 100644 index 0000000000..49441e4608 --- /dev/null +++ b/testing_and_setup/compass/ocean/drying_slope/marsh_flooding/idealized_transect/config_forward1.xml @@ -0,0 +1,22 @@ + + + + + + + +