Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |