Home

Context Navigation

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

Last change on this file since 4109 was 4109, checked in by suehring, 5 years ago
Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.
Property svn:keywords set to `Id`
File size: 67.7 KB

Line
1	!> @file prognostic_equations.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	! pass integer flag array to WS scalar advection routine which is now necessary
23	! as the flags may differ for scalars, e.g. pt can be cyclic while chemical
24	! species may be non-cyclic. Further, pass boundary flags.
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: prognostic_equations.f90 4109 2019-07-22 17:00:34Z suehring $
29	! Application of monotonic flux limiter for the vertical scalar advection
30	! up to the topography top (only for the cache-optimized version at the
31	! moment). Please note, at the moment the limiter is only applied for passive
32	! scalars.
33	!
34	! 4048 2019-06-21 21:00:21Z knoop
35	! Moved tcm_prognostic_equations to module_interface
36	!
37	! 3987 2019-05-22 09:52:13Z kanani
38	! Introduce alternative switch for debug output during timestepping
39	!
40	! 3956 2019-05-07 12:32:52Z monakurppa
41	! Removed salsa calls.
42	!
43	! 3931 2019-04-24 16:34:28Z schwenkel
44	! Correct/complete module_interface introduction for chemistry model
45	!
46	! 3899 2019-04-16 14:05:27Z monakurppa
47	! Corrections in the OpenMP version of salsa
48	!
49	! 3887 2019 -04-12 08:47:41Z schwenkel
50	! Implicit Bugfix for chemistry model, loop for non_transport_physics over
51	! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
52	!
53	! 3885 2019-04-11 11:29:34Z kanani
54	! Changes related to global restructuring of location messages and introduction
55	! of additional debug messages
56	!
57	! 3881 2019-04-10 09:31:22Z suehring
58	! Bugfix in OpenMP directive
59	!
60	! 3880 2019-04-08 21:43:02Z knoop
61	! Moved wtm_tendencies to module_interface_actions
62	!
63	! 3874 2019-04-08 16:53:48Z knoop
64	! Added non_transport_physics module interfaces and moved bcm code into it
65	!
66	! 3872 2019-04-08 15:03:06Z knoop
67	! Moving prognostic equations of bcm into bulk_cloud_model_mod
68	!
69	! 3864 2019-04-05 09:01:56Z monakurppa
70	! Modifications made for salsa:
71	! - salsa_prognostic_equations moved to salsa_mod (and the call to
72	! module_interface_mod)
73	! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
74	! ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
75	!
76	! 3840 2019-03-29 10:35:52Z knoop
77	! added USE chem_gasphase_mod for nspec, nspec and spc_names
78	!
79	! 3820 2019-03-27 11:53:41Z forkel
80	! renamed do_depo to deposition_dry (ecc)
81	!
82	! 3797 2019-03-15 11:15:38Z forkel
83	! Call chem_integegrate in OpenMP loop (ketelsen)
84	!
85	!
86	! 3771 2019-02-28 12:19:33Z raasch
87	! preprocessor directivs fro rrtmg added
88	!
89	! 3761 2019-02-25 15:31:42Z raasch
90	! unused variable removed
91	!
92	! 3719 2019-02-06 13:10:18Z kanani
93	! Cleaned up chemistry cpu measurements
94	!
95	! 3684 2019-01-20 20:20:58Z knoop
96	! OpenACC port for SPEC
97	!
98	! 3589 2018-11-30 15:09:51Z suehring
99	! Move the control parameter "salsa" from salsa_mod to control_parameters
100	! (M. Kurppa)
101	!
102	! 3582 2018-11-29 19:16:36Z suehring
103	! Implementation of a new aerosol module salsa.
104	! Remove cpu-logs from i,j loop in cache version.
105	!
106	! 3458 2018-10-30 14:51:23Z kanani
107	! remove duplicate USE chem_modules
108	! from chemistry branch r3443, banzhafs, basit:
109	! chem_depo call introduced
110	! code added for decycling chemistry
111	!
112	! 3386 2018-10-19 16:28:22Z gronemeier
113	! Renamed tcm_prognostic to tcm_prognostic_equations
114	!
115	! 3355 2018-10-16 14:03:34Z knoop
116	! (from branch resler)
117	! Fix for chemistry call
118	!
119	! 3302 2018-10-03 02:39:40Z raasch
120	! Stokes drift + wave breaking term added
121	!
122	! 3298 2018-10-02 12:21:11Z kanani
123	! Code added for decycling chemistry (basit)
124	!
125	! 3294 2018-10-01 02:37:10Z raasch
126	! changes concerning modularization of ocean option
127	!
128	! 3274 2018-09-24 15:42:55Z knoop
129	! Modularization of all bulk cloud physics code components
130	!
131	! 3241 2018-09-12 15:02:00Z raasch
132	! omp_get_thread_num now declared in omp directive
133	!
134	! 3183 2018-07-27 14:25:55Z suehring
135	! Remove unused variables from USE statements
136	!
137	! 3182 2018-07-27 13:36:03Z suehring
138	! Revise recent bugfix for nesting
139	!
140	! 3021 2018-05-16 08:14:20Z maronga
141	! Bugfix in IF clause for nesting
142	!
143	! 3014 2018-05-09 08:42:38Z maronga
144	! Fixed a bug in the IF condition to call pcm_tendency in case of
145	! potential temperature
146	!
147	! 2815 2018-02-19 11:29:57Z kanani
148	! Rename chem_tendency to chem_prognostic_equations,
149	! implement vector version for air chemistry
150	!
151	! 2766 2018-01-22 17:17:47Z kanani
152	! Removed preprocessor directive __chem
153	!
154	! 2746 2018-01-15 12:06:04Z suehring
155	! Move flag plant canopy to modules
156	!
157	! 2719 2018-01-02 09:02:06Z maronga
158	! Bugfix for last change.
159	!
160	! 2718 2018-01-02 08:49:38Z maronga
161	! Corrected "Former revisions" section
162	!
163	! 2696 2017-12-14 17:12:51Z kanani
164	! - Change in file header (GPL part)
165	! - Moved TKE equation to tcm_prognostic (TG)
166	! - Added switch for chemical reactions (RF, FK)
167	! - Implementation of chemistry module (RF, BK, FK)
168	!
169	! 2563 2017-10-19 15:36:10Z Giersch
170	! Variable wind_turbine moved to module control_parameters
171	!
172	! 2320 2017-07-21 12:47:43Z suehring
173	! Modularize large-scale forcing and nudging
174	!
175	! 2292 2017-06-20 09:51:42Z schwenkel
176	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
177	! includes two more prognostic equations for cloud drop concentration (nc)
178	! and cloud water content (qc).
179	!
180	! 2261 2017-06-08 14:25:57Z raasch
181	! bugfix for r2232: openmp directives removed
182	!
183	! 2233 2017-05-30 18:08:54Z suehring
184	!
185	! 2232 2017-05-30 17:47:52Z suehring
186	! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
187	! which is realized directly in diffusion_s now.
188	!
189	! 2192 2017-03-22 04:14:10Z raasch
190	! Bugfix for misplaced and missing openMP directives from r2155
191	!
192	! 2155 2017-02-21 09:57:40Z hoffmann
193	! Bugfix in the calculation of microphysical quantities on ghost points.
194	!
195	! 2118 2017-01-17 16:38:49Z raasch
196	! OpenACC version of subroutine removed
197	!
198	! 2031 2016-10-21 15:11:58Z knoop
199	! renamed variable rho to rho_ocean
200	!
201	! 2011 2016-09-19 17:29:57Z kanani
202	! Flag urban_surface is now defined in module control_parameters.
203	!
204	! 2007 2016-08-24 15:47:17Z kanani
205	! Added pt tendency calculation based on energy balance at urban surfaces
206	! (new urban surface model)
207	!
208	! 2000 2016-08-20 18:09:15Z knoop
209	! Forced header and separation lines into 80 columns
210	!
211	! 1976 2016-07-27 13:28:04Z maronga
212	! Simplied calls to radiation model
213	!
214	! 1960 2016-07-12 16:34:24Z suehring
215	! Separate humidity and passive scalar
216	!
217	! 1914 2016-05-26 14:44:07Z witha
218	! Added calls for wind turbine model
219	!
220	! 1873 2016-04-18 14:50:06Z maronga
221	! Module renamed (removed _mod)
222	!
223	! 1850 2016-04-08 13:29:27Z maronga
224	! Module renamed
225	!
226	! 1826 2016-04-07 12:01:39Z maronga
227	! Renamed canopy model calls.
228	!
229	! 1822 2016-04-07 07:49:42Z hoffmann
230	! Kessler microphysics scheme moved to microphysics.
231	!
232	! 1757 2016-02-22 15:49:32Z maronga
233	!
234	! 1691 2015-10-26 16:17:44Z maronga
235	! Added optional model spin-up without radiation / land surface model calls.
236	! Formatting corrections.
237	!
238	! 1682 2015-10-07 23:56:08Z knoop
239	! Code annotations made doxygen readable
240	!
241	! 1585 2015-04-30 07:05:52Z maronga
242	! Added call for temperature tendency calculation due to radiative flux divergence
243	!
244	! 1517 2015-01-07 19:12:25Z hoffmann
245	! advec_s_bc_mod addded, since advec_s_bc is now a module
246	!
247	! 1484 2014-10-21 10:53:05Z kanani
248	! Changes due to new module structure of the plant canopy model:
249	! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
250	! Removed double-listing of use_upstream_for_tke in ONLY-list of module
251	! control_parameters
252	!
253	! 1409 2014-05-23 12:11:32Z suehring
254	! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
255	! This ensures that left-hand side fluxes are also calculated for nxl in that
256	! case, even though the solution at nxl is overwritten in boundary_conds()
257	!
258	! 1398 2014-05-07 11:15:00Z heinze
259	! Rayleigh-damping for horizontal velocity components changed: instead of damping
260	! against ug and vg, damping against u_init and v_init is used to allow for a
261	! homogenized treatment in case of nudging
262	!
263	! 1380 2014-04-28 12:40:45Z heinze
264	! Change order of calls for scalar prognostic quantities:
265	! ls_advec -> nudging -> subsidence since initial profiles
266	!
267	! 1374 2014-04-25 12:55:07Z raasch
268	! missing variables added to ONLY lists
269	!
270	! 1365 2014-04-22 15:03:56Z boeske
271	! Calls of ls_advec for large scale advection added,
272	! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
273	! new argument ls_index added to the calls of subsidence
274	! +ls_index
275	!
276	! 1361 2014-04-16 15:17:48Z hoffmann
277	! Two-moment microphysics moved to the start of prognostic equations. This makes
278	! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
279	! Additionally, it is allowed to call the microphysics just once during the time
280	! step (not at each sub-time step).
281	!
282	! Two-moment cloud physics added for vector and accelerator optimization.
283	!
284	! 1353 2014-04-08 15:21:23Z heinze
285	! REAL constants provided with KIND-attribute
286	!
287	! 1337 2014-03-25 15:11:48Z heinze
288	! Bugfix: REAL constants provided with KIND-attribute
289	!
290	! 1332 2014-03-25 11:59:43Z suehring
291	! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
292	!
293	! 1330 2014-03-24 17:29:32Z suehring
294	! In case of SGS-particle velocity advection of TKE is also allowed with
295	! dissipative 5th-order scheme.
296	!
297	! 1320 2014-03-20 08:40:49Z raasch
298	! ONLY-attribute added to USE-statements,
299	! kind-parameters added to all INTEGER and REAL declaration statements,
300	! kinds are defined in new module kinds,
301	! old module precision_kind is removed,
302	! revision history before 2012 removed,
303	! comment fields (!:) to be used for variable explanations added to
304	! all variable declaration statements
305	!
306	! 1318 2014-03-17 13:35:16Z raasch
307	! module interfaces removed
308	!
309	! 1257 2013-11-08 15:18:40Z raasch
310	! openacc loop vector clauses removed, independent clauses added
311	!
312	! 1246 2013-11-01 08:59:45Z heinze
313	! enable nudging also for accelerator version
314	!
315	! 1241 2013-10-30 11:36:58Z heinze
316	! usage of nudging enabled (so far not implemented for accelerator version)
317	!
318	! 1179 2013-06-14 05:57:58Z raasch
319	! two arguments removed from routine buoyancy, ref_state updated on device
320	!
321	! 1128 2013-04-12 06:19:32Z raasch
322	! those parts requiring global communication moved to time_integration,
323	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
324	! j_north
325	!
326	! 1115 2013-03-26 18:16:16Z hoffmann
327	! optimized cloud physics: calculation of microphysical tendencies transfered
328	! to microphysics.f90; qr and nr are only calculated if precipitation is required
329	!
330	! 1111 2013-03-08 23:54:10Z raasch
331	! update directives for prognostic quantities removed
332	!
333	! 1106 2013-03-04 05:31:38Z raasch
334	! small changes in code formatting
335	!
336	! 1092 2013-02-02 11:24:22Z raasch
337	! unused variables removed
338	!
339	! 1053 2012-11-13 17:11:03Z hoffmann
340	! implementation of two new prognostic equations for rain drop concentration (nr)
341	! and rain water content (qr)
342	!
343	! currently, only available for cache loop optimization
344	!
345	! 1036 2012-10-22 13:43:42Z raasch
346	! code put under GPL (PALM 3.9)
347	!
348	! 1019 2012-09-28 06:46:45Z raasch
349	! non-optimized version of prognostic_equations removed
350	!
351	! 1015 2012-09-27 09:23:24Z raasch
352	! new branch prognostic_equations_acc
353	! OpenACC statements added + code changes required for GPU optimization
354	!
355	! 1001 2012-09-13 14:08:46Z raasch
356	! all actions concerning leapfrog- and upstream-spline-scheme removed
357	!
358	! 978 2012-08-09 08:28:32Z fricke
359	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
360	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
361	! boundary in case of non-cyclic lateral boundaries
362	! Bugfix: first thread index changes for WS-scheme at the inflow
363	!
364	! 940 2012-07-09 14:31:00Z raasch
365	! temperature equation can be switched off
366	!
367	! Revision 1.1 2000/04/13 14:56:27 schroeter
368	! Initial revision
369	!
370	!
371	! Description:
372	! ------------
373	!> Solving the prognostic equations.
374	!------------------------------------------------------------------------------!
375	MODULE prognostic_equations_mod
376
377	USE advec_s_bc_mod, &
378	ONLY: advec_s_bc
379
380	USE advec_s_pw_mod, &
381	ONLY: advec_s_pw
382
383	USE advec_s_up_mod, &
384	ONLY: advec_s_up
385
386	USE advec_u_pw_mod, &
387	ONLY: advec_u_pw
388
389	USE advec_u_up_mod, &
390	ONLY: advec_u_up
391
392	USE advec_v_pw_mod, &
393	ONLY: advec_v_pw
394
395	USE advec_v_up_mod, &
396	ONLY: advec_v_up
397
398	USE advec_w_pw_mod, &
399	ONLY: advec_w_pw
400
401	USE advec_w_up_mod, &
402	ONLY: advec_w_up
403
404	USE advec_ws, &
405	ONLY: advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
406
407	USE arrays_3d, &
408	ONLY: diss_l_e, diss_l_pt, diss_l_q, &
409	diss_l_s, diss_l_sa, diss_s_e, &
410	diss_s_pt, diss_s_q, diss_s_s, diss_s_sa, &
411	e, e_p, flux_s_e, flux_s_pt, flux_s_q, &
412	flux_s_s, flux_s_sa, flux_l_e, &
413	flux_l_pt, flux_l_q, flux_l_s, &
414	flux_l_sa, pt, ptdf_x, ptdf_y, pt_init, &
415	pt_p, prho, q, q_init, q_p, rdf, rdf_sc, &
416	ref_state, rho_ocean, s, s_init, s_p, tend, te_m, &
417	tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u, &
418	ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
419
420	USE buoyancy_mod, &
421	ONLY: buoyancy
422
423	USE control_parameters, &
424	ONLY: bc_dirichlet_l, &
425	bc_dirichlet_n, &
426	bc_dirichlet_r, &
427	bc_dirichlet_s, &
428	bc_radiation_l, &
429	bc_radiation_n, &
430	bc_radiation_r, &
431	bc_radiation_s, &
432	constant_diffusion, &
433	debug_output_timestep, &
434	dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, &
435	humidity, intermediate_timestep_count, &
436	intermediate_timestep_count_max, large_scale_forcing, &
437	large_scale_subsidence, &
438	monotonic_limiter_z, &
439	neutral, nudging, &
440	ocean_mode, passive_scalar, plant_canopy, pt_reference, &
441	scalar_advec, scalar_advec, simulated_time, sloping_surface, &
442	timestep_scheme, tsc, use_subsidence_tendencies, &
443	use_upstream_for_tke, wind_turbine, ws_scheme_mom, &
444	ws_scheme_sca, urban_surface, land_surface, &
445	time_since_reference_point, salsa
446
447	USE coriolis_mod, &
448	ONLY: coriolis
449
450	USE cpulog, &
451	ONLY: cpu_log, log_point, log_point_s
452
453	USE diffusion_s_mod, &
454	ONLY: diffusion_s
455
456	USE diffusion_u_mod, &
457	ONLY: diffusion_u
458
459	USE diffusion_v_mod, &
460	ONLY: diffusion_v
461
462	USE diffusion_w_mod, &
463	ONLY: diffusion_w
464
465	USE indices, &
466	ONLY: advc_flags_s, &
467	nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
468	nzb, nzt, wall_flags_0
469
470	USE kinds
471
472	USE lsf_nudging_mod, &
473	ONLY: ls_advec, nudge
474
475	USE module_interface, &
476	ONLY: module_interface_actions, &
477	module_interface_non_advective_processes, &
478	module_interface_exchange_horiz, &
479	module_interface_prognostic_equations
480
481	USE ocean_mod, &
482	ONLY: stokes_drift_terms, stokes_force, &
483	wave_breaking, wave_breaking_term
484
485	USE plant_canopy_model_mod, &
486	ONLY: cthf, pcm_tendency
487
488	#if defined( __rrtmg )
489	USE radiation_model_mod, &
490	ONLY: radiation, radiation_tendency, &
491	skip_time_do_radiation
492	#endif
493
494	USE statistics, &
495	ONLY: hom
496
497	USE subsidence_mod, &
498	ONLY: subsidence
499
500	USE surface_mod, &
501	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
502	surf_usm_v
503
504	IMPLICIT NONE
505
506	PRIVATE
507	PUBLIC prognostic_equations_cache, prognostic_equations_vector
508
509	INTERFACE prognostic_equations_cache
510	MODULE PROCEDURE prognostic_equations_cache
511	END INTERFACE prognostic_equations_cache
512
513	INTERFACE prognostic_equations_vector
514	MODULE PROCEDURE prognostic_equations_vector
515	END INTERFACE prognostic_equations_vector
516
517
518	CONTAINS
519
520
521	!------------------------------------------------------------------------------!
522	! Description:
523	! ------------
524	!> Version with one optimized loop over all equations. It is only allowed to
525	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
526	!>
527	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
528	!> so communication between CPUs is not allowed (does not make sense) within
529	!> these loops.
530	!>
531	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
532	!------------------------------------------------------------------------------!
533
534	SUBROUTINE prognostic_equations_cache
535
536
537	INTEGER(iwp) :: i !<
538	INTEGER(iwp) :: i_omp_start !<
539	INTEGER(iwp) :: j !<
540	INTEGER(iwp) :: k !<
541	!$ INTEGER(iwp) :: omp_get_thread_num !<
542	INTEGER(iwp) :: tn = 0 !<
543
544	LOGICAL :: loop_start !<
545
546
547	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'start' )
548	!
549	!-- Time measurement can only be performed for the whole set of equations
550	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
551
552	!$OMP PARALLEL PRIVATE (i,j)
553	!$OMP DO
554	DO i = nxl, nxr
555	DO j = nys, nyn
556	!
557	!-- Calculate non advective processes for all other modules
558	CALL module_interface_non_advective_processes( i, j )
559	ENDDO
560	ENDDO
561	!
562	!-- Module Inferface for exchange horiz after non_advective_processes but before
563	!-- advection. Therefore, non_advective_processes must not run for ghost points.
564	!$OMP END PARALLEL
565	CALL module_interface_exchange_horiz()
566	!
567	!-- Loop over all prognostic equations
568	!$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
569
570	!$ tn = omp_get_thread_num()
571	loop_start = .TRUE.
572
573	!$OMP DO
574	DO i = nxl, nxr
575
576	!
577	!-- Store the first loop index. It differs for each thread and is required
578	!-- later in advec_ws
579	IF ( loop_start ) THEN
580	loop_start = .FALSE.
581	i_omp_start = i
582	ENDIF
583
584	DO j = nys, nyn
585	!
586	!-- Tendency terms for u-velocity component. Please note, in case of
587	!-- non-cyclic boundary conditions the grid point i=0 is excluded from
588	!-- the prognostic equations for the u-component.
589	IF ( i >= nxlu ) THEN
590
591	tend(:,j,i) = 0.0_wp
592	IF ( timestep_scheme(1:5) == 'runge' ) THEN
593	IF ( ws_scheme_mom ) THEN
594	CALL advec_u_ws( i, j, i_omp_start, tn )
595	ELSE
596	CALL advec_u_pw( i, j )
597	ENDIF
598	ELSE
599	CALL advec_u_up( i, j )
600	ENDIF
601	CALL diffusion_u( i, j )
602	CALL coriolis( i, j, 1 )
603	IF ( sloping_surface .AND. .NOT. neutral ) THEN
604	CALL buoyancy( i, j, pt, 1 )
605	ENDIF
606
607	!
608	!-- Drag by plant canopy
609	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
610
611	!
612	!-- External pressure gradient
613	IF ( dp_external ) THEN
614	DO k = dp_level_ind_b+1, nzt
615	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
616	ENDDO
617	ENDIF
618
619	!
620	!-- Nudging
621	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
622
623	!
624	!-- Effect of Stokes drift (in ocean mode only)
625	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 )
626
627	CALL module_interface_actions( i, j, 'u-tendency' )
628	!
629	!-- Prognostic equation for u-velocity component
630	DO k = nzb+1, nzt
631
632	u_p(k,j,i) = u(k,j,i) + ( dt_3d * &
633	( tsc(2) * tend(k,j,i) + &
634	tsc(3) * tu_m(k,j,i) ) &
635	- tsc(5) * rdf(k) &
636	* ( u(k,j,i) - u_init(k) ) &
637	) * MERGE( 1.0_wp, 0.0_wp, &
638	BTEST( wall_flags_0(k,j,i), 1 )&
639	)
640	ENDDO
641
642	!
643	!-- Add turbulence generated by wave breaking (in ocean mode only)
644	IF ( wave_breaking .AND. &
645	intermediate_timestep_count == intermediate_timestep_count_max )&
646	THEN
647	CALL wave_breaking_term( i, j, 1 )
648	ENDIF
649
650	!
651	!-- Calculate tendencies for the next Runge-Kutta step
652	IF ( timestep_scheme(1:5) == 'runge' ) THEN
653	IF ( intermediate_timestep_count == 1 ) THEN
654	DO k = nzb+1, nzt
655	tu_m(k,j,i) = tend(k,j,i)
656	ENDDO
657	ELSEIF ( intermediate_timestep_count < &
658	intermediate_timestep_count_max ) THEN
659	DO k = nzb+1, nzt
660	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
661	+ 5.3125_wp * tu_m(k,j,i)
662	ENDDO
663	ENDIF
664	ENDIF
665
666	ENDIF
667	!
668	!-- Tendency terms for v-velocity component. Please note, in case of
669	!-- non-cyclic boundary conditions the grid point j=0 is excluded from
670	!-- the prognostic equations for the v-component. !--
671	IF ( j >= nysv ) THEN
672
673	tend(:,j,i) = 0.0_wp
674	IF ( timestep_scheme(1:5) == 'runge' ) THEN
675	IF ( ws_scheme_mom ) THEN
676	CALL advec_v_ws( i, j, i_omp_start, tn )
677	ELSE
678	CALL advec_v_pw( i, j )
679	ENDIF
680	ELSE
681	CALL advec_v_up( i, j )
682	ENDIF
683	CALL diffusion_v( i, j )
684	CALL coriolis( i, j, 2 )
685
686	!
687	!-- Drag by plant canopy
688	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
689
690	!
691	!-- External pressure gradient
692	IF ( dp_external ) THEN
693	DO k = dp_level_ind_b+1, nzt
694	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
695	ENDDO
696	ENDIF
697
698	!
699	!-- Nudging
700	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
701
702	!
703	!-- Effect of Stokes drift (in ocean mode only)
704	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 )
705
706	CALL module_interface_actions( i, j, 'v-tendency' )
707	!
708	!-- Prognostic equation for v-velocity component
709	DO k = nzb+1, nzt
710	v_p(k,j,i) = v(k,j,i) + ( dt_3d * &
711	( tsc(2) * tend(k,j,i) + &
712	tsc(3) * tv_m(k,j,i) ) &
713	- tsc(5) * rdf(k) &
714	* ( v(k,j,i) - v_init(k) )&
715	) * MERGE( 1.0_wp, 0.0_wp, &
716	BTEST( wall_flags_0(k,j,i), 2 )&
717	)
718	ENDDO
719
720	!
721	!-- Add turbulence generated by wave breaking (in ocean mode only)
722	IF ( wave_breaking .AND. &
723	intermediate_timestep_count == intermediate_timestep_count_max )&
724	THEN
725	CALL wave_breaking_term( i, j, 2 )
726	ENDIF
727
728	!
729	!-- Calculate tendencies for the next Runge-Kutta step
730	IF ( timestep_scheme(1:5) == 'runge' ) THEN
731	IF ( intermediate_timestep_count == 1 ) THEN
732	DO k = nzb+1, nzt
733	tv_m(k,j,i) = tend(k,j,i)
734	ENDDO
735	ELSEIF ( intermediate_timestep_count < &
736	intermediate_timestep_count_max ) THEN
737	DO k = nzb+1, nzt
738	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
739	+ 5.3125_wp * tv_m(k,j,i)
740	ENDDO
741	ENDIF
742	ENDIF
743
744	ENDIF
745
746	!
747	!-- Tendency terms for w-velocity component
748	tend(:,j,i) = 0.0_wp
749	IF ( timestep_scheme(1:5) == 'runge' ) THEN
750	IF ( ws_scheme_mom ) THEN
751	CALL advec_w_ws( i, j, i_omp_start, tn )
752	ELSE
753	CALL advec_w_pw( i, j )
754	END IF
755	ELSE
756	CALL advec_w_up( i, j )
757	ENDIF
758	CALL diffusion_w( i, j )
759	CALL coriolis( i, j, 3 )
760
761	IF ( .NOT. neutral ) THEN
762	IF ( ocean_mode ) THEN
763	CALL buoyancy( i, j, rho_ocean, 3 )
764	ELSE
765	IF ( .NOT. humidity ) THEN
766	CALL buoyancy( i, j, pt, 3 )
767	ELSE
768	CALL buoyancy( i, j, vpt, 3 )
769	ENDIF
770	ENDIF
771	ENDIF
772
773	!
774	!-- Drag by plant canopy
775	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
776
777	!
778	!-- Effect of Stokes drift (in ocean mode only)
779	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 3 )
780
781	CALL module_interface_actions( i, j, 'w-tendency' )
782	!
783	!-- Prognostic equation for w-velocity component
784	DO k = nzb+1, nzt-1
785	w_p(k,j,i) = w(k,j,i) + ( dt_3d * &
786	( tsc(2) * tend(k,j,i) + &
787	tsc(3) * tw_m(k,j,i) ) &
788	- tsc(5) * rdf(k) * w(k,j,i) &
789	) * MERGE( 1.0_wp, 0.0_wp, &
790	BTEST( wall_flags_0(k,j,i), 3 )&
791	)
792	ENDDO
793
794	!
795	!-- Calculate tendencies for the next Runge-Kutta step
796	IF ( timestep_scheme(1:5) == 'runge' ) THEN
797	IF ( intermediate_timestep_count == 1 ) THEN
798	DO k = nzb+1, nzt-1
799	tw_m(k,j,i) = tend(k,j,i)
800	ENDDO
801	ELSEIF ( intermediate_timestep_count < &
802	intermediate_timestep_count_max ) THEN
803	DO k = nzb+1, nzt-1
804	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
805	+ 5.3125_wp * tw_m(k,j,i)
806	ENDDO
807	ENDIF
808	ENDIF
809
810	!
811	!-- If required, compute prognostic equation for potential temperature
812	IF ( .NOT. neutral ) THEN
813	!
814	!-- Tendency terms for potential temperature
815	tend(:,j,i) = 0.0_wp
816	IF ( timestep_scheme(1:5) == 'runge' ) THEN
817	IF ( ws_scheme_sca ) THEN
818	CALL advec_s_ws( advc_flags_s, &
819	i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
820	flux_l_pt, diss_l_pt, i_omp_start, tn, &
821	bc_dirichlet_l .OR. bc_radiation_l, &
822	bc_dirichlet_n .OR. bc_radiation_n, &
823	bc_dirichlet_r .OR. bc_radiation_r, &
824	bc_dirichlet_s .OR. bc_radiation_s )
825	ELSE
826	CALL advec_s_pw( i, j, pt )
827	ENDIF
828	ELSE
829	CALL advec_s_up( i, j, pt )
830	ENDIF
831	CALL diffusion_s( i, j, pt, &
832	surf_def_h(0)%shf, surf_def_h(1)%shf, &
833	surf_def_h(2)%shf, &
834	surf_lsm_h%shf, surf_usm_h%shf, &
835	surf_def_v(0)%shf, surf_def_v(1)%shf, &
836	surf_def_v(2)%shf, surf_def_v(3)%shf, &
837	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
838	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
839	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
840	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
841
842	!
843	!-- Consideration of heat sources within the plant canopy
844	IF ( plant_canopy .AND. &
845	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
846	CALL pcm_tendency( i, j, 4 )
847	ENDIF
848
849	!
850	!-- Large scale advection
851	IF ( large_scale_forcing ) THEN
852	CALL ls_advec( i, j, simulated_time, 'pt' )
853	ENDIF
854
855	!
856	!-- Nudging
857	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
858
859	!
860	!-- If required, compute effect of large-scale subsidence/ascent
861	IF ( large_scale_subsidence .AND. &
862	.NOT. use_subsidence_tendencies ) THEN
863	CALL subsidence( i, j, tend, pt, pt_init, 2 )
864	ENDIF
865
866	#if defined( __rrtmg )
867	!
868	!-- If required, add tendency due to radiative heating/cooling
869	IF ( radiation .AND. &
870	simulated_time > skip_time_do_radiation ) THEN
871	CALL radiation_tendency ( i, j, tend )
872	ENDIF
873	#endif
874
875	CALL module_interface_actions( i, j, 'pt-tendency' )
876	!
877	!-- Prognostic equation for potential temperature
878	DO k = nzb+1, nzt
879	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * &
880	( tsc(2) * tend(k,j,i) + &
881	tsc(3) * tpt_m(k,j,i) ) &
882	- tsc(5) &
883	* ( pt(k,j,i) - pt_init(k) ) &
884	* ( rdf_sc(k) + ptdf_x(i) &
885	+ ptdf_y(j) ) &
886	) &
887	* MERGE( 1.0_wp, 0.0_wp, &
888	BTEST( wall_flags_0(k,j,i), 0 )&
889	)
890	ENDDO
891
892	!
893	!-- Calculate tendencies for the next Runge-Kutta step
894	IF ( timestep_scheme(1:5) == 'runge' ) THEN
895	IF ( intermediate_timestep_count == 1 ) THEN
896	DO k = nzb+1, nzt
897	tpt_m(k,j,i) = tend(k,j,i)
898	ENDDO
899	ELSEIF ( intermediate_timestep_count < &
900	intermediate_timestep_count_max ) THEN
901	DO k = nzb+1, nzt
902	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
903	5.3125_wp * tpt_m(k,j,i)
904	ENDDO
905	ENDIF
906	ENDIF
907
908	ENDIF
909
910	!
911	!-- If required, compute prognostic equation for total water content
912	IF ( humidity ) THEN
913
914	!
915	!-- Tendency-terms for total water content / scalar
916	tend(:,j,i) = 0.0_wp
917	IF ( timestep_scheme(1:5) == 'runge' ) &
918	THEN
919	IF ( ws_scheme_sca ) THEN
920	CALL advec_s_ws( advc_flags_s, &
921	i, j, q, 'q', flux_s_q, &
922	diss_s_q, flux_l_q, diss_l_q, &
923	i_omp_start, tn, &
924	bc_dirichlet_l .OR. bc_radiation_l, &
925	bc_dirichlet_n .OR. bc_radiation_n, &
926	bc_dirichlet_r .OR. bc_radiation_r, &
927	bc_dirichlet_s .OR. bc_radiation_s )
928	ELSE
929	CALL advec_s_pw( i, j, q )
930	ENDIF
931	ELSE
932	CALL advec_s_up( i, j, q )
933	ENDIF
934	CALL diffusion_s( i, j, q, &
935	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
936	surf_def_h(2)%qsws, &
937	surf_lsm_h%qsws, surf_usm_h%qsws, &
938	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
939	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
940	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
941	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
942	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
943	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
944
945	!
946	!-- Sink or source of humidity due to canopy elements
947	IF ( plant_canopy ) CALL pcm_tendency( i, j, 5 )
948
949	!
950	!-- Large scale advection
951	IF ( large_scale_forcing ) THEN
952	CALL ls_advec( i, j, simulated_time, 'q' )
953	ENDIF
954
955	!
956	!-- Nudging
957	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
958
959	!
960	!-- If required compute influence of large-scale subsidence/ascent
961	IF ( large_scale_subsidence .AND. &
962	.NOT. use_subsidence_tendencies ) THEN
963	CALL subsidence( i, j, tend, q, q_init, 3 )
964	ENDIF
965
966	CALL module_interface_actions( i, j, 'q-tendency' )
967
968	!
969	!-- Prognostic equation for total water content / scalar
970	DO k = nzb+1, nzt
971	q_p(k,j,i) = q(k,j,i) + ( dt_3d * &
972	( tsc(2) * tend(k,j,i) + &
973	tsc(3) * tq_m(k,j,i) ) &
974	- tsc(5) * rdf_sc(k) * &
975	( q(k,j,i) - q_init(k) ) &
976	) &
977	* MERGE( 1.0_wp, 0.0_wp, &
978	BTEST( wall_flags_0(k,j,i), 0 )&
979	)
980	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
981	ENDDO
982
983	!
984	!-- Calculate tendencies for the next Runge-Kutta step
985	IF ( timestep_scheme(1:5) == 'runge' ) THEN
986	IF ( intermediate_timestep_count == 1 ) THEN
987	DO k = nzb+1, nzt
988	tq_m(k,j,i) = tend(k,j,i)
989	ENDDO
990	ELSEIF ( intermediate_timestep_count < &
991	intermediate_timestep_count_max ) THEN
992	DO k = nzb+1, nzt
993	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
994	5.3125_wp * tq_m(k,j,i)
995	ENDDO
996	ENDIF
997	ENDIF
998
999	ENDIF
1000
1001	!
1002	!-- If required, compute prognostic equation for scalar
1003	IF ( passive_scalar ) THEN
1004	!
1005	!-- Tendency-terms for total water content / scalar
1006	tend(:,j,i) = 0.0_wp
1007	IF ( timestep_scheme(1:5) == 'runge' ) &
1008	THEN
1009	IF ( ws_scheme_sca ) THEN
1010	!
1011	!-- For scalar advection apply monotonic flux limiter near
1012	!-- topography.
1013	CALL advec_s_ws( advc_flags_s, &
1014	i, j, s, 's', flux_s_s, &
1015	diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
1016	tn, &
1017	bc_dirichlet_l .OR. bc_radiation_l, &
1018	bc_dirichlet_n .OR. bc_radiation_n, &
1019	bc_dirichlet_r .OR. bc_radiation_r, &
1020	bc_dirichlet_s .OR. bc_radiation_s, &
1021	monotonic_limiter_z )
1022	ELSE
1023	CALL advec_s_pw( i, j, s )
1024	ENDIF
1025	ELSE
1026	CALL advec_s_up( i, j, s )
1027	ENDIF
1028	CALL diffusion_s( i, j, s, &
1029	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1030	surf_def_h(2)%ssws, &
1031	surf_lsm_h%ssws, surf_usm_h%ssws, &
1032	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1033	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1034	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1035	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1036	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1037	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1038
1039	!
1040	!-- Sink or source of scalar concentration due to canopy elements
1041	IF ( plant_canopy ) CALL pcm_tendency( i, j, 7 )
1042
1043	!
1044	!-- Large scale advection, still need to be extended for scalars
1045	! IF ( large_scale_forcing ) THEN
1046	! CALL ls_advec( i, j, simulated_time, 's' )
1047	! ENDIF
1048
1049	!
1050	!-- Nudging, still need to be extended for scalars
1051	! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' )
1052
1053	!
1054	!-- If required compute influence of large-scale subsidence/ascent.
1055	!-- Note, the last argument is of no meaning in this case, as it is
1056	!-- only used in conjunction with large_scale_forcing, which is to
1057	!-- date not implemented for scalars.
1058	IF ( large_scale_subsidence .AND. &
1059	.NOT. use_subsidence_tendencies .AND. &
1060	.NOT. large_scale_forcing ) THEN
1061	CALL subsidence( i, j, tend, s, s_init, 3 )
1062	ENDIF
1063
1064	CALL module_interface_actions( i, j, 's-tendency' )
1065
1066	!
1067	!-- Prognostic equation for scalar
1068	DO k = nzb+1, nzt
1069	s_p(k,j,i) = s(k,j,i) + ( dt_3d * &
1070	( tsc(2) * tend(k,j,i) + &
1071	tsc(3) * ts_m(k,j,i) ) &
1072	- tsc(5) * rdf_sc(k) &
1073	* ( s(k,j,i) - s_init(k) )&
1074	) &
1075	* MERGE( 1.0_wp, 0.0_wp, &
1076	BTEST( wall_flags_0(k,j,i), 0 )&
1077	)
1078	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1079	ENDDO
1080
1081	!
1082	!-- Calculate tendencies for the next Runge-Kutta step
1083	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1084	IF ( intermediate_timestep_count == 1 ) THEN
1085	DO k = nzb+1, nzt
1086	ts_m(k,j,i) = tend(k,j,i)
1087	ENDDO
1088	ELSEIF ( intermediate_timestep_count < &
1089	intermediate_timestep_count_max ) THEN
1090	DO k = nzb+1, nzt
1091	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1092	5.3125_wp * ts_m(k,j,i)
1093	ENDDO
1094	ENDIF
1095	ENDIF
1096
1097	ENDIF
1098	!
1099	!-- Calculate prognostic equations for all other modules
1100	CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1101
1102	ENDDO ! loop over j
1103	ENDDO ! loop over i
1104	!$OMP END PARALLEL
1105
1106
1107	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1108
1109	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'end' )
1110
1111	END SUBROUTINE prognostic_equations_cache
1112
1113
1114	!------------------------------------------------------------------------------!
1115	! Description:
1116	! ------------
1117	!> Version for vector machines
1118	!------------------------------------------------------------------------------!
1119
1120	SUBROUTINE prognostic_equations_vector
1121
1122
1123	INTEGER(iwp) :: i !<
1124	INTEGER(iwp) :: j !<
1125	INTEGER(iwp) :: k !<
1126
1127	REAL(wp) :: sbt !<
1128
1129
1130	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'start' )
1131	!
1132	!-- Calculate non advective processes for all other modules
1133	CALL module_interface_non_advective_processes
1134	!
1135	!-- Module Inferface for exchange horiz after non_advective_processes but before
1136	!-- advection. Therefore, non_advective_processes must not run for ghost points.
1137	CALL module_interface_exchange_horiz()
1138	!
1139	!-- u-velocity component
1140	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1141
1142	!$ACC KERNELS PRESENT(tend)
1143	tend = 0.0_wp
1144	!$ACC END KERNELS
1145	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1146	IF ( ws_scheme_mom ) THEN
1147	CALL advec_u_ws
1148	ELSE
1149	CALL advec_u_pw
1150	ENDIF
1151	ELSE
1152	CALL advec_u_up
1153	ENDIF
1154	CALL diffusion_u
1155	CALL coriolis( 1 )
1156	IF ( sloping_surface .AND. .NOT. neutral ) THEN
1157	CALL buoyancy( pt, 1 )
1158	ENDIF
1159
1160	!
1161	!-- Drag by plant canopy
1162	IF ( plant_canopy ) CALL pcm_tendency( 1 )
1163
1164	!
1165	!-- External pressure gradient
1166	IF ( dp_external ) THEN
1167	DO i = nxlu, nxr
1168	DO j = nys, nyn
1169	DO k = dp_level_ind_b+1, nzt
1170	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1171	ENDDO
1172	ENDDO
1173	ENDDO
1174	ENDIF
1175
1176	!
1177	!-- Nudging
1178	IF ( nudging ) CALL nudge( simulated_time, 'u' )
1179
1180	!
1181	!-- Effect of Stokes drift (in ocean mode only)
1182	IF ( stokes_force ) CALL stokes_drift_terms( 1 )
1183
1184	CALL module_interface_actions( 'u-tendency' )
1185
1186	!
1187	!-- Prognostic equation for u-velocity component
1188	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1189	!$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1190	!$ACC PRESENT(tsc(2:5)) &
1191	!$ACC PRESENT(u_p)
1192	DO i = nxlu, nxr
1193	DO j = nys, nyn
1194	DO k = nzb+1, nzt
1195	u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1196	tsc(3) * tu_m(k,j,i) ) &
1197	- tsc(5) * rdf(k) * &
1198	( u(k,j,i) - u_init(k) ) &
1199	) * MERGE( 1.0_wp, 0.0_wp, &
1200	BTEST( wall_flags_0(k,j,i), 1 ) &
1201	)
1202	ENDDO
1203	ENDDO
1204	ENDDO
1205
1206	!
1207	!-- Add turbulence generated by wave breaking (in ocean mode only)
1208	IF ( wave_breaking .AND. &
1209	intermediate_timestep_count == intermediate_timestep_count_max ) &
1210	THEN
1211	CALL wave_breaking_term( 1 )
1212	ENDIF
1213
1214	!
1215	!-- Calculate tendencies for the next Runge-Kutta step
1216	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1217	IF ( intermediate_timestep_count == 1 ) THEN
1218	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1219	!$ACC PRESENT(tend, tu_m)
1220	DO i = nxlu, nxr
1221	DO j = nys, nyn
1222	DO k = nzb+1, nzt
1223	tu_m(k,j,i) = tend(k,j,i)
1224	ENDDO
1225	ENDDO
1226	ENDDO
1227	ELSEIF ( intermediate_timestep_count < &
1228	intermediate_timestep_count_max ) THEN
1229	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1230	!$ACC PRESENT(tend, tu_m)
1231	DO i = nxlu, nxr
1232	DO j = nys, nyn
1233	DO k = nzb+1, nzt
1234	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1235	+ 5.3125_wp * tu_m(k,j,i)
1236	ENDDO
1237	ENDDO
1238	ENDDO
1239	ENDIF
1240	ENDIF
1241
1242	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1243
1244	!
1245	!-- v-velocity component
1246	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1247
1248	!$ACC KERNELS PRESENT(tend)
1249	tend = 0.0_wp
1250	!$ACC END KERNELS
1251	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1252	IF ( ws_scheme_mom ) THEN
1253	CALL advec_v_ws
1254	ELSE
1255	CALL advec_v_pw
1256	END IF
1257	ELSE
1258	CALL advec_v_up
1259	ENDIF
1260	CALL diffusion_v
1261	CALL coriolis( 2 )
1262
1263	!
1264	!-- Drag by plant canopy
1265	IF ( plant_canopy ) CALL pcm_tendency( 2 )
1266
1267	!
1268	!-- External pressure gradient
1269	IF ( dp_external ) THEN
1270	DO i = nxl, nxr
1271	DO j = nysv, nyn
1272	DO k = dp_level_ind_b+1, nzt
1273	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1274	ENDDO
1275	ENDDO
1276	ENDDO
1277	ENDIF
1278
1279	!
1280	!-- Nudging
1281	IF ( nudging ) CALL nudge( simulated_time, 'v' )
1282
1283	!
1284	!-- Effect of Stokes drift (in ocean mode only)
1285	IF ( stokes_force ) CALL stokes_drift_terms( 2 )
1286
1287	CALL module_interface_actions( 'v-tendency' )
1288
1289	!
1290	!-- Prognostic equation for v-velocity component
1291	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1292	!$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1293	!$ACC PRESENT(tsc(2:5)) &
1294	!$ACC PRESENT(v_p)
1295	DO i = nxl, nxr
1296	DO j = nysv, nyn
1297	DO k = nzb+1, nzt
1298	v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1299	tsc(3) * tv_m(k,j,i) ) &
1300	- tsc(5) * rdf(k) * &
1301	( v(k,j,i) - v_init(k) ) &
1302	) * MERGE( 1.0_wp, 0.0_wp, &
1303	BTEST( wall_flags_0(k,j,i), 2 )&
1304	)
1305	ENDDO
1306	ENDDO
1307	ENDDO
1308
1309	!
1310	!-- Add turbulence generated by wave breaking (in ocean mode only)
1311	IF ( wave_breaking .AND. &
1312	intermediate_timestep_count == intermediate_timestep_count_max ) &
1313	THEN
1314	CALL wave_breaking_term( 2 )
1315	ENDIF
1316
1317	!
1318	!-- Calculate tendencies for the next Runge-Kutta step
1319	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1320	IF ( intermediate_timestep_count == 1 ) THEN
1321	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1322	!$ACC PRESENT(tend, tv_m)
1323	DO i = nxl, nxr
1324	DO j = nysv, nyn
1325	DO k = nzb+1, nzt
1326	tv_m(k,j,i) = tend(k,j,i)
1327	ENDDO
1328	ENDDO
1329	ENDDO
1330	ELSEIF ( intermediate_timestep_count < &
1331	intermediate_timestep_count_max ) THEN
1332	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1333	!$ACC PRESENT(tend, tv_m)
1334	DO i = nxl, nxr
1335	DO j = nysv, nyn
1336	DO k = nzb+1, nzt
1337	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1338	+ 5.3125_wp * tv_m(k,j,i)
1339	ENDDO
1340	ENDDO
1341	ENDDO
1342	ENDIF
1343	ENDIF
1344
1345	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1346
1347	!
1348	!-- w-velocity component
1349	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1350
1351	!$ACC KERNELS PRESENT(tend)
1352	tend = 0.0_wp
1353	!$ACC END KERNELS
1354	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1355	IF ( ws_scheme_mom ) THEN
1356	CALL advec_w_ws
1357	ELSE
1358	CALL advec_w_pw
1359	ENDIF
1360	ELSE
1361	CALL advec_w_up
1362	ENDIF
1363	CALL diffusion_w
1364	CALL coriolis( 3 )
1365
1366	IF ( .NOT. neutral ) THEN
1367	IF ( ocean_mode ) THEN
1368	CALL buoyancy( rho_ocean, 3 )
1369	ELSE
1370	IF ( .NOT. humidity ) THEN
1371	CALL buoyancy( pt, 3 )
1372	ELSE
1373	CALL buoyancy( vpt, 3 )
1374	ENDIF
1375	ENDIF
1376	ENDIF
1377
1378	!
1379	!-- Drag by plant canopy
1380	IF ( plant_canopy ) CALL pcm_tendency( 3 )
1381
1382	!
1383	!-- Effect of Stokes drift (in ocean mode only)
1384	IF ( stokes_force ) CALL stokes_drift_terms( 3 )
1385
1386	CALL module_interface_actions( 'w-tendency' )
1387
1388	!
1389	!-- Prognostic equation for w-velocity component
1390	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1391	!$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1392	!$ACC PRESENT(tsc(2:5)) &
1393	!$ACC PRESENT(w_p)
1394	DO i = nxl, nxr
1395	DO j = nys, nyn
1396	DO k = nzb+1, nzt-1
1397	w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1398	tsc(3) * tw_m(k,j,i) ) &
1399	- tsc(5) * rdf(k) * w(k,j,i) &
1400	) * MERGE( 1.0_wp, 0.0_wp, &
1401	BTEST( wall_flags_0(k,j,i), 3 )&
1402	)
1403	ENDDO
1404	ENDDO
1405	ENDDO
1406
1407	!
1408	!-- Calculate tendencies for the next Runge-Kutta step
1409	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1410	IF ( intermediate_timestep_count == 1 ) THEN
1411	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1412	!$ACC PRESENT(tend, tw_m)
1413	DO i = nxl, nxr
1414	DO j = nys, nyn
1415	DO k = nzb+1, nzt-1
1416	tw_m(k,j,i) = tend(k,j,i)
1417	ENDDO
1418	ENDDO
1419	ENDDO
1420	ELSEIF ( intermediate_timestep_count < &
1421	intermediate_timestep_count_max ) THEN
1422	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1423	!$ACC PRESENT(tend, tw_m)
1424	DO i = nxl, nxr
1425	DO j = nys, nyn
1426	DO k = nzb+1, nzt-1
1427	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1428	+ 5.3125_wp * tw_m(k,j,i)
1429	ENDDO
1430	ENDDO
1431	ENDDO
1432	ENDIF
1433	ENDIF
1434
1435	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1436
1437
1438	!
1439	!-- If required, compute prognostic equation for potential temperature
1440	IF ( .NOT. neutral ) THEN
1441
1442	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1443
1444	!
1445	!-- pt-tendency terms with communication
1446	sbt = tsc(2)
1447	IF ( scalar_advec == 'bc-scheme' ) THEN
1448
1449	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1450	!
1451	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1452	sbt = 1.0_wp
1453	ENDIF
1454	tend = 0.0_wp
1455	CALL advec_s_bc( pt, 'pt' )
1456
1457	ENDIF
1458
1459	!
1460	!-- pt-tendency terms with no communication
1461	IF ( scalar_advec /= 'bc-scheme' ) THEN
1462	!$ACC KERNELS PRESENT(tend)
1463	tend = 0.0_wp
1464	!$ACC END KERNELS
1465	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1466	IF ( ws_scheme_sca ) THEN
1467	CALL advec_s_ws( advc_flags_s, pt, 'pt', &
1468	bc_dirichlet_l .OR. bc_radiation_l, &
1469	bc_dirichlet_n .OR. bc_radiation_n, &
1470	bc_dirichlet_r .OR. bc_radiation_r, &
1471	bc_dirichlet_s .OR. bc_radiation_s )
1472	ELSE
1473	CALL advec_s_pw( pt )
1474	ENDIF
1475	ELSE
1476	CALL advec_s_up( pt )
1477	ENDIF
1478	ENDIF
1479
1480	CALL diffusion_s( pt, &
1481	surf_def_h(0)%shf, surf_def_h(1)%shf, &
1482	surf_def_h(2)%shf, &
1483	surf_lsm_h%shf, surf_usm_h%shf, &
1484	surf_def_v(0)%shf, surf_def_v(1)%shf, &
1485	surf_def_v(2)%shf, surf_def_v(3)%shf, &
1486	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
1487	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
1488	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
1489	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1490
1491	!
1492	!-- Consideration of heat sources within the plant canopy
1493	IF ( plant_canopy .AND. &
1494	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
1495	CALL pcm_tendency( 4 )
1496	ENDIF
1497
1498	!
1499	!-- Large scale advection
1500	IF ( large_scale_forcing ) THEN
1501	CALL ls_advec( simulated_time, 'pt' )
1502	ENDIF
1503
1504	!
1505	!-- Nudging
1506	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
1507
1508	!
1509	!-- If required compute influence of large-scale subsidence/ascent
1510	IF ( large_scale_subsidence .AND. &
1511	.NOT. use_subsidence_tendencies ) THEN
1512	CALL subsidence( tend, pt, pt_init, 2 )
1513	ENDIF
1514
1515	#if defined( __rrtmg )
1516	!
1517	!-- If required, add tendency due to radiative heating/cooling
1518	IF ( radiation .AND. &
1519	simulated_time > skip_time_do_radiation ) THEN
1520	CALL radiation_tendency ( tend )
1521	ENDIF
1522	#endif
1523
1524	CALL module_interface_actions( 'pt-tendency' )
1525
1526	!
1527	!-- Prognostic equation for potential temperature
1528	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1529	!$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1530	!$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1531	!$ACC PRESENT(tsc(3:5)) &
1532	!$ACC PRESENT(pt_p)
1533	DO i = nxl, nxr
1534	DO j = nys, nyn
1535	DO k = nzb+1, nzt
1536	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1537	tsc(3) * tpt_m(k,j,i) ) &
1538	- tsc(5) * &
1539	( pt(k,j,i) - pt_init(k) ) *&
1540	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1541	) &
1542	* MERGE( 1.0_wp, 0.0_wp, &
1543	BTEST( wall_flags_0(k,j,i), 0 ) &
1544	)
1545	ENDDO
1546	ENDDO
1547	ENDDO
1548	!
1549	!-- Calculate tendencies for the next Runge-Kutta step
1550	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1551	IF ( intermediate_timestep_count == 1 ) THEN
1552	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1553	!$ACC PRESENT(tend, tpt_m)
1554	DO i = nxl, nxr
1555	DO j = nys, nyn
1556	DO k = nzb+1, nzt
1557	tpt_m(k,j,i) = tend(k,j,i)
1558	ENDDO
1559	ENDDO
1560	ENDDO
1561	ELSEIF ( intermediate_timestep_count < &
1562	intermediate_timestep_count_max ) THEN
1563	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1564	!$ACC PRESENT(tend, tpt_m)
1565	DO i = nxl, nxr
1566	DO j = nys, nyn
1567	DO k = nzb+1, nzt
1568	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1569	5.3125_wp * tpt_m(k,j,i)
1570	ENDDO
1571	ENDDO
1572	ENDDO
1573	ENDIF
1574	ENDIF
1575
1576	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1577
1578	ENDIF
1579
1580	!
1581	!-- If required, compute prognostic equation for total water content
1582	IF ( humidity ) THEN
1583
1584	CALL cpu_log( log_point(29), 'q-equation', 'start' )
1585
1586	!
1587	!-- Scalar/q-tendency terms with communication
1588	sbt = tsc(2)
1589	IF ( scalar_advec == 'bc-scheme' ) THEN
1590
1591	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1592	!
1593	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1594	sbt = 1.0_wp
1595	ENDIF
1596	tend = 0.0_wp
1597	CALL advec_s_bc( q, 'q' )
1598
1599	ENDIF
1600
1601	!
1602	!-- Scalar/q-tendency terms with no communication
1603	IF ( scalar_advec /= 'bc-scheme' ) THEN
1604	tend = 0.0_wp
1605	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1606	IF ( ws_scheme_sca ) THEN
1607	CALL advec_s_ws( advc_flags_s, q, 'q', &
1608	bc_dirichlet_l .OR. bc_radiation_l, &
1609	bc_dirichlet_n .OR. bc_radiation_n, &
1610	bc_dirichlet_r .OR. bc_radiation_r, &
1611	bc_dirichlet_s .OR. bc_radiation_s )
1612	ELSE
1613	CALL advec_s_pw( q )
1614	ENDIF
1615	ELSE
1616	CALL advec_s_up( q )
1617	ENDIF
1618	ENDIF
1619
1620	CALL diffusion_s( q, &
1621	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
1622	surf_def_h(2)%qsws, &
1623	surf_lsm_h%qsws, surf_usm_h%qsws, &
1624	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
1625	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
1626	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
1627	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
1628	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
1629	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1630
1631	!
1632	!-- Sink or source of humidity due to canopy elements
1633	IF ( plant_canopy ) CALL pcm_tendency( 5 )
1634
1635	!
1636	!-- Large scale advection
1637	IF ( large_scale_forcing ) THEN
1638	CALL ls_advec( simulated_time, 'q' )
1639	ENDIF
1640
1641	!
1642	!-- Nudging
1643	IF ( nudging ) CALL nudge( simulated_time, 'q' )
1644
1645	!
1646	!-- If required compute influence of large-scale subsidence/ascent
1647	IF ( large_scale_subsidence .AND. &
1648	.NOT. use_subsidence_tendencies ) THEN
1649	CALL subsidence( tend, q, q_init, 3 )
1650	ENDIF
1651
1652	CALL module_interface_actions( 'q-tendency' )
1653
1654	!
1655	!-- Prognostic equation for total water content
1656	DO i = nxl, nxr
1657	DO j = nys, nyn
1658	DO k = nzb+1, nzt
1659	q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1660	tsc(3) * tq_m(k,j,i) ) &
1661	- tsc(5) * rdf_sc(k) * &
1662	( q(k,j,i) - q_init(k) ) &
1663	) * MERGE( 1.0_wp, 0.0_wp, &
1664	BTEST( wall_flags_0(k,j,i), 0 ) &
1665	)
1666	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
1667	ENDDO
1668	ENDDO
1669	ENDDO
1670
1671	!
1672	!-- Calculate tendencies for the next Runge-Kutta step
1673	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1674	IF ( intermediate_timestep_count == 1 ) THEN
1675	DO i = nxl, nxr
1676	DO j = nys, nyn
1677	DO k = nzb+1, nzt
1678	tq_m(k,j,i) = tend(k,j,i)
1679	ENDDO
1680	ENDDO
1681	ENDDO
1682	ELSEIF ( intermediate_timestep_count < &
1683	intermediate_timestep_count_max ) THEN
1684	DO i = nxl, nxr
1685	DO j = nys, nyn
1686	DO k = nzb+1, nzt
1687	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1688	+ 5.3125_wp * tq_m(k,j,i)
1689	ENDDO
1690	ENDDO
1691	ENDDO
1692	ENDIF
1693	ENDIF
1694
1695	CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1696
1697	ENDIF
1698	!
1699	!-- If required, compute prognostic equation for scalar
1700	IF ( passive_scalar ) THEN
1701
1702	CALL cpu_log( log_point(66), 's-equation', 'start' )
1703
1704	!
1705	!-- Scalar/q-tendency terms with communication
1706	sbt = tsc(2)
1707	IF ( scalar_advec == 'bc-scheme' ) THEN
1708
1709	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1710	!
1711	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1712	sbt = 1.0_wp
1713	ENDIF
1714	tend = 0.0_wp
1715	CALL advec_s_bc( s, 's' )
1716
1717	ENDIF
1718
1719	!
1720	!-- Scalar/q-tendency terms with no communication
1721	IF ( scalar_advec /= 'bc-scheme' ) THEN
1722	tend = 0.0_wp
1723	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1724	IF ( ws_scheme_sca ) THEN
1725	CALL advec_s_ws( advc_flags_s, s, 's', &
1726	bc_dirichlet_l .OR. bc_radiation_l, &
1727	bc_dirichlet_n .OR. bc_radiation_n, &
1728	bc_dirichlet_r .OR. bc_radiation_r, &
1729	bc_dirichlet_s .OR. bc_radiation_s )
1730	ELSE
1731	CALL advec_s_pw( s )
1732	ENDIF
1733	ELSE
1734	CALL advec_s_up( s )
1735	ENDIF
1736	ENDIF
1737
1738	CALL diffusion_s( s, &
1739	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1740	surf_def_h(2)%ssws, &
1741	surf_lsm_h%ssws, surf_usm_h%ssws, &
1742	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1743	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1744	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1745	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1746	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1747	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1748
1749	!
1750	!-- Sink or source of humidity due to canopy elements
1751	IF ( plant_canopy ) CALL pcm_tendency( 7 )
1752
1753	!
1754	!-- Large scale advection. Not implemented for scalars so far.
1755	! IF ( large_scale_forcing ) THEN
1756	! CALL ls_advec( simulated_time, 'q' )
1757	! ENDIF
1758
1759	!
1760	!-- Nudging. Not implemented for scalars so far.
1761	! IF ( nudging ) CALL nudge( simulated_time, 'q' )
1762
1763	!
1764	!-- If required compute influence of large-scale subsidence/ascent.
1765	!-- Not implemented for scalars so far.
1766	IF ( large_scale_subsidence .AND. &
1767	.NOT. use_subsidence_tendencies .AND. &
1768	.NOT. large_scale_forcing ) THEN
1769	CALL subsidence( tend, s, s_init, 3 )
1770	ENDIF
1771
1772	CALL module_interface_actions( 's-tendency' )
1773
1774	!
1775	!-- Prognostic equation for total water content
1776	DO i = nxl, nxr
1777	DO j = nys, nyn
1778	DO k = nzb+1, nzt
1779	s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1780	tsc(3) * ts_m(k,j,i) ) &
1781	- tsc(5) * rdf_sc(k) * &
1782	( s(k,j,i) - s_init(k) ) &
1783	) &
1784	* MERGE( 1.0_wp, 0.0_wp, &
1785	BTEST( wall_flags_0(k,j,i), 0 ) &
1786	)
1787	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1788	ENDDO
1789	ENDDO
1790	ENDDO
1791
1792	!
1793	!-- Calculate tendencies for the next Runge-Kutta step
1794	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1795	IF ( intermediate_timestep_count == 1 ) THEN
1796	DO i = nxl, nxr
1797	DO j = nys, nyn
1798	DO k = nzb+1, nzt
1799	ts_m(k,j,i) = tend(k,j,i)
1800	ENDDO
1801	ENDDO
1802	ENDDO
1803	ELSEIF ( intermediate_timestep_count < &
1804	intermediate_timestep_count_max ) THEN
1805	DO i = nxl, nxr
1806	DO j = nys, nyn
1807	DO k = nzb+1, nzt
1808	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1809	+ 5.3125_wp * ts_m(k,j,i)
1810	ENDDO
1811	ENDDO
1812	ENDDO
1813	ENDIF
1814	ENDIF
1815
1816	CALL cpu_log( log_point(66), 's-equation', 'stop' )
1817
1818	ENDIF
1819	!
1820	!-- Calculate prognostic equations for all other modules
1821	CALL module_interface_prognostic_equations()
1822
1823	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'end' )
1824
1825	END SUBROUTINE prognostic_equations_vector
1826
1827
1828	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |