Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |