Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |