Home

Context Navigation

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

Last change on this file since 138 was 138, checked in by letzel, 16 years ago
Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
Property svn:keywords set to `Id`
File size: 67.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |