Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |