Home

Context Navigation

prognostic_equations.f90 @ 132

Last change on this file since 132 was 129, checked in by letzel, 17 years ago

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

Bugfix: assignment of fluxes at walls

Updates to documentation:
chapter_4.2.html#mode_dvrp
chapter_3.5.4.html#time_series

Default for mrun options -q and -n is "sla3" for lctit. Queues bes1 and bes2
removed. DOC/misc/Tsubame.html updated.

Modified default paths (/work/...) for lctit in .mrun.config.default

Property svn:keywords set to Id

File size: 66.3 KB

Line
1	MODULE prognostic_equations_mod
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! add wall_*flux to parameter list of calls of diffusion_s
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: prognostic_equations.f90 129 2007-10-30 12:12:24Z letzel $
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, &
372	wall_heatflux, tend )
373	ELSE
374	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
375	tend(:,j,i) = 0.0
376	CALL advec_s_pw( i, j, pt )
377	ELSE
378	IF ( scalar_advec /= 'ups-scheme' ) THEN
379	tend(:,j,i) = 0.0
380	CALL advec_s_up( i, j, pt )
381	ENDIF
382	ENDIF
383	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
384	THEN
385	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
386	tswst_m, wall_heatflux, tend )
387	ELSE
388	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
389	wall_heatflux, tend )
390	ENDIF
391	ENDIF
392
393	!
394	!-- If required compute heating/cooling due to long wave radiation
395	!-- processes
396	IF ( radiation ) THEN
397	CALL calc_radiation( i, j )
398	ENDIF
399
400	!
401	!-- If required compute impact of latent heat due to precipitation
402	IF ( precipitation ) THEN
403	CALL impact_of_latent_heat( i, j )
404	ENDIF
405	CALL user_actions( i, j, 'pt-tendency' )
406
407	!
408	!-- Prognostic equation for potential temperature
409	DO k = nzb_s_inner(j,i)+1, nzt
410	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
411	dt_3d * ( &
412	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
413	) - &
414	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
415	ENDDO
416
417	!
418	!-- Calculate tendencies for the next Runge-Kutta step
419	IF ( timestep_scheme(1:5) == 'runge' ) THEN
420	IF ( intermediate_timestep_count == 1 ) THEN
421	DO k = nzb_s_inner(j,i)+1, nzt
422	tpt_m(k,j,i) = tend(k,j,i)
423	ENDDO
424	ELSEIF ( intermediate_timestep_count < &
425	intermediate_timestep_count_max ) THEN
426	DO k = nzb_s_inner(j,i)+1, nzt
427	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
428	ENDDO
429	ENDIF
430	ENDIF
431
432	ENDDO
433	ENDDO
434
435	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
436
437	!
438	!-- If required, compute prognostic equation for salinity
439	IF ( ocean ) THEN
440
441	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
442
443	!
444	!-- sa-tendency terms with communication
445	sat = tsc(1)
446	sbt = tsc(2)
447	IF ( scalar_advec == 'bc-scheme' ) THEN
448
449	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
450	!
451	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
452	!-- switched on. Thus:
453	sat = 1.0
454	sbt = 1.0
455	ENDIF
456	tend = 0.0
457	CALL advec_s_bc( sa, 'sa' )
458	ELSE
459	IF ( tsc(2) /= 2.0 ) THEN
460	IF ( scalar_advec == 'ups-scheme' ) THEN
461	tend = 0.0
462	CALL advec_s_ups( sa, 'sa' )
463	ENDIF
464	ENDIF
465	ENDIF
466
467	!
468	!-- sa terms with no communication
469	DO i = nxl, nxr
470	DO j = nys, nyn
471	!
472	!-- Tendency-terms
473	IF ( scalar_advec == 'bc-scheme' ) THEN
474	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
475	wall_salinityflux, tend )
476	ELSE
477	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
478	tend(:,j,i) = 0.0
479	CALL advec_s_pw( i, j, sa )
480	ELSE
481	IF ( scalar_advec /= 'ups-scheme' ) THEN
482	tend(:,j,i) = 0.0
483	CALL advec_s_up( i, j, sa )
484	ENDIF
485	ENDIF
486	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
487	wall_salinityflux, tend )
488	ENDIF
489
490	CALL user_actions( i, j, 'sa-tendency' )
491
492	!
493	!-- Prognostic equation for salinity
494	DO k = nzb_s_inner(j,i)+1, nzt
495	sa_p(k,j,i) = sat * sa(k,j,i) + &
496	dt_3d * ( &
497	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
498	) - &
499	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
500	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
501	ENDDO
502
503	!
504	!-- Calculate tendencies for the next Runge-Kutta step
505	IF ( timestep_scheme(1:5) == 'runge' ) THEN
506	IF ( intermediate_timestep_count == 1 ) THEN
507	DO k = nzb_s_inner(j,i)+1, nzt
508	tsa_m(k,j,i) = tend(k,j,i)
509	ENDDO
510	ELSEIF ( intermediate_timestep_count < &
511	intermediate_timestep_count_max ) THEN
512	DO k = nzb_s_inner(j,i)+1, nzt
513	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
514	5.3125 * tsa_m(k,j,i)
515	ENDDO
516	ENDIF
517	ENDIF
518
519	!
520	!-- Calculate density by the equation of state for seawater
521	CALL eqn_state_seawater( i, j )
522
523	ENDDO
524	ENDDO
525
526	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
527
528	ENDIF
529
530	!
531	!-- If required, compute prognostic equation for total water content / scalar
532	IF ( humidity .OR. passive_scalar ) THEN
533
534	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
535
536	!
537	!-- Scalar/q-tendency terms with communication
538	sat = tsc(1)
539	sbt = tsc(2)
540	IF ( scalar_advec == 'bc-scheme' ) THEN
541
542	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
543	!
544	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
545	!-- switched on. Thus:
546	sat = 1.0
547	sbt = 1.0
548	ENDIF
549	tend = 0.0
550	CALL advec_s_bc( q, 'q' )
551	ELSE
552	IF ( tsc(2) /= 2.0 ) THEN
553	IF ( scalar_advec == 'ups-scheme' ) THEN
554	tend = 0.0
555	CALL advec_s_ups( q, 'q' )
556	ENDIF
557	ENDIF
558	ENDIF
559
560	!
561	!-- Scalar/q-tendency terms with no communication
562	DO i = nxl, nxr
563	DO j = nys, nyn
564	!
565	!-- Tendency-terms
566	IF ( scalar_advec == 'bc-scheme' ) THEN
567	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
568	wall_qflux, tend )
569	ELSE
570	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
571	tend(:,j,i) = 0.0
572	CALL advec_s_pw( i, j, q )
573	ELSE
574	IF ( scalar_advec /= 'ups-scheme' ) THEN
575	tend(:,j,i) = 0.0
576	CALL advec_s_up( i, j, q )
577	ENDIF
578	ENDIF
579	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
580	THEN
581	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
582	qswst_m, wall_qflux, tend )
583	ELSE
584	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
585	wall_qflux, tend )
586	ENDIF
587	ENDIF
588
589	!
590	!-- If required compute decrease of total water content due to
591	!-- precipitation
592	IF ( precipitation ) THEN
593	CALL calc_precipitation( i, j )
594	ENDIF
595	CALL user_actions( i, j, 'q-tendency' )
596
597	!
598	!-- Prognostic equation for total water content / scalar
599	DO k = nzb_s_inner(j,i)+1, nzt
600	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
601	dt_3d * ( &
602	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
603	) - &
604	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
605	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
606	ENDDO
607
608	!
609	!-- Calculate tendencies for the next Runge-Kutta step
610	IF ( timestep_scheme(1:5) == 'runge' ) THEN
611	IF ( intermediate_timestep_count == 1 ) THEN
612	DO k = nzb_s_inner(j,i)+1, nzt
613	tq_m(k,j,i) = tend(k,j,i)
614	ENDDO
615	ELSEIF ( intermediate_timestep_count < &
616	intermediate_timestep_count_max ) THEN
617	DO k = nzb_s_inner(j,i)+1, nzt
618	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
619	ENDDO
620	ENDIF
621	ENDIF
622
623	ENDDO
624	ENDDO
625
626	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
627
628	ENDIF
629
630	!
631	!-- If required, compute prognostic equation for turbulent kinetic
632	!-- energy (TKE)
633	IF ( .NOT. constant_diffusion ) THEN
634
635	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
636
637	!
638	!-- TKE-tendency terms with communication
639	CALL production_e_init
640
641	sat = tsc(1)
642	sbt = tsc(2)
643	IF ( .NOT. use_upstream_for_tke ) THEN
644	IF ( scalar_advec == 'bc-scheme' ) THEN
645
646	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
647	!
648	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
649	!-- switched on. Thus:
650	sat = 1.0
651	sbt = 1.0
652	ENDIF
653	tend = 0.0
654	CALL advec_s_bc( e, 'e' )
655	ELSE
656	IF ( tsc(2) /= 2.0 ) THEN
657	IF ( scalar_advec == 'ups-scheme' ) THEN
658	tend = 0.0
659	CALL advec_s_ups( e, 'e' )
660	ENDIF
661	ENDIF
662	ENDIF
663	ENDIF
664
665	!
666	!-- TKE-tendency terms with no communication
667	DO i = nxl, nxr
668	DO j = nys, nyn
669	!
670	!-- Tendency-terms
671	IF ( scalar_advec == 'bc-scheme' .AND. &
672	.NOT. use_upstream_for_tke ) THEN
673	IF ( .NOT. humidity ) THEN
674	IF ( ocean ) THEN
675	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
676	l_grid, rho, prho_reference, rif, tend, &
677	zu, zw )
678	ELSE
679	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
680	l_grid, pt, pt_reference, rif, tend, &
681	zu, zw )
682	ENDIF
683	ELSE
684	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
685	l_grid, vpt, pt_reference, rif, tend, zu, &
686	zw )
687	ENDIF
688	ELSE
689	IF ( use_upstream_for_tke ) THEN
690	tend(:,j,i) = 0.0
691	CALL advec_s_up( i, j, e )
692	ELSE
693	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
694	THEN
695	tend(:,j,i) = 0.0
696	CALL advec_s_pw( i, j, e )
697	ELSE
698	IF ( scalar_advec /= 'ups-scheme' ) THEN
699	tend(:,j,i) = 0.0
700	CALL advec_s_up( i, j, e )
701	ENDIF
702	ENDIF
703	ENDIF
704	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
705	THEN
706	IF ( .NOT. humidity ) THEN
707	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
708	km_m, l_grid, pt_m, pt_reference, &
709	rif_m, tend, zu, zw )
710	ELSE
711	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
712	km_m, l_grid, vpt_m, pt_reference, &
713	rif_m, tend, zu, zw )
714	ENDIF
715	ELSE
716	IF ( .NOT. humidity ) THEN
717	IF ( ocean ) THEN
718	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
719	km, l_grid, rho, prho_reference, &
720	rif, tend, zu, zw )
721	ELSE
722	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
723	km, l_grid, pt, pt_reference, rif, &
724	tend, zu, zw )
725	ENDIF
726	ELSE
727	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
728	l_grid, vpt, pt_reference, rif, tend, &
729	zu, zw )
730	ENDIF
731	ENDIF
732	ENDIF
733	CALL production_e( i, j )
734	CALL user_actions( i, j, 'e-tendency' )
735
736	!
737	!-- Prognostic equation for TKE.
738	!-- Eliminate negative TKE values, which can occur due to numerical
739	!-- reasons in the course of the integration. In such cases the old TKE
740	!-- value is reduced by 90%.
741	DO k = nzb_s_inner(j,i)+1, nzt
742	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
743	dt_3d * ( &
744	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
745	)
746	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
747	ENDDO
748
749	!
750	!-- Calculate tendencies for the next Runge-Kutta step
751	IF ( timestep_scheme(1:5) == 'runge' ) THEN
752	IF ( intermediate_timestep_count == 1 ) THEN
753	DO k = nzb_s_inner(j,i)+1, nzt
754	te_m(k,j,i) = tend(k,j,i)
755	ENDDO
756	ELSEIF ( intermediate_timestep_count < &
757	intermediate_timestep_count_max ) THEN
758	DO k = nzb_s_inner(j,i)+1, nzt
759	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
760	ENDDO
761	ENDIF
762	ENDIF
763
764	ENDDO
765	ENDDO
766
767	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
768
769	ENDIF
770
771
772	END SUBROUTINE prognostic_equations_noopt
773
774
775	SUBROUTINE prognostic_equations_cache
776
777	!------------------------------------------------------------------------------!
778	! Version with one optimized loop over all equations. It is only allowed to
779	! be called for the standard Piascek-Williams advection scheme.
780	!
781	! Here the calls of most subroutines are embedded in two DO loops over i and j,
782	! so communication between CPUs is not allowed (does not make sense) within
783	! these loops.
784	!
785	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
786	!------------------------------------------------------------------------------!
787
788	IMPLICIT NONE
789
790	CHARACTER (LEN=9) :: time_to_string
791	INTEGER :: i, j, k
792
793
794	!
795	!-- Time measurement can only be performed for the whole set of equations
796	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
797
798
799	!
800	!-- Calculate those variables needed in the tendency terms which need
801	!-- global communication
802	CALL calc_mean_profile( pt, 4 )
803	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
804	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
805	IF ( .NOT. constant_diffusion ) CALL production_e_init
806
807
808	!
809	!-- Loop over all prognostic equations
810	!$OMP PARALLEL private (i,j,k)
811	!$OMP DO
812	DO i = nxl, nxr
813	DO j = nys, nyn
814	!
815	!-- Tendency terms for u-velocity component
816	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
817
818	tend(:,j,i) = 0.0
819	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
820	CALL advec_u_pw( i, j )
821	ELSE
822	CALL advec_u_up( i, j )
823	ENDIF
824	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
825	THEN
826	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
827	u_m, usws_m, uswst_m, v_m, w_m )
828	ELSE
829	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
830	usws, uswst, v, w )
831	ENDIF
832	CALL coriolis( i, j, 1 )
833	IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, &
834	4 )
835	CALL user_actions( i, j, 'u-tendency' )
836
837	!
838	!-- Prognostic equation for u-velocity component
839	DO k = nzb_u_inner(j,i)+1, nzt
840	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
841	dt_3d * ( &
842	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
843	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
844	) - &
845	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
846	ENDDO
847
848	!
849	!-- Calculate tendencies for the next Runge-Kutta step
850	IF ( timestep_scheme(1:5) == 'runge' ) THEN
851	IF ( intermediate_timestep_count == 1 ) THEN
852	DO k = nzb_u_inner(j,i)+1, nzt
853	tu_m(k,j,i) = tend(k,j,i)
854	ENDDO
855	ELSEIF ( intermediate_timestep_count < &
856	intermediate_timestep_count_max ) THEN
857	DO k = nzb_u_inner(j,i)+1, nzt
858	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
859	ENDDO
860	ENDIF
861	ENDIF
862
863	ENDIF
864
865	!
866	!-- Tendency terms for v-velocity component
867	IF ( .NOT. outflow_s .OR. j > nys ) THEN
868
869	tend(:,j,i) = 0.0
870	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
871	CALL advec_v_pw( i, j )
872	ELSE
873	CALL advec_v_up( i, j )
874	ENDIF
875	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
876	THEN
877	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
878	u_m, v_m, vsws_m, vswst_m, w_m )
879	ELSE
880	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
881	vsws, vswst, w )
882	ENDIF
883	CALL coriolis( i, j, 2 )
884	CALL user_actions( i, j, 'v-tendency' )
885
886	!
887	!-- Prognostic equation for v-velocity component
888	DO k = nzb_v_inner(j,i)+1, nzt
889	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
890	dt_3d * ( &
891	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
892	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
893	) - &
894	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
895	ENDDO
896
897	!
898	!-- Calculate tendencies for the next Runge-Kutta step
899	IF ( timestep_scheme(1:5) == 'runge' ) THEN
900	IF ( intermediate_timestep_count == 1 ) THEN
901	DO k = nzb_v_inner(j,i)+1, nzt
902	tv_m(k,j,i) = tend(k,j,i)
903	ENDDO
904	ELSEIF ( intermediate_timestep_count < &
905	intermediate_timestep_count_max ) THEN
906	DO k = nzb_v_inner(j,i)+1, nzt
907	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
908	ENDDO
909	ENDIF
910	ENDIF
911
912	ENDIF
913
914	!
915	!-- Tendency terms for w-velocity component
916	tend(:,j,i) = 0.0
917	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
918	CALL advec_w_pw( i, j )
919	ELSE
920	CALL advec_w_up( i, j )
921	ENDIF
922	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
923	THEN
924	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
925	km_damp_y, tend, u_m, v_m, w_m )
926	ELSE
927	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
928	tend, u, v, w )
929	ENDIF
930	CALL coriolis( i, j, 3 )
931	IF ( ocean ) THEN
932	CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
933	ELSE
934	IF ( .NOT. humidity ) THEN
935	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
936	ELSE
937	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
938	ENDIF
939	ENDIF
940	CALL user_actions( i, j, 'w-tendency' )
941
942	!
943	!-- Prognostic equation for w-velocity component
944	DO k = nzb_w_inner(j,i)+1, nzt-1
945	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
946	dt_3d * ( &
947	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
948	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
949	) - &
950	tsc(5) * rdf(k) * w(k,j,i)
951	ENDDO
952
953	!
954	!-- Calculate tendencies for the next Runge-Kutta step
955	IF ( timestep_scheme(1:5) == 'runge' ) THEN
956	IF ( intermediate_timestep_count == 1 ) THEN
957	DO k = nzb_w_inner(j,i)+1, nzt-1
958	tw_m(k,j,i) = tend(k,j,i)
959	ENDDO
960	ELSEIF ( intermediate_timestep_count < &
961	intermediate_timestep_count_max ) THEN
962	DO k = nzb_w_inner(j,i)+1, nzt-1
963	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
964	ENDDO
965	ENDIF
966	ENDIF
967
968	!
969	!-- Tendency terms for potential temperature
970	tend(:,j,i) = 0.0
971	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
972	CALL advec_s_pw( i, j, pt )
973	ELSE
974	CALL advec_s_up( i, j, pt )
975	ENDIF
976	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
977	THEN
978	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
979	tswst_m, wall_heatflux, tend )
980	ELSE
981	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
982	wall_heatflux, tend )
983	ENDIF
984
985	!
986	!-- If required compute heating/cooling due to long wave radiation
987	!-- processes
988	IF ( radiation ) THEN
989	CALL calc_radiation( i, j )
990	ENDIF
991
992	!
993	!-- If required compute impact of latent heat due to precipitation
994	IF ( precipitation ) THEN
995	CALL impact_of_latent_heat( i, j )
996	ENDIF
997	CALL user_actions( i, j, 'pt-tendency' )
998
999	!
1000	!-- Prognostic equation for potential temperature
1001	DO k = nzb_s_inner(j,i)+1, nzt
1002	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
1003	dt_3d * ( &
1004	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1005	) - &
1006	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1007	ENDDO
1008
1009	!
1010	!-- Calculate tendencies for the next Runge-Kutta step
1011	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1012	IF ( intermediate_timestep_count == 1 ) THEN
1013	DO k = nzb_s_inner(j,i)+1, nzt
1014	tpt_m(k,j,i) = tend(k,j,i)
1015	ENDDO
1016	ELSEIF ( intermediate_timestep_count < &
1017	intermediate_timestep_count_max ) THEN
1018	DO k = nzb_s_inner(j,i)+1, nzt
1019	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1020	5.3125 * tpt_m(k,j,i)
1021	ENDDO
1022	ENDIF
1023	ENDIF
1024
1025	!
1026	!-- If required, compute prognostic equation for salinity
1027	IF ( ocean ) THEN
1028
1029	!
1030	!-- Tendency-terms for salinity
1031	tend(:,j,i) = 0.0
1032	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1033	THEN
1034	CALL advec_s_pw( i, j, sa )
1035	ELSE
1036	CALL advec_s_up( i, j, sa )
1037	ENDIF
1038	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
1039	wall_salinityflux, tend )
1040
1041	CALL user_actions( i, j, 'sa-tendency' )
1042
1043	!
1044	!-- Prognostic equation for salinity
1045	DO k = nzb_s_inner(j,i)+1, nzt
1046	sa_p(k,j,i) = tsc(1) * sa(k,j,i) + &
1047	dt_3d * ( &
1048	tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1049	) - &
1050	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
1051	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1052	ENDDO
1053
1054	!
1055	!-- Calculate tendencies for the next Runge-Kutta step
1056	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1057	IF ( intermediate_timestep_count == 1 ) THEN
1058	DO k = nzb_s_inner(j,i)+1, nzt
1059	tsa_m(k,j,i) = tend(k,j,i)
1060	ENDDO
1061	ELSEIF ( intermediate_timestep_count < &
1062	intermediate_timestep_count_max ) THEN
1063	DO k = nzb_s_inner(j,i)+1, nzt
1064	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1065	5.3125 * tsa_m(k,j,i)
1066	ENDDO
1067	ENDIF
1068	ENDIF
1069
1070	!
1071	!-- Calculate density by the equation of state for seawater
1072	CALL eqn_state_seawater( i, j )
1073
1074	ENDIF
1075
1076	!
1077	!-- If required, compute prognostic equation for total water content /
1078	!-- scalar
1079	IF ( humidity .OR. passive_scalar ) THEN
1080
1081	!
1082	!-- Tendency-terms for total water content / scalar
1083	tend(:,j,i) = 0.0
1084	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1085	THEN
1086	CALL advec_s_pw( i, j, q )
1087	ELSE
1088	CALL advec_s_up( i, j, q )
1089	ENDIF
1090	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
1091	THEN
1092	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
1093	qswst_m, wall_qflux, tend )
1094	ELSE
1095	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
1096	wall_qflux, tend )
1097	ENDIF
1098
1099	!
1100	!-- If required compute decrease of total water content due to
1101	!-- precipitation
1102	IF ( precipitation ) THEN
1103	CALL calc_precipitation( i, j )
1104	ENDIF
1105	CALL user_actions( i, j, 'q-tendency' )
1106
1107	!
1108	!-- Prognostic equation for total water content / scalar
1109	DO k = nzb_s_inner(j,i)+1, nzt
1110	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
1111	dt_3d * ( &
1112	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1113	) - &
1114	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1115	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1116	ENDDO
1117
1118	!
1119	!-- Calculate tendencies for the next Runge-Kutta step
1120	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1121	IF ( intermediate_timestep_count == 1 ) THEN
1122	DO k = nzb_s_inner(j,i)+1, nzt
1123	tq_m(k,j,i) = tend(k,j,i)
1124	ENDDO
1125	ELSEIF ( intermediate_timestep_count < &
1126	intermediate_timestep_count_max ) THEN
1127	DO k = nzb_s_inner(j,i)+1, nzt
1128	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1129	5.3125 * tq_m(k,j,i)
1130	ENDDO
1131	ENDIF
1132	ENDIF
1133
1134	ENDIF
1135
1136	!
1137	!-- If required, compute prognostic equation for turbulent kinetic
1138	!-- energy (TKE)
1139	IF ( .NOT. constant_diffusion ) THEN
1140
1141	!
1142	!-- Tendency-terms for TKE
1143	tend(:,j,i) = 0.0
1144	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1145	.AND. .NOT. use_upstream_for_tke ) THEN
1146	CALL advec_s_pw( i, j, e )
1147	ELSE
1148	CALL advec_s_up( i, j, e )
1149	ENDIF
1150	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
1151	THEN
1152	IF ( .NOT. humidity ) THEN
1153	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1154	km_m, l_grid, pt_m, pt_reference, &
1155	rif_m, tend, zu, zw )
1156	ELSE
1157	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1158	km_m, l_grid, vpt_m, pt_reference, &
1159	rif_m, tend, zu, zw )
1160	ENDIF
1161	ELSE
1162	IF ( .NOT. humidity ) THEN
1163	IF ( ocean ) THEN
1164	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
1165	km, l_grid, rho, prho_reference, &
1166	rif, tend, zu, zw )
1167	ELSE
1168	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
1169	km, l_grid, pt, pt_reference, rif, &
1170	tend, zu, zw )
1171	ENDIF
1172	ELSE
1173	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
1174	l_grid, vpt, pt_reference, rif, tend, &
1175	zu, zw )
1176	ENDIF
1177	ENDIF
1178	CALL production_e( i, j )
1179	CALL user_actions( i, j, 'e-tendency' )
1180
1181	!
1182	!-- Prognostic equation for TKE.
1183	!-- Eliminate negative TKE values, which can occur due to numerical
1184	!-- reasons in the course of the integration. In such cases the old
1185	!-- TKE value is reduced by 90%.
1186	DO k = nzb_s_inner(j,i)+1, nzt
1187	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
1188	dt_3d * ( &
1189	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1190	)
1191	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1192	ENDDO
1193
1194	!
1195	!-- Calculate tendencies for the next Runge-Kutta step
1196	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1197	IF ( intermediate_timestep_count == 1 ) THEN
1198	DO k = nzb_s_inner(j,i)+1, nzt
1199	te_m(k,j,i) = tend(k,j,i)
1200	ENDDO
1201	ELSEIF ( intermediate_timestep_count < &
1202	intermediate_timestep_count_max ) THEN
1203	DO k = nzb_s_inner(j,i)+1, nzt
1204	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1205	5.3125 * te_m(k,j,i)
1206	ENDDO
1207	ENDIF
1208	ENDIF
1209
1210	ENDIF ! TKE equation
1211
1212	ENDDO
1213	ENDDO
1214	!$OMP END PARALLEL
1215
1216	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1217
1218
1219	END SUBROUTINE prognostic_equations_cache
1220
1221
1222	SUBROUTINE prognostic_equations_vector
1223
1224	!------------------------------------------------------------------------------!
1225	! Version for vector machines
1226	!------------------------------------------------------------------------------!
1227
1228	IMPLICIT NONE
1229
1230	CHARACTER (LEN=9) :: time_to_string
1231	INTEGER :: i, j, k
1232	REAL :: sat, sbt
1233
1234	!
1235	!-- Calculate those variables needed in the tendency terms which need
1236	!-- global communication
1237	CALL calc_mean_profile( pt, 4 )
1238	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
1239	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
1240
1241	!
1242	!-- u-velocity component
1243	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1244
1245	!
1246	!-- u-tendency terms with communication
1247	IF ( momentum_advec == 'ups-scheme' ) THEN
1248	tend = 0.0
1249	CALL advec_u_ups
1250	ENDIF
1251
1252	!
1253	!-- u-tendency terms with no communication
1254	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1255	tend = 0.0
1256	CALL advec_u_pw
1257	ELSE
1258	IF ( momentum_advec /= 'ups-scheme' ) THEN
1259	tend = 0.0
1260	CALL advec_u_up
1261	ENDIF
1262	ENDIF
1263	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1264	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
1265	uswst_m, v_m, w_m )
1266	ELSE
1267	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
1268	ENDIF
1269	CALL coriolis( 1 )
1270	IF ( sloping_surface ) CALL buoyancy( pt, pt_reference, 1, 4 )
1271	CALL user_actions( 'u-tendency' )
1272
1273	!
1274	!-- Prognostic equation for u-velocity component
1275	DO i = nxlu, nxr
1276	DO j = nys, nyn
1277	DO k = nzb_u_inner(j,i)+1, nzt
1278	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1279	dt_3d * ( &
1280	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1281	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
1282	) - &
1283	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1284	ENDDO
1285	ENDDO
1286	ENDDO
1287
1288	!
1289	!-- Calculate tendencies for the next Runge-Kutta step
1290	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1291	IF ( intermediate_timestep_count == 1 ) THEN
1292	DO i = nxlu, nxr
1293	DO j = nys, nyn
1294	DO k = nzb_u_inner(j,i)+1, nzt
1295	tu_m(k,j,i) = tend(k,j,i)
1296	ENDDO
1297	ENDDO
1298	ENDDO
1299	ELSEIF ( intermediate_timestep_count < &
1300	intermediate_timestep_count_max ) THEN
1301	DO i = nxlu, nxr
1302	DO j = nys, nyn
1303	DO k = nzb_u_inner(j,i)+1, nzt
1304	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1305	ENDDO
1306	ENDDO
1307	ENDDO
1308	ENDIF
1309	ENDIF
1310
1311	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1312
1313	!
1314	!-- v-velocity component
1315	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1316
1317	!
1318	!-- v-tendency terms with communication
1319	IF ( momentum_advec == 'ups-scheme' ) THEN
1320	tend = 0.0
1321	CALL advec_v_ups
1322	ENDIF
1323
1324	!
1325	!-- v-tendency terms with no communication
1326	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1327	tend = 0.0
1328	CALL advec_v_pw
1329	ELSE
1330	IF ( momentum_advec /= 'ups-scheme' ) THEN
1331	tend = 0.0
1332	CALL advec_v_up
1333	ENDIF
1334	ENDIF
1335	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1336	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1337	vswst_m, w_m )
1338	ELSE
1339	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
1340	ENDIF
1341	CALL coriolis( 2 )
1342	CALL user_actions( 'v-tendency' )
1343
1344	!
1345	!-- Prognostic equation for v-velocity component
1346	DO i = nxl, nxr
1347	DO j = nysv, nyn
1348	DO k = nzb_v_inner(j,i)+1, nzt
1349	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1350	dt_3d * ( &
1351	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1352	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
1353	) - &
1354	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1355	ENDDO
1356	ENDDO
1357	ENDDO
1358
1359	!
1360	!-- Calculate tendencies for the next Runge-Kutta step
1361	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1362	IF ( intermediate_timestep_count == 1 ) THEN
1363	DO i = nxl, nxr
1364	DO j = nysv, nyn
1365	DO k = nzb_v_inner(j,i)+1, nzt
1366	tv_m(k,j,i) = tend(k,j,i)
1367	ENDDO
1368	ENDDO
1369	ENDDO
1370	ELSEIF ( intermediate_timestep_count < &
1371	intermediate_timestep_count_max ) THEN
1372	DO i = nxl, nxr
1373	DO j = nysv, nyn
1374	DO k = nzb_v_inner(j,i)+1, nzt
1375	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1376	ENDDO
1377	ENDDO
1378	ENDDO
1379	ENDIF
1380	ENDIF
1381
1382	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1383
1384	!
1385	!-- w-velocity component
1386	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1387
1388	!
1389	!-- w-tendency terms with communication
1390	IF ( momentum_advec == 'ups-scheme' ) THEN
1391	tend = 0.0
1392	CALL advec_w_ups
1393	ENDIF
1394
1395	!
1396	!-- w-tendency terms with no communication
1397	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1398	tend = 0.0
1399	CALL advec_w_pw
1400	ELSE
1401	IF ( momentum_advec /= 'ups-scheme' ) THEN
1402	tend = 0.0
1403	CALL advec_w_up
1404	ENDIF
1405	ENDIF
1406	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1407	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1408	v_m, w_m )
1409	ELSE
1410	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1411	ENDIF
1412	CALL coriolis( 3 )
1413	IF ( ocean ) THEN
1414	CALL buoyancy( rho, prho_reference, 3, 64 )
1415	ELSE
1416	IF ( .NOT. humidity ) THEN
1417	CALL buoyancy( pt, pt_reference, 3, 4 )
1418	ELSE
1419	CALL buoyancy( vpt, pt_reference, 3, 44 )
1420	ENDIF
1421	ENDIF
1422	CALL user_actions( 'w-tendency' )
1423
1424	!
1425	!-- Prognostic equation for w-velocity component
1426	DO i = nxl, nxr
1427	DO j = nys, nyn
1428	DO k = nzb_w_inner(j,i)+1, nzt-1
1429	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1430	dt_3d * ( &
1431	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1432	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1433	) - &
1434	tsc(5) * rdf(k) * w(k,j,i)
1435	ENDDO
1436	ENDDO
1437	ENDDO
1438
1439	!
1440	!-- Calculate tendencies for the next Runge-Kutta step
1441	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1442	IF ( intermediate_timestep_count == 1 ) THEN
1443	DO i = nxl, nxr
1444	DO j = nys, nyn
1445	DO k = nzb_w_inner(j,i)+1, nzt-1
1446	tw_m(k,j,i) = tend(k,j,i)
1447	ENDDO
1448	ENDDO
1449	ENDDO
1450	ELSEIF ( intermediate_timestep_count < &
1451	intermediate_timestep_count_max ) THEN
1452	DO i = nxl, nxr
1453	DO j = nys, nyn
1454	DO k = nzb_w_inner(j,i)+1, nzt-1
1455	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1456	ENDDO
1457	ENDDO
1458	ENDDO
1459	ENDIF
1460	ENDIF
1461
1462	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1463
1464	!
1465	!-- potential temperature
1466	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1467
1468	!
1469	!-- pt-tendency terms with communication
1470	sat = tsc(1)
1471	sbt = tsc(2)
1472	IF ( scalar_advec == 'bc-scheme' ) THEN
1473
1474	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1475	!
1476	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1477	!-- switched on. Thus:
1478	sat = 1.0
1479	sbt = 1.0
1480	ENDIF
1481	tend = 0.0
1482	CALL advec_s_bc( pt, 'pt' )
1483	ELSE
1484	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
1485	tend = 0.0
1486	CALL advec_s_ups( pt, 'pt' )
1487	ENDIF
1488	ENDIF
1489
1490	!
1491	!-- pt-tendency terms with no communication
1492	IF ( scalar_advec == 'bc-scheme' ) THEN
1493	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1494	tend )
1495	ELSE
1496	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1497	tend = 0.0
1498	CALL advec_s_pw( pt )
1499	ELSE
1500	IF ( scalar_advec /= 'ups-scheme' ) THEN
1501	tend = 0.0
1502	CALL advec_s_up( pt )
1503	ENDIF
1504	ENDIF
1505	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1506	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
1507	wall_heatflux, tend )
1508	ELSE
1509	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1510	tend )
1511	ENDIF
1512	ENDIF
1513
1514	!
1515	!-- If required compute heating/cooling due to long wave radiation
1516	!-- processes
1517	IF ( radiation ) THEN
1518	CALL calc_radiation
1519	ENDIF
1520
1521	!
1522	!-- If required compute impact of latent heat due to precipitation
1523	IF ( precipitation ) THEN
1524	CALL impact_of_latent_heat
1525	ENDIF
1526	CALL user_actions( 'pt-tendency' )
1527
1528	!
1529	!-- Prognostic equation for potential temperature
1530	DO i = nxl, nxr
1531	DO j = nys, nyn
1532	DO k = nzb_s_inner(j,i)+1, nzt
1533	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
1534	dt_3d * ( &
1535	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1536	) - &
1537	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1538	ENDDO
1539	ENDDO
1540	ENDDO
1541
1542	!
1543	!-- Calculate tendencies for the next Runge-Kutta step
1544	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1545	IF ( intermediate_timestep_count == 1 ) THEN
1546	DO i = nxl, nxr
1547	DO j = nys, nyn
1548	DO k = nzb_s_inner(j,i)+1, nzt
1549	tpt_m(k,j,i) = tend(k,j,i)
1550	ENDDO
1551	ENDDO
1552	ENDDO
1553	ELSEIF ( intermediate_timestep_count < &
1554	intermediate_timestep_count_max ) THEN
1555	DO i = nxl, nxr
1556	DO j = nys, nyn
1557	DO k = nzb_s_inner(j,i)+1, nzt
1558	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1559	ENDDO
1560	ENDDO
1561	ENDDO
1562	ENDIF
1563	ENDIF
1564
1565	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1566
1567	!
1568	!-- If required, compute prognostic equation for salinity
1569	IF ( ocean ) THEN
1570
1571	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1572
1573	!
1574	!-- sa-tendency terms with communication
1575	sat = tsc(1)
1576	sbt = tsc(2)
1577	IF ( scalar_advec == 'bc-scheme' ) THEN
1578
1579	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1580	!
1581	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1582	!-- switched on. Thus:
1583	sat = 1.0
1584	sbt = 1.0
1585	ENDIF
1586	tend = 0.0
1587	CALL advec_s_bc( sa, 'sa' )
1588	ELSE
1589	IF ( tsc(2) /= 2.0 ) THEN
1590	IF ( scalar_advec == 'ups-scheme' ) THEN
1591	tend = 0.0
1592	CALL advec_s_ups( sa, 'sa' )
1593	ENDIF
1594	ENDIF
1595	ENDIF
1596
1597	!
1598	!-- sa-tendency terms with no communication
1599	IF ( scalar_advec == 'bc-scheme' ) THEN
1600	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1601	wall_salinityflux, tend )
1602	ELSE
1603	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1604	tend = 0.0
1605	CALL advec_s_pw( sa )
1606	ELSE
1607	IF ( scalar_advec /= 'ups-scheme' ) THEN
1608	tend = 0.0
1609	CALL advec_s_up( sa )
1610	ENDIF
1611	ENDIF
1612	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1613	wall_salinityflux, tend )
1614	ENDIF
1615
1616	CALL user_actions( 'sa-tendency' )
1617
1618	!
1619	!-- Prognostic equation for salinity
1620	DO i = nxl, nxr
1621	DO j = nys, nyn
1622	DO k = nzb_s_inner(j,i)+1, nzt
1623	sa_p(k,j,i) = sat * sa(k,j,i) + &
1624	dt_3d * ( &
1625	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1626	) - &
1627	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
1628	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1629	ENDDO
1630	ENDDO
1631	ENDDO
1632
1633	!
1634	!-- Calculate tendencies for the next Runge-Kutta step
1635	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1636	IF ( intermediate_timestep_count == 1 ) 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) = tend(k,j,i)
1641	ENDDO
1642	ENDDO
1643	ENDDO
1644	ELSEIF ( intermediate_timestep_count < &
1645	intermediate_timestep_count_max ) THEN
1646	DO i = nxl, nxr
1647	DO j = nys, nyn
1648	DO k = nzb_s_inner(j,i)+1, nzt
1649	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1650	5.3125 * tsa_m(k,j,i)
1651	ENDDO
1652	ENDDO
1653	ENDDO
1654	ENDIF
1655	ENDIF
1656
1657	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1658
1659	!
1660	!-- Calculate density by the equation of state for seawater
1661	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1662	CALL eqn_state_seawater
1663	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1664
1665	ENDIF
1666
1667	!
1668	!-- If required, compute prognostic equation for total water content / scalar
1669	IF ( humidity .OR. passive_scalar ) THEN
1670
1671	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1672
1673	!
1674	!-- Scalar/q-tendency terms with communication
1675	sat = tsc(1)
1676	sbt = tsc(2)
1677	IF ( scalar_advec == 'bc-scheme' ) THEN
1678
1679	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1680	!
1681	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1682	!-- switched on. Thus:
1683	sat = 1.0
1684	sbt = 1.0
1685	ENDIF
1686	tend = 0.0
1687	CALL advec_s_bc( q, 'q' )
1688	ELSE
1689	IF ( tsc(2) /= 2.0 ) THEN
1690	IF ( scalar_advec == 'ups-scheme' ) THEN
1691	tend = 0.0
1692	CALL advec_s_ups( q, 'q' )
1693	ENDIF
1694	ENDIF
1695	ENDIF
1696
1697	!
1698	!-- Scalar/q-tendency terms with no communication
1699	IF ( scalar_advec == 'bc-scheme' ) THEN
1700	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
1701	ELSE
1702	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1703	tend = 0.0
1704	CALL advec_s_pw( q )
1705	ELSE
1706	IF ( scalar_advec /= 'ups-scheme' ) THEN
1707	tend = 0.0
1708	CALL advec_s_up( q )
1709	ENDIF
1710	ENDIF
1711	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1712	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
1713	wall_qflux, tend )
1714	ELSE
1715	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
1716	wall_qflux, tend )
1717	ENDIF
1718	ENDIF
1719
1720	!
1721	!-- If required compute decrease of total water content due to
1722	!-- precipitation
1723	IF ( precipitation ) THEN
1724	CALL calc_precipitation
1725	ENDIF
1726	CALL user_actions( 'q-tendency' )
1727
1728	!
1729	!-- Prognostic equation for total water content / scalar
1730	DO i = nxl, nxr
1731	DO j = nys, nyn
1732	DO k = nzb_s_inner(j,i)+1, nzt
1733	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
1734	dt_3d * ( &
1735	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1736	) - &
1737	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1738	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1739	ENDDO
1740	ENDDO
1741	ENDDO
1742
1743	!
1744	!-- Calculate tendencies for the next Runge-Kutta step
1745	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1746	IF ( intermediate_timestep_count == 1 ) THEN
1747	DO i = nxl, nxr
1748	DO j = nys, nyn
1749	DO k = nzb_s_inner(j,i)+1, nzt
1750	tq_m(k,j,i) = tend(k,j,i)
1751	ENDDO
1752	ENDDO
1753	ENDDO
1754	ELSEIF ( intermediate_timestep_count < &
1755	intermediate_timestep_count_max ) THEN
1756	DO i = nxl, nxr
1757	DO j = nys, nyn
1758	DO k = nzb_s_inner(j,i)+1, nzt
1759	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1760	ENDDO
1761	ENDDO
1762	ENDDO
1763	ENDIF
1764	ENDIF
1765
1766	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1767
1768	ENDIF
1769
1770	!
1771	!-- If required, compute prognostic equation for turbulent kinetic
1772	!-- energy (TKE)
1773	IF ( .NOT. constant_diffusion ) THEN
1774
1775	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1776
1777	!
1778	!-- TKE-tendency terms with communication
1779	CALL production_e_init
1780
1781	sat = tsc(1)
1782	sbt = tsc(2)
1783	IF ( .NOT. use_upstream_for_tke ) THEN
1784	IF ( scalar_advec == 'bc-scheme' ) THEN
1785
1786	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1787	!
1788	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1789	!-- switched on. Thus:
1790	sat = 1.0
1791	sbt = 1.0
1792	ENDIF
1793	tend = 0.0
1794	CALL advec_s_bc( e, 'e' )
1795	ELSE
1796	IF ( tsc(2) /= 2.0 ) THEN
1797	IF ( scalar_advec == 'ups-scheme' ) THEN
1798	tend = 0.0
1799	CALL advec_s_ups( e, 'e' )
1800	ENDIF
1801	ENDIF
1802	ENDIF
1803	ENDIF
1804
1805	!
1806	!-- TKE-tendency terms with no communication
1807	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
1808	THEN
1809	IF ( .NOT. humidity ) THEN
1810	IF ( ocean ) THEN
1811	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
1812	prho_reference, rif, tend, zu, zw )
1813	ELSE
1814	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
1815	pt_reference, rif, tend, zu, zw )
1816	ENDIF
1817	ELSE
1818	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1819	pt_reference, rif, tend, zu, zw )
1820	ENDIF
1821	ELSE
1822	IF ( use_upstream_for_tke ) THEN
1823	tend = 0.0
1824	CALL advec_s_up( e )
1825	ELSE
1826	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1827	tend = 0.0
1828	CALL advec_s_pw( e )
1829	ELSE
1830	IF ( scalar_advec /= 'ups-scheme' ) THEN
1831	tend = 0.0
1832	CALL advec_s_up( e )
1833	ENDIF
1834	ENDIF
1835	ENDIF
1836	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1837	IF ( .NOT. humidity ) THEN
1838	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1839	pt_m, pt_reference, rif_m, tend, zu, zw )
1840	ELSE
1841	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
1842	vpt_m, pt_reference, rif_m, tend, zu, zw )
1843	ENDIF
1844	ELSE
1845	IF ( .NOT. humidity ) THEN
1846	IF ( ocean ) THEN
1847	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
1848	rho, prho_reference, rif, tend, zu, zw )
1849	ELSE
1850	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
1851	pt, pt_reference, rif, tend, zu, zw )
1852	ENDIF
1853	ELSE
1854	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
1855	pt_reference, rif, tend, zu, zw )
1856	ENDIF
1857	ENDIF
1858	ENDIF
1859	CALL production_e
1860	CALL user_actions( 'e-tendency' )
1861
1862	!
1863	!-- Prognostic equation for TKE.
1864	!-- Eliminate negative TKE values, which can occur due to numerical
1865	!-- reasons in the course of the integration. In such cases the old TKE
1866	!-- value is reduced by 90%.
1867	DO i = nxl, nxr
1868	DO j = nys, nyn
1869	DO k = nzb_s_inner(j,i)+1, nzt
1870	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
1871	dt_3d * ( &
1872	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1873	)
1874	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1875	ENDDO
1876	ENDDO
1877	ENDDO
1878
1879	!
1880	!-- Calculate tendencies for the next Runge-Kutta step
1881	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1882	IF ( intermediate_timestep_count == 1 ) THEN
1883	DO i = nxl, nxr
1884	DO j = nys, nyn
1885	DO k = nzb_s_inner(j,i)+1, nzt
1886	te_m(k,j,i) = tend(k,j,i)
1887	ENDDO
1888	ENDDO
1889	ENDDO
1890	ELSEIF ( intermediate_timestep_count < &
1891	intermediate_timestep_count_max ) THEN
1892	DO i = nxl, nxr
1893	DO j = nys, nyn
1894	DO k = nzb_s_inner(j,i)+1, nzt
1895	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1896	ENDDO
1897	ENDDO
1898	ENDDO
1899	ENDIF
1900	ENDIF
1901
1902	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1903
1904	ENDIF
1905
1906
1907	END SUBROUTINE prognostic_equations_vector
1908
1909
1910	END MODULE prognostic_equations_mod

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

Context Navigation

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

Download in other formats: