Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |