Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |