Home

Context Navigation

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

Last change on this file since 388 was 388, checked in by raasch, 15 years ago
in-situ AND potential density are calculated and used in the ocean version
Property svn:keywords set to `Id`
File size: 70.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |