Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |