Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |