{'fraction_non_alphanumeric': 0.03490541435134105, 'fraction_numerical': 0.01994021946277759, 'mean_word_length': 4.594837261503928, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Orographic‐Drag Parametrization Scheme Including Orographic Anisotropy for All Flow Directions

Orographic drag is an essential process for numerical weather predictions in complex terrain regions, which depends on the inflow direction. In this study, we define the orographic asymmetry vector (OAV) for a coarse grid as the normalized vector between the grid's center point and its center of mass, and the orographic asymmetry in a flow direction—which describes the inclination direction and the extent of the subgrid terrain—is calculated as the projection of OAV on this direction. Calculations of the effective orographic length (OL) and the model grid length λ are extended to all flow directions. A new orographic drag scheme, which considers the effect of orographic anisotropy in all directions, is then developed based on the OAV projection and the extended OL and λ for any given direction. Sensitivity tests of the orographic drag under the new scheme are conducted using a 5 m/s vertically uniform wind along different directions for four coarse grid points in typical mountain regions. The new scheme is shown to provide a more continuous transition of the orographic parameters and the resulting stress as a function of flow direction than piecewise transition of schemes with only eight directions. The predicted momentum flux profile of the new scheme is compared with mountain‐wave simulations obtained from the integrated modeling system IAP‐WRF (Institute of Atmospheric Physics‐Weather Research and Forecasting Model) for the Rocky Mountain. The new scheme is shown to predict an overall narrower stress scatter about the reference simulation than the old scheme.


Introduction
The orographic drag exerted by the Earth surface on the atmospheric flow plays an important role in the general circulation of the atmosphere, and the inclusion of the effects of orographic drag in numerical weather prediction and climate models is crucial for successful weather and climate predictions (e.g., Palmer et al., 1986;Miller et al., 1989;Kim et al., 2003). Subgrid-scale orographic drag parametrizations (SOPs) that describe the vertical transport of momentum induced by subgrid-scale orography on large-scale flow include representations of the following three processes: gravity wave drag (GWD), low-level wave breaking (LLWB), and flow blocking drag (FBD) (Kim & Doyle, 2005;Lott & Miller, 1997;Zadra et al., 2003). GWD is generated by flow over mountains, which in turn generates vertical propagating waves that produce drag at upper levels when breaking. The LLWB and trapped lee waves downstream can enhance the drag in the lower troposphere. The FBD is forced by flow blocked on the mountain flanks or flowing around the mountain under upstream stable conditions, providing drag near the surface, where the blocking occurs. All three processes are strongly related to the orientation of the terrain. Initially, studies about SOPs focused mainly on parameterizing the effects of the GWD using the two-dimensional (2-D) gravity wave theory for an idealized mountain (e.g., Boer et al., 1984;McFarlane, 1987). This theory is based on the assumption of an isotropic terrain and neglects the three-dimensional (3-D) spreading of the mountain waves and the low level flow that goes around the mountain below the FBD separation line (Gregory et al., 1998;Lott & Miller, 1997;Nappo & Chimonas, 1992), thereby leading to excessive drag for all flow directions that is commonly alleviated by imposing an upper limit to the drag (Miller et al., 1989;Palmer et al., 1986). Another potential adverse impact of this assumption is the bias it can generate on orographic precipitation, which has been shown to be highly sensitive to the inflow wind direction (Hughes et al., 2014;Neiman et al., 2011;Nuss & Miller, 2001;Picard & Mass, 2017).
With regard to the incorporation of the orographic anisotropy into the orographic drag schemes for SOPs under a 3-D framework, two types of schemes exist depending on the treatment of orographic specification. The first type assumes an elliptical shape of the orography based on an analytical function with several parameters to be derived using the subgrid orographic height data. The "best fit" elliptical shape is used to derive the relevant anisotropic parameters in the flow direction used to derive the orographic drag (Baines & Palmer, 1990;Lott & Miller, 1997;Phillips, 1984;Scinocca & McFarlane, 2000). This inherently assumes a symmetrical terrain shape on the windward and leeward side of the mountain in the inflow direction and hinders the input of high order details of the terrain shape specification into the orographic parameter derivation. The second type, on the other hand, derives the anisotropic parameters in bulk form, and hence does not require an a priori assumption of the terrain shape (Kim & Doyle, 2005;Hong et al., 2008;Choi & Hong, 2015). This brings more flexibility in introducing the higher-order moments of the mountain shape-such as orographic asymmetry-for the determination of the nonlinear orographic drag in the inflow direction (Kim & Doyle, 2005;Hong et al., 2008;Choi & Hong, 2015). However, the representation of orographic anisotropy in this line of treatment is restricted only to eight representative directions. This is due to the fact that partitioning of a coarse grid for orographic parameter derivation is relatively straightforward at the specific angles defined by multiples of π/4 (e.g., ±0°, ±45°, ±90°, ±135°; see, for example, Figures 2 and 6 in Kim & Doyle, 2005) but not as easy for other angles. Consequently, this treatment may induce sudden jumps in the values of the orographic statistics with small changes in the wind direction, which introduces potential bias in weather and climate prediction as discussed in Kim and Doyle (2005).
In this study, we generalize the representation of orographic parameters and develop a new subgrid orographic drag scheme that includes orographic anisotropy for all flow directions. This scheme is validated against high-resolution mesoscale simulations. The paper is organized as follows: Section 2 describes the subgrid orographic drag scheme; section 3 presents its sensitivity to the derived parameters and section 4 presents preliminary validation against high-resolution mesoscale simulation. Lastly, a summary and conclusions are provided in section 5.

Orographic Asymmetry, Effective Orographic Length, and Model Grid Length
In the original scheme of Kim and Doyle (2005), the direction-related parameters orographic asymmetry (OA), effective orographic length (OL), model grid length λ, and the derived orographic direction (OD) were applied in eight representative flow directions. This may cause a sudden jump near the boundary of neighboring directions (Kim & Doyle, 2005). In this study, we build on the derivation method of the parameters so they can be derived in all flow directions according to the procedure described below. 2.1.1. Orographic Asymmetry OA describes the inclination extent of the subgrid terrain in the windward direction. In Kim and Doyle (2005), the subgrid orography in a model grid-box was separated into upstream and downstream sides, and the OA was derived as the ratio between the number of points in the upstream and downstream regions that were higher than the grid box average.
Revisiting the original definition of the OA by Kim and Arakawa (1995), we find that the original skewness parameter in a 2-D space is a scalar value defined as the distance between the model's center grid point and the mode, normalized by the standard deviation of the scale in the given inflow direction [Appendix B in Kim & Arakawa, 1995]. If viewing the OA as the normalized distance between two points, the concept of this parameter can be easily generalized from 2-D to 3-D space. The scalar value in the 2-D space is then turned to a vector value in the 3-D space. This "OA inclination vector (OAV)" then describes the overall inclination direction and rate of change of the subgrid topography in the 3-D space. The OA value in the given inflow direction is then derived as the projection of the OAV onto this direction. Through this way, the method can easily be applied to all flow directions. An idealized example of the OAV is depicted in Figure 1. The idealized topography is a symmetric topography about thẽ 60°direction from the east. In the cross section of the 60°slice, an inflow that blows from the southwest to the northeast direction would encounter a mountain that leans to the downstream side; thus, having a negative OA in this direction. For wind directions other than along the~60°axis, however, the relief inclination is not as steep, and the OA is derived only as the projection of the OAV on this direction.
The derivation of the OAV is determined from the position of the model's center grid point and the mode. In a model grid, the position of the model's center grid point is fixed, and the problem falls only on the mode determination. Here, the original relative position of the mode defined in Kim and Arakawa (1995) is modified as where x j and h j are the x coordinate value and the height value of the high-resolution subgrid topography in the coarse grid, respectively. N B denotes the total number of subgrid topography points in the coarse grid. Through this derivation, we view the height as being similar to the density when determining the mass center position, and thus the mode represents the "center of mass" of the subgrid mountain heights. Prior to calculation of the OAV, points lower than the grid mean average height are masked in order to retain only the mountain structure. Under this new formula, the OA is conceptualized as the vector between the model grid's center point and the grid's center of mass normalized by the standard deviation of the scale.
Derivation of the OA in the commonly used Lat-Lon grid is not straightforward since calculation of the "center of mass" coordinate is written in the Cartesian coordinate form Baron (2004). For a more general derivation of the OA, we first show the mode and the OAV under the 3-D absolute Cartesian Coordinate (the z axis in this coordinate is aligned along the Earth's axis of rotation, and the origin is at the center of the Earth; Nair et al., 2005) and then transform the vectors onto Lat-Lon grid. The derivation of the OA during this process is seen as the normalized vector between two points on the Earth's sphere in 3-D space. It is worthwhile noting that for map projections (such as the Mercator projection) that are commonly used in regional models, the OA can be seen as the normalized vector between two points on the 2-D projected plane, and thus only the (x,y) coordinates are needed to derive the OA.
The absolute Cartesian coordinates of the mode are shown as follows: where X mode , Y mode ,Z mode are the absolute Cartesian coordinates of the center of mass on the global sphere, respectively, and x i ,y i ,z i ,h i are the absolute Cartesian coordinates and height of the fine-resolution grids, respectively. The absolute Cartesian coordinates x i ,y i ,z i can be derived from the correspondent grid point latitude and longitude using equation (24) in Nair et al. (2005). The OAV under the absolute Cartesian coordinate is then derived by the following formula: where X¯; Y¯; Z¯are the coordinate of the model center grid point on the sphere, and б x ,б y ,б z are the standard deviation of the distance scale in the x, y, z directions of the absolute Cartesian coordinate, respectively.
The OAV is then transformed onto Lat-Lon coordinate following: where OA lat , OA lon are the OA value in the latitudinal and longitudinal direction, respectively, and λ lon and φ lat are the corresponding latitude and longitude radius value of the coarse grid center. The OA value OA θ along the wind direction θ is derived as a point dot product of unit direction vector (cosθ, sinθ) and OA vector (OA lon , OA lat ) as follows: Under this scheme, the OA value is therefore derived as a continuous function of the inflow direction θ.

Effective Orographic Length
The OL estimates the bulk volume of the subgrid-scale orography associated with the nonlinearity of the flow in the low-level wind direction. The original 2-D definition of OL given in the Kim and Arakawa (1995) scheme is the fractional width of the area covered by subgrid orography higher than the critical height in the flow direction (e.g., ∑ x L h /Δx, where L h is the fractional width of the subgrid topographic heights and Δx is the domain interval). The 3-D generalization of OL given by Kim and Doyle (2005) is modified here to estimate the fractional area of the subgrid topographic heights that exceed the critical height in the center half of the grid box area. The derivation of OL in 3-D space using this definition in the eight representative directions is quite simple since it is quite easy to evenly separate the center half areas in these directions [ Figure 6 in Kim & Doyle, 2005]. When deriving for other directions, however, the main problem is turned into determining the positions of two parallel boundary lines with respect to the given direction that would effectively occupy the center half area of the coarse grid ( Figure 2a). We thus derive the OL for all directions by first determining the position of the two parallel boundary lines (e.g., the position of four points intersecting the grid box), and then the fractional area that exceeds the critical height ( Figure 2).
The typical expression for points on a plane is xsinθ − ycosθ = C, where C is a constant number to be determined, and (x,y) are the coordinates of a given point (Bronshtein & Semendyayev, 1973). The above expression can be used to represent a line on a plane once θ and C are given. For a given θ, the points (determined by C) between the two parallel lines (determined by the given C 1 and C 2 ) can be decided by min(C 1 ,C 2 ) < C < max(C 1 ,C 2 ). The derivation of OL is thus the following transformed formula: and where OL θ is the OL in the inflow direction θ, h is the subgrid topographic height, and h c is the coarse grid critical height. In Kim and Arakawa (1995), h c was derived following h c = 1,116.2 − 0.878б h , where б h is the standard deviation of the subgrid topographic heights. In this study, however, we use the definition of Kim and Doyle (2005), where h c is approximated by the mean grid height. C 1 and C 2 are the corresponding C values for the two parallel boundary lines. Determining C 1 and C 2 can be roughly divided into three categories depending on the angle: when two boundaries intersect the left and right sides of the grid ( Figure 2b); when two boundaries intersect four sides of the grid ( Figure 2c); and when two boundaries intersect the upper and lower sides of the grid (Figure 2d). These conditions are summarized in the following (see Appendix A for details of the derivation): where a and b are the meridional and zonal grid length, shown in Figure 2. For angles greater than 180°(or in the third and fourth quadrants), the OL is the same as the OL modulo π, that is,

Model Grid Length λ θ
In Kim and Doyle (2005), the model grid length λ θ was derived in the zonal, meridional, and diagonal directions. In this study, we extend the derivation of λ θ for all directions. For directions in the first quadrant, the derivation is categorized into two cases (Figures 3a and 3b) depending on the value of the angle: The λ θ value in the second quadrant can be derived using symmetric pattern, and those in the third and fourth quadrants are the same as their counterpart in the first and second quadrants minus 180°: This method allows representative model grid length to be obtained in the inflow direction, which provides a smooth transition of λ θ from one direction to another.

General Framework
The scheme devised in this study (hereafter referred as "new scheme") uses a SOP based on Kim and Arakawa (1995), which includes the effects of the GWD and the drag due to LLWB and nonhydrostatic wave trapping. The gravity wave stress (τ) at the reference level h ref (determined as max(2б h ,h pbl ) following Koo et al. (2018), where h pbl is the planetary boundary layer height), as follows: where ρ 0 is the low-level density, and E is the enhancement factor representing the nonlinear enhancement of drag due to the LLWB. This factor is controlled by the OA, which represents the shape and location of the subgrid-scale orography relative to the grid and the inverse Froude number for the reference level Fr 0 (normalized by the critical inverse Froude number Fr c =1). The parameter m is the number of subgrid-scale orography, and λ eff is the model effective grid length and can be used as a tunable parameter.
G is an asymptotic function that provides a smooth transition between the nonblocking and blocking cases and includes the effects of sharpness of the mountain through orographic convexity (OC)-a parameter that represents the sharpness of the mountain and corresponds to the vertical orographic aspect ratio (calculated where h¯denotes the average of the subgrid topography in the coarse grids).
U 0 is the horizontal wind speed projected to the direction of the low-level wind, and N 0 is the Brunt-Vaisala frequency. The coefficients C E and C G are set to 0.8 and 0.5, respectively, based on the mesoscale simulation Journal of Advances in Modeling Earth Systems results from Kim and Arakawa (1995). Subscript θ denotes the direction-related parameter in the low-level wind direction, and subscript o denotes the low-level average between the surface and the reference level href .
Fr 0 is calculated according to equation (4) in Kim and Doyle (2005) as follows: and OD is derived as the ratio between the crosswind effective orographic length (OLP) and the along-wind OL, (e.g., OD=OLP θ /OL θ ).
The FBD effects are also included in the new scheme. The FBDf stress is estimated as follows: where Δ 2 x is the grid box area, C d is the bulk drag coefficient defined by C d = max{2 − 1/OD,0}, Δ p x is the grid length in the crosswind direction, and h B is the height of the blocked layer.
The vertical stress profile above the reference level is determined according to linear instability theory (Lindzen, 1981), and the nonlinear resonant LLWB is adjusted according to the Scorer parameter (Kim & Arakawa, 1995;Kim & Doyle, 2005).

The minimum Richardson number for a layer is defined by
where Ri is the mean Richardson number Ri ¼ N 2 i = U z ð Þ 2 i , N i and (U z ) i are the Brunt-Vaisala frequency and vertical wind shear for the ith level, respectively, and i denotes the vertical layer index that decreases with height. The inverse Froude number for the ith level is A model layer is considered unstable if the minimum Richardson number Ri m is less than its critical number, that is, if Ri m < Ri c . The critical Richardson number Ri c is typically 0.25 according to linear stability theory (Howard, 1961;Miles, 1961). It was, however, shown that a stratified shear flow can be unstable when Ri c is of order unity in some observations and nonlinear cases (Abarbanel & Holm, 1986;Wood, 1969). Therefore, in the later sections, we set the value of Ri c following Lott and Miller (1997), so as to provide a better upper-level momentum profile. The vertical displacement h d is expressed by combining the expressions derived by Palmer et al. (1986) and Pierrehumbert (1986) as follows: The expression for the critical value of the vertical displacement (h d ) c is derived from the following equation: The vertical stress profile is determined from the previous equations according to the following steps: 1. Calculate τ GWD using equation ( (16). 4. τ i for the upper layer is calculated from equation (15) using replaced (h d ) c .
For drag at altitudes greater than the interface height h int (defined as the level where Ri first decreases with height above the reference level), and which meets the Ri criterion set in a model layer, linear instability theory is applied to the profile; for drag at altitudes greater than h ref and lower than h int , however, the gravity wave stress τ i is adjusted under the following method: where C l is an adjustment parameter and we currently set C l =1, and l 2 i is the Scorer parameter for layer i, which can be approximated by

Sensitivity of the Scheme to Resolution of the Inflow Direction
The sensitivity of the scheme to the resolution of the inflow direction's orographic parameter is investigated under the following design: An idealized vertically uniform 5 m/s wind of varying direction that covers the whole 360°on the 1.4°× 1.4°grid is used as input condition for the parametrization scheme for a comparison experiment. The first set uses the new scheme with direction-related parameters (OA, OL, λ θ ) for only eight representative directions (e.g., ±0°, ±45°, ±90°, ±135°, 8x scheme), and the other set uses the direction-related parameters for all the directions (ALLDIR scheme). Both experiments separately consist of 360 members. As a preliminary examination, four grid points were chosen in three mountain ranges to examine the behavior of the reference level GWD and the FBD stress as a function of inflow direction, and the subgrid topography of the coarse grid (1.4°× 1.4°) for the four grid points is shown in Figure 4. The grid points are taken from the Tibetan Plateau, (28°N, 93°E), (31°N, 80°E) (Figures 4a and 4b), the Andes (12°S, 73°E) (Figure 4c), and the Rocky Mountain (40°N, 105°W) ( Figure 4d). We compare the direction-related parameters (OA, OL, λ θ ) with the GWD and the FBD stress (calculated by equations (11) and (13)) from both experiments. In this study, we set the model effective grid length λ eff as triple the value of the root sum square of the model grid length λ θ in the along-wind and cross-wind directions as it demonstrated good agreement in the comparison between the parameterization scheme and the reference simulation discussed later in section 4.  (40°N, 105°W). The x and y axes represent the west-to-east and south-to-north directions, respectively, which have been scaled to (−0.7, 0.7), and the z axis denotes the height of the topography (m). The vector represents the OA inclination vector of the subgrid topography.
The first point is a typical high-altitude point at the southeastern boundary of the Tibetan Plateau. The high mountains in this subgrid topography are mainly located in the northwestern part of the coarse grid, and thus, the orientation follows a transition from high-to-low along a northwest-southeast axis (Figure 4a). This orientation is nicely captured by the OAV, with the vector pointing toward the southeast. As a result, the OA derived from the 8x scheme (black line in Figure 5a, first row) and the ALLDIR (red line in Figure 5a, first row) both show high values in the southeast direction and low values in the northwest direction. The difference between the two schemes, however, lies in that the 8x scheme only exhibits eight-section piecewise transitions in the 360°range while the ALLDIR scheme shows a more continuous transition. The derived OA values of the two schemes thus coincide along the original representative 8x directions, but generally exhibit normalized differences as high as 0.5 in the middle directions (MDs) of two nearby representative directions (e.g., ±22.5°, ±67.5°, ±112.5°, ±157.5°). These directions are indeed farthest from the 8x directions, and therefore, commonly present the largest difference with those 8x directions.
With regard to the OL/OLP, the highest mountains in this grid locate in the northwestern edge, and thus, the OL derived from the two schemes both exhibit a higher value in the northwest and southeast direction (and the OLP with a 90°shift). Similar to the OA, the OL/OLP is most similar in the 8x directions and most different in the MDs. The highest normal difference between the two schemes' OL/OLP derived is~0.1, and the ALLDIR scheme results generally exhibit a larger variability and more continuous transition than those of the 8x scheme (Figure 5a second row). The OD is also similar to that of the OL/OLP, and its variability generally follows that of the OLP (Figure 5a, third row). The results for the model effective grid length λ θ are similar to those of OA, OL, and OD, except that the variability is generally lower than that of OA, and thus, may cause minor differences in the stress comparing to other parameters (Figure 5a, fourth row). As a result of the large variability of the 8x direction OA as the wind direction changes, the predicted stress tends to show an overall large value in the southeast direction, and a low value in the northwest direction ( Figure 5a, fifth row). In some cases, however, high values of the OA are offset by low values of the OD for the GWD and OLP for the FBD that resulted in both stress values tending to local low values. The stress determined by the ALLDIR scheme tends to show a more continuous transition, and the largest difference between the two schemes for the GWD occurs in the range 157.5-180°, where large differences in OA exist. The approximation using the nearby 180°parameters would cause an overestimation of the OD, and thus, of the GWD. As for the FBD, the largest difference occurs near −67.5°to −90°, which is caused by overall large OLP and OA differences. Using the −90°parameters would cause an underestimation of the OLP and thus of the resulting FBD. The misestimation of the stress is as high as ±0.1 N/m 2 , which is approximately 30-40% of the value of the stress calculated with the 8x scheme.
The second point is in the western part of the southern boundary of the Tibetan Plateau (31°N, 80°E) ( Figure 4b). The orientation of the subgrid topography generally follows a transition from high-to-low along the north-south axis, which is also captured by the OAV. The parameters are overall relatively similar to those of the first site, except for the lower variability of OA, and higher values of OL ( Figure 5b). The largest difference in GWD and FBD between the two schemes are seen in the −112.5°to −135°range, as a result of the large OD and OLP differences. The resulting difference in the stress is as high as 0.05 and 0.1 N/m 2 for the GWD and the FBD, respectively, or approximately 30% and 25% of the value calculated by the 8x scheme, respectively. The results at the third and fourth points generally exhibit similar behaviors as those of the aforementioned two points (Figures 5c and 5d), and the largest difference between the two schemes also vary between 25% and 30%. Overall, compared to the 8x scheme, the ALLDIR scheme tends to show a more continuous transition, and the largest difference between the stress predicted by the two schemes ranges between 25% and 40% of the stress calculated by the old scheme.
To sketch a global map of the summed orographic stresses at the surface (GWD + FBD) and the differences between the two schemes, the two schemes are incorporated into the Institute of Atmospheric Physics-Atmospheric Global Climate Model (IAP-AGCM) version 4 . This model is the atmospheric component of the Chinese Academy of Sciences-Earth System Model (CAS-ESM). We initialize the IAP-AGCM with the real-time synoptic condition at 00 UTC 15 January 2012 using the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis ERA-Interim and run the model for 2 days. The resultant atmospheric condition (00 UTC 17 January 2015) is then used as input to the two schemes to simulate the orographic stress. It is shown that differences between the two schemes occur not only in idealized conditions but also for conditions that resemble reality ( Figure 6). The distribution of the ALLDIR shows All in all, the new scheme provides a more continuous variation of the derived topographic parameters than the simple piecewise transition of the 8x scheme. The largest difference occurring within the MDs range can be as high as 25-40% of the calculated stresses by the 8x scheme.

Evaluation of the Scheme
In this section, the performance of the scheme is evaluated in terms of its predicted surface stress and vertical momentum flux profile. Since only a limited amount of vertical momentum flux and surface stress observations are available, we use mesoscale mountain wave simulations similar to those of Kim and Doyle (2005) as reference. The ALLDIR scheme, along with the original scheme in Kim and Doyle (2005) (ORIG scheme), and the 8x scheme are compared with these reference simulations.

Mesoscale Model IAP-WRF
For the explicit simulations of mountain waves, we use the integrated modeling system IAP-Weather Research and Forecasting Model (IAP-WRF). IAP-AGCM version 4  undertakes the large-scale circulation simulations and generates the initial and lateral conditions for the system's mesoscale model WRF (version 3.2), which is used for regional high-resolution simulations. The coupling process between WRF and the IAP-AGCM is similar to what was depicted in He et al. (2013) where the WRF is initialized based on and forced by the atmospheric fields from the IAP-AGCM. The third-order Runge-Kutta time integration scheme is used for the WRF to integrate the equations in WRF. In the WRF physics, we use the WSM five-class microphysics scheme (Hong et al., 2004), and the Grell and Devenyi ensemble cumulus convection scheme (Grell & Devenyi, 2002). For turbulent processes, we use the (improved Mellor-Yamada) Nakanishi and Niino level 2.5 PBL scheme (Nakanishi & Niino, 2004). This scheme is a prognostic turbulence scheme similar to that used by Kim and Arakawa (1995) and Kim and Doyle (2005). The use of this scheme is essential since the GWD has to work together with the turbulence process to fully develop the gravity wave breaking process and thus the proper development of mountain waves with large amplitude (Kim & Doyle, 2005). For the land surface processes, we use the Monin-Obukhov surface layer scheme and the five-layer thermal diffusion scheme (Grell et al., 1994). The lateral boundary conditions of the WRF were produced by the IAP-AGCM, which was initialized everyday with the initial conditions ECMWF reanalysis ERA-Interim from December 1998 to January 1999 and December 2000 to January 2001 (Dee et al., 2011). The Rocky Mountain in Colorado, USA, which is renowned for frequent occurrences of strong downslope windstorms in its lee, is used as the simulation domain for mountain wave events. In WRF, nested-grid boundary conditions that are formulated using a one-way-interaction approach are used, with a horizontal resolution of each nested mesh equal to one third of the parent grid mesh. The innermost nested-grid domain of the model spans latitudes from 38-41.5°N and longitudes of 108-104°W (Figure 7) with the horizontal resolution of 2,000 m (174 × 174 grid points) and the vertical resolution varies with 60 levels between 50 m at the surface to~700 m at the top, which extends up to 25 hPa (~24 km). Lastly, 1-2-1 smoothing is applied to the topography to ensure model stability. Following Kim and Doyle (2005), we perform the "dry" integration (relative humidity set to 1%). To retain only the dynamical response of the mountain waves to perform a meaningful comparison, the surface heat-exchange coefficient is set to zero to suppress surface energy transfer and radiative forcing is not allowed. The integrations are hourly output for each initialized day, and only simulations after 3 hr are used for comparison.
The topographical parameters as a function of wind direction for the innermost domain is shown in Figure 8. The comparison is between the ORIG scheme, 8x scheme and the ALLDIR scheme. For the 8x scheme and the ALLDIR scheme, we use the method introduced in this paper, with the difference that for the 8x scheme, the parameter used for a specific angle is approximated by the topographic parameter of the closest angle in the eight representative directions. The ORIG scheme follows the setting of Kim and Doyle (2005) where OA = c (Nu − Nd)/(Nu + Nd), with c = 1, Nu and Nd being the number of points higher than the grid mean in the upstream and downstream area, respectively. The OL calculation follows equation (3) in that paper; λ θ is set to the same value as that of the 8x scheme (which approximates the λ θ in the inflow direction by the nearby angle in the eight representative direction).
The OA estimated from the ORIG scheme exhibits smaller values than those of the 8x scheme (Figure 8a). This is not surprising since our new method is consistent with c = 3, as indicated in the footnote of Kim and Doyle (2005), while for the ORIG scheme we only use c = 1 here. The results for the OL show strong similarities between ORIG and ALLDIR since the new method is equivalent to the ORIG method in the east and west direction. The diagonal wind direction presents small differences in the approximation of the area that cross the critical heights (Figure 8b). Lastly, the model effective grid length λ eff is set as triple the value of the root sum square of the model grid length λ θ in the along-wind and cross-wind directions as it has shown good agreement in the comparison between the parameterization scheme and the reference simulation.

Dry Mountain Wave Simulations
The simulated cases chosen to illustrate the results are shown as follows: (a) Case A simulation (08 UTC 12 January 1999). The low-level wind direction is 23.8°from the east direction at the mountain crest level. It generally goes along the direction of largest orographic asymmetry (Figure 8a) and tends to induce strong downslope windstorms (Figures 7 and 9a). A strong LLWB is evident in the vertical cross section of the wind profile, and isentropy that first occurs around 3-8 km over the steep lee-side slope of the highest peak, subsequently forms a lee wave breaking region that extends to the edge of the domain. This is clearly shown by the steepened isentrope in the lee of the mountain crest (Figure 9a). This downslope windstorm is accompanied by decelerated flow in the wave breaking region (shown by light blue colors).
(b) Case B simulation (20 UTC 13 January 2001). The low-level wind is−24.8°from the east direction at the mountain crest level. Unlike the topography of case A, the asymmetry of the mountain relief in this direction is smaller as the topography in the windward of the mountain is higher (Figure 7), and thus, presents a smaller OA (Figure 8a). Compared with the case A where there is a larger and wider wave-breaking region, the wave-breaking in case B is more centered in the lee (Figure 9b).
(c) Case C simulation (18 UTC 3 January 1999). The low-level wind direction is −67.2°from the east direction. Unlike the former two cases, the wind generally goes the along mountain ridge and does not go across the largest ridge of the Rocky Mountain that is apt for strong windstorm (Figures 7 Figure 8. Orographic parameter as a function of inflow direction for the inner simulation domain of IAP-WRF shown in Figure 7. The parameters are (a) OA, (b) OL, (c) OLP, and (d) model effective grid length λ eff (km), respectively. The green, red, and blue lines denote the parameter derived using the ORIG, the ALLDIR, and the 8x scheme, respectively.  and 8a). Consequently, the wave breaking in this case is generally weaker and more linear than in cases A and B, as potential temperature perturbation only exists at the mountain peak and none in the lee (Figure 9c).

Offline Evaluation of the New Scheme Using Mesoscale Simulations
We evaluate the new scheme following Kim and Doyle (2005) using the mesoscale simulations according to the following steps: 1. The horizontal directions of the simulated winds are domain-averaged over the depth of the surface to the height 2б h . 2. The zonal, meridional and vertical components of the simulated winds are decomposed into their area mean (overbar) and deviation (prime) from the mean. The horizontal and vertical wind deviations are high-pass filtered using the 2-D discrete Fourier transform (DFT) spatial filtering method described by Kruse and Smith (2015). The deviation fields are first transformed to spectral space using DFT, multiplied by a response function with the filter coefficients, and then inversed to produce filtered fields in physical space. A cutoff length of L = 200 km is used to remove the larger-scale variations in these fields that would obscure the gravity wave perturbations. 3. The horizontal momentum fluxes are calculated at each point, projected to the averaged low-level wind direction, then area-averaged and vertically interpolated to large-scale pressure levels of the IAP-AGCM to obtain the explicitly simulated momentum flux profile (ρu′w′ s ). 4. The area-averaged and vertically interpolated simulated variables, and the orographic statistics calculated from the topography used in the simulations are applied to the three schemes (following Broad, 1996 andDoyle, 2005) to obtain the parameterized momentum flux profile (ρu′w′ p ). 5. The parameterized momentum flux profile (ρu′w′ p ) is compared with the explicitly simulated momentum flux profile (ρu′w′ s ). 6. Below the reference level, the total parametrized surface stress (GWD plus FBD) is compared with that from the stress due to explicitly simulated surface pressure (computed following Smith et al., 2006), both of which are linearly interpolated in the vertical direction and connected to the flux profile at the reference level (Kim & Doyle, 2005;Lane et al., 2000;Webster et al., 2003).

10.1029/2019MS001921
Journal of Advances in Modeling Earth Systems 7. The parameterized flux profiles were shown to be sensitive to the choice of critical Richardson number in Kim and Doyle (2005). In the current study, the vertical stress profiles are tested under Ri c = 1 following Lott and Miller (1997) as this value is known to provide better results at the upper level, although choosing Ri c = 0.25 or Ri c = 0.5 would still yield a reasonable agreement.
Case A is first evaluated, where the low-level wind is in the 23.8°direction. Figure 10 shows the domain-averaged profiles of buoyancy frequency and projected wind. The atmospheric static stability ( Figure 10a) first decreases with height in the lower part of the domain (4 to 8 km), with a minimum value at 8 km. Above this height, stability increases to a high value of~0.022 at 13 km, and then varies between 0.017 to 0.022 between 13 and 24 km. The projected wind structure forms an approximate two-layered structure which is typically seen in strong resonant waves conditions, with a positive shear up to about 14 km and a negative shear above that (Figure 10b). The resulting Ri m diagnosed by the three schemes is thus two-layer structured, with potential wave-breaking region in the lower (~3-5 km) and upper (16-24 km) parts ( Figure 11d). The simulated vertical stress profile for Case A (Figure 11a) depicts a surface stress of about 0.47 N/m 2 , and a reference level stress of about 0.13 N/m 2 . The profile encounters the LLWB at the lower part at around 5-7 km, and then encounters an overall decreasing but wavy structure from near 7 to 16 km, and then decreases to 0 around 24 km height. The wavy structure may be due to small-scale noise in the pointwise calculation of the momentum flux, which was also noted by Kruse and Smith (2015). The strong low-level nonlinearity associated with this case is diagnosed by the ORIG and 8x scheme (e.g., 1.93 Figure 12. The parameterized stress versus the explicitly simulated domain-averaged stress at (a, d) the reference level (abscissa), (b, e) the surface, and (c, f) the weighted vertical average output for the whole 1-hr interval IAP-WRF simulation period, for which the wind direction is within ±5°of the middle of the 8x representative directions for the drag (MD). and 1.89; See Table 1), while the ALLDIR scheme depicts an even higher nonlinearity (Fr 0 = 2.62) as a result of its higher OD. This adjustment of the OD for the ALLDIR scheme helps to nearly double the value of the enhancement factor E compared with the other two schemes (Table 1). Thus, it yields a nearly double value of the reference level GWD (Figure 11a, red line) that pushes the surface stress closer to that of the reference simulation, thus alleviating potential underestimation of the surface stress (Figure 11a). In the upper level, the higher GWD in the ALLDIR profile also results in an overall lower Ri m (Figure 11d), and higher decrease in magnitude in the lower part of the profile that partly mimics the reference simulation (Figure 11a), which the other two schemes do not show.
Next, we evaluate case B, where the low-level wind is in the −24.8°direction. The wind structure and the atmospheric stability is similar to that of case A except for a low-level near-constant wind from the surface to about 7 km height, and an overall lower wind speed (Figure 10b). The low-level inverse Froude Number Fr 0 diagnosed by the ORIG scheme and the 8x scheme depicts a fairly nonlinear flow (1.07 and 1.09, Table 1), while the ALLDIR scheme depicts a strong nonlinear flow (1.61, Table 1) through the adjustment of OD. This led to a higher enhancement factor (2.72) for ALLDIR than that of the two other schemes (2.09 and 2.02), and thus a higher reference level GWD, which is closer to that of the reference simulation ( Figure 11b). The increase in reference level GWD is associated with the ALLDIR surface stress being Figure 13. Linear correlation (R), goodness of fit (R 2 , defined as the quotient of the regression sum of squares divided by the total sum squares), and the F ratio (defined as the quotient of the mean square due to regression divided by the mean square error) of the regression of the drag profile predicted by the ORIG scheme (blue), new scheme with 8x direction (green), and the ALLDIR scheme (red) to the reference simulation for (a) reference level stress, (b) surface stress, and (c) weighted vertical average stress and (d, e, f) the corresponding variables where only those wind directions within ±5°of the MD. For display purpose, the value of F ratio in the panels is divided by 4,000.
closer to that of the reference simulation, along with a lower Ri m in the near~5-6 km that caused decreased stress in the ALLDIR vertical stress profile similar to that of the reference simulation.
For the case C, the low-level wind is in the −67.2°direction and mostly crosses the higher plains of the domain (Figure 7). This case is shown in the previous section to be linear-unlike cases A and B-which is also consistent with the near-constant wind from 7 to 12 km height (Figure 10b). The diagnosis of the ORIG and 8x schemes indicates an overall nonlinear event (with Fr 0 as high as 1.7 and 1.75 for the ORIG and 8x scheme, respectively). The ALLDIR scheme, however, indicates a near linear event (Fr 0 = 1.12), and thus, its reference level GWD is adjusted to a closer value to that of the reference simulation ( Figure 11c) compared to the overestimated reference level GWD in the ORIG and 8x schemes. This causes a perfect match of the surface stress and an upper level stress profile closer to the reference profile than that of the other two schemes.
Overall, the ALLDIR scheme shown in these three cases tends to perform better than the other two schemes when compared to the reference simulation. In addition to this case study comparison, we also compare the predicted reference level stress, the surface stress, and the weighted vertical average stress of the three schemes against the hourly simulated stress from the IAP-WRF for all flow directions (Figures 12a-12c), and for direction within the eight ranges of the MDs (Figures 12d and 12e). In all, there are 2,018 and 480 points for all directions and within the range of the MDs used for the scatter plot, respectively. While the three schemes all show a reasonable shape about the reference simulation, the ALLDIR tends to show a narrower scatter than that of the other two schemes. This better performance is also demonstrated by the evaluation metrics of the scatter using linear correlation, the goodness of fit (defined as the quotient of the regression sum of squares divided by the total sum of squares), and the F ratio (defined as the quotient of the mean square due to the regression divided by the mean square error) (Wilks, 2006) shown in Figures 13a-13e. For practically every comparison, the ALLDIR scheme tends to perform better. This indicates the potential benefit in further considering the 3-D anisotropy of the subgrid terrain in the GWD and FBD schemes. Note that while the reference GWD and the weighted vertical average stress of the three schemes tend to show good correspondence with the reference simulation, the surface pressure stress tends to show higher uncertainty and a potential underestimation. Since the reference level GWD predicted by the three schemes exhibits a closer agreement to the reference simulation, this indicates potential underestimation of the FBD for the three schemes. It might be tempting to retune the C d to pull the surface stress higher. However, this may lead to the overestimation of the surface stress for the few points that are equal to or higher than the reference simulation's stress. Therefore, the FBD is currently kept as it is, but the improvement of the FBD estimate is a subject which merits further study.

Flow Regime Diagram
Another measure of the overall behavior of the drag parameterization schemes is the "regime diagram" mentioned by Kim and Doyle (2005). This serves to check the relationship among the various sources of drag caused by the subgrid-scale orography (such as GWD due to hydrostatic waves, LLWB, trapped lee waves, internal wave reflections, and FBD due to low-level flow blocking) considered with the drag parametrizations. In this section, we follow Kim and Doyle (2005) and plot the low-level stress normalized by its linear counterpart (kρ 0 N 0 U 0 б 2 h with k = 8 × 10 −6 , following McFarlane, 1987) as a function of the inverse Froude number divided by OD. The ALLDIR scheme is then compared with the values obtained by the ORIG scheme (note that the 8x scheme results are not shown since they yield similar results to that of the ORIG scheme). As shown in Figure 14, while ORIG yields a two-group scatter that corresponds to the high and low values of OD, the ALLDIR scheme tends to show a large scattering, with values all well within the bound of the ORIG scheme. This is not surprising since the application of the parameters to all directions would Figure 14. Flow diagram of the ratio of parametrized GWD at the reference level and/or FBD to the linear drag as a function of the domain-averaged low-level inverse Froude number divided by OD. Total of 2,018 "snapshots" have been taken at various times with 1-hr intervals from the mesoscale simulations performed with IAP-WRF. The dots denote ORIG GWD (gray), FBD (blue), GWD + FBD (black), ALLDIR GWD (green), FBD (yellow), and GWD + FBD (red). allow a larger range of OD than that of the four values in the ORIG scheme. Another trait to note is that despite the higher scattering, the ALLDIR scheme shows a dominance of GWD in low Fr 0 /OD state, and comparable contribution of the GWD and FBD under high Fr 0 /OD state like that of the ORIG scheme.

Conclusions and Discussion
The present study introduced a new orographic drag parameterization that considers orographic anisotropy for all flow directions. This was done by refining the direction-related parameters as follows: An OAV was first defined for a coarse grid as the normalized vector between the grid center point and the center of mass in the grid, from which the OA in a given flow direction was derived as the projection of the OAV in this direction; the OL for a direction was calculated as the percentage of the areal coverage that exceeds the mean grid height in the half area along the given flow directions; and the calculation of the model grid length λ was extended to all directions. Sensitivity tests of the orographic drag were conducted under the new scheme using a 5 m/s vertically uniform wind along different directions at four coarse grid points in typical mountain regions. The new (ALLDIR) scheme is shown to provide a more continuous transition of the derived topographic parameters than the simple piecewise transition of the 8x scheme, where the largest difference that occurs within the MDs range can be as high as 25-40% of the calculated stress by the scheme with only eight directions. Global simulation initialized with realistic initial conditions showed that these differences exist in some of the major mountain regions. The predicted stress and momentum flux profile by the various schemes were compared with the realistic mountain-wave simulations run under the integrated modeling system IAP-WRF over the Rocky Mountain for 4 months of daily simulations initialized with realistic initial conditions. It was shown that the new scheme predicted an overall narrower scatter about the reference simulation for the predicted reference level stress, surface stress and weighted vertical average stress than the old scheme, with a higher value of linear correlation, goodness of fit, and the F ratio, thus indicating its potential benefit over the other two schemes.
Aside from the improvements brought by the new scheme, issues remain to be investigated in future studies. First, the new scheme follows the original scheme in that it only treats drag divergence and does not provide means for treating drag convergence (Kim & Arakawa, 1995). In addition, the scheme's performance in the large-scale model needs further evaluation.