Home

Context Navigation

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

Last change on this file since 1324 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 69.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |