Home

Context Navigation

prognostic_equations.f90 @ 3864

Last change on this file since 3864 was 3864, checked in by monakurppa, 5 years ago

major changes in salsa: data input, format and performance

Time-dependent emissions enabled: lod=1 for yearly PM emissions that are normalised depending on the time, and lod=2 for preprocessed emissions (similar to the chemistry module).
Additionally, 'uniform' emissions allowed. This emission is set constant on all horisontal upward facing surfaces and it is created based on parameters surface_aerosol_flux, aerosol_flux_dpg/sigmag/mass_fracs_a/mass_fracs_b.
All emissions are now implemented as surface fluxes! No 3D sources anymore.
Update the emission information by calling salsa_emission_update if skip_time_do_salsa >= time_since_reference_point and next_aero_emission_update <= time_since_reference_point
Aerosol background concentrations read from PIDS_DYNAMIC. The vertical grid must match the one applied in the model.
Gas emissions and background concentrations can be also read in in salsa_mod if the chemistry module is not applied.
In deposition, information on the land use type can be now imported from the land use model
Use SI units in PARIN, i.e. n_lognorm given in #/m3 and dpg in metres.
Apply 100 character line limit
Change all variable names from capital to lowercase letter
Change real exponents to integer if possible. If not, precalculate the value of exponent
Rename in1a to start_subrange_1a, fn2a to end_subrange_1a etc.
Rename nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and ngast --> ngases_salsa
Rename ibc to index_bc, idu to index_du etc.
Renamed loop indices b, c and sg to ib, ic and ig
run_salsa subroutine removed
Corrected a bud in salsa_driver: falsely applied ino instead of inh
Call salsa_tendency within salsa_prognostic_equations which is called in module_interface_mod instead of prognostic_equations_mod
Removed tailing white spaces and unused variables
Change error message to start by PA instead of SA

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: