Home

Context Navigation

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

Last change on this file since 70 was 70, checked in by raasch, 18 years ago
bugs fixed for particle code and bc-scheme
Property svn:keywords set to `Id`
File size: 54.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |