Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |