Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |