Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |