Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |