Home

Context Navigation

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

Last change on this file since 1330 was 1330, checked in by suehring, 10 years ago
In case of SGS-particle velocity advection of TKE is also allowed with 5th-order scheme.
Property svn:keywords set to `Id`
File size: 69.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |