Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |