Home

Context Navigation

prognostic_equations.f90 @ 97

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

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:

inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:

Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

Property svn:keywords set to Id

File size: 66.0 KB

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

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

Context Navigation

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

Download in other formats: