Home

Context Navigation

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

Last change on this file since 2 was 1, checked in by raasch, 17 years ago
Initial repository layout and content
File size: 56.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |