Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |