Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |