Content of presentation

1. Design of ALADIN NH with Vertical Finite Element Discretisation Jozef Vivoda, SHMÚPierre Bénard, METEO France Karim Yessad, METEO France Emeline Larrieu Rosina, SHMÚ 5-6 December, 2006

2. Content of presentation • Introduction • VFE – basic features • VFE derivative and integration operator • Discretization of linear model • Analysis of stability of 2TL ICI NESC scheme • First 2D idealized tests • Conclusions and future plans 5-6 December, 2006

3. Introduction Motivation: • VFE scheme successfully implemented into HY model (Untch and Hortal,2004) • To extent the VFE scheme into NH model with HY model as a limit case • VFE has two times higher accuracy than FD method on the same stencil (eight order for cubic elements in the case of evenly spaced mesh) • More accurate vertical velocities for SL scheme Work done so far: • Benard – study of accuracy of vertical integral operators and of compatibility of VFE with existing model choices • Vivoda – linear analysis of stability of VFE scheme • Vivoda, Yessad –start of implementation of VFE into ALADIN code 5-6 December, 2006

4. VFE scheme basic features Starting point – Untch and Hortal,2004: • Basis functions – cubic B-splines with compact support • No staggering – all variables are defined on model full levels • Only integration/derivation is performed in FE space, the products of variables are done in physical space • In SL version of HY model (ECMWF or ARPEGE) only non-local operations in the vertical are integrations. In NH version derivatives plays crucial role (structure equation contains vertical laplacian). 5-6 December, 2006

5. FE derivative operator We expand F and f in terms of chosen set of functions: Truncation error R is orthogonalized (required to vanish in weighted integral sense) on interval (0,1): Finally we incorporate the transforms from physical to FE space and back: 5-6 December, 2006

6. L+1 L+3 top L-1 L+4 L L+2 surf FE derivative operator basis functions Example for 35 unevenly spaced levels • In order to cover the physical domain (eta=<0,1>) with the full set of basis functions the 4 additional functions must be determined, the two at the bottom and the two at the top (analogical approach as Untch and Hortal, 2004). 5-6 December, 2006

7. FD FD VFE VFE FE derivative operator accuracy The 4 additional assumptions: 5-6 December, 2006

8. FE integration operator We expand F and f in terms of chosen set of functions: Truncation error R is orthogonalized (required to vanish in weighted integral sense) on interval (0,1): Finally we incorporate the transforms from physical to FE space and back: In: N values of f Out: N+1 values of integral N on full levels + integral from Model top to surface 5-6 December, 2006

9. FE integration operator basis functions Untch and Hortal, 2004 -3 -1 0 1 L L+1 L+2 L+3 L+4 5-6 December, 2006

10. Vertical operators Linear system Prognostic var. Mass weighted vertical integral operators (already in HY model) Mass weighted vertical laplacian operator (NH specific) VFE scheme – linear isothermal NH model System describing small aplitude waves around reference state 5-6 December, 2006

11. Constraint C1 on vertical operators • When eliminating the linear system to obtain single variable Helmholtz solver the following two equations for (d4,D) are obtained: Constrain C1 is 0 in continuous case: • In discrete case C1 is NOT automatically 0. The FD vertical discretization is designed in order to satisfy C1. If C1 is not satisfied further elimination is not feasible. • In order to be consistent with existing HY and NH code the C1 constrain shall be satisfied. 5-6 December, 2006

12. Constraint C2 on vertical operators When C1 is satisfied we could eliminate D from the system. It gives the structure equation With C2 constraint: • In continuous case T* is an identity matrix. In FD case T* is the tridiagonal matrix. • For stability reasons of SI time stepping algorithm we require: • L* to have real and negative eigenvalues • T* to have positive and real eigenvalues • In the case of C1 unsatisfied above conditions could be insufficient, we we still require them assuming that C1 is close to zero 5-6 December, 2006

13. Laplacian term FE treatment In linear and nonlinear model the laplacian has the same form (thanks to fact that we use vertical divergence related prognostic variable): V2: V1: • Boundary conditions the same as in FD model version 5-6 December, 2006

14. Imaginary part of eigenvalues Real part of eigenvalues Laplacian term FE treatment V2: V1: Eigenvalues of laplacian for 100 levels Green – eigenvalues of V1 Red – eigenvalues of V2 Imaginary part multiplied by 100 5-6 December, 2006

15. Stability of 2TL Non-Extrapolating Iterative Centered Implicit scheme Predictor: Nth-corrector: For the purpose of analysis of stability we linearise the nonlinear residual R around resting, isothermal, hydrostatically balanced reference state: We analyse the temperature instability due to fact that The following analyses were computed for the wave with dx=5km, dt=600s, 35 vertical levels, T*=350K, T*a=100K, p*s=101325Pa (perfect pressure) 5-6 December, 2006

16. Solid – VFE Dotted - FD SI N=1 N=2 N=3 C1 constrain – modification of G* • from C1 constrain we redefine the operator G*. It must be computed via iterative procedure. In linear and nonlinear model In linear model only 5-6 December, 2006

17. Stability – C1 satisfied for NL model as well Stability – C1 satisfied in linear model only Solid – VFE Dotted - FD SI N=1 N=2 N=3 C1 constrain – modification of S* • Directly from C1 constrain we can define operator S* The time stepping is unstable, because eigenvalues of L*A2 are complex. 5-6 December, 2006

18. Eigenvalues of C2 constraint Black – no C1 Red - S* redefinition in lin and nonlin Green – G* redefinition in lin and nonlin Imaginary part multiplied by 10 Circle radius 0.4 Solid – VFE Dotted - FD SI N=1 N=2 N=3 C1 constrain unsatisfied Dt=600s Dt=999999s • Solution of 2Lx2L soler for each spectral component – noncompatibility with existing model (LxL solver), problem with efficient introduction of map factor horizontal dependence (LSIDG resp. LESIDG key) • Yessad – proposal of iterative procedure, build from existing pieces of code, study of convergence required 5-6 December, 2006

19. VFE scheme – first tests • Eta coordinate • NLNH flow regime • dx=200m • V=10ms-1 • Agnesi hill, h=1000m, L=1000m • CY30T1 • Tested in 2D • 2TL ICI NESC n=1 • sponge Ref VFE, 100lev, dt=1s (X-term, Z-term FD) Ref FD, 100lev, dt=5s 5-6 December, 2006

20. VFE scheme – first tests Ref VFE, 35lev, dt=5s X-term, Z-term-FD Ref FD, 35, dt=5s 5-6 December, 2006

21. VFE scheme – first tests Ref VFE, 35lev, dt=5s X-term FD, Z-term-VFE Ref VFE, 35lev, dt=5s X-term VFE, Z-term-FD 5-6 December, 2006

22. VFE scheme – conclusions • Analysis of stability in linear framework suggests that it is possible to extend HY model VFE scheme to NH model • From stability point of view the crucial is the definition of vertical laplacian operator (all eigenvalues must be real and negative) • The VFE scheme is stable in linear context when FD boundary conditions are used in VFE (we can adopt BBC and TBC from FD model), but oscillatory behaviour near model bottom is observed (this is not true for large number of levels !) • So far no acceptable stable solutions that satisfy C1 was found 5-6 December, 2006

23. VFE scheme – future work • Non-oscillatory discretization of vertical laplacian • Extension of just described VFE discretization concept into full NL model • Coding of “draft” of VFE scheme in NH ALADIN with 2Lx2L solver (finished) • testing 5-6 December, 2006

24. Thank you for your attention ! Ďakujem za pozornosť ! Merci de votre attention ! Hvala ! 5-6 December, 2006