Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |