Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |