• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Iterative Methods for Solving the Nonlinear Balance Equation with Optimal Truncation

    2021-04-20 00:41:56QinXUandJieCAO
    Advances in Atmospheric Sciences 2021年5期

    Qin XU and Jie CAO

    1NOAA/National Severe Storms Laboratory, Norman, Oklahoma 73069, USA

    2Cooperative Institute for Mesoscale Meteorological Studies, University of Oklahoma, Norman, Oklahoma 73072, USA

    3Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China

    (Received 27 August 2020; revised 16 November 2020; accepted 20 December 2020)

    ABSTRACT Two types of existing iterative methods for solving the nonlinear balance equation (NBE) are revisited. In the first type, the NBE is rearranged into a linearized equation for a presumably small correction to the initial guess or the subsequent updated solution. In the second type, the NBE is rearranged into a quadratic form of the absolute vorticity with the positive root of this quadratic form used in the form of a Poisson equation to solve NBE iteratively. The two methods are rederived by expanding the solution asymptotically upon a small Rossby number, and a criterion for optimally truncating the asymptotic expansion is proposed to obtain the super-asymptotic approximation of the solution. For each rederived method, two iterative procedures are designed using the integral-form Poisson solver versus the over-relaxation scheme to solve the boundary value problem in each iteration. Upon testing with analytically formulated wavering jet flows on the synoptic, sub-synoptic and meso-α scales, the iterative procedure designed for the first method with the Poisson solver, named M1a, is found to be the most accurate and efficient. For the synoptic wavering jet flow in which the NBE is entirely elliptic, M1a is extremely accurate. For the sub-synoptic wavering jet flow in which the NBE is mostly elliptic,M1a is sufficiently accurate. For the meso-α wavering jet flow in which the NBE is partially hyperbolic so its boundary value problem becomes seriously ill-posed, M1a can effectively reduce the solution error for the cyclonically curved part of the wavering jet flow, but not for the anti-cyclonically curved part.

    Key words: nonlinear balance, iterative method, optimal truncation

    1. Introduction

    For flows of synoptic and sub-synoptic scales in the middle and high latitudes, the nonlinear balance equation(NBE) links the streamfunction field with the geopotential field more accurately than the geostrophic balance (Bolin,1955; Charney, 1955). However, solving the streamfunction from the NBE for a given geopotential field can be very challenging due to complicated issues on the existence of solution in conjunction with difficulties caused by nonlinearity(Courant and Hilbert, 1962). It is well known mathematically that the NBE is a special case of the Monge-Ampere differential equation for the streamfunction (Charney, 1955). If the geostrophic vorticity (that is, the vorticiy of geostrophic flow associated with the given geopotential field) is larger than ?f/2 for a constant f where f is the Coriolis parameter,then the NBE is of the elliptic type and its associated boundary value problem can have no more than two solutions (see Section 6.3 in Chapter 4 of Courant and Hilbert, 1962). If the geostrophic vorticity is smaller than ?f/2 in a local area,then the NBE becomes locally hyperbolic. In this case, the boundary value problem becomes ill-posed and thus may have no solution although the NBE can be integrated along the characteristic lines within the locally hyperbolic area(see Section 3 of Appendix I in Chapter 5 of Courant and Hilbert, 1962).

    To avoid the complication and difficulties caused by the local non-ellipticity in solving the NBE, one can simply enforce the ellipticity condition to a certain extent by slightly smoothing or adjusting the given geopotential field.This type of treatment has been commonly used in iterative methods to solve the NBE as a boundary value problem(Bolin, 1955, 1956; Shuman, 1955, 1957; Bushby and Huckle, 1956; Miyakoda, 1956; Arnason, 1958; Bring and Charasch, 1958; Liao and Chow, 1962; Asselin, 1967;Bijlsma and Hoogendoorn, 1983). However, regardless of the above treatment, the convergence properties of these or any other iterative methods can be not only scale-dependent but also flow-dependent and thus very difficult to study theoretically and rigorously.

    The above reviewed iterative methods can be classified into two types. In the first type (originally proposed by Bolin, 1955), the NBE is transformed into a linearized equation for a presumably small correction to the initial guess or to the subsequent updated solution when this linearized equation is solved iteratively. In the second type (originally proposed by Shuman, 1955, 1957; Miyakoda, 1956), the NBE is rearranged into a quadratic form of the absolute vorticity and the positive root of this quadratic form is used in the form of a Poisson equation to solve for the streamfunction iteratively. The initial guess for both types is the geostrophic streamfunction. Their convergence properties have been analyzed theoretically, but the analysis lacks rigor and generality, because the coefficients of the linearized differential operator for the first type and the forcing terms on the right-hand side of the iterative form of the linearized equation for the second type were functions of space but treated as constants (Arnason, 1958; Bijlsma and Hoogendoorn,1983). Therefore, the convergence properties of the previously iterative methods were examined mainly through numerical experiments. Besides, due to the very limited computer memories and speed at that time, these earlier iterative methods employed the memory-saving sequential relaxation scheme based on the classical Liebmann-type iteration algorithm (Southwell, 1946) and applied to coarse resolution grids for large-scale flows. The sequential relaxation and successive over-relaxation (SOR) schemes have been used in the second type of iterative method (Shuman, 1955,1957) to solve the NBE for hurricane flows (Zhu et al.,2002). However, using these iterative methods to solve the NBE still faces various difficulties especially when the spatial scale is smaller than the sub-synoptic scale. In particular, there are unaddressed challenging issues concerning whether and how the solutions can be obtained approximately and efficiently through a limited numbers of iterations,especially when the NBE becomes locally hyperbolic (due mainly to reduced spatial scales) and thus the iterative methods fail to converge.

    This paper aims to address the above challenging issues. In particular, we will rederive the above two types of iterative methods formally and systematically by expanding the solution asymptotically upon a small Rossby number and substituting it into the NBE. Since the asymptotic expansion is not ensured to converge, especially when the Rossby number is not sufficiently small, the concept of optimal truncation of asymptotic expansion is employed and a criterion is proposed to obtain the super-asymptotic approximation of the solution based on the heuristic theory of asymptotic analysis (Boyd, 1999). As will be seen in this paper, by employing the optimal truncation, the issue of non-convergence of the iterative methods caused by the increase of Rossby number can be addressed to a certain extent. Besides, the recently developed Poisson solver based on integral formulas (Cao and Xu, 2011; Xu et al., 2011) will be used in comparison with the aforementioned classical SOR scheme to solve the boundary value problem in each iterative step. In particular, for flows of sub-synoptic scale or meso-

    α

    scale,the NBE can become locally hyperbolic and the solution will be checked in this paper via the proposed optimal truncation under certain conditions.

    The paper is organized as follows. The next section presents formal and systematical derivations of the above reviewed two iterative methods. Section 2 formulates the criterion for optimal truncation, and section 3 constructs four different iterative procedures with optimal truncation and designs numerical experiments for testing the iterative procedures. Section 4 examines and compares the results of experiments performed with the four iterative procedures, followed by conclusions in section 5.

    2. Derivations of two iterative methods

    2.1. Scaling and asymptotic expansion based on small Rossby number

    The NBE can be expressed in the following form (Charney, 1955):

    where ?

    D

    denotes the domain boundary, and ψ≡?/

    f

    is the global geostrophic streamfunction (Kuo, 1959; Charney and Stern, 1962; Schubert et al., 2009).Formally, we can scale x and y by L, scale f = f+ f

    by f, and scale

    ψ

    and

    ?

    by UL and fUL, respectively, where U is the horizontal velocity scale, L is the horizontal length scale, fis a constant reference value of f which can be the value of f at the domain center. The scaled variables are still denoted by their respectively original symbols, so the NBE can have the following non-dimensional form:

    where F = f

    /(fRo) ≤ O(1) and O() is the “order-of-magnitude” symbol. Thus,

    ψ

    can have the following asymptotic expansion:

    where

    ψ

    =

    ψ

    and Σdenotes the summation over k from 1 to ∞. The kth order truncation of the asymptotic expansion of

    ψ

    in Eq. (4) is given by ψ≡ψ+Σεδψ′ , where Σdenotes the summation over k' from 1 to k. Formally,

    ψ

    =

    ψ

    +O(

    ε

    ), so

    ψ

    is accurate up to O(

    ε

    ) as an approximation of

    ψ

    .By substituting Eq. (4) into Eq. (3) and Eq. (1b), and then collecting terms of the same order of

    ε

    , we obtain

    Here, Eq. (5) gives a formal series of linearized equations for computing

    δψ

    consecutively from

    δψ

    to increasingly higher-order term in the expansion of

    ψ

    with the boundary conditions of

    δψ

    = 0 (on ?

    D

    for k = 1, 2, 3, …). The equations in Eq. (5), however, are inconvenient to use, because the equation at each given order becomes increasingly complex as the order k increases. It is thus desirable to modify Eq. (5) into a recursive form, and this can be done nonuniquely by first combining the equations in Eq. (5) with?( f

    ψ

    ) = ?

    ?

    into a series of equations for

    ψ

    (instead of

    δψ

    ) and then adding properly selected higher-order terms to the equation for

    ψ

    at each order without affecting the order of accuracy of the equation. In particular, two different modifications will be made in the next two subsections.From these two modifications, the two types of iterative methods reviewed in the introduction for solving the NBE can be derived formally and systematically via the asymptotic expansion of

    ψ

    in Eq. (4).

    2.2. Derivation of method-1

    The equations in Eq. (5) can be combined with ?(f

    ψ

    ) =?

    ?

    at O(

    ε

    ) into a series of equations for

    ψ

    defined in Eq.(4) as shown blow:

    Formally

    ψ

    is accurate up to O(

    ε

    ) and so is ?(f

    ψ

    ) on the left-hand side of the above kequation. This implies that the kequation is accurate only up to O(

    ε

    ), so the last term O(

    ε

    ) (that represents all the high-order terms) on the right-hand side can be neglected without degrading the order of accuracy of the equation. This leads to the following recursive form of equation and boundary condition for solving the NBE iteratively:

    If

    ε

    is sufficiently small to ensure the convergence of the asymptotic expansion in Eq. (4), then

    ψ

    ψ

    gives the solution of the NBE in the limit of k → ∞.Substituting ε?

    F

    =?

    f

    /

    f

    and

    ε

    = Ro ≡ U/fL into Eq.(7) gives the dimensional form of Eq. (7):

    For f = constant, Eq. (8a) recovers Eq. (5) of Bushby and Huckle (1956), but this recursive form of equation is derived here formally and systematically via the asymptotic expansion of the solution in Eq. (4). Substituting the dimensional form of

    ψ

    =

    ψ

    +

    ε

    δψ

    , that is,

    ψ

    =

    ψ

    +

    δψ

    into Eq. (8) gives

    where N() is the nonlinear differential operator defined in Eq. (1a). Analytically, Eq. (9a) is identical to Eq. (8a) but expressed in an incremental form. Numerically, however,solving

    δψ

    from Eq. (9) and updating

    ψ

    to

    ψ

    =

    ψ

    +

    δψ

    iteratively does not give exactly the same solution as that obtained by solving

    ψ

    from Eq. (8) iteratively. According to our additional numerical experiments (not shown),the solutions obtained from Eq. (8) are less accurate (by about an order of magnitude for the case of Ro = 0.1) than their counterpart solutions obtained from Eq. (9), so the non-incremental form of boundary value problem in Eq. (8)will not be considered in this paper.

    2.3. Derivation of method-2

    The equation for

    ψ

    in Eq. (6.4) can be multiplied by 2 and rewritten as

    where

    ψ

    =

    ψ

    +

    ε

    δψ

    =

    ψ

    + O(

    ε

    ) and ?(f

    ψ

    ) = f

    ζ

    +( ?f)·

    ?

    ψ

    + ?·(

    ψ

    ? f ) = f

    ζ

    +

    ε

    ( ?F)·

    ?

    ψ

    +

    ε

    ?·(

    ψ

    ? F) = f

    ζ

    +

    ε

    ( ?F)·(

    ?

    ψ

    ) +

    ε

    ?· (

    ψ

    ? F) + O(

    ε

    ) are used. One can verify that ?4

    ε

    J(?ψ, ?ψ) =

    ε

    (

    ζ

    ? A? B) =

    εζ

    ?

    ε

    (A+ B) + O(

    ε

    ) where

    ζ

    = ?

    ψ

    , A≡(???)

    ψ

    , B≡ 2??ψ, and A= (???)(

    ψ

    +

    ε

    δψ

    ) = A+ O(

    ε

    ) and B= 2 ??(

    ψ

    +

    ε

    δψ

    ) = B+O(

    ε

    ) are used. Substituting these into Eq. (10) gives

    This leads to the following recursive form of equation that is accurate up to O(

    ε

    ):

    Substituting ε ?

    F

    =?

    f

    /

    f

    and

    ε

    = Ro ≡ U/(fL) into Eq.(12a) gives its dimensional form which can be rewritten as

    The non-negative condition of (f +

    ζ

    )≥ 0 requires M≥ 0 on the right-hand side of Eq. (12b). Also, as a quadratic equation of f +

    ζ

    for given

    ?

    and

    ψ

    , Eq. (12b) has two roots, but only the positive root, given by f +

    ζ

    =M, is physically acceptable (because f +

    ζ

    ≥ 0 is required for stably balanced flow). This leads to the following recursive form of equation and boundary condition for solving the NBE iteratively:

    where M≥ 0 is ensured by setting M= 0 when the computed Mfrom the previous step becomes negative. Here,Eq. (13a) gives essentially the same recursive form of equation as that in Eq. (8) of Shuman (1957) for solving the NBE iteratively, but this recursive form of equation is derived here via the asymptotic expansion of the solution in Eq. (4).

    3. Iterative procedures with optimal truncation and experiment design

    3.1. Criterion for optimal truncation

    When the Rossby number is not sufficiently small to ensure the convergence of the asymptotic expansion, the optimal truncation of the asymptotic expansion of

    ψ

    in Eq.(4) can be determined (Boyd, 1999) by an empirical criterion in the following dimensional form:

    where N() is the function form defined in Eq. (1a), K is the number of optimal truncation, E[N(

    ψ

    )] ≡ ||

    ε

    [N(

    ψ

    )]||′, || ||′denotes the root-mean-square (RMS) of discretized field of the variable inside || ||′ computed over all the interior grid points (excluding the boundary points) of domain D, and

    ε

    [N(

    ψ

    )] ≡ [N(

    ψ

    ) ? N(

    ψ

    )]/||N(

    ψ

    )||′ = [N(

    ψ

    ) ? ?

    ?

    ]/||

    ?

    ?

    ||′is the relative error of N(

    ψ

    ) with respect to N(

    ψ

    ) which is also the normalized (by ||

    ?

    ?

    ||′) residual error of the NBE caused by the approximation of

    ψ

    ψ

    , and

    ψ

    denotes the true solution. Here, E[N(

    ψ

    )] is expected to be the global minimum of E[N(

    ψ

    )]. If E[N(

    ψ

    )] does not oscillate as k increases, then it is sufficient to set m = 1 in Eq. (14). Otherwise, m should be sufficiently large to ensure E[N(

    ψ

    )] be the global minimum of E[N(

    ψ

    )].

    3.2. Iterative procedures

    The iterative procedure for method-1 performs the following steps:

    1. Start from k = 0 and set

    ψ

    =

    ψ

    ?

    /f in D and ?

    D

    .2. Substitute

    ψ

    (=

    ψ

    for k = 1) into N(

    ψ

    ) to compute the right-hand-side of Eq. (9a), and then solve the boundary value problem in Eq. (9) for

    δψ

    .3. Substitute

    ψ

    =

    ψ

    +

    αδψ

    into ||N(

    ψ

    ) ? ?

    ?

    ||′ and save the computed ||N(

    ψ

    ) ? ?

    ?

    ||′ where

    α

    is an adjustable parameter in the range of 0 <

    α

    ≤ 1.4. If k ≥ 2m, then find min||N(

    ψ

    ) ? ?

    ?

    ||′, say at k′ =K′, for k′ = k, k ? 1, … k ? 2m. If K′ < k ? m, then K = K′and

    ψ

    gives the optimally truncated solution–the final solution that ends the iteration. Otherwise, go back to step 2.When the Poisson solver (or SOR scheme) is used to solve the boundary value problem in the above step 2, the iterative procedure designed for method-1 is named M1a (or M1b). For the Poisson solver used in this paper, the internally induced solution is obtained by using the scheme S2 described in section 2.1 of Cao and Xu (2011) and the externally induced solution obtained by using the Cauchy integral method described in section 4.1 of Cao and Xu (2011). For M1a with Ro < 0.4 (or Ro = 0.4), it is sufficient to set m = 1 and

    α

    = 1 (or 1/2). For M1b, it is sufficient to set m = 3 and

    α

    = 1.

    The iterative procedure for method-2 performs the following steps:

    1. Start from k = 0 and set

    ψ

    =

    ψ

    ?

    /f in D and ?

    D

    .2. Substitute

    ψ

    into Mdefined in Eq. (12b) to compute the right-hand-side of Eq. (13a), and then solve the boundary value problem in Eq. (13) for

    ψ

    .3. Compute and save ||N(

    ψ

    ) ? ?

    ?

    ||′.

    4. Perform this step as described above for step 4 of method-1.

    When the Poisson solver (or SOR scheme) is used to solve boundary value problem in the above step 2, the iterative procedure designed for method-2 is named M2a (or M2b). For M2a and M2b, it is sufficient to set m = 1 and

    α

    = 1.

    3.3. Experiment design

    To examine and compare the accuracies and computational efficiencies of the four iterative procedures, the true streamfunction field is formulated for a wavering jet flow by

    and the associated velocity components are given by

    and

    where U = 20 m sis the maximum zonal speed of the wavering jet flow, y = ?0.25Lcos(πx′/L) is the longitudinal location (in y-coordinate) of the wavering jet axis as a function of x′ = x ? x, and xis the zonal location of wave ridge. By setting the half-wavelength L to 2000, 1000 and 500 km, the flow fields formulated in Eqs. (15) and (16) resemble wavering westerly jet flows on the synoptic, sub-synoptic and meso-

    α

    scales, respectively (as often observed on northernhemisphere mid-latitude 500 hPa weather maps).Four sets of experiments are designed to test and compare the iterative procedures with

    ψ

    given in Eq. (15) over a square domain of D ≡ [?L ≤ x ≤ L, ?L ≤ y ≤ L]. The first set consists of four experiments to test the four iterative procedures (that is, M1a, M1b, M2a and M2b) on the synoptic scale by setting L = 2000 km and x= 0 for

    ψ

    in Eq. (15).The second set also consists of four experiments but to test the four iterative procedures on the sub-synoptic scale by setting L = 1000 km and x= 0 for

    ψ

    in Eq. (15). The third(or fourth) set still consists of four experiments to test the four iterative procedures on the meso-

    α

    scale by setting L =500 km and x= 0 (or L) for

    ψ

    in Eq. (15). Note that setting x= 0 (or L) places the ridge (or trough) of the wavering jet in the middle of domain D, so the nonlinearly balanced flow used for the tests in the third (or fourth) set is curved anti-cyclonically (or cyclonically) in the middle of domain D. For simplicity, the Coriolis parameter f is assumed to be constant and set to f = f= 10sin all the experiments. The Rossby number, defined by Ro = U/fL, is thus 0.1, 0.2 and 0.4 for L = 2000, 1000 and 500 km, respectively.The true geopotential field,

    ?

    , is obtained by solving the Poisson equation, ?

    ?

    = N(

    ψ

    ), numerically on a 51×51 grid over domain D with the boundary condition given by

    ?

    =f

    ψ

    . In this case,

    ψ

    in Eq. (15) is also discretized on the same 51×51 grid over the same square domain, and is used to compute the right-hand side of ?

    ?

    = N(

    ψ

    ) via standard finite-differencing. Then,

    ?

    is solved numerically by using the Poisson solver of Cao and Xu (2011). The SOR scheme can be also used to solve for

    ?

    , but the solution is generally less accurate than that obtained by using the Poisson solver.The NBE discretization error (scaled by ||

    ?

    ?

    ||′) can be denoted and defined by

    This error is 3.25×10(or 4.33×10) for

    ?

    obtained by using the Poisson solver with L = 2000 (or 1000) km but increases to 5.58×10(or 5.78×10) for

    ?

    obtained by using the SOR scheme. Thus, the solution obtained by using the Poisson solver is used as the input field of

    ?

    in the NBE to test the iterative procedures in each set of experiments.

    4. Results of experiments

    4.1. Results from first set of experiments

    For this set of experiments,

    ψ

    and (u, v) are plotted in Fig. 1a,

    ψ

    and (u, v) ≡ (? ?ψ, ??ψ) are plotted in Fig.1b, the vorticity

    ζ

    ≡ ?

    ψ

    is plotted in Fig. 1c, and the geostrophic vorticity

    ζ

    ≡ ?

    ψ

    is plotted in Fig. 1d. Figure 1c shows that the absolute vorticity, defined by f +

    ζ

    , is positive everywhere so the nonlinearly balanced wavering jet flow is inertially stable over the entire domain (see the proof in Appendix C of Xu, 1994). Figure 1c also shows that the geostrophic vorticity

    ζ

    is larger than ?f/2 (= ?f/2)everywhere, so the NBE is elliptic over the entire domain and its associated boundary value problem in Eq. (1) is well posed.The relative error of

    ψ

    with respect to

    ψ

    can be denoted and defined by

    Fig. 1. (a) ψt plotted by color contours every 4.0 in the unit of 106 m2 s?1 and (ut, vt) plotted by black arrows over domain D ≡ [?L ≤ x ≤ L, ?L ≤ y ≤ L] with L = 2000 km for the first set of experiments. (b) As in (a) but for ψg and(ug, vg) with ψg ≡ ?/f and ? computed from ψt by setting f = f0 = 10?4 s?1 as described in section 3.3. (c) Vorticity ζt ≡?2ψtplotted by color contours every 0.1 in the unit of 10?4 s?1 over domain D. (d) As in (c) but for geostrophic vorticity ζg ≡?2ψg. The wavering jet axis is along the green contour of ψt = 0 in (a) with its ridge at x = 0 and two troughs at x = ±L on the west and east boundaries of domain D.

    where || || denotes the RMS of discretized field of the variable inside || || computed over all the grid points (including the boundary points) of domain D. The accuracy of the solution

    ψ

    obtained during the iterative process in each experiment can be evaluated by the RMS of

    ε

    (

    ψ

    ), denoted and defined by

    where || || is defined in Eq. (18). The accuracy to which the NBE is satisfied by

    ψ

    can be measured by E[N(

    ψ

    )] defined in Eq. (14).Table 1 lists the values of E(

    ψ

    ) and E[N(

    ψ

    )] for the initial guess

    ψ

    (=

    ψ

    ) in row 1 and the optimally truncated solutions

    ψ

    from the four experiments in rows 2?5. As shown in row 2 versus row 1 of Table 1, M1a reaches the optimal truncation at k = K = 6 where E[N(

    ψ

    )] is reduced (from 0.120 at k = 0) to its minimum [= 2.411×10< E(

    ?

    ?

    ) =3.25×10–the NBE discretization error defined in Eq.(17)] with E(

    ψ

    ) reduced (from 2.43×10at k = 0) to 4.87×10. Figure 2a shows that E(

    ψ

    ) reaches its minimum(= 4.79×10) at k = 10. This minimum is slightly below E(

    ψ

    ) = 4.87×10but undetectable in real-case applications.

    Table 1. Values of E(ψk) and E[N(ψk)] listed in row 1 for the initial guess ψ0 (= ψg) with k = 0 and in rows 2?5 for ψK from the four iterative procedures in the first set of experiments (with Ro =0.1). Here, E(ψk) is defined in Eq. (19), E[N(ψk)] is defined in Eq.(14), k is the iteration number, and ψK is the optimally truncated solution at k = K.

    Fig. 2. (a) E[N(ψk)] and E(ψk) from M1a in the first set of experiments plotted by red and blue curves, respectively,as functions of k over the range of 1 ≤ k ≤ 20. (b) As in (a) but from M1b plotted over the range of 1 ≤ k ≤ 4×104.(c) As in (a) but from M2a plotted over the range of 1 ≤ k ≤ 60. (d) As in (c) but from M2b. In each panel, the ordinate of E[N(ψk)] is on the left side labeled in red and the ordinate of E(ψk) is on the right side labeled in blue.

    On the contrary, as shown in row 3 of Table 1 and Fig.2b, M1b reaches the optimal truncation very slowly at k =K = 38493 where E[N(

    ψ

    )] is reduced to its global minimum (= 1.81×10) with E(

    ψ

    ) reduced to 1.68×10. Here,E[N(

    ψ

    )] has three extremely shallow and small local minima (at k = 32408, 38490 and 38497) not visible in Fig. 2b.These local minima are detected and passed by setting m =3 in Eq. (14) for M1b. Clearly M1b is less accurate and much less efficient than M1a.Figure 2c (or 2d) shows that M2a (or M2b) reaches the optimal truncation at k = K = 19 (or 26) where E[N(

    ψ

    )] is reduced to its global minimum [= 3.55×10(or 2.66×10)]with E(

    ψ

    ) reduced to 4.55×10(or 2.69×10), and E(

    ψ

    )decreases continuously toward its minimum [= 2.45×10(or 1.62×10)] as k increases beyond K. Thus, M2a and M2b are less efficient and much less accurate than M1a for Ro = 0.1.

    4.2. Results from second set of experiments

    For this set of experiments,

    ψ

    and (u, v) have the same patterns as those in Fig. 1a, and

    ψ

    and (u, v) are similar to those in Fig. 1b, but the contour intervals of

    ψ

    and

    ψ

    are reduced by 50% as L is reduced from 2000 to 1000 km with Ro increased to 0.2, so the wavering jet flow is on the subsynoptic scale. In this case, the nonlinearly balanced jet flow is still inertially stable over the entire domain since

    ζ

    >?f everywhere as shown in Fig. 3a, but

    ζ

    < ?f/2 in the two small yellow colored areas as shown in Fig. 3b, so the NBE becomes hyperbolic locally in this small area and the boundary value problem in Eq. (1) is not fully well posed.In this case, as shown in row 2 versus row 1 of Table 2,M1a reaches the optimal truncation at k = K = 13 where E[N(

    ψ

    )] is reduced (from 0.243 at k = 0) to its minimum [=5.23×10close to E(

    ?

    ?

    ) = 4.33×10] with E(

    ψ

    ) reduced(from 4.86×10at k = 0) to 1.24×10. The rapid descending processes of E(

    ψ

    ) and E[N(

    ψ

    )] (not shown) are similar to those in Fig. 2a for M1a in the first set of experiments.As shown in row 3 of Table 2, M1b takes K = 48057 iterations to reach the optimal truncation and the values of E[N(

    ψ

    )] and E(

    ψ

    ) at k = K are about four times larger than those from M1a. The extremely slow descending processes of E(

    ψ

    ) and E[N(

    ψ

    )] (not shown) are similar to those in Fig. 2b for M1b in the first set of experiments. As shown in row 4 (or 5) of Table 2, M2a (or M2b) reaches the optimal truncation at k = K = 26 (or 35) and the values of E[N(

    ψ

    )]and E(

    ψ

    ) are more than (or about) 4 times of those from M1a. Thus, M1a is still more accurate and much more efficient than M1b and is more efficient and much more accurate than M2a and M2b for Ro = 0.2, although the boundary value problem in Eq. (1) in this case is not fully (but nearly)well posed.

    4.3. Results from third set of experiments

    For this set of experiments,

    ψ

    and (u, v) have the same patterns as those in Fig. 1a but the contour interval of

    ψ

    is reduced 4 times as L is reduced from 2000 to 500 km with Ro increased to 0.4, so the wavering jet flow is on the meso-

    α

    scale. Figure 4a shows the fields of

    ψ

    and (u, v)for the nonlinearly balanced jet flow. This nonlinearly balanced jet flow is inertially unstable in the yellow colored area south of the ridge of wavering jet axis in the middle of domain D where

    ζ

    < ?f as shown in Fig. 4c. Figure 4d shows that

    ζ

    < ?f/2 in the long and broad yellow colored area along and around the wavering jet, so the NBE is hyperbolic in this area and the boundary value problem in Eq. (1)becomes seriously ill-posed.

    Fig. 3. (a) ζt plotted by color contours every 0.25 in the unit of 10?4 s?1 in domain D with L = 1000 km and Ro = 0.2 for the second set of experiments. (b) As in (a) but for ζg. As shown in (b), ζg < ?f/2 (= ?f0/2) in the two small yellow colored areas where the NBE becomes locally hyperbolic.

    Table 2. As in Table 1 but for the second set of experiments(with Ro = 0.2).

    In this case, as shown in row 2 of Table 3 and Fig. 5a,M1a reaches the optimal truncation at k = K = 2 where E[N(

    ψ

    )] is decreased (from 0.57 at k = 0) to its minimum(= 0.13), while E(

    ψ

    ) decreases from 9.72×10at k = 0 to 8.20×10at k = K = 2 and then to its minimum (=7.38×10) at k = 6. As k increases beyond 6, M1a diverges.Its optimally truncated solution

    ψ

    is merely slightly more accurate than the initial guess

    ψ

    . As shown in row 3 of Table 3 and Fig. 5b, M1b reaches the optimal truncation at k = K = 10325 where E[N(

    ψ

    )] is decreased to its global minimum (= 0.15), while E(

    ψ

    ) decreases to 8.31×10at k = K and then to its minimum (= 7.68×10) at k = 23515. Thus,M1b is still less accurate and much efficient than M1a.

    Fig. 4. (a) ψg plotted by color contours every 1.0 in the unit of 106 m2 s?1 and (ug, vg) plotted by black arrows over domain D with L = 500 km and Ro = 0.4 for the third set of experiments. (b) As in (a) but for ε(ψ0) = ε(ψg) plotted by color contours every 5.0 in the unit of 10?2. (c) As in (a) but for ζt plotted by color contours every 0.5 in the unit of 10?4 s?1 in domain D. (d)As in (c) but for ζg. As shown in (c), ζt < ?f in the yellow colored area south of the ridge of wavering jet axis where the jet flow becomes inertially unstable. As shown in (c), ζg < ?f/2 (= ?f0/2) in the long and broad yellow colored area (along and around the wavering jet) where the NBE becomes hyperbolic.

    Table 3. As in Table 1 but for the third set of experiments (with Ro = 0.4 and x0 = 0).

    Figure 5c (or 5d) shows that M2a (or M2b) reaches the optimal truncation at k = K = 26 (or 29) where E[N(

    ψ

    )] is reduced to its minimum [= 0.11 (or 0.10)], while E(

    ψ

    ) is reduced to its minimum [= 8.24×10(or 8.24×10)] at k =25 (or 26) and then increases slightly to 8.25×10(or 8.26×10) at k = K = 26 (or 29). As shown in row 4 (or 5)versus row 2 of Table 3, E(

    ψ

    ) from M2a (or M2b) is larger than that from M1a, so M2a (or M2b) is still less accurate than M1a in this case.

    Fig. 5. (a) E[N(ψk)] and E(ψk) from M1a in the third set of experiments plotted by red and blue curves, respectively, as functions of k over the range of 1 ≤ k ≤ 8, (b) As in (a) but from M1b plotted over the range of 1 ≤ k ≤ 3×104. (c) As in (a)but from M2a plotted over the range of 1 ≤ k ≤ 60. (d) As in (c) but from M2b. In each panel, the ordinates of E[N(ψk)] and E(ψk) are placed and labeled as in Fig. 2.

    Fig. 6. ε(ψK) plotted by color contours every 0.5 in the unit of 10?2 for ψK from (a) M1a, (b) M1b, (c) M2a and (d) M2b in the third set of experiments.

    Figure 6a (or 6b) shows that

    ε

    (

    ψ

    ) from M1a (or M1b)peaks positively and negatively in the middle of domain D as

    ε

    (

    ψ

    ) does in Fig. 4b but with slightly reduced amplitudes. Figure 6c (or 6d) shows that

    ε

    (

    ψ

    ) from M2a (or M2b) has a broad negative peak south of the ridge of wavering jet axis similar to that of

    ε

    (

    ψ

    ) in Fig. 4b but with a slightly enhanced amplitude. In this case, M1a is still slightly more accurate than the other three iterative procedures, but it cannot effectively reduce the solution error in the central part of the domain where not only is NBE hyperbolic (with

    ζ

    < ?f/2 as shown in Fig. 4d), but also the jet flow is strongly anti-cyclonically curved and subject to inertial instability (with

    ζ

    < ?f as shown in Fig. 4c).

    4.4. Results from fourth set of experiments

    For this set of experiments,

    ψ

    and (u, v) are plotted in Fig. 7a. These fields represent the same nonlinearly balanced wavering westerly jet flow as that in the third set of experiments except that the wave fields are shifted by a half-wavelength so the jet flow is curved cyclonically in the middle of domain D. In this case,

    ψ

    and (u, v) are nearly the same as the half-wavelength shifted fields (not shown)from Fig. 4a but with small differences, mainly along and around the trough and ridge lines due to the boundary condition,

    ?

    ≡ f

    ψ

    = f

    ψ

    , used here along the two trough lines(instead of the two ridge lines in Fig. 4a) for solving

    ?

    from?

    ?

    = N(

    ψ

    ). Figure 7c shows that the jet flow becomes inertially unstable in the two yellow colored areas (where

    ζ

    ζ

    < ?f/2) that is nearly the same as the yellow colored area in Fig. 4d but half-wavelength shifted, so the area of

    ζ

    < ?f (that is, the area of

    ζ

    + f < 0 in which the initial guess field is inertially unstable) in Fig. 4d is moved with the ridge line to the west and east boundaries in Fig. 7d. As the area of

    ζ

    < ?f and area of

    ζ

    < ?f are moved away from the domain center to the domain boundaries where

    ψ

    is known and given by

    ?

    /f,the NBE becomes less difficult to solve in this fourth set of experiments than in the third set.

    Fig. 7. (a) As in Fig. 4a but for ψt and (ut, vt) in the fourth set of experiments with L = 500 km and x0 = L (instead of x0 = 0).(b) As in (a) but for ε(ψ0) = ε(ψg) plotted by color contours every 6.0 in the unit of 10?2. (c) As in (a) but for ζt plotted by color contours every 0.5 in the unit of 10?4 s?1 in domain D. (d) As in (c) but for ζg.

    In this case, as shown in row 2 of Table 4 and Fig. 8a,M1a reaches the optimal truncation at k = K = 7 where E[N(

    ψ

    )] is decreased (from 0.76 at k = 0) to its minimum(= 3.81×10), while E(

    ψ

    ) decreases from 9.71×10at k =0 to 2.29×10at k = K = 7 and then to its flat minimum(= 2.25×10) at k = 12, so

    ψ

    is significantly more accurate than

    ψ

    and slightly less accurate than

    ψ

    at k = 12(which is undetectable in real-case applications). As shown in row 3 of Table 4 and Fig. 8b, M1b reaches the optimal truncation at k = K = 31830 where E[N(

    ψ

    )] is decreased to its global minimum (= 4.54×10), while E(

    ψ

    ) decreases to 2.37×10at k = K and then to its minimum (= 2.21×10)at k = 57 586. Thus, M1b is still much less efficient and less accurate than M1a.

    Table 4. As in Table 1 but for the fourth set of experiments (with Ro = 0.4 and x0 = L).

    Fig. 8. (a) E[N(ψk)] and E(ψk) from M1a in the fourth set of experiments plotted by red and blue curves, respectively, as functions of k over the range of 1 ≤ k ≤ 24, (b) As in (a) but from M1b plotted over the range of 1 ≤ k ≤ 6×104. (c) As in (a)but from M2a plotted over the range of 0 ≤ k ≤ 60. (d) As in (c) but from M2b. In each panel, the ordinates of E[N(ψk)] and E(ψk) are placed and labeled as in Fig. 2.

    Figure 8c (or 8d) shows that M2a (or M2b) reaches the optimal truncation at k = K = 27 (or 32) where E[N(

    ψ

    )] is reduced to its minimum [= 5.42×10(or 4.66×10)], E(

    ψ

    )reduces to 3.03×10(or 2.64×10) at k = K and then to its minimum [= 2.72×10(or 2.43×10)] at k = 36 (or 44), so M2a (or M2b) is still less efficient and less accurate than M1a in this case.Figure 7b shows that

    ε

    (

    ψ

    ) has a broad positive (or negative) peak south (or north) of the trough of wavering jet axis in the middle of domain D. These broad peaks are mostly reduced by M1a as shown by

    ε

    (

    ψ

    ) in Fig. 9a but slightly less reduced by M1b as shown in Fig. 9b and less reduced by M2a (or M2b) as shown in Fig. 9c (or 9d). However, the small secondary negative peak of

    ε

    (

    ψ

    ) near the west or east boundary in Fig. 7b is reduced only about 30% by M1a (or M1b) as shown by

    ε

    (

    ψ

    ) in Fig. 9a (or 9b) and even less reduced by M2a (or M2b) as shown in Fig. 9c (or 9d). Thus,all the four iterative procedures have difficulties in reducing the errors of their optimally truncated solutions near the west and east boundaries where not only is the NBE hyperbolic (with

    ζ

    < ?f/2 as shown in Fig. 7d), but also the jet flow is subject to inertial instability (with

    ζ

    < ?f as shown in Fig. 7c). Nevertheless, since the area of

    ζ

    < ?f is moved with the ridge of wavering jet axis to the domain boundaries in Fig. 7c, all of the four iterative procedures perform significantly better in this set of experiments than in the previous third set, as shown in Fig. 9 and Table 4 versus Fig. 6 and Table 3. In this case, M1a is still most accurate and M1b is still least efficient among the four iterative procedures.

    Fig. 9. ε(ψK) plotted by color contours every 2.0 in the unit of 10?2 for ψK from (a) M1a, (b) M1b, (c) M2a and (d)M2b in the fourth set of experiments.

    5. Conclusions

    In this paper, two types of existing iterative methods for solving the NBE are reviewed and revisited. The first type was originally proposed by Bolin (1955), in which the NBE is transformed into a linearized equation for a presumably small correction to the initial guess or the subsequently updated solution. The second type was originally proposed by Shuman (1955, 1957) and Miyakoda (1956), in which the NBE is rearranged into a quadratic form of the absolute vorticity and the positive root of this quadratic form is used in the form of a Poisson equation to obtain the solution iteratively. These two types of methods are rederived formally by expanding the solution asymptotically upon a small Rossby number (see section 2), and the rederived methods are called method-1 and method-2, respectively.

    Since the rearranged asymptotic expansion is not ensured to converge, especially when the Rossby number is not sufficiently small, a criterion for optimal truncation of asymptotic expansion is proposed [see Eq. (14)] to obtain the super-asymptotic approximation of the solution based on the heuristic theory of asymptotic analysis (Boyd, 1999).In addition, the Poisson solver based on the integral formulas (Cao and Xu, 2011; Xu et al., 2011) is used versus the SOR scheme to solve the boundary value problem in each iterative step.

    The four iterative procedures are tested with analytically formulated wavering jet flows on different spatial scales in four sets of experiments. The computational domain covers one full wavelength and is centered at the ridge of the wavering jet in the first three sets of experiments but centered at the trough in the last set. In the first set of experiments, the wavering jet flow is formulated on the synoptic scale [with the half wavelength L = 2000 km and the associated Rossby number Ro = 0.1]. In this case,the NBE is of the elliptic type over the entire domain and therefore its boundary value problem is well posed. In the second set of experiments, the wavering jet flow is formulated on the sub-synoptic scale [with L = 1000 km and Ro =0.2]. In this case, the NBE is of the elliptic type over nearly the entire domain so that its boundary value problem is nearly well posed. In the third (and fourth) sets of experiments, the wavering jet flow is formulated on the meso-

    α

    scale with Ro = 0.4, and the wavering jet flow is curved anti-cyclonically (or cyclonically) in the middle of the domain where the absolute vorticity is locally negative (or strongly positive). In this case, the NBE becomes hyperbolic broadly along and around the wavering jet so that its boundary value problem is seriously ill-posed.The test results can be summarized as follows: For wavering jet flows on the synoptic and sub-synoptic scales, all four iterative procedures can reach their respective optimal truncations and the solution error (originally from the initial guess–the geostrophic streamfunction) can be reduced at the optimal truncation by an order of magnitude or nearly so even when the NBE is not entirely elliptic. Among the four iterative procedures, M1a is most accurate and efficient while M1b is least efficient. The results for wavering jet flows on the synoptic and sub-synoptic scales are insensitive to the location of the wavering jet in the computational domain. In particular, according to our additional experiments (not shown in this paper), when the wavering jet is shifted zonally by a half of wavelength (with the trough moved to the domain center), the solution errors become slightly smaller and the optimal truncation numbers for M1a and M1b (or M2a and M2b) become slightly smaller (or larger)than those listed in Tables 1 and 2. For wavering jet flows on the meso-

    α

    scale in which the NBE’s boundary value problem is seriously ill-posed, the four iterative procedures still can reach their respective optimal truncations with the solution error reduced effectively for cyclonically curved part of the wavering jet flow but not for the anti-cyclonically curved part. In this case, M1a is still most accurate and efficient while M1b is least efficient.

    In comparison with M1b, the high accuracy and efficiency of M1a can be explained by the fact that the solution obtained by the Poisson solver based on the integral formulas is not only more accurate but also smoother than the solution obtained by the SOR scheme in each step of nonlinear iteration. Consequently, in each next step, the nonlinear differential term on the right-hand side of the incremental-form iteration equation [see Eq. (9a)] is computed more accurately in M1a than in M1b and so is the entire right-hand side. This is especially true and important when the entire right-hand side becomes very small (toward zero) in the late stage of iterations, as it also explains why M1b reaches the optimal truncation much slower than M1a (see Tables 1?4). In comparison with M2a and M2b, the high accuracy and efficiency of M1a can be explained by the fact that the solution in M1a is updated incrementally and the increment is small relative to the entire solution, and so is the error of the increment computed in each step of nonlinear iteration. On the other hand,the solution in M2a or M2b is updated entirely and the entire solution is large relative to the increment and so is the error of the entire solution computed in each step of nonlinear iteration. Moreover, the recursive form of equation [see Eq. (13)] used by M2a and M2b contains a square root term on its right-hand side, so it cannot be converted into an incremental form. Furthermore, this square root term must set to zero when the term inside the square root becomes negative,although the term inside the square root corresponds to the squared absolute vorticity. This problem is caused by the non-negative absolute vorticity assumed in the derivation of the recursive form of equation for M2a and M2b.

    Cyclonically curved meso-

    α

    scale jet flows in the middle and upper troposphere are often precursors of severe weather especially when the curved jet flow evolves into a cut-off cyclone atop a meso-

    α

    scale low pressure system in the lower troposphere. In this case, M1a can be potentially and particularly useful for severe weather analyses in the context of semi-balanced dynamics (Xu, 1994; Xu and Cao,2012). In addition, since the mass fields can be estimated from Advanced Microwave Sounding Unit (AMSU) observations, using the NBE to retrieve the horizontal winds in and around tropical cyclones (TC) from the estimated mass fields have potentially important applications for TC warnings and improving TC initial conditions in numerical predictions (Velden and Smith, 1983; Bessho et al, 2006). Applications of M1a in the aforementioned directions deserve continued studies. In particular, the gradient wind can be easily computed for the axisymmetric part of a cutoff cyclone (or TC)and used to improve the initial guess for the iterative procedure. This use of gradient wind can be somewhat similar to the use of gradient wind associated with the axisymmetric part of a hurricane to improve the basic-state potential vorticity (PV) construction for hurricane PV diagnoses (Wang and Zhang, 2003; Kieu and Zhang, 2010). Furthermore,either the gradient wind or the optimal truncated solution from M1a can be used as a new improved initial guess. In this case, the asymptotic expansion can be reformulated upon a new small parameter associated with the reduced error of the new initial guess and this new small parameter can be smaller or much smaller than the Rossby number used for the asymptotic expansion in this paper. The reformulated asymptotic expansion may be truncated to yield a more accurate “hyper-asymptotic” approximation of the solution according to the heuristic theory of asymptotic analysis(see section 5 of Boyd, 1999). This approach deserves further explorations.The authors are thankful to Dr. Ming XUE for reviewing the original manuscript and to the anonymous reviewers for their constructive comments and suggestions. This work was supported by the NSF of China Grants 91937301 and 41675060, the National Key Scientific and Technological Infrastructure Project "EarthLab", and the ONR Grants N000141712375 and N000142012449 to the University of Oklahoma (OU). The numerical experiments were performed at the OU supercomputer Schooner. Funding was also provided to CIMMS by NOAA/Office of Oceanic and Atmospheric Research under NOAA-OU Cooperative Agreement #NA110AR4320072, U.S.Department of Commerce.

    老司机福利观看| 国产中年淑女户外野战色| 综合色av麻豆| 日韩av不卡免费在线播放| 国产午夜精品一二区理论片| 麻豆乱淫一区二区| 91在线精品国自产拍蜜月| 日韩精品有码人妻一区| 波野结衣二区三区在线| 国产又黄又爽又无遮挡在线| 老司机影院成人| 亚洲无线在线观看| 国产午夜精品久久久久久一区二区三区| 内射极品少妇av片p| 成人亚洲精品av一区二区| 免费人成视频x8x8入口观看| 久久精品国产清高在天天线| 成年女人看的毛片在线观看| 国内精品久久久久精免费| 国产精品久久久久久精品电影小说 | 国产一区二区在线观看日韩| 白带黄色成豆腐渣| 又黄又爽又刺激的免费视频.| 人妻少妇偷人精品九色| or卡值多少钱| 国产 一区精品| 国产老妇女一区| 亚洲在线自拍视频| 久久精品国产亚洲网站| 日韩av在线大香蕉| 美女黄网站色视频| 国产三级中文精品| 亚洲欧美日韩东京热| 国产综合懂色| 美女 人体艺术 gogo| 成人午夜高清在线视频| 国产精品一及| 日韩一本色道免费dvd| 国产精品乱码一区二三区的特点| 日本在线视频免费播放| av女优亚洲男人天堂| 免费大片18禁| 99热网站在线观看| av在线蜜桃| 国产亚洲精品av在线| 精品99又大又爽又粗少妇毛片| 丝袜喷水一区| 熟妇人妻久久中文字幕3abv| 亚洲最大成人av| 大型黄色视频在线免费观看| 亚洲精品自拍成人| 日韩精品有码人妻一区| 成人特级av手机在线观看| 欧美一区二区国产精品久久精品| 深夜精品福利| 搡女人真爽免费视频火全软件| 精品少妇黑人巨大在线播放 | 男人的好看免费观看在线视频| 久久久久久久午夜电影| 男女下面进入的视频免费午夜| 日本-黄色视频高清免费观看| av在线亚洲专区| 久久久a久久爽久久v久久| 免费在线观看成人毛片| 国产成人a∨麻豆精品| 久久精品综合一区二区三区| 亚洲美女视频黄频| 波多野结衣高清无吗| 国产高清不卡午夜福利| 亚洲欧美成人精品一区二区| 高清在线视频一区二区三区 | 久久人人爽人人爽人人片va| 美女大奶头视频| 我的老师免费观看完整版| 国内精品一区二区在线观看| 欧美最黄视频在线播放免费| 成熟少妇高潮喷水视频| 国产真实乱freesex| 99久久无色码亚洲精品果冻| 精品久久久久久久久av| 搡女人真爽免费视频火全软件| 日韩精品青青久久久久久| 国产精品.久久久| 在线国产一区二区在线| 精品人妻一区二区三区麻豆| 青春草国产在线视频 | 内地一区二区视频在线| 国产精品美女特级片免费视频播放器| 日韩av在线大香蕉| 我要看日韩黄色一级片| 爱豆传媒免费全集在线观看| 97人妻精品一区二区三区麻豆| 人妻夜夜爽99麻豆av| h日本视频在线播放| 亚洲国产日韩欧美精品在线观看| 国产 一区 欧美 日韩| 蜜桃亚洲精品一区二区三区| 乱系列少妇在线播放| 日韩av在线大香蕉| 嫩草影院精品99| 精品午夜福利在线看| 少妇丰满av| 2021天堂中文幕一二区在线观| 精品一区二区三区人妻视频| 国语自产精品视频在线第100页| 一级毛片aaaaaa免费看小| 国产精品永久免费网站| 国产单亲对白刺激| 少妇人妻精品综合一区二区 | 中文字幕人妻熟人妻熟丝袜美| 干丝袜人妻中文字幕| 麻豆国产av国片精品| 丝袜美腿在线中文| 欧美一区二区亚洲| 国产一区二区亚洲精品在线观看| 在线免费观看的www视频| 日韩精品有码人妻一区| 亚洲成人精品中文字幕电影| 日韩人妻高清精品专区| av.在线天堂| 国产精品.久久久| 国产免费男女视频| 18+在线观看网站| 成人欧美大片| 中文字幕久久专区| 国产在视频线在精品| 成年av动漫网址| 嘟嘟电影网在线观看| 日本黄色片子视频| 最近2019中文字幕mv第一页| 美女cb高潮喷水在线观看| 亚洲五月天丁香| 91久久精品国产一区二区成人| 舔av片在线| a级毛片免费高清观看在线播放| 波多野结衣高清作品| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 欧美潮喷喷水| 亚洲欧美日韩高清在线视频| 欧洲精品卡2卡3卡4卡5卡区| 久久草成人影院| 日韩成人av中文字幕在线观看| 精品久久久久久久久av| 亚洲熟妇中文字幕五十中出| 日韩欧美三级三区| 久久久久久久久大av| 国产伦精品一区二区三区视频9| 国产午夜精品久久久久久一区二区三区| 波多野结衣高清无吗| ponron亚洲| 亚洲中文字幕一区二区三区有码在线看| 在线观看66精品国产| 亚洲经典国产精华液单| 少妇熟女aⅴ在线视频| 搞女人的毛片| 欧美激情国产日韩精品一区| 国产成人一区二区在线| 国产激情偷乱视频一区二区| 可以在线观看毛片的网站| 青春草视频在线免费观看| 黄色欧美视频在线观看| 日韩大尺度精品在线看网址| 欧美一级a爱片免费观看看| 国产色爽女视频免费观看| 天堂影院成人在线观看| 99久久人妻综合| 直男gayav资源| 日本在线视频免费播放| 高清在线视频一区二区三区 | 亚洲av熟女| 色播亚洲综合网| 男女视频在线观看网站免费| 免费搜索国产男女视频| 亚洲精品影视一区二区三区av| 夜夜爽天天搞| 99久久无色码亚洲精品果冻| 99久久无色码亚洲精品果冻| 嫩草影院入口| 性色avwww在线观看| 久久久久九九精品影院| 一级毛片我不卡| 禁无遮挡网站| 国产片特级美女逼逼视频| 麻豆国产av国片精品| 人妻少妇偷人精品九色| 欧美xxxx性猛交bbbb| 少妇熟女欧美另类| 色视频www国产| 成人鲁丝片一二三区免费| 国产单亲对白刺激| 中文字幕免费在线视频6| 成年av动漫网址| 国产在视频线在精品| 精品久久久久久久久久免费视频| 欧美一级a爱片免费观看看| 99精品在免费线老司机午夜| 国产在线男女| 亚洲性久久影院| 少妇熟女aⅴ在线视频| 婷婷精品国产亚洲av| 亚洲人成网站在线播放欧美日韩| 成年女人看的毛片在线观看| 国产高潮美女av| 人人妻人人澡欧美一区二区| 日韩欧美精品v在线| 日本色播在线视频| 日韩制服骚丝袜av| 99热只有精品国产| 国产黄a三级三级三级人| 久久久精品欧美日韩精品| 校园人妻丝袜中文字幕| 欧美激情在线99| 一区二区三区免费毛片| 高清毛片免费观看视频网站| 高清毛片免费看| 亚洲无线观看免费| 天天一区二区日本电影三级| 成年女人看的毛片在线观看| 成人毛片60女人毛片免费| 又粗又爽又猛毛片免费看| 亚洲国产精品sss在线观看| 国内揄拍国产精品人妻在线| 久久中文看片网| 久久久久久久午夜电影| 中国国产av一级| 99久久无色码亚洲精品果冻| 午夜福利视频1000在线观看| av女优亚洲男人天堂| 成人二区视频| 欧洲精品卡2卡3卡4卡5卡区| 国产色爽女视频免费观看| 少妇的逼水好多| 99久久九九国产精品国产免费| 成人特级黄色片久久久久久久| 22中文网久久字幕| 深夜a级毛片| 色综合亚洲欧美另类图片| 亚洲熟妇中文字幕五十中出| 18禁在线播放成人免费| 久久午夜福利片| 亚洲国产色片| 国产69精品久久久久777片| 日本三级黄在线观看| 国产又黄又爽又无遮挡在线| 长腿黑丝高跟| 久久久久久久午夜电影| 欧美激情在线99| 老司机影院成人| 婷婷色av中文字幕| 最近最新中文字幕大全电影3| 久久久欧美国产精品| 欧美区成人在线视频| 国产精品福利在线免费观看| 99热网站在线观看| 综合色av麻豆| 久久韩国三级中文字幕| 精品欧美国产一区二区三| 久久久国产成人精品二区| 中文字幕人妻熟人妻熟丝袜美| 午夜视频国产福利| www.av在线官网国产| 国产免费一级a男人的天堂| www日本黄色视频网| 长腿黑丝高跟| 一区二区三区四区激情视频 | 亚洲,欧美,日韩| 精品熟女少妇av免费看| 中文欧美无线码| 免费不卡的大黄色大毛片视频在线观看 | 99久久无色码亚洲精品果冻| 人妻少妇偷人精品九色| 亚洲久久久久久中文字幕| 国产激情偷乱视频一区二区| 亚洲成人av在线免费| 免费看a级黄色片| 欧美精品国产亚洲| а√天堂www在线а√下载| 1024手机看黄色片| 精品99又大又爽又粗少妇毛片| 亚洲丝袜综合中文字幕| 国产白丝娇喘喷水9色精品| 日韩一本色道免费dvd| 九九久久精品国产亚洲av麻豆| 淫秽高清视频在线观看| 在线观看午夜福利视频| 亚洲成a人片在线一区二区| 丝袜喷水一区| 色综合色国产| 午夜福利高清视频| 丰满的人妻完整版| 少妇裸体淫交视频免费看高清| 久久午夜亚洲精品久久| 校园春色视频在线观看| 久久久欧美国产精品| 亚洲国产高清在线一区二区三| 亚洲国产精品sss在线观看| 国产国拍精品亚洲av在线观看| 少妇猛男粗大的猛烈进出视频 | 国产人妻一区二区三区在| 特大巨黑吊av在线直播| 麻豆一二三区av精品| 免费av不卡在线播放| 久久人人精品亚洲av| 亚洲内射少妇av| 久久韩国三级中文字幕| 国产日本99.免费观看| 国产成人精品一,二区 | 免费大片18禁| 国产淫片久久久久久久久| 亚洲精品乱码久久久v下载方式| 久久亚洲国产成人精品v| 国内精品美女久久久久久| 成人三级黄色视频| 成人特级av手机在线观看| 最近手机中文字幕大全| 中文字幕制服av| 欧美日本视频| 国产午夜福利久久久久久| 久久精品夜夜夜夜夜久久蜜豆| 欧美性猛交╳xxx乱大交人| av女优亚洲男人天堂| 国产又黄又爽又无遮挡在线| www.av在线官网国产| 国产亚洲5aaaaa淫片| 久久精品影院6| 午夜免费激情av| 高清毛片免费看| 国产精品99久久久久久久久| 欧美一区二区精品小视频在线| 我的女老师完整版在线观看| 99国产精品一区二区蜜桃av| 久久久久性生活片| 可以在线观看毛片的网站| 全区人妻精品视频| 最新中文字幕久久久久| 免费看美女性在线毛片视频| 51国产日韩欧美| 欧美最黄视频在线播放免费| 99视频精品全部免费 在线| 欧美日韩精品成人综合77777| 听说在线观看完整版免费高清| 国产精品美女特级片免费视频播放器| 最近最新中文字幕大全电影3| h日本视频在线播放| 亚洲精品成人久久久久久| 黄色视频,在线免费观看| 国产美女午夜福利| 国产伦精品一区二区三区四那| 国产一级毛片七仙女欲春2| 国产成人91sexporn| 久久韩国三级中文字幕| 国产精品蜜桃在线观看 | 悠悠久久av| 亚洲国产欧洲综合997久久,| 性色avwww在线观看| 日本在线视频免费播放| 国语自产精品视频在线第100页| 欧美日韩精品成人综合77777| 国内少妇人妻偷人精品xxx网站| 国产伦精品一区二区三区视频9| 最近中文字幕高清免费大全6| 女人被狂操c到高潮| 欧美日本视频| 日日啪夜夜撸| 久久久久久大精品| 简卡轻食公司| 99精品在免费线老司机午夜| 色哟哟哟哟哟哟| 国产成人精品婷婷| 亚洲人与动物交配视频| kizo精华| 国产精品,欧美在线| 日本黄色片子视频| av福利片在线观看| 亚洲综合色惰| 美女国产视频在线观看| 69av精品久久久久久| 又爽又黄无遮挡网站| 久久鲁丝午夜福利片| 亚洲美女视频黄频| 午夜免费激情av| 99精品在免费线老司机午夜| 亚洲激情五月婷婷啪啪| 一级毛片久久久久久久久女| 欧美成人免费av一区二区三区| 亚洲av二区三区四区| 亚洲av免费在线观看| 午夜老司机福利剧场| 日韩欧美国产在线观看| 欧美最新免费一区二区三区| 国产私拍福利视频在线观看| 久久久久网色| 日韩欧美在线乱码| 91aial.com中文字幕在线观看| 69av精品久久久久久| 老师上课跳d突然被开到最大视频| 亚洲一区二区三区色噜噜| 综合色av麻豆| 欧美成人精品欧美一级黄| 在线a可以看的网站| 亚洲av成人av| a级毛色黄片| 国产精品.久久久| av又黄又爽大尺度在线免费看 | av女优亚洲男人天堂| 三级经典国产精品| 爱豆传媒免费全集在线观看| 国产单亲对白刺激| av.在线天堂| 久久久久久大精品| 乱码一卡2卡4卡精品| 国产熟女欧美一区二区| 欧美日韩精品成人综合77777| 中文字幕精品亚洲无线码一区| 亚洲中文字幕一区二区三区有码在线看| 久久久精品94久久精品| 久久久精品大字幕| 国产不卡一卡二| 一本久久中文字幕| 欧美潮喷喷水| 亚洲精品久久久久久婷婷小说 | 人妻夜夜爽99麻豆av| 欧美一区二区国产精品久久精品| a级毛色黄片| 国产精品一区二区三区四区免费观看| 色综合亚洲欧美另类图片| 日本欧美国产在线视频| 国产激情偷乱视频一区二区| 2021天堂中文幕一二区在线观| 亚洲欧美日韩卡通动漫| 又爽又黄a免费视频| 1000部很黄的大片| 中国国产av一级| 99热这里只有精品一区| 在线a可以看的网站| АⅤ资源中文在线天堂| 日本撒尿小便嘘嘘汇集6| 黄色一级大片看看| 在线免费观看的www视频| 亚洲人与动物交配视频| 一个人看视频在线观看www免费| 国产蜜桃级精品一区二区三区| 春色校园在线视频观看| 成人欧美大片| 日本色播在线视频| 久久久久九九精品影院| 赤兔流量卡办理| 精品久久久久久久久亚洲| 一级毛片我不卡| 国产午夜精品一二区理论片| 欧美日韩精品成人综合77777| 男插女下体视频免费在线播放| 国产在线男女| 少妇丰满av| 老司机影院成人| 欧美+日韩+精品| 亚洲av男天堂| 热99re8久久精品国产| 亚洲欧美日韩东京热| 97超视频在线观看视频| 日本撒尿小便嘘嘘汇集6| 三级经典国产精品| 久久精品影院6| 一级av片app| 成年女人看的毛片在线观看| 中文字幕熟女人妻在线| 一级毛片aaaaaa免费看小| 欧美又色又爽又黄视频| 国产成人a∨麻豆精品| 免费人成在线观看视频色| 免费观看a级毛片全部| 国产精华一区二区三区| 国产伦理片在线播放av一区 | 亚洲欧美日韩卡通动漫| 久久久精品大字幕| 亚洲成人精品中文字幕电影| 男插女下体视频免费在线播放| 日韩在线高清观看一区二区三区| 婷婷六月久久综合丁香| 狂野欧美白嫩少妇大欣赏| 亚洲av二区三区四区| 久久久久久国产a免费观看| 午夜精品国产一区二区电影 | 国产精品电影一区二区三区| av免费观看日本| 国内精品一区二区在线观看| 免费不卡的大黄色大毛片视频在线观看 | 老女人水多毛片| 久久99精品国语久久久| 性色avwww在线观看| 91精品一卡2卡3卡4卡| 久久久久久伊人网av| a级毛色黄片| 久久精品91蜜桃| 校园人妻丝袜中文字幕| 成人美女网站在线观看视频| 美女大奶头视频| 99在线视频只有这里精品首页| 熟妇人妻久久中文字幕3abv| 五月玫瑰六月丁香| 晚上一个人看的免费电影| 亚洲成人久久性| 一个人观看的视频www高清免费观看| a级一级毛片免费在线观看| 九九久久精品国产亚洲av麻豆| 三级毛片av免费| 熟女电影av网| 91在线精品国自产拍蜜月| 成人欧美大片| 九九久久精品国产亚洲av麻豆| 日韩一区二区三区影片| av.在线天堂| 亚洲欧美日韩卡通动漫| 午夜视频国产福利| 99久久中文字幕三级久久日本| 久久久国产成人精品二区| av在线蜜桃| 亚洲婷婷狠狠爱综合网| 麻豆成人av视频| 变态另类丝袜制服| 国产久久久一区二区三区| 免费看a级黄色片| 久久久久网色| 成人亚洲欧美一区二区av| 青青草视频在线视频观看| 大又大粗又爽又黄少妇毛片口| 九九久久精品国产亚洲av麻豆| 久久99热6这里只有精品| av在线蜜桃| 在线观看午夜福利视频| 伊人久久精品亚洲午夜| 一级毛片电影观看 | 亚洲婷婷狠狠爱综合网| 国产成人一区二区在线| 色哟哟·www| av在线老鸭窝| 亚洲美女视频黄频| 久久久国产成人精品二区| 日本免费一区二区三区高清不卡| 最后的刺客免费高清国语| 久久鲁丝午夜福利片| 国产视频首页在线观看| 欧美一级a爱片免费观看看| 国产伦精品一区二区三区视频9| 色视频www国产| 大型黄色视频在线免费观看| 日本一二三区视频观看| 三级经典国产精品| 一夜夜www| 欧美日本视频| 国产一区二区激情短视频| 男女那种视频在线观看| 久久精品国产亚洲av香蕉五月| 蜜桃亚洲精品一区二区三区| 亚洲精品久久久久久婷婷小说 | 爱豆传媒免费全集在线观看| 内地一区二区视频在线| 深夜a级毛片| a级毛色黄片| 男女啪啪激烈高潮av片| 久久久久网色| 天堂中文最新版在线下载 | 99久久精品热视频| 晚上一个人看的免费电影| 人人妻人人澡欧美一区二区| 国产美女午夜福利| 悠悠久久av| 国产精品嫩草影院av在线观看| 小说图片视频综合网站| 久久韩国三级中文字幕| 九草在线视频观看| 搡老妇女老女人老熟妇| 亚洲无线观看免费| 我要看日韩黄色一级片| 中文字幕av在线有码专区| 亚洲内射少妇av| 亚洲三级黄色毛片| 亚洲国产欧洲综合997久久,| 国产高清三级在线| 91麻豆精品激情在线观看国产| 日本三级黄在线观看| 男女下面进入的视频免费午夜| 国产精品蜜桃在线观看 | 成人特级av手机在线观看| 可以在线观看毛片的网站| АⅤ资源中文在线天堂| 午夜精品一区二区三区免费看| 深夜精品福利| 天堂网av新在线| 国产成人精品婷婷| 啦啦啦观看免费观看视频高清| 国产高清视频在线观看网站| 日本一本二区三区精品| 黑人高潮一二区| 精品久久久噜噜| 久久人人爽人人爽人人片va| 日韩视频在线欧美| 亚洲av成人av| 成人欧美大片| 成人特级av手机在线观看| 亚洲人成网站在线观看播放| 国产精品三级大全| 亚洲av中文字字幕乱码综合| 九九热线精品视视频播放| 国产精品人妻久久久久久| 听说在线观看完整版免费高清| 亚洲无线观看免费| 观看美女的网站| 亚洲国产精品合色在线| 亚洲国产精品成人综合色| 亚洲av二区三区四区| 一区二区三区免费毛片| 国产精品一区二区性色av| 国产美女午夜福利| 欧美日韩在线观看h|