Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |