Home

Context Navigation

source: palm/trunk/SOURCE/prognostic_equations_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 83.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |