Home

Context Navigation

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

Last change on this file since 63 was 63, checked in by raasch, 17 years ago
preliminary changes concerning update of BC-scheme
Property svn:keywords set to `Id`
File size: 53.9 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 63 2007-03-13 03:52:49Z 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	IF ( scalar_advec == 'bc-scheme' ) THEN
319	!
320	!-- Bott-Chlond scheme always uses Euler time step. Thus:
321	sat = 1.0
322	sbt = 1.0
323	tend = 0.0
324	CALL advec_s_bc( pt, 'pt' )
325	ELSE
326	sat = tsc(1)
327	sbt = tsc(2)
328	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
329	tend = 0.0
330	CALL advec_s_ups( pt, 'pt' )
331	ENDIF
332	ENDIF
333
334	!
335	!-- pt-tendency terms with no communication
336	DO i = nxl, nxr
337	DO j = nys, nyn
338	!
339	!-- Tendency terms
340	IF ( scalar_advec == 'bc-scheme' ) THEN
341	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
342	ELSE
343	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
344	tend(:,j,i) = 0.0
345	CALL advec_s_pw( i, j, pt )
346	ELSE
347	IF ( scalar_advec /= 'ups-scheme' ) THEN
348	tend(:,j,i) = 0.0
349	CALL advec_s_up( i, j, pt )
350	ENDIF
351	ENDIF
352	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
353	THEN
354	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
355	tswst_m, tend )
356	ELSE
357	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
358	ENDIF
359	ENDIF
360
361	!
362	!-- If required compute heating/cooling due to long wave radiation
363	!-- processes
364	IF ( radiation ) THEN
365	CALL calc_radiation( i, j )
366	ENDIF
367
368	!
369	!-- If required compute impact of latent heat due to precipitation
370	IF ( precipitation ) THEN
371	CALL impact_of_latent_heat( i, j )
372	ENDIF
373	CALL user_actions( i, j, 'pt-tendency' )
374
375	!
376	!-- Prognostic equation for potential temperature
377	DO k = nzb_s_inner(j,i)+1, nzt
378	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
379	dt_3d * ( &
380	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
381	) - &
382	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
383	ENDDO
384
385	!
386	!-- Calculate tendencies for the next Runge-Kutta step
387	IF ( timestep_scheme(1:5) == 'runge' ) THEN
388	IF ( intermediate_timestep_count == 1 ) THEN
389	DO k = nzb_s_inner(j,i)+1, nzt
390	tpt_m(k,j,i) = tend(k,j,i)
391	ENDDO
392	ELSEIF ( intermediate_timestep_count < &
393	intermediate_timestep_count_max ) THEN
394	DO k = nzb_s_inner(j,i)+1, nzt
395	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
396	ENDDO
397	ENDIF
398	ENDIF
399
400	ENDDO
401	ENDDO
402
403	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
404
405	!
406	!-- If required, compute prognostic equation for total water content / scalar
407	IF ( moisture .OR. passive_scalar ) THEN
408
409	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
410
411	!
412	!-- Scalar/q-tendency terms with communication
413	IF ( scalar_advec == 'bc-scheme' ) THEN
414	!
415	!-- Bott-Chlond scheme always uses Euler time step. Thus:
416	sat = 1.0
417	sbt = 1.0
418	tend = 0.0
419	CALL advec_s_bc( q, 'q' )
420	ELSE
421	sat = tsc(1)
422	sbt = tsc(2)
423	IF ( tsc(2) /= 2.0 ) THEN
424	IF ( scalar_advec == 'ups-scheme' ) THEN
425	tend = 0.0
426	CALL advec_s_ups( q, 'q' )
427	ENDIF
428	ENDIF
429	ENDIF
430
431	!
432	!-- Scalar/q-tendency terms with no communication
433	DO i = nxl, nxr
434	DO j = nys, nyn
435	!
436	!-- Tendency-terms
437	IF ( scalar_advec == 'bc-scheme' ) THEN
438	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, tend )
439	ELSE
440	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
441	tend(:,j,i) = 0.0
442	CALL advec_s_pw( i, j, q )
443	ELSE
444	IF ( scalar_advec /= 'ups-scheme' ) THEN
445	tend(:,j,i) = 0.0
446	CALL advec_s_up( i, j, q )
447	ENDIF
448	ENDIF
449	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
450	THEN
451	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
452	qswst_m, tend )
453	ELSE
454	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
455	tend )
456	ENDIF
457	ENDIF
458
459	!
460	!-- If required compute decrease of total water content due to
461	!-- precipitation
462	IF ( precipitation ) THEN
463	CALL calc_precipitation( i, j )
464	ENDIF
465	CALL user_actions( i, j, 'q-tendency' )
466
467	!
468	!-- Prognostic equation for total water content / scalar
469	DO k = nzb_s_inner(j,i)+1, nzt
470	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
471	dt_3d * ( &
472	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
473	) - &
474	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
475	ENDDO
476
477	!
478	!-- Calculate tendencies for the next Runge-Kutta step
479	IF ( timestep_scheme(1:5) == 'runge' ) THEN
480	IF ( intermediate_timestep_count == 1 ) THEN
481	DO k = nzb_s_inner(j,i)+1, nzt
482	tq_m(k,j,i) = tend(k,j,i)
483	ENDDO
484	ELSEIF ( intermediate_timestep_count < &
485	intermediate_timestep_count_max ) THEN
486	DO k = nzb_s_inner(j,i)+1, nzt
487	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
488	ENDDO
489	ENDIF
490	ENDIF
491
492	ENDDO
493	ENDDO
494
495	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
496
497	ENDIF
498
499	!
500	!-- If required, compute prognostic equation for turbulent kinetic
501	!-- energy (TKE)
502	IF ( .NOT. constant_diffusion ) THEN
503
504	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
505
506	!
507	!-- TKE-tendency terms with communication
508	CALL production_e_init
509	IF ( .NOT. use_upstream_for_tke ) THEN
510	IF ( scalar_advec == 'bc-scheme' ) THEN
511	!
512	!-- Bott-Chlond scheme always uses Euler time step. Thus:
513	sat = 1.0
514	sbt = 1.0
515	tend = 0.0
516	CALL advec_s_bc( e, 'e' )
517	ELSE
518	sat = tsc(1)
519	sbt = tsc(2)
520	IF ( tsc(2) /= 2.0 ) THEN
521	IF ( scalar_advec == 'ups-scheme' ) THEN
522	tend = 0.0
523	CALL advec_s_ups( e, 'e' )
524	ENDIF
525	ENDIF
526	ENDIF
527	ENDIF
528
529	!
530	!-- TKE-tendency terms with no communication
531	DO i = nxl, nxr
532	DO j = nys, nyn
533	!
534	!-- Tendency-terms
535	IF ( scalar_advec == 'bc-scheme' .AND. &
536	.NOT. use_upstream_for_tke ) THEN
537	IF ( .NOT. moisture ) THEN
538	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
539	l_grid, pt, rif, tend, zu )
540	ELSE
541	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
542	l_grid, vpt, rif, tend, zu )
543	ENDIF
544	ELSE
545	IF ( use_upstream_for_tke ) THEN
546	tend(:,j,i) = 0.0
547	CALL advec_s_up( i, j, e )
548	ELSE
549	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
550	THEN
551	tend(:,j,i) = 0.0
552	CALL advec_s_pw( i, j, e )
553	ELSE
554	IF ( scalar_advec /= 'ups-scheme' ) THEN
555	tend(:,j,i) = 0.0
556	CALL advec_s_up( i, j, e )
557	ENDIF
558	ENDIF
559	ENDIF
560	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
561	THEN
562	IF ( .NOT. moisture ) THEN
563	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
564	km_m, l_grid, pt_m, rif_m, tend, zu )
565	ELSE
566	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
567	km_m, l_grid, vpt_m, rif_m, tend, zu )
568	ENDIF
569	ELSE
570	IF ( .NOT. moisture ) THEN
571	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
572	l_grid, pt, rif, tend, zu )
573	ELSE
574	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
575	l_grid, vpt, rif, tend, zu )
576	ENDIF
577	ENDIF
578	ENDIF
579	CALL production_e( i, j )
580	CALL user_actions( i, j, 'e-tendency' )
581
582	!
583	!-- Prognostic equation for TKE.
584	!-- Eliminate negative TKE values, which can occur due to numerical
585	!-- reasons in the course of the integration. In such cases the old TKE
586	!-- value is reduced by 90%.
587	DO k = nzb_s_inner(j,i)+1, nzt
588	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
589	dt_3d * ( &
590	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
591	)
592	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
593	ENDDO
594
595	!
596	!-- Calculate tendencies for the next Runge-Kutta step
597	IF ( timestep_scheme(1:5) == 'runge' ) THEN
598	IF ( intermediate_timestep_count == 1 ) THEN
599	DO k = nzb_s_inner(j,i)+1, nzt
600	te_m(k,j,i) = tend(k,j,i)
601	ENDDO
602	ELSEIF ( intermediate_timestep_count < &
603	intermediate_timestep_count_max ) THEN
604	DO k = nzb_s_inner(j,i)+1, nzt
605	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
606	ENDDO
607	ENDIF
608	ENDIF
609
610	ENDDO
611	ENDDO
612
613	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
614
615	ENDIF
616
617
618	END SUBROUTINE prognostic_equations_noopt
619
620
621	SUBROUTINE prognostic_equations_cache
622
623	!------------------------------------------------------------------------------!
624	! Version with one optimized loop over all equations. It is only allowed to
625	! be called for the standard Piascek-Williams advection scheme.
626	!
627	! The call of this subroutine is embedded in two DO loops over i and j, thus
628	! communication between CPUs is not allowed in this subroutine.
629	!
630	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
631	!------------------------------------------------------------------------------!
632
633	IMPLICIT NONE
634
635	CHARACTER (LEN=9) :: time_to_string
636	INTEGER :: i, j, k
637
638
639	!
640	!-- Time measurement can only be performed for the whole set of equations
641	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
642
643
644	!
645	!-- Calculate those variables needed in the tendency terms which need
646	!-- global communication
647	CALL calc_mean_pt_profile( pt, 4 )
648	IF ( moisture ) CALL calc_mean_pt_profile( vpt, 44 )
649	IF ( .NOT. constant_diffusion ) CALL production_e_init
650
651
652	!
653	!-- Loop over all prognostic equations
654	!$OMP PARALLEL private (i,j,k)
655	!$OMP DO
656	DO i = nxl, nxr+uxrp ! Additional levels for non cyclic boundary
657	DO j = nys, nyn+vynp ! conditions are included
658	!
659	!-- Tendency terms for u-velocity component
660	IF ( j < nyn+1 ) THEN
661
662	tend(:,j,i) = 0.0
663	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
664	CALL advec_u_pw( i, j )
665	ELSE
666	CALL advec_u_up( i, j )
667	ENDIF
668	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
669	THEN
670	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
671	u_m, usws_m, v_m, w_m )
672	ELSE
673	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
674	usws, v, w )
675	ENDIF
676	CALL coriolis( i, j, 1 )
677	IF ( sloping_surface ) CALL buoyancy( i, j, pt, 1, 4 )
678	CALL user_actions( i, j, 'u-tendency' )
679
680	!
681	!-- Prognostic equation for u-velocity component
682	DO k = nzb_u_inner(j,i)+1, nzt
683	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
684	dt_3d * ( &
685	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
686	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
687	) - &
688	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
689	ENDDO
690
691	!
692	!-- Calculate tendencies for the next Runge-Kutta step
693	IF ( timestep_scheme(1:5) == 'runge' ) THEN
694	IF ( intermediate_timestep_count == 1 ) THEN
695	DO k = nzb_u_inner(j,i)+1, nzt
696	tu_m(k,j,i) = tend(k,j,i)
697	ENDDO
698	ELSEIF ( intermediate_timestep_count < &
699	intermediate_timestep_count_max ) THEN
700	DO k = nzb_u_inner(j,i)+1, nzt
701	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
702	ENDDO
703	ENDIF
704	ENDIF
705
706	ENDIF
707
708	!
709	!-- Tendency terms for v-velocity component
710	IF ( i < nxr+1 ) THEN
711
712	tend(:,j,i) = 0.0
713	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
714	CALL advec_v_pw( i, j )
715	ELSE
716	CALL advec_v_up( i, j )
717	ENDIF
718	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
719	THEN
720	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
721	u_m, v_m, vsws_m, w_m )
722	ELSE
723	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
724	vsws, w )
725	ENDIF
726	CALL coriolis( i, j, 2 )
727	CALL user_actions( i, j, 'v-tendency' )
728
729	!
730	!-- Prognostic equation for v-velocity component
731	DO k = nzb_v_inner(j,i)+1, nzt
732	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
733	dt_3d * ( &
734	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
735	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
736	) - &
737	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
738	ENDDO
739
740	!
741	!-- Calculate tendencies for the next Runge-Kutta step
742	IF ( timestep_scheme(1:5) == 'runge' ) THEN
743	IF ( intermediate_timestep_count == 1 ) THEN
744	DO k = nzb_v_inner(j,i)+1, nzt
745	tv_m(k,j,i) = tend(k,j,i)
746	ENDDO
747	ELSEIF ( intermediate_timestep_count < &
748	intermediate_timestep_count_max ) THEN
749	DO k = nzb_v_inner(j,i)+1, nzt
750	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
751	ENDDO
752	ENDIF
753	ENDIF
754
755	ENDIF
756
757	!
758	!-- Tendency terms for w-velocity component
759	IF ( i < nxr+1 .AND. j < nyn+1 ) THEN
760
761	tend(:,j,i) = 0.0
762	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
763	CALL advec_w_pw( i, j )
764	ELSE
765	CALL advec_w_up( i, j )
766	ENDIF
767	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
768	THEN
769	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
770	km_damp_y, tend, u_m, v_m, w_m )
771	ELSE
772	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
773	tend, u, v, w )
774	ENDIF
775	CALL coriolis( i, j, 3 )
776	IF ( .NOT. moisture ) THEN
777	CALL buoyancy( i, j, pt, 3, 4 )
778	ELSE
779	CALL buoyancy( i, j, vpt, 3, 44 )
780	ENDIF
781	CALL user_actions( i, j, 'w-tendency' )
782
783	!
784	!-- Prognostic equation for w-velocity component
785	DO k = nzb_w_inner(j,i)+1, nzt-1
786	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
787	dt_3d * ( &
788	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
789	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
790	) - &
791	tsc(5) * rdf(k) * w(k,j,i)
792	ENDDO
793
794	!
795	!-- Calculate tendencies for the next Runge-Kutta step
796	IF ( timestep_scheme(1:5) == 'runge' ) THEN
797	IF ( intermediate_timestep_count == 1 ) THEN
798	DO k = nzb_w_inner(j,i)+1, nzt-1
799	tw_m(k,j,i) = tend(k,j,i)
800	ENDDO
801	ELSEIF ( intermediate_timestep_count < &
802	intermediate_timestep_count_max ) THEN
803	DO k = nzb_w_inner(j,i)+1, nzt-1
804	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
805	ENDDO
806	ENDIF
807	ENDIF
808
809	!
810	!-- Tendency terms for potential temperature
811	tend(:,j,i) = 0.0
812	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
813	CALL advec_s_pw( i, j, pt )
814	ELSE
815	CALL advec_s_up( i, j, pt )
816	ENDIF
817	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
818	THEN
819	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
820	tswst_m, tend )
821	ELSE
822	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, tend )
823	ENDIF
824
825	!
826	!-- If required compute heating/cooling due to long wave radiation
827	!-- processes
828	IF ( radiation ) THEN
829	CALL calc_radiation( i, j )
830	ENDIF
831
832	!
833	!-- If required compute impact of latent heat due to precipitation
834	IF ( precipitation ) THEN
835	CALL impact_of_latent_heat( i, j )
836	ENDIF
837	CALL user_actions( i, j, 'pt-tendency' )
838
839	!
840	!-- Prognostic equation for potential temperature
841	DO k = nzb_s_inner(j,i)+1, nzt
842	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
843	dt_3d * ( &
844	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
845	) - &
846	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
847	ENDDO
848
849	!
850	!-- Calculate tendencies for the next Runge-Kutta step
851	IF ( timestep_scheme(1:5) == 'runge' ) THEN
852	IF ( intermediate_timestep_count == 1 ) THEN
853	DO k = nzb_s_inner(j,i)+1, nzt
854	tpt_m(k,j,i) = tend(k,j,i)
855	ENDDO
856	ELSEIF ( intermediate_timestep_count < &
857	intermediate_timestep_count_max ) THEN
858	DO k = nzb_s_inner(j,i)+1, nzt
859	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
860	5.3125 * tpt_m(k,j,i)
861	ENDDO
862	ENDIF
863	ENDIF
864
865	!
866	!-- If required, compute prognostic equation for total water content /
867	!-- scalar
868	IF ( moisture .OR. passive_scalar ) THEN
869
870	!
871	!-- Tendency-terms for total water content / scalar
872	tend(:,j,i) = 0.0
873	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
874	THEN
875	CALL advec_s_pw( i, j, q )
876	ELSE
877	CALL advec_s_up( i, j, q )
878	ENDIF
879	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
880	THEN
881	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
882	qswst_m, tend )
883	ELSE
884	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
885	tend )
886	ENDIF
887
888	!
889	!-- If required compute decrease of total water content due to
890	!-- precipitation
891	IF ( precipitation ) THEN
892	CALL calc_precipitation( i, j )
893	ENDIF
894	CALL user_actions( i, j, 'q-tendency' )
895
896	!
897	!-- Prognostic equation for total water content / scalar
898	DO k = nzb_s_inner(j,i)+1, nzt
899	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
900	dt_3d * ( &
901	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
902	) - &
903	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
904	ENDDO
905
906	!
907	!-- Calculate tendencies for the next Runge-Kutta step
908	IF ( timestep_scheme(1:5) == 'runge' ) THEN
909	IF ( intermediate_timestep_count == 1 ) THEN
910	DO k = nzb_s_inner(j,i)+1, nzt
911	tq_m(k,j,i) = tend(k,j,i)
912	ENDDO
913	ELSEIF ( intermediate_timestep_count < &
914	intermediate_timestep_count_max ) THEN
915	DO k = nzb_s_inner(j,i)+1, nzt
916	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
917	5.3125 * tq_m(k,j,i)
918	ENDDO
919	ENDIF
920	ENDIF
921
922	ENDIF
923
924	!
925	!-- If required, compute prognostic equation for turbulent kinetic
926	!-- energy (TKE)
927	IF ( .NOT. constant_diffusion ) THEN
928
929	!
930	!-- Tendency-terms for TKE
931	tend(:,j,i) = 0.0
932	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
933	.AND. .NOT. use_upstream_for_tke ) THEN
934	CALL advec_s_pw( i, j, e )
935	ELSE
936	CALL advec_s_up( i, j, e )
937	ENDIF
938	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
939	THEN
940	IF ( .NOT. moisture ) THEN
941	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
942	km_m, l_grid, pt_m, rif_m, tend, zu )
943	ELSE
944	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
945	km_m, l_grid, vpt_m, rif_m, tend, zu )
946	ENDIF
947	ELSE
948	IF ( .NOT. moisture ) THEN
949	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
950	l_grid, pt, rif, tend, zu )
951	ELSE
952	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
953	l_grid, vpt, rif, tend, zu )
954	ENDIF
955	ENDIF
956	CALL production_e( i, j )
957	CALL user_actions( i, j, 'e-tendency' )
958
959	!
960	!-- Prognostic equation for TKE.
961	!-- Eliminate negative TKE values, which can occur due to numerical
962	!-- reasons in the course of the integration. In such cases the old
963	!-- TKE value is reduced by 90%.
964	DO k = nzb_s_inner(j,i)+1, nzt
965	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
966	dt_3d * ( &
967	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
968	)
969	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
970	ENDDO
971
972	!
973	!-- Calculate tendencies for the next Runge-Kutta step
974	IF ( timestep_scheme(1:5) == 'runge' ) THEN
975	IF ( intermediate_timestep_count == 1 ) THEN
976	DO k = nzb_s_inner(j,i)+1, nzt
977	te_m(k,j,i) = tend(k,j,i)
978	ENDDO
979	ELSEIF ( intermediate_timestep_count < &
980	intermediate_timestep_count_max ) THEN
981	DO k = nzb_s_inner(j,i)+1, nzt
982	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
983	5.3125 * te_m(k,j,i)
984	ENDDO
985	ENDIF
986	ENDIF
987
988	ENDIF ! TKE equation
989
990	ENDIF ! Gridpoints excluding the non-cyclic wall
991
992	ENDDO
993	ENDDO
994	!$OMP END PARALLEL
995
996	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
997
998
999	END SUBROUTINE prognostic_equations_cache
1000
1001
1002	SUBROUTINE prognostic_equations_vector
1003
1004	!------------------------------------------------------------------------------!
1005	! Version for vector machines
1006	!------------------------------------------------------------------------------!
1007
1008	IMPLICIT NONE
1009
1010	CHARACTER (LEN=9) :: time_to_string
1011	INTEGER :: i, j, k
1012	REAL :: sat, sbt
1013
1014	!
1015	!-- Calculate those variables needed in the tendency terms which need
1016	!-- global communication
1017	CALL calc_mean_pt_profile( pt, 4 )
1018	IF ( moisture ) CALL calc_mean_pt_profile( vpt, 44 )
1019
1020	!
1021	!-- u-velocity component
1022	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1023
1024	!
1025	!-- u-tendency terms with communication
1026	IF ( momentum_advec == 'ups-scheme' ) THEN
1027	tend = 0.0
1028	CALL advec_u_ups
1029	ENDIF
1030
1031	!
1032	!-- u-tendency terms with no communication
1033	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1034	tend = 0.0
1035	CALL advec_u_pw
1036	ELSE
1037	IF ( momentum_advec /= 'ups-scheme' ) THEN
1038	tend = 0.0
1039	CALL advec_u_up
1040	ENDIF
1041	ENDIF
1042	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1043	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
1044	w_m )
1045	ELSE
1046	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
1047	ENDIF
1048	CALL coriolis( 1 )
1049	IF ( sloping_surface ) CALL buoyancy( pt, 1, 4 )
1050	CALL user_actions( 'u-tendency' )
1051
1052	!
1053	!-- Prognostic equation for u-velocity component
1054	DO i = nxl, nxr+uxrp
1055	DO j = nys, nyn
1056	DO k = nzb_u_inner(j,i)+1, nzt
1057	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1058	dt_3d * ( &
1059	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1060	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
1061	) - &
1062	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1063	ENDDO
1064	ENDDO
1065	ENDDO
1066
1067	!
1068	!-- Calculate tendencies for the next Runge-Kutta step
1069	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1070	IF ( intermediate_timestep_count == 1 ) THEN
1071	DO i = nxl, nxr+uxrp
1072	DO j = nys, nyn
1073	DO k = nzb_u_inner(j,i)+1, nzt
1074	tu_m(k,j,i) = tend(k,j,i)
1075	ENDDO
1076	ENDDO
1077	ENDDO
1078	ELSEIF ( intermediate_timestep_count < &
1079	intermediate_timestep_count_max ) THEN
1080	DO i = nxl, nxr+uxrp
1081	DO j = nys, nyn
1082	DO k = nzb_u_inner(j,i)+1, nzt
1083	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1084	ENDDO
1085	ENDDO
1086	ENDDO
1087	ENDIF
1088	ENDIF
1089
1090	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1091
1092	!
1093	!-- v-velocity component
1094	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1095
1096	!
1097	!-- v-tendency terms with communication
1098	IF ( momentum_advec == 'ups-scheme' ) THEN
1099	tend = 0.0
1100	CALL advec_v_ups
1101	ENDIF
1102
1103	!
1104	!-- v-tendency terms with no communication
1105	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1106	tend = 0.0
1107	CALL advec_v_pw
1108	ELSE
1109	IF ( momentum_advec /= 'ups-scheme' ) THEN
1110	tend = 0.0
1111	CALL advec_v_up
1112	ENDIF
1113	ENDIF
1114	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1115	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1116	w_m )
1117	ELSE
1118	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
1119	ENDIF
1120	CALL coriolis( 2 )
1121	CALL user_actions( 'v-tendency' )
1122
1123	!
1124	!-- Prognostic equation for v-velocity component
1125	DO i = nxl, nxr
1126	DO j = nys, nyn+vynp
1127	DO k = nzb_v_inner(j,i)+1, nzt
1128	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1129	dt_3d * ( &
1130	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1131	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
1132	) - &
1133	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1134	ENDDO
1135	ENDDO
1136	ENDDO
1137
1138	!
1139	!-- Calculate tendencies for the next Runge-Kutta step
1140	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1141	IF ( intermediate_timestep_count == 1 ) THEN
1142	DO i = nxl, nxr
1143	DO j = nys, nyn+vynp
1144	DO k = nzb_v_inner(j,i)+1, nzt
1145	tv_m(k,j,i) = tend(k,j,i)
1146	ENDDO
1147	ENDDO
1148	ENDDO
1149	ELSEIF ( intermediate_timestep_count < &
1150	intermediate_timestep_count_max ) THEN
1151	DO i = nxl, nxr
1152	DO j = nys, nyn+vynp
1153	DO k = nzb_v_inner(j,i)+1, nzt
1154	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1155	ENDDO
1156	ENDDO
1157	ENDDO
1158	ENDIF
1159	ENDIF
1160
1161	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1162
1163	!
1164	!-- w-velocity component
1165	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1166
1167	!
1168	!-- w-tendency terms with communication
1169	IF ( momentum_advec == 'ups-scheme' ) THEN
1170	tend = 0.0
1171	CALL advec_w_ups
1172	ENDIF
1173
1174	!
1175	!-- w-tendency terms with no communication
1176	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1177	tend = 0.0
1178	CALL advec_w_pw
1179	ELSE
1180	IF ( momentum_advec /= 'ups-scheme' ) THEN
1181	tend = 0.0
1182	CALL advec_w_up
1183	ENDIF
1184	ENDIF
1185	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1186	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1187	v_m, w_m )
1188	ELSE
1189	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1190	ENDIF
1191	CALL coriolis( 3 )
1192	IF ( .NOT. moisture ) THEN
1193	CALL buoyancy( pt, 3, 4 )
1194	ELSE
1195	CALL buoyancy( vpt, 3, 44 )
1196	ENDIF
1197	CALL user_actions( 'w-tendency' )
1198
1199	!
1200	!-- Prognostic equation for w-velocity component
1201	DO i = nxl, nxr
1202	DO j = nys, nyn
1203	DO k = nzb_w_inner(j,i)+1, nzt-1
1204	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1205	dt_3d * ( &
1206	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1207	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1208	) - &
1209	tsc(5) * rdf(k) * w(k,j,i)
1210	ENDDO
1211	ENDDO
1212	ENDDO
1213
1214	!
1215	!-- Calculate tendencies for the next Runge-Kutta step
1216	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1217	IF ( intermediate_timestep_count == 1 ) THEN
1218	DO i = nxl, nxr
1219	DO j = nys, nyn
1220	DO k = nzb_w_inner(j,i)+1, nzt-1
1221	tw_m(k,j,i) = tend(k,j,i)
1222	ENDDO
1223	ENDDO
1224	ENDDO
1225	ELSEIF ( intermediate_timestep_count < &
1226	intermediate_timestep_count_max ) THEN
1227	DO i = nxl, nxr
1228	DO j = nys, nyn
1229	DO k = nzb_w_inner(j,i)+1, nzt-1
1230	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1231	ENDDO
1232	ENDDO
1233	ENDDO
1234	ENDIF
1235	ENDIF
1236
1237	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1238
1239	!
1240	!-- potential temperature
1241	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1242
1243	!
1244	!-- pt-tendency terms with communication
1245	sat = tsc(1)
1246	sbt = tsc(2)
1247	IF ( scalar_advec == 'bc-scheme' ) THEN
1248
1249	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1250	!
1251	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1252	!-- switched on. Thus:
1253	sat = 1.0
1254	sbt = 1.0
1255	ENDIF
1256	tend = 0.0
1257	CALL advec_s_bc( pt, 'pt' )
1258	ELSE
1259	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
1260	tend = 0.0
1261	CALL advec_s_ups( pt, 'pt' )
1262	ENDIF
1263	ENDIF
1264
1265	!
1266	!-- pt-tendency terms with no communication
1267	IF ( scalar_advec == 'bc-scheme' ) THEN
1268	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1269	ELSE
1270	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1271	tend = 0.0
1272	CALL advec_s_pw( pt )
1273	ELSE
1274	IF ( scalar_advec /= 'ups-scheme' ) THEN
1275	tend = 0.0
1276	CALL advec_s_up( pt )
1277	ENDIF
1278	ENDIF
1279	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1280	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, tend )
1281	ELSE
1282	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, tend )
1283	ENDIF
1284	ENDIF
1285
1286	!
1287	!-- If required compute heating/cooling due to long wave radiation
1288	!-- processes
1289	IF ( radiation ) THEN
1290	CALL calc_radiation
1291	ENDIF
1292
1293	!
1294	!-- If required compute impact of latent heat due to precipitation
1295	IF ( precipitation ) THEN
1296	CALL impact_of_latent_heat
1297	ENDIF
1298	CALL user_actions( 'pt-tendency' )
1299
1300	!
1301	!-- Prognostic equation for potential temperature
1302	DO i = nxl, nxr
1303	DO j = nys, nyn
1304	DO k = nzb_s_inner(j,i)+1, nzt
1305	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
1306	dt_3d * ( &
1307	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1308	) - &
1309	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1310	ENDDO
1311	ENDDO
1312	ENDDO
1313
1314	!
1315	!-- Calculate tendencies for the next Runge-Kutta step
1316	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1317	IF ( intermediate_timestep_count == 1 ) THEN
1318	DO i = nxl, nxr
1319	DO j = nys, nyn
1320	DO k = nzb_s_inner(j,i)+1, nzt
1321	tpt_m(k,j,i) = tend(k,j,i)
1322	ENDDO
1323	ENDDO
1324	ENDDO
1325	ELSEIF ( intermediate_timestep_count < &
1326	intermediate_timestep_count_max ) THEN
1327	DO i = nxl, nxr
1328	DO j = nys, nyn
1329	DO k = nzb_s_inner(j,i)+1, nzt
1330	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1331	ENDDO
1332	ENDDO
1333	ENDDO
1334	ENDIF
1335	ENDIF
1336
1337	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1338
1339	!
1340	!-- If required, compute prognostic equation for total water content / scalar
1341	IF ( moisture .OR. passive_scalar ) THEN
1342
1343	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1344
1345	!
1346	!-- Scalar/q-tendency terms with communication
1347	sat = tsc(1)
1348	sbt = tsc(2)
1349	IF ( scalar_advec == 'bc-scheme' ) THEN
1350
1351	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1352	!
1353	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1354	!-- switched on. Thus:
1355	sat = 1.0
1356	sbt = 1.0
1357	ENDIF
1358	tend = 0.0
1359	CALL advec_s_bc( q, 'q' )
1360	ELSE
1361	IF ( tsc(2) /= 2.0 ) THEN
1362	IF ( scalar_advec == 'ups-scheme' ) THEN
1363	tend = 0.0
1364	CALL advec_s_ups( q, 'q' )
1365	ENDIF
1366	ENDIF
1367	ENDIF
1368
1369	!
1370	!-- Scalar/q-tendency terms with no communication
1371	IF ( scalar_advec == 'bc-scheme' ) THEN
1372	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1373	ELSE
1374	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1375	tend = 0.0
1376	CALL advec_s_pw( q )
1377	ELSE
1378	IF ( scalar_advec /= 'ups-scheme' ) THEN
1379	tend = 0.0
1380	CALL advec_s_up( q )
1381	ENDIF
1382	ENDIF
1383	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1384	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, tend )
1385	ELSE
1386	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, tend )
1387	ENDIF
1388	ENDIF
1389
1390	!
1391	!-- If required compute decrease of total water content due to
1392	!-- precipitation
1393	IF ( precipitation ) THEN
1394	CALL calc_precipitation
1395	ENDIF
1396	CALL user_actions( 'q-tendency' )
1397
1398	!
1399	!-- Prognostic equation for total water content / scalar
1400	DO i = nxl, nxr
1401	DO j = nys, nyn
1402	DO k = nzb_s_inner(j,i)+1, nzt
1403	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
1404	dt_3d * ( &
1405	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1406	) - &
1407	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1408	ENDDO
1409	ENDDO
1410	ENDDO
1411
1412	!
1413	!-- Calculate tendencies for the next Runge-Kutta step
1414	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1415	IF ( intermediate_timestep_count == 1 ) THEN
1416	DO i = nxl, nxr
1417	DO j = nys, nyn
1418	DO k = nzb_s_inner(j,i)+1, nzt
1419	tq_m(k,j,i) = tend(k,j,i)
1420	ENDDO
1421	ENDDO
1422	ENDDO
1423	ELSEIF ( intermediate_timestep_count < &
1424	intermediate_timestep_count_max ) THEN
1425	DO i = nxl, nxr
1426	DO j = nys, nyn
1427	DO k = nzb_s_inner(j,i)+1, nzt
1428	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1429	ENDDO
1430	ENDDO
1431	ENDDO
1432	ENDIF
1433	ENDIF
1434
1435	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1436
1437	ENDIF
1438
1439	!
1440	!-- If required, compute prognostic equation for turbulent kinetic
1441	!-- energy (TKE)
1442	IF ( .NOT. constant_diffusion ) THEN
1443
1444	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1445
1446	!
1447	!-- TKE-tendency terms with communication
1448	CALL production_e_init
1449
1450	sat = tsc(1)
1451	sbt = tsc(2)
1452	IF ( .NOT. use_upstream_for_tke ) THEN
1453	IF ( scalar_advec == 'bc-scheme' ) THEN
1454
1455	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1456	!
1457	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1458	!-- switched on. Thus:
1459	sat = 1.0
1460	sbt = 1.0
1461	ENDIF
1462	tend = 0.0
1463	CALL advec_s_bc( e, 'e' )
1464	ELSE
1465	IF ( tsc(2) /= 2.0 ) THEN
1466	IF ( scalar_advec == 'ups-scheme' ) THEN
1467	tend = 0.0
1468	CALL advec_s_ups( e, 'e' )
1469	ENDIF
1470	ENDIF
1471	ENDIF
1472	ENDIF
1473
1474	!
1475	!-- TKE-tendency terms with no communication
1476	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
1477	THEN
1478	IF ( .NOT. moisture ) THEN
1479	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1480	rif, tend, zu )
1481	ELSE
1482	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1483	rif, tend, zu )
1484	ENDIF
1485	ELSE
1486	IF ( use_upstream_for_tke ) THEN
1487	tend = 0.0
1488	CALL advec_s_up( e )
1489	ELSE
1490	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1491	tend = 0.0
1492	CALL advec_s_pw( e )
1493	ELSE
1494	IF ( scalar_advec /= 'ups-scheme' ) THEN
1495	tend = 0.0
1496	CALL advec_s_up( e )
1497	ENDIF
1498	ENDIF
1499	ENDIF
1500	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1501	IF ( .NOT. moisture ) THEN
1502	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1503	pt_m, rif_m, tend, zu )
1504	ELSE
1505	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1506	vpt_m, rif_m, tend, zu )
1507	ENDIF
1508	ELSE
1509	IF ( .NOT. moisture ) THEN
1510	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1511	rif, tend, zu )
1512	ELSE
1513	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1514	rif, tend, zu )
1515	ENDIF
1516	ENDIF
1517	ENDIF
1518	CALL production_e
1519	CALL user_actions( 'e-tendency' )
1520
1521	!
1522	!-- Prognostic equation for TKE.
1523	!-- Eliminate negative TKE values, which can occur due to numerical
1524	!-- reasons in the course of the integration. In such cases the old TKE
1525	!-- value is reduced by 90%.
1526	DO i = nxl, nxr
1527	DO j = nys, nyn
1528	DO k = nzb_s_inner(j,i)+1, nzt
1529	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
1530	dt_3d * ( &
1531	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1532	)
1533	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1534	ENDDO
1535	ENDDO
1536	ENDDO
1537
1538	!
1539	!-- Calculate tendencies for the next Runge-Kutta step
1540	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1541	IF ( intermediate_timestep_count == 1 ) THEN
1542	DO i = nxl, nxr
1543	DO j = nys, nyn
1544	DO k = nzb_s_inner(j,i)+1, nzt
1545	te_m(k,j,i) = tend(k,j,i)
1546	ENDDO
1547	ENDDO
1548	ENDDO
1549	ELSEIF ( intermediate_timestep_count < &
1550	intermediate_timestep_count_max ) THEN
1551	DO i = nxl, nxr
1552	DO j = nys, nyn
1553	DO k = nzb_s_inner(j,i)+1, nzt
1554	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1555	ENDDO
1556	ENDDO
1557	ENDDO
1558	ENDIF
1559	ENDIF
1560
1561	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1562
1563	ENDIF
1564
1565
1566	END SUBROUTINE prognostic_equations_vector
1567
1568
1569	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |