Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |