Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |