Home

Context Navigation

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

Last change on this file since 669 was 669, checked in by gryschka, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 75.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |