Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |