Home

Context Navigation

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

Last change on this file since 95 was 95, checked in by raasch, 17 years ago
further preliminary uncomplete changes for ocean version
Property svn:keywords set to `Id`
File size: 62.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |