Home

Context Navigation

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

Last change on this file since 3836 was 3833, checked in by forkel, 6 years ago
removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead
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 3833 2019-03-28 15:04:04Z raasch $
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
436	USE ocean_mod, &
437	ONLY: ocean_prognostic_equations, stokes_drift_terms, stokes_force, &
438	wave_breaking, wave_breaking_term
439
440	USE plant_canopy_model_mod, &
441	ONLY: cthf, pcm_tendency
442
443	#if defined( __rrtmg )
444	USE radiation_model_mod, &
445	ONLY: radiation, radiation_tendency, &
446	skip_time_do_radiation
447	#endif
448
449	USE salsa_mod, &
450	ONLY: aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, nbins, &
451	ncc_tot, ngast, salsa_boundary_conds, salsa_diagnostics, &
452	salsa_driver, salsa_gas, salsa_gases_from_chem, salsa_tendency, &
453	skip_time_do_salsa
454
455	USE salsa_util_mod, &
456	ONLY: sums_salsa_ws_l
457
458	USE statistics, &
459	ONLY: hom
460
461	USE subsidence_mod, &
462	ONLY: subsidence
463
464	USE surface_mod, &
465	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
466	surf_usm_v
467
468	USE turbulence_closure_mod, &
469	ONLY: tcm_prognostic_equations
470
471	USE wind_turbine_model_mod, &
472	ONLY: wtm_tendencies
473
474
475	PRIVATE
476	PUBLIC prognostic_equations_cache, prognostic_equations_vector
477
478	INTERFACE prognostic_equations_cache
479	MODULE PROCEDURE prognostic_equations_cache
480	END INTERFACE prognostic_equations_cache
481
482	INTERFACE prognostic_equations_vector
483	MODULE PROCEDURE prognostic_equations_vector
484	END INTERFACE prognostic_equations_vector
485
486
487	CONTAINS
488
489
490	!------------------------------------------------------------------------------!
491	! Description:
492	! ------------
493	!> Version with one optimized loop over all equations. It is only allowed to
494	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
495	!>
496	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
497	!> so communication between CPUs is not allowed (does not make sense) within
498	!> these loops.
499	!>
500	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
501	!------------------------------------------------------------------------------!
502
503	SUBROUTINE prognostic_equations_cache
504
505
506	IMPLICIT NONE
507
508	INTEGER(iwp) :: b !< index for aerosol size bins (salsa)
509	INTEGER(iwp) :: c !< index for chemical compounds (salsa)
510	INTEGER(iwp) :: g !< index for gaseous compounds (salsa)
511	INTEGER(iwp) :: i !<
512	INTEGER(iwp) :: i_omp_start !<
513	INTEGER(iwp) :: j !<
514	INTEGER(iwp) :: k !<
515	!$ INTEGER(iwp) :: omp_get_thread_num !<
516	INTEGER(iwp) :: tn = 0 !<
517
518	LOGICAL :: loop_start !<
519	INTEGER(iwp) :: lsp
520	INTEGER(iwp) :: lsp_usr !< lsp running index for chem spcs
521
522
523	!
524	!-- Time measurement can only be performed for the whole set of equations
525	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
526
527	!
528	!-- Calculation of chemical reactions. This is done outside of main loop,
529	!-- since exchange of ghost points is required after this update of the
530	!-- concentrations of chemical species
531	IF ( air_chemistry ) THEN
532	lsp_usr = 1
533	!
534	!-- Chemical reactions and deposition
535	IF ( chem_gasphase_on ) THEN
536	!
537	!-- If required, calculate photolysis frequencies -
538	!-- UNFINISHED: Why not before the intermediate timestep loop?
539	IF ( intermediate_timestep_count == 1 ) THEN
540	CALL photolysis_control
541	ENDIF
542
543	IF ( intermediate_timestep_count == 1 .OR. &
544	call_chem_at_all_substeps ) THEN
545
546	CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
547	!$OMP PARALLEL PRIVATE (i,j)
548	!$OMP DO schedule(static,1)
549	DO i = nxl, nxr
550	DO j = nys, nyn
551	CALL chem_integrate (i,j)
552	ENDDO
553	ENDDO
554	!$OMP END PARALLEL
555	CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
556
557	IF ( deposition_dry ) THEN
558	CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
559	DO i = nxl, nxr
560	DO j = nys, nyn
561	CALL chem_depo(i,j)
562	ENDDO
563	ENDDO
564	CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
565	ENDIF
566	ENDIF
567	ENDIF
568	!
569	!-- Loop over chemical species
570	CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
571	DO lsp = 1, nspec
572	CALL exchange_horiz( chem_species(lsp)%conc, nbgp )
573	lsp_usr = 1
574	DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
575	IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) ) THEN
576
577	CALL chem_boundary_conds( chem_species(lsp)%conc_p, &
578	chem_species(lsp)%conc_pr_init )
579
580	ENDIF
581	lsp_usr = lsp_usr +1
582	ENDDO
583
584
585	ENDDO
586	CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
587
588	ENDIF
589	!
590	!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
591	!-- step. The exchange of ghost points is required after this update of the
592	!-- concentrations of aerosol number and mass
593	IF ( salsa ) THEN
594
595	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
596	IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa ) &
597	THEN
598	CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
599	!$OMP PARALLEL PRIVATE (i,j,b,c,g)
600	!$OMP DO
601	!
602	!-- Call salsa processes
603	DO i = nxl, nxr
604	DO j = nys, nyn
605	CALL salsa_diagnostics( i, j )
606	CALL salsa_driver( i, j, 3 )
607	CALL salsa_diagnostics( i, j )
608	ENDDO
609	ENDDO
610
611	CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
612
613	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
614	!
615	!-- Exchange ghost points and decycle if needed.
616	DO b = 1, nbins
617	CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
618	CALL salsa_boundary_conds( aerosol_number(b)%conc_p, &
619	aerosol_number(b)%init )
620	DO c = 1, ncc_tot
621	CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
622	CALL salsa_boundary_conds( &
623	aerosol_mass((c-1)*nbins+b)%conc_p, &
624	aerosol_mass((c-1)*nbins+b)%init )
625	ENDDO
626	ENDDO
627
628	IF ( .NOT. salsa_gases_from_chem ) THEN
629	DO g = 1, ngast
630	CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
631	CALL salsa_boundary_conds( salsa_gas(g)%conc_p, &
632	salsa_gas(g)%init )
633	ENDDO
634	ENDIF
635	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
636
637	!$OMP END PARALLEL
638	last_salsa_time = time_since_reference_point
639
640	ENDIF
641
642	ENDIF
643
644	ENDIF
645
646	!
647	!-- If required, calculate cloud microphysics
648	IF ( bulk_cloud_model .AND. .NOT. microphysics_sat_adjust .AND. &
649	( intermediate_timestep_count == 1 .OR. &
650	call_microphysics_at_all_substeps ) ) &
651	THEN
652	!$OMP PARALLEL PRIVATE (i,j)
653	!$OMP DO
654	DO i = nxlg, nxrg
655	DO j = nysg, nyng
656	CALL bcm_actions( i, j )
657	ENDDO
658	ENDDO
659	!$OMP END PARALLEL
660	ENDIF
661
662	!
663	!-- Loop over all prognostic equations
664	!-- b, c ang g added for SALSA
665	!$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
666
667	!$ tn = omp_get_thread_num()
668	loop_start = .TRUE.
669
670	!$OMP DO
671	DO i = nxl, nxr
672
673	!
674	!-- Store the first loop index. It differs for each thread and is required
675	!-- later in advec_ws
676	IF ( loop_start ) THEN
677	loop_start = .FALSE.
678	i_omp_start = i
679	ENDIF
680
681	DO j = nys, nyn
682	!
683	!-- Tendency terms for u-velocity component. Please note, in case of
684	!-- non-cyclic boundary conditions the grid point i=0 is excluded from
685	!-- the prognostic equations for the u-component.
686	IF ( i >= nxlu ) THEN
687
688	tend(:,j,i) = 0.0_wp
689	IF ( timestep_scheme(1:5) == 'runge' ) THEN
690	IF ( ws_scheme_mom ) THEN
691	CALL advec_u_ws( i, j, i_omp_start, tn )
692	ELSE
693	CALL advec_u_pw( i, j )
694	ENDIF
695	ELSE
696	CALL advec_u_up( i, j )
697	ENDIF
698	CALL diffusion_u( i, j )
699	CALL coriolis( i, j, 1 )
700	IF ( sloping_surface .AND. .NOT. neutral ) THEN
701	CALL buoyancy( i, j, pt, 1 )
702	ENDIF
703
704	!
705	!-- Drag by plant canopy
706	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
707
708	!
709	!-- External pressure gradient
710	IF ( dp_external ) THEN
711	DO k = dp_level_ind_b+1, nzt
712	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
713	ENDDO
714	ENDIF
715
716	!
717	!-- Nudging
718	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
719
720	!
721	!-- Effect of Stokes drift (in ocean mode only)
722	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 )
723
724	!
725	!-- Forces by wind turbines
726	IF ( wind_turbine ) CALL wtm_tendencies( i, j, 1 )
727
728	CALL module_interface_actions( i, j, 'u-tendency' )
729	!
730	!-- Prognostic equation for u-velocity component
731	DO k = nzb+1, nzt
732
733	u_p(k,j,i) = u(k,j,i) + ( dt_3d * &
734	( tsc(2) * tend(k,j,i) + &
735	tsc(3) * tu_m(k,j,i) ) &
736	- tsc(5) * rdf(k) &
737	* ( u(k,j,i) - u_init(k) ) &
738	) * MERGE( 1.0_wp, 0.0_wp, &
739	BTEST( wall_flags_0(k,j,i), 1 )&
740	)
741	ENDDO
742
743	!
744	!-- Add turbulence generated by wave breaking (in ocean mode only)
745	IF ( wave_breaking .AND. &
746	intermediate_timestep_count == intermediate_timestep_count_max )&
747	THEN
748	CALL wave_breaking_term( i, j, 1 )
749	ENDIF
750
751	!
752	!-- Calculate tendencies for the next Runge-Kutta step
753	IF ( timestep_scheme(1:5) == 'runge' ) THEN
754	IF ( intermediate_timestep_count == 1 ) THEN
755	DO k = nzb+1, nzt
756	tu_m(k,j,i) = tend(k,j,i)
757	ENDDO
758	ELSEIF ( intermediate_timestep_count < &
759	intermediate_timestep_count_max ) THEN
760	DO k = nzb+1, nzt
761	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
762	+ 5.3125_wp * tu_m(k,j,i)
763	ENDDO
764	ENDIF
765	ENDIF
766
767	ENDIF
768	!
769	!-- Tendency terms for v-velocity component. Please note, in case of
770	!-- non-cyclic boundary conditions the grid point j=0 is excluded from
771	!-- the prognostic equations for the v-component. !--
772	IF ( j >= nysv ) THEN
773
774	tend(:,j,i) = 0.0_wp
775	IF ( timestep_scheme(1:5) == 'runge' ) THEN
776	IF ( ws_scheme_mom ) THEN
777	CALL advec_v_ws( i, j, i_omp_start, tn )
778	ELSE
779	CALL advec_v_pw( i, j )
780	ENDIF
781	ELSE
782	CALL advec_v_up( i, j )
783	ENDIF
784	CALL diffusion_v( i, j )
785	CALL coriolis( i, j, 2 )
786
787	!
788	!-- Drag by plant canopy
789	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
790
791	!
792	!-- External pressure gradient
793	IF ( dp_external ) THEN
794	DO k = dp_level_ind_b+1, nzt
795	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
796	ENDDO
797	ENDIF
798
799	!
800	!-- Nudging
801	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
802
803	!
804	!-- Effect of Stokes drift (in ocean mode only)
805	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 )
806
807	!
808	!-- Forces by wind turbines
809	IF ( wind_turbine ) CALL wtm_tendencies( i, j, 2 )
810
811	CALL module_interface_actions( i, j, 'v-tendency' )
812	!
813	!-- Prognostic equation for v-velocity component
814	DO k = nzb+1, nzt
815	v_p(k,j,i) = v(k,j,i) + ( dt_3d * &
816	( tsc(2) * tend(k,j,i) + &
817	tsc(3) * tv_m(k,j,i) ) &
818	- tsc(5) * rdf(k) &
819	* ( v(k,j,i) - v_init(k) )&
820	) * MERGE( 1.0_wp, 0.0_wp, &
821	BTEST( wall_flags_0(k,j,i), 2 )&
822	)
823	ENDDO
824
825	!
826	!-- Add turbulence generated by wave breaking (in ocean mode only)
827	IF ( wave_breaking .AND. &
828	intermediate_timestep_count == intermediate_timestep_count_max )&
829	THEN
830	CALL wave_breaking_term( i, j, 2 )
831	ENDIF
832
833	!
834	!-- Calculate tendencies for the next Runge-Kutta step
835	IF ( timestep_scheme(1:5) == 'runge' ) THEN
836	IF ( intermediate_timestep_count == 1 ) THEN
837	DO k = nzb+1, nzt
838	tv_m(k,j,i) = tend(k,j,i)
839	ENDDO
840	ELSEIF ( intermediate_timestep_count < &
841	intermediate_timestep_count_max ) THEN
842	DO k = nzb+1, nzt
843	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
844	+ 5.3125_wp * tv_m(k,j,i)
845	ENDDO
846	ENDIF
847	ENDIF
848
849	ENDIF
850
851	!
852	!-- Tendency terms for w-velocity component
853	tend(:,j,i) = 0.0_wp
854	IF ( timestep_scheme(1:5) == 'runge' ) THEN
855	IF ( ws_scheme_mom ) THEN
856	CALL advec_w_ws( i, j, i_omp_start, tn )
857	ELSE
858	CALL advec_w_pw( i, j )
859	END IF
860	ELSE
861	CALL advec_w_up( i, j )
862	ENDIF
863	CALL diffusion_w( i, j )
864	CALL coriolis( i, j, 3 )
865
866	IF ( .NOT. neutral ) THEN
867	IF ( ocean_mode ) THEN
868	CALL buoyancy( i, j, rho_ocean, 3 )
869	ELSE
870	IF ( .NOT. humidity ) THEN
871	CALL buoyancy( i, j, pt, 3 )
872	ELSE
873	CALL buoyancy( i, j, vpt, 3 )
874	ENDIF
875	ENDIF
876	ENDIF
877
878	!
879	!-- Drag by plant canopy
880	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
881
882	!
883	!-- Effect of Stokes drift (in ocean mode only)
884	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 3 )
885
886	!
887	!-- Forces by wind turbines
888	IF ( wind_turbine ) CALL wtm_tendencies( i, j, 3 )
889
890	CALL module_interface_actions( i, j, 'w-tendency' )
891	!
892	!-- Prognostic equation for w-velocity component
893	DO k = nzb+1, nzt-1
894	w_p(k,j,i) = w(k,j,i) + ( dt_3d * &
895	( tsc(2) * tend(k,j,i) + &
896	tsc(3) * tw_m(k,j,i) ) &
897	- tsc(5) * rdf(k) * w(k,j,i) &
898	) * MERGE( 1.0_wp, 0.0_wp, &
899	BTEST( wall_flags_0(k,j,i), 3 )&
900	)
901	ENDDO
902
903	!
904	!-- Calculate tendencies for the next Runge-Kutta step
905	IF ( timestep_scheme(1:5) == 'runge' ) THEN
906	IF ( intermediate_timestep_count == 1 ) THEN
907	DO k = nzb+1, nzt-1
908	tw_m(k,j,i) = tend(k,j,i)
909	ENDDO
910	ELSEIF ( intermediate_timestep_count < &
911	intermediate_timestep_count_max ) THEN
912	DO k = nzb+1, nzt-1
913	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
914	+ 5.3125_wp * tw_m(k,j,i)
915	ENDDO
916	ENDIF
917	ENDIF
918
919	!
920	!-- If required, compute prognostic equation for potential temperature
921	IF ( .NOT. neutral ) THEN
922	!
923	!-- Tendency terms for potential temperature
924	tend(:,j,i) = 0.0_wp
925	IF ( timestep_scheme(1:5) == 'runge' ) THEN
926	IF ( ws_scheme_sca ) THEN
927	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
928	flux_l_pt, diss_l_pt, i_omp_start, tn )
929	ELSE
930	CALL advec_s_pw( i, j, pt )
931	ENDIF
932	ELSE
933	CALL advec_s_up( i, j, pt )
934	ENDIF
935	CALL diffusion_s( i, j, pt, &
936	surf_def_h(0)%shf, surf_def_h(1)%shf, &
937	surf_def_h(2)%shf, &
938	surf_lsm_h%shf, surf_usm_h%shf, &
939	surf_def_v(0)%shf, surf_def_v(1)%shf, &
940	surf_def_v(2)%shf, surf_def_v(3)%shf, &
941	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
942	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
943	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
944	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
945
946	!
947	!-- Consideration of heat sources within the plant canopy
948	IF ( plant_canopy .AND. &
949	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
950	CALL pcm_tendency( i, j, 4 )
951	ENDIF
952
953	!
954	!-- Large scale advection
955	IF ( large_scale_forcing ) THEN
956	CALL ls_advec( i, j, simulated_time, 'pt' )
957	ENDIF
958
959	!
960	!-- Nudging
961	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
962
963	!
964	!-- If required, compute effect of large-scale subsidence/ascent
965	IF ( large_scale_subsidence .AND. &
966	.NOT. use_subsidence_tendencies ) THEN
967	CALL subsidence( i, j, tend, pt, pt_init, 2 )
968	ENDIF
969
970	#if defined( __rrtmg )
971	!
972	!-- If required, add tendency due to radiative heating/cooling
973	IF ( radiation .AND. &
974	simulated_time > skip_time_do_radiation ) THEN
975	CALL radiation_tendency ( i, j, tend )
976	ENDIF
977	#endif
978
979	CALL module_interface_actions( i, j, 'pt-tendency' )
980	!
981	!-- Prognostic equation for potential temperature
982	DO k = nzb+1, nzt
983	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * &
984	( tsc(2) * tend(k,j,i) + &
985	tsc(3) * tpt_m(k,j,i) ) &
986	- tsc(5) &
987	* ( pt(k,j,i) - pt_init(k) ) &
988	* ( rdf_sc(k) + ptdf_x(i) &
989	+ ptdf_y(j) ) &
990	) &
991	* MERGE( 1.0_wp, 0.0_wp, &
992	BTEST( wall_flags_0(k,j,i), 0 )&
993	)
994	ENDDO
995
996	!
997	!-- Calculate tendencies for the next Runge-Kutta step
998	IF ( timestep_scheme(1:5) == 'runge' ) THEN
999	IF ( intermediate_timestep_count == 1 ) THEN
1000	DO k = nzb+1, nzt
1001	tpt_m(k,j,i) = tend(k,j,i)
1002	ENDDO
1003	ELSEIF ( intermediate_timestep_count < &
1004	intermediate_timestep_count_max ) THEN
1005	DO k = nzb+1, nzt
1006	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1007	5.3125_wp * tpt_m(k,j,i)
1008	ENDDO
1009	ENDIF
1010	ENDIF
1011
1012	ENDIF
1013
1014	!
1015	!-- If required, compute prognostic equation for total water content
1016	IF ( humidity ) THEN
1017
1018	!
1019	!-- Tendency-terms for total water content / scalar
1020	tend(:,j,i) = 0.0_wp
1021	IF ( timestep_scheme(1:5) == 'runge' ) &
1022	THEN
1023	IF ( ws_scheme_sca ) THEN
1024	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
1025	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1026	ELSE
1027	CALL advec_s_pw( i, j, q )
1028	ENDIF
1029	ELSE
1030	CALL advec_s_up( i, j, q )
1031	ENDIF
1032	CALL diffusion_s( i, j, q, &
1033	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
1034	surf_def_h(2)%qsws, &
1035	surf_lsm_h%qsws, surf_usm_h%qsws, &
1036	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
1037	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
1038	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
1039	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
1040	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
1041	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1042
1043	!
1044	!-- Sink or source of humidity due to canopy elements
1045	IF ( plant_canopy ) CALL pcm_tendency( i, j, 5 )
1046
1047	!
1048	!-- Large scale advection
1049	IF ( large_scale_forcing ) THEN
1050	CALL ls_advec( i, j, simulated_time, 'q' )
1051	ENDIF
1052
1053	!
1054	!-- Nudging
1055	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
1056
1057	!
1058	!-- If required compute influence of large-scale subsidence/ascent
1059	IF ( large_scale_subsidence .AND. &
1060	.NOT. use_subsidence_tendencies ) THEN
1061	CALL subsidence( i, j, tend, q, q_init, 3 )
1062	ENDIF
1063
1064	CALL module_interface_actions( i, j, 'q-tendency' )
1065
1066	!
1067	!-- Prognostic equation for total water content / scalar
1068	DO k = nzb+1, nzt
1069	q_p(k,j,i) = q(k,j,i) + ( dt_3d * &
1070	( tsc(2) * tend(k,j,i) + &
1071	tsc(3) * tq_m(k,j,i) ) &
1072	- tsc(5) * rdf_sc(k) * &
1073	( q(k,j,i) - q_init(k) ) &
1074	) &
1075	* MERGE( 1.0_wp, 0.0_wp, &
1076	BTEST( wall_flags_0(k,j,i), 0 )&
1077	)
1078	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
1079	ENDDO
1080
1081	!
1082	!-- Calculate tendencies for the next Runge-Kutta step
1083	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1084	IF ( intermediate_timestep_count == 1 ) THEN
1085	DO k = nzb+1, nzt
1086	tq_m(k,j,i) = tend(k,j,i)
1087	ENDDO
1088	ELSEIF ( intermediate_timestep_count < &
1089	intermediate_timestep_count_max ) THEN
1090	DO k = nzb+1, nzt
1091	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1092	5.3125_wp * tq_m(k,j,i)
1093	ENDDO
1094	ENDIF
1095	ENDIF
1096
1097	!
1098	!-- If required, calculate prognostic equations for cloud water content
1099	!-- and cloud drop concentration
1100	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
1101	!
1102	!-- Calculate prognostic equation for cloud water content
1103	tend(:,j,i) = 0.0_wp
1104	IF ( timestep_scheme(1:5) == 'runge' ) &
1105	THEN
1106	IF ( ws_scheme_sca ) THEN
1107	CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc, &
1108	diss_s_qc, flux_l_qc, diss_l_qc, &
1109	i_omp_start, tn )
1110	ELSE
1111	CALL advec_s_pw( i, j, qc )
1112	ENDIF
1113	ELSE
1114	CALL advec_s_up( i, j, qc )
1115	ENDIF
1116	CALL diffusion_s( i, j, qc, &
1117	surf_def_h(0)%qcsws, surf_def_h(1)%qcsws, &
1118	surf_def_h(2)%qcsws, &
1119	surf_lsm_h%qcsws, surf_usm_h%qcsws, &
1120	surf_def_v(0)%qcsws, surf_def_v(1)%qcsws, &
1121	surf_def_v(2)%qcsws, surf_def_v(3)%qcsws, &
1122	surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws, &
1123	surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws, &
1124	surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws, &
1125	surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
1126
1127	!
1128	!-- Prognostic equation for cloud water content
1129	DO k = nzb+1, nzt
1130	qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * &
1131	( tsc(2) * tend(k,j,i) + &
1132	tsc(3) * tqc_m(k,j,i) )&
1133	- tsc(5) * rdf_sc(k) &
1134	* qc(k,j,i) &
1135	) &
1136	* MERGE( 1.0_wp, 0.0_wp, &
1137	BTEST( wall_flags_0(k,j,i), 0 )&
1138	)
1139	IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp
1140	ENDDO
1141	!
1142	!-- Calculate tendencies for the next Runge-Kutta step
1143	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1144	IF ( intermediate_timestep_count == 1 ) THEN
1145	DO k = nzb+1, nzt
1146	tqc_m(k,j,i) = tend(k,j,i)
1147	ENDDO
1148	ELSEIF ( intermediate_timestep_count < &
1149	intermediate_timestep_count_max ) THEN
1150	DO k = nzb+1, nzt
1151	tqc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1152	5.3125_wp * tqc_m(k,j,i)
1153	ENDDO
1154	ENDIF
1155	ENDIF
1156
1157	!
1158	!-- Calculate prognostic equation for cloud drop concentration.
1159	tend(:,j,i) = 0.0_wp
1160	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1161	IF ( ws_scheme_sca ) THEN
1162	CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc, &
1163	diss_s_nc, flux_l_nc, diss_l_nc, &
1164	i_omp_start, tn )
1165	ELSE
1166	CALL advec_s_pw( i, j, nc )
1167	ENDIF
1168	ELSE
1169	CALL advec_s_up( i, j, nc )
1170	ENDIF
1171	CALL diffusion_s( i, j, nc, &
1172	surf_def_h(0)%ncsws, surf_def_h(1)%ncsws, &
1173	surf_def_h(2)%ncsws, &
1174	surf_lsm_h%ncsws, surf_usm_h%ncsws, &
1175	surf_def_v(0)%ncsws, surf_def_v(1)%ncsws, &
1176	surf_def_v(2)%ncsws, surf_def_v(3)%ncsws, &
1177	surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws, &
1178	surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws, &
1179	surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws, &
1180	surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
1181
1182	!
1183	!-- Prognostic equation for cloud drop concentration
1184	DO k = nzb+1, nzt
1185	nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * &
1186	( tsc(2) * tend(k,j,i) + &
1187	tsc(3) * tnc_m(k,j,i) )&
1188	- tsc(5) * rdf_sc(k) &
1189	* nc(k,j,i) &
1190	) &
1191	* MERGE( 1.0_wp, 0.0_wp, &
1192	BTEST( wall_flags_0(k,j,i), 0 )&
1193	)
1194	IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp
1195	ENDDO
1196	!
1197	!-- Calculate tendencies for the next Runge-Kutta step
1198	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1199	IF ( intermediate_timestep_count == 1 ) THEN
1200	DO k = nzb+1, nzt
1201	tnc_m(k,j,i) = tend(k,j,i)
1202	ENDDO
1203	ELSEIF ( intermediate_timestep_count < &
1204	intermediate_timestep_count_max ) THEN
1205	DO k = nzb+1, nzt
1206	tnc_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1207	5.3125_wp * tnc_m(k,j,i)
1208	ENDDO
1209	ENDIF
1210	ENDIF
1211
1212	ENDIF
1213	!
1214	!-- If required, calculate prognostic equations for rain water content
1215	!-- and rain drop concentration
1216	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
1217	!
1218	!-- Calculate prognostic equation for rain water content
1219	tend(:,j,i) = 0.0_wp
1220	IF ( timestep_scheme(1:5) == 'runge' ) &
1221	THEN
1222	IF ( ws_scheme_sca ) THEN
1223	CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr, &
1224	diss_s_qr, flux_l_qr, diss_l_qr, &
1225	i_omp_start, tn )
1226	ELSE
1227	CALL advec_s_pw( i, j, qr )
1228	ENDIF
1229	ELSE
1230	CALL advec_s_up( i, j, qr )
1231	ENDIF
1232	CALL diffusion_s( i, j, qr, &
1233	surf_def_h(0)%qrsws, surf_def_h(1)%qrsws, &
1234	surf_def_h(2)%qrsws, &
1235	surf_lsm_h%qrsws, surf_usm_h%qrsws, &
1236	surf_def_v(0)%qrsws, surf_def_v(1)%qrsws, &
1237	surf_def_v(2)%qrsws, surf_def_v(3)%qrsws, &
1238	surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws, &
1239	surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws, &
1240	surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws, &
1241	surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1242
1243	!
1244	!-- Prognostic equation for rain water content
1245	DO k = nzb+1, nzt
1246	qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * &
1247	( tsc(2) * tend(k,j,i) + &
1248	tsc(3) * tqr_m(k,j,i) )&
1249	- tsc(5) * rdf_sc(k) &
1250	* qr(k,j,i) &
1251	) &
1252	* MERGE( 1.0_wp, 0.0_wp, &
1253	BTEST( wall_flags_0(k,j,i), 0 )&
1254	)
1255	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
1256	ENDDO
1257	!
1258	!-- Calculate tendencies for the next Runge-Kutta step
1259	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1260	IF ( intermediate_timestep_count == 1 ) THEN
1261	DO k = nzb+1, nzt
1262	tqr_m(k,j,i) = tend(k,j,i)
1263	ENDDO
1264	ELSEIF ( intermediate_timestep_count < &
1265	intermediate_timestep_count_max ) THEN
1266	DO k = nzb+1, nzt
1267	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1268	5.3125_wp * tqr_m(k,j,i)
1269	ENDDO
1270	ENDIF
1271	ENDIF
1272
1273	!
1274	!-- Calculate prognostic equation for rain drop concentration.
1275	tend(:,j,i) = 0.0_wp
1276	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1277	IF ( ws_scheme_sca ) THEN
1278	CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr, &
1279	diss_s_nr, flux_l_nr, diss_l_nr, &
1280	i_omp_start, tn )
1281	ELSE
1282	CALL advec_s_pw( i, j, nr )
1283	ENDIF
1284	ELSE
1285	CALL advec_s_up( i, j, nr )
1286	ENDIF
1287	CALL diffusion_s( i, j, nr, &
1288	surf_def_h(0)%nrsws, surf_def_h(1)%nrsws, &
1289	surf_def_h(2)%nrsws, &
1290	surf_lsm_h%nrsws, surf_usm_h%nrsws, &
1291	surf_def_v(0)%nrsws, surf_def_v(1)%nrsws, &
1292	surf_def_v(2)%nrsws, surf_def_v(3)%nrsws, &
1293	surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws, &
1294	surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws, &
1295	surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws, &
1296	surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1297
1298	!
1299	!-- Prognostic equation for rain drop concentration
1300	DO k = nzb+1, nzt
1301	nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * &
1302	( tsc(2) * tend(k,j,i) + &
1303	tsc(3) * tnr_m(k,j,i) )&
1304	- tsc(5) * rdf_sc(k) &
1305	* nr(k,j,i) &
1306	) &
1307	* MERGE( 1.0_wp, 0.0_wp, &
1308	BTEST( wall_flags_0(k,j,i), 0 )&
1309	)
1310	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
1311	ENDDO
1312	!
1313	!-- Calculate tendencies for the next Runge-Kutta step
1314	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1315	IF ( intermediate_timestep_count == 1 ) THEN
1316	DO k = nzb+1, nzt
1317	tnr_m(k,j,i) = tend(k,j,i)
1318	ENDDO
1319	ELSEIF ( intermediate_timestep_count < &
1320	intermediate_timestep_count_max ) THEN
1321	DO k = nzb+1, nzt
1322	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1323	5.3125_wp * tnr_m(k,j,i)
1324	ENDDO
1325	ENDIF
1326	ENDIF
1327
1328	ENDIF
1329
1330	ENDIF
1331
1332	!
1333	!-- If required, compute prognostic equation for scalar
1334	IF ( passive_scalar ) THEN
1335	!
1336	!-- Tendency-terms for total water content / scalar
1337	tend(:,j,i) = 0.0_wp
1338	IF ( timestep_scheme(1:5) == 'runge' ) &
1339	THEN
1340	IF ( ws_scheme_sca ) THEN
1341	CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1342	diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1343	ELSE
1344	CALL advec_s_pw( i, j, s )
1345	ENDIF
1346	ELSE
1347	CALL advec_s_up( i, j, s )
1348	ENDIF
1349	CALL diffusion_s( i, j, s, &
1350	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1351	surf_def_h(2)%ssws, &
1352	surf_lsm_h%ssws, surf_usm_h%ssws, &
1353	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1354	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1355	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1356	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1357	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1358	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1359
1360	!
1361	!-- Sink or source of scalar concentration due to canopy elements
1362	IF ( plant_canopy ) CALL pcm_tendency( i, j, 7 )
1363
1364	!
1365	!-- Large scale advection, still need to be extended for scalars
1366	! IF ( large_scale_forcing ) THEN
1367	! CALL ls_advec( i, j, simulated_time, 's' )
1368	! ENDIF
1369
1370	!
1371	!-- Nudging, still need to be extended for scalars
1372	! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' )
1373
1374	!
1375	!-- If required compute influence of large-scale subsidence/ascent.
1376	!-- Note, the last argument is of no meaning in this case, as it is
1377	!-- only used in conjunction with large_scale_forcing, which is to
1378	!-- date not implemented for scalars.
1379	IF ( large_scale_subsidence .AND. &
1380	.NOT. use_subsidence_tendencies .AND. &
1381	.NOT. large_scale_forcing ) THEN
1382	CALL subsidence( i, j, tend, s, s_init, 3 )
1383	ENDIF
1384
1385	CALL module_interface_actions( i, j, 's-tendency' )
1386
1387	!
1388	!-- Prognostic equation for scalar
1389	DO k = nzb+1, nzt
1390	s_p(k,j,i) = s(k,j,i) + ( dt_3d * &
1391	( tsc(2) * tend(k,j,i) + &
1392	tsc(3) * ts_m(k,j,i) ) &
1393	- tsc(5) * rdf_sc(k) &
1394	* ( s(k,j,i) - s_init(k) )&
1395	) &
1396	* MERGE( 1.0_wp, 0.0_wp, &
1397	BTEST( wall_flags_0(k,j,i), 0 )&
1398	)
1399	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1400	ENDDO
1401
1402	!
1403	!-- Calculate tendencies for the next Runge-Kutta step
1404	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1405	IF ( intermediate_timestep_count == 1 ) THEN
1406	DO k = nzb+1, nzt
1407	ts_m(k,j,i) = tend(k,j,i)
1408	ENDDO
1409	ELSEIF ( intermediate_timestep_count < &
1410	intermediate_timestep_count_max ) THEN
1411	DO k = nzb+1, nzt
1412	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1413	5.3125_wp * ts_m(k,j,i)
1414	ENDDO
1415	ENDIF
1416	ENDIF
1417
1418	ENDIF
1419	!
1420	!-- Calculate prognostic equations for turbulence closure
1421	CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1422
1423	!
1424	!-- Calculate prognostic equation for chemical quantites
1425	IF ( air_chemistry ) THEN
1426	!> TODO: remove time measurement since it slows down performance because it will be called extremely often
1427	!
1428	!-- Loop over chemical species
1429	DO lsp = 1, nvar
1430	CALL chem_prognostic_equations( chem_species(lsp)%conc_p, &
1431	chem_species(lsp)%conc, &
1432	chem_species(lsp)%tconc_m, &
1433	chem_species(lsp)%conc_pr_init, &
1434	i, j, i_omp_start, tn, lsp, &
1435	chem_species(lsp)%flux_s_cs, &
1436	chem_species(lsp)%diss_s_cs, &
1437	chem_species(lsp)%flux_l_cs, &
1438	chem_species(lsp)%diss_l_cs )
1439	ENDDO
1440
1441	ENDIF ! Chemical equations
1442	!
1443	!-- Calculate prognostic equations for the ocean
1444	IF ( ocean_mode ) THEN
1445	CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
1446	ENDIF
1447
1448	IF ( salsa ) THEN
1449	!
1450	!-- Loop over aerosol size bins: number and mass bins
1451	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
1452
1453	DO b = 1, nbins
1454	sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
1455	CALL salsa_tendency( 'aerosol_number', &
1456	aerosol_number(b)%conc_p, &
1457	aerosol_number(b)%conc, &
1458	aerosol_number(b)%tconc_m, &
1459	i, j, i_omp_start, tn, b, b, &
1460	aerosol_number(b)%flux_s, &
1461	aerosol_number(b)%diss_s, &
1462	aerosol_number(b)%flux_l, &
1463	aerosol_number(b)%diss_l, &
1464	aerosol_number(b)%init )
1465	aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
1466	DO c = 1, ncc_tot
1467	sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
1468	CALL salsa_tendency( 'aerosol_mass', &
1469	aerosol_mass((c-1)*nbins+b)%conc_p,&
1470	aerosol_mass((c-1)*nbins+b)%conc, &
1471	aerosol_mass((c-1)*nbins+b)%tconc_m,&
1472	i, j, i_omp_start, tn, b, c, &
1473	aerosol_mass((c-1)*nbins+b)%flux_s,&
1474	aerosol_mass((c-1)*nbins+b)%diss_s,&
1475	aerosol_mass((c-1)*nbins+b)%flux_l,&
1476	aerosol_mass((c-1)*nbins+b)%diss_l,&
1477	aerosol_mass((c-1)*nbins+b)%init )
1478	aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
1479	ENDDO
1480	ENDDO
1481	IF ( .NOT. salsa_gases_from_chem ) THEN
1482	DO g = 1, ngast
1483	sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
1484	CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p, &
1485	salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
1486	i, j, i_omp_start, tn, g, g, &
1487	salsa_gas(g)%flux_s, salsa_gas(g)%diss_s,&
1488	salsa_gas(g)%flux_l, salsa_gas(g)%diss_l,&
1489	salsa_gas(g)%init )
1490	salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
1491	ENDDO
1492	ENDIF
1493
1494	ENDIF
1495
1496	ENDIF
1497
1498	ENDDO ! loop over j
1499	ENDDO ! loop over i
1500	!$OMP END PARALLEL
1501
1502
1503	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1504
1505
1506	END SUBROUTINE prognostic_equations_cache
1507
1508
1509	!------------------------------------------------------------------------------!
1510	! Description:
1511	! ------------
1512	!> Version for vector machines
1513	!------------------------------------------------------------------------------!
1514
1515	SUBROUTINE prognostic_equations_vector
1516
1517
1518	IMPLICIT NONE
1519
1520	INTEGER(iwp) :: b !< index for aerosol size bins (salsa)
1521	INTEGER(iwp) :: c !< index for chemical compounds (salsa)
1522	INTEGER(iwp) :: g !< index for gaseous compounds (salsa)
1523	INTEGER(iwp) :: i !<
1524	INTEGER(iwp) :: j !<
1525	INTEGER(iwp) :: k !<
1526	INTEGER(iwp) :: lsp !< running index for chemical species
1527
1528	REAL(wp) :: sbt !<
1529
1530	!
1531	!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1532	!-- step. The exchange of ghost points is required after this update of the
1533	!-- concentrations of aerosol number and mass
1534	IF ( salsa ) THEN
1535
1536	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
1537
1538	IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa ) &
1539	THEN
1540
1541	CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1542	!$OMP PARALLEL PRIVATE (i,j,b,c,g)
1543	!$OMP DO
1544	!
1545	!-- Call salsa processes
1546	DO i = nxl, nxr
1547	DO j = nys, nyn
1548	CALL salsa_diagnostics( i, j )
1549	CALL salsa_driver( i, j, 3 )
1550	CALL salsa_diagnostics( i, j )
1551	ENDDO
1552	ENDDO
1553
1554	CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1555
1556	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1557	!
1558	!-- Exchange ghost points and decycle if needed.
1559	DO b = 1, nbins
1560	CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
1561	CALL salsa_boundary_conds( aerosol_number(b)%conc_p, &
1562	aerosol_number(b)%init )
1563	DO c = 1, ncc_tot
1564	CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
1565	CALL salsa_boundary_conds( &
1566	aerosol_mass((c-1)*nbins+b)%conc_p, &
1567	aerosol_mass((c-1)*nbins+b)%init )
1568	ENDDO
1569	ENDDO
1570
1571	IF ( .NOT. salsa_gases_from_chem ) THEN
1572	DO g = 1, ngast
1573	CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
1574	CALL salsa_boundary_conds( salsa_gas(g)%conc_p, &
1575	salsa_gas(g)%init )
1576	ENDDO
1577	ENDIF
1578	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1579
1580	!$OMP END PARALLEL
1581	last_salsa_time = time_since_reference_point
1582
1583	ENDIF
1584
1585	ENDIF
1586
1587	ENDIF
1588
1589	!
1590	!-- If required, calculate cloud microphysical impacts
1591	IF ( bulk_cloud_model .AND. .NOT. microphysics_sat_adjust .AND. &
1592	( intermediate_timestep_count == 1 .OR. &
1593	call_microphysics_at_all_substeps ) &
1594	) THEN
1595	CALL cpu_log( log_point(51), 'microphysics', 'start' )
1596	CALL bcm_actions
1597	CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1598	ENDIF
1599
1600	!
1601	!-- u-velocity component
1602	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1603
1604	!$ACC KERNELS PRESENT(tend)
1605	tend = 0.0_wp
1606	!$ACC END KERNELS
1607	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1608	IF ( ws_scheme_mom ) THEN
1609	CALL advec_u_ws
1610	ELSE
1611	CALL advec_u_pw
1612	ENDIF
1613	ELSE
1614	CALL advec_u_up
1615	ENDIF
1616	CALL diffusion_u
1617	CALL coriolis( 1 )
1618	IF ( sloping_surface .AND. .NOT. neutral ) THEN
1619	CALL buoyancy( pt, 1 )
1620	ENDIF
1621
1622	!
1623	!-- Drag by plant canopy
1624	IF ( plant_canopy ) CALL pcm_tendency( 1 )
1625
1626	!
1627	!-- External pressure gradient
1628	IF ( dp_external ) THEN
1629	DO i = nxlu, nxr
1630	DO j = nys, nyn
1631	DO k = dp_level_ind_b+1, nzt
1632	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1633	ENDDO
1634	ENDDO
1635	ENDDO
1636	ENDIF
1637
1638	!
1639	!-- Nudging
1640	IF ( nudging ) CALL nudge( simulated_time, 'u' )
1641
1642	!
1643	!-- Effect of Stokes drift (in ocean mode only)
1644	IF ( stokes_force ) CALL stokes_drift_terms( 1 )
1645
1646	!
1647	!-- Forces by wind turbines
1648	IF ( wind_turbine ) CALL wtm_tendencies( 1 )
1649
1650	CALL module_interface_actions( 'u-tendency' )
1651
1652	!
1653	!-- Prognostic equation for u-velocity component
1654	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1655	!$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1656	!$ACC PRESENT(tsc(2:5)) &
1657	!$ACC PRESENT(u_p)
1658	DO i = nxlu, nxr
1659	DO j = nys, nyn
1660	DO k = nzb+1, nzt
1661	u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1662	tsc(3) * tu_m(k,j,i) ) &
1663	- tsc(5) * rdf(k) * &
1664	( u(k,j,i) - u_init(k) ) &
1665	) * MERGE( 1.0_wp, 0.0_wp, &
1666	BTEST( wall_flags_0(k,j,i), 1 ) &
1667	)
1668	ENDDO
1669	ENDDO
1670	ENDDO
1671
1672	!
1673	!-- Add turbulence generated by wave breaking (in ocean mode only)
1674	IF ( wave_breaking .AND. &
1675	intermediate_timestep_count == intermediate_timestep_count_max ) &
1676	THEN
1677	CALL wave_breaking_term( 1 )
1678	ENDIF
1679
1680	!
1681	!-- Calculate tendencies for the next Runge-Kutta step
1682	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1683	IF ( intermediate_timestep_count == 1 ) THEN
1684	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1685	!$ACC PRESENT(tend, tu_m)
1686	DO i = nxlu, nxr
1687	DO j = nys, nyn
1688	DO k = nzb+1, nzt
1689	tu_m(k,j,i) = tend(k,j,i)
1690	ENDDO
1691	ENDDO
1692	ENDDO
1693	ELSEIF ( intermediate_timestep_count < &
1694	intermediate_timestep_count_max ) THEN
1695	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1696	!$ACC PRESENT(tend, tu_m)
1697	DO i = nxlu, nxr
1698	DO j = nys, nyn
1699	DO k = nzb+1, nzt
1700	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1701	+ 5.3125_wp * tu_m(k,j,i)
1702	ENDDO
1703	ENDDO
1704	ENDDO
1705	ENDIF
1706	ENDIF
1707
1708	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1709
1710	!
1711	!-- v-velocity component
1712	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1713
1714	!$ACC KERNELS PRESENT(tend)
1715	tend = 0.0_wp
1716	!$ACC END KERNELS
1717	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1718	IF ( ws_scheme_mom ) THEN
1719	CALL advec_v_ws
1720	ELSE
1721	CALL advec_v_pw
1722	END IF
1723	ELSE
1724	CALL advec_v_up
1725	ENDIF
1726	CALL diffusion_v
1727	CALL coriolis( 2 )
1728
1729	!
1730	!-- Drag by plant canopy
1731	IF ( plant_canopy ) CALL pcm_tendency( 2 )
1732
1733	!
1734	!-- External pressure gradient
1735	IF ( dp_external ) THEN
1736	DO i = nxl, nxr
1737	DO j = nysv, nyn
1738	DO k = dp_level_ind_b+1, nzt
1739	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1740	ENDDO
1741	ENDDO
1742	ENDDO
1743	ENDIF
1744
1745	!
1746	!-- Nudging
1747	IF ( nudging ) CALL nudge( simulated_time, 'v' )
1748
1749	!
1750	!-- Effect of Stokes drift (in ocean mode only)
1751	IF ( stokes_force ) CALL stokes_drift_terms( 2 )
1752
1753	!
1754	!-- Forces by wind turbines
1755	IF ( wind_turbine ) CALL wtm_tendencies( 2 )
1756
1757	CALL module_interface_actions( 'v-tendency' )
1758
1759	!
1760	!-- Prognostic equation for v-velocity component
1761	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1762	!$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1763	!$ACC PRESENT(tsc(2:5)) &
1764	!$ACC PRESENT(v_p)
1765	DO i = nxl, nxr
1766	DO j = nysv, nyn
1767	DO k = nzb+1, nzt
1768	v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1769	tsc(3) * tv_m(k,j,i) ) &
1770	- tsc(5) * rdf(k) * &
1771	( v(k,j,i) - v_init(k) ) &
1772	) * MERGE( 1.0_wp, 0.0_wp, &
1773	BTEST( wall_flags_0(k,j,i), 2 )&
1774	)
1775	ENDDO
1776	ENDDO
1777	ENDDO
1778
1779	!
1780	!-- Add turbulence generated by wave breaking (in ocean mode only)
1781	IF ( wave_breaking .AND. &
1782	intermediate_timestep_count == intermediate_timestep_count_max ) &
1783	THEN
1784	CALL wave_breaking_term( 2 )
1785	ENDIF
1786
1787	!
1788	!-- Calculate tendencies for the next Runge-Kutta step
1789	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1790	IF ( intermediate_timestep_count == 1 ) THEN
1791	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1792	!$ACC PRESENT(tend, tv_m)
1793	DO i = nxl, nxr
1794	DO j = nysv, nyn
1795	DO k = nzb+1, nzt
1796	tv_m(k,j,i) = tend(k,j,i)
1797	ENDDO
1798	ENDDO
1799	ENDDO
1800	ELSEIF ( intermediate_timestep_count < &
1801	intermediate_timestep_count_max ) THEN
1802	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1803	!$ACC PRESENT(tend, tv_m)
1804	DO i = nxl, nxr
1805	DO j = nysv, nyn
1806	DO k = nzb+1, nzt
1807	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1808	+ 5.3125_wp * tv_m(k,j,i)
1809	ENDDO
1810	ENDDO
1811	ENDDO
1812	ENDIF
1813	ENDIF
1814
1815	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1816
1817	!
1818	!-- w-velocity component
1819	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1820
1821	!$ACC KERNELS PRESENT(tend)
1822	tend = 0.0_wp
1823	!$ACC END KERNELS
1824	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1825	IF ( ws_scheme_mom ) THEN
1826	CALL advec_w_ws
1827	ELSE
1828	CALL advec_w_pw
1829	ENDIF
1830	ELSE
1831	CALL advec_w_up
1832	ENDIF
1833	CALL diffusion_w
1834	CALL coriolis( 3 )
1835
1836	IF ( .NOT. neutral ) THEN
1837	IF ( ocean_mode ) THEN
1838	CALL buoyancy( rho_ocean, 3 )
1839	ELSE
1840	IF ( .NOT. humidity ) THEN
1841	CALL buoyancy( pt, 3 )
1842	ELSE
1843	CALL buoyancy( vpt, 3 )
1844	ENDIF
1845	ENDIF
1846	ENDIF
1847
1848	!
1849	!-- Drag by plant canopy
1850	IF ( plant_canopy ) CALL pcm_tendency( 3 )
1851
1852	!
1853	!-- Effect of Stokes drift (in ocean mode only)
1854	IF ( stokes_force ) CALL stokes_drift_terms( 3 )
1855
1856	!
1857	!-- Forces by wind turbines
1858	IF ( wind_turbine ) CALL wtm_tendencies( 3 )
1859
1860	CALL module_interface_actions( 'w-tendency' )
1861
1862	!
1863	!-- Prognostic equation for w-velocity component
1864	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1865	!$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1866	!$ACC PRESENT(tsc(2:5)) &
1867	!$ACC PRESENT(w_p)
1868	DO i = nxl, nxr
1869	DO j = nys, nyn
1870	DO k = nzb+1, nzt-1
1871	w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1872	tsc(3) * tw_m(k,j,i) ) &
1873	- tsc(5) * rdf(k) * w(k,j,i) &
1874	) * MERGE( 1.0_wp, 0.0_wp, &
1875	BTEST( wall_flags_0(k,j,i), 3 )&
1876	)
1877	ENDDO
1878	ENDDO
1879	ENDDO
1880
1881	!
1882	!-- Calculate tendencies for the next Runge-Kutta step
1883	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1884	IF ( intermediate_timestep_count == 1 ) THEN
1885	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1886	!$ACC PRESENT(tend, tw_m)
1887	DO i = nxl, nxr
1888	DO j = nys, nyn
1889	DO k = nzb+1, nzt-1
1890	tw_m(k,j,i) = tend(k,j,i)
1891	ENDDO
1892	ENDDO
1893	ENDDO
1894	ELSEIF ( intermediate_timestep_count < &
1895	intermediate_timestep_count_max ) THEN
1896	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1897	!$ACC PRESENT(tend, tw_m)
1898	DO i = nxl, nxr
1899	DO j = nys, nyn
1900	DO k = nzb+1, nzt-1
1901	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1902	+ 5.3125_wp * tw_m(k,j,i)
1903	ENDDO
1904	ENDDO
1905	ENDDO
1906	ENDIF
1907	ENDIF
1908
1909	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1910
1911
1912	!
1913	!-- If required, compute prognostic equation for potential temperature
1914	IF ( .NOT. neutral ) THEN
1915
1916	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1917
1918	!
1919	!-- pt-tendency terms with communication
1920	sbt = tsc(2)
1921	IF ( scalar_advec == 'bc-scheme' ) THEN
1922
1923	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1924	!
1925	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1926	sbt = 1.0_wp
1927	ENDIF
1928	tend = 0.0_wp
1929	CALL advec_s_bc( pt, 'pt' )
1930
1931	ENDIF
1932
1933	!
1934	!-- pt-tendency terms with no communication
1935	IF ( scalar_advec /= 'bc-scheme' ) THEN
1936	!$ACC KERNELS PRESENT(tend)
1937	tend = 0.0_wp
1938	!$ACC END KERNELS
1939	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1940	IF ( ws_scheme_sca ) THEN
1941	CALL advec_s_ws( pt, 'pt' )
1942	ELSE
1943	CALL advec_s_pw( pt )
1944	ENDIF
1945	ELSE
1946	CALL advec_s_up( pt )
1947	ENDIF
1948	ENDIF
1949
1950	CALL diffusion_s( pt, &
1951	surf_def_h(0)%shf, surf_def_h(1)%shf, &
1952	surf_def_h(2)%shf, &
1953	surf_lsm_h%shf, surf_usm_h%shf, &
1954	surf_def_v(0)%shf, surf_def_v(1)%shf, &
1955	surf_def_v(2)%shf, surf_def_v(3)%shf, &
1956	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
1957	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
1958	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
1959	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1960
1961	!
1962	!-- Consideration of heat sources within the plant canopy
1963	IF ( plant_canopy .AND. &
1964	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
1965	CALL pcm_tendency( 4 )
1966	ENDIF
1967
1968	!
1969	!-- Large scale advection
1970	IF ( large_scale_forcing ) THEN
1971	CALL ls_advec( simulated_time, 'pt' )
1972	ENDIF
1973
1974	!
1975	!-- Nudging
1976	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
1977
1978	!
1979	!-- If required compute influence of large-scale subsidence/ascent
1980	IF ( large_scale_subsidence .AND. &
1981	.NOT. use_subsidence_tendencies ) THEN
1982	CALL subsidence( tend, pt, pt_init, 2 )
1983	ENDIF
1984
1985	#if defined( __rrtmg )
1986	!
1987	!-- If required, add tendency due to radiative heating/cooling
1988	IF ( radiation .AND. &
1989	simulated_time > skip_time_do_radiation ) THEN
1990	CALL radiation_tendency ( tend )
1991	ENDIF
1992	#endif
1993
1994	CALL module_interface_actions( 'pt-tendency' )
1995
1996	!
1997	!-- Prognostic equation for potential temperature
1998	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1999	!$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
2000	!$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
2001	!$ACC PRESENT(tsc(3:5)) &
2002	!$ACC PRESENT(pt_p)
2003	DO i = nxl, nxr
2004	DO j = nys, nyn
2005	DO k = nzb+1, nzt
2006	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2007	tsc(3) * tpt_m(k,j,i) ) &
2008	- tsc(5) * &
2009	( pt(k,j,i) - pt_init(k) ) *&
2010	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
2011	) &
2012	* MERGE( 1.0_wp, 0.0_wp, &
2013	BTEST( wall_flags_0(k,j,i), 0 ) &
2014	)
2015	ENDDO
2016	ENDDO
2017	ENDDO
2018	!
2019	!-- Calculate tendencies for the next Runge-Kutta step
2020	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2021	IF ( intermediate_timestep_count == 1 ) THEN
2022	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2023	!$ACC PRESENT(tend, tpt_m)
2024	DO i = nxl, nxr
2025	DO j = nys, nyn
2026	DO k = nzb+1, nzt
2027	tpt_m(k,j,i) = tend(k,j,i)
2028	ENDDO
2029	ENDDO
2030	ENDDO
2031	ELSEIF ( intermediate_timestep_count < &
2032	intermediate_timestep_count_max ) THEN
2033	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2034	!$ACC PRESENT(tend, tpt_m)
2035	DO i = nxl, nxr
2036	DO j = nys, nyn
2037	DO k = nzb+1, nzt
2038	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2039	5.3125_wp * tpt_m(k,j,i)
2040	ENDDO
2041	ENDDO
2042	ENDDO
2043	ENDIF
2044	ENDIF
2045
2046	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
2047
2048	ENDIF
2049
2050	!
2051	!-- If required, compute prognostic equation for total water content
2052	IF ( humidity ) THEN
2053
2054	CALL cpu_log( log_point(29), 'q-equation', 'start' )
2055
2056	!
2057	!-- Scalar/q-tendency terms with communication
2058	sbt = tsc(2)
2059	IF ( scalar_advec == 'bc-scheme' ) THEN
2060
2061	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2062	!
2063	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2064	sbt = 1.0_wp
2065	ENDIF
2066	tend = 0.0_wp
2067	CALL advec_s_bc( q, 'q' )
2068
2069	ENDIF
2070
2071	!
2072	!-- Scalar/q-tendency terms with no communication
2073	IF ( scalar_advec /= 'bc-scheme' ) THEN
2074	tend = 0.0_wp
2075	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2076	IF ( ws_scheme_sca ) THEN
2077	CALL advec_s_ws( q, 'q' )
2078	ELSE
2079	CALL advec_s_pw( q )
2080	ENDIF
2081	ELSE
2082	CALL advec_s_up( q )
2083	ENDIF
2084	ENDIF
2085
2086	CALL diffusion_s( q, &
2087	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
2088	surf_def_h(2)%qsws, &
2089	surf_lsm_h%qsws, surf_usm_h%qsws, &
2090	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
2091	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
2092	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
2093	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
2094	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
2095	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
2096
2097	!
2098	!-- Sink or source of humidity due to canopy elements
2099	IF ( plant_canopy ) CALL pcm_tendency( 5 )
2100
2101	!
2102	!-- Large scale advection
2103	IF ( large_scale_forcing ) THEN
2104	CALL ls_advec( simulated_time, 'q' )
2105	ENDIF
2106
2107	!
2108	!-- Nudging
2109	IF ( nudging ) CALL nudge( simulated_time, 'q' )
2110
2111	!
2112	!-- If required compute influence of large-scale subsidence/ascent
2113	IF ( large_scale_subsidence .AND. &
2114	.NOT. use_subsidence_tendencies ) THEN
2115	CALL subsidence( tend, q, q_init, 3 )
2116	ENDIF
2117
2118	CALL module_interface_actions( 'q-tendency' )
2119
2120	!
2121	!-- Prognostic equation for total water content
2122	DO i = nxl, nxr
2123	DO j = nys, nyn
2124	DO k = nzb+1, nzt
2125	q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2126	tsc(3) * tq_m(k,j,i) ) &
2127	- tsc(5) * rdf_sc(k) * &
2128	( q(k,j,i) - q_init(k) ) &
2129	) * MERGE( 1.0_wp, 0.0_wp, &
2130	BTEST( wall_flags_0(k,j,i), 0 ) &
2131	)
2132	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
2133	ENDDO
2134	ENDDO
2135	ENDDO
2136
2137	!
2138	!-- Calculate tendencies for the next Runge-Kutta step
2139	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2140	IF ( intermediate_timestep_count == 1 ) THEN
2141	DO i = nxl, nxr
2142	DO j = nys, nyn
2143	DO k = nzb+1, nzt
2144	tq_m(k,j,i) = tend(k,j,i)
2145	ENDDO
2146	ENDDO
2147	ENDDO
2148	ELSEIF ( intermediate_timestep_count < &
2149	intermediate_timestep_count_max ) THEN
2150	DO i = nxl, nxr
2151	DO j = nys, nyn
2152	DO k = nzb+1, nzt
2153	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2154	+ 5.3125_wp * tq_m(k,j,i)
2155	ENDDO
2156	ENDDO
2157	ENDDO
2158	ENDIF
2159	ENDIF
2160
2161	CALL cpu_log( log_point(29), 'q-equation', 'stop' )
2162
2163	!
2164	!-- If required, calculate prognostic equations for cloud water content
2165	!-- and cloud drop concentration
2166	IF ( bulk_cloud_model .AND. microphysics_morrison ) THEN
2167
2168	CALL cpu_log( log_point(67), 'qc-equation', 'start' )
2169
2170	!
2171	!-- Calculate prognostic equation for cloud water content
2172	sbt = tsc(2)
2173	IF ( scalar_advec == 'bc-scheme' ) THEN
2174
2175	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2176	!
2177	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2178	sbt = 1.0_wp
2179	ENDIF
2180	tend = 0.0_wp
2181	CALL advec_s_bc( qc, 'qc' )
2182
2183	ENDIF
2184
2185	!
2186	!-- qc-tendency terms with no communication
2187	IF ( scalar_advec /= 'bc-scheme' ) THEN
2188	tend = 0.0_wp
2189	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2190	IF ( ws_scheme_sca ) THEN
2191	CALL advec_s_ws( qc, 'qc' )
2192	ELSE
2193	CALL advec_s_pw( qc )
2194	ENDIF
2195	ELSE
2196	CALL advec_s_up( qc )
2197	ENDIF
2198	ENDIF
2199
2200	CALL diffusion_s( qc, &
2201	surf_def_h(0)%qcsws, surf_def_h(1)%qcsws, &
2202	surf_def_h(2)%qcsws, &
2203	surf_lsm_h%qcsws, surf_usm_h%qcsws, &
2204	surf_def_v(0)%qcsws, surf_def_v(1)%qcsws, &
2205	surf_def_v(2)%qcsws, surf_def_v(3)%qcsws, &
2206	surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws, &
2207	surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws, &
2208	surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws, &
2209	surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
2210
2211	!
2212	!-- Prognostic equation for cloud water content
2213	DO i = nxl, nxr
2214	DO j = nys, nyn
2215	DO k = nzb+1, nzt
2216	qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2217	tsc(3) * tqc_m(k,j,i) ) &
2218	- tsc(5) * rdf_sc(k) * &
2219	qc(k,j,i) &
2220	) &
2221	* MERGE( 1.0_wp, 0.0_wp, &
2222	BTEST( wall_flags_0(k,j,i), 0 ) &
2223	)
2224	IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp
2225	ENDDO
2226	ENDDO
2227	ENDDO
2228
2229	!
2230	!-- Calculate tendencies for the next Runge-Kutta step
2231	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2232	IF ( intermediate_timestep_count == 1 ) THEN
2233	DO i = nxl, nxr
2234	DO j = nys, nyn
2235	DO k = nzb+1, nzt
2236	tqc_m(k,j,i) = tend(k,j,i)
2237	ENDDO
2238	ENDDO
2239	ENDDO
2240	ELSEIF ( intermediate_timestep_count < &
2241	intermediate_timestep_count_max ) THEN
2242	DO i = nxl, nxr
2243	DO j = nys, nyn
2244	DO k = nzb+1, nzt
2245	tqc_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2246	+ 5.3125_wp * tqc_m(k,j,i)
2247	ENDDO
2248	ENDDO
2249	ENDDO
2250	ENDIF
2251	ENDIF
2252
2253	CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
2254
2255	CALL cpu_log( log_point(68), 'nc-equation', 'start' )
2256	!
2257	!-- Calculate prognostic equation for cloud drop concentration
2258	sbt = tsc(2)
2259	IF ( scalar_advec == 'bc-scheme' ) THEN
2260
2261	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2262	!
2263	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2264	sbt = 1.0_wp
2265	ENDIF
2266	tend = 0.0_wp
2267	CALL advec_s_bc( nc, 'nc' )
2268
2269	ENDIF
2270
2271	!
2272	!-- nc-tendency terms with no communication
2273	IF ( scalar_advec /= 'bc-scheme' ) THEN
2274	tend = 0.0_wp
2275	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2276	IF ( ws_scheme_sca ) THEN
2277	CALL advec_s_ws( nc, 'nc' )
2278	ELSE
2279	CALL advec_s_pw( nc )
2280	ENDIF
2281	ELSE
2282	CALL advec_s_up( nc )
2283	ENDIF
2284	ENDIF
2285
2286	CALL diffusion_s( nc, &
2287	surf_def_h(0)%ncsws, surf_def_h(1)%ncsws, &
2288	surf_def_h(2)%ncsws, &
2289	surf_lsm_h%ncsws, surf_usm_h%ncsws, &
2290	surf_def_v(0)%ncsws, surf_def_v(1)%ncsws, &
2291	surf_def_v(2)%ncsws, surf_def_v(3)%ncsws, &
2292	surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws, &
2293	surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws, &
2294	surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws, &
2295	surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
2296
2297	!
2298	!-- Prognostic equation for cloud drop concentration
2299	DO i = nxl, nxr
2300	DO j = nys, nyn
2301	DO k = nzb+1, nzt
2302	nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2303	tsc(3) * tnc_m(k,j,i) ) &
2304	- tsc(5) * rdf_sc(k) * &
2305	nc(k,j,i) &
2306	) &
2307	* MERGE( 1.0_wp, 0.0_wp, &
2308	BTEST( wall_flags_0(k,j,i), 0 ) &
2309	)
2310	IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp
2311	ENDDO
2312	ENDDO
2313	ENDDO
2314
2315	!
2316	!-- Calculate tendencies for the next Runge-Kutta step
2317	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2318	IF ( intermediate_timestep_count == 1 ) THEN
2319	DO i = nxl, nxr
2320	DO j = nys, nyn
2321	DO k = nzb+1, nzt
2322	tnc_m(k,j,i) = tend(k,j,i)
2323	ENDDO
2324	ENDDO
2325	ENDDO
2326	ELSEIF ( intermediate_timestep_count < &
2327	intermediate_timestep_count_max ) THEN
2328	DO i = nxl, nxr
2329	DO j = nys, nyn
2330	DO k = nzb+1, nzt
2331	tnc_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2332	+ 5.3125_wp * tnc_m(k,j,i)
2333	ENDDO
2334	ENDDO
2335	ENDDO
2336	ENDIF
2337	ENDIF
2338
2339	CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
2340
2341	ENDIF
2342	!
2343	!-- If required, calculate prognostic equations for rain water content
2344	!-- and rain drop concentration
2345	IF ( bulk_cloud_model .AND. microphysics_seifert ) THEN
2346
2347	CALL cpu_log( log_point(52), 'qr-equation', 'start' )
2348
2349	!
2350	!-- Calculate prognostic equation for rain water content
2351	sbt = tsc(2)
2352	IF ( scalar_advec == 'bc-scheme' ) THEN
2353
2354	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2355	!
2356	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2357	sbt = 1.0_wp
2358	ENDIF
2359	tend = 0.0_wp
2360	CALL advec_s_bc( qr, 'qr' )
2361
2362	ENDIF
2363
2364	!
2365	!-- qr-tendency terms with no communication
2366	IF ( scalar_advec /= 'bc-scheme' ) THEN
2367	tend = 0.0_wp
2368	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2369	IF ( ws_scheme_sca ) THEN
2370	CALL advec_s_ws( qr, 'qr' )
2371	ELSE
2372	CALL advec_s_pw( qr )
2373	ENDIF
2374	ELSE
2375	CALL advec_s_up( qr )
2376	ENDIF
2377	ENDIF
2378
2379	CALL diffusion_s( qr, &
2380	surf_def_h(0)%qrsws, surf_def_h(1)%qrsws, &
2381	surf_def_h(2)%qrsws, &
2382	surf_lsm_h%qrsws, surf_usm_h%qrsws, &
2383	surf_def_v(0)%qrsws, surf_def_v(1)%qrsws, &
2384	surf_def_v(2)%qrsws, surf_def_v(3)%qrsws, &
2385	surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws, &
2386	surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws, &
2387	surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws, &
2388	surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
2389
2390	!
2391	!-- Prognostic equation for rain water content
2392	DO i = nxl, nxr
2393	DO j = nys, nyn
2394	DO k = nzb+1, nzt
2395	qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2396	tsc(3) * tqr_m(k,j,i) ) &
2397	- tsc(5) * rdf_sc(k) * &
2398	qr(k,j,i) &
2399	) &
2400	* MERGE( 1.0_wp, 0.0_wp, &
2401	BTEST( wall_flags_0(k,j,i), 0 ) &
2402	)
2403	IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp
2404	ENDDO
2405	ENDDO
2406	ENDDO
2407
2408	!
2409	!-- Calculate tendencies for the next Runge-Kutta step
2410	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2411	IF ( intermediate_timestep_count == 1 ) THEN
2412	DO i = nxl, nxr
2413	DO j = nys, nyn
2414	DO k = nzb+1, nzt
2415	tqr_m(k,j,i) = tend(k,j,i)
2416	ENDDO
2417	ENDDO
2418	ENDDO
2419	ELSEIF ( intermediate_timestep_count < &
2420	intermediate_timestep_count_max ) THEN
2421	DO i = nxl, nxr
2422	DO j = nys, nyn
2423	DO k = nzb+1, nzt
2424	tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2425	+ 5.3125_wp * tqr_m(k,j,i)
2426	ENDDO
2427	ENDDO
2428	ENDDO
2429	ENDIF
2430	ENDIF
2431
2432	CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
2433	CALL cpu_log( log_point(53), 'nr-equation', 'start' )
2434
2435	!
2436	!-- Calculate prognostic equation for rain drop concentration
2437	sbt = tsc(2)
2438	IF ( scalar_advec == 'bc-scheme' ) THEN
2439
2440	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2441	!
2442	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2443	sbt = 1.0_wp
2444	ENDIF
2445	tend = 0.0_wp
2446	CALL advec_s_bc( nr, 'nr' )
2447
2448	ENDIF
2449
2450	!
2451	!-- nr-tendency terms with no communication
2452	IF ( scalar_advec /= 'bc-scheme' ) THEN
2453	tend = 0.0_wp
2454	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2455	IF ( ws_scheme_sca ) THEN
2456	CALL advec_s_ws( nr, 'nr' )
2457	ELSE
2458	CALL advec_s_pw( nr )
2459	ENDIF
2460	ELSE
2461	CALL advec_s_up( nr )
2462	ENDIF
2463	ENDIF
2464
2465	CALL diffusion_s( nr, &
2466	surf_def_h(0)%nrsws, surf_def_h(1)%nrsws, &
2467	surf_def_h(2)%nrsws, &
2468	surf_lsm_h%nrsws, surf_usm_h%nrsws, &
2469	surf_def_v(0)%nrsws, surf_def_v(1)%nrsws, &
2470	surf_def_v(2)%nrsws, surf_def_v(3)%nrsws, &
2471	surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws, &
2472	surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws, &
2473	surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws, &
2474	surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
2475
2476	!
2477	!-- Prognostic equation for rain drop concentration
2478	DO i = nxl, nxr
2479	DO j = nys, nyn
2480	DO k = nzb+1, nzt
2481	nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2482	tsc(3) * tnr_m(k,j,i) ) &
2483	- tsc(5) * rdf_sc(k) * &
2484	nr(k,j,i) &
2485	) &
2486	* MERGE( 1.0_wp, 0.0_wp, &
2487	BTEST( wall_flags_0(k,j,i), 0 ) &
2488	)
2489	IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp
2490	ENDDO
2491	ENDDO
2492	ENDDO
2493
2494	!
2495	!-- Calculate tendencies for the next Runge-Kutta step
2496	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2497	IF ( intermediate_timestep_count == 1 ) THEN
2498	DO i = nxl, nxr
2499	DO j = nys, nyn
2500	DO k = nzb+1, nzt
2501	tnr_m(k,j,i) = tend(k,j,i)
2502	ENDDO
2503	ENDDO
2504	ENDDO
2505	ELSEIF ( intermediate_timestep_count < &
2506	intermediate_timestep_count_max ) THEN
2507	DO i = nxl, nxr
2508	DO j = nys, nyn
2509	DO k = nzb+1, nzt
2510	tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2511	+ 5.3125_wp * tnr_m(k,j,i)
2512	ENDDO
2513	ENDDO
2514	ENDDO
2515	ENDIF
2516	ENDIF
2517
2518	CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
2519
2520	ENDIF
2521
2522	ENDIF
2523	!
2524	!-- If required, compute prognostic equation for scalar
2525	IF ( passive_scalar ) THEN
2526
2527	CALL cpu_log( log_point(66), 's-equation', 'start' )
2528
2529	!
2530	!-- Scalar/q-tendency terms with communication
2531	sbt = tsc(2)
2532	IF ( scalar_advec == 'bc-scheme' ) THEN
2533
2534	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2535	!
2536	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2537	sbt = 1.0_wp
2538	ENDIF
2539	tend = 0.0_wp
2540	CALL advec_s_bc( s, 's' )
2541
2542	ENDIF
2543
2544	!
2545	!-- Scalar/q-tendency terms with no communication
2546	IF ( scalar_advec /= 'bc-scheme' ) THEN
2547	tend = 0.0_wp
2548	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2549	IF ( ws_scheme_sca ) THEN
2550	CALL advec_s_ws( s, 's' )
2551	ELSE
2552	CALL advec_s_pw( s )
2553	ENDIF
2554	ELSE
2555	CALL advec_s_up( s )
2556	ENDIF
2557	ENDIF
2558
2559	CALL diffusion_s( s, &
2560	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
2561	surf_def_h(2)%ssws, &
2562	surf_lsm_h%ssws, surf_usm_h%ssws, &
2563	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
2564	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
2565	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
2566	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
2567	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
2568	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2569
2570	!
2571	!-- Sink or source of humidity due to canopy elements
2572	IF ( plant_canopy ) CALL pcm_tendency( 7 )
2573
2574	!
2575	!-- Large scale advection. Not implemented for scalars so far.
2576	! IF ( large_scale_forcing ) THEN
2577	! CALL ls_advec( simulated_time, 'q' )
2578	! ENDIF
2579
2580	!
2581	!-- Nudging. Not implemented for scalars so far.
2582	! IF ( nudging ) CALL nudge( simulated_time, 'q' )
2583
2584	!
2585	!-- If required compute influence of large-scale subsidence/ascent.
2586	!-- Not implemented for scalars so far.
2587	IF ( large_scale_subsidence .AND. &
2588	.NOT. use_subsidence_tendencies .AND. &
2589	.NOT. large_scale_forcing ) THEN
2590	CALL subsidence( tend, s, s_init, 3 )
2591	ENDIF
2592
2593	CALL module_interface_actions( 's-tendency' )
2594
2595	!
2596	!-- Prognostic equation for total water content
2597	DO i = nxl, nxr
2598	DO j = nys, nyn
2599	DO k = nzb+1, nzt
2600	s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2601	tsc(3) * ts_m(k,j,i) ) &
2602	- tsc(5) * rdf_sc(k) * &
2603	( s(k,j,i) - s_init(k) ) &
2604	) &
2605	* MERGE( 1.0_wp, 0.0_wp, &
2606	BTEST( wall_flags_0(k,j,i), 0 ) &
2607	)
2608	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
2609	ENDDO
2610	ENDDO
2611	ENDDO
2612
2613	!
2614	!-- Calculate tendencies for the next Runge-Kutta step
2615	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2616	IF ( intermediate_timestep_count == 1 ) THEN
2617	DO i = nxl, nxr
2618	DO j = nys, nyn
2619	DO k = nzb+1, nzt
2620	ts_m(k,j,i) = tend(k,j,i)
2621	ENDDO
2622	ENDDO
2623	ENDDO
2624	ELSEIF ( intermediate_timestep_count < &
2625	intermediate_timestep_count_max ) THEN
2626	DO i = nxl, nxr
2627	DO j = nys, nyn
2628	DO k = nzb+1, nzt
2629	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2630	+ 5.3125_wp * ts_m(k,j,i)
2631	ENDDO
2632	ENDDO
2633	ENDDO
2634	ENDIF
2635	ENDIF
2636
2637	CALL cpu_log( log_point(66), 's-equation', 'stop' )
2638
2639	ENDIF
2640
2641	!
2642	!-- Calculate prognostic equations for turbulence closure
2643	CALL tcm_prognostic_equations()
2644
2645	!
2646	!-- Calculate prognostic equation for chemical quantites
2647	IF ( air_chemistry ) THEN
2648	CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2649	!
2650	!-- Loop over chemical species
2651	DO lsp = 1, nvar
2652	CALL chem_prognostic_equations( chem_species(lsp)%conc_p, &
2653	chem_species(lsp)%conc, &
2654	chem_species(lsp)%tconc_m, &
2655	chem_species(lsp)%conc_pr_init, &
2656	lsp )
2657	ENDDO
2658
2659	CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )
2660	ENDIF ! Chemicals equations
2661
2662	IF ( salsa ) THEN
2663	CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'start' )
2664	!
2665	!-- Loop over aerosol size bins: number and mass bins
2666	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
2667
2668	DO b = 1, nbins
2669	sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
2670	CALL salsa_tendency( 'aerosol_number', aerosol_number(b)%conc_p, &
2671	aerosol_number(b)%conc, &
2672	aerosol_number(b)%tconc_m, &
2673	b, b, aerosol_number(b)%init )
2674	aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
2675	DO c = 1, ncc_tot
2676	sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
2677	CALL salsa_tendency( 'aerosol_mass', &
2678	aerosol_mass((c-1)*nbins+b)%conc_p, &
2679	aerosol_mass((c-1)*nbins+b)%conc, &
2680	aerosol_mass((c-1)*nbins+b)%tconc_m, &
2681	b, c, aerosol_mass((c-1)*nbins+b)%init )
2682	aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
2683	ENDDO
2684	ENDDO
2685	IF ( .NOT. salsa_gases_from_chem ) THEN
2686	DO g = 1, ngast
2687	sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
2688	CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p, &
2689	salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
2690	g, g, salsa_gas(g)%init )
2691	salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
2692	ENDDO
2693	ENDIF
2694
2695	ENDIF
2696
2697	CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'stop' )
2698	ENDIF
2699
2700	!
2701	!-- Calculate prognostic equations for the ocean
2702	IF ( ocean_mode ) CALL ocean_prognostic_equations()
2703
2704	END SUBROUTINE prognostic_equations_vector
2705
2706
2707	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |