Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |