Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |