Home

Context Navigation

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

Last change on this file since 505 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 71.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |