Application of Digital Filtering Initialization in A Limited Area Model over Indian Region

In order to efficiently and self-consistently describe the transport ofhigh-latitude ionospheric plasma into the magnetosphere,we have developeda Dynamic coupled Fluid-Kinetic(DyFK)model for describing plasmaflow along a magnetic flux tube.The collision-dominated ionospheric plasmais treated with a low-speed fluid approach for altitudes between 120-1100km,while a generalized semikinetic approach is used for the topside andhigher altitudes,starting at 800km.This paper presents a description ofthe new DyFK model,along with illustrative results obtained from modelingthe effects of soft(100eV)auroral electron precipitation using the newmodel.The results demonstrate the accuracy and efficiency of the newDyFK model and illustrate how soft electron precipitation provides a mechanismto drive observed upflows of high-latitude ionospheric plasma alonggeomagnetic field lines.


INTRODUCTION
It is well known that if the mass and wind fields in the initial state are not dynamically balanced, spurious high-frequency large-amplitude inertial gravity-wave oscillations occur in numerical weather prediction based on primitive equations. It is therefore necessary to remove or reduce the high-frequency components from the initial data so as to bring the mass and wind fields into a proper state of balance. Although the high-frequency noise can be controlled to some extent by choosing suitable dissipation schemes and time filters, the best way to elimi nate them is to enforce a balance between the mass and wind fields in the model's initial state through a process of adjustment commonly known as "initialization". Many initialization schemes have been proposed viz., static initialization (Phillips, 1960), dynamic initialization (Miyakoda and Moyer, 1968), linear and non-linear normal-mode initialization (Dickenson andWilliamson, 1972. Machenhauer, 1977;Baer and Tribbia, 1977) and dynamic normal mode initialization (Sugi, 1986). More recently, Lynch and Huang (1992) have put forward the concept of digital filtering initialization (DFI) technique based on digital filters which are of two types -non-recursive and recursive. In non-recursive filters, input is processed to pro duce output, whereas in recursive filters, the feedback of the output is involved. Recursive filters are more powerful than non-recursive ones but are more problematic as feedback of output can give rise to instability. The DFI technique can be used adiabatically as well as diabatically. The filter operates on a time series of dependent variables which are produced by adiabatic backward and adiabatic or diabatic forward time integration of the forecast model. The digital filter removes the high-frequency noise from the time series. A balanced initial state is thus achieved through an adjustment of all dependent variables of the numerical model. In this technique, there is no need for simplifying assumptions such as constant Coriolis parameter and linearization about a basic state which are generally made in normal-mode initialization (NM!). Lynch and Huang (1992) have applied a non-recursive adiabatic digital filtering initialization (ADFI) procedure to initialize the input data for a high-resolution lim ited area model (HIRLAM). Huang and Lynch (1993) have extended ADF1 to include diabatic processes and formulated a non-recursive diabatic digital filtering initialization technique (DDFI). Their results from applications of DFI to HIRLAM have shown that both ADFI and DDFI can yield results which are at least as good as those of NMI. Both ADFI and DDFI produce more realistic and synoptically better-organized initial fields than NMI. Moreover, the primary advantage of DFI over NMI is its simplicity, both in physical understanding and practical application. Another significant and practical advantage of DF1 is that the initializa tion procedure involves the model integration which ensures compatibility between initializa tion and forecast. Any changes to the model formulation automatically apply to the filtering process which may easily be incorporated in the model code. In this paper, the authors have made an attempt to implement DFI in a limited area model. Both the adiabatic (ADFI) and diabatic (DDFI) versions are used separately in the model and their impact on the forecast performance of the model for short range prediction over Indian monsoon region is examined and discussed with reference to the uninitialized forecast.

THEMODEL
The fo recast model used in the present study is a limited-area primitive equation grid point model (Singh et al. 1990) in terrain following a sigma (er) coordinate system. The gov erning equations are : the zonal and meridional momentum equations, thermodynamic equa tion, moisture continuity equation, surface-pressure tendency equation, hydrostatic equation and gas law. The equations are : (3) (8) u and v are the zonal and meridional component of horizontal wind vector respectively. fJ is the potential temperature , q is the specific humidity, ¢ is the geopotential, T the temperature, f the coriolis parameter, m the map scale factor on Mercator projection, u* = u/m, v* =vim, P = (p/lOOQ)Rfcp, tr = Psp T , p T is the pressure at the top of model atmosphere, Ps is the surface-pressure, p is the pr � ssure at any level, R the gas constant, CP is the specific heat of air at constant pressure. F 1s the eddy stress of momentum, D u, D v, D and D9 are the horizontal diffusion ofu,v, fJ and q respectively.Hand E are the vertical eddy fluxes of heat and moisture respectively. Q is the diabatic heating rate per unit mass and M is the change of moisture per unit mass. The system of equations (1) to (6) form a complete set for the prog nostic Variables U, V, fJ, q and tr, with</> and Gare determined diagnostically from (6) and (7).
All the variables except for vertical velocity a , are defined at the middle oflayers. a is specified at the interfaces of the model layers. a is set zero both at the top (a = 0) and bottom ( a=l) of the model atmosphere.
The original version (Singh et al., 1990) of the model had six layers in vertical and a was p-pT defined as a = with p T as 100 hPa. In the present study, the number of vertical p, -pT layers is increased to ten and the top of the atmosphere is set at zero pressure level (p T = 0) so that the definitions of a and n are changed to The other features of the present version remain the same as in original version of the model. The horizontal distribution of prognostic variables follows Arakawa's B-grid on the sigma (er) coordinate (Arakawa, 1972). The forecast area spans from 10°S to 40°N and from 40°E to l 10°E and contains 51 x 39 grid points with a horizontal resolution of 150 km on Mercator projection. Following Arakawa and Mintz (1974) and Arakawa and Lamb (1977), a fourth order finite difference scheme for horizontal advection terms in momentum equations and a second order scheme for other terms in the governing equations are used.
The physical processes included in the model are cumulus convection, large-scale con densation, dry-convective adjustment, shallow-convection, horizontal and vertical diffusion, surface fluxes of momentum, heat and moisture, and short-wave and long-wave radia tion. Cumulus convection is parameterized following the scheme of Kuo (1974). The param eterization of large-scale condensation is based on grid-point saturation. Whenever, grid point super-saturation occurs, excess moisture is condensed as precipitation, latent heat is released and the temperature and specific humidity are adjusted to saturation using the scheme of Manabe et al. (1965). The evaporation of falling rains in the unsaturated layers below is not allowed for. Whenever temperature lapse-rate exceeds the adiabatic lapse-rate, the atmosphere be comes hydrodynamically unstable to the small-scale convective motions which rapidly trans port momentum, heat and moisture vertically. These transports are represented in the model by a dry-convective adjustment scheme (Kanamitsu, 1975). The parameterization of horizon tal and vertical diffusion and surface fluxes computation are as in Singh et al. (1990). The horizontal diffusion is based on fourth-order Laplacian with constant diffusion coefficients for all the grid points in the domain. The surface-fluxes of momentum, heat and moisture are evaluated through Bulk-aerodynamic formulas in which the computations of Bulk coefficients are based on Kondo (1975). The surface-fluxes are considered over sea only. The long-and short-wave radiation-flux computations are based on the band-model of Harshvardhan and Corsetti (1984) and Lacis and Hansen (1974) respectively. Other details of parameterization are available in Singh (1985) and Singh et al. (1990).
The model is integrated by an explicit time integration scheme (Tatsumi, 1983;Bandyopadhyay and Singh, 1994). During forecasting, boundary values are replaced by ob served values. A five-point relaxation zone is used to ensure smooth transition between the boundary and interior values. To do this, the model predicted value of any variable Xis modi fied in the five point relaxation zone by using the following scheme. Xm = aXp + q-a)X R where a=O, .1, .4 .7 and .9 at the boundary, one, two, three and four points inside the boundary. For all other points a=l. X P stands for the model predicted value of X, X R is the real (observed) value, and Xm is the modified value.

Digital Filter
A filter is any device designed to carry out the separation of waves based on frequencies.
The process of filtering involves selection of those components of an assemblage of waves having some particular properties and removal of wave-components which lack it. The selec tion is generally based on frequency of the component waves.
Filtering theory originated from the need to design electronic circuits with precise fr e quency selection characteristics like in radio, telecommunications etc. These are analogue filters and act on continuous time-signals. Digital-filters developed in the context of discrete time-signal processing are recently being used in the analysis and solution of problems of geophysics and NWP. In the field ofNWP, digital-filters have been applied to the problems of initialization and objective analysis and they have considerable potential for other applications in NWP. In NWP, input consists of low-frequency waves contaminated by high-frequency waves (Noise). A filter is, therefore, used to carry-out the selection oflow and high-frequency waves and to filter-out the noise.
It is mentioned in Section 1 that digital-filters are of two types -recursive and non-recur sive. For a discrete function f of time t, a non-recursive digita� filter may be defined math ematically as N Yn = L a k f,, -k (9) k=-N and a recursive digital filter as (10) From the above two definitions, it is clear that for a non-recursive filter, out-put depends both on past and future values off but not on the other out-put values. For a recursive filter, out-put depends on past and present values of input and also on the previous out-put values. Since in the present study our interest is on the application of non-recursive filter, the design of non-recursive digital filter is summarised here briefly.
Let us consider a function f(t) of time t consisting of low and high-frequency components. We further assume that f(t) is known only at some discrete moments t = nAt at a time- where h is the discrete notation of h(t) = h(t ) = h(n At)= h n n . n.
From the well-known procedure for finding the Fourier coefficients, hn can be found as from the combination of two relations ( 11) and ( 13 ), the expression for hn reduces to hn = sin ( n e )l(n 1t') If fn * denotes the low-f requency parts of (f " ) from which all the high-frequency components with frequencies greater than e c are eliminated, then Now the convolution theory of Fourier series states that H( 8)*F( 8) is the transform of convo lution of h with f , that is, where, F stands for the Fourier transform operator. So, the above convolution theory implies that f 11 * = 0 * J 1 = th kj n-k k�o: once we find out the weights hn and compute the sum in the right-hand side of (18), the result gives the filtered out-put f" * from which the high-frequency components are removed. In prac tice, however, a finite number of coefficients h can be used, and so the summation in (18) n must be truncated at some finite value of k say, N. We finally, therefore, compute In (19), the sequence (hn) is called the filtered response and fn * is the filtered output we are interested in. However, the use of (19) has associated problems since the truncated Fourier series gives rise to spurious oscillations . To reduce the oscillations and to improve the effi ciency of the filter, the coefficients given by (15) are multiplied by a window-function. Lynch and Huang (1992), has used a Lanczos window defined as The same window is used in the present study. So the co-efficients h0 are finally computed as More details on digital-filters are available in any book on the subject.

Data
The source of data is the global atmospheric research program (GARP). The First GARP Global experiment level Illb data of 12 GMT, 7 July, 1979 are used in this study. The prevail ing synoptic situation governing this data over the domain of interest was a monsoon depres sion. The globally analysed data of 1.875° resolution and at 15 standard pressure levels are interpolated to the model grid-points and ten sigma levels.

Initialization
The digital filter defined by (19) is used to initialize the impact date. Two sets of initial ized input are prepared for the forecast model. One set is obtained by applying ADFI tech nique and the other set is produced by D D FI scheme. The general procedure of DFI is to apply the filter to a time series of dependent variables produced by the forecast model. In ADFI two adiabatic short range integrations, one forward (from time t=O to t =N /:1 t hrs.) and the other backward (from t=O to t=-N /:1 t hrs.), are carried out in time. As the integration proceeds, the model variable is multiplied by a filter coefficient and accumulated to produce a balanced initial state. This is done for all the prognostic variables u, v, 8, q and p s' at each grid point and at each model level.
Formally, the ADFI procedure can be written as where XADF is the initialized field, h, is the filter-coefficient, i is the time-step, N is the number of integration steps for half of the filter-span and X\ is the time-series generated by adiabatic integration of the model. In DDFI, first the model is integrated backward adiabatically from time t=O to t=-N /:1 t and then a forward diabatic integration is carried out for the full filter span from t =-N /:1 t to t=+N /:1 t hours. The backward integration is done to ensure that the filter output is centered around the initial time. Since the diabatic processes are irreversible, backward integration is done adiabatically. During forward integration, each variable is multiplied by a filter coeffi cient at the end of each time step and summed up to ach ieve a balanced initial state for smooth time evolution of the forecast.
The DDFI method can be expressed as The equation (23) has the same form as (22) except that the time-series x n i is yielded by the diabatic integration of the model with all the diabatic processes included. Thus �DF is the initialized field. Lynch (1990) in his experiment with a barotropic model noticed that a cut off period of 3 hrs produce similar results to those with 6 hrs cut-off period. But a cut off period of 12 hrs showed damping effect on the meteorological mode. This is our reason for choosing the cut off period of 6 hrs in the present study.
In both ADFI and DDFI, the total filter span and cut-off period are chosen to be 6 hrs (Lynch and Huang, 1992). For time integration,an explicit leap-frog scheme with a time step ( /:1 t) of 3 mins is used. The number of time steps (N) required to cover half the filter span (3 hr) is determined from the relation 2N /:1 t = 6 hr. The filter coefficients are computed follow ing Lynch and Huang (1992). The coefficients obtained from (21) are functions of time only and remain the same for both ADFI and DDFI methods. The boundary values of the variables are held constant during initialization.

RESULTS AND DISCUSSION
The results of the impact of digital filtering initialization on the forecast performance of the model are discussed in this section. Three 24-hr forecasts are produced starting from ADFI, DDFI and uninitialized analysis (hereafter referred to as NOi). The performance of ADFI and DDFI on the model forecasts are evaluated by comparing the ADFI and DDFI forecasts with those obtained from NOi input. The two initialized forecasts are also compared with each other. The individual performance of ADFI and DDFI are assessed in the light of the three essential requirements for any initialization scheme. These are, i) the changes made in the initial data due to initialization should be reasonably small, ii) the forecast produced from initialized fields should not be degraded as compared to the forecast obtained from uninitialized data (NOi), and iii) spurious high-frequency noise is eliminated or reduced from the forecast.

Changes Made in the Initial Fields by Initialization
To examine the changes made in the original analysis due to initialization, the root mean square difference (RMS) and maximum difference (MAXD) in the wind components (u, v), temperature (T), and surface pressure (ps) between the initialized and uninitialized fields are computed and presented in Table 1. It can be seen that the changes made in the original analy sis by both ADFI and DDFI are reasonably small. The maximum changes as well as RMS changes made by DDFI are slightly smaller than the changes made by ADFI.
The mean sea level (MSL) pressure maps for all the three sets of input ADFI, DDFI and NOi are shown in Figurel, which shows that all the three maps closely resemble each other; hardly any difference between the three maps is discemable.
The uninitialized (NOi) and the two initialized (ADFI and DDFI) wind fields are shown in Figure 2. The centre of the cyclonic circulation, flow pattern, and the strength of the wind are basically similar in all three initial fields. This confirms that the changes made in the original wind field by both the initialization schemes are very small.

Changes in the Forecast
To examine the changes made in the forecasts by ADFI and DDFI as compared to NOi, the RMS and MAXD values are computed for u, v, T and p,, between the 24-hr forecast fields obtained from two initialized input and the uninitialized analysis. The statistics presented in Table 2 show that, the RMS and MAXD values for all four parameters between ADFI and NOi forecast and also between DDFI and NOi forecasts are reasonably small. The changes made in the forecast fields by DDFI are smaller than the changes made by ADFI. As compared to the values of Table 1, the RMS and MAXD values in the forecasts for u, v and T for both ADFI and DDFI are somewhat larger than the respective changes induced by ADFI and DDFI in the original analysis. In the case of surface pressure however, the forecast changes induced by both the initialization schemes are substantially reduced compared to the corresponding changes in the original analysis. In particular, the forecast RMS values are reduced by nearly 50%. The surface-pressure tendency is related to the vertically integrated divergence field. In other words, surface-pressure is sensitive to noise in a vertically integrated sense. Thus the better forecast changes regarding surface-pressure is primarily due to the reduction of noise in a vertically integrated sense.
To further assess the changes made by ADFI and DDFI on the model forecast as com pared to NO! forecast, spatial distribution of the 24-hr forecast MSL pressure obtained from MSL Pressure (hPa) T=OO AOFI (lsotachs in rns-1).  the two initialized data sets of ADFI and DDFI and from NOi data are examined and presented in Figure 3 along with the verifying analysis (12 GMT 8 July 1979). The forecast RMS errors (RMSE) for u,v, T and p,, between the observed analysis and each of the three forecasts from NOi, ADFI and DDFI data are presented in Table 3. It can be seen that for all variables forecast RMSE starting from the two initialized data ADFI and DDFI compare well with each other at almost all levels. Moreover, the initialized forecast errors for u and v are, in general, slightly smaller than the uninitialized forecast errors in the lower troposphere. In the upper troposphere, however, initialized forecast errors are margin ally higher than or comparable to uninitialized errors. The RMSE for temperature (T) in all the forecasts are seen to be comparable with each other at all levels. For surface pressure ps, however, the errors in NOi forecast are marginally smaller than the errors in the initialized runs.

Surface pressure tendency
The basic purpose of any initialization scheme is to eliminate or reduce the high-fre quency noise from the forecast fields. A common practice for measuring the noise level in the forecast is to examine the oscillations of the mean absolute surface pressure tendency (ASPT) as the time integration proceeds. The parameter ASPT is defined as ASPT =ti\ a�; t

485
where M is the total number of horizontal grid points. The time variation of ASPT in the first 6 hrs of model forecast obtained from each of the three sets of inputs ADFI, DDFI and NOi are depicted in Figure 5 which clearly shows that ASPT oscillates noisily during uninitialized forecast and this indicates the presence of imbalance between the mass and wind fields in the uninitialized input. The initial value of this quantity (ASPT) in the NOi curve is of the order of 10 hPa (3 hr)-1• The ADFI and DDFI curves show that the high-frequency oscillations in NOi forecast is effectively controlled when the model is integrated with either ADFI or DDFI data as input. The initial values of ASPT in ADFI and DDFI data are respectively of the order of 2 hPa (3 hr)-1 and 3.6 hPa (3 hr)-1• Figure 5 also reveals that the large-amplitude high frequency oscillations in NOi forecast decreases as the forecast progresses. This is perhaps due to the inherent damping processes in the model. From the two initialized forecast curves, it is seen that at the beginning of the forecast, ADFI curve is less noisy than the DDFI curve. It is also found that after about 12 hrs of forecast, all the three curves almost coincide (figure not shown) and thereafter display very close fairly steady runs.

Surface pressure
The time evolution of surface pressure for the first 6 hrs of forecast at a model grid point (27,16) is depicted in Figure 6. This point is chosen because it is close to the closed low pressure system depicted in forecast MSL pressure chart (Figure 3). Severe contamination of forecast by high frequency noise in the uninitialized run is clear in the figure. This once again indicates the imbalance between the mass and wind fields in the uninitialized data. Although the residue of high-frequency noise is still there in the two initialized forecasts, the noise level is considerably less in both initialized runs. The difference in the amplitude of oscillations shown in ADFI and DDFI curves are very small (less than 0.5 hPa).    The reduction of noise is to a great extent, clearly seen in the two initialized forecasts. This indicates the efficiency of digital filtering in reducing or controlling the high frequency noise in the forecast. In the first few hours DDFI forecast is less noisy than ADFI run.

Precipitation
The area-averaged forecast rainfall rates (mm/day) starting from NOi, ADFI and DDFI input are shown in Figure 8. Rainfall rates in each forecast increases with forecast time. No significant differences are evident in the three forecast curves. Figure 9 shows the spatial distribution of 24-hr accumulated rainfall amount for each of the three forecasts. The figure shows that the forecast rainfall patterns from both initialized runs are in good agreement with that of the uninitialized one. There is virtually no difference in the rainfall patterns. The areas as well as the rainfall amounts of maximum rainfall are nearly identical in all the runs.

5.SUMMARY
Non-recursive digital filtering initialization technique is used to initialize the input data for a limited area model for short range prediction over Indian region. Two sets of initialized  input are prepared. One set is produced by ADFI process of Lynch and Huang (1992) and the other set is through DDFI method of Huang and Lynch (1993). The model variables at the boundaries are held constant during initialization. It would be more realistic if the variables were allowed to vary in time. However, since the time-dependent boundary conditions nor mally vary slowly, and as the integration for the initialization purpose is done for a very short period, it is unlikely that the use of temporal variation at the boundaries would make any significant change on the initialized fields and thereby improve the forecast. The performance of DFI in the forecasting is assessed by comparing it with the forecast produced from uninitialized input (NO!). There is marginal improvement in the DFI forecast over the NO! forecast because the model has built-in damping processes and smoothers. These smoothers/filters suppress the growth of high-frequency noise in the uninitialized forecast and at the end of 24-hr integration a smooth forecast field appears even without initialization, and as such the NO! forecast field is found quite comparable with those of DFI forecast. However, it should be noted that DFI shows considerable impact at the beginning of the forecast. This impact can be seen from Figs. 5, 6 and 7 which suggest clearly that DFI could effectively control the high-frequency oscillations in the forecast at least in the first six hours.
Further the individual performance of ADFI and DDFI are evaluated in the light of three essential characteristics (mentioned in section 4) any initialization technique should have. Both the ADFI and DDFI are found to produce a fairly well balanced model initial state. Both the maximum (MAXD) and RMS changes induced in the original analysis by DDFI are somewhat smaller. The changes made by DDFI in the 24-hr forecast fields are also smaller than the changes made by ADFI.
It is expected that an initialization scheme including physical processes (DDFI) should produce more ·realistic and better organized initial fields and consequently better forecasts than its adiabatic version (ADFI). The reason that the performance of DDFI is not as good as expected remains, probably, in the inherent disadvantage of DDFI. In non-cursive DDFI a backward integration is needed in order to ensure that the filter output is valid around the initial time, and this implies that the subsequent forward integration will not pass through the original analysis at t=O. In an ideal case, what is needed is a balanced initial state deduced from a forward diabatic integration starting from the initial data. This could be achieved by using a recursive filter which allows for feedback of previous output. Nevertheless, the impact of DFI on the model has been found promising and encouraging to implement DFI with a recursive filter.