Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |