Home

Context Navigation

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

Last change on this file since 1398 was 1398, checked in by heinze, 10 years ago
adjustments in case of nudging + bugfix for KIND in CMPLX function
Property svn:keywords set to `Id`
File size: 83.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |