Home

Context Navigation

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

Last change on this file since 1777 was 1758, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 85.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |