Home

Context Navigation

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

Last change on this file since 4343 was 4329, checked in by motisi, 5 years ago
Renamed wall_flags_0 to wall_flags_static_0
Property svn:keywords set to `Id`
File size: 59.2 KB

Line
1	!> @file prognostic_equations.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: prognostic_equations.f90 4329 2019-12-10 15:46:36Z oliver.maas $
27	! Renamed wall_flags_0 to wall_flags_static_0
28	!
29	! 4182 2019-08-22 15:20:23Z scharf
30	! Corrected "Former revisions" section
31	!
32	! 4110 2019-07-22 17:05:21Z suehring
33	! pass integer flag array to WS scalar advection routine which is now necessary
34	! as the flags may differ for scalars, e.g. pt can be cyclic while chemical
35	! species may be non-cyclic. Further, pass boundary flags.
36	!
37	! 4109 2019-07-22 17:00:34Z suehring
38	! Application of monotonic flux limiter for the vertical scalar advection
39	! up to the topography top (only for the cache-optimized version at the
40	! moment). Please note, at the moment the limiter is only applied for passive
41	! scalars.
42	!
43	! 4048 2019-06-21 21:00:21Z knoop
44	! Moved tcm_prognostic_equations to module_interface
45	!
46	! 3987 2019-05-22 09:52:13Z kanani
47	! Introduce alternative switch for debug output during timestepping
48	!
49	! 3956 2019-05-07 12:32:52Z monakurppa
50	! Removed salsa calls.
51	!
52	! 3931 2019-04-24 16:34:28Z schwenkel
53	! Correct/complete module_interface introduction for chemistry model
54	!
55	! 3899 2019-04-16 14:05:27Z monakurppa
56	! Corrections in the OpenMP version of salsa
57	!
58	! 3887 2019 -04-12 08:47:41Z schwenkel
59	! Implicit Bugfix for chemistry model, loop for non_transport_physics over
60	! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
61	!
62	! 3885 2019-04-11 11:29:34Z kanani
63	! Changes related to global restructuring of location messages and introduction
64	! of additional debug messages
65	!
66	! 3881 2019-04-10 09:31:22Z suehring
67	! Bugfix in OpenMP directive
68	!
69	! 3880 2019-04-08 21:43:02Z knoop
70	! Moved wtm_tendencies to module_interface_actions
71	!
72	! 3874 2019-04-08 16:53:48Z knoop
73	! Added non_transport_physics module interfaces and moved bcm code into it
74	!
75	! 3872 2019-04-08 15:03:06Z knoop
76	! Moving prognostic equations of bcm into bulk_cloud_model_mod
77	!
78	! 3864 2019-04-05 09:01:56Z monakurppa
79	! Modifications made for salsa:
80	! - salsa_prognostic_equations moved to salsa_mod (and the call to
81	! module_interface_mod)
82	! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
83	! ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
84	!
85	! 3840 2019-03-29 10:35:52Z knoop
86	! added USE chem_gasphase_mod for nspec, nspec and spc_names
87	!
88	! 3820 2019-03-27 11:53:41Z forkel
89	! renamed do_depo to deposition_dry (ecc)
90	!
91	! 3797 2019-03-15 11:15:38Z forkel
92	! Call chem_integegrate in OpenMP loop (ketelsen)
93	!
94	!
95	! 3771 2019-02-28 12:19:33Z raasch
96	! preprocessor directivs fro rrtmg added
97	!
98	! 3761 2019-02-25 15:31:42Z raasch
99	! unused variable removed
100	!
101	! 3719 2019-02-06 13:10:18Z kanani
102	! Cleaned up chemistry cpu measurements
103	!
104	! 3684 2019-01-20 20:20:58Z knoop
105	! OpenACC port for SPEC
106	!
107	! Revision 1.1 2000/04/13 14:56:27 schroeter
108	! Initial revision
109	!
110	!
111	! Description:
112	! ------------
113	!> Solving the prognostic equations.
114	!------------------------------------------------------------------------------!
115	MODULE prognostic_equations_mod
116
117	USE advec_s_bc_mod, &
118	ONLY: advec_s_bc
119
120	USE advec_s_pw_mod, &
121	ONLY: advec_s_pw
122
123	USE advec_s_up_mod, &
124	ONLY: advec_s_up
125
126	USE advec_u_pw_mod, &
127	ONLY: advec_u_pw
128
129	USE advec_u_up_mod, &
130	ONLY: advec_u_up
131
132	USE advec_v_pw_mod, &
133	ONLY: advec_v_pw
134
135	USE advec_v_up_mod, &
136	ONLY: advec_v_up
137
138	USE advec_w_pw_mod, &
139	ONLY: advec_w_pw
140
141	USE advec_w_up_mod, &
142	ONLY: advec_w_up
143
144	USE advec_ws, &
145	ONLY: advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
146
147	USE arrays_3d, &
148	ONLY: diss_l_e, diss_l_pt, diss_l_q, &
149	diss_l_s, diss_l_sa, diss_s_e, &
150	diss_s_pt, diss_s_q, diss_s_s, diss_s_sa, &
151	e, e_p, flux_s_e, flux_s_pt, flux_s_q, &
152	flux_s_s, flux_s_sa, flux_l_e, &
153	flux_l_pt, flux_l_q, flux_l_s, &
154	flux_l_sa, pt, ptdf_x, ptdf_y, pt_init, &
155	pt_p, prho, q, q_init, q_p, rdf, rdf_sc, &
156	ref_state, rho_ocean, s, s_init, s_p, tend, te_m, &
157	tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u, &
158	ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
159
160	USE buoyancy_mod, &
161	ONLY: buoyancy
162
163	USE control_parameters, &
164	ONLY: bc_dirichlet_l, &
165	bc_dirichlet_n, &
166	bc_dirichlet_r, &
167	bc_dirichlet_s, &
168	bc_radiation_l, &
169	bc_radiation_n, &
170	bc_radiation_r, &
171	bc_radiation_s, &
172	constant_diffusion, &
173	debug_output_timestep, &
174	dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, &
175	humidity, intermediate_timestep_count, &
176	intermediate_timestep_count_max, large_scale_forcing, &
177	large_scale_subsidence, &
178	monotonic_limiter_z, &
179	neutral, nudging, &
180	ocean_mode, passive_scalar, plant_canopy, pt_reference, &
181	scalar_advec, scalar_advec, simulated_time, sloping_surface, &
182	timestep_scheme, tsc, use_subsidence_tendencies, &
183	use_upstream_for_tke, wind_turbine, ws_scheme_mom, &
184	ws_scheme_sca, urban_surface, land_surface, &
185	time_since_reference_point, salsa
186
187	USE coriolis_mod, &
188	ONLY: coriolis
189
190	USE cpulog, &
191	ONLY: cpu_log, log_point, log_point_s
192
193	USE diffusion_s_mod, &
194	ONLY: diffusion_s
195
196	USE diffusion_u_mod, &
197	ONLY: diffusion_u
198
199	USE diffusion_v_mod, &
200	ONLY: diffusion_v
201
202	USE diffusion_w_mod, &
203	ONLY: diffusion_w
204
205	USE indices, &
206	ONLY: advc_flags_s, &
207	nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
208	nzb, nzt, wall_flags_static_0
209
210	USE kinds
211
212	USE lsf_nudging_mod, &
213	ONLY: ls_advec, nudge
214
215	USE module_interface, &
216	ONLY: module_interface_actions, &
217	module_interface_non_advective_processes, &
218	module_interface_exchange_horiz, &
219	module_interface_prognostic_equations
220
221	USE ocean_mod, &
222	ONLY: stokes_drift_terms, stokes_force, &
223	wave_breaking, wave_breaking_term
224
225	USE plant_canopy_model_mod, &
226	ONLY: cthf, pcm_tendency
227
228	#if defined( __rrtmg )
229	USE radiation_model_mod, &
230	ONLY: radiation, radiation_tendency, &
231	skip_time_do_radiation
232	#endif
233
234	USE statistics, &
235	ONLY: hom
236
237	USE subsidence_mod, &
238	ONLY: subsidence
239
240	USE surface_mod, &
241	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
242	surf_usm_v
243
244	IMPLICIT NONE
245
246	PRIVATE
247	PUBLIC prognostic_equations_cache, prognostic_equations_vector
248
249	INTERFACE prognostic_equations_cache
250	MODULE PROCEDURE prognostic_equations_cache
251	END INTERFACE prognostic_equations_cache
252
253	INTERFACE prognostic_equations_vector
254	MODULE PROCEDURE prognostic_equations_vector
255	END INTERFACE prognostic_equations_vector
256
257
258	CONTAINS
259
260
261	!------------------------------------------------------------------------------!
262	! Description:
263	! ------------
264	!> Version with one optimized loop over all equations. It is only allowed to
265	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
266	!>
267	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
268	!> so communication between CPUs is not allowed (does not make sense) within
269	!> these loops.
270	!>
271	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
272	!------------------------------------------------------------------------------!
273
274	SUBROUTINE prognostic_equations_cache
275
276
277	INTEGER(iwp) :: i !<
278	INTEGER(iwp) :: i_omp_start !<
279	INTEGER(iwp) :: j !<
280	INTEGER(iwp) :: k !<
281	!$ INTEGER(iwp) :: omp_get_thread_num !<
282	INTEGER(iwp) :: tn = 0 !<
283
284	LOGICAL :: loop_start !<
285
286
287	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'start' )
288	!
289	!-- Time measurement can only be performed for the whole set of equations
290	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
291
292	!$OMP PARALLEL PRIVATE (i,j)
293	!$OMP DO
294	DO i = nxl, nxr
295	DO j = nys, nyn
296	!
297	!-- Calculate non advective processes for all other modules
298	CALL module_interface_non_advective_processes( i, j )
299	ENDDO
300	ENDDO
301	!
302	!-- Module Inferface for exchange horiz after non_advective_processes but before
303	!-- advection. Therefore, non_advective_processes must not run for ghost points.
304	!$OMP END PARALLEL
305	CALL module_interface_exchange_horiz()
306	!
307	!-- Loop over all prognostic equations
308	!$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
309
310	!$ tn = omp_get_thread_num()
311	loop_start = .TRUE.
312
313	!$OMP DO
314	DO i = nxl, nxr
315
316	!
317	!-- Store the first loop index. It differs for each thread and is required
318	!-- later in advec_ws
319	IF ( loop_start ) THEN
320	loop_start = .FALSE.
321	i_omp_start = i
322	ENDIF
323
324	DO j = nys, nyn
325	!
326	!-- Tendency terms for u-velocity component. Please note, in case of
327	!-- non-cyclic boundary conditions the grid point i=0 is excluded from
328	!-- the prognostic equations for the u-component.
329	IF ( i >= nxlu ) THEN
330
331	tend(:,j,i) = 0.0_wp
332	IF ( timestep_scheme(1:5) == 'runge' ) THEN
333	IF ( ws_scheme_mom ) THEN
334	CALL advec_u_ws( i, j, i_omp_start, tn )
335	ELSE
336	CALL advec_u_pw( i, j )
337	ENDIF
338	ELSE
339	CALL advec_u_up( i, j )
340	ENDIF
341	CALL diffusion_u( i, j )
342	CALL coriolis( i, j, 1 )
343	IF ( sloping_surface .AND. .NOT. neutral ) THEN
344	CALL buoyancy( i, j, pt, 1 )
345	ENDIF
346
347	!
348	!-- Drag by plant canopy
349	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
350
351	!
352	!-- External pressure gradient
353	IF ( dp_external ) THEN
354	DO k = dp_level_ind_b+1, nzt
355	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
356	ENDDO
357	ENDIF
358
359	!
360	!-- Nudging
361	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
362
363	!
364	!-- Effect of Stokes drift (in ocean mode only)
365	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 )
366
367	CALL module_interface_actions( i, j, 'u-tendency' )
368	!
369	!-- Prognostic equation for u-velocity component
370	DO k = nzb+1, nzt
371
372	u_p(k,j,i) = u(k,j,i) + ( dt_3d * &
373	( tsc(2) * tend(k,j,i) + &
374	tsc(3) * tu_m(k,j,i) ) &
375	- tsc(5) * rdf(k) &
376	* ( u(k,j,i) - u_init(k) ) &
377	) * MERGE( 1.0_wp, 0.0_wp, &
378	BTEST( wall_flags_static_0(k,j,i), 1 )&
379	)
380	ENDDO
381
382	!
383	!-- Add turbulence generated by wave breaking (in ocean mode only)
384	IF ( wave_breaking .AND. &
385	intermediate_timestep_count == intermediate_timestep_count_max )&
386	THEN
387	CALL wave_breaking_term( i, j, 1 )
388	ENDIF
389
390	!
391	!-- Calculate tendencies for the next Runge-Kutta step
392	IF ( timestep_scheme(1:5) == 'runge' ) THEN
393	IF ( intermediate_timestep_count == 1 ) THEN
394	DO k = nzb+1, nzt
395	tu_m(k,j,i) = tend(k,j,i)
396	ENDDO
397	ELSEIF ( intermediate_timestep_count < &
398	intermediate_timestep_count_max ) THEN
399	DO k = nzb+1, nzt
400	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
401	+ 5.3125_wp * tu_m(k,j,i)
402	ENDDO
403	ENDIF
404	ENDIF
405
406	ENDIF
407	!
408	!-- Tendency terms for v-velocity component. Please note, in case of
409	!-- non-cyclic boundary conditions the grid point j=0 is excluded from
410	!-- the prognostic equations for the v-component. !--
411	IF ( j >= nysv ) THEN
412
413	tend(:,j,i) = 0.0_wp
414	IF ( timestep_scheme(1:5) == 'runge' ) THEN
415	IF ( ws_scheme_mom ) THEN
416	CALL advec_v_ws( i, j, i_omp_start, tn )
417	ELSE
418	CALL advec_v_pw( i, j )
419	ENDIF
420	ELSE
421	CALL advec_v_up( i, j )
422	ENDIF
423	CALL diffusion_v( i, j )
424	CALL coriolis( i, j, 2 )
425
426	!
427	!-- Drag by plant canopy
428	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
429
430	!
431	!-- External pressure gradient
432	IF ( dp_external ) THEN
433	DO k = dp_level_ind_b+1, nzt
434	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
435	ENDDO
436	ENDIF
437
438	!
439	!-- Nudging
440	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
441
442	!
443	!-- Effect of Stokes drift (in ocean mode only)
444	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 )
445
446	CALL module_interface_actions( i, j, 'v-tendency' )
447	!
448	!-- Prognostic equation for v-velocity component
449	DO k = nzb+1, nzt
450	v_p(k,j,i) = v(k,j,i) + ( dt_3d * &
451	( tsc(2) * tend(k,j,i) + &
452	tsc(3) * tv_m(k,j,i) ) &
453	- tsc(5) * rdf(k) &
454	* ( v(k,j,i) - v_init(k) )&
455	) * MERGE( 1.0_wp, 0.0_wp, &
456	BTEST( wall_flags_static_0(k,j,i), 2 )&
457	)
458	ENDDO
459
460	!
461	!-- Add turbulence generated by wave breaking (in ocean mode only)
462	IF ( wave_breaking .AND. &
463	intermediate_timestep_count == intermediate_timestep_count_max )&
464	THEN
465	CALL wave_breaking_term( i, j, 2 )
466	ENDIF
467
468	!
469	!-- Calculate tendencies for the next Runge-Kutta step
470	IF ( timestep_scheme(1:5) == 'runge' ) THEN
471	IF ( intermediate_timestep_count == 1 ) THEN
472	DO k = nzb+1, nzt
473	tv_m(k,j,i) = tend(k,j,i)
474	ENDDO
475	ELSEIF ( intermediate_timestep_count < &
476	intermediate_timestep_count_max ) THEN
477	DO k = nzb+1, nzt
478	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
479	+ 5.3125_wp * tv_m(k,j,i)
480	ENDDO
481	ENDIF
482	ENDIF
483
484	ENDIF
485
486	!
487	!-- Tendency terms for w-velocity component
488	tend(:,j,i) = 0.0_wp
489	IF ( timestep_scheme(1:5) == 'runge' ) THEN
490	IF ( ws_scheme_mom ) THEN
491	CALL advec_w_ws( i, j, i_omp_start, tn )
492	ELSE
493	CALL advec_w_pw( i, j )
494	END IF
495	ELSE
496	CALL advec_w_up( i, j )
497	ENDIF
498	CALL diffusion_w( i, j )
499	CALL coriolis( i, j, 3 )
500
501	IF ( .NOT. neutral ) THEN
502	IF ( ocean_mode ) THEN
503	CALL buoyancy( i, j, rho_ocean, 3 )
504	ELSE
505	IF ( .NOT. humidity ) THEN
506	CALL buoyancy( i, j, pt, 3 )
507	ELSE
508	CALL buoyancy( i, j, vpt, 3 )
509	ENDIF
510	ENDIF
511	ENDIF
512
513	!
514	!-- Drag by plant canopy
515	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
516
517	!
518	!-- Effect of Stokes drift (in ocean mode only)
519	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 3 )
520
521	CALL module_interface_actions( i, j, 'w-tendency' )
522	!
523	!-- Prognostic equation for w-velocity component
524	DO k = nzb+1, nzt-1
525	w_p(k,j,i) = w(k,j,i) + ( dt_3d * &
526	( tsc(2) * tend(k,j,i) + &
527	tsc(3) * tw_m(k,j,i) ) &
528	- tsc(5) * rdf(k) * w(k,j,i) &
529	) * MERGE( 1.0_wp, 0.0_wp, &
530	BTEST( wall_flags_static_0(k,j,i), 3 )&
531	)
532	ENDDO
533
534	!
535	!-- Calculate tendencies for the next Runge-Kutta step
536	IF ( timestep_scheme(1:5) == 'runge' ) THEN
537	IF ( intermediate_timestep_count == 1 ) THEN
538	DO k = nzb+1, nzt-1
539	tw_m(k,j,i) = tend(k,j,i)
540	ENDDO
541	ELSEIF ( intermediate_timestep_count < &
542	intermediate_timestep_count_max ) THEN
543	DO k = nzb+1, nzt-1
544	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
545	+ 5.3125_wp * tw_m(k,j,i)
546	ENDDO
547	ENDIF
548	ENDIF
549
550	!
551	!-- If required, compute prognostic equation for potential temperature
552	IF ( .NOT. neutral ) THEN
553	!
554	!-- Tendency terms for potential temperature
555	tend(:,j,i) = 0.0_wp
556	IF ( timestep_scheme(1:5) == 'runge' ) THEN
557	IF ( ws_scheme_sca ) THEN
558	CALL advec_s_ws( advc_flags_s, &
559	i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
560	flux_l_pt, diss_l_pt, i_omp_start, tn, &
561	bc_dirichlet_l .OR. bc_radiation_l, &
562	bc_dirichlet_n .OR. bc_radiation_n, &
563	bc_dirichlet_r .OR. bc_radiation_r, &
564	bc_dirichlet_s .OR. bc_radiation_s )
565	ELSE
566	CALL advec_s_pw( i, j, pt )
567	ENDIF
568	ELSE
569	CALL advec_s_up( i, j, pt )
570	ENDIF
571	CALL diffusion_s( i, j, pt, &
572	surf_def_h(0)%shf, surf_def_h(1)%shf, &
573	surf_def_h(2)%shf, &
574	surf_lsm_h%shf, surf_usm_h%shf, &
575	surf_def_v(0)%shf, surf_def_v(1)%shf, &
576	surf_def_v(2)%shf, surf_def_v(3)%shf, &
577	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
578	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
579	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
580	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
581
582	!
583	!-- Consideration of heat sources within the plant canopy
584	IF ( plant_canopy .AND. &
585	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
586	CALL pcm_tendency( i, j, 4 )
587	ENDIF
588
589	!
590	!-- Large scale advection
591	IF ( large_scale_forcing ) THEN
592	CALL ls_advec( i, j, simulated_time, 'pt' )
593	ENDIF
594
595	!
596	!-- Nudging
597	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
598
599	!
600	!-- If required, compute effect of large-scale subsidence/ascent
601	IF ( large_scale_subsidence .AND. &
602	.NOT. use_subsidence_tendencies ) THEN
603	CALL subsidence( i, j, tend, pt, pt_init, 2 )
604	ENDIF
605
606	#if defined( __rrtmg )
607	!
608	!-- If required, add tendency due to radiative heating/cooling
609	IF ( radiation .AND. &
610	simulated_time > skip_time_do_radiation ) THEN
611	CALL radiation_tendency ( i, j, tend )
612	ENDIF
613	#endif
614
615	CALL module_interface_actions( i, j, 'pt-tendency' )
616	!
617	!-- Prognostic equation for potential temperature
618	DO k = nzb+1, nzt
619	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * &
620	( tsc(2) * tend(k,j,i) + &
621	tsc(3) * tpt_m(k,j,i) ) &
622	- tsc(5) &
623	* ( pt(k,j,i) - pt_init(k) ) &
624	* ( rdf_sc(k) + ptdf_x(i) &
625	+ ptdf_y(j) ) &
626	) &
627	* MERGE( 1.0_wp, 0.0_wp, &
628	BTEST( wall_flags_static_0(k,j,i), 0 )&
629	)
630	ENDDO
631
632	!
633	!-- Calculate tendencies for the next Runge-Kutta step
634	IF ( timestep_scheme(1:5) == 'runge' ) THEN
635	IF ( intermediate_timestep_count == 1 ) THEN
636	DO k = nzb+1, nzt
637	tpt_m(k,j,i) = tend(k,j,i)
638	ENDDO
639	ELSEIF ( intermediate_timestep_count < &
640	intermediate_timestep_count_max ) THEN
641	DO k = nzb+1, nzt
642	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
643	5.3125_wp * tpt_m(k,j,i)
644	ENDDO
645	ENDIF
646	ENDIF
647
648	ENDIF
649
650	!
651	!-- If required, compute prognostic equation for total water content
652	IF ( humidity ) THEN
653
654	!
655	!-- Tendency-terms for total water content / scalar
656	tend(:,j,i) = 0.0_wp
657	IF ( timestep_scheme(1:5) == 'runge' ) &
658	THEN
659	IF ( ws_scheme_sca ) THEN
660	CALL advec_s_ws( advc_flags_s, &
661	i, j, q, 'q', flux_s_q, &
662	diss_s_q, flux_l_q, diss_l_q, &
663	i_omp_start, tn, &
664	bc_dirichlet_l .OR. bc_radiation_l, &
665	bc_dirichlet_n .OR. bc_radiation_n, &
666	bc_dirichlet_r .OR. bc_radiation_r, &
667	bc_dirichlet_s .OR. bc_radiation_s )
668	ELSE
669	CALL advec_s_pw( i, j, q )
670	ENDIF
671	ELSE
672	CALL advec_s_up( i, j, q )
673	ENDIF
674	CALL diffusion_s( i, j, q, &
675	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
676	surf_def_h(2)%qsws, &
677	surf_lsm_h%qsws, surf_usm_h%qsws, &
678	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
679	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
680	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
681	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
682	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
683	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
684
685	!
686	!-- Sink or source of humidity due to canopy elements
687	IF ( plant_canopy ) CALL pcm_tendency( i, j, 5 )
688
689	!
690	!-- Large scale advection
691	IF ( large_scale_forcing ) THEN
692	CALL ls_advec( i, j, simulated_time, 'q' )
693	ENDIF
694
695	!
696	!-- Nudging
697	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
698
699	!
700	!-- If required compute influence of large-scale subsidence/ascent
701	IF ( large_scale_subsidence .AND. &
702	.NOT. use_subsidence_tendencies ) THEN
703	CALL subsidence( i, j, tend, q, q_init, 3 )
704	ENDIF
705
706	CALL module_interface_actions( i, j, 'q-tendency' )
707
708	!
709	!-- Prognostic equation for total water content / scalar
710	DO k = nzb+1, nzt
711	q_p(k,j,i) = q(k,j,i) + ( dt_3d * &
712	( tsc(2) * tend(k,j,i) + &
713	tsc(3) * tq_m(k,j,i) ) &
714	- tsc(5) * rdf_sc(k) * &
715	( q(k,j,i) - q_init(k) ) &
716	) &
717	* MERGE( 1.0_wp, 0.0_wp, &
718	BTEST( wall_flags_static_0(k,j,i), 0 )&
719	)
720	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
721	ENDDO
722
723	!
724	!-- Calculate tendencies for the next Runge-Kutta step
725	IF ( timestep_scheme(1:5) == 'runge' ) THEN
726	IF ( intermediate_timestep_count == 1 ) THEN
727	DO k = nzb+1, nzt
728	tq_m(k,j,i) = tend(k,j,i)
729	ENDDO
730	ELSEIF ( intermediate_timestep_count < &
731	intermediate_timestep_count_max ) THEN
732	DO k = nzb+1, nzt
733	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
734	5.3125_wp * tq_m(k,j,i)
735	ENDDO
736	ENDIF
737	ENDIF
738
739	ENDIF
740
741	!
742	!-- If required, compute prognostic equation for scalar
743	IF ( passive_scalar ) THEN
744	!
745	!-- Tendency-terms for total water content / scalar
746	tend(:,j,i) = 0.0_wp
747	IF ( timestep_scheme(1:5) == 'runge' ) &
748	THEN
749	IF ( ws_scheme_sca ) THEN
750	!
751	!-- For scalar advection apply monotonic flux limiter near
752	!-- topography.
753	CALL advec_s_ws( advc_flags_s, &
754	i, j, s, 's', flux_s_s, &
755	diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
756	tn, &
757	bc_dirichlet_l .OR. bc_radiation_l, &
758	bc_dirichlet_n .OR. bc_radiation_n, &
759	bc_dirichlet_r .OR. bc_radiation_r, &
760	bc_dirichlet_s .OR. bc_radiation_s, &
761	monotonic_limiter_z )
762	ELSE
763	CALL advec_s_pw( i, j, s )
764	ENDIF
765	ELSE
766	CALL advec_s_up( i, j, s )
767	ENDIF
768	CALL diffusion_s( i, j, s, &
769	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
770	surf_def_h(2)%ssws, &
771	surf_lsm_h%ssws, surf_usm_h%ssws, &
772	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
773	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
774	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
775	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
776	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
777	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
778
779	!
780	!-- Sink or source of scalar concentration due to canopy elements
781	IF ( plant_canopy ) CALL pcm_tendency( i, j, 7 )
782
783	!
784	!-- Large scale advection, still need to be extended for scalars
785	! IF ( large_scale_forcing ) THEN
786	! CALL ls_advec( i, j, simulated_time, 's' )
787	! ENDIF
788
789	!
790	!-- Nudging, still need to be extended for scalars
791	! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' )
792
793	!
794	!-- If required compute influence of large-scale subsidence/ascent.
795	!-- Note, the last argument is of no meaning in this case, as it is
796	!-- only used in conjunction with large_scale_forcing, which is to
797	!-- date not implemented for scalars.
798	IF ( large_scale_subsidence .AND. &
799	.NOT. use_subsidence_tendencies .AND. &
800	.NOT. large_scale_forcing ) THEN
801	CALL subsidence( i, j, tend, s, s_init, 3 )
802	ENDIF
803
804	CALL module_interface_actions( i, j, 's-tendency' )
805
806	!
807	!-- Prognostic equation for scalar
808	DO k = nzb+1, nzt
809	s_p(k,j,i) = s(k,j,i) + ( dt_3d * &
810	( tsc(2) * tend(k,j,i) + &
811	tsc(3) * ts_m(k,j,i) ) &
812	- tsc(5) * rdf_sc(k) &
813	* ( s(k,j,i) - s_init(k) )&
814	) &
815	* MERGE( 1.0_wp, 0.0_wp, &
816	BTEST( wall_flags_static_0(k,j,i), 0 )&
817	)
818	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
819	ENDDO
820
821	!
822	!-- Calculate tendencies for the next Runge-Kutta step
823	IF ( timestep_scheme(1:5) == 'runge' ) THEN
824	IF ( intermediate_timestep_count == 1 ) THEN
825	DO k = nzb+1, nzt
826	ts_m(k,j,i) = tend(k,j,i)
827	ENDDO
828	ELSEIF ( intermediate_timestep_count < &
829	intermediate_timestep_count_max ) THEN
830	DO k = nzb+1, nzt
831	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
832	5.3125_wp * ts_m(k,j,i)
833	ENDDO
834	ENDIF
835	ENDIF
836
837	ENDIF
838	!
839	!-- Calculate prognostic equations for all other modules
840	CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
841
842	ENDDO ! loop over j
843	ENDDO ! loop over i
844	!$OMP END PARALLEL
845
846
847	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
848
849	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_cache', 'end' )
850
851	END SUBROUTINE prognostic_equations_cache
852
853
854	!------------------------------------------------------------------------------!
855	! Description:
856	! ------------
857	!> Version for vector machines
858	!------------------------------------------------------------------------------!
859
860	SUBROUTINE prognostic_equations_vector
861
862
863	INTEGER(iwp) :: i !<
864	INTEGER(iwp) :: j !<
865	INTEGER(iwp) :: k !<
866
867	REAL(wp) :: sbt !<
868
869
870	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'start' )
871	!
872	!-- Calculate non advective processes for all other modules
873	CALL module_interface_non_advective_processes
874	!
875	!-- Module Inferface for exchange horiz after non_advective_processes but before
876	!-- advection. Therefore, non_advective_processes must not run for ghost points.
877	CALL module_interface_exchange_horiz()
878	!
879	!-- u-velocity component
880	CALL cpu_log( log_point(5), 'u-equation', 'start' )
881
882	!$ACC KERNELS PRESENT(tend)
883	tend = 0.0_wp
884	!$ACC END KERNELS
885	IF ( timestep_scheme(1:5) == 'runge' ) THEN
886	IF ( ws_scheme_mom ) THEN
887	CALL advec_u_ws
888	ELSE
889	CALL advec_u_pw
890	ENDIF
891	ELSE
892	CALL advec_u_up
893	ENDIF
894	CALL diffusion_u
895	CALL coriolis( 1 )
896	IF ( sloping_surface .AND. .NOT. neutral ) THEN
897	CALL buoyancy( pt, 1 )
898	ENDIF
899
900	!
901	!-- Drag by plant canopy
902	IF ( plant_canopy ) CALL pcm_tendency( 1 )
903
904	!
905	!-- External pressure gradient
906	IF ( dp_external ) THEN
907	DO i = nxlu, nxr
908	DO j = nys, nyn
909	DO k = dp_level_ind_b+1, nzt
910	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
911	ENDDO
912	ENDDO
913	ENDDO
914	ENDIF
915
916	!
917	!-- Nudging
918	IF ( nudging ) CALL nudge( simulated_time, 'u' )
919
920	!
921	!-- Effect of Stokes drift (in ocean mode only)
922	IF ( stokes_force ) CALL stokes_drift_terms( 1 )
923
924	CALL module_interface_actions( 'u-tendency' )
925
926	!
927	!-- Prognostic equation for u-velocity component
928	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
929	!$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_static_0) &
930	!$ACC PRESENT(tsc(2:5)) &
931	!$ACC PRESENT(u_p)
932	DO i = nxlu, nxr
933	DO j = nys, nyn
934	DO k = nzb+1, nzt
935	u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
936	tsc(3) * tu_m(k,j,i) ) &
937	- tsc(5) * rdf(k) * &
938	( u(k,j,i) - u_init(k) ) &
939	) * MERGE( 1.0_wp, 0.0_wp, &
940	BTEST( wall_flags_static_0(k,j,i), 1 ) &
941	)
942	ENDDO
943	ENDDO
944	ENDDO
945
946	!
947	!-- Add turbulence generated by wave breaking (in ocean mode only)
948	IF ( wave_breaking .AND. &
949	intermediate_timestep_count == intermediate_timestep_count_max ) &
950	THEN
951	CALL wave_breaking_term( 1 )
952	ENDIF
953
954	!
955	!-- Calculate tendencies for the next Runge-Kutta step
956	IF ( timestep_scheme(1:5) == 'runge' ) THEN
957	IF ( intermediate_timestep_count == 1 ) THEN
958	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
959	!$ACC PRESENT(tend, tu_m)
960	DO i = nxlu, nxr
961	DO j = nys, nyn
962	DO k = nzb+1, nzt
963	tu_m(k,j,i) = tend(k,j,i)
964	ENDDO
965	ENDDO
966	ENDDO
967	ELSEIF ( intermediate_timestep_count < &
968	intermediate_timestep_count_max ) THEN
969	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
970	!$ACC PRESENT(tend, tu_m)
971	DO i = nxlu, nxr
972	DO j = nys, nyn
973	DO k = nzb+1, nzt
974	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
975	+ 5.3125_wp * tu_m(k,j,i)
976	ENDDO
977	ENDDO
978	ENDDO
979	ENDIF
980	ENDIF
981
982	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
983
984	!
985	!-- v-velocity component
986	CALL cpu_log( log_point(6), 'v-equation', 'start' )
987
988	!$ACC KERNELS PRESENT(tend)
989	tend = 0.0_wp
990	!$ACC END KERNELS
991	IF ( timestep_scheme(1:5) == 'runge' ) THEN
992	IF ( ws_scheme_mom ) THEN
993	CALL advec_v_ws
994	ELSE
995	CALL advec_v_pw
996	END IF
997	ELSE
998	CALL advec_v_up
999	ENDIF
1000	CALL diffusion_v
1001	CALL coriolis( 2 )
1002
1003	!
1004	!-- Drag by plant canopy
1005	IF ( plant_canopy ) CALL pcm_tendency( 2 )
1006
1007	!
1008	!-- External pressure gradient
1009	IF ( dp_external ) THEN
1010	DO i = nxl, nxr
1011	DO j = nysv, nyn
1012	DO k = dp_level_ind_b+1, nzt
1013	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1014	ENDDO
1015	ENDDO
1016	ENDDO
1017	ENDIF
1018
1019	!
1020	!-- Nudging
1021	IF ( nudging ) CALL nudge( simulated_time, 'v' )
1022
1023	!
1024	!-- Effect of Stokes drift (in ocean mode only)
1025	IF ( stokes_force ) CALL stokes_drift_terms( 2 )
1026
1027	CALL module_interface_actions( 'v-tendency' )
1028
1029	!
1030	!-- Prognostic equation for v-velocity component
1031	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1032	!$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_static_0) &
1033	!$ACC PRESENT(tsc(2:5)) &
1034	!$ACC PRESENT(v_p)
1035	DO i = nxl, nxr
1036	DO j = nysv, nyn
1037	DO k = nzb+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) * &
1041	( v(k,j,i) - v_init(k) ) &
1042	) * MERGE( 1.0_wp, 0.0_wp, &
1043	BTEST( wall_flags_static_0(k,j,i), 2 )&
1044	)
1045	ENDDO
1046	ENDDO
1047	ENDDO
1048
1049	!
1050	!-- Add turbulence generated by wave breaking (in ocean mode only)
1051	IF ( wave_breaking .AND. &
1052	intermediate_timestep_count == intermediate_timestep_count_max ) &
1053	THEN
1054	CALL wave_breaking_term( 2 )
1055	ENDIF
1056
1057	!
1058	!-- Calculate tendencies for the next Runge-Kutta step
1059	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1060	IF ( intermediate_timestep_count == 1 ) THEN
1061	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1062	!$ACC PRESENT(tend, tv_m)
1063	DO i = nxl, nxr
1064	DO j = nysv, nyn
1065	DO k = nzb+1, nzt
1066	tv_m(k,j,i) = tend(k,j,i)
1067	ENDDO
1068	ENDDO
1069	ENDDO
1070	ELSEIF ( intermediate_timestep_count < &
1071	intermediate_timestep_count_max ) THEN
1072	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1073	!$ACC PRESENT(tend, tv_m)
1074	DO i = nxl, nxr
1075	DO j = nysv, nyn
1076	DO k = nzb+1, nzt
1077	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1078	+ 5.3125_wp * tv_m(k,j,i)
1079	ENDDO
1080	ENDDO
1081	ENDDO
1082	ENDIF
1083	ENDIF
1084
1085	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1086
1087	!
1088	!-- w-velocity component
1089	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1090
1091	!$ACC KERNELS PRESENT(tend)
1092	tend = 0.0_wp
1093	!$ACC END KERNELS
1094	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1095	IF ( ws_scheme_mom ) THEN
1096	CALL advec_w_ws
1097	ELSE
1098	CALL advec_w_pw
1099	ENDIF
1100	ELSE
1101	CALL advec_w_up
1102	ENDIF
1103	CALL diffusion_w
1104	CALL coriolis( 3 )
1105
1106	IF ( .NOT. neutral ) THEN
1107	IF ( ocean_mode ) THEN
1108	CALL buoyancy( rho_ocean, 3 )
1109	ELSE
1110	IF ( .NOT. humidity ) THEN
1111	CALL buoyancy( pt, 3 )
1112	ELSE
1113	CALL buoyancy( vpt, 3 )
1114	ENDIF
1115	ENDIF
1116	ENDIF
1117
1118	!
1119	!-- Drag by plant canopy
1120	IF ( plant_canopy ) CALL pcm_tendency( 3 )
1121
1122	!
1123	!-- Effect of Stokes drift (in ocean mode only)
1124	IF ( stokes_force ) CALL stokes_drift_terms( 3 )
1125
1126	CALL module_interface_actions( 'w-tendency' )
1127
1128	!
1129	!-- Prognostic equation for w-velocity component
1130	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1131	!$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_static_0) &
1132	!$ACC PRESENT(tsc(2:5)) &
1133	!$ACC PRESENT(w_p)
1134	DO i = nxl, nxr
1135	DO j = nys, nyn
1136	DO k = nzb+1, nzt-1
1137	w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1138	tsc(3) * tw_m(k,j,i) ) &
1139	- tsc(5) * rdf(k) * w(k,j,i) &
1140	) * MERGE( 1.0_wp, 0.0_wp, &
1141	BTEST( wall_flags_static_0(k,j,i), 3 )&
1142	)
1143	ENDDO
1144	ENDDO
1145	ENDDO
1146
1147	!
1148	!-- Calculate tendencies for the next Runge-Kutta step
1149	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1150	IF ( intermediate_timestep_count == 1 ) THEN
1151	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1152	!$ACC PRESENT(tend, tw_m)
1153	DO i = nxl, nxr
1154	DO j = nys, nyn
1155	DO k = nzb+1, nzt-1
1156	tw_m(k,j,i) = tend(k,j,i)
1157	ENDDO
1158	ENDDO
1159	ENDDO
1160	ELSEIF ( intermediate_timestep_count < &
1161	intermediate_timestep_count_max ) THEN
1162	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1163	!$ACC PRESENT(tend, tw_m)
1164	DO i = nxl, nxr
1165	DO j = nys, nyn
1166	DO k = nzb+1, nzt-1
1167	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1168	+ 5.3125_wp * tw_m(k,j,i)
1169	ENDDO
1170	ENDDO
1171	ENDDO
1172	ENDIF
1173	ENDIF
1174
1175	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1176
1177
1178	!
1179	!-- If required, compute prognostic equation for potential temperature
1180	IF ( .NOT. neutral ) THEN
1181
1182	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1183
1184	!
1185	!-- pt-tendency terms with communication
1186	sbt = tsc(2)
1187	IF ( scalar_advec == 'bc-scheme' ) THEN
1188
1189	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1190	!
1191	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1192	sbt = 1.0_wp
1193	ENDIF
1194	tend = 0.0_wp
1195	CALL advec_s_bc( pt, 'pt' )
1196
1197	ENDIF
1198
1199	!
1200	!-- pt-tendency terms with no communication
1201	IF ( scalar_advec /= 'bc-scheme' ) THEN
1202	!$ACC KERNELS PRESENT(tend)
1203	tend = 0.0_wp
1204	!$ACC END KERNELS
1205	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1206	IF ( ws_scheme_sca ) THEN
1207	CALL advec_s_ws( advc_flags_s, pt, 'pt', &
1208	bc_dirichlet_l .OR. bc_radiation_l, &
1209	bc_dirichlet_n .OR. bc_radiation_n, &
1210	bc_dirichlet_r .OR. bc_radiation_r, &
1211	bc_dirichlet_s .OR. bc_radiation_s )
1212	ELSE
1213	CALL advec_s_pw( pt )
1214	ENDIF
1215	ELSE
1216	CALL advec_s_up( pt )
1217	ENDIF
1218	ENDIF
1219
1220	CALL diffusion_s( pt, &
1221	surf_def_h(0)%shf, surf_def_h(1)%shf, &
1222	surf_def_h(2)%shf, &
1223	surf_lsm_h%shf, surf_usm_h%shf, &
1224	surf_def_v(0)%shf, surf_def_v(1)%shf, &
1225	surf_def_v(2)%shf, surf_def_v(3)%shf, &
1226	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
1227	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
1228	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
1229	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1230
1231	!
1232	!-- Consideration of heat sources within the plant canopy
1233	IF ( plant_canopy .AND. &
1234	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
1235	CALL pcm_tendency( 4 )
1236	ENDIF
1237
1238	!
1239	!-- Large scale advection
1240	IF ( large_scale_forcing ) THEN
1241	CALL ls_advec( simulated_time, 'pt' )
1242	ENDIF
1243
1244	!
1245	!-- Nudging
1246	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
1247
1248	!
1249	!-- If required compute influence of large-scale subsidence/ascent
1250	IF ( large_scale_subsidence .AND. &
1251	.NOT. use_subsidence_tendencies ) THEN
1252	CALL subsidence( tend, pt, pt_init, 2 )
1253	ENDIF
1254
1255	#if defined( __rrtmg )
1256	!
1257	!-- If required, add tendency due to radiative heating/cooling
1258	IF ( radiation .AND. &
1259	simulated_time > skip_time_do_radiation ) THEN
1260	CALL radiation_tendency ( tend )
1261	ENDIF
1262	#endif
1263
1264	CALL module_interface_actions( 'pt-tendency' )
1265
1266	!
1267	!-- Prognostic equation for potential temperature
1268	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1269	!$ACC PRESENT(pt, tend, tpt_m, wall_flags_static_0) &
1270	!$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1271	!$ACC PRESENT(tsc(3:5)) &
1272	!$ACC PRESENT(pt_p)
1273	DO i = nxl, nxr
1274	DO j = nys, nyn
1275	DO k = nzb+1, nzt
1276	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1277	tsc(3) * tpt_m(k,j,i) ) &
1278	- tsc(5) * &
1279	( pt(k,j,i) - pt_init(k) ) *&
1280	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1281	) &
1282	* MERGE( 1.0_wp, 0.0_wp, &
1283	BTEST( wall_flags_static_0(k,j,i), 0 ) &
1284	)
1285	ENDDO
1286	ENDDO
1287	ENDDO
1288	!
1289	!-- Calculate tendencies for the next Runge-Kutta step
1290	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1291	IF ( intermediate_timestep_count == 1 ) THEN
1292	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1293	!$ACC PRESENT(tend, tpt_m)
1294	DO i = nxl, nxr
1295	DO j = nys, nyn
1296	DO k = nzb+1, nzt
1297	tpt_m(k,j,i) = tend(k,j,i)
1298	ENDDO
1299	ENDDO
1300	ENDDO
1301	ELSEIF ( intermediate_timestep_count < &
1302	intermediate_timestep_count_max ) THEN
1303	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1304	!$ACC PRESENT(tend, tpt_m)
1305	DO i = nxl, nxr
1306	DO j = nys, nyn
1307	DO k = nzb+1, nzt
1308	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1309	5.3125_wp * tpt_m(k,j,i)
1310	ENDDO
1311	ENDDO
1312	ENDDO
1313	ENDIF
1314	ENDIF
1315
1316	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1317
1318	ENDIF
1319
1320	!
1321	!-- If required, compute prognostic equation for total water content
1322	IF ( humidity ) THEN
1323
1324	CALL cpu_log( log_point(29), 'q-equation', 'start' )
1325
1326	!
1327	!-- Scalar/q-tendency terms with communication
1328	sbt = tsc(2)
1329	IF ( scalar_advec == 'bc-scheme' ) THEN
1330
1331	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1332	!
1333	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1334	sbt = 1.0_wp
1335	ENDIF
1336	tend = 0.0_wp
1337	CALL advec_s_bc( q, 'q' )
1338
1339	ENDIF
1340
1341	!
1342	!-- Scalar/q-tendency terms with no communication
1343	IF ( scalar_advec /= 'bc-scheme' ) THEN
1344	tend = 0.0_wp
1345	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1346	IF ( ws_scheme_sca ) THEN
1347	CALL advec_s_ws( advc_flags_s, q, 'q', &
1348	bc_dirichlet_l .OR. bc_radiation_l, &
1349	bc_dirichlet_n .OR. bc_radiation_n, &
1350	bc_dirichlet_r .OR. bc_radiation_r, &
1351	bc_dirichlet_s .OR. bc_radiation_s )
1352	ELSE
1353	CALL advec_s_pw( q )
1354	ENDIF
1355	ELSE
1356	CALL advec_s_up( q )
1357	ENDIF
1358	ENDIF
1359
1360	CALL diffusion_s( q, &
1361	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
1362	surf_def_h(2)%qsws, &
1363	surf_lsm_h%qsws, surf_usm_h%qsws, &
1364	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
1365	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
1366	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
1367	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
1368	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
1369	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1370
1371	!
1372	!-- Sink or source of humidity due to canopy elements
1373	IF ( plant_canopy ) CALL pcm_tendency( 5 )
1374
1375	!
1376	!-- Large scale advection
1377	IF ( large_scale_forcing ) THEN
1378	CALL ls_advec( simulated_time, 'q' )
1379	ENDIF
1380
1381	!
1382	!-- Nudging
1383	IF ( nudging ) CALL nudge( simulated_time, 'q' )
1384
1385	!
1386	!-- If required compute influence of large-scale subsidence/ascent
1387	IF ( large_scale_subsidence .AND. &
1388	.NOT. use_subsidence_tendencies ) THEN
1389	CALL subsidence( tend, q, q_init, 3 )
1390	ENDIF
1391
1392	CALL module_interface_actions( 'q-tendency' )
1393
1394	!
1395	!-- Prognostic equation for total water content
1396	DO i = nxl, nxr
1397	DO j = nys, nyn
1398	DO k = nzb+1, nzt
1399	q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1400	tsc(3) * tq_m(k,j,i) ) &
1401	- tsc(5) * rdf_sc(k) * &
1402	( q(k,j,i) - q_init(k) ) &
1403	) * MERGE( 1.0_wp, 0.0_wp, &
1404	BTEST( wall_flags_static_0(k,j,i), 0 ) &
1405	)
1406	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
1407	ENDDO
1408	ENDDO
1409	ENDDO
1410
1411	!
1412	!-- Calculate tendencies for the next Runge-Kutta step
1413	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1414	IF ( intermediate_timestep_count == 1 ) THEN
1415	DO i = nxl, nxr
1416	DO j = nys, nyn
1417	DO k = nzb+1, nzt
1418	tq_m(k,j,i) = tend(k,j,i)
1419	ENDDO
1420	ENDDO
1421	ENDDO
1422	ELSEIF ( intermediate_timestep_count < &
1423	intermediate_timestep_count_max ) THEN
1424	DO i = nxl, nxr
1425	DO j = nys, nyn
1426	DO k = nzb+1, nzt
1427	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1428	+ 5.3125_wp * tq_m(k,j,i)
1429	ENDDO
1430	ENDDO
1431	ENDDO
1432	ENDIF
1433	ENDIF
1434
1435	CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1436
1437	ENDIF
1438	!
1439	!-- If required, compute prognostic equation for scalar
1440	IF ( passive_scalar ) THEN
1441
1442	CALL cpu_log( log_point(66), 's-equation', 'start' )
1443
1444	!
1445	!-- Scalar/q-tendency terms with communication
1446	sbt = tsc(2)
1447	IF ( scalar_advec == 'bc-scheme' ) THEN
1448
1449	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1450	!
1451	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1452	sbt = 1.0_wp
1453	ENDIF
1454	tend = 0.0_wp
1455	CALL advec_s_bc( s, 's' )
1456
1457	ENDIF
1458
1459	!
1460	!-- Scalar/q-tendency terms with no communication
1461	IF ( scalar_advec /= 'bc-scheme' ) THEN
1462	tend = 0.0_wp
1463	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1464	IF ( ws_scheme_sca ) THEN
1465	CALL advec_s_ws( advc_flags_s, s, 's', &
1466	bc_dirichlet_l .OR. bc_radiation_l, &
1467	bc_dirichlet_n .OR. bc_radiation_n, &
1468	bc_dirichlet_r .OR. bc_radiation_r, &
1469	bc_dirichlet_s .OR. bc_radiation_s )
1470	ELSE
1471	CALL advec_s_pw( s )
1472	ENDIF
1473	ELSE
1474	CALL advec_s_up( s )
1475	ENDIF
1476	ENDIF
1477
1478	CALL diffusion_s( s, &
1479	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1480	surf_def_h(2)%ssws, &
1481	surf_lsm_h%ssws, surf_usm_h%ssws, &
1482	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1483	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1484	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1485	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1486	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1487	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1488
1489	!
1490	!-- Sink or source of humidity due to canopy elements
1491	IF ( plant_canopy ) CALL pcm_tendency( 7 )
1492
1493	!
1494	!-- Large scale advection. Not implemented for scalars so far.
1495	! IF ( large_scale_forcing ) THEN
1496	! CALL ls_advec( simulated_time, 'q' )
1497	! ENDIF
1498
1499	!
1500	!-- Nudging. Not implemented for scalars so far.
1501	! IF ( nudging ) CALL nudge( simulated_time, 'q' )
1502
1503	!
1504	!-- If required compute influence of large-scale subsidence/ascent.
1505	!-- Not implemented for scalars so far.
1506	IF ( large_scale_subsidence .AND. &
1507	.NOT. use_subsidence_tendencies .AND. &
1508	.NOT. large_scale_forcing ) THEN
1509	CALL subsidence( tend, s, s_init, 3 )
1510	ENDIF
1511
1512	CALL module_interface_actions( 's-tendency' )
1513
1514	!
1515	!-- Prognostic equation for total water content
1516	DO i = nxl, nxr
1517	DO j = nys, nyn
1518	DO k = nzb+1, nzt
1519	s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1520	tsc(3) * ts_m(k,j,i) ) &
1521	- tsc(5) * rdf_sc(k) * &
1522	( s(k,j,i) - s_init(k) ) &
1523	) &
1524	* MERGE( 1.0_wp, 0.0_wp, &
1525	BTEST( wall_flags_static_0(k,j,i), 0 ) &
1526	)
1527	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1528	ENDDO
1529	ENDDO
1530	ENDDO
1531
1532	!
1533	!-- Calculate tendencies for the next Runge-Kutta step
1534	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1535	IF ( intermediate_timestep_count == 1 ) THEN
1536	DO i = nxl, nxr
1537	DO j = nys, nyn
1538	DO k = nzb+1, nzt
1539	ts_m(k,j,i) = tend(k,j,i)
1540	ENDDO
1541	ENDDO
1542	ENDDO
1543	ELSEIF ( intermediate_timestep_count < &
1544	intermediate_timestep_count_max ) THEN
1545	DO i = nxl, nxr
1546	DO j = nys, nyn
1547	DO k = nzb+1, nzt
1548	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1549	+ 5.3125_wp * ts_m(k,j,i)
1550	ENDDO
1551	ENDDO
1552	ENDDO
1553	ENDIF
1554	ENDIF
1555
1556	CALL cpu_log( log_point(66), 's-equation', 'stop' )
1557
1558	ENDIF
1559	!
1560	!-- Calculate prognostic equations for all other modules
1561	CALL module_interface_prognostic_equations()
1562
1563	IF ( debug_output_timestep ) CALL debug_message( 'prognostic_equations_vector', 'end' )
1564
1565	END SUBROUTINE prognostic_equations_vector
1566
1567
1568	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |