Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |