Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |