Home

Context Navigation

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

Last change on this file since 1466 was 1410, checked in by suehring, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 83.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |