Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |