Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |