Home

Context Navigation

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

Last change on this file since 106 was 106, checked in by raasch, 17 years ago
preliminary update of bugfixes and extensions for non-cyclic BCs
Property svn:keywords set to `Id`
File size: 65.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |