Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |