Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |