Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |