Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |