Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |