Home

Context Navigation

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

Last change on this file since 1128 was 1128, checked in by raasch, 11 years ago
asynchronous transfer of ghost point data for acc-optimized version
Property svn:keywords set to `Id`
File size: 65.0 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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	! those parts requiring global communication moved to time_integration,
23	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
24	! j_north
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: prognostic_equations.f90 1128 2013-04-12 06:19:32Z raasch $
29	!
30	! 1115 2013-03-26 18:16:16Z hoffmann
31	! optimized cloud physics: calculation of microphysical tendencies transfered
32	! to microphysics.f90; qr and nr are only calculated if precipitation is required
33	!
34	! 1111 2013-03-08 23:54:10Z raasch
35	! update directives for prognostic quantities removed
36	!
37	! 1106 2013-03-04 05:31:38Z raasch
38	! small changes in code formatting
39	!
40	! 1092 2013-02-02 11:24:22Z raasch
41	! unused variables removed
42	!
43	! 1053 2012-11-13 17:11:03Z hoffmann
44	! implementation of two new prognostic equations for rain drop concentration (nr)
45	! and rain water content (qr)
46	!
47	! currently, only available for cache loop optimization
48	!
49	! 1036 2012-10-22 13:43:42Z raasch
50	! code put under GPL (PALM 3.9)
51	!
52	! 1019 2012-09-28 06:46:45Z raasch
53	! non-optimized version of prognostic_equations removed
54	!
55	! 1015 2012-09-27 09:23:24Z raasch
56	! new branch prognostic_equations_acc
57	! OpenACC statements added + code changes required for GPU optimization
58	!
59	! 1001 2012-09-13 14:08:46Z raasch
60	! all actions concerning leapfrog- and upstream-spline-scheme removed
61	!
62	! 978 2012-08-09 08:28:32Z fricke
63	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
64	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
65	! boundary in case of non-cyclic lateral boundaries
66	! Bugfix: first thread index changes for WS-scheme at the inflow
67	!
68	! 940 2012-07-09 14:31:00Z raasch
69	! temperature equation can be switched off
70	!
71	! 785 2011-11-28 09:47:19Z raasch
72	! new factor rdf_sc allows separate Rayleigh damping of scalars
73	!
74	! 736 2011-08-17 14:13:26Z suehring
75	! Bugfix: determination of first thread index i for WS-scheme
76	!
77	! 709 2011-03-30 09:31:40Z raasch
78	! formatting adjustments
79	!
80	! 673 2011-01-18 16:19:48Z suehring
81	! Consideration of the pressure gradient (steered by tsc(4)) during the time
82	! integration removed.
83	!
84	! 667 2010-12-23 12:06:00Z suehring/gryschka
85	! Calls of the advection routines with WS5 added.
86	! Calls of ws_statistics added to set the statistical arrays to zero after each
87	! time step.
88	!
89	! 531 2010-04-21 06:47:21Z heinze
90	! add call of subsidence in the equation for humidity / passive scalar
91	!
92	! 411 2009-12-11 14:15:58Z heinze
93	! add call of subsidence in the equation for potential temperature
94	!
95	! 388 2009-09-23 09:40:33Z raasch
96	! prho is used instead of rho in diffusion_e,
97	! external pressure gradient
98	!
99	! 153 2008-03-19 09:41:30Z steinfeld
100	! add call of plant_canopy_model in the prognostic equation for
101	! the potential temperature and for the passive scalar
102	!
103	! 138 2007-11-28 10:03:58Z letzel
104	! add call of subroutines that evaluate the canopy drag terms,
105	! add wall_*flux to parameter list of calls of diffusion_s
106	!
107	! 106 2007-08-16 14:30:26Z raasch
108	! +uswst, vswst as arguments in calls of diffusion_u\|v,
109	! loops for u and v are starting from index nxlu, nysv, respectively (needed
110	! for non-cyclic boundary conditions)
111	!
112	! 97 2007-06-21 08:23:15Z raasch
113	! prognostic equation for salinity, density is calculated from equation of
114	! state for seawater and is used for calculation of buoyancy,
115	! +eqn_state_seawater_mod
116	! diffusion_e is called with argument rho in case of ocean runs,
117	! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
118	! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
119	! calc_mean_profile
120	!
121	! 75 2007-03-22 09:54:05Z raasch
122	! checking for negative q and limiting for positive values,
123	! z0 removed from arguments in calls of diffusion_u/v/w, uxrp, vynp eliminated,
124	! subroutine names changed to .._noopt, .._cache, and .._vector,
125	! moisture renamed humidity, Bott-Chlond-scheme can be used in the
126	! _vector-version
127	!
128	! 19 2007-02-23 04:53:48Z raasch
129	! Calculation of e, q, and pt extended for gridpoint nzt,
130	! handling of given temperature/humidity/scalar fluxes at top surface
131	!
132	! RCS Log replace by Id keyword, revision history cleaned up
133	!
134	! Revision 1.21 2006/08/04 15:01:07 raasch
135	! upstream scheme can be forced to be used for tke (use_upstream_for_tke)
136	! regardless of the timestep scheme used for the other quantities,
137	! new argument diss in call of diffusion_e
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	USE control_parameters
150	USE cpulog
151	USE eqn_state_seawater_mod
152	USE grid_variables
153	USE indices
154	USE interfaces
155	USE pegrid
156	USE pointer_interfaces
157	USE statistics
158	USE advec_ws
159	USE advec_s_pw_mod
160	USE advec_s_up_mod
161	USE advec_u_pw_mod
162	USE advec_u_up_mod
163	USE advec_v_pw_mod
164	USE advec_v_up_mod
165	USE advec_w_pw_mod
166	USE advec_w_up_mod
167	USE buoyancy_mod
168	USE calc_precipitation_mod
169	USE calc_radiation_mod
170	USE coriolis_mod
171	USE diffusion_e_mod
172	USE diffusion_s_mod
173	USE diffusion_u_mod
174	USE diffusion_v_mod
175	USE diffusion_w_mod
176	USE impact_of_latent_heat_mod
177	USE microphysics_mod
178	USE plant_canopy_model_mod
179	USE production_e_mod
180	USE subsidence_mod
181	USE user_actions_mod
182
183
184	PRIVATE
185	PUBLIC prognostic_equations_cache, prognostic_equations_vector, &
186	prognostic_equations_acc
187
188	INTERFACE prognostic_equations_cache
189	MODULE PROCEDURE prognostic_equations_cache
190	END INTERFACE prognostic_equations_cache
191
192	INTERFACE prognostic_equations_vector
193	MODULE PROCEDURE prognostic_equations_vector
194	END INTERFACE prognostic_equations_vector
195
196	INTERFACE prognostic_equations_acc
197	MODULE PROCEDURE prognostic_equations_acc
198	END INTERFACE prognostic_equations_acc
199
200
201	CONTAINS
202
203
204	SUBROUTINE prognostic_equations_cache
205
206	!------------------------------------------------------------------------------!
207	! Version with one optimized loop over all equations. It is only allowed to
208	! be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
209	!
210	! Here the calls of most subroutines are embedded in two DO loops over i and j,
211	! so communication between CPUs is not allowed (does not make sense) within
212	! these loops.
213	!
214	! (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
215	!------------------------------------------------------------------------------!
216
217	IMPLICIT NONE
218
219	INTEGER :: i, i_omp_start, j, k, omp_get_thread_num, tn = 0
220	LOGICAL :: loop_start
221
222
223	!
224	!-- Time measurement can only be performed for the whole set of equations
225	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
226
227	!
228	!-- Loop over all prognostic equations
229	!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
230
231	!$ tn = omp_get_thread_num()
232	loop_start = .TRUE.
233	!$OMP DO
234	DO i = nxl, nxr
235
236	!
237	!-- Store the first loop index. It differs for each thread and is required
238	!-- later in advec_ws
239	IF ( loop_start ) THEN
240	loop_start = .FALSE.
241	i_omp_start = i
242	ENDIF
243
244	DO j = nys, nyn
245	!
246	!-- Tendency terms for u-velocity component
247	IF ( .NOT. outflow_l .OR. i > nxl ) THEN
248
249	tend(:,j,i) = 0.0
250	IF ( timestep_scheme(1:5) == 'runge' ) THEN
251	IF ( ws_scheme_mom ) THEN
252	IF ( ( inflow_l .OR. outflow_l ) .AND. i_omp_start == nxl ) THEN
253	CALL advec_u_ws( i, j, i_omp_start + 1, tn )
254	ELSE
255	CALL advec_u_ws( i, j, i_omp_start, tn )
256	ENDIF
257	ELSE
258	CALL advec_u_pw( i, j )
259	ENDIF
260	ELSE
261	CALL advec_u_up( i, j )
262	ENDIF
263	CALL diffusion_u( i, j )
264	CALL coriolis( i, j, 1 )
265	IF ( sloping_surface .AND. .NOT. neutral ) THEN
266	CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
267	ENDIF
268
269	!
270	!-- Drag by plant canopy
271	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 1 )
272
273	!
274	!-- External pressure gradient
275	IF ( dp_external ) THEN
276	DO k = dp_level_ind_b+1, nzt
277	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
278	ENDDO
279	ENDIF
280
281	CALL user_actions( i, j, 'u-tendency' )
282	!
283	!-- Prognostic equation for u-velocity component
284	DO k = nzb_u_inner(j,i)+1, nzt
285	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
286	tsc(3) * tu_m(k,j,i) ) &
287	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
288	ENDDO
289
290	!
291	!-- Calculate tendencies for the next Runge-Kutta step
292	IF ( timestep_scheme(1:5) == 'runge' ) THEN
293	IF ( intermediate_timestep_count == 1 ) THEN
294	DO k = nzb_u_inner(j,i)+1, nzt
295	tu_m(k,j,i) = tend(k,j,i)
296	ENDDO
297	ELSEIF ( intermediate_timestep_count < &
298	intermediate_timestep_count_max ) THEN
299	DO k = nzb_u_inner(j,i)+1, nzt
300	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
301	ENDDO
302	ENDIF
303	ENDIF
304
305	ENDIF
306
307	!
308	!-- Tendency terms for v-velocity component
309	IF ( .NOT. outflow_s .OR. j > nys ) THEN
310
311	tend(:,j,i) = 0.0
312	IF ( timestep_scheme(1:5) == 'runge' ) THEN
313	IF ( ws_scheme_mom ) THEN
314	CALL advec_v_ws( i, j, i_omp_start, tn )
315	ELSE
316	CALL advec_v_pw( i, j )
317	ENDIF
318	ELSE
319	CALL advec_v_up( i, j )
320	ENDIF
321	CALL diffusion_v( i, j )
322	CALL coriolis( i, j, 2 )
323
324	!
325	!-- Drag by plant canopy
326	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 2 )
327
328	!
329	!-- External pressure gradient
330	IF ( dp_external ) THEN
331	DO k = dp_level_ind_b+1, nzt
332	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
333	ENDDO
334	ENDIF
335
336	CALL user_actions( i, j, 'v-tendency' )
337	!
338	!-- Prognostic equation for v-velocity component
339	DO k = nzb_v_inner(j,i)+1, nzt
340	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
341	tsc(3) * tv_m(k,j,i) ) &
342	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
343	ENDDO
344
345	!
346	!-- Calculate tendencies for the next Runge-Kutta step
347	IF ( timestep_scheme(1:5) == 'runge' ) THEN
348	IF ( intermediate_timestep_count == 1 ) THEN
349	DO k = nzb_v_inner(j,i)+1, nzt
350	tv_m(k,j,i) = tend(k,j,i)
351	ENDDO
352	ELSEIF ( intermediate_timestep_count < &
353	intermediate_timestep_count_max ) THEN
354	DO k = nzb_v_inner(j,i)+1, nzt
355	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
356	ENDDO
357	ENDIF
358	ENDIF
359
360	ENDIF
361
362	!
363	!-- Tendency terms for w-velocity component
364	tend(:,j,i) = 0.0
365	IF ( timestep_scheme(1:5) == 'runge' ) THEN
366	IF ( ws_scheme_mom ) THEN
367	CALL advec_w_ws( i, j, i_omp_start, tn )
368	ELSE
369	CALL advec_w_pw( i, j )
370	END IF
371	ELSE
372	CALL advec_w_up( i, j )
373	ENDIF
374	CALL diffusion_w( i, j )
375	CALL coriolis( i, j, 3 )
376
377	IF ( .NOT. neutral ) THEN
378	IF ( ocean ) THEN
379	CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
380	ELSE
381	IF ( .NOT. humidity ) THEN
382	CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
383	ELSE
384	CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
385	ENDIF
386	ENDIF
387	ENDIF
388
389	!
390	!-- Drag by plant canopy
391	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 3 )
392
393	CALL user_actions( i, j, 'w-tendency' )
394
395	!
396	!-- Prognostic equation for w-velocity component
397	DO k = nzb_w_inner(j,i)+1, nzt-1
398	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
399	tsc(3) * tw_m(k,j,i) ) &
400	- tsc(5) * rdf(k) * w(k,j,i)
401	ENDDO
402
403	!
404	!-- Calculate tendencies for the next Runge-Kutta step
405	IF ( timestep_scheme(1:5) == 'runge' ) THEN
406	IF ( intermediate_timestep_count == 1 ) THEN
407	DO k = nzb_w_inner(j,i)+1, nzt-1
408	tw_m(k,j,i) = tend(k,j,i)
409	ENDDO
410	ELSEIF ( intermediate_timestep_count < &
411	intermediate_timestep_count_max ) THEN
412	DO k = nzb_w_inner(j,i)+1, nzt-1
413	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
414	ENDDO
415	ENDIF
416	ENDIF
417	!
418	!-- If required, calculate tendencies for total water content, liquid water
419	!-- potential temperature, rain water content and rain drop concentration
420	IF ( cloud_physics .AND. icloud_scheme == 0 ) CALL microphysics_control( i, j )
421	!
422	!-- If required, compute prognostic equation for potential temperature
423	IF ( .NOT. neutral ) THEN
424	!
425	!-- Tendency terms for potential temperature
426	tend(:,j,i) = 0.0
427	IF ( timestep_scheme(1:5) == 'runge' ) THEN
428	IF ( ws_scheme_sca ) THEN
429	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
430	flux_l_pt, diss_l_pt, i_omp_start, tn )
431	ELSE
432	CALL advec_s_pw( i, j, pt )
433	ENDIF
434	ELSE
435	CALL advec_s_up( i, j, pt )
436	ENDIF
437	CALL diffusion_s( i, j, pt, shf, tswst, wall_heatflux )
438
439	!
440	!-- If required compute heating/cooling due to long wave radiation
441	!-- processes
442	IF ( radiation ) THEN
443	CALL calc_radiation( i, j )
444	ENDIF
445
446	!
447	!-- Using microphysical tendencies (latent heat)
448	IF ( cloud_physics ) THEN
449	IF ( icloud_scheme == 0 ) THEN
450	tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i)
451	ELSEIF ( icloud_scheme == 1 .AND. precipitation ) THEN
452	CALL impact_of_latent_heat( i, j )
453	ENDIF
454	ENDIF
455
456	!
457	!-- Consideration of heat sources within the plant canopy
458	IF ( plant_canopy .AND. cthf /= 0.0 ) THEN
459	CALL plant_canopy_model( i, j, 4 )
460	ENDIF
461
462	!
463	!-- If required, compute effect of large-scale subsidence/ascent
464	IF ( large_scale_subsidence ) THEN
465	CALL subsidence( i, j, tend, pt, pt_init )
466	ENDIF
467
468	CALL user_actions( i, j, 'pt-tendency' )
469
470	!
471	!-- Prognostic equation for potential temperature
472	DO k = nzb_s_inner(j,i)+1, nzt
473	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
474	tsc(3) * tpt_m(k,j,i) ) &
475	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
476	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
477	ENDDO
478
479	!
480	!-- Calculate tendencies for the next Runge-Kutta step
481	IF ( timestep_scheme(1:5) == 'runge' ) THEN
482	IF ( intermediate_timestep_count == 1 ) THEN
483	DO k = nzb_s_inner(j,i)+1, nzt
484	tpt_m(k,j,i) = tend(k,j,i)
485	ENDDO
486	ELSEIF ( intermediate_timestep_count < &
487	intermediate_timestep_count_max ) THEN
488	DO k = nzb_s_inner(j,i)+1, nzt
489	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
490	5.3125 * tpt_m(k,j,i)
491	ENDDO
492	ENDIF
493	ENDIF
494
495	ENDIF
496
497	!
498	!-- If required, compute prognostic equation for salinity
499	IF ( ocean ) THEN
500
501	!
502	!-- Tendency-terms for salinity
503	tend(:,j,i) = 0.0
504	IF ( timestep_scheme(1:5) == 'runge' ) &
505	THEN
506	IF ( ws_scheme_sca ) THEN
507	CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, &
508	diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn )
509	ELSE
510	CALL advec_s_pw( i, j, sa )
511	ENDIF
512	ELSE
513	CALL advec_s_up( i, j, sa )
514	ENDIF
515	CALL diffusion_s( i, j, sa, saswsb, saswst, wall_salinityflux )
516
517	CALL user_actions( i, j, 'sa-tendency' )
518
519	!
520	!-- Prognostic equation for salinity
521	DO k = nzb_s_inner(j,i)+1, nzt
522	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
523	tsc(3) * tsa_m(k,j,i) ) &
524	- tsc(5) * rdf_sc(k) * &
525	( sa(k,j,i) - sa_init(k) )
526	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
527	ENDDO
528
529	!
530	!-- Calculate tendencies for the next Runge-Kutta step
531	IF ( timestep_scheme(1:5) == 'runge' ) THEN
532	IF ( intermediate_timestep_count == 1 ) THEN
533	DO k = nzb_s_inner(j,i)+1, nzt
534	tsa_m(k,j,i) = tend(k,j,i)
535	ENDDO
536	ELSEIF ( intermediate_timestep_count < &
537	intermediate_timestep_count_max ) THEN
538	DO k = nzb_s_inner(j,i)+1, nzt
539	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
540	5.3125 * tsa_m(k,j,i)
541	ENDDO
542	ENDIF
543	ENDIF
544
545	!
546	!-- Calculate density by the equation of state for seawater
547	CALL eqn_state_seawater( i, j )
548
549	ENDIF
550
551	!
552	!-- If required, compute prognostic equation for total water content /
553	!-- scalar
554	IF ( humidity .OR. passive_scalar ) THEN
555
556	!
557	!-- Tendency-terms for total water content / scalar
558	tend(:,j,i) = 0.0
559	IF ( timestep_scheme(1:5) == 'runge' ) &
560	THEN
561	IF ( ws_scheme_sca ) THEN
562	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
563	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
564	ELSE
565	CALL advec_s_pw( i, j, q )
566	ENDIF
567	ELSE
568	CALL advec_s_up( i, j, q )
569	ENDIF
570	CALL diffusion_s( i, j, q, qsws, qswst, wall_qflux )
571
572	!
573	!-- Using microphysical tendencies
574	IF ( cloud_physics ) THEN
575	IF ( icloud_scheme == 0 ) THEN
576	tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i)
577	ELSEIF ( icloud_scheme == 1 .AND. precipitation ) THEN
578	CALL calc_precipitation( i, j )
579	ENDIF
580	ENDIF
581	!
582	!-- Sink or source of scalar concentration due to canopy elements
583	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
584
585	!
586	!-- If required compute influence of large-scale subsidence/ascent
587	IF ( large_scale_subsidence ) THEN
588	CALL subsidence( i, j, tend, q, q_init )
589	ENDIF
590
591	CALL user_actions( i, j, 'q-tendency' )
592
593	!
594	!-- Prognostic equation for total water content / scalar
595	DO k = nzb_s_inner(j,i)+1, nzt
596	q_p(k,j,i) = q(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
597	tsc(3) * tq_m(k,j,i) ) &
598	- tsc(5) * rdf_sc(k) * &
599	( q(k,j,i) - q_init(k) )
600	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
601	ENDDO
602
603	!
604	!-- Calculate tendencies for the next Runge-Kutta step
605	IF ( timestep_scheme(1:5) == 'runge' ) THEN
606	IF ( intermediate_timestep_count == 1 ) THEN
607	DO k = nzb_s_inner(j,i)+1, nzt
608	tq_m(k,j,i) = tend(k,j,i)
609	ENDDO
610	ELSEIF ( intermediate_timestep_count < &
611	intermediate_timestep_count_max ) THEN
612	DO k = nzb_s_inner(j,i)+1, nzt
613	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + &
614	5.3125 * tq_m(k,j,i)
615	ENDDO
616	ENDIF
617	ENDIF
618
619	!
620	!-- If required, calculate prognostic equations for rain water content
621	!-- and rain drop concentration
622	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
623	precipitation ) THEN
624	!
625	!-- Calculate prognostic equation for rain water content
626	tend(:,j,i) = 0.0
627	IF ( timestep_scheme(1:5) == 'runge' ) &
628	THEN
629	IF ( ws_scheme_sca ) THEN
630	CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr, &
631	diss_s_qr, flux_l_qr, diss_l_qr, &
632	i_omp_start, tn )
633	ELSE
634	CALL advec_s_pw( i, j, qr )
635	ENDIF
636	ELSE
637	CALL advec_s_up( i, j, qr )
638	ENDIF
639	CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux )
640	!
641	!-- Using microphysical tendencies (autoconversion, accretion,
642	!-- evaporation; if required: sedimentation)
643	tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)
644
645	CALL user_actions( i, j, 'qr-tendency' )
646	!
647	!-- Prognostic equation for rain water content
648	DO k = nzb_s_inner(j,i)+1, nzt
649	qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
650	tsc(3) * tqr_m(k,j,i) ) &
651	- tsc(5) * rdf_sc(k) * qr(k,j,i)
652	IF ( qr_p(k,j,i) < 0.0 ) qr_p(k,j,i) = 0.0
653	ENDDO
654	!
655	!-- Calculate tendencies for the next Runge-Kutta step
656	IF ( timestep_scheme(1:5) == 'runge' ) THEN
657	IF ( intermediate_timestep_count == 1 ) THEN
658	DO k = nzb_s_inner(j,i)+1, nzt
659	tqr_m(k,j,i) = tend(k,j,i)
660	ENDDO
661	ELSEIF ( intermediate_timestep_count < &
662	intermediate_timestep_count_max ) THEN
663	DO k = nzb_s_inner(j,i)+1, nzt
664	tqr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
665	5.3125 * tqr_m(k,j,i)
666	ENDDO
667	ENDIF
668	ENDIF
669
670	!
671	!-- Calculate prognostic equation for rain drop concentration.
672	tend(:,j,i) = 0.0
673	IF ( timestep_scheme(1:5) == 'runge' ) THEN
674	IF ( ws_scheme_sca ) THEN
675	CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr, &
676	diss_s_nr, flux_l_nr, diss_l_nr, &
677	i_omp_start, tn )
678	ELSE
679	CALL advec_s_pw( i, j, nr )
680	ENDIF
681	ELSE
682	CALL advec_s_up( i, j, nr )
683	ENDIF
684	CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux )
685	!
686	!-- Using microphysical tendencies (autoconversion, accretion,
687	!-- selfcollection, breakup, evaporation;
688	!-- if required: sedimentation)
689	tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)
690
691	CALL user_actions( i, j, 'nr-tendency' )
692	!
693	!-- Prognostic equation for rain drop concentration
694	DO k = nzb_s_inner(j,i)+1, nzt
695	nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
696	tsc(3) * tnr_m(k,j,i) ) &
697	- tsc(5) * rdf_sc(k) * nr(k,j,i)
698	IF ( nr_p(k,j,i) < 0.0 ) nr_p(k,j,i) = 0.0
699	ENDDO
700	!
701	!-- Calculate tendencies for the next Runge-Kutta step
702	IF ( timestep_scheme(1:5) == 'runge' ) THEN
703	IF ( intermediate_timestep_count == 1 ) THEN
704	DO k = nzb_s_inner(j,i)+1, nzt
705	tnr_m(k,j,i) = tend(k,j,i)
706	ENDDO
707	ELSEIF ( intermediate_timestep_count < &
708	intermediate_timestep_count_max ) THEN
709	DO k = nzb_s_inner(j,i)+1, nzt
710	tnr_m(k,j,i) = -9.5625 * tend(k,j,i) + &
711	5.3125 * tnr_m(k,j,i)
712	ENDDO
713	ENDIF
714	ENDIF
715
716	ENDIF
717
718	ENDIF
719
720	!
721	!-- If required, compute prognostic equation for turbulent kinetic
722	!-- energy (TKE)
723	IF ( .NOT. constant_diffusion ) THEN
724
725	!
726	!-- Tendency-terms for TKE
727	tend(:,j,i) = 0.0
728	IF ( timestep_scheme(1:5) == 'runge' &
729	.AND. .NOT. use_upstream_for_tke ) THEN
730	IF ( ws_scheme_sca ) THEN
731	CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
732	flux_l_e, diss_l_e , i_omp_start, tn )
733	ELSE
734	CALL advec_s_pw( i, j, e )
735	ENDIF
736	ELSE
737	CALL advec_s_up( i, j, e )
738	ENDIF
739	IF ( .NOT. humidity ) THEN
740	IF ( ocean ) THEN
741	CALL diffusion_e( i, j, prho, prho_reference )
742	ELSE
743	CALL diffusion_e( i, j, pt, pt_reference )
744	ENDIF
745	ELSE
746	CALL diffusion_e( i, j, vpt, pt_reference )
747	ENDIF
748	CALL production_e( i, j )
749
750	!
751	!-- Additional sink term for flows through plant canopies
752	IF ( plant_canopy ) CALL plant_canopy_model( i, j, 6 )
753
754	CALL user_actions( i, j, 'e-tendency' )
755
756	!
757	!-- Prognostic equation for TKE.
758	!-- Eliminate negative TKE values, which can occur due to numerical
759	!-- reasons in the course of the integration. In such cases the old
760	!-- TKE value is reduced by 90%.
761	DO k = nzb_s_inner(j,i)+1, nzt
762	e_p(k,j,i) = e(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
763	tsc(3) * te_m(k,j,i) )
764	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
765	ENDDO
766
767	!
768	!-- Calculate tendencies for the next Runge-Kutta step
769	IF ( timestep_scheme(1:5) == 'runge' ) THEN
770	IF ( intermediate_timestep_count == 1 ) THEN
771	DO k = nzb_s_inner(j,i)+1, nzt
772	te_m(k,j,i) = tend(k,j,i)
773	ENDDO
774	ELSEIF ( intermediate_timestep_count < &
775	intermediate_timestep_count_max ) THEN
776	DO k = nzb_s_inner(j,i)+1, nzt
777	te_m(k,j,i) = -9.5625 * tend(k,j,i) + &
778	5.3125 * te_m(k,j,i)
779	ENDDO
780	ENDIF
781	ENDIF
782
783	ENDIF ! TKE equation
784
785	ENDDO
786	ENDDO
787	!$OMP END PARALLEL
788
789	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
790
791
792	END SUBROUTINE prognostic_equations_cache
793
794
795	SUBROUTINE prognostic_equations_vector
796
797	!------------------------------------------------------------------------------!
798	! Version for vector machines
799	!------------------------------------------------------------------------------!
800
801	IMPLICIT NONE
802
803	INTEGER :: i, j, k
804	REAL :: sbt
805
806
807	!
808	!-- u-velocity component
809	CALL cpu_log( log_point(5), 'u-equation', 'start' )
810
811	tend = 0.0
812	IF ( timestep_scheme(1:5) == 'runge' ) THEN
813	IF ( ws_scheme_mom ) THEN
814	CALL advec_u_ws
815	ELSE
816	CALL advec_u_pw
817	ENDIF
818	ELSE
819	CALL advec_u_up
820	ENDIF
821	CALL diffusion_u
822	CALL coriolis( 1 )
823	IF ( sloping_surface .AND. .NOT. neutral ) THEN
824	CALL buoyancy( pt, pt_reference, 1, 4 )
825	ENDIF
826
827	!
828	!-- Drag by plant canopy
829	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
830
831	!
832	!-- External pressure gradient
833	IF ( dp_external ) THEN
834	DO i = nxlu, nxr
835	DO j = nys, nyn
836	DO k = dp_level_ind_b+1, nzt
837	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
838	ENDDO
839	ENDDO
840	ENDDO
841	ENDIF
842
843	CALL user_actions( 'u-tendency' )
844
845	!
846	!-- Prognostic equation for u-velocity component
847	DO i = nxlu, nxr
848	DO j = nys, nyn
849	DO k = nzb_u_inner(j,i)+1, nzt
850	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
851	tsc(3) * tu_m(k,j,i) ) &
852	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
853	ENDDO
854	ENDDO
855	ENDDO
856
857	!
858	!-- Calculate tendencies for the next Runge-Kutta step
859	IF ( timestep_scheme(1:5) == 'runge' ) THEN
860	IF ( intermediate_timestep_count == 1 ) THEN
861	DO i = nxlu, nxr
862	DO j = nys, nyn
863	DO k = nzb_u_inner(j,i)+1, nzt
864	tu_m(k,j,i) = tend(k,j,i)
865	ENDDO
866	ENDDO
867	ENDDO
868	ELSEIF ( intermediate_timestep_count < &
869	intermediate_timestep_count_max ) THEN
870	DO i = nxlu, nxr
871	DO j = nys, nyn
872	DO k = nzb_u_inner(j,i)+1, nzt
873	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
874	ENDDO
875	ENDDO
876	ENDDO
877	ENDIF
878	ENDIF
879
880	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
881
882	!
883	!-- v-velocity component
884	CALL cpu_log( log_point(6), 'v-equation', 'start' )
885
886	tend = 0.0
887	IF ( timestep_scheme(1:5) == 'runge' ) THEN
888	IF ( ws_scheme_mom ) THEN
889	CALL advec_v_ws
890	ELSE
891	CALL advec_v_pw
892	END IF
893	ELSE
894	CALL advec_v_up
895	ENDIF
896	CALL diffusion_v
897	CALL coriolis( 2 )
898
899	!
900	!-- Drag by plant canopy
901	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
902
903	!
904	!-- External pressure gradient
905	IF ( dp_external ) THEN
906	DO i = nxl, nxr
907	DO j = nysv, nyn
908	DO k = dp_level_ind_b+1, nzt
909	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
910	ENDDO
911	ENDDO
912	ENDDO
913	ENDIF
914
915	CALL user_actions( 'v-tendency' )
916
917	!
918	!-- Prognostic equation for v-velocity component
919	DO i = nxl, nxr
920	DO j = nysv, nyn
921	DO k = nzb_v_inner(j,i)+1, nzt
922	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
923	tsc(3) * tv_m(k,j,i) ) &
924	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
925	ENDDO
926	ENDDO
927	ENDDO
928
929	!
930	!-- Calculate tendencies for the next Runge-Kutta step
931	IF ( timestep_scheme(1:5) == 'runge' ) THEN
932	IF ( intermediate_timestep_count == 1 ) THEN
933	DO i = nxl, nxr
934	DO j = nysv, nyn
935	DO k = nzb_v_inner(j,i)+1, nzt
936	tv_m(k,j,i) = tend(k,j,i)
937	ENDDO
938	ENDDO
939	ENDDO
940	ELSEIF ( intermediate_timestep_count < &
941	intermediate_timestep_count_max ) THEN
942	DO i = nxl, nxr
943	DO j = nysv, nyn
944	DO k = nzb_v_inner(j,i)+1, nzt
945	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
946	ENDDO
947	ENDDO
948	ENDDO
949	ENDIF
950	ENDIF
951
952	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
953
954	!
955	!-- w-velocity component
956	CALL cpu_log( log_point(7), 'w-equation', 'start' )
957
958	tend = 0.0
959	IF ( timestep_scheme(1:5) == 'runge' ) THEN
960	IF ( ws_scheme_mom ) THEN
961	CALL advec_w_ws
962	ELSE
963	CALL advec_w_pw
964	ENDIF
965	ELSE
966	CALL advec_w_up
967	ENDIF
968	CALL diffusion_w
969	CALL coriolis( 3 )
970
971	IF ( .NOT. neutral ) THEN
972	IF ( ocean ) THEN
973	CALL buoyancy( rho, rho_reference, 3, 64 )
974	ELSE
975	IF ( .NOT. humidity ) THEN
976	CALL buoyancy( pt, pt_reference, 3, 4 )
977	ELSE
978	CALL buoyancy( vpt, pt_reference, 3, 44 )
979	ENDIF
980	ENDIF
981	ENDIF
982
983	!
984	!-- Drag by plant canopy
985	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
986
987	CALL user_actions( 'w-tendency' )
988
989	!
990	!-- Prognostic equation for w-velocity component
991	DO i = nxl, nxr
992	DO j = nys, nyn
993	DO k = nzb_w_inner(j,i)+1, nzt-1
994	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
995	tsc(3) * tw_m(k,j,i) ) &
996	- tsc(5) * rdf(k) * w(k,j,i)
997	ENDDO
998	ENDDO
999	ENDDO
1000
1001	!
1002	!-- Calculate tendencies for the next Runge-Kutta step
1003	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1004	IF ( intermediate_timestep_count == 1 ) THEN
1005	DO i = nxl, nxr
1006	DO j = nys, nyn
1007	DO k = nzb_w_inner(j,i)+1, nzt-1
1008	tw_m(k,j,i) = tend(k,j,i)
1009	ENDDO
1010	ENDDO
1011	ENDDO
1012	ELSEIF ( intermediate_timestep_count < &
1013	intermediate_timestep_count_max ) THEN
1014	DO i = nxl, nxr
1015	DO j = nys, nyn
1016	DO k = nzb_w_inner(j,i)+1, nzt-1
1017	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1018	ENDDO
1019	ENDDO
1020	ENDDO
1021	ENDIF
1022	ENDIF
1023
1024	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1025
1026
1027	!
1028	!-- If required, compute prognostic equation for potential temperature
1029	IF ( .NOT. neutral ) THEN
1030
1031	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1032
1033	!
1034	!-- pt-tendency terms with communication
1035	sbt = tsc(2)
1036	IF ( scalar_advec == 'bc-scheme' ) THEN
1037
1038	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1039	!
1040	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1041	sbt = 1.0
1042	ENDIF
1043	tend = 0.0
1044	CALL advec_s_bc( pt, 'pt' )
1045
1046	ENDIF
1047
1048	!
1049	!-- pt-tendency terms with no communication
1050	IF ( scalar_advec /= 'bc-scheme' ) THEN
1051	tend = 0.0
1052	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1053	IF ( ws_scheme_sca ) THEN
1054	CALL advec_s_ws( pt, 'pt' )
1055	ELSE
1056	CALL advec_s_pw( pt )
1057	ENDIF
1058	ELSE
1059	CALL advec_s_up( pt )
1060	ENDIF
1061	ENDIF
1062
1063	CALL diffusion_s( pt, shf, tswst, wall_heatflux )
1064
1065	!
1066	!-- If required compute heating/cooling due to long wave radiation processes
1067	IF ( radiation ) THEN
1068	CALL calc_radiation
1069	ENDIF
1070
1071	!
1072	!-- If required compute impact of latent heat due to precipitation
1073	IF ( precipitation ) THEN
1074	CALL impact_of_latent_heat
1075	ENDIF
1076
1077	!
1078	!-- Consideration of heat sources within the plant canopy
1079	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1080	CALL plant_canopy_model( 4 )
1081	ENDIF
1082
1083	!
1084	!-- If required compute influence of large-scale subsidence/ascent
1085	IF ( large_scale_subsidence ) THEN
1086	CALL subsidence( tend, pt, pt_init )
1087	ENDIF
1088
1089	CALL user_actions( 'pt-tendency' )
1090
1091	!
1092	!-- Prognostic equation for potential temperature
1093	DO i = nxl, nxr
1094	DO j = nys, nyn
1095	DO k = nzb_s_inner(j,i)+1, nzt
1096	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1097	tsc(3) * tpt_m(k,j,i) ) &
1098	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
1099	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
1100	ENDDO
1101	ENDDO
1102	ENDDO
1103
1104	!
1105	!-- Calculate tendencies for the next Runge-Kutta step
1106	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1107	IF ( intermediate_timestep_count == 1 ) THEN
1108	DO i = nxl, nxr
1109	DO j = nys, nyn
1110	DO k = nzb_s_inner(j,i)+1, nzt
1111	tpt_m(k,j,i) = tend(k,j,i)
1112	ENDDO
1113	ENDDO
1114	ENDDO
1115	ELSEIF ( intermediate_timestep_count < &
1116	intermediate_timestep_count_max ) THEN
1117	DO i = nxl, nxr
1118	DO j = nys, nyn
1119	DO k = nzb_s_inner(j,i)+1, nzt
1120	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1121	5.3125 * tpt_m(k,j,i)
1122	ENDDO
1123	ENDDO
1124	ENDDO
1125	ENDIF
1126	ENDIF
1127
1128	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1129
1130	ENDIF
1131
1132	!
1133	!-- If required, compute prognostic equation for salinity
1134	IF ( ocean ) THEN
1135
1136	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1137
1138	!
1139	!-- sa-tendency terms with communication
1140	sbt = tsc(2)
1141	IF ( scalar_advec == 'bc-scheme' ) THEN
1142
1143	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1144	!
1145	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1146	sbt = 1.0
1147	ENDIF
1148	tend = 0.0
1149	CALL advec_s_bc( sa, 'sa' )
1150
1151	ENDIF
1152
1153	!
1154	!-- sa-tendency terms with no communication
1155	IF ( scalar_advec /= 'bc-scheme' ) THEN
1156	tend = 0.0
1157	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1158	IF ( ws_scheme_sca ) THEN
1159	CALL advec_s_ws( sa, 'sa' )
1160	ELSE
1161	CALL advec_s_pw( sa )
1162	ENDIF
1163	ELSE
1164	CALL advec_s_up( sa )
1165	ENDIF
1166	ENDIF
1167
1168	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
1169
1170	CALL user_actions( 'sa-tendency' )
1171
1172	!
1173	!-- Prognostic equation for salinity
1174	DO i = nxl, nxr
1175	DO j = nys, nyn
1176	DO k = nzb_s_inner(j,i)+1, nzt
1177	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1178	tsc(3) * tsa_m(k,j,i) ) &
1179	- tsc(5) * rdf_sc(k) * &
1180	( sa(k,j,i) - sa_init(k) )
1181	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1182	ENDDO
1183	ENDDO
1184	ENDDO
1185
1186	!
1187	!-- Calculate tendencies for the next Runge-Kutta step
1188	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1189	IF ( intermediate_timestep_count == 1 ) THEN
1190	DO i = nxl, nxr
1191	DO j = nys, nyn
1192	DO k = nzb_s_inner(j,i)+1, nzt
1193	tsa_m(k,j,i) = tend(k,j,i)
1194	ENDDO
1195	ENDDO
1196	ENDDO
1197	ELSEIF ( intermediate_timestep_count < &
1198	intermediate_timestep_count_max ) THEN
1199	DO i = nxl, nxr
1200	DO j = nys, nyn
1201	DO k = nzb_s_inner(j,i)+1, nzt
1202	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + &
1203	5.3125 * tsa_m(k,j,i)
1204	ENDDO
1205	ENDDO
1206	ENDDO
1207	ENDIF
1208	ENDIF
1209
1210	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1211
1212	!
1213	!-- Calculate density by the equation of state for seawater
1214	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1215	CALL eqn_state_seawater
1216	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1217
1218	ENDIF
1219
1220	!
1221	!-- If required, compute prognostic equation for total water content / scalar
1222	IF ( humidity .OR. passive_scalar ) THEN
1223
1224	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1225
1226	!
1227	!-- Scalar/q-tendency terms with communication
1228	sbt = tsc(2)
1229	IF ( scalar_advec == 'bc-scheme' ) THEN
1230
1231	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1232	!
1233	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1234	sbt = 1.0
1235	ENDIF
1236	tend = 0.0
1237	CALL advec_s_bc( q, 'q' )
1238
1239	ENDIF
1240
1241	!
1242	!-- Scalar/q-tendency terms with no communication
1243	IF ( scalar_advec /= 'bc-scheme' ) THEN
1244	tend = 0.0
1245	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1246	IF ( ws_scheme_sca ) THEN
1247	CALL advec_s_ws( q, 'q' )
1248	ELSE
1249	CALL advec_s_pw( q )
1250	ENDIF
1251	ELSE
1252	CALL advec_s_up( q )
1253	ENDIF
1254	ENDIF
1255
1256	CALL diffusion_s( q, qsws, qswst, wall_qflux )
1257
1258	!
1259	!-- If required compute decrease of total water content due to
1260	!-- precipitation
1261	IF ( precipitation ) THEN
1262	CALL calc_precipitation
1263	ENDIF
1264
1265	!
1266	!-- Sink or source of scalar concentration due to canopy elements
1267	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
1268
1269	!
1270	!-- If required compute influence of large-scale subsidence/ascent
1271	IF ( large_scale_subsidence ) THEN
1272	CALL subsidence( tend, q, q_init )
1273	ENDIF
1274
1275	CALL user_actions( 'q-tendency' )
1276
1277	!
1278	!-- Prognostic equation for total water content / scalar
1279	DO i = nxl, nxr
1280	DO j = nys, nyn
1281	DO k = nzb_s_inner(j,i)+1, nzt
1282	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1283	tsc(3) * tq_m(k,j,i) ) &
1284	- tsc(5) * rdf_sc(k) * &
1285	( q(k,j,i) - q_init(k) )
1286	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1287	ENDDO
1288	ENDDO
1289	ENDDO
1290
1291	!
1292	!-- Calculate tendencies for the next Runge-Kutta step
1293	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1294	IF ( intermediate_timestep_count == 1 ) THEN
1295	DO i = nxl, nxr
1296	DO j = nys, nyn
1297	DO k = nzb_s_inner(j,i)+1, nzt
1298	tq_m(k,j,i) = tend(k,j,i)
1299	ENDDO
1300	ENDDO
1301	ENDDO
1302	ELSEIF ( intermediate_timestep_count < &
1303	intermediate_timestep_count_max ) THEN
1304	DO i = nxl, nxr
1305	DO j = nys, nyn
1306	DO k = nzb_s_inner(j,i)+1, nzt
1307	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1308	ENDDO
1309	ENDDO
1310	ENDDO
1311	ENDIF
1312	ENDIF
1313
1314	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1315
1316	ENDIF
1317
1318	!
1319	!-- If required, compute prognostic equation for turbulent kinetic
1320	!-- energy (TKE)
1321	IF ( .NOT. constant_diffusion ) THEN
1322
1323	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1324
1325	sbt = tsc(2)
1326	IF ( .NOT. use_upstream_for_tke ) THEN
1327	IF ( scalar_advec == 'bc-scheme' ) THEN
1328
1329	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1330	!
1331	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1332	sbt = 1.0
1333	ENDIF
1334	tend = 0.0
1335	CALL advec_s_bc( e, 'e' )
1336
1337	ENDIF
1338	ENDIF
1339
1340	!
1341	!-- TKE-tendency terms with no communication
1342	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
1343	IF ( use_upstream_for_tke ) THEN
1344	tend = 0.0
1345	CALL advec_s_up( e )
1346	ELSE
1347	tend = 0.0
1348	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1349	IF ( ws_scheme_sca ) THEN
1350	CALL advec_s_ws( e, 'e' )
1351	ELSE
1352	CALL advec_s_pw( e )
1353	ENDIF
1354	ELSE
1355	CALL advec_s_up( e )
1356	ENDIF
1357	ENDIF
1358	ENDIF
1359
1360	IF ( .NOT. humidity ) THEN
1361	IF ( ocean ) THEN
1362	CALL diffusion_e( prho, prho_reference )
1363	ELSE
1364	CALL diffusion_e( pt, pt_reference )
1365	ENDIF
1366	ELSE
1367	CALL diffusion_e( vpt, pt_reference )
1368	ENDIF
1369
1370	CALL production_e
1371
1372	!
1373	!-- Additional sink term for flows through plant canopies
1374	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
1375	CALL user_actions( 'e-tendency' )
1376
1377	!
1378	!-- Prognostic equation for TKE.
1379	!-- Eliminate negative TKE values, which can occur due to numerical
1380	!-- reasons in the course of the integration. In such cases the old TKE
1381	!-- value is reduced by 90%.
1382	DO i = nxl, nxr
1383	DO j = nys, nyn
1384	DO k = nzb_s_inner(j,i)+1, nzt
1385	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1386	tsc(3) * te_m(k,j,i) )
1387	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1388	ENDDO
1389	ENDDO
1390	ENDDO
1391
1392	!
1393	!-- Calculate tendencies for the next Runge-Kutta step
1394	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1395	IF ( intermediate_timestep_count == 1 ) THEN
1396	DO i = nxl, nxr
1397	DO j = nys, nyn
1398	DO k = nzb_s_inner(j,i)+1, nzt
1399	te_m(k,j,i) = tend(k,j,i)
1400	ENDDO
1401	ENDDO
1402	ENDDO
1403	ELSEIF ( intermediate_timestep_count < &
1404	intermediate_timestep_count_max ) THEN
1405	DO i = nxl, nxr
1406	DO j = nys, nyn
1407	DO k = nzb_s_inner(j,i)+1, nzt
1408	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1409	ENDDO
1410	ENDDO
1411	ENDDO
1412	ENDIF
1413	ENDIF
1414
1415	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1416
1417	ENDIF
1418
1419
1420	END SUBROUTINE prognostic_equations_vector
1421
1422
1423	SUBROUTINE prognostic_equations_acc
1424
1425	!------------------------------------------------------------------------------!
1426	! Version for accelerator boards
1427	!------------------------------------------------------------------------------!
1428
1429	IMPLICIT NONE
1430
1431	INTEGER :: i, j, k, runge_step
1432	REAL :: sbt
1433
1434	!
1435	!-- Set switch for intermediate Runge-Kutta step
1436	runge_step = 0
1437	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1438	IF ( intermediate_timestep_count == 1 ) THEN
1439	runge_step = 1
1440	ELSEIF ( intermediate_timestep_count < &
1441	intermediate_timestep_count_max ) THEN
1442	runge_step = 2
1443	ENDIF
1444	ENDIF
1445
1446	!
1447	!-- u-velocity component
1448	!++ Statistics still not ported to accelerators
1449	!$acc update device( hom )
1450	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1451
1452	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1453	IF ( ws_scheme_mom ) THEN
1454	CALL advec_u_ws_acc
1455	ELSE
1456	tend = 0.0 ! to be removed later??
1457	CALL advec_u_pw
1458	ENDIF
1459	ELSE
1460	CALL advec_u_up
1461	ENDIF
1462	CALL diffusion_u_acc
1463	CALL coriolis_acc( 1 )
1464	IF ( sloping_surface .AND. .NOT. neutral ) THEN
1465	CALL buoyancy( pt, pt_reference, 1, 4 )
1466	ENDIF
1467
1468	!
1469	!-- Drag by plant canopy
1470	IF ( plant_canopy ) CALL plant_canopy_model( 1 )
1471
1472	!
1473	!-- External pressure gradient
1474	IF ( dp_external ) THEN
1475	DO i = i_left, i_right
1476	DO j = j_south, j_north
1477	DO k = dp_level_ind_b+1, nzt
1478	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1479	ENDDO
1480	ENDDO
1481	ENDDO
1482	ENDIF
1483
1484	CALL user_actions( 'u-tendency' )
1485
1486	!
1487	!-- Prognostic equation for u-velocity component
1488	!$acc kernels present( nzb_u_inner, rdf, tend, tu_m, u, ug, u_p )
1489	!$acc loop
1490	DO i = i_left, i_right
1491	DO j = j_south, j_north
1492	!$acc loop vector( 32 )
1493	DO k = 1, nzt
1494	IF ( k > nzb_u_inner(j,i) ) THEN
1495	u_p(k,j,i) = u(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
1496	tsc(3) * tu_m(k,j,i) ) &
1497	- tsc(5) * rdf(k) * ( u(k,j,i) - ug(k) )
1498	!
1499	!-- Tendencies for the next Runge-Kutta step
1500	IF ( runge_step == 1 ) THEN
1501	tu_m(k,j,i) = tend(k,j,i)
1502	ELSEIF ( runge_step == 2 ) THEN
1503	tu_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tu_m(k,j,i)
1504	ENDIF
1505	ENDIF
1506	ENDDO
1507	ENDDO
1508	ENDDO
1509	!$acc end kernels
1510
1511	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1512
1513	!
1514	!-- v-velocity component
1515	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1516
1517	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1518	IF ( ws_scheme_mom ) THEN
1519	CALL advec_v_ws_acc
1520	ELSE
1521	tend = 0.0 ! to be removed later??
1522	CALL advec_v_pw
1523	END IF
1524	ELSE
1525	CALL advec_v_up
1526	ENDIF
1527	CALL diffusion_v_acc
1528	CALL coriolis_acc( 2 )
1529
1530	!
1531	!-- Drag by plant canopy
1532	IF ( plant_canopy ) CALL plant_canopy_model( 2 )
1533
1534	!
1535	!-- External pressure gradient
1536	IF ( dp_external ) THEN
1537	DO i = i_left, i_right
1538	DO j = j_south, j_north
1539	DO k = dp_level_ind_b+1, nzt
1540	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1541	ENDDO
1542	ENDDO
1543	ENDDO
1544	ENDIF
1545
1546	CALL user_actions( 'v-tendency' )
1547
1548	!
1549	!-- Prognostic equation for v-velocity component
1550	!$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
1551	!$acc loop
1552	DO i = i_left, i_right
1553	DO j = j_south, j_north
1554	!$acc loop vector( 32 )
1555	DO k = 1, nzt
1556	IF ( k > nzb_v_inner(j,i) ) THEN
1557	v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
1558	tsc(3) * tv_m(k,j,i) ) &
1559	- tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1560	!
1561	!-- Tendencies for the next Runge-Kutta step
1562	IF ( runge_step == 1 ) THEN
1563	tv_m(k,j,i) = tend(k,j,i)
1564	ELSEIF ( runge_step == 2 ) THEN
1565	tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1566	ENDIF
1567	ENDIF
1568	ENDDO
1569	ENDDO
1570	ENDDO
1571	!$acc end kernels
1572
1573	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1574
1575	!
1576	!-- w-velocity component
1577	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1578
1579	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1580	IF ( ws_scheme_mom ) THEN
1581	CALL advec_w_ws_acc
1582	ELSE
1583	tend = 0.0 ! to be removed later??
1584	CALL advec_w_pw
1585	ENDIF
1586	ELSE
1587	CALL advec_w_up
1588	ENDIF
1589	CALL diffusion_w_acc
1590	CALL coriolis_acc( 3 )
1591
1592	IF ( .NOT. neutral ) THEN
1593	IF ( ocean ) THEN
1594	CALL buoyancy( rho, rho_reference, 3, 64 )
1595	ELSE
1596	IF ( .NOT. humidity ) THEN
1597	CALL buoyancy_acc( pt, pt_reference, 3, 4 )
1598	ELSE
1599	CALL buoyancy( vpt, pt_reference, 3, 44 )
1600	ENDIF
1601	ENDIF
1602	ENDIF
1603
1604	!
1605	!-- Drag by plant canopy
1606	IF ( plant_canopy ) CALL plant_canopy_model( 3 )
1607
1608	CALL user_actions( 'w-tendency' )
1609
1610	!
1611	!-- Prognostic equation for w-velocity component
1612	!$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
1613	!$acc loop
1614	DO i = i_left, i_right
1615	DO j = j_south, j_north
1616	!$acc loop vector( 32 )
1617	DO k = 1, nzt-1
1618	IF ( k > nzb_w_inner(j,i) ) THEN
1619	w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) + &
1620	tsc(3) * tw_m(k,j,i) ) &
1621	- tsc(5) * rdf(k) * w(k,j,i)
1622	!
1623	!-- Tendencies for the next Runge-Kutta step
1624	IF ( runge_step == 1 ) THEN
1625	tw_m(k,j,i) = tend(k,j,i)
1626	ELSEIF ( runge_step == 2 ) THEN
1627	tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1628	ENDIF
1629	ENDIF
1630	ENDDO
1631	ENDDO
1632	ENDDO
1633	!$acc end kernels
1634
1635	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1636
1637
1638	!
1639	!-- If required, compute prognostic equation for potential temperature
1640	IF ( .NOT. neutral ) THEN
1641
1642	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1643
1644	!
1645	!-- pt-tendency terms with communication
1646	sbt = tsc(2)
1647	IF ( scalar_advec == 'bc-scheme' ) THEN
1648
1649	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1650	!
1651	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1652	sbt = 1.0
1653	ENDIF
1654	tend = 0.0
1655	CALL advec_s_bc( pt, 'pt' )
1656
1657	ENDIF
1658
1659	!
1660	!-- pt-tendency terms with no communication
1661	IF ( scalar_advec /= 'bc-scheme' ) THEN
1662	tend = 0.0
1663	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1664	IF ( ws_scheme_sca ) THEN
1665	CALL advec_s_ws_acc( pt, 'pt' )
1666	ELSE
1667	tend = 0.0 ! to be removed later??
1668	CALL advec_s_pw( pt )
1669	ENDIF
1670	ELSE
1671	CALL advec_s_up( pt )
1672	ENDIF
1673	ENDIF
1674
1675	CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
1676
1677	!
1678	!-- If required compute heating/cooling due to long wave radiation processes
1679	IF ( radiation ) THEN
1680	CALL calc_radiation
1681	ENDIF
1682
1683	!
1684	!-- If required compute impact of latent heat due to precipitation
1685	IF ( precipitation ) THEN
1686	CALL impact_of_latent_heat
1687	ENDIF
1688
1689	!
1690	!-- Consideration of heat sources within the plant canopy
1691	IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1692	CALL plant_canopy_model( 4 )
1693	ENDIF
1694
1695	!
1696	!-- If required compute influence of large-scale subsidence/ascent
1697	IF ( large_scale_subsidence ) THEN
1698	CALL subsidence( tend, pt, pt_init )
1699	ENDIF
1700
1701	CALL user_actions( 'pt-tendency' )
1702
1703	!
1704	!-- Prognostic equation for potential temperature
1705	!$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
1706	!$acc present( tend, tpt_m, pt, pt_p )
1707	!$acc loop
1708	DO i = i_left, i_right
1709	DO j = j_south, j_north
1710	!$acc loop vector( 32 )
1711	DO k = 1, nzt
1712	IF ( k > nzb_s_inner(j,i) ) THEN
1713	pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1714	tsc(3) * tpt_m(k,j,i) ) &
1715	- tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
1716	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
1717	!
1718	!-- Tendencies for the next Runge-Kutta step
1719	IF ( runge_step == 1 ) THEN
1720	tpt_m(k,j,i) = tend(k,j,i)
1721	ELSEIF ( runge_step == 2 ) THEN
1722	tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1723	ENDIF
1724	ENDIF
1725	ENDDO
1726	ENDDO
1727	ENDDO
1728	!$acc end kernels
1729
1730	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1731
1732	ENDIF
1733
1734	!
1735	!-- If required, compute prognostic equation for salinity
1736	IF ( ocean ) THEN
1737
1738	CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1739
1740	!
1741	!-- sa-tendency terms with communication
1742	sbt = tsc(2)
1743	IF ( scalar_advec == 'bc-scheme' ) THEN
1744
1745	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1746	!
1747	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1748	sbt = 1.0
1749	ENDIF
1750	tend = 0.0
1751	CALL advec_s_bc( sa, 'sa' )
1752
1753	ENDIF
1754
1755	!
1756	!-- sa-tendency terms with no communication
1757	IF ( scalar_advec /= 'bc-scheme' ) THEN
1758	tend = 0.0
1759	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1760	IF ( ws_scheme_sca ) THEN
1761	CALL advec_s_ws( sa, 'sa' )
1762	ELSE
1763	CALL advec_s_pw( sa )
1764	ENDIF
1765	ELSE
1766	CALL advec_s_up( sa )
1767	ENDIF
1768	ENDIF
1769
1770	CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
1771
1772	CALL user_actions( 'sa-tendency' )
1773
1774	!
1775	!-- Prognostic equation for salinity
1776	DO i = i_left, i_right
1777	DO j = j_south, j_north
1778	DO k = nzb_s_inner(j,i)+1, nzt
1779	sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1780	tsc(3) * tsa_m(k,j,i) ) &
1781	- tsc(5) * rdf_sc(k) * &
1782	( sa(k,j,i) - sa_init(k) )
1783	IF ( sa_p(k,j,i) < 0.0 ) sa_p(k,j,i) = 0.1 * sa(k,j,i)
1784	!
1785	!-- Tendencies for the next Runge-Kutta step
1786	IF ( runge_step == 1 ) THEN
1787	tsa_m(k,j,i) = tend(k,j,i)
1788	ELSEIF ( runge_step == 2 ) THEN
1789	tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i)
1790	ENDIF
1791	ENDDO
1792	ENDDO
1793	ENDDO
1794
1795	CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1796
1797	!
1798	!-- Calculate density by the equation of state for seawater
1799	CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1800	CALL eqn_state_seawater
1801	CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1802
1803	ENDIF
1804
1805	!
1806	!-- If required, compute prognostic equation for total water content / scalar
1807	IF ( humidity .OR. passive_scalar ) THEN
1808
1809	CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1810
1811	!
1812	!-- Scalar/q-tendency terms with communication
1813	sbt = tsc(2)
1814	IF ( scalar_advec == 'bc-scheme' ) THEN
1815
1816	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1817	!
1818	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1819	sbt = 1.0
1820	ENDIF
1821	tend = 0.0
1822	CALL advec_s_bc( q, 'q' )
1823
1824	ENDIF
1825
1826	!
1827	!-- Scalar/q-tendency terms with no communication
1828	IF ( scalar_advec /= 'bc-scheme' ) THEN
1829	tend = 0.0
1830	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1831	IF ( ws_scheme_sca ) THEN
1832	CALL advec_s_ws( q, 'q' )
1833	ELSE
1834	CALL advec_s_pw( q )
1835	ENDIF
1836	ELSE
1837	CALL advec_s_up( q )
1838	ENDIF
1839	ENDIF
1840
1841	CALL diffusion_s( q, qsws, qswst, wall_qflux )
1842
1843	!
1844	!-- If required compute decrease of total water content due to
1845	!-- precipitation
1846	IF ( precipitation ) THEN
1847	CALL calc_precipitation
1848	ENDIF
1849
1850	!
1851	!-- Sink or source of scalar concentration due to canopy elements
1852	IF ( plant_canopy ) CALL plant_canopy_model( 5 )
1853
1854	!
1855	!-- If required compute influence of large-scale subsidence/ascent
1856	IF ( large_scale_subsidence ) THEN
1857	CALL subsidence( tend, q, q_init )
1858	ENDIF
1859
1860	CALL user_actions( 'q-tendency' )
1861
1862	!
1863	!-- Prognostic equation for total water content / scalar
1864	DO i = i_left, i_right
1865	DO j = j_south, j_north
1866	DO k = nzb_s_inner(j,i)+1, nzt
1867	q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1868	tsc(3) * tq_m(k,j,i) ) &
1869	- tsc(5) * rdf_sc(k) * &
1870	( q(k,j,i) - q_init(k) )
1871	IF ( q_p(k,j,i) < 0.0 ) q_p(k,j,i) = 0.1 * q(k,j,i)
1872	!
1873	!-- Tendencies for the next Runge-Kutta step
1874	IF ( runge_step == 1 ) THEN
1875	tq_m(k,j,i) = tend(k,j,i)
1876	ELSEIF ( runge_step == 2 ) THEN
1877	tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1878	ENDIF
1879	ENDDO
1880	ENDDO
1881	ENDDO
1882
1883	CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1884
1885	ENDIF
1886
1887	!
1888	!-- If required, compute prognostic equation for turbulent kinetic
1889	!-- energy (TKE)
1890	IF ( .NOT. constant_diffusion ) THEN
1891
1892	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1893
1894	sbt = tsc(2)
1895	IF ( .NOT. use_upstream_for_tke ) THEN
1896	IF ( scalar_advec == 'bc-scheme' ) THEN
1897
1898	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1899	!
1900	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1901	sbt = 1.0
1902	ENDIF
1903	tend = 0.0
1904	CALL advec_s_bc( e, 'e' )
1905
1906	ENDIF
1907	ENDIF
1908
1909	!
1910	!-- TKE-tendency terms with no communication
1911	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
1912	IF ( use_upstream_for_tke ) THEN
1913	tend = 0.0
1914	CALL advec_s_up( e )
1915	ELSE
1916	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1917	IF ( ws_scheme_sca ) THEN
1918	CALL advec_s_ws_acc( e, 'e' )
1919	ELSE
1920	tend = 0.0 ! to be removed later??
1921	CALL advec_s_pw( e )
1922	ENDIF
1923	ELSE
1924	tend = 0.0 ! to be removed later??
1925	CALL advec_s_up( e )
1926	ENDIF
1927	ENDIF
1928	ENDIF
1929
1930	IF ( .NOT. humidity ) THEN
1931	IF ( ocean ) THEN
1932	CALL diffusion_e( prho, prho_reference )
1933	ELSE
1934	CALL diffusion_e_acc( pt, pt_reference )
1935	ENDIF
1936	ELSE
1937	CALL diffusion_e( vpt, pt_reference )
1938	ENDIF
1939
1940	CALL production_e_acc
1941
1942	!
1943	!-- Additional sink term for flows through plant canopies
1944	IF ( plant_canopy ) CALL plant_canopy_model( 6 )
1945	CALL user_actions( 'e-tendency' )
1946
1947	!
1948	!-- Prognostic equation for TKE.
1949	!-- Eliminate negative TKE values, which can occur due to numerical
1950	!-- reasons in the course of the integration. In such cases the old TKE
1951	!-- value is reduced by 90%.
1952	!$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
1953	!$acc loop
1954	DO i = i_left, i_right
1955	DO j = j_south, j_north
1956	!$acc loop vector( 32 )
1957	DO k = 1, nzt
1958	IF ( k > nzb_s_inner(j,i) ) THEN
1959	e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + &
1960	tsc(3) * te_m(k,j,i) )
1961	IF ( e_p(k,j,i) < 0.0 ) e_p(k,j,i) = 0.1 * e(k,j,i)
1962	!
1963	!-- Tendencies for the next Runge-Kutta step
1964	IF ( runge_step == 1 ) THEN
1965	te_m(k,j,i) = tend(k,j,i)
1966	ELSEIF ( runge_step == 2 ) THEN
1967	te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
1968	ENDIF
1969	ENDIF
1970	ENDDO
1971	ENDDO
1972	ENDDO
1973	!$acc end kernels
1974
1975	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
1976
1977	ENDIF
1978
1979
1980	END SUBROUTINE prognostic_equations_acc
1981
1982
1983	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |