Home

Context Navigation

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

Last change on this file since 2723 was 2719, checked in by maronga, 7 years ago
bugfix for last commit
Property svn:keywords set to `Id`
File size: 93.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |