Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |