Home

Context Navigation

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

Last change on this file since 1409 was 1409, checked in by suehring, 11 years ago
Bugfix: set inflow boundary conditions also if no humidity or passive_scalar is used
Property svn:keywords set to `Id`
File size: 83.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |