Home

Context Navigation

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

Last change on this file since 1504 was 1497, checked in by maronga, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 83.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |