Home

Context Navigation

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

Last change on this file since 1074 was 1054, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 67.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |