Home

Context Navigation

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

Last change on this file since 16 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 53.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |