Numerical Modelling of Weather and Climate
Thomas Kuster and Andr´ e Welti 3. Juli 2007
1 Einf¨ uhrung
Mit einem simplen isentropischen Model (isentrop.f) werden verschiedene Effekte studiert. Das Model wird in Fortran implementiert, die berechneten Daten werden als ASCII-Dateien ausgeschrieben und k¨ onnen in Matlab mit zur Veerf¨ ugung gestellten Funktionen visualisiert werden.
2 Modell
2.1 Vereinfachungen
Im Modell werden einige Vereinfachungen angenommen. Die Rotation der Erde wird vernachl¨ assigt, d. h. es gibt keine Corioliskraft:
f = 0
Nur adiabatische Fl¨ usse, d. h. die potentizelle Temperatur (θ) bleibt erhalten:
θ ˙ = Dθ Dt = 0
Nur einen zweidimensionalen Fluss in der x-z-Ebenene:
∂
∂y = 0 v = 0
Die untere Grenzschicht ist eine isentrope Fl¨ ache, d. h. die potentielle Tem- peratur der untersten Schicht ist konstant.
θ(z = z s ) = θ s = konstant
2.2 Gleichungen und ihre Diskretisierungen
Die dem Modell zu Grunde liegenden Gleichungen werden diskretisiert und als Fortran Code implementiert. F¨ ur die Integration werden zentrierte Diffe- renzen verwendet. Die Gleichung wird dann nach dem n¨ achsten Zeitschritt gel¨ ost und als Fortran Code implementiert
2.2.1 Gleichungsset
Horizontale Impulsgleichung in x-Richtung:
Du
Dt = − ∂M
∂x θ
(1) mit dem vereinfachten advektions Operator:
D
Dt = ∂
∂t + u ∂
∂x
!
θ
(2) und dem Montgomerypotential:
M = φ + c p T = gz + c p T (3) Die zweidimensionale Kontinuit¨ atsgleichung:
∂σ
∂t + ∂σu
∂x
!
θ
= 0 (4)
mit der isentropen Dichte:
σ = − 1 g
∂p
∂θ (5)
Die hydrostatische Beziehung:
π = ∂M
∂θ (6)
mit (Exner-Funktion):
π = c p p p ref
! R/c
p(7)
2.3 Diskretisierungen in Fortan Code
variablenew steht f¨ ur die Variable im n¨ achsten Zeitschritt (n+1), variablenow f¨ ur diejenige im aktuellen Zeitschritt (n).
2.3.1 Isentropische Massen Dichte, Gleichung (4)
1 ! time step for i s e n t r o p i c mass d e n s i t y
2 do k =1 , nz
3 do i =1 , nx
4 snew ( i , k ) = sold ( i , k ) - dtdx *(
5 & 0 . 5 * ( unow ( i +2 , k ) + unow ( i +1 , k ) ) * snow ( i +1 , k )
6 & -0.5*( unow ( i , k ) + unow ( i -1 , k ) ) *
snow ( i -1 , k ) )
7 e n d d o
8 e n d d o
2.3.2 Geschwindigkeit im n¨ achsten Zeitschritt
Die Ortsindices sind alle genau zwischen dem Gitter und daher Halbzahlig um trotzdem ganze Zahlen zu haben (Arrayindices m¨ ussen ganzzahlig sein), werden diese Indices alle auf die n¨ achste ganze Zahl aufgerundet.
1 ! time step for v e l o c i t y
2 do k =1 , nz
3 do i =1 , nx +1
4 unew ( i , k ) = uold ( i , k )
5 & - dtdx * unow ( i , k ) *( unow ( i +1 , k ) - unow ( i -1 , k ) )
6 & -2.* dtdx *( mtg ( i , k ) - mtg ( i -1 , k ) )
7 e n d d o
8 e n d d o
2.3.3 Druckberechnung
Der Druck wird von oben nach unten berechnet. Es wird angenommen dass
der Druck an der oberen und unteren Randen am Anfang konstant ist. Die
Druckwerte werden mit der Gleichung (5) von oben nach unten berechnet:
1 ! P r e s s u r e d i a g n o s t i c s ( u p p e r b o u n d a r y c o n d i t i o n
2 ! and i n t e g r a t i o n d o w n w a r d s )
3 do i = -1 , nx +2
4 prs ( i , nz +1) = prs0 ( nz +1)
5 do k = nz ,1 , -1
6 prs ( i , k ) = prs ( i , k +1) + g * dth * snow ( i , k )
7 e n d d o
8 e n d d o
2.3.4 Montgomery Potential
Um die das Montgomery Potential zu berechnen wird zuerst die Exner- Funktion, Gleichung (7), diskretisiert:
1 ! C a l c u l a t e E x n e r f u n c t i o n
2 do i = -1 , nx +2
3 do k =1 , nz +1
4 exn ( i , k ) = cp *( prs ( i , k ) / pref ) ** rdcp
5 e n d d o
6 e n d d o
Die unteren Randwerte sind gegeben, daher kann das Montgomery Potential, Gleichung (6), aufw¨ arts integriert werden:
1 ! M o n t g o m e r y p o t e n t i a l d i a g n o s t i c s ( l o w e r b o u n d a r y c o n d i t i o n
2 ! and i n t e g r a t i o n u p w a r d s )
3 do i = -1 , nx +2
4 mtg ( i ,1) = g * topo ( i ) * t o p o f a c t + th0 (1) * exn ( i ,1)
5 mtg ( i ,1) = mtg ( i ,1) + dth * ( 3 . * exn ( i ,1) + exn ( i ,2) ) /8.
6 do k =2 , nz
7 mtg ( i , k ) = mtg ( i , k -1) + dth * exn ( i , k )
8 e n d d o
9 e n d d o
2.4 Anfangsbedingungen
Die Geschwindigkeit zum Zeitpunkt 0 ist u 0 , bzw. ein vertikales Windprofil kann angegeben werden.
Die Schichtdicke, Gleichung (5) ist Konstant:
σ = − 1 g
∂p
∂θ = 0 (8)
Die Topographie wird wie folgt beschrieben:
a · exp −
x i b
2 !
(9) a ist H¨ ohe des Berges.
Um die vertikale Schichtung zu definieren wird die Brunt-V¨ ais¨ al¨ a-Frequenz benutzt:
N 2 = g θ
∂θ
∂z (10)
Die Anfangsschichtung ist gleichm¨ assig (N (x, z, t = 0) = N 0 ).
3 Visualisierung in Matlab
F¨ ur die Visualisierung in Matlab stehen bereits erstellte Funktion zur Verf¨ ugung.
3.1 Standard Werte
Die f¨ ur die Simulation verwendeten Standardwerte, in Tabelle 1, sind typisch f¨ ur synoptische Systeme der mittleren Breite und h¨ aufigen orographischen Hindernissen. Wird isentrop.f mit den Standardwerten aufgerufen ergibt sich f¨ ur die Zeit t=4 h der Plot in Abbildung 1.
Bei der ¨ Uberstr¨ omung von Hindernis z. B. einem Berg wird bei stabilen Bedingungen die Bildung von Wellen angeregt. Die Wellen propagieren in vertikaler Richtung und reduzieren die Stabilit¨ at der Schichtung durch Auf- wind und erzeugen Windsch¨ arung.
Die Isotropenfl¨ ache ist in Regionen wo Wellen entstehen geneigt. Bei der
Simulation der Str¨ omung w¨ ahrend einer Dauer von 10 Stunden zeigt sich ab
der neunten Stunde eine deutliche Entwicklung von einer station¨ aren Gravity
Wave ¨ uber dem Berg (Abbildung 2).
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 1: Plot zum Zeitpunkt t = 4 h mit Standardwerten. Die d¨ unnen Linien sind isentropen Linien (Orte mit der selben potentieller Temperatur, θ). Die dickeren Linien sind die horizontalen Geschwindigkeiten (gepunktet:
negativ, ausgezogen: positiv)
Tabelle 1: Standardwerte
Bezeichung Wert Einheit
Anfangsgeschwindigkeit 10 m/s
Brunt-V¨ ais¨ al¨ a Frequenz 0.01 Hz
Bergh¨ ohe 500 m
Bergbreite auf halber H¨ ohe 50 km
Bergwachstumsrate 1800 s
Laterale Grenzen periodisch Boolean
Anzahl der Gitterpunkte im Absorber 30
Maximale Absorbtion 1
Simulationsbreite 500 km
Simulationsh¨ ohe 60 K
Zeitschritt 10 s
Integrationszeit 4 h
Diffusionskoeffizient 0.02
Ausgabeintervall NTS/8
0 50 100 150 200 250 300 350 400 450 500
0 1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 2: Plot zum Zeitpunkt t = 9 h mit Standardwerten. Station¨ are
Gravity Wave ¨ uber dem Berg
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 3: Periodische St¨ orung pflanzt sich ¨ uber den Rand fort.
3.2 Probleme
Die x-Werte sind periodisch, dadurch k¨ onnen sich St¨ orungen bei l¨ angeren In- tegrationszeiten ¨ uber die rechte Seite fortpflanzen und erscheinen links wie- der, wie dies in Abbildung 3 der Fall ist.
3.3 Absorption
Durch eine D¨ ampfung der Werte am lateralen Rand, wird eine Propagation von St¨ orungen ¨ uber die R¨ andern verhindert. Dies wird in der Funktion relax implementiert. In der Abbildung 4 ist die deutliche Verbesserung zu sehen gegen¨ uber der Model mit periodischen R¨ andern zu sehen (Abbildung 3).
3.4 Stabiler Zustand
Werden die horizontalen Geschwindigkeiten ¨ uber dem Berg, gegen die Zeit
geplottet (Abbildung 5) ist ersichtlich, dass sich nach einiger Zeit ein stabi-
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 4: D¨ ampfung der Werte an den R¨ andern, keine periodische
St¨ orung welche sich ¨ uber den Rand fortpflanzt.
5 4 5
5
5 5
6 6
6
6
6 6
7
7 7
7 7
7
8
8 8
8
8
8 8
9
9
9 9
9
9 9
9 9
10
10
10 10
10 10
10
10
10 11 10 11
11 11
11
11 11
11 11
12 12
12 12
12 12
12
12
12
12 13
13
13 13
13
13 13
13
14
14 14
14
14 14
14
15 15
15
15
15 15
16 16 Velocity at x = 50
Time [h]
Height [km]
10 20 30 40 50 60 70 80 90 100
2 4 6 8 10 12 14 16
Abbildung 5: Ver¨ anderung der Geschwindigkeit ¨ uber dem Berg mit der Zeit, bei ged¨ ampften R¨ andern
ler Zustand einstellt (keine Ver¨ anderungen der Geschwindigkeit). Wird die
” mountain groth time
” erh¨ oht, wird der stabile Zustand viel schneller er- reicht, da die Abweichung zum stabilen Windfeld jeweils viel kleiner ist.
Bei periodischen R¨ andern kann sich kein stabiler Zustand einstellen (Ab- bildung 6).
Ebenso bildet sich bei sehr hohen Windgeschwindigkeiten keine Gravita- tionswelle aus, da ein laminares Feld entsteht, kleine St¨ orungen werden sofort mit dem Wind wegtransportiert.
3.5 Scherwinde
Gravitationswellen k¨ onnen eine Scherzone nicht durchdringen dadurch kommt
es zu Resonanzerscheinungen und die Windgeschwindigkeit nimmt stark zu
4
4 4
4
5
5 5 5 55
5
6
6 6 6 6
6 7 6
7 7
7
7 7
8
8 8
8 8
8 8
9
9 9
9
9
9 9
9 9
9 9
9 9
10 10
10 10
10
10 10
10 10
10 10
11 11
11 11
11 11
11
11 11
11 11
12 12 12 12
12 12
12 12
12
12 12
12
13
13 13
13 13
13 13 13
14
14 14
14
14
14
15
15 15
15
16 16
Velocity at x = 50
Time [h]
Height [km]
10 20 30 40 50 60 70 80 90 100
2 4 6 8 10 12 14 16
Abbildung 6: Ver¨ anderung der Geschwindigkeit ¨ uber dem Berg mit der Zeit,
bei periodischen R¨ andern
17 18 19
19
19 19 19
20 20
20
20 20
20 20
20
20 20
21
21 21
21
22
22 22
22
23 23 23
23
24
24
24
24
25
25 25
25
26 26 26
26 27 27
27
27 28 28
28
28
29 29
29
29
30
30
30
30 31 31
31
31
32
32 32
32
33
33
33
33
34 34 34
34
35
35 35
35
36
36 36
36
37
37
37
38
38
38
39
39
39 40
40
40
Velocity at zlev = 8
x [km]
Time [h]
50 100 150 200 250 300 350 400 450 500
10 20 30 40 50 60 70 80 90 100
Abbildung 7: Schnitt auf Layer 8: Entwicklung der Windgeschwindigkeit mit der Zeit.
(Abbildung 7).
In den Modellen wird die Windgeschwindigkeit oberhalb der Scherzone gleich 0 m/s gesetzt und unterhalb gleich 20 m/s (Abbildung 8 und 9). Zum Vergleich ist in Abbildung 10 das selbe Modell ohne Scherzone dargestellt.
A Quellcode
Die ben¨ otigen Programme sind in diesem PDF angeh¨ angt.
Fortran (Modell)
Matlab (Matlab-Skripte zur Visualisierung) isentrop.f
matlab.zip
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 8: Scherzone von 5 km bis 7 km, Bergh¨ ohe 1.4 km, Geschwindig-
keit oberhalb 0 m/s unterhalb 20 m/s zum Zeitpunkt 0. Auf der Leeseite des
Berges bildet sich ein Abw¨ arts Windsturm mit sehr hohen Windgeschwin-
digkeiten aus. Die Gravitationswelle wird an der Scherzone geblockt und es
kommt zu Resonanzen mit hohen Windgeschwindigkeiten
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]
Abbildung 9: Scherzone von 9 km bis 11 km, Bergh¨ ohe 1.4 km, Geschwindig-
keit oberhalb 0 m/s unterhalb 20 m/s zum Zeitpunkt 0. Die vertikale Gravi-
tationswelle wird an der Scherzone abgeblockt.
0 50 100 150 200 250 300 350 400 450 500 0
1 2 3 4 5 6 7 8 9 10
x [km]
z [km]