Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |