Home

Context Navigation

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

Last change on this file since 28 was 19, checked in by raasch, 17 years ago
preliminary version of modified boundary conditions at top
Property svn:keywords set to `Id`
File size: 53.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |