Home

Context Navigation

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

Last change on this file since 2116 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 92.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |