Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |