Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |