Home

Context Navigation

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

Last change on this file since 75 was 75, checked in by raasch, 17 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |