Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |