Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |