Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |