Home

Context Navigation

source: palm/trunk/SOURCE/prognostic_equations.f90 @ 94

Last change on this file since 94 was 94, checked in by raasch, 17 years ago
preliminary uncomplete changes for ocean version
Property svn:keywords set to `Id`
File size: 54.8 KB

Line
1	MODULE prognostic_equations_mod
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! new argument zw in calls of diffusion_e
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: prognostic_equations.f90 94 2007-06-01 15:25:22Z raasch $
11	!
12	! 75 2007-03-22 09:54:05Z raasch
13	! checking for negative q and limiting for positive values,
14	! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
15	! subroutine names changed to .._noopt, .._cache, and .._vector,
16	! moisture renamed humidity, Bott-Chlond-scheme can be used in the
17	! _vector-version
18	!
19	! 19 2007-02-23 04:53:48Z raasch
20	! Calculation of e, q, and pt extended for gridpoint nzt,
21	! handling of given temperature/humidity/scalar fluxes at top surface
22	!
23	! RCS Log replace by Id keyword, revision history cleaned up
24	!
25	! Revision 1.21 2006/08/04 15:01:07 raasch
26	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
27	! regardless of the timestep scheme used for the other quantities,
28	! new argument diss in call of diffusion_e
29	!
30	! Revision 1.1 2000/04/13 14:56:27 schroeter
31	! Initial revision
32	!
33	!
34	! Description:
35	! ------------
36	! Solving the prognostic equations.
37	!------------------------------------------------------------------------------!
38
39	USE arrays_3d
40	USE control_parameters
41	USE cpulog
42	USE grid_variables
43	USE indices
44	USE interfaces
45	USE pegrid
46	USE pointer_interfaces
47	USE statistics
48
49	USE advec_s_pw_mod
50	USE advec_s_up_mod
51	USE advec_u_pw_mod
52	USE advec_u_up_mod
53	USE advec_v_pw_mod
54	USE advec_v_up_mod
55	USE advec_w_pw_mod
56	USE advec_w_up_mod
57	USE buoyancy_mod
58	USE calc_precipitation_mod
59	USE calc_radiation_mod
60	USE coriolis_mod
61	USE diffusion_e_mod
62	USE diffusion_s_mod
63	USE diffusion_u_mod
64	USE diffusion_v_mod
65	USE diffusion_w_mod
66	USE impact_of_latent_heat_mod
67	USE production_e_mod
68	USE user_actions_mod
69
70
71	PRIVATE
72	PUBLIC prognostic_equations_noopt, prognostic_equations_cache, &
73	prognostic_equations_vector
74
75	INTERFACE prognostic_equations_noopt
76	MODULE PROCEDURE prognostic_equations_noopt
77	END INTERFACE prognostic_equations_noopt
78
79	INTERFACE prognostic_equations_cache
80	MODULE PROCEDURE prognostic_equations_cache
81	END INTERFACE prognostic_equations_cache
82
83	INTERFACE prognostic_equations_vector
84	MODULE PROCEDURE prognostic_equations_vector
85	END INTERFACE prognostic_equations_vector
86
87
88	CONTAINS
89
90
91	SUBROUTINE prognostic_equations_noopt
92
93	!------------------------------------------------------------------------------!
94	! Version with single loop optimization
95	!
96	! (Optimized over each single prognostic equation.)
97	!------------------------------------------------------------------------------!
98
99	IMPLICIT NONE
100
101	CHARACTER (LEN=9) :: time_to_string
102	INTEGER :: i, j, k
103	REAL :: sat, sbt
104
105	!
106	!-- Calculate those variables needed in the tendency terms which need
107	!-- global communication
108	CALL calc_mean_pt_profile( pt, 4 )
109	IF ( humidity ) CALL calc_mean_pt_profile( vpt, 44 )
110
111	!
112	!-- u-velocity component
113	CALL cpu_log( log_point(5), 'u-equation', 'start' )
114
115	!
116	!-- u-tendency terms with communication
117	IF ( momentum_advec == 'ups-scheme' ) THEN
118	tend = 0.0
119	CALL advec_u_ups
120	ENDIF
121
122	!
123	!-- u-tendency terms with no communication
124	DO i = nxl, nxr
125	DO j = nys, nyn
126	!
127	!-- Tendency terms
128	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
129	tend(:,j,i) = 0.0
130	CALL advec_u_pw( i, j )
131	ELSE
132	IF ( momentum_advec /= 'ups-scheme' ) THEN
133	tend(:,j,i) = 0.0
134	CALL advec_u_up( i, j )
135	ENDIF
136	ENDIF
137	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
138	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m, &
139	usws_m, v_m, w_m )
140	ELSE
141	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
142	v, w )
143	ENDIF
144	CALL coriolis( i, j, 1 )
145	IF ( sloping_surface ) CALL buoyancy( i, j, pt, 1, 4 )
146	CALL user_actions( i, j, 'u-tendency' )
147
148	!
149	!-- Prognostic equation for u-velocity component
150	DO k = nzb_u_inner(j,i)+1, nzt
151	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
152	dt_3d * ( &
153	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
154	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
155	) - &
156	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
157	ENDDO
158
159	!
160	!-- Calculate tendencies for the next Runge-Kutta step
161	IF ( timestep_scheme(1:5) == 'runge' ) THEN
162	IF ( intermediate_timestep_count == 1 ) THEN
163	DO k = nzb_u_inner(j,i)+1, nzt
164	tu_m(k,j,i) = tend(k,j,i)
165	ENDDO
166	ELSEIF ( intermediate_timestep_count < &
167	intermediate_timestep_count_max ) THEN
168	DO k = nzb_u_inner(j,i)+1, nzt
169	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
170	ENDDO
171	ENDIF
172	ENDIF
173
174	ENDDO
175	ENDDO
176
177	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
178
179	!
180	!-- v-velocity component
181	CALL cpu_log( log_point(6), 'v-equation', 'start' )
182
183	!
184	!-- v-tendency terms with communication
185	IF ( momentum_advec == 'ups-scheme' ) THEN
186	tend = 0.0
187	CALL advec_v_ups
188	ENDIF
189
190	!
191	!-- v-tendency terms with no communication
192	DO i = nxl, nxr
193	DO j = nys, nyn
194	!
195	!-- Tendency terms
196	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
197	tend(:,j,i) = 0.0
198	CALL advec_v_pw( i, j )
199	ELSE
200	IF ( momentum_advec /= 'ups-scheme' ) THEN
201	tend(:,j,i) = 0.0
202	CALL advec_v_up( i, j )
203	ENDIF
204	ENDIF
205	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
206	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
207	v_m, vsws_m, w_m )
208	ELSE
209	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
210	vsws, w )
211	ENDIF
212	CALL coriolis( i, j, 2 )
213	CALL user_actions( i, j, 'v-tendency' )
214
215	!
216	!-- Prognostic equation for v-velocity component
217	DO k = nzb_v_inner(j,i)+1, nzt
218	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
219	dt_3d * ( &
220	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
221	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
222	) - &
223	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
224	ENDDO
225
226	!
227	!-- Calculate tendencies for the next Runge-Kutta step
228	IF ( timestep_scheme(1:5) == 'runge' ) THEN
229	IF ( intermediate_timestep_count == 1 ) THEN
230	DO k = nzb_v_inner(j,i)+1, nzt
231	tv_m(k,j,i) = tend(k,j,i)
232	ENDDO
233	ELSEIF ( intermediate_timestep_count < &
234	intermediate_timestep_count_max ) THEN
235	DO k = nzb_v_inner(j,i)+1, nzt
236	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
237	ENDDO
238	ENDIF
239	ENDIF
240
241	ENDDO
242	ENDDO
243
244	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
245
246	!
247	!-- w-velocity component
248	CALL cpu_log( log_point(7), 'w-equation', 'start' )
249
250	!
251	!-- w-tendency terms with communication
252	IF ( momentum_advec == 'ups-scheme' ) THEN
253	tend = 0.0
254	CALL advec_w_ups
255	ENDIF
256
257	!
258	!-- w-tendency terms with no communication
259	DO i = nxl, nxr
260	DO j = nys, nyn
261	!
262	!-- Tendency terms
263	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
264	tend(:,j,i) = 0.0
265	CALL advec_w_pw( i, j )
266	ELSE
267	IF ( momentum_advec /= 'ups-scheme' ) THEN
268	tend(:,j,i) = 0.0
269	CALL advec_w_up( i, j )
270	ENDIF
271	ENDIF
272	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
273	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, km_damp_y, &
274	tend, u_m, v_m, w_m )
275	ELSE
276	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
277	tend, u, v, w )
278	ENDIF
279	CALL coriolis( i, j, 3 )
280	IF ( .NOT. humidity ) THEN
281	CALL buoyancy( i, j, pt, 3, 4 )
282	ELSE
283	CALL buoyancy( i, j, vpt, 3, 44 )
284	ENDIF
285	CALL user_actions( i, j, 'w-tendency' )
286
287	!
288	!-- Prognostic equation for w-velocity component
289	DO k = nzb_w_inner(j,i)+1, nzt-1
290	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
291	dt_3d * ( &
292	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
293	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
294	) - &
295	tsc(5) * rdf(k) * w(k,j,i)
296	ENDDO
297
298	!
299	!-- Calculate tendencies for the next Runge-Kutta step
300	IF ( timestep_scheme(1:5) == 'runge' ) THEN
301	IF ( intermediate_timestep_count == 1 ) THEN
302	DO k = nzb_w_inner(j,i)+1, nzt-1
303	tw_m(k,j,i) = tend(k,j,i)
304	ENDDO
305	ELSEIF ( intermediate_timestep_count < &
306	intermediate_timestep_count_max ) THEN
307	DO k = nzb_w_inner(j,i)+1, nzt-1
308	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
309	ENDDO
310	ENDIF
311	ENDIF
312
313	ENDDO
314	ENDDO
315
316	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
317
318	!
319	!-- potential temperature
320	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
321
322	!
323	!-- pt-tendency terms with communication
324	sat = tsc(1)
325	sbt = tsc(2)
326	IF ( scalar_advec == 'bc-scheme' ) THEN
327
328	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
329	!
330	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
331	!-- switched on. Thus:
332	sat = 1.0
333	sbt = 1.0
334	ENDIF
335	tend = 0.0
336	CALL advec_s_bc( pt, 'pt' )
337	ELSE
338	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
339	tend = 0.0
340	CALL advec_s_ups( pt, 'pt' )
341	ENDIF
342	ENDIF
343
344	!
345	!-- pt-tendency terms with no communication
346	DO i = nxl, nxr
347	DO j = nys, nyn
348	!
349	!-- Tendency terms
350	IF ( scalar_advec == 'bc-scheme' ) THEN
351	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
352	ELSE
353	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
354	tend(:,j,i) = 0.0
355	CALL advec_s_pw( i, j, pt )
356	ELSE
357	IF ( scalar_advec /= 'ups-scheme' ) THEN
358	tend(:,j,i) = 0.0
359	CALL advec_s_up( i, j, pt )
360	ENDIF
361	ENDIF
362	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
363	THEN
364	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
365	tswst_m, tend )
366	ELSE
367	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
368	ENDIF
369	ENDIF
370
371	!
372	!-- If required compute heating/cooling due to long wave radiation
373	!-- processes
374	IF ( radiation ) THEN
375	CALL calc_radiation( i, j )
376	ENDIF
377
378	!
379	!-- If required compute impact of latent heat due to precipitation
380	IF ( precipitation ) THEN
381	CALL impact_of_latent_heat( i, j )
382	ENDIF
383	CALL user_actions( i, j, 'pt-tendency' )
384
385	!
386	!-- Prognostic equation for potential temperature
387	DO k = nzb_s_inner(j,i)+1, nzt
388	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
389	dt_3d * ( &
390	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
391	) - &
392	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
393	ENDDO
394
395	!
396	!-- Calculate tendencies for the next Runge-Kutta step
397	IF ( timestep_scheme(1:5) == 'runge' ) THEN
398	IF ( intermediate_timestep_count == 1 ) THEN
399	DO k = nzb_s_inner(j,i)+1, nzt
400	tpt_m(k,j,i) = tend(k,j,i)
401	ENDDO
402	ELSEIF ( intermediate_timestep_count < &
403	intermediate_timestep_count_max ) THEN
404	DO k = nzb_s_inner(j,i)+1, nzt
405	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
406	ENDDO
407	ENDIF
408	ENDIF
409
410	ENDDO
411	ENDDO
412
413	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
414
415	!
416	!-- If required, compute prognostic equation for total water content / scalar
417	IF ( humidity .OR. passive_scalar ) THEN
418
419	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
420
421	!
422	!-- Scalar/q-tendency terms with communication
423	sat = tsc(1)
424	sbt = tsc(2)
425	IF ( scalar_advec == 'bc-scheme' ) THEN
426
427	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
428	!
429	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
430	!-- switched on. Thus:
431	sat = 1.0
432	sbt = 1.0
433	ENDIF
434	tend = 0.0
435	CALL advec_s_bc( q, 'q' )
436	ELSE
437	IF ( tsc(2) /= 2.0 ) THEN
438	IF ( scalar_advec == 'ups-scheme' ) THEN
439	tend = 0.0
440	CALL advec_s_ups( q, 'q' )
441	ENDIF
442	ENDIF
443	ENDIF
444
445	!
446	!-- Scalar/q-tendency terms with no communication
447	DO i = nxl, nxr
448	DO j = nys, nyn
449	!
450	!-- Tendency-terms
451	IF ( scalar_advec == 'bc-scheme' ) THEN
452	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
453	ELSE
454	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
455	tend(:,j,i) = 0.0
456	CALL advec_s_pw( i, j, q )
457	ELSE
458	IF ( scalar_advec /= 'ups-scheme' ) THEN
459	tend(:,j,i) = 0.0
460	CALL advec_s_up( i, j, q )
461	ENDIF
462	ENDIF
463	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
464	THEN
465	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
466	qswst_m, tend )
467	ELSE
468	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
469	tend )
470	ENDIF
471	ENDIF
472
473	!
474	!-- If required compute decrease of total water content due to
475	!-- precipitation
476	IF ( precipitation ) THEN
477	CALL calc_precipitation( i, j )
478	ENDIF
479	CALL user_actions( i, j, 'q-tendency' )
480
481	!
482	!-- Prognostic equation for total water content / scalar
483	DO k = nzb_s_inner(j,i)+1, nzt
484	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
485	dt_3d * ( &
486	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
487	) - &
488	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
489	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
490	ENDDO
491
492	!
493	!-- Calculate tendencies for the next Runge-Kutta step
494	IF ( timestep_scheme(1:5) == 'runge' ) THEN
495	IF ( intermediate_timestep_count == 1 ) THEN
496	DO k = nzb_s_inner(j,i)+1, nzt
497	tq_m(k,j,i) = tend(k,j,i)
498	ENDDO
499	ELSEIF ( intermediate_timestep_count < &
500	intermediate_timestep_count_max ) THEN
501	DO k = nzb_s_inner(j,i)+1, nzt
502	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
503	ENDDO
504	ENDIF
505	ENDIF
506
507	ENDDO
508	ENDDO
509
510	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
511
512	ENDIF
513
514	!
515	!-- If required, compute prognostic equation for turbulent kinetic
516	!-- energy (TKE)
517	IF ( .NOT. constant_diffusion ) THEN
518
519	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
520
521	!
522	!-- TKE-tendency terms with communication
523	CALL production_e_init
524
525	sat = tsc(1)
526	sbt = tsc(2)
527	IF ( .NOT. use_upstream_for_tke ) THEN
528	IF ( scalar_advec == 'bc-scheme' ) THEN
529
530	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
531	!
532	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
533	!-- switched on. Thus:
534	sat = 1.0
535	sbt = 1.0
536	ENDIF
537	tend = 0.0
538	CALL advec_s_bc( e, 'e' )
539	ELSE
540	IF ( tsc(2) /= 2.0 ) THEN
541	IF ( scalar_advec == 'ups-scheme' ) THEN
542	tend = 0.0
543	CALL advec_s_ups( e, 'e' )
544	ENDIF
545	ENDIF
546	ENDIF
547	ENDIF
548
549	!
550	!-- TKE-tendency terms with no communication
551	DO i = nxl, nxr
552	DO j = nys, nyn
553	!
554	!-- Tendency-terms
555	IF ( scalar_advec == 'bc-scheme' .AND. &
556	.NOT. use_upstream_for_tke ) THEN
557	IF ( .NOT. humidity ) THEN
558	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
559	l_grid, pt, rif, tend, zu, zw )
560	ELSE
561	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
562	l_grid, vpt, rif, tend, zu, zw )
563	ENDIF
564	ELSE
565	IF ( use_upstream_for_tke ) THEN
566	tend(:,j,i) = 0.0
567	CALL advec_s_up( i, j, e )
568	ELSE
569	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
570	THEN
571	tend(:,j,i) = 0.0
572	CALL advec_s_pw( i, j, e )
573	ELSE
574	IF ( scalar_advec /= 'ups-scheme' ) THEN
575	tend(:,j,i) = 0.0
576	CALL advec_s_up( i, j, e )
577	ENDIF
578	ENDIF
579	ENDIF
580	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
581	THEN
582	IF ( .NOT. humidity ) THEN
583	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
584	km_m, l_grid, pt_m, rif_m, tend, zu, &
585	zw )
586	ELSE
587	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
588	km_m, l_grid, vpt_m, rif_m, tend, zu, &
589	zw )
590	ENDIF
591	ELSE
592	IF ( .NOT. humidity ) THEN
593	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
594	l_grid, pt, rif, tend, zu, zw )
595	ELSE
596	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
597	l_grid, vpt, rif, tend, zu, zw )
598	ENDIF
599	ENDIF
600	ENDIF
601	CALL production_e( i, j )
602	CALL user_actions( i, j, 'e-tendency' )
603
604	!
605	!-- Prognostic equation for TKE.
606	!-- Eliminate negative TKE values, which can occur due to numerical
607	!-- reasons in the course of the integration. In such cases the old TKE
608	!-- value is reduced by 90%.
609	DO k = nzb_s_inner(j,i)+1, nzt
610	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
611	dt_3d * ( &
612	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
613	)
614	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
615	ENDDO
616
617	!
618	!-- Calculate tendencies for the next Runge-Kutta step
619	IF ( timestep_scheme(1:5) == 'runge' ) THEN
620	IF ( intermediate_timestep_count == 1 ) THEN
621	DO k = nzb_s_inner(j,i)+1, nzt
622	te_m(k,j,i) = tend(k,j,i)
623	ENDDO
624	ELSEIF ( intermediate_timestep_count < &
625	intermediate_timestep_count_max ) THEN
626	DO k = nzb_s_inner(j,i)+1, nzt
627	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
628	ENDDO
629	ENDIF
630	ENDIF
631
632	ENDDO
633	ENDDO
634
635	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
636
637	ENDIF
638
639
640	END SUBROUTINE prognostic_equations_noopt
641
642
643	SUBROUTINE prognostic_equations_cache
644
645	!------------------------------------------------------------------------------!
646	! Version with one optimized loop over all equations. It is only allowed to
647	! be called for the standard Piascek-Williams advection scheme.
648	!
649	! The call of this subroutine is embedded in two DO loops over i and j, thus
650	! communication between CPUs is not allowed in this subroutine.
651	!
652	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
653	!------------------------------------------------------------------------------!
654
655	IMPLICIT NONE
656
657	CHARACTER (LEN=9) :: time_to_string
658	INTEGER :: i, j, k
659
660
661	!
662	!-- Time measurement can only be performed for the whole set of equations
663	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
664
665
666	!
667	!-- Calculate those variables needed in the tendency terms which need
668	!-- global communication
669	CALL calc_mean_pt_profile( pt, 4 )
670	IF ( humidity ) CALL calc_mean_pt_profile( vpt, 44 )
671	IF ( .NOT. constant_diffusion ) CALL production_e_init
672
673
674	!
675	!-- Loop over all prognostic equations
676	!$OMP PARALLEL private (i,j,k)
677	!$OMP DO
678	DO i = nxl, nxr
679	DO j = nys, nyn
680	!
681	!-- Tendency terms for u-velocity component
682	IF ( j < nyn+1 ) THEN
683
684	tend(:,j,i) = 0.0
685	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
686	CALL advec_u_pw( i, j )
687	ELSE
688	CALL advec_u_up( i, j )
689	ENDIF
690	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
691	THEN
692	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
693	u_m, usws_m, v_m, w_m )
694	ELSE
695	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
696	usws, v, w )
697	ENDIF
698	CALL coriolis( i, j, 1 )
699	IF ( sloping_surface ) CALL buoyancy( i, j, pt, 1, 4 )
700	CALL user_actions( i, j, 'u-tendency' )
701
702	!
703	!-- Prognostic equation for u-velocity component
704	DO k = nzb_u_inner(j,i)+1, nzt
705	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
706	dt_3d * ( &
707	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
708	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
709	) - &
710	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
711	ENDDO
712
713	!
714	!-- Calculate tendencies for the next Runge-Kutta step
715	IF ( timestep_scheme(1:5) == 'runge' ) THEN
716	IF ( intermediate_timestep_count == 1 ) THEN
717	DO k = nzb_u_inner(j,i)+1, nzt
718	tu_m(k,j,i) = tend(k,j,i)
719	ENDDO
720	ELSEIF ( intermediate_timestep_count < &
721	intermediate_timestep_count_max ) THEN
722	DO k = nzb_u_inner(j,i)+1, nzt
723	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
724	ENDDO
725	ENDIF
726	ENDIF
727
728	ENDIF
729
730	!
731	!-- Tendency terms for v-velocity component
732	IF ( i < nxr+1 ) THEN
733
734	tend(:,j,i) = 0.0
735	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
736	CALL advec_v_pw( i, j )
737	ELSE
738	CALL advec_v_up( i, j )
739	ENDIF
740	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
741	THEN
742	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
743	u_m, v_m, vsws_m, w_m )
744	ELSE
745	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
746	vsws, w )
747	ENDIF
748	CALL coriolis( i, j, 2 )
749	CALL user_actions( i, j, 'v-tendency' )
750
751	!
752	!-- Prognostic equation for v-velocity component
753	DO k = nzb_v_inner(j,i)+1, nzt
754	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
755	dt_3d * ( &
756	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
757	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
758	) - &
759	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
760	ENDDO
761
762	!
763	!-- Calculate tendencies for the next Runge-Kutta step
764	IF ( timestep_scheme(1:5) == 'runge' ) THEN
765	IF ( intermediate_timestep_count == 1 ) THEN
766	DO k = nzb_v_inner(j,i)+1, nzt
767	tv_m(k,j,i) = tend(k,j,i)
768	ENDDO
769	ELSEIF ( intermediate_timestep_count < &
770	intermediate_timestep_count_max ) THEN
771	DO k = nzb_v_inner(j,i)+1, nzt
772	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
773	ENDDO
774	ENDIF
775	ENDIF
776
777	ENDIF
778
779	!
780	!-- Tendency terms for w-velocity component
781	IF ( i < nxr+1 .AND. j < nyn+1 ) THEN
782
783	tend(:,j,i) = 0.0
784	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
785	CALL advec_w_pw( i, j )
786	ELSE
787	CALL advec_w_up( i, j )
788	ENDIF
789	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
790	THEN
791	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
792	km_damp_y, tend, u_m, v_m, w_m )
793	ELSE
794	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
795	tend, u, v, w )
796	ENDIF
797	CALL coriolis( i, j, 3 )
798	IF ( .NOT. humidity ) THEN
799	CALL buoyancy( i, j, pt, 3, 4 )
800	ELSE
801	CALL buoyancy( i, j, vpt, 3, 44 )
802	ENDIF
803	CALL user_actions( i, j, 'w-tendency' )
804
805	!
806	!-- Prognostic equation for w-velocity component
807	DO k = nzb_w_inner(j,i)+1, nzt-1
808	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
809	dt_3d * ( &
810	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
811	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
812	) - &
813	tsc(5) * rdf(k) * w(k,j,i)
814	ENDDO
815
816	!
817	!-- Calculate tendencies for the next Runge-Kutta step
818	IF ( timestep_scheme(1:5) == 'runge' ) THEN
819	IF ( intermediate_timestep_count == 1 ) THEN
820	DO k = nzb_w_inner(j,i)+1, nzt-1
821	tw_m(k,j,i) = tend(k,j,i)
822	ENDDO
823	ELSEIF ( intermediate_timestep_count < &
824	intermediate_timestep_count_max ) THEN
825	DO k = nzb_w_inner(j,i)+1, nzt-1
826	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
827	ENDDO
828	ENDIF
829	ENDIF
830
831	!
832	!-- Tendency terms for potential temperature
833	tend(:,j,i) = 0.0
834	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
835	CALL advec_s_pw( i, j, pt )
836	ELSE
837	CALL advec_s_up( i, j, pt )
838	ENDIF
839	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
840	THEN
841	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
842	tswst_m, tend )
843	ELSE
844	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
845	ENDIF
846
847	!
848	!-- If required compute heating/cooling due to long wave radiation
849	!-- processes
850	IF ( radiation ) THEN
851	CALL calc_radiation( i, j )
852	ENDIF
853
854	!
855	!-- If required compute impact of latent heat due to precipitation
856	IF ( precipitation ) THEN
857	CALL impact_of_latent_heat( i, j )
858	ENDIF
859	CALL user_actions( i, j, 'pt-tendency' )
860
861	!
862	!-- Prognostic equation for potential temperature
863	DO k = nzb_s_inner(j,i)+1, nzt
864	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
865	dt_3d * ( &
866	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
867	) - &
868	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
869	ENDDO
870
871	!
872	!-- Calculate tendencies for the next Runge-Kutta step
873	IF ( timestep_scheme(1:5) == 'runge' ) THEN
874	IF ( intermediate_timestep_count == 1 ) THEN
875	DO k = nzb_s_inner(j,i)+1, nzt
876	tpt_m(k,j,i) = tend(k,j,i)
877	ENDDO
878	ELSEIF ( intermediate_timestep_count < &
879	intermediate_timestep_count_max ) THEN
880	DO k = nzb_s_inner(j,i)+1, nzt
881	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
882	5.3125 * tpt_m(k,j,i)
883	ENDDO
884	ENDIF
885	ENDIF
886
887	!
888	!-- If required, compute prognostic equation for total water content /
889	!-- scalar
890	IF ( humidity .OR. passive_scalar ) THEN
891
892	!
893	!-- Tendency-terms for total water content / scalar
894	tend(:,j,i) = 0.0
895	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
896	THEN
897	CALL advec_s_pw( i, j, q )
898	ELSE
899	CALL advec_s_up( i, j, q )
900	ENDIF
901	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
902	THEN
903	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
904	qswst_m, tend )
905	ELSE
906	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
907	tend )
908	ENDIF
909
910	!
911	!-- If required compute decrease of total water content due to
912	!-- precipitation
913	IF ( precipitation ) THEN
914	CALL calc_precipitation( i, j )
915	ENDIF
916	CALL user_actions( i, j, 'q-tendency' )
917
918	!
919	!-- Prognostic equation for total water content / scalar
920	DO k = nzb_s_inner(j,i)+1, nzt
921	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
922	dt_3d * ( &
923	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
924	) - &
925	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
926	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
927	ENDDO
928
929	!
930	!-- Calculate tendencies for the next Runge-Kutta step
931	IF ( timestep_scheme(1:5) == 'runge' ) THEN
932	IF ( intermediate_timestep_count == 1 ) THEN
933	DO k = nzb_s_inner(j,i)+1, nzt
934	tq_m(k,j,i) = tend(k,j,i)
935	ENDDO
936	ELSEIF ( intermediate_timestep_count < &
937	intermediate_timestep_count_max ) THEN
938	DO k = nzb_s_inner(j,i)+1, nzt
939	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
940	5.3125 * tq_m(k,j,i)
941	ENDDO
942	ENDIF
943	ENDIF
944
945	ENDIF
946
947	!
948	!-- If required, compute prognostic equation for turbulent kinetic
949	!-- energy (TKE)
950	IF ( .NOT. constant_diffusion ) THEN
951
952	!
953	!-- Tendency-terms for TKE
954	tend(:,j,i) = 0.0
955	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
956	.AND. .NOT. use_upstream_for_tke ) THEN
957	CALL advec_s_pw( i, j, e )
958	ELSE
959	CALL advec_s_up( i, j, e )
960	ENDIF
961	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
962	THEN
963	IF ( .NOT. humidity ) THEN
964	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
965	km_m, l_grid, pt_m, rif_m, tend, zu, &
966	zw )
967	ELSE
968	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
969	km_m, l_grid, vpt_m, rif_m, tend, zu, &
970	zw )
971	ENDIF
972	ELSE
973	IF ( .NOT. humidity ) THEN
974	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
975	l_grid, pt, rif, tend, zu, zw )
976	ELSE
977	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
978	l_grid, vpt, rif, tend, zu, zw )
979	ENDIF
980	ENDIF
981	CALL production_e( i, j )
982	CALL user_actions( i, j, 'e-tendency' )
983
984	!
985	!-- Prognostic equation for TKE.
986	!-- Eliminate negative TKE values, which can occur due to numerical
987	!-- reasons in the course of the integration. In such cases the old
988	!-- TKE value is reduced by 90%.
989	DO k = nzb_s_inner(j,i)+1, nzt
990	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
991	dt_3d * ( &
992	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
993	)
994	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
995	ENDDO
996
997	!
998	!-- Calculate tendencies for the next Runge-Kutta step
999	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1000	IF ( intermediate_timestep_count == 1 ) THEN
1001	DO k = nzb_s_inner(j,i)+1, nzt
1002	te_m(k,j,i) = tend(k,j,i)
1003	ENDDO
1004	ELSEIF ( intermediate_timestep_count < &
1005	intermediate_timestep_count_max ) THEN
1006	DO k = nzb_s_inner(j,i)+1, nzt
1007	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1008	5.3125 * te_m(k,j,i)
1009	ENDDO
1010	ENDIF
1011	ENDIF
1012
1013	ENDIF ! TKE equation
1014
1015	ENDIF ! Gridpoints excluding the non-cyclic wall
1016
1017	ENDDO
1018	ENDDO
1019	!$OMP END PARALLEL
1020
1021	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1022
1023
1024	END SUBROUTINE prognostic_equations_cache
1025
1026
1027	SUBROUTINE prognostic_equations_vector
1028
1029	!------------------------------------------------------------------------------!
1030	! Version for vector machines
1031	!------------------------------------------------------------------------------!
1032
1033	IMPLICIT NONE
1034
1035	CHARACTER (LEN=9) :: time_to_string
1036	INTEGER :: i, j, k
1037	REAL :: sat, sbt
1038
1039	!
1040	!-- Calculate those variables needed in the tendency terms which need
1041	!-- global communication
1042	CALL calc_mean_pt_profile( pt, 4 )
1043	IF ( humidity ) CALL calc_mean_pt_profile( vpt, 44 )
1044
1045	!
1046	!-- u-velocity component
1047	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1048
1049	!
1050	!-- u-tendency terms with communication
1051	IF ( momentum_advec == 'ups-scheme' ) THEN
1052	tend = 0.0
1053	CALL advec_u_ups
1054	ENDIF
1055
1056	!
1057	!-- u-tendency terms with no communication
1058	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1059	tend = 0.0
1060	CALL advec_u_pw
1061	ELSE
1062	IF ( momentum_advec /= 'ups-scheme' ) THEN
1063	tend = 0.0
1064	CALL advec_u_up
1065	ENDIF
1066	ENDIF
1067	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1068	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
1069	w_m )
1070	ELSE
1071	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
1072	ENDIF
1073	CALL coriolis( 1 )
1074	IF ( sloping_surface ) CALL buoyancy( pt, 1, 4 )
1075	CALL user_actions( 'u-tendency' )
1076
1077	!
1078	!-- Prognostic equation for u-velocity component
1079	DO i = nxl, nxr
1080	DO j = nys, nyn
1081	DO k = nzb_u_inner(j,i)+1, nzt
1082	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1083	dt_3d * ( &
1084	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1085	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
1086	) - &
1087	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1088	ENDDO
1089	ENDDO
1090	ENDDO
1091
1092	!
1093	!-- Calculate tendencies for the next Runge-Kutta step
1094	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1095	IF ( intermediate_timestep_count == 1 ) THEN
1096	DO i = nxl, nxr
1097	DO j = nys, nyn
1098	DO k = nzb_u_inner(j,i)+1, nzt
1099	tu_m(k,j,i) = tend(k,j,i)
1100	ENDDO
1101	ENDDO
1102	ENDDO
1103	ELSEIF ( intermediate_timestep_count < &
1104	intermediate_timestep_count_max ) THEN
1105	DO i = nxl, nxr
1106	DO j = nys, nyn
1107	DO k = nzb_u_inner(j,i)+1, nzt
1108	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1109	ENDDO
1110	ENDDO
1111	ENDDO
1112	ENDIF
1113	ENDIF
1114
1115	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1116
1117	!
1118	!-- v-velocity component
1119	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1120
1121	!
1122	!-- v-tendency terms with communication
1123	IF ( momentum_advec == 'ups-scheme' ) THEN
1124	tend = 0.0
1125	CALL advec_v_ups
1126	ENDIF
1127
1128	!
1129	!-- v-tendency terms with no communication
1130	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1131	tend = 0.0
1132	CALL advec_v_pw
1133	ELSE
1134	IF ( momentum_advec /= 'ups-scheme' ) THEN
1135	tend = 0.0
1136	CALL advec_v_up
1137	ENDIF
1138	ENDIF
1139	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1140	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1141	w_m )
1142	ELSE
1143	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
1144	ENDIF
1145	CALL coriolis( 2 )
1146	CALL user_actions( 'v-tendency' )
1147
1148	!
1149	!-- Prognostic equation for v-velocity component
1150	DO i = nxl, nxr
1151	DO j = nys, nyn
1152	DO k = nzb_v_inner(j,i)+1, nzt
1153	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1154	dt_3d * ( &
1155	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1156	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
1157	) - &
1158	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1159	ENDDO
1160	ENDDO
1161	ENDDO
1162
1163	!
1164	!-- Calculate tendencies for the next Runge-Kutta step
1165	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1166	IF ( intermediate_timestep_count == 1 ) THEN
1167	DO i = nxl, nxr
1168	DO j = nys, nyn
1169	DO k = nzb_v_inner(j,i)+1, nzt
1170	tv_m(k,j,i) = tend(k,j,i)
1171	ENDDO
1172	ENDDO
1173	ENDDO
1174	ELSEIF ( intermediate_timestep_count < &
1175	intermediate_timestep_count_max ) THEN
1176	DO i = nxl, nxr
1177	DO j = nys, nyn
1178	DO k = nzb_v_inner(j,i)+1, nzt
1179	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1180	ENDDO
1181	ENDDO
1182	ENDDO
1183	ENDIF
1184	ENDIF
1185
1186	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1187
1188	!
1189	!-- w-velocity component
1190	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1191
1192	!
1193	!-- w-tendency terms with communication
1194	IF ( momentum_advec == 'ups-scheme' ) THEN
1195	tend = 0.0
1196	CALL advec_w_ups
1197	ENDIF
1198
1199	!
1200	!-- w-tendency terms with no communication
1201	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1202	tend = 0.0
1203	CALL advec_w_pw
1204	ELSE
1205	IF ( momentum_advec /= 'ups-scheme' ) THEN
1206	tend = 0.0
1207	CALL advec_w_up
1208	ENDIF
1209	ENDIF
1210	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1211	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1212	v_m, w_m )
1213	ELSE
1214	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1215	ENDIF
1216	CALL coriolis( 3 )
1217	IF ( .NOT. humidity ) THEN
1218	CALL buoyancy( pt, 3, 4 )
1219	ELSE
1220	CALL buoyancy( vpt, 3, 44 )
1221	ENDIF
1222	CALL user_actions( 'w-tendency' )
1223
1224	!
1225	!-- Prognostic equation for w-velocity component
1226	DO i = nxl, nxr
1227	DO j = nys, nyn
1228	DO k = nzb_w_inner(j,i)+1, nzt-1
1229	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1230	dt_3d * ( &
1231	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1232	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1233	) - &
1234	tsc(5) * rdf(k) * w(k,j,i)
1235	ENDDO
1236	ENDDO
1237	ENDDO
1238
1239	!
1240	!-- Calculate tendencies for the next Runge-Kutta step
1241	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1242	IF ( intermediate_timestep_count == 1 ) THEN
1243	DO i = nxl, nxr
1244	DO j = nys, nyn
1245	DO k = nzb_w_inner(j,i)+1, nzt-1
1246	tw_m(k,j,i) = tend(k,j,i)
1247	ENDDO
1248	ENDDO
1249	ENDDO
1250	ELSEIF ( intermediate_timestep_count < &
1251	intermediate_timestep_count_max ) THEN
1252	DO i = nxl, nxr
1253	DO j = nys, nyn
1254	DO k = nzb_w_inner(j,i)+1, nzt-1
1255	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1256	ENDDO
1257	ENDDO
1258	ENDDO
1259	ENDIF
1260	ENDIF
1261
1262	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1263
1264	!
1265	!-- potential temperature
1266	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1267
1268	!
1269	!-- pt-tendency terms with communication
1270	sat = tsc(1)
1271	sbt = tsc(2)
1272	IF ( scalar_advec == 'bc-scheme' ) THEN
1273
1274	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1275	!
1276	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1277	!-- switched on. Thus:
1278	sat = 1.0
1279	sbt = 1.0
1280	ENDIF
1281	tend = 0.0
1282	CALL advec_s_bc( pt, 'pt' )
1283	ELSE
1284	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
1285	tend = 0.0
1286	CALL advec_s_ups( pt, 'pt' )
1287	ENDIF
1288	ENDIF
1289
1290	!
1291	!-- pt-tendency terms with no communication
1292	IF ( scalar_advec == 'bc-scheme' ) THEN
1293	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1294	ELSE
1295	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1296	tend = 0.0
1297	CALL advec_s_pw( pt )
1298	ELSE
1299	IF ( scalar_advec /= 'ups-scheme' ) THEN
1300	tend = 0.0
1301	CALL advec_s_up( pt )
1302	ENDIF
1303	ENDIF
1304	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1305	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
1306	ELSE
1307	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1308	ENDIF
1309	ENDIF
1310
1311	!
1312	!-- If required compute heating/cooling due to long wave radiation
1313	!-- processes
1314	IF ( radiation ) THEN
1315	CALL calc_radiation
1316	ENDIF
1317
1318	!
1319	!-- If required compute impact of latent heat due to precipitation
1320	IF ( precipitation ) THEN
1321	CALL impact_of_latent_heat
1322	ENDIF
1323	CALL user_actions( 'pt-tendency' )
1324
1325	!
1326	!-- Prognostic equation for potential temperature
1327	DO i = nxl, nxr
1328	DO j = nys, nyn
1329	DO k = nzb_s_inner(j,i)+1, nzt
1330	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
1331	dt_3d * ( &
1332	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1333	) - &
1334	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1335	ENDDO
1336	ENDDO
1337	ENDDO
1338
1339	!
1340	!-- Calculate tendencies for the next Runge-Kutta step
1341	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1342	IF ( intermediate_timestep_count == 1 ) THEN
1343	DO i = nxl, nxr
1344	DO j = nys, nyn
1345	DO k = nzb_s_inner(j,i)+1, nzt
1346	tpt_m(k,j,i) = tend(k,j,i)
1347	ENDDO
1348	ENDDO
1349	ENDDO
1350	ELSEIF ( intermediate_timestep_count < &
1351	intermediate_timestep_count_max ) THEN
1352	DO i = nxl, nxr
1353	DO j = nys, nyn
1354	DO k = nzb_s_inner(j,i)+1, nzt
1355	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1356	ENDDO
1357	ENDDO
1358	ENDDO
1359	ENDIF
1360	ENDIF
1361
1362	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1363
1364	!
1365	!-- If required, compute prognostic equation for total water content / scalar
1366	IF ( humidity .OR. passive_scalar ) THEN
1367
1368	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1369
1370	!
1371	!-- Scalar/q-tendency terms with communication
1372	sat = tsc(1)
1373	sbt = tsc(2)
1374	IF ( scalar_advec == 'bc-scheme' ) THEN
1375
1376	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1377	!
1378	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1379	!-- switched on. Thus:
1380	sat = 1.0
1381	sbt = 1.0
1382	ENDIF
1383	tend = 0.0
1384	CALL advec_s_bc( q, 'q' )
1385	ELSE
1386	IF ( tsc(2) /= 2.0 ) THEN
1387	IF ( scalar_advec == 'ups-scheme' ) THEN
1388	tend = 0.0
1389	CALL advec_s_ups( q, 'q' )
1390	ENDIF
1391	ENDIF
1392	ENDIF
1393
1394	!
1395	!-- Scalar/q-tendency terms with no communication
1396	IF ( scalar_advec == 'bc-scheme' ) THEN
1397	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1398	ELSE
1399	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1400	tend = 0.0
1401	CALL advec_s_pw( q )
1402	ELSE
1403	IF ( scalar_advec /= 'ups-scheme' ) THEN
1404	tend = 0.0
1405	CALL advec_s_up( q )
1406	ENDIF
1407	ENDIF
1408	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1409	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
1410	ELSE
1411	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1412	ENDIF
1413	ENDIF
1414
1415	!
1416	!-- If required compute decrease of total water content due to
1417	!-- precipitation
1418	IF ( precipitation ) THEN
1419	CALL calc_precipitation
1420	ENDIF
1421	CALL user_actions( 'q-tendency' )
1422
1423	!
1424	!-- Prognostic equation for total water content / scalar
1425	DO i = nxl, nxr
1426	DO j = nys, nyn
1427	DO k = nzb_s_inner(j,i)+1, nzt
1428	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
1429	dt_3d * ( &
1430	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1431	) - &
1432	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1433	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1434	ENDDO
1435	ENDDO
1436	ENDDO
1437
1438	!
1439	!-- Calculate tendencies for the next Runge-Kutta step
1440	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1441	IF ( intermediate_timestep_count == 1 ) THEN
1442	DO i = nxl, nxr
1443	DO j = nys, nyn
1444	DO k = nzb_s_inner(j,i)+1, nzt
1445	tq_m(k,j,i) = tend(k,j,i)
1446	ENDDO
1447	ENDDO
1448	ENDDO
1449	ELSEIF ( intermediate_timestep_count < &
1450	intermediate_timestep_count_max ) THEN
1451	DO i = nxl, nxr
1452	DO j = nys, nyn
1453	DO k = nzb_s_inner(j,i)+1, nzt
1454	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1455	ENDDO
1456	ENDDO
1457	ENDDO
1458	ENDIF
1459	ENDIF
1460
1461	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1462
1463	ENDIF
1464
1465	!
1466	!-- If required, compute prognostic equation for turbulent kinetic
1467	!-- energy (TKE)
1468	IF ( .NOT. constant_diffusion ) THEN
1469
1470	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1471
1472	!
1473	!-- TKE-tendency terms with communication
1474	CALL production_e_init
1475
1476	sat = tsc(1)
1477	sbt = tsc(2)
1478	IF ( .NOT. use_upstream_for_tke ) THEN
1479	IF ( scalar_advec == 'bc-scheme' ) THEN
1480
1481	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1482	!
1483	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1484	!-- switched on. Thus:
1485	sat = 1.0
1486	sbt = 1.0
1487	ENDIF
1488	tend = 0.0
1489	CALL advec_s_bc( e, 'e' )
1490	ELSE
1491	IF ( tsc(2) /= 2.0 ) THEN
1492	IF ( scalar_advec == 'ups-scheme' ) THEN
1493	tend = 0.0
1494	CALL advec_s_ups( e, 'e' )
1495	ENDIF
1496	ENDIF
1497	ENDIF
1498	ENDIF
1499
1500	!
1501	!-- TKE-tendency terms with no communication
1502	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
1503	THEN
1504	IF ( .NOT. humidity ) THEN
1505	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1506	rif, tend, zu, zw )
1507	ELSE
1508	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1509	rif, tend, zu, zw )
1510	ENDIF
1511	ELSE
1512	IF ( use_upstream_for_tke ) THEN
1513	tend = 0.0
1514	CALL advec_s_up( e )
1515	ELSE
1516	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1517	tend = 0.0
1518	CALL advec_s_pw( e )
1519	ELSE
1520	IF ( scalar_advec /= 'ups-scheme' ) THEN
1521	tend = 0.0
1522	CALL advec_s_up( e )
1523	ENDIF
1524	ENDIF
1525	ENDIF
1526	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1527	IF ( .NOT. humidity ) THEN
1528	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1529	pt_m, rif_m, tend, zu, zw )
1530	ELSE
1531	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1532	vpt_m, rif_m, tend, zu, zw )
1533	ENDIF
1534	ELSE
1535	IF ( .NOT. humidity ) THEN
1536	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1537	rif, tend, zu, zw )
1538	ELSE
1539	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1540	rif, tend, zu, zw )
1541	ENDIF
1542	ENDIF
1543	ENDIF
1544	CALL production_e
1545	CALL user_actions( 'e-tendency' )
1546
1547	!
1548	!-- Prognostic equation for TKE.
1549	!-- Eliminate negative TKE values, which can occur due to numerical
1550	!-- reasons in the course of the integration. In such cases the old TKE
1551	!-- value is reduced by 90%.
1552	DO i = nxl, nxr
1553	DO j = nys, nyn
1554	DO k = nzb_s_inner(j,i)+1, nzt
1555	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
1556	dt_3d * ( &
1557	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1558	)
1559	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1560	ENDDO
1561	ENDDO
1562	ENDDO
1563
1564	!
1565	!-- Calculate tendencies for the next Runge-Kutta step
1566	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1567	IF ( intermediate_timestep_count == 1 ) THEN
1568	DO i = nxl, nxr
1569	DO j = nys, nyn
1570	DO k = nzb_s_inner(j,i)+1, nzt
1571	te_m(k,j,i) = tend(k,j,i)
1572	ENDDO
1573	ENDDO
1574	ENDDO
1575	ELSEIF ( intermediate_timestep_count < &
1576	intermediate_timestep_count_max ) THEN
1577	DO i = nxl, nxr
1578	DO j = nys, nyn
1579	DO k = nzb_s_inner(j,i)+1, nzt
1580	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1581	ENDDO
1582	ENDDO
1583	ENDDO
1584	ENDIF
1585	ENDIF
1586
1587	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1588
1589	ENDIF
1590
1591
1592	END SUBROUTINE prognostic_equations_vector
1593
1594
1595	END MODULE prognostic_equations_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |