Home

Context Navigation

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

Last change on this file since 3754 was 3719, checked in by kanani, 5 years ago
Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)
Property svn:keywords set to `Id`
File size: 104.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |