Home

Context Navigation

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

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