A Method of Up grad ing a Hy dro static Model to a Nonhydrostatic Model

As the sigma-p co or di nate un der hy dro static ap prox i ma tion can be in ter preted as the mass co or di nate with out the hy dro static ap prox i ma tion, we pro pose a method that up grades a hy dro static model to a nonhydrostatic model with rel a tively less ef fort. The method adds to the prim i tive equa tions the ex tra terms omit ted by the hy dro static ap prox i ma tion and two prog nos tic equa tions for ver ti cal speed w and nonhydrostatic part pres sure p’. With prop erly for mu lated gov ern ing equa tions, at each time step, the dy namic part of the model is first in te grated as that for the orig i nal hy dro static model and then nonhydrostatic con tri bu tions are added as cor rec tions to the hy dro static so lu tions. In ap ply ing phys i cal parameterizations af ter the dy namic part in te gra tion, all phys ics pack ages of the orig i nal hy dro static model can be di rectly used in the nonhydrostatic model, since the up graded nonhydrostatic model shares the same ver ti cal co or di nates with the orig i nal hy dro static model. In this way, the ma jor ity codes of the nonhydrostatic model come from the orig i nal hy dro static model. The ex tra codes are only needed for the cal cu la tion ad di tional to the prim i tive equa tions. In or der to han dle sound waves, we use smaller time steps in the nonhydrostatic part dy namic time in te gra tion with a split-ex plicit scheme for hor i zon tal mo men tum and tem per a ture and a semi-im plicit scheme for w and p’. Sim u la tions of 2-di men sional moun tain waves and den sity flows as so ci ated with a cold bub ble have been used to test the method. The ide al ized case tests dem on strate that the pro posed method re al is ti cally sim u lates the nonhydrostatic ef fects on dif fer ent at mo spheric cir cu la tions that are re vealed in the o ret i cal so lu tions and simulations from other nonhydrostatic models. This method can be used in upgrading any global or mesoscale models from a hydrostatic to nonhydrostatic model.


IN TRO DUC TION
The prim i tive equa tions used as gov ern ing equa tions in at mo spheric hy dro static mod els in clude hy dro static approximation that re places the ver ti cal mo men tum equa tion by the hy dro static bal ance equa tion.The ap prox i ma tion is very ac cu rate when ver ti cal ve loc ity, and hence ver ti cal acce leration, is very small.The hy dro static ap prox i ma tion not only sim pli fies the gov ern ing equa tions but also elim i nates sound waves from the so lu tion to avoid the ne ces sity in handling those fast mov ing waves in time in te gra tion.However, if one needs to sim u late at mo spheric cir cu la tion with strong ver ti cal mo tions, such as thun der storms or hur ricanes, a nonhydrostatic model may be de sired.In this pa per, we de scribe a method that ex tends a hy dro static model to a nonhydrostatic model by add ing the terms and equa tions neglected by the hy dro static ap prox i ma tion to the prim i tive equa tions as cor rec tions to the dy namic part so lu tion at each time step.With this ap proach, the model can be run as a hydro static or nonhydrostatic model de pend ing upon the user's choice for ap pli ca tions.The only ap prox i ma tion made in this method is to linearize the pres sure ten dency equa tion to for mu late a sta ble scheme for the non hy dro static part time in te gra tion.Since sound wave pre dic tion may be dis torted by linearization, nonhydrostatic mod els de vel oped by this method are not suit able to sim u late sound wave dom i nated phe nom ena.How ever, they should be per fectly ad e quate for nu mer i cal weather pre dic tion (NWP) since sound waves play an in sig nif i cant role on weather changes.
For a given tem per a ture pro file, the hy dro static approximation makes height and pres sure monotonically de pend on each other so that pres sure can be used as a ver ti cal co ordi nate.To better han dle bot tom bound ary con di tions, the terrain-fol low ing sigma-p co or di nate, de fined as s = (p -p Top ) / (p Sur face -p Top ), is usu ally used as the ver ti cal co or di nate in hy dro static mod els.Laprise (1992) in tro duced a new concept that if we in ter pret the hy dro stati cally bal anced part of pres sure as the air mass weight above the area, the sigma-p co or di nate can still be used in nonhydrostatic mod els and it is called the mass co or di nate.When in ter pret ing the hy drostatic pres sure as air mass weight per area, we can for mu late the gov ern ing equa tions for a fully com press ible at mo sphere with out the hy dro static ap prox i ma tion in the mass co or dinates.Con se quently, the gov ern ing equa tions, and hence the whole codes, of a hy dro static model can be used as the ma jor part of a nonhydrostatic model ex tended from the hydrostatic model.In the dy namic in te gra tion of the nonhydrostatic model, af ter the hy dro static ten dency has been cal cu lated the nonhydrostatic ten dency from terms and equa tions omit ted by the hy dro static ap prox i ma tion can then be added to com plete the time step.This cor rec tion type time in te gra tion for the dy namic part works as the time split ting in te gra tion, with which dif fer ent terms are in te grated by differ ent nu mer i cal meth ods.In the phys ics part time in te gration af ter the dy namic in te gra tion at each time step, all physical parameterization pack ages of the orig i nal hy drostatic model can be di rectly ap plied for the nonhydrostatic model as they share the same ver ti cal co or di nate.The challenge in this ap proach is to or ga nize the gov ern ing equa tions into a form that hy dro static and nonhydrostatic terms are sep a ra ble, and to se lect proper prog nos tic vari ables and numerical schemes for which the time in te gra tion is sta ble and ef fi cient.The Weather Re search Fore cast ing (WRF) model also uses mass co or di nates in its two dy namic cores: NMM de vel oped at the Na tional Cen ters for En vi ron men tal Pre dic tion (NCEP, Janjic et al. 2001) and ARW de vel oped at the Na tional Cen ter for At mo spheric Re search (NCAR, Klemp et al. 2007).How ever, we choose dif fer ent pro gnostic vari ables and dif fer ent nu mer i cal meth ods in handling sound waves.These dif fer ences will be dis cussed in sec tion 4.
There are other ver ti cal co or di nates used in non hy drostatic mod els, such as the sigma-z co or di nate in the Cou pled Ocean/At mo sphere Mesoscale Pre dic tion Sys tem (COAMPS, Hodur 1997) and the sigma-p* co or di nate in the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5, Dudhia 1993) and the Purdue Mesoscale Model (PMM, Hsu and Sun 2001).The sigma-p* co or di nate de fines the ver ti cal co or di nate by a ref er ence sur face pres sure p* that is con stant in time and var ies in height only.There fore, the sigma-p* coor di nate is equiv a lent to the sigma-z co or di nate.The sigma-z co or di nate has an ad van tage in that the co or di nate is not changed with time.How ever, the sigma-z co or di nate involves with a rigid up per lid that can ar ti fi cially in crease pres sure or re duce mass (Klemp et al. 2007).In the MM5, the non-con serv ing prob lem is min i mized by ig nor ing diabatic terms in the pres sure pre dic tion equa tion (Dudhia 1993).The mass co or di nate is fa vored over the sigma-z coor di nate in the WRF-ARW dy namic core be cause of mass con ser va tion and con ve nience in switch ing to hy dro static in te gra tion (Klemp et al. 2007).

GOV ERN ING EQUA TIONS
For a hy dro static model, the prim i tive equa tions in the sigma-p co or di nates can be writ ten as: (1) where is hor i zon tal wind ve loc ity, f is the Coriolis pa ram e ter, f = 2Wsinj, f is the geopotential of a sigma-p level, f = gz s - F r is mo men tum fric tion, and S is mois ture source.For sim plic ity we have cho sen p Top = 0 in de fin ing the sigma-p co or di nate.Equation (3) is a form of the con ti nu ity equa tion de rived from the trans formed con ti nu ity equa tion for any gen eral co or dinate h (Kasahara 1974): The hy dro static bal ance Eq. ( 5) is the ver ti cal mo men tum equa tion with the hy dro static ap prox i ma tion, i.e., dw/dt » 0. The pres sure ten dency equa tion be comes a triv ial equation with the hy dro static ap prox i ma tion (see Ap pen dix) and the pres sure can be di ag nosed from the sur face pressure p s as p = p s s.
To add nonhydrostatic com po nents to the prim i tive equa tions, we se lect ver ti cal ve loc ity (w) and the nonhydrostatic part of pres sure (p¢), de fined as the dif fer ence be tween to tal pres sure and hy dro static pres sure, as two ex tra prog nos tic vari ables.Fol low ing Laprise (1992), we can write the gov ern ing equa tions with out the hy dro static approx i ma tion in the mass co or di nates as: The pres sure ten dency Eq. ( 10) in volves the very del i cate bal ance be tween heat ing and di ver gence.We can re write this equa tion in terms of hy dro static and nonhydrostatic pres sure as: (12) Since p is much greater than p¢ for cir cu la tion re lated to weather pre dic tion and dp¢/dt is im por tant mainly for sound waves, we sim plify this equa tion by first as sum ing the nonhydrostatic part of pres sure plays no sig nif i cant role on the hy dro static part of pres sure (or mass) ten dency cal cu lation.As dem on strated in the ap pen dix, in the mass co or dinates, if only the hy dro static part of pres sure is used in cal -cu lat ing D 3 by Eq. (A7), a triv ial equa tion for the hy drostatic pres sure ten dency can be ob tained by add ing the conti nu ity and ther mo dy namic equa tions to gether as: (13) By sub tract ing Eq. ( 13) from Eq. ( 12), a nonhydrostatic pres sure ten dency equa tion can be ob tained as: (14) Since Eq. ( 13) is a triv ial equa tion, it will not be in cluded in the time in te gra tion and the pres sure ten dency Eq. ( 12) is sim pli fied to Eq. ( 14) un der this as sump tion.We fur ther linearize Eq. ( 14) by re plac ing p¢ in the sec ond term with a time-in de pend ent ref er ence pres sure p as: (15) The first sim pli fi ca tion is equiv a lent to as sum ing the hydro static and nonhydrostatic pres sure ten den cies are balanced in de pend ently.The fur ther linearization made in Eq. ( 15) is to for mu late a sta ble semi-im plicit scheme for in te grat ing dw/dt and dp¢/dt equa tions.By sub sti tut ing the pres sure-ten dency Eqs. ( 13) and ( 15 (1) to ( 5) form a com plete set of gov ern ing equa tions for a fully com press ible at mosphere with the only ap prox i ma tion be ing to linearize the pres sure ten dency equa tion, which is im por tant to sound waves only.We found that a choice of p to be 10% of the stan dard at mo sphere pres sure gives very ac cept able re sults both in ide al ized case sim u la tions (see below) and real data fore casts.

TIME IN TE GRA TION
For the time in te gra tion of the nonhydrostatic model developed in this ap proach, the dry dy namic part of the hydrostatic Eqs.(1) to ( 4) is first in te grated as they are done in the orig i nal hy dro static model.In our orig i nal hy dro static model, the gov ern ing Eqs. ( 1) to ( 4) are nu mer i cally solved by fourth or der cen tral dif fer enc ing in the hor i zon tal on staggered C grids and sec ond or der cen tral dif fer enc ing in the ver ti cal on stag gered sigma lev els.The ver ti cal stag ger ing car ries ver ti cal ve loc ity at the full lev els and all other variables at the half lev els.The time in te gra tion is con ducted by the leap frog scheme to gether with an ex plicit time-split scheme (Madala 1981) to han dle fast mov ing grav ity waves.The time-split scheme de com poses grav ity wave re lated vari ables, hor i zon tal di ver gence and pseudo geopotential, into ver ti cal modes.It then uses the same leap frog scheme but smaller time steps to com pute cor rec tions to the first few fast mov ing modes to main tain their sta bil ity.Rather than damp ing those fast mov ing grav ity wave modes by an implicit scheme, as is nor mally done in a semi-im plicit scheme, the Madala scheme in cludes those fast mov ing grav ity  20) are then in te grated with a suit able method to han dle fast mov ing sound waves.In our de vel opment, we use an ex plicit-split scheme to in te grate Eqs. ( 17), (18), and the hor i zon tal part Eq. ( 20) and a semi-im plicit scheme to in te grate Eq. ( 19) and the ver ti cal part Eq. ( 20).
Since our orig i nal hy dro static model is in te grated by a flux form of Eqs. ( 1) to ( 4), we will now ex press Eqs. ( 17) and ( 18) in a flux form as well.We re write the Eqs.( 17) to (20) as: where and v is the an gle be tween Y and north di rec tions.
The right hand side terms of the above equa tions are in tegrated with the time scheme used in the hy dro static in te gration, as they are in volved with Rossby and grav ity waves.How ever, the left hand side terms are in te grated with the fol low ing meth ods to han dle sound waves.A for wardback ward scheme with a smaller time step is used in in tegrat ing the sec ond and third terms of Eqs. ( 22) and ( 23) and the third term of Eq. ( 26).In this time in te gra tion, the second and third terms of Eqs. ( 22) and ( 23) are in te grated for ward and then the third term of Eq. ( 26) is in te grated back ward af ter the winds have been up dated.An im plicit scheme with the same small time step of the for wardbackward scheme is used in in te grat ing the sec ond terms of Eqs. ( 25) and ( 26).The im plicit scheme in te grates these terms by solv ing a tri-di ag o nal ma trix.Once p¢ has been up dated to the next reg u lar time step, Eq. ( 24 27) and ( 28).Once ũn + 1 and ṽn + 1 are com puted, the m term of Eq. ( 30) can be eval u ated and w n + 1 and p¢ n + 1 can be cal cu lated by a tri-di ag o nal matrix solver af ter one or the other of them is elim i nated through sub sti tu tion.We choose to elim i nate w for solv ing p¢ n + 1 with bound ary con di tions p¢ n + 1 = 0 at the model top and w = v V × Ñz s at the sur face.The linearization made in Eq. ( 15) or (20) al lows us to for mu late Eq. ( 30) with out p¢ in volved in the co ef fi cient .

COM PAR I SON WITH THE WRF MODEL
Both dy namic cores of the WRF model use the mass coor di nates but with dif fer ent prog nos tic vari ables and nu meri cal meth ods.Sim i lar to our choice, the NMM core chooses tem per a ture, rather than po ten tial tem per a ture, as a prognostic vari able for con ve nience in com put ing the gas law and ap ply ing phys ics pack ages.On the other hand, the ARW core chooses the po ten tial tem per a ture for con ve nience in en ergy con ser va tion.How ever, the ma jor dif fer ences among the two WRF dy namic cores and our method are in add ing the nonhydrostatic con tri bu tions to the prim i tive equa tions.We add w and p¢ as two ad di tional prog nos tic vari ables, and add the third mo men tum Eq. ( 19) and nonhydrostatic pressure ten dency Eq. ( 20 to accom mo date a com pli cated time in te gra tion method that splits the dy namic in te gra tion into two en ergy con serv ing sub sys tems.The method uses an im plicit scheme to solve for w, p, and T, the Ad ams-Bashforth scheme for the advection terms, the back ward scheme for the pres sure gra di ent force term, and the trap e zoidal scheme for the Coriolis term.The WRF-ARW core also adds w and to tal pres sure p as two additional prog nos tic vari ables, while it adds the third momentum equa tion and wg = d dt F to gether with the gas law to close the sys tem.It in te grates the gov ern ing equa tions by a time-split scheme sim i lar to our method de scribed in the previous sec tion ex cept: (i) it uses a third-or der Runge-Kutta scheme in the large time step in te gra tion; (ii) it considers mois ture in the den sity cal cu la tion; and (iii) it includes the pre dic tion of mass change in the small time step in te gra tion for sound waves.
As sug gested by Laprise (1992), nonhydrostatic equations for mu lated with the mass co or di nates can be solved by ei ther add ing a prog nos tic equa tion for pres sure or a prognos tic equa tion for the geopotential.Two dy namic cores of the WRF choose to add the prog nos tic equa tion for the geopotential, while we choose to add the prog nos tic equa tion for pres sure.The choice of add ing the prog nos tic equa tion for pres sure makes the terms in volved with sound waves easy to iden tify for for mu lat ing the time-split scheme.On the other hand, spe cial care must be made to en sure con sistency be tween pres sure and mass pre dic tion.We solve this prob lem by sep a rat ing the pres sure pre dic tion into mass conser va tion and p¢ pre dic tion.All three meth ods dis cussed here in clude some kind of linearization in their nu mer i cal schemes.We linearize the prog nos tic equa tion for p¢, the WRF-ARW core linearizes the gas law when it com putes pres sure per tur ba tion in the small time step in te gra tion, and the WRF-NMM core linearizes the co ef fi cient of a sec ondor der dif fer en tial equa tion solved iteratively to close the implicit scheme for w and p.The WRF-NMM core and our method use ef fi cient sec ond-or der time schemes in the time in te gra tion.Only the WRF-ARW core uses a more ex pensive higher or der time scheme in the time in te gra tion.The third-or der Runge-Kutta scheme is very ac cu rate in cal cu lating the advection to pre serve the shape of the ma te rial transported, which is im por tant for cloud and aero sol mod el ing.How ever, for NWP ap pli ca tions, as our model is de signed for, a sec ond or der time scheme is suf fi cient since ma jor numer i cal er rors come from spa tial dif fer enc ing rather than the time scheme and a rel a tively small time step is usu ally required by the CFL lin ear sta bil ity con di tion.

IDE AL IZED CASE TESTS
We have tested this method of ex tend ing a hy dro static model to a nonhydrostatic model by 2-di men sional (2D) moun tain waves and den sity flows as so ci ated with a cold bub ble.In these ide al ized tests, the nonhydrostatic model is mod i fied to be two di men sional with 481 grid points in the hor i zon tal and 87 lev els (for moun tain wave sim u la tion) or 77 lev els (for cold bub ble sim u la tion) in the ver ti cal.The model is in te grated with dry dy nam ics only.To con trol bound ary re flec tion at the model top, we treat the top 15 model lay ers as sponge lay ers where sec ond-or der ver ti cal dif fu sion is ap plied to u and T, Ray leigh type damp ing is applied to w and p¢, and sec ond-or der hor i zon tal dif fu sion is ap plied to ver ti cal dif fer ence du and dT and p¢.Ra di a tive type bound ary con di tions (Orlanski 1976) are used at lat eral bound ary points where no gra di ent is as sumed at the in flow bound ary and ex trap o la tion is made at the out flow bound ary.We treat 25 grid points next to each lat eral bound ary as a lateral bound ary zone where sec ond-or der hor i zon tal diffusion is ap plied to all prog nos tic vari ables to con trol boundary noise.In ad di tion to these bound ary treat ments, weak dif fu sion of the fourth-or der in the hor i zon tal and sec ondorder in the ver ti cal are ap plied at in te rior grid points.
In moun tain wave tests, we first re peat the 2D moun tain wave sim u la tions con ducted by Dudhia (1993) for test ing the MM5.In these nu mer i cal sim u la tions, a bell-shaped moun tain h s (x) = h x a is used.In this ex pression, a rep re sents the half width and h 0 rep re sents the max imum height of the moun tain.Five dif fer ent half-width moun tains, 100 m, 1, 10, 100, and 1000 km to gether with h 0 = 400 m are se lected to test the ac cu racy of the model in simu lat ing five dif fer ent char ac ters of lin ear moun tain waves dis cussed in Queney (1948) and Dudhia (1993).As in the MM5 sim u la tions, the grid res o lu tion Dx is 1/5 of the half width and the time step is pro por tional to the grid res o lu tion with Dt = 1 s for Dx = 200 m.Coriolis force is in cluded in all five si mulations with Coriolis pa ram e ter f = 1.0 e -4 .The initial conditions are in geostrophic and hy dro static bal ance with U 0 = 10 m s -1 zonal wind and zero me rid i o nal wind every where.Tem per a ture and height are in i tial ized with a constant strat i fi ca tion of N 0 = g z ¶q q ¶ = 0.01 s -1 , 300 K sur face tem per a ture, and 1000 hPa pres sure at ground (z = 0).These con di tions give N h U 0 0 0 = 0.4, which cor re sponds to lin ear moun tain waves (Laprise and Peltier 1989).The only difference in our sim u la tion setup is that we use 481 grid points in the hor i zon tal di rec tion and 87 lev els in the vertical.The dispersive na ture of the nonhydrostatic moun tain waves causes large am pli tude waves to prop a gate into model bound aries.As very crude bound ary con di tion treat ments ap plied at the top and lat eral bound aries, we need more grid points in the model do main to en sure the bound ary noise will be kept far away from the re gion where we are look ing for steady so lu tions.The ver ti cal res o lu tion is 13 hPa in the bottom 60 lev els and is grad u ally re duced to 7 hPa in top 16 levels.With this se lec tion the sponge layer is above 105 hPa (ap prox i mately 12.2 km.) For con ve nience in com par i son, we plot our sim u la tion re sults at 2160 time steps in the same way as that pre sented in the MM5 sim u la tions, i.e., 50 grid points in each side of the moun tain peak and up to 10 km from the ground with the 400 m moun tain at the sur face.Fig ure 1 shows the ver ti cal ve loc ity of the five moun tain wave sim u la tions from our nonhydrostatic model.The sim u la tions cap ture the char acter is tics of the lin ear moun tain waves from the o ret i cal and nu mer i cal stud ies (Queney 1948;Dudhia 1993); that is, that the wave is ev a nes cent for a = 100 m, down stream propagation for a = 1 km, up right for a = 10 km, in flu enced by the inertial force for a = 100 km, and non-sta tion ary for a = 1000 km.The flow pat tern, wave length, and wave am plitude all agree well with the MM5 sim u la tions.To ex am ine the nonhydrostatic ef fects on these steady moun tain waves, we have con ducted the same ex per i ments but with out non hy dro static con tri bu tions (Fig. 2).The re sults clearly show that nonhydrostatic ef fects are sig nif i cant mainly for cases with N a U 0 0 £ 1, as pre dicted by early the o ret i cal studies (Laprise and Peltier 1989).To check the im pact of the free-slip lower bound ary con di tion on the non hy dro static moun tain waves, we con duct an other 1 km half-width moun tain sim u la tion with vis cos ity.The vis cos ity is modeled by the sur face stress t Vï, with C d = 0.01, and free at mo sphere stress t ¶ ¶s = k V , with k = 0.001, at the bottom ¼ sigma lev els.The vis cos ity does not change the main char ac ter is tics of the moun tain wave, ex cept to sig nif icantly re duce the wave am pli tude (Fig. 3) and make the so lution con verge very slowly to a quasi-steady state (not shown).The re sult is con sis tent with the the o ret i cal stud ies in that char ac ter is tics of lin ear moun tain waves are de termined by U 0 , N 0 , a, and f, but not vis cos ity.
In the sec ond moun tain wave test, we com pare the surface pres sure per tur ba tion of sim u lated moun tain waves with that com puted by a lin ear the ory (Queney 1948;Hsu and Sun 2001).Fol low ing Hsu and Sun (2001), we use a 10 m high moun tain to in te grate our 2D nonhydrostatic model to a steady state and mul ti ply the pres sure per tur bation by 100 to com pare it with the an a lytic lin ear so lu tion for a 1000-m high moun tain.In this in te gra tion, the Coriolis force is not in cluded and the same model setup for 1-km half width with out vis cos ity is used.The in te gra tion is made for 12960 time steps to en sure a steady state so lu tion is reached.Since pres sure is a prog nos tic vari able in our model, we can eas ily cal cu late sur face pres sure per tur ba tion as the sum of the nonhydrostatic and sur face hy dro static pres sure per turba tions, i.e., dp s = (p s -p s0 + p¢) where p s0 is the ini tial hydro static pres sure at the sur face.The sim u lated sur face pressure per tur ba tion matches well with the an a lytic so lu tion com puted by Hsu and Sun (2001) and the ver ti cal ve loc ity of the lin ear wave sim u la tion is very sim i lar to the 400-m height so lu tion in Fig. 1b, but the am pli tude is about 40 times less (Fig. 4).
To sim u late den sity flows as so ci ated with a cold bub ble, we set up a nu mer i cal ex per i ment fol low ing Straka et al. (1993) and Janjic et al. (2001).The ini tial con di tions are spec i fied as an at mo sphere at rest with 1000 hPa pres sure at the ground and hy dro static bal ance above.The ini tial temperature is spec i fied with a neu tral mean strat i fi ca tion of 300 K po ten tial tem per a ture ev ery where and tem per a ture per tur ba tions of a cold bub ble, DT(x, z) = -15cos 2 (pL / 2), added at the do main cen tral area of L < 1, where x c = 0 m, z c = 3000 m, x i = 4000 m, and z i = 2000 m.The 77 ver ti cal lev els for this test are se lected to have 15 sponge lay ers above s = 0.442 (approx i mately 6.5 km) with a ver ti cal res o lu tion vary ing from 5 to 18 hPa.We use 100-m hor i zon tal res o lu tion and 0.3 s time steps in this nonhydrostatic sim u la tion and ap ply second-or der dif fu sion in both hor i zon tal and ver ti cal di rections us ing the dif fu sion co ef fi cient K = 75 m 2 s -1 as in Straka et al. (1993) and Janjic et al. (2001).The cold po ten tial temper a ture per tur ba tions from our 2D nonhydrostatic si mulations at ini tial time, 300, 600, and 900 s are dis played in Fig. 5.The area shown ex tends from the model do main center to 19.2 km to the right, and from the ground to 4.8 km.
The mi nor com pu ta tional noise on zero con tours is due to nu mer i cal er rors as so ci ated with the non-pos i tive def i nite schemes used in our model for NWP ap pli ca tions.The hor izon tal and ver ti cal wind com po nents af ter 900 s are displayed in Fig. 6 with the same set ting.The nonhydrostatic sim u la tion agrees well with the re sults from the two pre vi ous stud ies that the den sity flow ex hib its one ro tor af ter 600 s and three ro tors af ter 900 s gen er ated by the Kel vin-Helmholtz in sta bil ity.We also check the nonhydrostatic effects on this fall ing cold bub ble by con duct ing a hy drostatic sim u la tion with out the nonhydrostatic com po nents.In this hy dro static sim u la tion, we have to use a shorter time step of 0.1 s for the same grid res o lu tion due to very large ver ti cal mo tions gen er ated at ini tial pe ri ods by as sum ing the hy dro static bal ance.Fig ures 7 and 8 show the hy dro static sim u la tions that cor re spond to Figs. 5 and 6, re spec tively.
No ro tors are gen er ated in the model sim u la tion with out the nonhydrostatic com po nents.The Rich ard son num ber of the hy dro static sim u la tion af ter 300 s is about 0.7, which is too large for the Kel vin-Helmholtz in sta bil ity.Hy dro static approx i ma tion se verely dis torts the cold bub ble sim u la tion at the ini tial pe riod of fall ing when the down ward ver ti cal veloc ity is ac cel er ated by neg a tive buoy ancy.In this early period, the hy dro static sim u la tion gives too fast fall ing and too noisy den sity flows (Fig. 9), as is ex pected when the hy drostatic as sump tion is vi o lated.

SUM MARY
By in ter pret ing the hy dro static part of pres sure as air mass weight above an area, the prim i tive equa tions in the sigma-p co or di nates can be di rectly used to form the ma jor       (Leou 2004).The NFS is cur rently run ning twice (or 4 times when a ty phoon is ap proach ing Tai wan) a day with tri ple nests of 45/15/5 km res o lu tions for mesoscale and ty phoon fore casts.This method can also be used in up grad ing hy dro static global mod els to non hy drostatic global mod els.

Ac knowl edge ments
The ma jor part of this work was done be tween 1999 and 2003 dur ing many vis its by the first author to the CWB.The au thors ap pre ci ate the sup port by the for mer di rec tor gen eral, Mr. Shinn-Liang Shieh.
ver ti cal speed, p is hy dro static pres sure, p/p s , p is to tal pres sure, p = p + p¢, F = gz sa unit vec tor de fined as the cross prod uct of the unit vec tor in y-di rec tion and unit vec tor in north-di rec tion, and D 3 is three di men sional di - equa tion for en ergy con ser va tion af ter linear izing the nonhydrostatic pres sure ten dency equation.With the sim pli fi ca tion, the ex tra terms and equa tions added to the prim i tive equa tions for a nonhydrostatic model can be sum ma rized as: note the hy dro static ten den cies con tri buted by the dry dy nam ics and f¢ º F -f de notes the correction to the geopotential due to nonhydrostatic ef fects.Equa tions (17) to (21) to gether with Eqs.
waves in the model so lu tion dur ing time in te gra tion.If model in te gra tion is se lected to in clude non hy drostatic com po nents, dry dy namic hy dro static ten den cies dV dt H cal cu lated by the dif fer ence af ter and before the dy namic in te gra tion di vided by the length of the time step.Equa tions (17) to ( ) can be eas ily in te grated with dp¢/dt com puted from the p¢ change in this time step.Two or three smaller time steps per reg u lar time step are usu ally enough to han dle the sound wave in te gration.The smaller time step in te gra tion for the left hand side terms is il lus trated with the fol low ing equa tions: and where ũ rep re sents (p s u), ṽ rep re sents (p s v), h u and h v represent the sec ond and third terms of Eqs.(22) and (23), respectively, H u and H v rep re sent the right hand side terms of Eqs.(22) and (23), a and b rep re sent the co ef fi cients of ¶/ ¶s de riv a tives in the sec ond terms of Eqs.(25) and (26), I and J rep re sent the right hand side terms of Eqs.(25) and (26), and m rep re sents the third term of Eq. (26).For a given reg u lar time step N, the (n + 1) small time step in tegra tion at level k can be writ ten as: are a and b weighted by layer thick ness factors as so ci ated with the ver ti cal dif fer enc ing, and p¢ n* and w n* are weighted av er ages be tween small time steps n and (n + 1), i.e., p¢ n* = np¢ n + 1 + (1 -n) p¢ n with 0.5 £ n £ 1 rep re sent ing the de gree of im plic it ness.A de fault value of n = 0.8 is typ i cally used in the in te gra tion.As with all variables known at the time steps N and n, ũn + 1 and ṽn + 1 can be easily cal cu lated by Eqs. ( ) as two ad di tional equa tions to close the sys tem.The WRF-NMM core adds w and to tal pres sure p as two ad di tional prog nos tic vari ables.How ever, it adds three prog nos tic equa tions: the third mo men tum equa tion, pres sure ten dency equa tion, and the def i ni tion of wg = d dt F together with the con straint of = s + p s RT

Fig. 2 .
Fig.2.Same as Fig.1, ex cept for hy dro static sim u la tions with out the nonhydrostatic con tri bu tions.

Fig. 3 .
Fig.3.Moun tain wave sim u la tions with vis cos ity for 1-km half-width hill at 2160 time steps: (a) ver ti cal ve loc ity from nonhydrostatic, (b) zonal wind per tur ba tion from nonhydrostatic, (c) ver ti cal ve loc ity from hy dro static, and (d) zonal wind per tur ba tion from hy dro static runs.The zonal wind per tur ba tion is the wind de vi a tion from ini tial 10 m s -1 mean wind.Con tour in ter vals are 30 cm s -1 for ver ti cal ve loc ity and 50 cm s -1 for zonal wind per tur ba tion.

Fig. 5 .
Fig. 5. Cold po ten tial tem per a ture per tur ba tion of the cold bub ble sim u la tion with all nonhydrostatic com po nents in cluded at: (a) 0, (b) 300, (c) 600, and (d) 900 s.The con tour in ter val is 1 K.

Fig. 6 .
Fig. 6.(a) Zonal wind and (b) ver ti cal wind of the cold bub ble sim u la tion with all nonhydrostatic com po nents in cluded at 900 s.The con tour in ter val is 2 m s -1 with dashed con tours for neg a tive val ues.
gov ern ing equa tions in mass co or di nates for a fully com press ible at mo sphere with out the hy dro static approx i ma tion.Af ter add ing ex tra terms to the prim i tive equations and two more equa tions for w and p¢ pre dic tions, a hydro static model can be eas ily up graded to a nonhydrostatic model.The only sim pli fi ca tion made in this nu mer i cal method is the linearization of the pres sure ten dency equa tion, which may dis tort sound wave sim u la tions but has ne gligible im pact on weather fore casts.With this ap proach, almost all codes of a hy dro static model can be used in up grading the model to a nonhydrostatic model.In our ef forts to develop a nonhydrostatic model at the Cen tral Weather Bureau (CWB) of Tai wan, in ad di tional to mod i fy ing the model driver and out put rou tines to ac count for ex tra variables and rou tine calls, we only need to add two more routines to com pute the right hand side terms of RH ( ) for Eqs.(22) to (26) and per form smaller time step in te gra tion for Eqs.(22) to (26).The ide al ized case tests show that the pro posed method can re al is ti cally sim u late nonhydrostatic ef fects on dif fer ent at mo spheric cir cu la tions, which are revealed in the o ret i cal so lu tions and sim u la tions from other nonhydrostatic mod els.The nonhydrostatic mesoscale mo -del de vel oped by this method at the CWB, Nonhydrostatic Fore cast Sys tem (NFS), has been op er a tional since Decem ber 2003

Fig. 9 .
Fig. 9. Cold po ten tial tem per a ture per tur ba tion of the cold bub ble simulation at 75 s for (a) nonhydrostatic and (b) hy dro static sim u lations.The con tour in ter val is 1 K.

Fig. 8 .
Fig. 8. Same as Fig. 6a, ex cept for the zonal wind of the hy dro static sim u la tion with out nonhydrostatic com po nents.