Home

Context Navigation

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

Last change on this file since 1365 was 1365, checked in by boeske, 10 years ago
large scale forcing enabled
Property svn:keywords set to `Id`
File size: 82.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |