Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |