Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |