Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |