Home

Context Navigation

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

Last change on this file since 1213 was 1182, checked in by raasch, 11 years ago
last commit documented, rc-file for example run updated
Property svn:keywords set to `Id`
File size: 65.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |