Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |