Home

Context Navigation

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

Last change on this file since 1493 was 1485, checked in by kanani, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 83.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |