Updating docs

Docs/sphinx_doc/index.rst

@@ -49,6 +49,7 @@ In addition to this documentation, there is API documentation for ERF generated
 
    theory/DryEquations.rst
    theory/WetEquations.rst
+   theory/Initialization.rst
    theory/Buoyancy.rst
    theory/Microphysics.rst
    theory/DNSvsLES.rst

Docs/sphinx_doc/theory/Buoyancy.rst

@@ -109,3 +109,60 @@ which removes the need to compute horizontal averages of these quantities.   Thi
 
 We note that this version of the buoyancy force matches that given in Marat F. Khairoutdinov and David A. Randall's paper (J. Atm Sciences, 607, 1983)
 if we neglect :math:`\frac{p^\prime}{\bar{p_0}}`.
+
+Type 4
+------
+This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002benchmark`_.
+
+.. _`khairoutdinov2003cloud`: https://journals.ametsoc.org/view/journals/atsc/60/4/1520-0469_2003_060_0607_crmota_2.0.co_2.xml
+.. _`bryan2002benchmark`: https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml
+
+.. math::
+
+    \begin{equation}
+    \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg),
+    \end{equation}
+
+The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as
+
+.. math::
+
+    \rho = \frac{p (1 + q_v + q_c + q_p)}{R_dT\Bigg(1 + \cfrac{R_v}{R_d}q_v\Bigg)}
+
+This can be written using binomial expansion as
+
+.. math::
+
+    \begin{align*}
+    \rho &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 + \frac{R_v}{R_d}q_v\Bigg)^{-1} \\
+    &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 - \frac{R_v}{R_d}q_v + O(q_v^2)\Bigg) \\
+    &= \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v +  \text{H.O.T. such as } O(q_v^2) + O(q_vq_c)\Bigg) \\
+    &\approx \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v\Bigg)
+    \end{align*}
+
+Taking log on both sides, we get
+
+.. math::
+
+    \log{\rho} = \log{p} - \log{R_d} - \log{T} + \log(1 - 0.61 q_v + q_c + q_p)
+
+Taking derivative gives
+
+.. math::
+
+    \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + \frac{(-0.61 q_v' + q_c' + q_p')}{(1 - 0.61 q_v + q_c + q_p)}
+
+Using :math:`- 0.61 q_v + q_c + q_p \ll 1`, we have
+
+.. math::
+
+    \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + (-0.61 q_v' + q_c' + q_p')
+
+Since the background values of cloud water and precipitate mass mixing ratios -- :math:`q_c` and :math:`q_p` are zero, we have :math:`q_c' = q_c` and :math:`q_p' = q_p`. Hence, we have
+
+.. math::
+
+    \begin{equation}
+    \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg),
+    \end{equation}
+

Docs/sphinx_doc/theory/Initialization.rst

@@ -0,0 +1,149 @@
+
+ .. role:: cpp(code)
+    :language: c++
+
+ .. role:: f(code)
+    :language: fortran
+
+.. _Initialization:
+
+Initialization
+==================
+
+To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient
+and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the
+cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density
+:math:`\rho = \rho_d(1 + q_t)`, where :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water, instead
+of the dry density :math:`\rho_d` for cases without moisture.
+
+Computation of the dry density
+-------------------------------
+We express the total pressure :math:`p` as:
+
+.. math::
+   p = \rho_d R_d T_v.
+
+By definition, we have:
+
+.. math::
+   T_v = \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p},
+
+.. math::
+   T = \theta_d\left(\frac{p}{p_0}\right)^{R_d/C_p}.
+
+This gives:
+
+.. math::
+   p = \rho_d R_d \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p},
+
+which, on rearranging, yields:
+
+.. math::
+   p = p_0\left(\frac{\rho_d R_d \theta_m}{p_0}\right)^\gamma.
+
+To obtain :math:`\theta_m`, consider the density of dry air:
+
+.. math::
+   \rho_d = \frac{p - p_v}{R_d T}.
+
+Substituting for :math:`\rho_d` from the above equation, we get:
+
+.. math::
+   \frac{p}{T_v} = \frac{p - p_v}{T},
+
+which implies:
+
+.. math::
+   \frac{T_v}{T} = \frac{p}{p - p_v} = \frac{1}{1-\left(\cfrac{p_v}{p}\right)}.
+
+We also have:
+
+.. math::
+   q_v = \frac{\rho_v}{\rho_d} = \frac{p_v}{R_v T}\frac{R_d T}{p-p_v} = \frac{r p_v}{p - p_v},
+
+where :math:`p_v` is the partial pressure of water vapor, :math:`r = R_d/R_v \approx 0.622`, and :math:`q_v` is the vapor mass mixing ratio—the ratio of the
+mass of vapor to the mass of dry air. Rearranging and using :math:`q_v \ll r`, we get:
+
+.. math::
+   \frac{p_v}{p} = \frac{1}{1 + \left(\cfrac{r}{q_v}\right)} \approx \frac{q_v}{r},
+
+which, on substitution in the equation for :math:`\frac{T_v}{T}`, gives:
+
+.. math::
+   \frac{T_v}{T} = \frac{1}{1 - \left(\cfrac{q_v}{r}\right)}.
+
+As :math:`q_v \ll 1`, a binomial expansion, ignoring higher-order terms, gives:
+
+.. math::
+   T_v = T\left(1 + \frac{R_v}{R_d}q_v\right).
+
+Hence, the density of dry air is given by:
+
+.. math::
+   \rho_d = \frac{p}{R_d T_v} = \frac{p}{R_d T\left(1 + \cfrac{R_v}{R_d}q_v\right)}.
+
+
+Initialization with a second-order integration of the hydrostatic equation
+----------------------------------------------------------------------------
+
+We have the hydrostatic equation given by
+
+.. math::
+
+    \frac{\partial p}{\partial z} = -\rho g,
+
+where :math:`\rho = \rho_d(1 + q_t)`, :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water. Using an average value of :math:`\rho` for the integration from :math:`z = z(k-1)` to :math:`z(k)`, we get
+
+.. math::
+
+    p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z.
+
+The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write :math:`\rho(k) = f(p(k))`. Hence
+
+.. math::
+
+    p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z.
+
+Now, we define
+
+.. math::
+
+    F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0.
+
+This is a non-linear equation in :math:`p(k)`. Consider a Newton-Raphson iteration (where :math:`n` denotes the iteration number) procedure
+
+.. math::
+
+    F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0,
+
+which implies
+
+.. math::
+
+    \delta p = -\frac{F}{F'},
+
+with the gradient being evaluated as
+
+.. math::
+
+    F' = \frac{F(p+\epsilon) - F(p)}{\epsilon},
+
+and the iteration update is given by
+
+.. math::
+
+    p^{n+1} = p^n + \delta p.
+
+For the first cell (:math:`k=0`), which is at a height of :math:`z = \frac{\Delta z}{2}` from the base, we have
+
+.. math::
+
+    p(0) = p_0 - \rho(0)g\frac{\Delta z}{2},
+
+where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we define
+
+.. math::
+
+    F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2},
+
+and the Newton-Raphson procedure is the same.

Docs/sphinx_doc/theory/Microphysics.rst

@@ -7,8 +7,12 @@
 
 .. _Microphysics:
 
+Microphysics model
+====================
+
 Kessler Microphysics model
-===========================
+---------------------------
+The Kessler microphysics model is a simple version of cloud microphysics which has precipitation only in the form of rain. Hence :math:`q_p = q_r`.
 Governing equations for the microphysical quantities for Kessler microphysics from `gabervsek2012dry`_ are
 
 .. math::
@@ -24,14 +28,78 @@ Governing equations for the microphysical quantities for Kessler microphysics fr
 where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor,
 :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of
 rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given
-below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`.
+below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. The parametrization
+of the source terms are given below.
 
 .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml
+
+.. raw:: latex
+
+   \newgeometry{left=2cm,right=2cm,top=2cm,bottom=2cm}
+
+Rate of condensation of water vapor/evaporation of cloud water
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+From `klemp1978simulation`_, we have the following expressions.
+
 .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml
 
+If the air is not saturated, i.e. :math:`q_v > q_{vs}`
+
+.. math::
+    C_c = \frac{q_v - q_{vs}}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}}
+
+If the air is not saturated, i.e. :math:`q_v < q_{vs}`, then cloud water evaporates to water vapor at the rate
+
+.. math::
+    E_c = \frac{q_{vs} - q_v}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}}
+
+Rate of autoconversion of cloud water into rain
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The rate of autoconversion of cloud water into rain is given by
+
+.. math::
+    A_c = k_1(q_c - a)
+
+where :math:`k_1 = 0.001` s\ :sup:`-1`, :math:`a = 0.001` kg kg\ :sup:`-1`.
+
+Rate of accretion of cloud water into rain water drops
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The rate of accretion of cloud water into rain water drops is given by
+
+.. math::
+    K_c = k_2q_cq_r^{0.875}
+
+where :math:`k_2= 2.2` s\ :sup:`-1`.
+
+The rate of evaporation of rain into water vapor
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The rate of evaporation of rain into water vapor is given by
+
+.. math::
+    E_r = \cfrac{1}{\overline{\rho}}\cfrac{(1- q_v/q_s)C(\overline{\rho}q_r)^{0.525}}{5.4\times10^5 + 2.55\times10^6/(\overline{p}q_s)},
+
+where the ventilation factor :math:`C = 1.6 + 124.9(\overline{\rho}q_r)^{0.2046}`.
+
+Terminal fall velocity of rain
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The terminal fall velocity of rain is given by
+
+.. math::
+    V = 3634(\overline{\rho}q_r)^{0.1346}\Bigg(\cfrac{\overline{\rho}}{\rho_0}\Bigg)^{-\frac{1}{2}}~\text{[cm/s]}
+
+.. raw:: latex
+
+   \restoregeometry
+
+
 
 Single Moment Microphysics Model
-===================================
+----------------------------------
 The conversion rates among the moist hydrometeors are parameterized assuming that
 
 .. math::
@@ -54,7 +122,7 @@ Marat F. Khairoutdinov and David A. Randall's (J. Atm Sciences, 607, 1983).
 The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/~marat/SAM.html)
 
 Accretion
-------------------
+~~~~~~~~~~~~
 There are several different type of accretional growth mechanisms that need to be included; these describe
 the interaction of water vapor and cloud water with rain water.
 
@@ -69,21 +137,21 @@ The accretion of raindrops by accretion of cloud water is
    Q_{racw} = \frac{\pi E_{RW}n_{0R}\alpha q_{c}\Gamma(3+b)}{4\lambda_{R}^{3+b}}(\frac{\rho_{0}}{\rho})^{1/2}
 
 The bergeron Process
-------------------------
+~~~~~~~~~~~~~~~~~~~~~~
 The cloud water transform to snow by deposition and rimming can be written as
 
 .. math::
    Q_{sfw} = N_{150}\left(\alpha_{1}m_{150}^{\alpha_{2}}+\pi E_{iw}\rho q_{c}R_{150}^{2}U_{150}\right)
 
 Autoconversion
-------------------------
+~~~~~~~~~~~~~~~~~~~~~~
 The collision and coalescence of cloud water to from raindrops is parameterized as following
 
 .. math::
    Q_{raut} = \rho\left(q_{c}-q_{c0}\right)^{2}\left[1.2 \times 10^{-4}+{1.569 \times 10^{-12}N_{1}/[D_{0}(q_{c}-q_{c0})]}\right]^{-1}
 
 Evaporation
-------------------------
+~~~~~~~~~~~~~~~~~~~~~~
 The evaporation rate of rain is
 
 .. math::

Docs/sphinx_doc/theory/WetEquations.rst

@@ -138,4 +138,3 @@ In this set of equations, the subgrid turbulent parameterization effects are inc
 :math:`\mathbf{F}` stands for the external force, and :math:`Q` and :math:`F_Q` represent the mass and energy transformation
 of water vapor to/from water through condensation/evaporation, which is determined by the microphysics parameterization processes.
 :math:`\mathbf{B}` is the buoyancy force, which is defined in :ref:`Buoyancy <Buoyancy>`.
-

Exec/SquallLine_2D/inputs_moisture

@@ -70,7 +70,7 @@ erf.molec_diff_type = "ConstantAlpha"
 erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
 erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
 erf.alpha_T = 00.0 # [m^2/s]
-erf.alpha_C = 25.0
+erf.alpha_C = 100.0
 
 erf.moisture_model = "Kessler"
 erf.use_moist_background = true