Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |