Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |