Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |