Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |