Home

Context Navigation

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

Last change on this file since 531 was 531, checked in by heinze, 14 years ago
call of subsidence added in prognostic equations for humidity/passive scalar, some bugfixes
Property svn:keywords set to `Id`
File size: 71.6 KB

Line
1	MODULE prognostic_equations_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! add call of subsidence in the equation for humidity / passive scalar
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: prognostic_equations.f90 531 2010-04-21 06:47:21Z heinze $
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	!
662	!-- If required compute influence of large-scale subsidence/ascent
663	IF ( large_scale_subsidence ) THEN
664	CALL subsidence ( i, j, tend, q, q_init )
665	ENDIF
666
667	CALL user_actions( i, j, 'q-tendency' )
668
669	!
670	!-- Prognostic equation for total water content / scalar
671	DO k = nzb_s_inner(j,i)+1, nzt
672	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
673	dt_3d * ( &
674	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
675	) - &
676	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
677	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
678	ENDDO
679
680	!
681	!-- Calculate tendencies for the next Runge-Kutta step
682	IF ( timestep_scheme(1:5) == 'runge' ) THEN
683	IF ( intermediate_timestep_count == 1 ) THEN
684	DO k = nzb_s_inner(j,i)+1, nzt
685	tq_m(k,j,i) = tend(k,j,i)
686	ENDDO
687	ELSEIF ( intermediate_timestep_count < &
688	intermediate_timestep_count_max ) THEN
689	DO k = nzb_s_inner(j,i)+1, nzt
690	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
691	ENDDO
692	ENDIF
693	ENDIF
694
695	ENDDO
696	ENDDO
697
698	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
699
700	ENDIF
701
702	!
703	!-- If required, compute prognostic equation for turbulent kinetic
704	!-- energy (TKE)
705	IF ( .NOT. constant_diffusion ) THEN
706
707	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
708
709	!
710	!-- TKE-tendency terms with communication
711	CALL production_e_init
712
713	sat = tsc(1)
714	sbt = tsc(2)
715	IF ( .NOT. use_upstream_for_tke ) THEN
716	IF ( scalar_advec == 'bc-scheme' ) THEN
717
718	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
719	!
720	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
721	!-- switched on. Thus:
722	sat = 1.0
723	sbt = 1.0
724	ENDIF
725	tend = 0.0
726	CALL advec_s_bc( e, 'e' )
727	ELSE
728	IF ( tsc(2) /= 2.0 ) THEN
729	IF ( scalar_advec == 'ups-scheme' ) THEN
730	tend = 0.0
731	CALL advec_s_ups( e, 'e' )
732	ENDIF
733	ENDIF
734	ENDIF
735	ENDIF
736
737	!
738	!-- TKE-tendency terms with no communication
739	DO i = nxl, nxr
740	DO j = nys, nyn
741	!
742	!-- Tendency-terms
743	IF ( scalar_advec == 'bc-scheme' .AND. &
744	.NOT. use_upstream_for_tke ) THEN
745	IF ( .NOT. humidity ) THEN
746	IF ( ocean ) THEN
747	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
748	l_grid, prho, prho_reference, rif, &
749	tend, zu, zw )
750	ELSE
751	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
752	l_grid, pt, pt_reference, rif, tend, &
753	zu, zw )
754	ENDIF
755	ELSE
756	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
757	l_grid, vpt, pt_reference, rif, tend, zu, &
758	zw )
759	ENDIF
760	ELSE
761	IF ( use_upstream_for_tke ) THEN
762	tend(:,j,i) = 0.0
763	CALL advec_s_up( i, j, e )
764	ELSE
765	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
766	THEN
767	tend(:,j,i) = 0.0
768	CALL advec_s_pw( i, j, e )
769	ELSE
770	IF ( scalar_advec /= 'ups-scheme' ) THEN
771	tend(:,j,i) = 0.0
772	CALL advec_s_up( i, j, e )
773	ENDIF
774	ENDIF
775	ENDIF
776	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
777	THEN
778	IF ( .NOT. humidity ) THEN
779	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
780	km_m, l_grid, pt_m, pt_reference, &
781	rif_m, tend, zu, zw )
782	ELSE
783	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
784	km_m, l_grid, vpt_m, pt_reference, &
785	rif_m, tend, zu, zw )
786	ENDIF
787	ELSE
788	IF ( .NOT. humidity ) THEN
789	IF ( ocean ) THEN
790	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
791	km, l_grid, prho, prho_reference, &
792	rif, tend, zu, zw )
793	ELSE
794	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
795	km, l_grid, pt, pt_reference, rif, &
796	tend, zu, zw )
797	ENDIF
798	ELSE
799	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
800	l_grid, vpt, pt_reference, rif, tend, &
801	zu, zw )
802	ENDIF
803	ENDIF
804	ENDIF
805	CALL production_e( i, j )
806
807	!
808	!-- Additional sink term for flows through plant canopies
809	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
810
811	CALL user_actions( i, j, 'e-tendency' )
812
813	!
814	!-- Prognostic equation for TKE.
815	!-- Eliminate negative TKE values, which can occur due to numerical
816	!-- reasons in the course of the integration. In such cases the old TKE
817	!-- value is reduced by 90%.
818	DO k = nzb_s_inner(j,i)+1, nzt
819	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
820	dt_3d * ( &
821	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
822	)
823	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
824	ENDDO
825
826	!
827	!-- Calculate tendencies for the next Runge-Kutta step
828	IF ( timestep_scheme(1:5) == 'runge' ) THEN
829	IF ( intermediate_timestep_count == 1 ) THEN
830	DO k = nzb_s_inner(j,i)+1, nzt
831	te_m(k,j,i) = tend(k,j,i)
832	ENDDO
833	ELSEIF ( intermediate_timestep_count < &
834	intermediate_timestep_count_max ) THEN
835	DO k = nzb_s_inner(j,i)+1, nzt
836	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
837	ENDDO
838	ENDIF
839	ENDIF
840
841	ENDDO
842	ENDDO
843
844	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
845
846	ENDIF
847
848
849	END SUBROUTINE prognostic_equations_noopt
850
851
852	SUBROUTINE prognostic_equations_cache
853
854	!------------------------------------------------------------------------------!
855	! Version with one optimized loop over all equations. It is only allowed to
856	! be called for the standard Piascek-Williams advection scheme.
857	!
858	! Here the calls of most subroutines are embedded in two DO loops over i and j,
859	! so communication between CPUs is not allowed (does not make sense) within
860	! these loops.
861	!
862	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
863	!------------------------------------------------------------------------------!
864
865	IMPLICIT NONE
866
867	CHARACTER (LEN=9) :: time_to_string
868	INTEGER :: i, j, k
869
870
871	!
872	!-- Time measurement can only be performed for the whole set of equations
873	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
874
875
876	!
877	!-- Calculate those variables needed in the tendency terms which need
878	!-- global communication
879	CALL calc_mean_profile( pt, 4 )
880	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
881	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
882	IF ( .NOT. constant_diffusion ) CALL production_e_init
883
884
885	!
886	!-- Loop over all prognostic equations
887	!$OMP PARALLEL private (i,j,k)
888	!$OMP DO
889	DO i = nxl, nxr
890	DO j = nys, nyn
891	!
892	!-- Tendency terms for u-velocity component
893	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
894
895	tend(:,j,i) = 0.0
896	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
897	CALL advec_u_pw( i, j )
898	ELSE
899	CALL advec_u_up( i, j )
900	ENDIF
901	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
902	THEN
903	CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, &
904	u_m, usws_m, uswst_m, v_m, w_m )
905	ELSE
906	CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
907	usws, uswst, v, w )
908	ENDIF
909	CALL coriolis( i, j, 1 )
910	IF ( sloping_surface ) CALL buoyancy( i, j, pt, pt_reference, 1, &
911	4 )
912
913	!
914	!-- Drag by plant canopy
915	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
916
917	!
918	!-- External pressure gradient
919	IF ( dp_external ) THEN
920	DO k = dp_level_ind_b+1, nzt
921	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
922	ENDDO
923	ENDIF
924
925	CALL user_actions( i, j, 'u-tendency' )
926
927	!
928	!-- Prognostic equation for u-velocity component
929	DO k = nzb_u_inner(j,i)+1, nzt
930	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
931	dt_3d * ( &
932	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
933	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
934	) - &
935	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
936	ENDDO
937
938	!
939	!-- Calculate tendencies for the next Runge-Kutta step
940	IF ( timestep_scheme(1:5) == 'runge' ) THEN
941	IF ( intermediate_timestep_count == 1 ) THEN
942	DO k = nzb_u_inner(j,i)+1, nzt
943	tu_m(k,j,i) = tend(k,j,i)
944	ENDDO
945	ELSEIF ( intermediate_timestep_count < &
946	intermediate_timestep_count_max ) THEN
947	DO k = nzb_u_inner(j,i)+1, nzt
948	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
949	ENDDO
950	ENDIF
951	ENDIF
952
953	ENDIF
954
955	!
956	!-- Tendency terms for v-velocity component
957	IF ( .NOT. outflow_s .OR. j > nys ) THEN
958
959	tend(:,j,i) = 0.0
960	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
961	CALL advec_v_pw( i, j )
962	ELSE
963	CALL advec_v_up( i, j )
964	ENDIF
965	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
966	THEN
967	CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, &
968	u_m, v_m, vsws_m, vswst_m, w_m )
969	ELSE
970	CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
971	vsws, vswst, w )
972	ENDIF
973	CALL coriolis( i, j, 2 )
974
975	!
976	!-- Drag by plant canopy
977	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
978
979	!
980	!-- External pressure gradient
981	IF ( dp_external ) THEN
982	DO k = dp_level_ind_b+1, nzt
983	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
984	ENDDO
985	ENDIF
986
987	CALL user_actions( i, j, 'v-tendency' )
988
989	!
990	!-- Prognostic equation for v-velocity component
991	DO k = nzb_v_inner(j,i)+1, nzt
992	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
993	dt_3d * ( &
994	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
995	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
996	) - &
997	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
998	ENDDO
999
1000	!
1001	!-- Calculate tendencies for the next Runge-Kutta step
1002	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1003	IF ( intermediate_timestep_count == 1 ) THEN
1004	DO k = nzb_v_inner(j,i)+1, nzt
1005	tv_m(k,j,i) = tend(k,j,i)
1006	ENDDO
1007	ELSEIF ( intermediate_timestep_count < &
1008	intermediate_timestep_count_max ) THEN
1009	DO k = nzb_v_inner(j,i)+1, nzt
1010	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1011	ENDDO
1012	ENDIF
1013	ENDIF
1014
1015	ENDIF
1016
1017	!
1018	!-- Tendency terms for w-velocity component
1019	tend(:,j,i) = 0.0
1020	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1021	CALL advec_w_pw( i, j )
1022	ELSE
1023	CALL advec_w_up( i, j )
1024	ENDIF
1025	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
1026	THEN
1027	CALL diffusion_w( i, j, ddzu, ddzw, km_m, km_damp_x, &
1028	km_damp_y, tend, u_m, v_m, w_m )
1029	ELSE
1030	CALL diffusion_w( i, j, ddzu, ddzw, km, km_damp_x, km_damp_y, &
1031	tend, u, v, w )
1032	ENDIF
1033	CALL coriolis( i, j, 3 )
1034	IF ( ocean ) THEN
1035	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
1036	ELSE
1037	IF ( .NOT. humidity ) THEN
1038	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
1039	ELSE
1040	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
1041	ENDIF
1042	ENDIF
1043
1044	!
1045	!-- Drag by plant canopy
1046	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
1047
1048	CALL user_actions( i, j, 'w-tendency' )
1049
1050	!
1051	!-- Prognostic equation for w-velocity component
1052	DO k = nzb_w_inner(j,i)+1, nzt-1
1053	w_p(k,j,i) = ( 1.0-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1054	dt_3d * ( &
1055	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1056	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1057	) - &
1058	tsc(5) * rdf(k) * w(k,j,i)
1059	ENDDO
1060
1061	!
1062	!-- Calculate tendencies for the next Runge-Kutta step
1063	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1064	IF ( intermediate_timestep_count == 1 ) THEN
1065	DO k = nzb_w_inner(j,i)+1, nzt-1
1066	tw_m(k,j,i) = tend(k,j,i)
1067	ENDDO
1068	ELSEIF ( intermediate_timestep_count < &
1069	intermediate_timestep_count_max ) THEN
1070	DO k = nzb_w_inner(j,i)+1, nzt-1
1071	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1072	ENDDO
1073	ENDIF
1074	ENDIF
1075
1076	!
1077	!-- Tendency terms for potential temperature
1078	tend(:,j,i) = 0.0
1079	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1080	CALL advec_s_pw( i, j, pt )
1081	ELSE
1082	CALL advec_s_up( i, j, pt )
1083	ENDIF
1084	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) &
1085	THEN
1086	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, pt_m, shf_m, &
1087	tswst_m, wall_heatflux, tend )
1088	ELSE
1089	CALL diffusion_s( i, j, ddzu, ddzw, kh, pt, shf, tswst, &
1090	wall_heatflux, tend )
1091	ENDIF
1092
1093	!
1094	!-- If required compute heating/cooling due to long wave radiation
1095	!-- processes
1096	IF ( radiation ) THEN
1097	CALL calc_radiation( i, j )
1098	ENDIF
1099
1100	!
1101	!-- If required compute impact of latent heat due to precipitation
1102	IF ( precipitation ) THEN
1103	CALL impact_of_latent_heat( i, j )
1104	ENDIF
1105
1106	!
1107	!-- Consideration of heat sources within the plant canopy
1108	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1109	CALL plant_canopy_model( i, j, 4 )
1110	ENDIF
1111
1112
1113	!-- If required compute influence of large-scale subsidence/ascent
1114	IF ( large_scale_subsidence ) THEN
1115	CALL subsidence ( i, j, tend, pt, pt_init )
1116	ENDIF
1117
1118	CALL user_actions( i, j, 'pt-tendency' )
1119
1120	!
1121	!-- Prognostic equation for potential temperature
1122	DO k = nzb_s_inner(j,i)+1, nzt
1123	pt_p(k,j,i) = ( 1.0-tsc(1) ) * pt_m(k,j,i) + tsc(1)*pt(k,j,i) +&
1124	dt_3d * ( &
1125	tsc(2) * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1126	) - &
1127	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1128	ENDDO
1129
1130	!
1131	!-- Calculate tendencies for the next Runge-Kutta step
1132	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1133	IF ( intermediate_timestep_count == 1 ) THEN
1134	DO k = nzb_s_inner(j,i)+1, nzt
1135	tpt_m(k,j,i) = tend(k,j,i)
1136	ENDDO
1137	ELSEIF ( intermediate_timestep_count < &
1138	intermediate_timestep_count_max ) THEN
1139	DO k = nzb_s_inner(j,i)+1, nzt
1140	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1141	5.3125 * tpt_m(k,j,i)
1142	ENDDO
1143	ENDIF
1144	ENDIF
1145
1146	!
1147	!-- If required, compute prognostic equation for salinity
1148	IF ( ocean ) THEN
1149
1150	!
1151	!-- Tendency-terms for salinity
1152	tend(:,j,i) = 0.0
1153	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1154	THEN
1155	CALL advec_s_pw( i, j, sa )
1156	ELSE
1157	CALL advec_s_up( i, j, sa )
1158	ENDIF
1159	CALL diffusion_s( i, j, ddzu, ddzw, kh, sa, saswsb, saswst, &
1160	wall_salinityflux, tend )
1161
1162	CALL user_actions( i, j, 'sa-tendency' )
1163
1164	!
1165	!-- Prognostic equation for salinity
1166	DO k = nzb_s_inner(j,i)+1, nzt
1167	sa_p(k,j,i) = tsc(1) * sa(k,j,i) + &
1168	dt_3d * ( &
1169	tsc(2) * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1170	) - &
1171	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
1172	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1173	ENDDO
1174
1175	!
1176	!-- Calculate tendencies for the next Runge-Kutta step
1177	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1178	IF ( intermediate_timestep_count == 1 ) THEN
1179	DO k = nzb_s_inner(j,i)+1, nzt
1180	tsa_m(k,j,i) = tend(k,j,i)
1181	ENDDO
1182	ELSEIF ( intermediate_timestep_count < &
1183	intermediate_timestep_count_max ) THEN
1184	DO k = nzb_s_inner(j,i)+1, nzt
1185	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1186	5.3125 * tsa_m(k,j,i)
1187	ENDDO
1188	ENDIF
1189	ENDIF
1190
1191	!
1192	!-- Calculate density by the equation of state for seawater
1193	CALL eqn_state_seawater( i, j )
1194
1195	ENDIF
1196
1197	!
1198	!-- If required, compute prognostic equation for total water content /
1199	!-- scalar
1200	IF ( humidity .OR. passive_scalar ) THEN
1201
1202	!
1203	!-- Tendency-terms for total water content / scalar
1204	tend(:,j,i) = 0.0
1205	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1206	THEN
1207	CALL advec_s_pw( i, j, q )
1208	ELSE
1209	CALL advec_s_up( i, j, q )
1210	ENDIF
1211	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
1212	THEN
1213	CALL diffusion_s( i, j, ddzu, ddzw, kh_m, q_m, qsws_m, &
1214	qswst_m, wall_qflux, tend )
1215	ELSE
1216	CALL diffusion_s( i, j, ddzu, ddzw, kh, q, qsws, qswst, &
1217	wall_qflux, tend )
1218	ENDIF
1219
1220	!
1221	!-- If required compute decrease of total water content due to
1222	!-- precipitation
1223	IF ( precipitation ) THEN
1224	CALL calc_precipitation( i, j )
1225	ENDIF
1226
1227	!
1228	!-- Sink or source of scalar concentration due to canopy elements
1229	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
1230
1231
1232	!-- If required compute influence of large-scale subsidence/ascent
1233	IF ( large_scale_subsidence ) THEN
1234	CALL subsidence ( i, j, tend, q, q_init )
1235	ENDIF
1236
1237
1238	CALL user_actions( i, j, 'q-tendency' )
1239
1240	!
1241	!-- Prognostic equation for total water content / scalar
1242	DO k = nzb_s_inner(j,i)+1, nzt
1243	q_p(k,j,i) = ( 1.0-tsc(1) ) * q_m(k,j,i) + tsc(1)*q(k,j,i) +&
1244	dt_3d * ( &
1245	tsc(2) * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1246	) - &
1247	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1248	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1249	ENDDO
1250
1251	!
1252	!-- Calculate tendencies for the next Runge-Kutta step
1253	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1254	IF ( intermediate_timestep_count == 1 ) THEN
1255	DO k = nzb_s_inner(j,i)+1, nzt
1256	tq_m(k,j,i) = tend(k,j,i)
1257	ENDDO
1258	ELSEIF ( intermediate_timestep_count < &
1259	intermediate_timestep_count_max ) THEN
1260	DO k = nzb_s_inner(j,i)+1, nzt
1261	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1262	5.3125 * tq_m(k,j,i)
1263	ENDDO
1264	ENDIF
1265	ENDIF
1266
1267	ENDIF
1268
1269	!
1270	!-- If required, compute prognostic equation for turbulent kinetic
1271	!-- energy (TKE)
1272	IF ( .NOT. constant_diffusion ) THEN
1273
1274	!
1275	!-- Tendency-terms for TKE
1276	tend(:,j,i) = 0.0
1277	IF ( ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) &
1278	.AND. .NOT. use_upstream_for_tke ) THEN
1279	CALL advec_s_pw( i, j, e )
1280	ELSE
1281	CALL advec_s_up( i, j, e )
1282	ENDIF
1283	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' )&
1284	THEN
1285	IF ( .NOT. humidity ) THEN
1286	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1287	km_m, l_grid, pt_m, pt_reference, &
1288	rif_m, tend, zu, zw )
1289	ELSE
1290	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
1291	km_m, l_grid, vpt_m, pt_reference, &
1292	rif_m, tend, zu, zw )
1293	ENDIF
1294	ELSE
1295	IF ( .NOT. humidity ) THEN
1296	IF ( ocean ) THEN
1297	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
1298	km, l_grid, prho, prho_reference, &
1299	rif, tend, zu, zw )
1300	ELSE
1301	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
1302	km, l_grid, pt, pt_reference, rif, &
1303	tend, zu, zw )
1304	ENDIF
1305	ELSE
1306	CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
1307	l_grid, vpt, pt_reference, rif, tend, &
1308	zu, zw )
1309	ENDIF
1310	ENDIF
1311	CALL production_e( i, j )
1312
1313	!
1314	!-- Additional sink term for flows through plant canopies
1315	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
1316
1317	CALL user_actions( i, j, 'e-tendency' )
1318
1319	!
1320	!-- Prognostic equation for TKE.
1321	!-- Eliminate negative TKE values, which can occur due to numerical
1322	!-- reasons in the course of the integration. In such cases the old
1323	!-- TKE value is reduced by 90%.
1324	DO k = nzb_s_inner(j,i)+1, nzt
1325	e_p(k,j,i) = ( 1.0-tsc(1) ) * e_m(k,j,i) + tsc(1)*e(k,j,i) +&
1326	dt_3d * ( &
1327	tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
1328	)
1329	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1330	ENDDO
1331
1332	!
1333	!-- Calculate tendencies for the next Runge-Kutta step
1334	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1335	IF ( intermediate_timestep_count == 1 ) THEN
1336	DO k = nzb_s_inner(j,i)+1, nzt
1337	te_m(k,j,i) = tend(k,j,i)
1338	ENDDO
1339	ELSEIF ( intermediate_timestep_count < &
1340	intermediate_timestep_count_max ) THEN
1341	DO k = nzb_s_inner(j,i)+1, nzt
1342	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1343	5.3125 * te_m(k,j,i)
1344	ENDDO
1345	ENDIF
1346	ENDIF
1347
1348	ENDIF ! TKE equation
1349
1350	ENDDO
1351	ENDDO
1352	!$OMP END PARALLEL
1353
1354	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1355
1356
1357	END SUBROUTINE prognostic_equations_cache
1358
1359
1360	SUBROUTINE prognostic_equations_vector
1361
1362	!------------------------------------------------------------------------------!
1363	! Version for vector machines
1364	!------------------------------------------------------------------------------!
1365
1366	IMPLICIT NONE
1367
1368	CHARACTER (LEN=9) :: time_to_string
1369	INTEGER :: i, j, k
1370	REAL :: sat, sbt
1371
1372	!
1373	!-- Calculate those variables needed in the tendency terms which need
1374	!-- global communication
1375	CALL calc_mean_profile( pt, 4 )
1376	IF ( ocean ) CALL calc_mean_profile( rho, 64 )
1377	IF ( humidity ) CALL calc_mean_profile( vpt, 44 )
1378
1379	!
1380	!-- u-velocity component
1381	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1382
1383	!
1384	!-- u-tendency terms with communication
1385	IF ( momentum_advec == 'ups-scheme' ) THEN
1386	tend = 0.0
1387	CALL advec_u_ups
1388	ENDIF
1389
1390	!
1391	!-- u-tendency terms with no communication
1392	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1393	tend = 0.0
1394	CALL advec_u_pw
1395	ELSE
1396	IF ( momentum_advec /= 'ups-scheme' ) THEN
1397	tend = 0.0
1398	CALL advec_u_up
1399	ENDIF
1400	ENDIF
1401	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1402	CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
1403	uswst_m, v_m, w_m )
1404	ELSE
1405	CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
1406	ENDIF
1407	CALL coriolis( 1 )
1408	IF ( sloping_surface ) CALL buoyancy( pt, pt_reference, 1, 4 )
1409
1410	!
1411	!-- Drag by plant canopy
1412	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
1413
1414	!
1415	!-- External pressure gradient
1416	IF ( dp_external ) THEN
1417	DO i = nxlu, nxr
1418	DO j = nys, nyn
1419	DO k = dp_level_ind_b+1, nzt
1420	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1421	ENDDO
1422	ENDDO
1423	ENDDO
1424	ENDIF
1425
1426	CALL user_actions( 'u-tendency' )
1427
1428	!
1429	!-- Prognostic equation for u-velocity component
1430	DO i = nxlu, nxr
1431	DO j = nys, nyn
1432	DO k = nzb_u_inner(j,i)+1, nzt
1433	u_p(k,j,i) = ( 1.0-tsc(1) ) * u_m(k,j,i) + tsc(1) * u(k,j,i) + &
1434	dt_3d * ( &
1435	tsc(2) * tend(k,j,i) + tsc(3) * tu_m(k,j,i) &
1436	- tsc(4) * ( p(k,j,i) - p(k,j,i-1) ) * ddx &
1437	) - &
1438	tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1439	ENDDO
1440	ENDDO
1441	ENDDO
1442
1443	!
1444	!-- Calculate tendencies for the next Runge-Kutta step
1445	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1446	IF ( intermediate_timestep_count == 1 ) THEN
1447	DO i = nxlu, nxr
1448	DO j = nys, nyn
1449	DO k = nzb_u_inner(j,i)+1, nzt
1450	tu_m(k,j,i) = tend(k,j,i)
1451	ENDDO
1452	ENDDO
1453	ENDDO
1454	ELSEIF ( intermediate_timestep_count < &
1455	intermediate_timestep_count_max ) THEN
1456	DO i = nxlu, nxr
1457	DO j = nys, nyn
1458	DO k = nzb_u_inner(j,i)+1, nzt
1459	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1460	ENDDO
1461	ENDDO
1462	ENDDO
1463	ENDIF
1464	ENDIF
1465
1466	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1467
1468	!
1469	!-- v-velocity component
1470	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1471
1472	!
1473	!-- v-tendency terms with communication
1474	IF ( momentum_advec == 'ups-scheme' ) THEN
1475	tend = 0.0
1476	CALL advec_v_ups
1477	ENDIF
1478
1479	!
1480	!-- v-tendency terms with no communication
1481	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1482	tend = 0.0
1483	CALL advec_v_pw
1484	ELSE
1485	IF ( momentum_advec /= 'ups-scheme' ) THEN
1486	tend = 0.0
1487	CALL advec_v_up
1488	ENDIF
1489	ENDIF
1490	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1491	CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
1492	vswst_m, w_m )
1493	ELSE
1494	CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
1495	ENDIF
1496	CALL coriolis( 2 )
1497
1498	!
1499	!-- Drag by plant canopy
1500	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
1501
1502	!
1503	!-- External pressure gradient
1504	IF ( dp_external ) THEN
1505	DO i = nxl, nxr
1506	DO j = nysv, nyn
1507	DO k = dp_level_ind_b+1, nzt
1508	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1509	ENDDO
1510	ENDDO
1511	ENDDO
1512	ENDIF
1513
1514	CALL user_actions( 'v-tendency' )
1515
1516	!
1517	!-- Prognostic equation for v-velocity component
1518	DO i = nxl, nxr
1519	DO j = nysv, nyn
1520	DO k = nzb_v_inner(j,i)+1, nzt
1521	v_p(k,j,i) = ( 1.0-tsc(1) ) * v_m(k,j,i) + tsc(1) * v(k,j,i) + &
1522	dt_3d * ( &
1523	tsc(2) * tend(k,j,i) + tsc(3) * tv_m(k,j,i) &
1524	- tsc(4) * ( p(k,j,i) - p(k,j-1,i) ) * ddy &
1525	) - &
1526	tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1527	ENDDO
1528	ENDDO
1529	ENDDO
1530
1531	!
1532	!-- Calculate tendencies for the next Runge-Kutta step
1533	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1534	IF ( intermediate_timestep_count == 1 ) THEN
1535	DO i = nxl, nxr
1536	DO j = nysv, nyn
1537	DO k = nzb_v_inner(j,i)+1, nzt
1538	tv_m(k,j,i) = tend(k,j,i)
1539	ENDDO
1540	ENDDO
1541	ENDDO
1542	ELSEIF ( intermediate_timestep_count < &
1543	intermediate_timestep_count_max ) THEN
1544	DO i = nxl, nxr
1545	DO j = nysv, nyn
1546	DO k = nzb_v_inner(j,i)+1, nzt
1547	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1548	ENDDO
1549	ENDDO
1550	ENDDO
1551	ENDIF
1552	ENDIF
1553
1554	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1555
1556	!
1557	!-- w-velocity component
1558	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1559
1560	!
1561	!-- w-tendency terms with communication
1562	IF ( momentum_advec == 'ups-scheme' ) THEN
1563	tend = 0.0
1564	CALL advec_w_ups
1565	ENDIF
1566
1567	!
1568	!-- w-tendency terms with no communication
1569	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1570	tend = 0.0
1571	CALL advec_w_pw
1572	ELSE
1573	IF ( momentum_advec /= 'ups-scheme' ) THEN
1574	tend = 0.0
1575	CALL advec_w_up
1576	ENDIF
1577	ENDIF
1578	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1579	CALL diffusion_w( ddzu, ddzw, km_m, km_damp_x, km_damp_y, tend, u_m, &
1580	v_m, w_m )
1581	ELSE
1582	CALL diffusion_w( ddzu, ddzw, km, km_damp_x, km_damp_y, tend, u, v, w )
1583	ENDIF
1584	CALL coriolis( 3 )
1585	IF ( ocean ) THEN
1586	CALL buoyancy( rho, rho_reference, 3, 64 )
1587	ELSE
1588	IF ( .NOT. humidity ) THEN
1589	CALL buoyancy( pt, pt_reference, 3, 4 )
1590	ELSE
1591	CALL buoyancy( vpt, pt_reference, 3, 44 )
1592	ENDIF
1593	ENDIF
1594
1595	!
1596	!-- Drag by plant canopy
1597	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
1598
1599	CALL user_actions( 'w-tendency' )
1600
1601	!
1602	!-- Prognostic equation for w-velocity component
1603	DO i = nxl, nxr
1604	DO j = nys, nyn
1605	DO k = nzb_w_inner(j,i)+1, nzt-1
1606	w_p(k,j,i) = ( 1-tsc(1) ) * w_m(k,j,i) + tsc(1) * w(k,j,i) + &
1607	dt_3d * ( &
1608	tsc(2) * tend(k,j,i) + tsc(3) * tw_m(k,j,i) &
1609	- tsc(4) * ( p(k+1,j,i) - p(k,j,i) ) * ddzu(k+1) &
1610	) - &
1611	tsc(5) * rdf(k) * w(k,j,i)
1612	ENDDO
1613	ENDDO
1614	ENDDO
1615
1616	!
1617	!-- Calculate tendencies for the next Runge-Kutta step
1618	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1619	IF ( intermediate_timestep_count == 1 ) THEN
1620	DO i = nxl, nxr
1621	DO j = nys, nyn
1622	DO k = nzb_w_inner(j,i)+1, nzt-1
1623	tw_m(k,j,i) = tend(k,j,i)
1624	ENDDO
1625	ENDDO
1626	ENDDO
1627	ELSEIF ( intermediate_timestep_count < &
1628	intermediate_timestep_count_max ) THEN
1629	DO i = nxl, nxr
1630	DO j = nys, nyn
1631	DO k = nzb_w_inner(j,i)+1, nzt-1
1632	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1633	ENDDO
1634	ENDDO
1635	ENDDO
1636	ENDIF
1637	ENDIF
1638
1639	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1640
1641	!
1642	!-- potential temperature
1643	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1644
1645	!
1646	!-- pt-tendency terms with communication
1647	sat = tsc(1)
1648	sbt = tsc(2)
1649	IF ( scalar_advec == 'bc-scheme' ) THEN
1650
1651	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1652	!
1653	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1654	!-- switched on. Thus:
1655	sat = 1.0
1656	sbt = 1.0
1657	ENDIF
1658	tend = 0.0
1659	CALL advec_s_bc( pt, 'pt' )
1660	ELSE
1661	IF ( tsc(2) /= 2.0 .AND. scalar_advec == 'ups-scheme' ) THEN
1662	tend = 0.0
1663	CALL advec_s_ups( pt, 'pt' )
1664	ENDIF
1665	ENDIF
1666
1667	!
1668	!-- pt-tendency terms with no communication
1669	IF ( scalar_advec == 'bc-scheme' ) THEN
1670	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1671	tend )
1672	ELSE
1673	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1674	tend = 0.0
1675	CALL advec_s_pw( pt )
1676	ELSE
1677	IF ( scalar_advec /= 'ups-scheme' ) THEN
1678	tend = 0.0
1679	CALL advec_s_up( pt )
1680	ENDIF
1681	ENDIF
1682	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1683	CALL diffusion_s( ddzu, ddzw, kh_m, pt_m, shf_m, tswst_m, &
1684	wall_heatflux, tend )
1685	ELSE
1686	CALL diffusion_s( ddzu, ddzw, kh, pt, shf, tswst, wall_heatflux, &
1687	tend )
1688	ENDIF
1689	ENDIF
1690
1691	!
1692	!-- If required compute heating/cooling due to long wave radiation
1693	!-- processes
1694	IF ( radiation ) THEN
1695	CALL calc_radiation
1696	ENDIF
1697
1698	!
1699	!-- If required compute impact of latent heat due to precipitation
1700	IF ( precipitation ) THEN
1701	CALL impact_of_latent_heat
1702	ENDIF
1703
1704	!
1705	!-- Consideration of heat sources within the plant canopy
1706	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1707	CALL plant_canopy_model( 4 )
1708	ENDIF
1709
1710	!--If required compute influence of large-scale subsidence/ascent
1711	IF ( large_scale_subsidence ) THEN
1712	CALL subsidence ( tend, pt, pt_init )
1713	ENDIF
1714
1715	CALL user_actions( 'pt-tendency' )
1716
1717	!
1718	!-- Prognostic equation for potential temperature
1719	DO i = nxl, nxr
1720	DO j = nys, nyn
1721	DO k = nzb_s_inner(j,i)+1, nzt
1722	pt_p(k,j,i) = ( 1 - sat ) * pt_m(k,j,i) + sat * pt(k,j,i) + &
1723	dt_3d * ( &
1724	sbt * tend(k,j,i) + tsc(3) * tpt_m(k,j,i) &
1725	) - &
1726	tsc(5) * rdf(k) * ( pt(k,j,i) - pt_init(k) )
1727	ENDDO
1728	ENDDO
1729	ENDDO
1730
1731	!
1732	!-- Calculate tendencies for the next Runge-Kutta step
1733	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1734	IF ( intermediate_timestep_count == 1 ) THEN
1735	DO i = nxl, nxr
1736	DO j = nys, nyn
1737	DO k = nzb_s_inner(j,i)+1, nzt
1738	tpt_m(k,j,i) = tend(k,j,i)
1739	ENDDO
1740	ENDDO
1741	ENDDO
1742	ELSEIF ( intermediate_timestep_count < &
1743	intermediate_timestep_count_max ) THEN
1744	DO i = nxl, nxr
1745	DO j = nys, nyn
1746	DO k = nzb_s_inner(j,i)+1, nzt
1747	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1748	ENDDO
1749	ENDDO
1750	ENDDO
1751	ENDIF
1752	ENDIF
1753
1754	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1755
1756	!
1757	!-- If required, compute prognostic equation for salinity
1758	IF ( ocean ) THEN
1759
1760	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1761
1762	!
1763	!-- sa-tendency terms with communication
1764	sat = tsc(1)
1765	sbt = tsc(2)
1766	IF ( scalar_advec == 'bc-scheme' ) THEN
1767
1768	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1769	!
1770	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1771	!-- switched on. Thus:
1772	sat = 1.0
1773	sbt = 1.0
1774	ENDIF
1775	tend = 0.0
1776	CALL advec_s_bc( sa, 'sa' )
1777	ELSE
1778	IF ( tsc(2) /= 2.0 ) THEN
1779	IF ( scalar_advec == 'ups-scheme' ) THEN
1780	tend = 0.0
1781	CALL advec_s_ups( sa, 'sa' )
1782	ENDIF
1783	ENDIF
1784	ENDIF
1785
1786	!
1787	!-- sa-tendency terms with no communication
1788	IF ( scalar_advec == 'bc-scheme' ) THEN
1789	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1790	wall_salinityflux, tend )
1791	ELSE
1792	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1793	tend = 0.0
1794	CALL advec_s_pw( sa )
1795	ELSE
1796	IF ( scalar_advec /= 'ups-scheme' ) THEN
1797	tend = 0.0
1798	CALL advec_s_up( sa )
1799	ENDIF
1800	ENDIF
1801	CALL diffusion_s( ddzu, ddzw, kh, sa, saswsb, saswst, &
1802	wall_salinityflux, tend )
1803	ENDIF
1804
1805	CALL user_actions( 'sa-tendency' )
1806
1807	!
1808	!-- Prognostic equation for salinity
1809	DO i = nxl, nxr
1810	DO j = nys, nyn
1811	DO k = nzb_s_inner(j,i)+1, nzt
1812	sa_p(k,j,i) = sat * sa(k,j,i) + &
1813	dt_3d * ( &
1814	sbt * tend(k,j,i) + tsc(3) * tsa_m(k,j,i) &
1815	) - &
1816	tsc(5) * rdf(k) * ( sa(k,j,i) - sa_init(k) )
1817	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1818	ENDDO
1819	ENDDO
1820	ENDDO
1821
1822	!
1823	!-- Calculate tendencies for the next Runge-Kutta step
1824	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1825	IF ( intermediate_timestep_count == 1 ) THEN
1826	DO i = nxl, nxr
1827	DO j = nys, nyn
1828	DO k = nzb_s_inner(j,i)+1, nzt
1829	tsa_m(k,j,i) = tend(k,j,i)
1830	ENDDO
1831	ENDDO
1832	ENDDO
1833	ELSEIF ( intermediate_timestep_count < &
1834	intermediate_timestep_count_max ) THEN
1835	DO i = nxl, nxr
1836	DO j = nys, nyn
1837	DO k = nzb_s_inner(j,i)+1, nzt
1838	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1839	5.3125 * tsa_m(k,j,i)
1840	ENDDO
1841	ENDDO
1842	ENDDO
1843	ENDIF
1844	ENDIF
1845
1846	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1847
1848	!
1849	!-- Calculate density by the equation of state for seawater
1850	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1851	CALL eqn_state_seawater
1852	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1853
1854	ENDIF
1855
1856	!
1857	!-- If required, compute prognostic equation for total water content / scalar
1858	IF ( humidity .OR. passive_scalar ) THEN
1859
1860	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1861
1862	!
1863	!-- Scalar/q-tendency terms with communication
1864	sat = tsc(1)
1865	sbt = tsc(2)
1866	IF ( scalar_advec == 'bc-scheme' ) THEN
1867
1868	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1869	!
1870	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1871	!-- switched on. Thus:
1872	sat = 1.0
1873	sbt = 1.0
1874	ENDIF
1875	tend = 0.0
1876	CALL advec_s_bc( q, 'q' )
1877	ELSE
1878	IF ( tsc(2) /= 2.0 ) THEN
1879	IF ( scalar_advec == 'ups-scheme' ) THEN
1880	tend = 0.0
1881	CALL advec_s_ups( q, 'q' )
1882	ENDIF
1883	ENDIF
1884	ENDIF
1885
1886	!
1887	!-- Scalar/q-tendency terms with no communication
1888	IF ( scalar_advec == 'bc-scheme' ) THEN
1889	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, wall_qflux, tend )
1890	ELSE
1891	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
1892	tend = 0.0
1893	CALL advec_s_pw( q )
1894	ELSE
1895	IF ( scalar_advec /= 'ups-scheme' ) THEN
1896	tend = 0.0
1897	CALL advec_s_up( q )
1898	ENDIF
1899	ENDIF
1900	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
1901	CALL diffusion_s( ddzu, ddzw, kh_m, q_m, qsws_m, qswst_m, &
1902	wall_qflux, tend )
1903	ELSE
1904	CALL diffusion_s( ddzu, ddzw, kh, q, qsws, qswst, &
1905	wall_qflux, tend )
1906	ENDIF
1907	ENDIF
1908
1909	!
1910	!-- If required compute decrease of total water content due to
1911	!-- precipitation
1912	IF ( precipitation ) THEN
1913	CALL calc_precipitation
1914	ENDIF
1915
1916	!
1917	!-- Sink or source of scalar concentration due to canopy elements
1918	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
1919
1920	!
1921	!-- If required compute influence of large-scale subsidence/ascent
1922	IF ( large_scale_subsidence ) THEN
1923	CALL subsidence ( tend, q, q_init )
1924	ENDIF
1925
1926	CALL user_actions( 'q-tendency' )
1927
1928	!
1929	!-- Prognostic equation for total water content / scalar
1930	DO i = nxl, nxr
1931	DO j = nys, nyn
1932	DO k = nzb_s_inner(j,i)+1, nzt
1933	q_p(k,j,i) = ( 1 - sat ) * q_m(k,j,i) + sat * q(k,j,i) + &
1934	dt_3d * ( &
1935	sbt * tend(k,j,i) + tsc(3) * tq_m(k,j,i) &
1936	) - &
1937	tsc(5) * rdf(k) * ( q(k,j,i) - q_init(k) )
1938	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1939	ENDDO
1940	ENDDO
1941	ENDDO
1942
1943	!
1944	!-- Calculate tendencies for the next Runge-Kutta step
1945	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1946	IF ( intermediate_timestep_count == 1 ) THEN
1947	DO i = nxl, nxr
1948	DO j = nys, nyn
1949	DO k = nzb_s_inner(j,i)+1, nzt
1950	tq_m(k,j,i) = tend(k,j,i)
1951	ENDDO
1952	ENDDO
1953	ENDDO
1954	ELSEIF ( intermediate_timestep_count < &
1955	intermediate_timestep_count_max ) THEN
1956	DO i = nxl, nxr
1957	DO j = nys, nyn
1958	DO k = nzb_s_inner(j,i)+1, nzt
1959	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1960	ENDDO
1961	ENDDO
1962	ENDDO
1963	ENDIF
1964	ENDIF
1965
1966	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1967
1968	ENDIF
1969
1970	!
1971	!-- If required, compute prognostic equation for turbulent kinetic
1972	!-- energy (TKE)
1973	IF ( .NOT. constant_diffusion ) THEN
1974
1975	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1976
1977	!
1978	!-- TKE-tendency terms with communication
1979	CALL production_e_init
1980
1981	sat = tsc(1)
1982	sbt = tsc(2)
1983	IF ( .NOT. use_upstream_for_tke ) THEN
1984	IF ( scalar_advec == 'bc-scheme' ) THEN
1985
1986	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1987	!
1988	!-- Bott-Chlond scheme always uses Euler time step when leapfrog is
1989	!-- switched on. Thus:
1990	sat = 1.0
1991	sbt = 1.0
1992	ENDIF
1993	tend = 0.0
1994	CALL advec_s_bc( e, 'e' )
1995	ELSE
1996	IF ( tsc(2) /= 2.0 ) THEN
1997	IF ( scalar_advec == 'ups-scheme' ) THEN
1998	tend = 0.0
1999	CALL advec_s_ups( e, 'e' )
2000	ENDIF
2001	ENDIF
2002	ENDIF
2003	ENDIF
2004
2005	!
2006	!-- TKE-tendency terms with no communication
2007	IF ( scalar_advec == 'bc-scheme' .AND. .NOT. use_upstream_for_tke ) &
2008	THEN
2009	IF ( .NOT. humidity ) THEN
2010	IF ( ocean ) THEN
2011	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2012	prho, prho_reference, rif, tend, zu, zw )
2013	ELSE
2014	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
2015	pt_reference, rif, tend, zu, zw )
2016	ENDIF
2017	ELSE
2018	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
2019	pt_reference, rif, tend, zu, zw )
2020	ENDIF
2021	ELSE
2022	IF ( use_upstream_for_tke ) THEN
2023	tend = 0.0
2024	CALL advec_s_up( e )
2025	ELSE
2026	IF ( tsc(2) == 2.0 .OR. timestep_scheme(1:5) == 'runge' ) THEN
2027	tend = 0.0
2028	CALL advec_s_pw( e )
2029	ELSE
2030	IF ( scalar_advec /= 'ups-scheme' ) THEN
2031	tend = 0.0
2032	CALL advec_s_up( e )
2033	ENDIF
2034	ENDIF
2035	ENDIF
2036	IF ( tsc(2) == 2.0 .AND. timestep_scheme(1:8) == 'leapfrog' ) THEN
2037	IF ( .NOT. humidity ) THEN
2038	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2039	pt_m, pt_reference, rif_m, tend, zu, zw )
2040	ELSE
2041	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
2042	vpt_m, pt_reference, rif_m, tend, zu, zw )
2043	ENDIF
2044	ELSE
2045	IF ( .NOT. humidity ) THEN
2046	IF ( ocean ) THEN
2047	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2048	prho, prho_reference, rif, tend, zu, zw )
2049	ELSE
2050	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
2051	pt, pt_reference, rif, tend, zu, zw )
2052	ENDIF
2053	ELSE
2054	CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
2055	pt_reference, rif, tend, zu, zw )
2056	ENDIF
2057	ENDIF
2058	ENDIF
2059	CALL production_e
2060
2061	!
2062	!-- Additional sink term for flows through plant canopies
2063	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
2064	CALL user_actions( 'e-tendency' )
2065
2066	!
2067	!-- Prognostic equation for TKE.
2068	!-- Eliminate negative TKE values, which can occur due to numerical
2069	!-- reasons in the course of the integration. In such cases the old TKE
2070	!-- value is reduced by 90%.
2071	DO i = nxl, nxr
2072	DO j = nys, nyn
2073	DO k = nzb_s_inner(j,i)+1, nzt
2074	e_p(k,j,i) = ( 1 - sat ) * e_m(k,j,i) + sat * e(k,j,i) + &
2075	dt_3d * ( &
2076	sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) &
2077	)
2078	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
2079	ENDDO
2080	ENDDO
2081	ENDDO
2082
2083	!
2084	!-- Calculate tendencies for the next Runge-Kutta step
2085	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2086	IF ( intermediate_timestep_count == 1 ) THEN
2087	DO i = nxl, nxr
2088	DO j = nys, nyn
2089	DO k = nzb_s_inner(j,i)+1, nzt
2090	te_m(k,j,i) = tend(k,j,i)
2091	ENDDO
2092	ENDDO
2093	ENDDO
2094	ELSEIF ( intermediate_timestep_count < &
2095	intermediate_timestep_count_max ) THEN
2096	DO i = nxl, nxr
2097	DO j = nys, nyn
2098	DO k = nzb_s_inner(j,i)+1, nzt
2099	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
2100	ENDDO
2101	ENDDO
2102	ENDDO
2103	ENDIF
2104	ENDIF
2105
2106	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2107
2108	ENDIF
2109
2110
2111	END SUBROUTINE prognostic_equations_vector
2112
2113
2114	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |