Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |