Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |