Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |