Home

Context Navigation

prognostic_equations_mod.f90 @ 1851

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

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

Context Navigation

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

Download in other formats: