Home

Context Navigation

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

Last change on this file since 558 was 532, checked in by heinze, 15 years ago
NCL scripts allow the output of png files
Property svn:keywords set to `Id`
File size: 71.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |