Home

Context Navigation

prognostic_equations.f90 @ 113

Last change on this file since 113 was 110, checked in by raasch, 17 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:

Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:

Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 65.6 KB

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

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

Context Navigation

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

Download in other formats: