Home

Context Navigation

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

Last change on this file since 1986 was 1977, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 91.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |