Home

Context Navigation

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

Last change on this file since 2219 was 2193, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 66.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |