Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |