Home

Context Navigation

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

Last change on this file since 100 was 98, checked in by raasch, 17 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 66.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |