1 / 49

Numeric developments in COSMO SRNWP / EWGLAM-Meeting Dubrovnik, 08.-12.10.2007

Numeric developments in COSMO SRNWP / EWGLAM-Meeting Dubrovnik, 08.-12.10.2007. Michael Baldauf 1 , Jochen Förstner 1 , Uli Schättler 1 Pier Luigi Vitagliano 2 , Gabriella Ceci 2 , Lucio Torrisi 3 , Ronny Petrik 4

baakir
Download Presentation

Numeric developments in COSMO SRNWP / EWGLAM-Meeting Dubrovnik, 08.-12.10.2007

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Numeric developments in COSMOSRNWP / EWGLAM-MeetingDubrovnik, 08.-12.10.2007 Michael Baldauf1, Jochen Förstner1, Uli Schättler1 Pier Luigi Vitagliano2, Gabriella Ceci2 , Lucio Torrisi3, Ronny Petrik4 1Deutscher Wetterdienst, 2CIRA-Institute, 3USMA (Rome), 4Max-Plank-Institut Hamburg

  2. Outlook Dynamical cores in the COSMO model: • Terrain following • Leapfrog time integration • Runge-Kutta time integration (COSMO Priority Project) • ‚operational version‘ • Stability considerations (Winter storm ‚Kyrill‘, ...) • p'T'-dynamics • Moisture advection • Deep / shallow atmosphere • Physics/Dynamics coupling • alternatives (A. Gassmann) • Semi-implicit (S. Thomas, ...) • LM-Z (COSMO Priority Project)

  3. The operational Model Chain of DWD: GME, COSMO-EU and -DE(since 16. April 2007) COSMO-EU (LME) COSMO-DE (LMK) GME hydrostatic parameterised convection x  40 km 368642 * 40 GP t = 133 sec., T = 7 days non-hydrostatic parameterised convection x = 7 km 665 * 657 * 40 GP t = 40 sec., T= 78 h non-hydrostatic resolved convection x = 2.8 km 421 * 461 * 50 GP t = 25 sec., T = 21 h

  4. COSMO - Working Group 2 (Numerics) • COSMO Priority Project 'LM-Z' • several improvements on the code: • prevent decoupling of z-grid (dynamics) and tf-grid (physics) by 'nudging' • implicit vertical advection  increase in time step • tendencies of data assimilation are now also transformed to the z-grid • Comparison of LM-Z and an older version of LM (COSMO-model) (e.g. without prognostic precipitation) • --> report: end 2007 Collaboration with Univ. of Leeds started

  5. COSMO - Working Group 2 (Numerics) • COSMO-Priority Project ‚Runge-Kutta‘: • New Developments • NEW: Divergence damping in a 3D-(isotropic) version • NEW: DFI for RK • Advection of moisture quantities in conservation form • Higher order discretization in the vertical • Physics coupling scheme • Testing of alternative fast wave scheme • Development of a more conservative dynamics (planned) • Development of an efficient semi-implicit solver in combination with RK time integration scheme (planned) • Developing diagnostic tools • Conservation inspection tool (finished) • Investigation of convergence • Known problems • Looking at pressure bias • Deep valleys • (Different filter options for orography) (finished)

  6. Numerics and Dynamics - COSMO-DE developments Grid structure horizontal: Arakawa C, vertical: Lorenz Prognostic Var. cartesian components u, v, w, p’,T’ (LME: T) Time integration time-splitting between fast and slow modes - 3-timelevels: Leapfrog (+centered diff.) (Klemp, Wilhelmson, 1978) - 2-timelevels: Runge-Kutta: 2. order, 3. order(Wicker, Skamarock, 1998, 2002) Fast modes (=sound waves, buoyancy, divergence filtering) centered diff. 2. order, vertical implicit, (p’T’-Dyn.) Advection for u,v,w,p’,T’: horizontal. adv.: upwind 3., 5. order / centered diff. 4. 6. order vertical adv.: implicit 2. orderfor qv, qc, qi, qr, qs, qg, TKE:LME: qv, qc: centered diff. 2nd order qi: 2nd ord. flux-form advection scheme qr, qs: semi-lagrange (tri-linear interpol.) Courant-number-independent (CNI)-advection: - Bott (1989) (2., 4. order), in conservation form - Semi-Lagrange (tricubic interpol.) Other slow modes(optional: complete Coriolis terms) Smoothing 3D divergence damping horizontal diffusion 4. order applied only in the boundary relaxation zone slope dependent orographic filtering

  7. Stability considerations Winter storm ‚Kyrill‘, 18.01.2007 crash of all COSMO-DE (2.8 km)-runs from 03, 06, 09, ... UTC • two measures necessary: • timestep: • old: t = 30 sec. (winter storm ‚Lothar' could be simulated) • new: t = 25 sec • time integration scheme: • old: TVD-RK3 (Shu, Osher, 1988)  • new: 3-stage 2nd order RK3 (Wicker, Skamarock 2002) 

  8. COSMO-DE (2.8 km), 18.01.2007

  9. COSMO-DE (2.8 km), 18.01.2007

  10. Von-Neumann stability analysis of a 2-dim., linearised Advection-Sound-Buoyancy-system

  11. RK3-scheme(WS2002) • upwind 5th order • Sound: =0.6 • x/ z=10 • T/ t=6 Crank-Nicholson-parameter for buoyancy terms in the p‘T‘-dynamics =0.6 =0.5 (‚pure‘ Crank-Nic.) =0.7 amplification factor =0.8 =1.0 (pure implicit) =0.9 Cadv = u T / x Csnd = cst / x  choose =0.7 as the best value

  12. What is the influence of divergence filtering ? • fast processes (operatorsplitting): • sound (Crank-Nic., =0.6), • divergence damping (vertical implicit) • no buoyancy • slow process: upwind 5. order • time splitting RK 3. order (WS2002-Version) • aspect ratio: x / z=10 • T / t=6 --> Divergence damping is needed in this dynamical core!

  13. Influence of Cdiv stability limit by long waves (k0) Cdiv = div t/x2 in COSMO-model: Cdiv = xkd * (cs * t/ x)2 Cdiv=0 Cdiv=0.025 ~0.35 amplification factor Cadv = u T / x Cdiv=0.05 Cdiv=0.1 Cdiv=0.15 Csnd = cst / x

  14. Advantages of p'T'-dynamics over p'T-dynamics 1. Improved representation of T-advection in terrain-following coordinates Terms (a) and (b) cancel analytically, but not numerically 2. Better representation of buoyancy term in fast waves solver Buoyancy term alone generates an oscillation equation:  = g/cs using T:  = a = acoustic cut-off frequency using T':

  15. point 1.): 'improved T-advection' ... Idealised test case: Steady atmosphere with mountain base state: T0, p0 deviations from base state: T', p'  0  introduces spurious circulations!

  16. Runge-Kuttaold p*-T-dynamics Leapfrog contours: vertical velocity w isolines: potential temperature 

  17. Runge-Kuttanew p*-T*-Dynamik Runge-Kuttaold p*-T-Dynamik contours: vertical velocity w isolines: potential temperature 

  18. Climate simulations • start: 1. july 1979 + 324 h (~2 weeks) • results: accumulated precipitation (TOT_PREC) and PMSL Problems: • unrealistic prediction of pressure and precipitation distribution • strong dependency from the time step These problems occur in the Leapfrog and the (old) Runge-Kutta-Version (both p'T-dynamics) but not in the semi-implicit solver or the RK-p'T'-dynamics. assumption: point 2.) 'treatment of the buoyancy term' improves this case (simulations: U. Schättler, in cooperation with the CLM-community)

  19. Leapfrog – Dt = 90s Leapfrog – Dt = 75s RR (mm/h)

  20. RK (p*-T) – Dt = 180s RK (p*-T) – Dt = 150s RR (mm/h)

  21. LF (semi-implizit) – Dt = 90s LF (semi-implizit) – Dt = 75s RR (mm/h)

  22. RK (p*-T*) – Dt = 180s RK (p*-T*) – Dt = 150s RR (mm/h)

  23. Advection of moisture quantities qx • implementation of the Bott (1989)-scheme into the Courant-number independent advection algorithm for moisture densities (Easter, 1993, Skamarock, 2004, 2006) • ‚classical‘ semi-Lagrange advection with 2nd order backtrajectory and tri-cubic interpolation (using 64 points) (Staniforth, Coté, 1991)

  24. Problems found with Bott (1989)-scheme in the meanwhile: 1.) Directional splitting of the scheme: Parallel Marchuk-splitting of conservation equation for density can lead to a complete evacuation of cellsSolution: Easter (1993), Skamarock (2004, 2006), mass-consistent splitting 2.) Strang-splitting ( 'x-y-z' and 'z-y-x' in 2 time steps) produces 2*dt oscillations Solution: proper Strang-Splitting ('x-y-2z-y-x') in every time step solves the problem, but nearly doubles the computation time 3.) metric terms prevent the scheme to be properly mass conserving <-- Schär–test case of an unconfined jet and ‚tracer=1‘ initialisation (remark: exact mass conservation is already violated by the 'flux-shifting' to make the Bott-scheme Courant-number independent)

  25. COSMO-ITA 2.8 km: comparisonRK+Bott / RK+Semi-Lagrange RK+SL for light precipitation:TS is larger, whereas FBI is smaller than that for RK+Bott. Moreover, RK+SL has slightly less domain-averaged precipitation and larger maximum prec. values than RK. L. Torrisi

  26. Semi-Lagrangian advection in COSMO-model ‚classical‘ semi-Lagrange advection (Staniforth, Coté, 1991) with 2nd order backtrajectory and tri-cubic interpolation (using 64 points) SL is not positive definite  clipping necessary 'multiplicative filling' (Rood, 1987) combines clipping with global conservation problem: global summation is not ‚reproducible‘ (dependent from domain decomposition) -> solution: REAL -> INTEGER mapping Moisture transport in COSMO model: DWD: COSMO-DE: Bott-scheme used COSMO-EU: SL scheme planned operationally MeteoCH: COSMO-S2 and COSMO-S7: SL scheme used pre-operationally CNMCA: COSMO-ITA 2.8: SL-scheme used pre-operationally

  27. Deep / shallow atmosphere Momentum equations for deep atmosphere (spherical coordinates): • additionally: • introduce a hydrostatic, steady base state • transformation to terrain following coordinates • shallow atmosphere approximation: • r ~ a • neglect terms in advection and Coriolis force • deep atmosphere terms are implemented in COSMO 3.21 • diploma thesis R. Petrik, Univ. Leipzig • White, Bromley (1995), QJRMS • Davies et al. (2005), QJRMS

  28. wmax Test case Weisman, Klemp (1982): warm bubble in a base flow with vertical velocity shear + Coriolis force dx= 2 km RR precipitation distribution ‚deep‘ (shaded), ‚shallow‘ (isolines) RRdeep- RRshallow (shaded)

  29. Case study ‚12.08.2004‘ (Diploma thesis R. Petrik) • summary for precipitation forecast in ‚deep‘, convection resolving models: • additional advection terms: not relevant • additional Coriolis terms: • have a certain influence, but don't seem to be important for COSMO-DE application • could be important for simulations near the equator

  30. Physics coupling scheme original idea: problems with reduced precipitation could be due to a nonadequate coupling between physics scheme and dynamics Problems in new physics-dynamics coupling (NPDC): • Negative feedback between NPDC and operational moistturbulence parameterization (not present in dry turbulence parameterization) • 2-z - structures in the specific cloud water field (qc) • 2-z - structures in the TKE field, unrealistic high values, where qc > 0 • Work to do: • what are the reasons for the failure of the WRF-PD-scheme in LM? (turbulence scheme?) • Test different sequences of dynamics and physics (especially physics after dynamics)  test tool (Bryan-Fritsch-case) is developed in PP ‚QPF‘, task 4.1

  31. Physics-Dynamics-Coupling n = (u, v, w, pp, T, ...)n Descr. of Advanced Research WRF Ver. 2 (2005) • Physics (I) • Radiation • Shallow Convection • Coriolis force • Turbulence ‚Physics (I)‘-Tendencies: n(phys I) + n-1(phys II) Dynamics Runge-Kutta [ (phys) +(adv)  fast waves ] * = (u, v, w, pp, T, ...)* - n-1(phys II) • Physics (II) • Cloud Microphysics ‚Physics (II)‘-Tendencies: n(phys II) n+1 = (u, v, w, pp, T, ...)n+1

  32. Plans (RK-core, short, medium range) • 3D- (isotropic) divergence filtering in fast waves solver • implicit advection of 3. order in the verticalbut: implicit adv. 3. order in every RK-substep needs ~ 30% of total computational time!  planned: use outside of RK-scheme (splitting-error?, stability with fast waves?) • Efficiency gains by using RK4? • Development of a more conservative dynamics (rho’-Theta’-dynamics?) • diabatic terms in the pressure equation (up to now neglected, e.g. Dhudia, 1991) • radiation upper boundary condition (non-local in time ) • continue diagnostics: • convergence (mountain flows) • conservation: mass, moisture variables, energy

  33. Stability limit of the ‚effective Courant-number‘ for advection schemes Ceff := C / s, s= stage of RK-scheme up1 cd2 up3 cd4 up5 cd6 Euler 1 0 0 0 0 0 LC-RK2 0.5 0 0.437 0 0 0 LC-RK3 0.419 0.577 0.542 0.421 0.478 0.364 LC-RK4 0.348 0.707 0.436 0.515 0.433 0.446 LC-RK5 0.322 0 0.391 0 0.329 0 LC-RK6 0.296 0 0.385 0 0.311 0 LC-RK7 0.282 0.252 0.369 0.184 0.323 0.159 Baldauf (2007), submitted to J. Comput. Phys.

  34. Plans (long range) Higher order discretization on unstructured grids using Discontinuous Galerkin methods Univ. Freiburg: Kröner, Dedner, NN.,   DWD: Baldauf In the DFG priority program 'METSTROEM' a new dynamical core for the COSMO-model will be developed. It will use Discontinuous Galerkin methods to achieve higher order, conservative discretizations. Currently the building of an adequate library is under development. The work with the COSMO-model will start probably at the end of 2009. This is therefore base research especially to clarify, if these methods can lead to efficient solvers for NWP. start: 2007, start of implementation into COSMO: 2009  

  35. Investigation of convergence solution with a damping layer of 85 levels and nRΔt=200. Analytical solution (Klemp-Lilly (1978) JAS)

  36. CONVERGENCE OF VERTICAL VELOCITY w L1 = average of errors L = maximum error Convergence slightly less than 2. order. (2. order at smaller scales?)

  37. NON LINEAR HYDROSTATIC FLOW Convergence of vertical velocity w L1 = average of absolute errors L = maximum error Stable and stationary solution of this non-linear case!

  38. Operational timetable of the DWD model suite GME, COSMO-EU, COSMO-DE and WAVE

  39. Equation system of LM/LMK in spherical coordinates • additionally: • introduce a hydrostatic, steady base state • Transformation to terrain-following coordinates • shallow/deep atmosphere

  40. (from spatial discretization of advection operator)

  41. RK3-scheme • slow process: upwind 5. order • aspect ratio: dx/dz=10 • dT/dt=6 • How to handle the fast processes with buoyancy? • with buoyancy (Cbuoy = adt = 0.15, standard atmosphere) • different fast processes: • operatorsplitting (Marchuk-Splitting): ‘Sound -> Div. -> Buoyancy‘ • partial adding of tendencies: ‘(Sound+Buoyancy) -> Div.') • adding of all fast tendencies: ‘Sound+Div.+Buoyancy‘ • different Crank-Nicholson-weights for buoyancy: • =0.6 / 0.7

  42. ‘Sound -> Div. -> Buoyancy‘ ‘(Sound+Buoyancy) -> Div.') ‘Sound+Div.+Buoyancy' =0.6 =0.7 Cadv = u T / x Csnd = cst / x amplification factor curious result: operator splitting of the fast processes is not the best choice, better: simple addition of tendencies.

  43. Task 3: Conservation (Baldauf) Tool for inspection of conservation properties will be developed. balance equation for scalar : temporal change sources / sinks flux divergence integration area = arbitrarily chosen cuboid (in the transformed grid, i.e. terrain-following) • Status: available in LM 3.23: • Subr. init_integral_3D: define cuboid (in the transformed grid!), prepare domain decomp. • Function integral_3D_total: calc. volume integral Vijk Vijk • Subr. surface_integral_total: calc. surface integrals V jijk * Aijk • prelimineary idealised tests were carried out • report finished; will be published in the next COSMO-Newsletter Nr. 7 (2007) Task is finished (Study of conservation properties will be continued in collaboration with MPI-Hamburg, see WG2 Task 2.10.1)

  44. Task 3: Weisman-Klemp (1982)-test case without physical parameterisation (only advection & condensation/evaporation) Semi-Lagrange-Adv. for qx with multiplicative filling x :=  (qv + qc ) total moisture mass M =  x dV (Mn-Mn-1) / t Res. violation in moisture conservation (?) total surface flux timestep

  45. Task 3: Weisman-Klemp (1982)-test case with warmer bubble (10 K) without physical parameterisation, without Condensation/Evap. Semi-Lagrange-Adv. for qx with multiplicative filling x :=  (qv + qc ) Residuum  0  advection seems to be ‚conservative enough‘ total moisture mass M =  x dV (Mn-Mn-1) / t possible reasons for conservation violation: saturation adjustment conserves specific mass (and specific energy) but not mass (and energy) itself ! Res. total surface flux Baldauf (2007), COSMO-Newsletter Nr. 7 timestep

  46. COSMO-ITA: RK+SL / RK+new Bott SL Bott RK+new Bott has a larger FBI for all precipitation thresholds than RK+SL(= COSMO-ITA operational run). Moreover, RK+new Bott has a deterioration in MSLP bias and RMSE after T+12h.

  47. Verbesserte Vertikaladvektion für dynamische Var. u, v, w, T, p‘ analytic sol. implicit 2. order implicit 3. order implicit 4. order C=1.5 80 timesteps Idealized 1D advection test C=2.5 48 timesteps

More Related