1 / 48

Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany

COSMO Priority Project: Further developments of the Runge-Kutta Time Integration Scheme report ‘Oct. 2007 – Sept. 2008’ / final report COSMO General Meeting, Cracow 15.-19.09.2008. Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany. Tasks of the Priority Project ‚Runge-Kutta‘:

bianca-atkins
Download Presentation

Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany

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. COSMO Priority Project:Further developments of the Runge-Kutta Time Integration Scheme report ‘Oct. 2007 – Sept. 2008’ / final reportCOSMO General Meeting, Cracow15.-19.09.2008 Michael Baldauf Deutscher Wetterdienst, Offenbach, Germany

  2. Tasks of the Priority Project ‚Runge-Kutta‘: Repair detected model deficiencies: 1. Looking at pressure bias 4. Advection of moisture quantities in conservation form 6. Deep valleys 7. Different filter options for orography 14. DFI for RK New developments: 8. Higher order discretization in the vertical for Runge Kutta scheme 9. Physics coupling scheme 10. Testing of alternative fast wave scheme 13. Divergence damping in a truely 3D-version Tool development: 2. Continue RK case studies 3. Conservation tool 5. Investigation of convergence 11. Development of a more conservative dynamics (planned) 12. Development of an efficient semi-implicit solver in combination with RK time integration scheme (planned)

  3. Task 1: Looking at pressure bias (Torrisi, Zängl) Goals: verifications of LM 7 km runs showed a higher positive pressure bias for the RK core than for the Leapfrog core, whereas other variables show comparable behaviour. Reasons and solutions? Talk by Lucio Torrisi RK starting point of the task: Leapfrog

  4. Task 2: Continue RK case studies (deMorsier,Torrisi) Identify problems of the RK scheme Several unstable cases found in previous winter periods (e.g. ‚13. Jan. 2004‘) most of them could be simulated with Semi-Lagrange Adv. for moisture variables Winter storms Kyrill ('18.01.2007') and Lothar ('26.12.1999') simulated with MeteoCH new pre-operational model chain (2.2 km and 6.6 km):  new configuration: 1.) WRF-like RK3 used (instead of TVD-RK3) (as found at DWD for the Kyrill case)2.) Semi-Lagrange-Adv. for moisture (instead of Bott-scheme)3.) new level distribution especially in boundary layer (cures problems with TKE scheme)

  5. TP BL Vertical level distribution • Test-chains in July 2007 using operational (L60.2) and a new (L60.1) vertical level distributions • Three test cases • 12.7.2006 (convection) • 23.12.2006 (fog) • 18.1.2007 (Kyrill)

  6. TKE Instability 23.12.2006 (fog) gridpoint (60,180) L60 v2 L60 v1 • Strong checkerboard instability in TKE-field • Same effect also 12.7.2006 (convection) and 18.1.2007 (Kyrill)

  7. Stability of TKE-diffusion O. Fuhrer Alternative: vertically implicit (Crank-Nicholson) scheme was implementedthis cures the most problems; some artefacts remain (stability functions?) If equation is solved explicitly, stability constraints apply In COSMO, the diffusion constant is limited Default value for securi = 0.85 is wrong!!!

  8. 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)

  9. Task 4: Advection of moisture quantities in conservation form (Förstner, Baldauf) • Status: two schemes available • implementation of the Bott (1989)-scheme into the Courant-number independent advection algorithm for the moisture densities with mass consistency (Easter, 1993, Skamarock, 2004, 2006) • Task was finished in Sept. 2006 because implemented • schemes (Bott-2, Bott-4) behaved well • But in the meanwhile: stability problems occured in some cases (steep orography!) • revival of the task necessary! • Semi-Lagrange-scheme (backtraj. 2nd order, tri-cubic interpolation)multiplicative filling algorithm for global conservation

  10. Task 5: investigation of convergence • (Ceci, Vitagliano) • Goals: determination of the spatial and temporal order of convergence of the RK-scheme in combination with advection schemes of higher order. • Test cases: • linear, 2D, hydrostatic mountain flow (h=10 m, a=10 km) • linear, 2D, non-hydrostatic mountain flow (h=10 m, a=500 m) • nonlinear, 2D mountain flows (dry case) (h=500 m, a=10 km) • linear, 3D mountain flow • nonlinear mountain flows with precipitation

  11. Deliver a program to calculate linear analytic solutions (e.g. for convergence tests) (M. Baldauf) • Starting point: compressible Euler equations • Preconditions: • no friction • only adiabatic processes (in particular no phase changes) • ideal gas law • cp=const., cV=const., R=const. • no Coriolis force • all movements take place on a plane (no earth curvature) • these preconditions can easily be fulfilled by a dynamical core (‚switches‘) • Only 2 approximations will be made: • linearisation (1/Fr<<<1  very flat mountains; not too small U) • the assumption that kz=const (see below; not absolutely necessary) • confidence into the accuracy of the linear solution for comparison with numerical models

  12. Stationary case (=0) From perturbation equations: express u', v', ' and p' by w'  equation 2nd order for w'(kx, ky, z): with coefficient functions: The only approximation so far is linearisation!

  13. Example 2: 2D-test case from Schaer et al (2002) w [m/s]

  14. w [m/s] Example 3: 3D Gaussian Hill

  15. G. Zängl Initialization of the perturbation pressure field • The present initialization of the perturbation pressure field (executed in src_artifdata for idealized simulations; otherwise in int2LM) is not exactly consistent with the discretized buoyancy term in the vertical momentum equation • The error is too small to be noticeable in real-case applications; however, it becomes evident in idealized simulations with constant flow and a very low mountain (or no mountain at all) • To fix the problem, a new initialization procedure has been developed by solving the discretized vertical wind equation (for dw/dt = 0) for p‘; ideally, this would ensure strict absence of buoyancy at the lateral model boundaries

  16. Simulation with flat surface, u = 10m/s, and fixed relaxation b.c.‘s, t = 12 h Fields: θ (contour interval 2 K), w (colours) Old p‘ initialization Error amplitude: 1 mm/s New p‘ initialization Error amplitude: 10-4 mm/s

  17. Gaussian mountain height=750 m size=10 km Horizontal resolution 4 km 3D TEST CASES: HYDROSTATIC FLOW P. L. Vitagliano, G. Ceci 3th order upwind 5th order upwind

  18. Gaussian mountain height=750 m size=10 km Horizontal resolution 8 km 3D TEST CASES: HYDROSTATIC FLOW 3th order upwind 5th order upwind

  19. Gaussian mountain height=750 m size=10 km Horizontal resolution 16 km 3D TEST CASES: HYDROSTATIC FLOW 3th order upwind 5th order upwind

  20. P. L. Vitagliano, G. Ceci CONVERGENCE OF KINETIC ENERGY

  21. CONVERGENCE OF KINETIC ENERGY

  22. Conclusions from Convergence Tests: 2D mountain flow • slightly less than 2nd order spatial convergence (fast waves scheme dominates) • TVD and non-TVD 3 stages RungeKutta show similar behaviour • time step has minor effect (if any) on spatial convergence • important issues with upper and lateral boundary condition • difficult to compare with analytical solutions, due to b.c. 3D mountain flow • smaller influence of damping layer on 3d mountain waves and drag • optimal damping parameter dt*nrdtau increases to 1000 s • with poor resolution different scheme can give different solutions • with poor resolution higher order upwind can improve results

  23. Task 6: deep valleys Goal:detection of the reason for the unrealistic ‚cold pools‘ in Alpine valleys + Task 7: Different filter options for orography (Förstner, Torrisi, Reinhardt, deMorsier) Status:The reason for the cold pools was identified: metric terms of the pressure gradient Dynamical Bottom boundary condition (DBBC) (A. Gassmann (2004), COSMO-Newsl.) and a slope-dependent orography-filtering cures the problem to a certain extent. example: For the COSMO-DE first the orography is filtered globally to remove scales approximately smaller than 4-. In a second step a stronger filter (5-) is used for all points with a step of the orography still bigger than 625 m.

  24. (G. Zängl) Spurious noise over mountains in a resting atmosphere • Tests reveal a 2Δz structure in the horizontal and vertical wind field • Depending on the difference between base state and actual temperature profile, it can take more than 12 h until the noise reaches a significant amplitude • Afterwards, it rapidly grows within a time scale of a few hours until some sort of saturation is reached • Tests indicate that a modified discretization of the dw/dz term in the pressure tendency equation may damp the noise Setup of test experiments: mountain with h = 1500 m, a = 5 km; Δx = 1 km, no ambient winds; results are shown for t = 24 h

  25. Results with implicit 2nd-order vertical advection θ (contour interval 1 K), u (colours) standard discretization with damping discretization

  26. Results with implicit 2nd-order vertical advection θ (contour interval 1 K), w (colours) standard discretization with damping discretization

  27. Results for quasi-linear flow over a mountain, h = 300 m, u = 10 m/s θ (contour interval 1 K), u (colours) standard discretization with damping discretization

  28. Spurious noise over mountains in a resting atmosphere • In the modified version, the term is not only evaluated between half-levels but also between full-levels (which damps 2Δz waves), followed by a weighting of both terms • A weight of 0.05 of the damping discretization turned out to suffice for eliminating the noise • Normally very small impact on flow dynamics, butstability problems over steep topography in the presence of strong winds

  29. Task 8: Higher order discretization in the vertical for RK-scheme (Baldauf) Improved vertical advection for the dynamic var. u, v, w, T (or T‘), p‘ motivation: resolved convection • vertical advection has increased importance => use scheme of higher order (compare: horizontal adv. from 2. order to 5. order) • => bigger w (~20 m/s) => Courant-crit. is violated =>implicit scheme or CNI-explicit scheme up to now: implicit (Crank-Nicholson) advection 2. order (centered differences) new: implicit (Crank-N.) advektion 3. order  LES with 5-banddiagonal-matrix but: 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?) Work to do: best combination with time integration scheme?

  30. Comparison of the two implicit vertical advection schemes Test with constant vertical velocity; initial cone distribution implicit cent. diff. 2nd order implicit cent. diff. 3rd order

  31. Task 9: Physics coupling scheme (deMorsier, Förstner) 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) (=WRF-like coupling): • 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: Is the problem cured now also in the moist turbulence case with the improvement of the TKE-Diffusion (solution of O. Fuhrer)?

  32. Task 10: Testing of alternative fast wave scheme • (Torrisi, Gassmann) • Goals: • p‘T‘-RK-scheme • ‚shortened-RK2‘-scheme (Gassmann) • this allows the use of the ‚radiative upper boundary condition‘ (RUBC) • Properties of A. Gassmann dyn. core: • Splitting up of vertical advection of p*/T into fast/slow mode equations and consistent boundary conditions • Vertical average to half levels: mass weighted mean (in RK simple mean) and base-state consistent formulation of the discrete w-equation • Different horizontal pressure gradient discretization • Divergence in conservative flux-form • Slightly different buoyancy term • No artificial divergence damping

  33. Status: • The fast waves part (Gassmann) is combined with the Leapfrog scheme in LM 3.21 • Original Gassmann dynamical core poses stability problems in several cases! • Gassmann fast waves part in RK3 worked in only 1 case • ‚shortened RK2‘-scheme (Gassmann (2002), Gassmann and Herzog (2007)) is implemented into LM 3.21 using the fast waves solver of RK3 and the RK3 advection/physics subroutines • Preliminary investigation of this dynamical core (L. Torrisi)tested in real cases for a five days period: similar results to the RK3 splitting method • Separate inspection of divergence in conservation form and vertical staggering • Implementation questions pointed out: • Splitting of contravariant vertical velocityposes problems in formulation of lower boundary conditions Overall assessment: needs too much work to bring to operational use

  34. Task 13: Divergence damping in a true 3D-version (Baldauf) Description: Cases occured, where the up to now used 'quasi-2D' divergence filtering lead to unstable results. But a complete abandoning of the divergence filtering (as proposed by A. Gassmann for her dynamical core) also leads to several instabilities. This was also shown by stability analyses of the RK-core by M. Baldauf. P. Prohl (DWD) could demonstrate, that the Bryan-Fritsch- test case of a rising warm bubble is unstable with 'quasi-2D' divergence damping but becomes stable only with a full 3D (=isotropic) version (realised with a preliminary explicit formulation). For operational use an implicit version of 3D divergence damping is necessary.

  35. Task 14: DFI for RK (L. Torrisi) An initialization scheme Twice Digital Filter Initialization Adiabatic backward integration Diabatic forward integration *

  36. Digital Filter Initialization in RK core (L. Torrisi) • Some modifications (mostly in the adv functions ) are needed to run DFI with RK core. • They are in: • - dfi_initialization.f90: add initialization of rho_snow • - src_runge_kutta.f90: correction in wind Rayleigh damping • - src_advection_rk.f90: changes in cfl control and changes in adv function interfaces • - fast_waves_rk.f90: changes in adv function interfaces • - numeric_utilities_rk.f90: changes in adv function interfaces and correction to run with DFI • All the odd order advection operators are changed to run in the backward integration of the DFI. The odd order advection operators implicitly contain a dissipative term that needs a special treatment in the backward integration of the DFI. • The dissipative terms are treated as the horizontal diffusion operator in the backward integration of DFI (when dt<0 , -1 is multiplied to the dissipative term). • For example: • 5th order velo*ds/dx operator = 6th order velo*ds/dx operator + dissipative term * SIGN(1.,dt) Sign(1.,dt) *

  37. Digital Filter Initialization in RK core • DFI seems to work well using a 7 km grid spacing 2h 2h

  38. COSMO-IT (2.8km) Digital Filter Initialization in RK core • Using a 2.8km grid spacing DFI works only with explicit vertical advection 1h 1h

  39. (M. Baldauf) Optimization of horizontal advection: up to COSMO 4.3: 'advection operators' = a subroutine acting on every single grid point  compiler has problems to optimize loops since COSMO 4.4: advection routines using 'field operations' (and additionally the DFI modifications of Lucio Torrisi) • Efficiency gain for routine COSMO-DE at DWD (IBM): • speedup of the horizontal advection alone: ~ 3 times faster • overall reduction of model run time:~ 1 Min. / 20 Min. ~ 5% Furthermore, some inconsistencies using metrical factors could be repaired: acrlat(j,1)acrlat(j,2) lent to an error of ~ -0.05% in the term v dw/dy

  40. Introduction of RK-scheme into operational models DWD: COSMO-DE (2.8km): since 16.04.2007 COSMO-EU (7km): planned for ~Q4/2008 (if pressure bias problems removed) (weak artificial horizontal diffusion, SL-scheme, new aver. reference pressure) MeteoCH: COSMO-S2: operational since April 2008 COSMO-S7: CNMCA: COSMO-IT (2.8km): since Oct. 2007 COSMO-ME (7km): in next future

  41. Publications of the PP 'Runge-Kutta' • Reviewed articles • M. Baldauf (2008): Stability analysis for linear discretisations of the advection equation with Runge-Kutta time integration, J. Comput. Phys. 227, 6638-6659 • Other articles • M. Baldauf (2008): A Tool for Testing Conservation Properties in the COSMO-Model (LM), COSMO-Newsletter 7, 7-17 • J. Förstner, M. Baldauf, A. Seifert (2006), Courant Number Independent Advection of the Moisture Quantities for the LMK, COSMO-Newsletter 6, 51-64 • L. Torrisi (2006): Sensitivity experiments with the Runge-Kutta time integration scheme, COSMO-Newsletter No. 6 Final report: Draft version (12.09.2008) available

  42. Thanks to all contributing scientists (in alphabetical order): Michael Baldauf1, Gabriella Ceci2, Jochen Förstner1, Oliver Fuhrer4, Almut Gassmann5, Hans-Joachim Herzog1, Guy deMorsier4, Thorsten Reinhardt 7, Gdaly Rivin6, Lucio Torrisi3, Pier Luigi Vitagliano2, Günther Zängl1 1Deutscher Wetterdienst (DWD), Germany 2 Centro Italiano Ricerche Aerospaziali (CIRA), Italy 3 Centro Nazionale di Meteorologia e Climatologia Aeronautica (CNMCA), Italy 4 MeteoSchweiz, Switzerland 5 Max-Planck-Insitut, Hamburg, Germany 6 Federal Service for Hydrometeorology and Environmental Monitoring, Russia 7 Universität Köln, Germany

  43. List of people contributing to the project during Oct. 2007 - Sept. 2008: • (alphabetical order) • Michael Baldauf (DWD, D) • Gabriella Ceci (CIRA, I) • Oliver Fuhrer (MeteoCH, CH) • Lucio Torrisi (CNMCA, I) • Pier Luigi Vitagliano (CIRA, I) • Gdaly Rivin (Roshydromet, RU) • Günther Zängl (DWD,D) Additional meeting of PP-RK-Group during the LM-User-Workshop, Langen, 05.03.2008

  44. 04.09.2007 Task 3: Weisman-Klemp (1982)-test case without physical parameterisation (only advection & Condensation/Evap.) 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. Convergence tests G. Ceci, P. L. Vitagliano HYDROSTATIC LINEAR / NON LINEAR a = 10 km H = 10 m / 500 m Time = 60 h / 100 h dt = 2.5” Domain size 500x19.5 km2 Horizontal resolution = 4km, 2km, 1km, 500m, 250m, 125m NON HYDROSTATIC a = 500 m H = 10 m Time = 100 h dt = 2.5” Domain size 250x19.5 km2 Horizontal resolution = 1km, 500m, 250m, 125m, 62.5m • All test cases runned again with constant time step = 2.5” • Test cases repeated with non-TVD 3-stage RK

  46. P. L. Vitagliano, G. Ceci CONVERGENCE OF VERTICAL VELOCITY w

  47. CONVERGENCE OF VERTICAL VELOCITY w

More Related