Home

Context Navigation

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

Last change on this file since 1526 was 1518, checked in by hoffmann, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 84.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |