Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |