Home

Context Navigation

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

Last change on this file since 3910 was 3899, checked in by monakurppa, 6 years ago
corrected the OpenMP implementation for salsa and some minor bugs in salsa_mod
Property svn:keywords set to `Id`
File size: 69.4 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 3899 2019-04-16 14:05:27Z suehring $
27	! Corrections in the OpenMP version of salsa
28	!
29	! 3887 2019-04-12 08:47:41Z schwenkel
30	! Implicit Bugfix for chemistry model, loop for non_transport_physics over
31	! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
32	!
33	! 3885 2019-04-11 11:29:34Z kanani
34	! Changes related to global restructuring of location messages and introduction
35	! of additional debug messages
36	!
37	! 3881 2019-04-10 09:31:22Z suehring
38	! Bugfix in OpenMP directive
39	!
40	! 3880 2019-04-08 21:43:02Z knoop
41	! Moved wtm_tendencies to module_interface_actions
42	!
43	! 3874 2019-04-08 16:53:48Z knoop
44	! Added non_transport_physics module interfaces and moved bcm code into it
45	!
46	! 3872 2019-04-08 15:03:06Z knoop
47	! Moving prognostic equations of bcm into bulk_cloud_model_mod
48	!
49	! 3864 2019-04-05 09:01:56Z monakurppa
50	! Modifications made for salsa:
51	! - salsa_prognostic_equations moved to salsa_mod (and the call to
52	! module_interface_mod)
53	! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
54	! ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
55	!
56	! 3840 2019-03-29 10:35:52Z knoop
57	! added USE chem_gasphase_mod for nspec, nspec and spc_names
58	!
59	! 3820 2019-03-27 11:53:41Z forkel
60	! renamed do_depo to deposition_dry (ecc)
61	!
62	! 3797 2019-03-15 11:15:38Z forkel
63	! Call chem_integegrate in OpenMP loop (ketelsen)
64	!
65	!
66	! 3771 2019-02-28 12:19:33Z raasch
67	! preprocessor directivs fro rrtmg added
68	!
69	! 3761 2019-02-25 15:31:42Z raasch
70	! unused variable removed
71	!
72	! 3719 2019-02-06 13:10:18Z kanani
73	! Cleaned up chemistry cpu measurements
74	!
75	! 3684 2019-01-20 20:20:58Z knoop
76	! OpenACC port for SPEC
77	!
78	! 3589 2018-11-30 15:09:51Z suehring
79	! Move the control parameter "salsa" from salsa_mod to control_parameters
80	! (M. Kurppa)
81	!
82	! 3582 2018-11-29 19:16:36Z suehring
83	! Implementation of a new aerosol module salsa.
84	! Remove cpu-logs from i,j loop in cache version.
85	!
86	! 3458 2018-10-30 14:51:23Z kanani
87	! remove duplicate USE chem_modules
88	! from chemistry branch r3443, banzhafs, basit:
89	! chem_depo call introduced
90	! code added for decycling chemistry
91	!
92	! 3386 2018-10-19 16:28:22Z gronemeier
93	! Renamed tcm_prognostic to tcm_prognostic_equations
94	!
95	! 3355 2018-10-16 14:03:34Z knoop
96	! (from branch resler)
97	! Fix for chemistry call
98	!
99	! 3302 2018-10-03 02:39:40Z raasch
100	! Stokes drift + wave breaking term added
101	!
102	! 3298 2018-10-02 12:21:11Z kanani
103	! Code added for decycling chemistry (basit)
104	!
105	! 3294 2018-10-01 02:37:10Z raasch
106	! changes concerning modularization of ocean option
107	!
108	! 3274 2018-09-24 15:42:55Z knoop
109	! Modularization of all bulk cloud physics code components
110	!
111	! 3241 2018-09-12 15:02:00Z raasch
112	! omp_get_thread_num now declared in omp directive
113	!
114	! 3183 2018-07-27 14:25:55Z suehring
115	! Remove unused variables from USE statements
116	!
117	! 3182 2018-07-27 13:36:03Z suehring
118	! Revise recent bugfix for nesting
119	!
120	! 3021 2018-05-16 08:14:20Z maronga
121	! Bugfix in IF clause for nesting
122	!
123	! 3014 2018-05-09 08:42:38Z maronga
124	! Fixed a bug in the IF condition to call pcm_tendency in case of
125	! potential temperature
126	!
127	! 2815 2018-02-19 11:29:57Z kanani
128	! Rename chem_tendency to chem_prognostic_equations,
129	! implement vector version for air chemistry
130	!
131	! 2766 2018-01-22 17:17:47Z kanani
132	! Removed preprocessor directive __chem
133	!
134	! 2746 2018-01-15 12:06:04Z suehring
135	! Move flag plant canopy to modules
136	!
137	! 2719 2018-01-02 09:02:06Z maronga
138	! Bugfix for last change.
139	!
140	! 2718 2018-01-02 08:49:38Z maronga
141	! Corrected "Former revisions" section
142	!
143	! 2696 2017-12-14 17:12:51Z kanani
144	! - Change in file header (GPL part)
145	! - Moved TKE equation to tcm_prognostic (TG)
146	! - Added switch for chemical reactions (RF, FK)
147	! - Implementation of chemistry module (RF, BK, FK)
148	!
149	! 2563 2017-10-19 15:36:10Z Giersch
150	! Variable wind_turbine moved to module control_parameters
151	!
152	! 2320 2017-07-21 12:47:43Z suehring
153	! Modularize large-scale forcing and nudging
154	!
155	! 2292 2017-06-20 09:51:42Z schwenkel
156	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
157	! includes two more prognostic equations for cloud drop concentration (nc)
158	! and cloud water content (qc).
159	!
160	! 2261 2017-06-08 14:25:57Z raasch
161	! bugfix for r2232: openmp directives removed
162	!
163	! 2233 2017-05-30 18:08:54Z suehring
164	!
165	! 2232 2017-05-30 17:47:52Z suehring
166	! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
167	! which is realized directly in diffusion_s now.
168	!
169	! 2192 2017-03-22 04:14:10Z raasch
170	! Bugfix for misplaced and missing openMP directives from r2155
171	!
172	! 2155 2017-02-21 09:57:40Z hoffmann
173	! Bugfix in the calculation of microphysical quantities on ghost points.
174	!
175	! 2118 2017-01-17 16:38:49Z raasch
176	! OpenACC version of subroutine removed
177	!
178	! 2031 2016-10-21 15:11:58Z knoop
179	! renamed variable rho to rho_ocean
180	!
181	! 2011 2016-09-19 17:29:57Z kanani
182	! Flag urban_surface is now defined in module control_parameters.
183	!
184	! 2007 2016-08-24 15:47:17Z kanani
185	! Added pt tendency calculation based on energy balance at urban surfaces
186	! (new urban surface model)
187	!
188	! 2000 2016-08-20 18:09:15Z knoop
189	! Forced header and separation lines into 80 columns
190	!
191	! 1976 2016-07-27 13:28:04Z maronga
192	! Simplied calls to radiation model
193	!
194	! 1960 2016-07-12 16:34:24Z suehring
195	! Separate humidity and passive scalar
196	!
197	! 1914 2016-05-26 14:44:07Z witha
198	! Added calls for wind turbine model
199	!
200	! 1873 2016-04-18 14:50:06Z maronga
201	! Module renamed (removed _mod)
202	!
203	! 1850 2016-04-08 13:29:27Z maronga
204	! Module renamed
205	!
206	! 1826 2016-04-07 12:01:39Z maronga
207	! Renamed canopy model calls.
208	!
209	! 1822 2016-04-07 07:49:42Z hoffmann
210	! Kessler microphysics scheme moved to microphysics.
211	!
212	! 1757 2016-02-22 15:49:32Z maronga
213	!
214	! 1691 2015-10-26 16:17:44Z maronga
215	! Added optional model spin-up without radiation / land surface model calls.
216	! Formatting corrections.
217	!
218	! 1682 2015-10-07 23:56:08Z knoop
219	! Code annotations made doxygen readable
220	!
221	! 1585 2015-04-30 07:05:52Z maronga
222	! Added call for temperature tendency calculation due to radiative flux divergence
223	!
224	! 1517 2015-01-07 19:12:25Z hoffmann
225	! advec_s_bc_mod addded, since advec_s_bc is now a module
226	!
227	! 1484 2014-10-21 10:53:05Z kanani
228	! Changes due to new module structure of the plant canopy model:
229	! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
230	! Removed double-listing of use_upstream_for_tke in ONLY-list of module
231	! control_parameters
232	!
233	! 1409 2014-05-23 12:11:32Z suehring
234	! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
235	! This ensures that left-hand side fluxes are also calculated for nxl in that
236	! case, even though the solution at nxl is overwritten in boundary_conds()
237	!
238	! 1398 2014-05-07 11:15:00Z heinze
239	! Rayleigh-damping for horizontal velocity components changed: instead of damping
240	! against ug and vg, damping against u_init and v_init is used to allow for a
241	! homogenized treatment in case of nudging
242	!
243	! 1380 2014-04-28 12:40:45Z heinze
244	! Change order of calls for scalar prognostic quantities:
245	! ls_advec -> nudging -> subsidence since initial profiles
246	!
247	! 1374 2014-04-25 12:55:07Z raasch
248	! missing variables added to ONLY lists
249	!
250	! 1365 2014-04-22 15:03:56Z boeske
251	! Calls of ls_advec for large scale advection added,
252	! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
253	! new argument ls_index added to the calls of subsidence
254	! +ls_index
255	!
256	! 1361 2014-04-16 15:17:48Z hoffmann
257	! Two-moment microphysics moved to the start of prognostic equations. This makes
258	! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
259	! Additionally, it is allowed to call the microphysics just once during the time
260	! step (not at each sub-time step).
261	!
262	! Two-moment cloud physics added for vector and accelerator optimization.
263	!
264	! 1353 2014-04-08 15:21:23Z heinze
265	! REAL constants provided with KIND-attribute
266	!
267	! 1337 2014-03-25 15:11:48Z heinze
268	! Bugfix: REAL constants provided with KIND-attribute
269	!
270	! 1332 2014-03-25 11:59:43Z suehring
271	! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
272	!
273	! 1330 2014-03-24 17:29:32Z suehring
274	! In case of SGS-particle velocity advection of TKE is also allowed with
275	! dissipative 5th-order scheme.
276	!
277	! 1320 2014-03-20 08:40:49Z raasch
278	! ONLY-attribute added to USE-statements,
279	! kind-parameters added to all INTEGER and REAL declaration statements,
280	! kinds are defined in new module kinds,
281	! old module precision_kind is removed,
282	! revision history before 2012 removed,
283	! comment fields (!:) to be used for variable explanations added to
284	! all variable declaration statements
285	!
286	! 1318 2014-03-17 13:35:16Z raasch
287	! module interfaces removed
288	!
289	! 1257 2013-11-08 15:18:40Z raasch
290	! openacc loop vector clauses removed, independent clauses added
291	!
292	! 1246 2013-11-01 08:59:45Z heinze
293	! enable nudging also for accelerator version
294	!
295	! 1241 2013-10-30 11:36:58Z heinze
296	! usage of nudging enabled (so far not implemented for accelerator version)
297	!
298	! 1179 2013-06-14 05:57:58Z raasch
299	! two arguments removed from routine buoyancy, ref_state updated on device
300	!
301	! 1128 2013-04-12 06:19:32Z raasch
302	! those parts requiring global communication moved to time_integration,
303	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
304	! j_north
305	!
306	! 1115 2013-03-26 18:16:16Z hoffmann
307	! optimized cloud physics: calculation of microphysical tendencies transfered
308	! to microphysics.f90; qr and nr are only calculated if precipitation is required
309	!
310	! 1111 2013-03-08 23:54:10Z raasch
311	! update directives for prognostic quantities removed
312	!
313	! 1106 2013-03-04 05:31:38Z raasch
314	! small changes in code formatting
315	!
316	! 1092 2013-02-02 11:24:22Z raasch
317	! unused variables removed
318	!
319	! 1053 2012-11-13 17:11:03Z hoffmann
320	! implementation of two new prognostic equations for rain drop concentration (nr)
321	! and rain water content (qr)
322	!
323	! currently, only available for cache loop optimization
324	!
325	! 1036 2012-10-22 13:43:42Z raasch
326	! code put under GPL (PALM 3.9)
327	!
328	! 1019 2012-09-28 06:46:45Z raasch
329	! non-optimized version of prognostic_equations removed
330	!
331	! 1015 2012-09-27 09:23:24Z raasch
332	! new branch prognostic_equations_acc
333	! OpenACC statements added + code changes required for GPU optimization
334	!
335	! 1001 2012-09-13 14:08:46Z raasch
336	! all actions concerning leapfrog- and upstream-spline-scheme removed
337	!
338	! 978 2012-08-09 08:28:32Z fricke
339	! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
340	! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
341	! boundary in case of non-cyclic lateral boundaries
342	! Bugfix: first thread index changes for WS-scheme at the inflow
343	!
344	! 940 2012-07-09 14:31:00Z raasch
345	! temperature equation can be switched off
346	!
347	! Revision 1.1 2000/04/13 14:56:27 schroeter
348	! Initial revision
349	!
350	!
351	! Description:
352	! ------------
353	!> Solving the prognostic equations.
354	!------------------------------------------------------------------------------!
355	MODULE prognostic_equations_mod
356
357	USE advec_s_bc_mod, &
358	ONLY: advec_s_bc
359
360	USE advec_s_pw_mod, &
361	ONLY: advec_s_pw
362
363	USE advec_s_up_mod, &
364	ONLY: advec_s_up
365
366	USE advec_u_pw_mod, &
367	ONLY: advec_u_pw
368
369	USE advec_u_up_mod, &
370	ONLY: advec_u_up
371
372	USE advec_v_pw_mod, &
373	ONLY: advec_v_pw
374
375	USE advec_v_up_mod, &
376	ONLY: advec_v_up
377
378	USE advec_w_pw_mod, &
379	ONLY: advec_w_pw
380
381	USE advec_w_up_mod, &
382	ONLY: advec_w_up
383
384	USE advec_ws, &
385	ONLY: advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
386
387	USE arrays_3d, &
388	ONLY: diss_l_e, diss_l_pt, diss_l_q, &
389	diss_l_s, diss_l_sa, diss_s_e, &
390	diss_s_pt, diss_s_q, diss_s_s, diss_s_sa, &
391	e, e_p, flux_s_e, flux_s_pt, flux_s_q, &
392	flux_s_s, flux_s_sa, flux_l_e, &
393	flux_l_pt, flux_l_q, flux_l_s, &
394	flux_l_sa, pt, ptdf_x, ptdf_y, pt_init, &
395	pt_p, prho, q, q_init, q_p, rdf, rdf_sc, &
396	ref_state, rho_ocean, s, s_init, s_p, tend, te_m, &
397	tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u, &
398	ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
399
400	USE buoyancy_mod, &
401	ONLY: buoyancy
402
403	USE chemistry_model_mod, &
404	ONLY: chem_boundary_conds_decycle
405
406	USE control_parameters, &
407	ONLY: air_chemistry, constant_diffusion, &
408	debug_output, debug_string, &
409	dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, &
410	humidity, intermediate_timestep_count, &
411	intermediate_timestep_count_max, large_scale_forcing, &
412	large_scale_subsidence, neutral, nudging, &
413	ocean_mode, passive_scalar, plant_canopy, pt_reference, &
414	scalar_advec, scalar_advec, simulated_time, sloping_surface, &
415	timestep_scheme, tsc, use_subsidence_tendencies, &
416	use_upstream_for_tke, wind_turbine, ws_scheme_mom, &
417	ws_scheme_sca, urban_surface, land_surface, &
418	time_since_reference_point, salsa
419
420	USE coriolis_mod, &
421	ONLY: coriolis
422
423	USE cpulog, &
424	ONLY: cpu_log, log_point, log_point_s
425
426	USE diffusion_s_mod, &
427	ONLY: diffusion_s
428
429	USE diffusion_u_mod, &
430	ONLY: diffusion_u
431
432	USE diffusion_v_mod, &
433	ONLY: diffusion_v
434
435	USE diffusion_w_mod, &
436	ONLY: diffusion_w
437
438	USE indices, &
439	ONLY: nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, &
440	nzb, nzt, wall_flags_0
441
442	USE kinds
443
444	USE lsf_nudging_mod, &
445	ONLY: ls_advec, nudge
446
447	USE module_interface, &
448	ONLY: module_interface_actions, &
449	module_interface_non_transport_physics, &
450	module_interface_exchange_horiz, &
451	module_interface_prognostic_equations
452
453	USE ocean_mod, &
454	ONLY: stokes_drift_terms, stokes_force, &
455	wave_breaking, wave_breaking_term
456
457	USE plant_canopy_model_mod, &
458	ONLY: cthf, pcm_tendency
459
460	#if defined( __rrtmg )
461	USE radiation_model_mod, &
462	ONLY: radiation, radiation_tendency, &
463	skip_time_do_radiation
464	#endif
465
466	USE salsa_mod, &
467	ONLY: aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, &
468	nbins_aerosol, ncomponents_mass, ngases_salsa, &
469	salsa_boundary_conds, salsa_diagnostics, salsa_driver, &
470	salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
471
472	USE statistics, &
473	ONLY: hom
474
475	USE subsidence_mod, &
476	ONLY: subsidence
477
478	USE surface_mod, &
479	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
480	surf_usm_v
481
482	USE turbulence_closure_mod, &
483	ONLY: tcm_prognostic_equations
484
485	IMPLICIT NONE
486
487	PRIVATE
488	PUBLIC prognostic_equations_cache, prognostic_equations_vector
489
490	INTERFACE prognostic_equations_cache
491	MODULE PROCEDURE prognostic_equations_cache
492	END INTERFACE prognostic_equations_cache
493
494	INTERFACE prognostic_equations_vector
495	MODULE PROCEDURE prognostic_equations_vector
496	END INTERFACE prognostic_equations_vector
497
498
499	CONTAINS
500
501
502	!------------------------------------------------------------------------------!
503	! Description:
504	! ------------
505	!> Version with one optimized loop over all equations. It is only allowed to
506	!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
507	!>
508	!> Here the calls of most subroutines are embedded in two DO loops over i and j,
509	!> so communication between CPUs is not allowed (does not make sense) within
510	!> these loops.
511	!>
512	!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
513	!------------------------------------------------------------------------------!
514
515	SUBROUTINE prognostic_equations_cache
516
517
518	INTEGER(iwp) :: i !<
519	INTEGER(iwp) :: i_omp_start !<
520	INTEGER(iwp) :: ib !< index for aerosol size bins (salsa)
521	INTEGER(iwp) :: ic !< index for chemical compounds (salsa)
522	INTEGER(iwp) :: icc !< additional index for chemical compounds (salsa)
523	INTEGER(iwp) :: ig !< index for gaseous compounds (salsa)
524	INTEGER(iwp) :: j !<
525	INTEGER(iwp) :: k !<
526	!$ INTEGER(iwp) :: omp_get_thread_num !<
527	INTEGER(iwp) :: tn = 0 !<
528
529	LOGICAL :: loop_start !<
530
531
532
533	IF ( debug_output ) THEN
534	WRITE( debug_string, * ) 'prognostic_equations_cache'
535	CALL debug_message( debug_string, 'start' )
536	ENDIF
537	!
538	!-- Time measurement can only be performed for the whole set of equations
539	CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
540
541	!$OMP PARALLEL PRIVATE (i,j)
542	!$OMP DO
543	DO i = nxl, nxr
544	DO j = nys, nyn
545	!
546	!-- Calculate non transport physics for all other modules
547	CALL module_interface_non_transport_physics( i, j )
548	ENDDO
549	ENDDO
550	!
551	!-- Module Inferface for exchange horiz after non_transport_physics but before
552	!-- advection. Therefore, non_transport_physics must not run for ghost points.
553	CALL module_interface_exchange_horiz()
554	!$OMP END PARALLEL
555
556	IF ( air_chemistry ) CALL chem_boundary_conds_decycle
557	!
558	!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
559	!-- step. The exchange of ghost points is required after this update of the
560	!-- concentrations of aerosol number and mass
561	IF ( salsa ) THEN
562
563	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
564	IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa ) &
565	THEN
566	CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
567	!$OMP PARALLEL PRIVATE (i,j)
568	!
569	!-- Call salsa processes
570	!$OMP DO
571	DO i = nxl, nxr
572	DO j = nys, nyn
573	CALL salsa_diagnostics( i, j )
574	CALL salsa_driver( i, j, 3 )
575	CALL salsa_diagnostics( i, j )
576	ENDDO
577	ENDDO
578	!$OMP END PARALLEL
579
580	CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
581
582	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
583	!
584	!-- Exchange ghost points and decycle if needed.
585	DO ib = 1, nbins_aerosol
586	CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
587	CALL salsa_boundary_conds( aerosol_number(ib)%conc, &
588	aerosol_number(ib)%init )
589	DO ic = 1, ncomponents_mass
590	icc = ( ic - 1 ) * nbins_aerosol + ib
591	CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
592	CALL salsa_boundary_conds( aerosol_mass(icc)%conc, &
593	aerosol_mass(icc)%init )
594	ENDDO
595	ENDDO
596
597	IF ( .NOT. salsa_gases_from_chem ) THEN
598	DO ig = 1, ngases_salsa
599	CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
600	CALL salsa_boundary_conds( salsa_gas(ig)%conc, &
601	salsa_gas(ig)%init )
602	ENDDO
603	ENDIF
604	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
605	last_salsa_time = time_since_reference_point
606
607	ENDIF
608
609	ENDIF
610
611	ENDIF
612
613	!
614	!-- Loop over all prognostic equations
615	!$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
616
617	!$ tn = omp_get_thread_num()
618	loop_start = .TRUE.
619
620	!$OMP DO
621	DO i = nxl, nxr
622
623	!
624	!-- Store the first loop index. It differs for each thread and is required
625	!-- later in advec_ws
626	IF ( loop_start ) THEN
627	loop_start = .FALSE.
628	i_omp_start = i
629	ENDIF
630
631	DO j = nys, nyn
632	!
633	!-- Tendency terms for u-velocity component. Please note, in case of
634	!-- non-cyclic boundary conditions the grid point i=0 is excluded from
635	!-- the prognostic equations for the u-component.
636	IF ( i >= nxlu ) THEN
637
638	tend(:,j,i) = 0.0_wp
639	IF ( timestep_scheme(1:5) == 'runge' ) THEN
640	IF ( ws_scheme_mom ) THEN
641	CALL advec_u_ws( i, j, i_omp_start, tn )
642	ELSE
643	CALL advec_u_pw( i, j )
644	ENDIF
645	ELSE
646	CALL advec_u_up( i, j )
647	ENDIF
648	CALL diffusion_u( i, j )
649	CALL coriolis( i, j, 1 )
650	IF ( sloping_surface .AND. .NOT. neutral ) THEN
651	CALL buoyancy( i, j, pt, 1 )
652	ENDIF
653
654	!
655	!-- Drag by plant canopy
656	IF ( plant_canopy ) CALL pcm_tendency( i, j, 1 )
657
658	!
659	!-- External pressure gradient
660	IF ( dp_external ) THEN
661	DO k = dp_level_ind_b+1, nzt
662	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
663	ENDDO
664	ENDIF
665
666	!
667	!-- Nudging
668	IF ( nudging ) CALL nudge( i, j, simulated_time, 'u' )
669
670	!
671	!-- Effect of Stokes drift (in ocean mode only)
672	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 1 )
673
674	CALL module_interface_actions( i, j, 'u-tendency' )
675	!
676	!-- Prognostic equation for u-velocity component
677	DO k = nzb+1, nzt
678
679	u_p(k,j,i) = u(k,j,i) + ( dt_3d * &
680	( tsc(2) * tend(k,j,i) + &
681	tsc(3) * tu_m(k,j,i) ) &
682	- tsc(5) * rdf(k) &
683	* ( u(k,j,i) - u_init(k) ) &
684	) * MERGE( 1.0_wp, 0.0_wp, &
685	BTEST( wall_flags_0(k,j,i), 1 )&
686	)
687	ENDDO
688
689	!
690	!-- Add turbulence generated by wave breaking (in ocean mode only)
691	IF ( wave_breaking .AND. &
692	intermediate_timestep_count == intermediate_timestep_count_max )&
693	THEN
694	CALL wave_breaking_term( i, j, 1 )
695	ENDIF
696
697	!
698	!-- Calculate tendencies for the next Runge-Kutta step
699	IF ( timestep_scheme(1:5) == 'runge' ) THEN
700	IF ( intermediate_timestep_count == 1 ) THEN
701	DO k = nzb+1, nzt
702	tu_m(k,j,i) = tend(k,j,i)
703	ENDDO
704	ELSEIF ( intermediate_timestep_count < &
705	intermediate_timestep_count_max ) THEN
706	DO k = nzb+1, nzt
707	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
708	+ 5.3125_wp * tu_m(k,j,i)
709	ENDDO
710	ENDIF
711	ENDIF
712
713	ENDIF
714	!
715	!-- Tendency terms for v-velocity component. Please note, in case of
716	!-- non-cyclic boundary conditions the grid point j=0 is excluded from
717	!-- the prognostic equations for the v-component. !--
718	IF ( j >= nysv ) THEN
719
720	tend(:,j,i) = 0.0_wp
721	IF ( timestep_scheme(1:5) == 'runge' ) THEN
722	IF ( ws_scheme_mom ) THEN
723	CALL advec_v_ws( i, j, i_omp_start, tn )
724	ELSE
725	CALL advec_v_pw( i, j )
726	ENDIF
727	ELSE
728	CALL advec_v_up( i, j )
729	ENDIF
730	CALL diffusion_v( i, j )
731	CALL coriolis( i, j, 2 )
732
733	!
734	!-- Drag by plant canopy
735	IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 )
736
737	!
738	!-- External pressure gradient
739	IF ( dp_external ) THEN
740	DO k = dp_level_ind_b+1, nzt
741	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
742	ENDDO
743	ENDIF
744
745	!
746	!-- Nudging
747	IF ( nudging ) CALL nudge( i, j, simulated_time, 'v' )
748
749	!
750	!-- Effect of Stokes drift (in ocean mode only)
751	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 2 )
752
753	CALL module_interface_actions( i, j, 'v-tendency' )
754	!
755	!-- Prognostic equation for v-velocity component
756	DO k = nzb+1, nzt
757	v_p(k,j,i) = v(k,j,i) + ( dt_3d * &
758	( tsc(2) * tend(k,j,i) + &
759	tsc(3) * tv_m(k,j,i) ) &
760	- tsc(5) * rdf(k) &
761	* ( v(k,j,i) - v_init(k) )&
762	) * MERGE( 1.0_wp, 0.0_wp, &
763	BTEST( wall_flags_0(k,j,i), 2 )&
764	)
765	ENDDO
766
767	!
768	!-- Add turbulence generated by wave breaking (in ocean mode only)
769	IF ( wave_breaking .AND. &
770	intermediate_timestep_count == intermediate_timestep_count_max )&
771	THEN
772	CALL wave_breaking_term( i, j, 2 )
773	ENDIF
774
775	!
776	!-- Calculate tendencies for the next Runge-Kutta step
777	IF ( timestep_scheme(1:5) == 'runge' ) THEN
778	IF ( intermediate_timestep_count == 1 ) THEN
779	DO k = nzb+1, nzt
780	tv_m(k,j,i) = tend(k,j,i)
781	ENDDO
782	ELSEIF ( intermediate_timestep_count < &
783	intermediate_timestep_count_max ) THEN
784	DO k = nzb+1, nzt
785	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
786	+ 5.3125_wp * tv_m(k,j,i)
787	ENDDO
788	ENDIF
789	ENDIF
790
791	ENDIF
792
793	!
794	!-- Tendency terms for w-velocity component
795	tend(:,j,i) = 0.0_wp
796	IF ( timestep_scheme(1:5) == 'runge' ) THEN
797	IF ( ws_scheme_mom ) THEN
798	CALL advec_w_ws( i, j, i_omp_start, tn )
799	ELSE
800	CALL advec_w_pw( i, j )
801	END IF
802	ELSE
803	CALL advec_w_up( i, j )
804	ENDIF
805	CALL diffusion_w( i, j )
806	CALL coriolis( i, j, 3 )
807
808	IF ( .NOT. neutral ) THEN
809	IF ( ocean_mode ) THEN
810	CALL buoyancy( i, j, rho_ocean, 3 )
811	ELSE
812	IF ( .NOT. humidity ) THEN
813	CALL buoyancy( i, j, pt, 3 )
814	ELSE
815	CALL buoyancy( i, j, vpt, 3 )
816	ENDIF
817	ENDIF
818	ENDIF
819
820	!
821	!-- Drag by plant canopy
822	IF ( plant_canopy ) CALL pcm_tendency( i, j, 3 )
823
824	!
825	!-- Effect of Stokes drift (in ocean mode only)
826	IF ( stokes_force ) CALL stokes_drift_terms( i, j, 3 )
827
828	CALL module_interface_actions( i, j, 'w-tendency' )
829	!
830	!-- Prognostic equation for w-velocity component
831	DO k = nzb+1, nzt-1
832	w_p(k,j,i) = w(k,j,i) + ( dt_3d * &
833	( tsc(2) * tend(k,j,i) + &
834	tsc(3) * tw_m(k,j,i) ) &
835	- tsc(5) * rdf(k) * w(k,j,i) &
836	) * MERGE( 1.0_wp, 0.0_wp, &
837	BTEST( wall_flags_0(k,j,i), 3 )&
838	)
839	ENDDO
840
841	!
842	!-- Calculate tendencies for the next Runge-Kutta step
843	IF ( timestep_scheme(1:5) == 'runge' ) THEN
844	IF ( intermediate_timestep_count == 1 ) THEN
845	DO k = nzb+1, nzt-1
846	tw_m(k,j,i) = tend(k,j,i)
847	ENDDO
848	ELSEIF ( intermediate_timestep_count < &
849	intermediate_timestep_count_max ) THEN
850	DO k = nzb+1, nzt-1
851	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
852	+ 5.3125_wp * tw_m(k,j,i)
853	ENDDO
854	ENDIF
855	ENDIF
856
857	!
858	!-- If required, compute prognostic equation for potential temperature
859	IF ( .NOT. neutral ) THEN
860	!
861	!-- Tendency terms for potential temperature
862	tend(:,j,i) = 0.0_wp
863	IF ( timestep_scheme(1:5) == 'runge' ) THEN
864	IF ( ws_scheme_sca ) THEN
865	CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt, &
866	flux_l_pt, diss_l_pt, i_omp_start, tn )
867	ELSE
868	CALL advec_s_pw( i, j, pt )
869	ENDIF
870	ELSE
871	CALL advec_s_up( i, j, pt )
872	ENDIF
873	CALL diffusion_s( i, j, pt, &
874	surf_def_h(0)%shf, surf_def_h(1)%shf, &
875	surf_def_h(2)%shf, &
876	surf_lsm_h%shf, surf_usm_h%shf, &
877	surf_def_v(0)%shf, surf_def_v(1)%shf, &
878	surf_def_v(2)%shf, surf_def_v(3)%shf, &
879	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
880	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
881	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
882	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
883
884	!
885	!-- Consideration of heat sources within the plant canopy
886	IF ( plant_canopy .AND. &
887	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
888	CALL pcm_tendency( i, j, 4 )
889	ENDIF
890
891	!
892	!-- Large scale advection
893	IF ( large_scale_forcing ) THEN
894	CALL ls_advec( i, j, simulated_time, 'pt' )
895	ENDIF
896
897	!
898	!-- Nudging
899	IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' )
900
901	!
902	!-- If required, compute effect of large-scale subsidence/ascent
903	IF ( large_scale_subsidence .AND. &
904	.NOT. use_subsidence_tendencies ) THEN
905	CALL subsidence( i, j, tend, pt, pt_init, 2 )
906	ENDIF
907
908	#if defined( __rrtmg )
909	!
910	!-- If required, add tendency due to radiative heating/cooling
911	IF ( radiation .AND. &
912	simulated_time > skip_time_do_radiation ) THEN
913	CALL radiation_tendency ( i, j, tend )
914	ENDIF
915	#endif
916
917	CALL module_interface_actions( i, j, 'pt-tendency' )
918	!
919	!-- Prognostic equation for potential temperature
920	DO k = nzb+1, nzt
921	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * &
922	( tsc(2) * tend(k,j,i) + &
923	tsc(3) * tpt_m(k,j,i) ) &
924	- tsc(5) &
925	* ( pt(k,j,i) - pt_init(k) ) &
926	* ( rdf_sc(k) + ptdf_x(i) &
927	+ ptdf_y(j) ) &
928	) &
929	* MERGE( 1.0_wp, 0.0_wp, &
930	BTEST( wall_flags_0(k,j,i), 0 )&
931	)
932	ENDDO
933
934	!
935	!-- Calculate tendencies for the next Runge-Kutta step
936	IF ( timestep_scheme(1:5) == 'runge' ) THEN
937	IF ( intermediate_timestep_count == 1 ) THEN
938	DO k = nzb+1, nzt
939	tpt_m(k,j,i) = tend(k,j,i)
940	ENDDO
941	ELSEIF ( intermediate_timestep_count < &
942	intermediate_timestep_count_max ) THEN
943	DO k = nzb+1, nzt
944	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
945	5.3125_wp * tpt_m(k,j,i)
946	ENDDO
947	ENDIF
948	ENDIF
949
950	ENDIF
951
952	!
953	!-- If required, compute prognostic equation for total water content
954	IF ( humidity ) THEN
955
956	!
957	!-- Tendency-terms for total water content / scalar
958	tend(:,j,i) = 0.0_wp
959	IF ( timestep_scheme(1:5) == 'runge' ) &
960	THEN
961	IF ( ws_scheme_sca ) THEN
962	CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
963	diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
964	ELSE
965	CALL advec_s_pw( i, j, q )
966	ENDIF
967	ELSE
968	CALL advec_s_up( i, j, q )
969	ENDIF
970	CALL diffusion_s( i, j, q, &
971	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
972	surf_def_h(2)%qsws, &
973	surf_lsm_h%qsws, surf_usm_h%qsws, &
974	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
975	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
976	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
977	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
978	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
979	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
980
981	!
982	!-- Sink or source of humidity due to canopy elements
983	IF ( plant_canopy ) CALL pcm_tendency( i, j, 5 )
984
985	!
986	!-- Large scale advection
987	IF ( large_scale_forcing ) THEN
988	CALL ls_advec( i, j, simulated_time, 'q' )
989	ENDIF
990
991	!
992	!-- Nudging
993	IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' )
994
995	!
996	!-- If required compute influence of large-scale subsidence/ascent
997	IF ( large_scale_subsidence .AND. &
998	.NOT. use_subsidence_tendencies ) THEN
999	CALL subsidence( i, j, tend, q, q_init, 3 )
1000	ENDIF
1001
1002	CALL module_interface_actions( i, j, 'q-tendency' )
1003
1004	!
1005	!-- Prognostic equation for total water content / scalar
1006	DO k = nzb+1, nzt
1007	q_p(k,j,i) = q(k,j,i) + ( dt_3d * &
1008	( tsc(2) * tend(k,j,i) + &
1009	tsc(3) * tq_m(k,j,i) ) &
1010	- tsc(5) * rdf_sc(k) * &
1011	( q(k,j,i) - q_init(k) ) &
1012	) &
1013	* MERGE( 1.0_wp, 0.0_wp, &
1014	BTEST( wall_flags_0(k,j,i), 0 )&
1015	)
1016	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
1017	ENDDO
1018
1019	!
1020	!-- Calculate tendencies for the next Runge-Kutta step
1021	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1022	IF ( intermediate_timestep_count == 1 ) THEN
1023	DO k = nzb+1, nzt
1024	tq_m(k,j,i) = tend(k,j,i)
1025	ENDDO
1026	ELSEIF ( intermediate_timestep_count < &
1027	intermediate_timestep_count_max ) THEN
1028	DO k = nzb+1, nzt
1029	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1030	5.3125_wp * tq_m(k,j,i)
1031	ENDDO
1032	ENDIF
1033	ENDIF
1034
1035	ENDIF
1036
1037	!
1038	!-- If required, compute prognostic equation for scalar
1039	IF ( passive_scalar ) THEN
1040	!
1041	!-- Tendency-terms for total water content / scalar
1042	tend(:,j,i) = 0.0_wp
1043	IF ( timestep_scheme(1:5) == 'runge' ) &
1044	THEN
1045	IF ( ws_scheme_sca ) THEN
1046	CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1047	diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1048	ELSE
1049	CALL advec_s_pw( i, j, s )
1050	ENDIF
1051	ELSE
1052	CALL advec_s_up( i, j, s )
1053	ENDIF
1054	CALL diffusion_s( i, j, s, &
1055	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1056	surf_def_h(2)%ssws, &
1057	surf_lsm_h%ssws, surf_usm_h%ssws, &
1058	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1059	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1060	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1061	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1062	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1063	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1064
1065	!
1066	!-- Sink or source of scalar concentration due to canopy elements
1067	IF ( plant_canopy ) CALL pcm_tendency( i, j, 7 )
1068
1069	!
1070	!-- Large scale advection, still need to be extended for scalars
1071	! IF ( large_scale_forcing ) THEN
1072	! CALL ls_advec( i, j, simulated_time, 's' )
1073	! ENDIF
1074
1075	!
1076	!-- Nudging, still need to be extended for scalars
1077	! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' )
1078
1079	!
1080	!-- If required compute influence of large-scale subsidence/ascent.
1081	!-- Note, the last argument is of no meaning in this case, as it is
1082	!-- only used in conjunction with large_scale_forcing, which is to
1083	!-- date not implemented for scalars.
1084	IF ( large_scale_subsidence .AND. &
1085	.NOT. use_subsidence_tendencies .AND. &
1086	.NOT. large_scale_forcing ) THEN
1087	CALL subsidence( i, j, tend, s, s_init, 3 )
1088	ENDIF
1089
1090	CALL module_interface_actions( i, j, 's-tendency' )
1091
1092	!
1093	!-- Prognostic equation for scalar
1094	DO k = nzb+1, nzt
1095	s_p(k,j,i) = s(k,j,i) + ( dt_3d * &
1096	( tsc(2) * tend(k,j,i) + &
1097	tsc(3) * ts_m(k,j,i) ) &
1098	- tsc(5) * rdf_sc(k) &
1099	* ( s(k,j,i) - s_init(k) )&
1100	) &
1101	* MERGE( 1.0_wp, 0.0_wp, &
1102	BTEST( wall_flags_0(k,j,i), 0 )&
1103	)
1104	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1105	ENDDO
1106
1107	!
1108	!-- Calculate tendencies for the next Runge-Kutta step
1109	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1110	IF ( intermediate_timestep_count == 1 ) THEN
1111	DO k = nzb+1, nzt
1112	ts_m(k,j,i) = tend(k,j,i)
1113	ENDDO
1114	ELSEIF ( intermediate_timestep_count < &
1115	intermediate_timestep_count_max ) THEN
1116	DO k = nzb+1, nzt
1117	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1118	5.3125_wp * ts_m(k,j,i)
1119	ENDDO
1120	ENDIF
1121	ENDIF
1122
1123	ENDIF
1124	!
1125	!-- Calculate prognostic equations for turbulence closure
1126	CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1127	!
1128	!-- Calculate prognostic equations for all other modules
1129	CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1130
1131	ENDDO ! loop over j
1132	ENDDO ! loop over i
1133	!$OMP END PARALLEL
1134
1135
1136	CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1137
1138	IF ( debug_output ) THEN
1139	WRITE( debug_string, * ) 'prognostic_equations_cache'
1140	CALL debug_message( debug_string, 'end' )
1141	ENDIF
1142
1143	END SUBROUTINE prognostic_equations_cache
1144
1145
1146	!------------------------------------------------------------------------------!
1147	! Description:
1148	! ------------
1149	!> Version for vector machines
1150	!------------------------------------------------------------------------------!
1151
1152	SUBROUTINE prognostic_equations_vector
1153
1154
1155	INTEGER(iwp) :: i !<
1156	INTEGER(iwp) :: ib !< index for aerosol size bins (salsa)
1157	INTEGER(iwp) :: ic !< index for chemical compounds (salsa)
1158	INTEGER(iwp) :: icc !< additional index for chemical compounds (salsa)
1159	INTEGER(iwp) :: ig !< index for gaseous compounds (salsa)
1160	INTEGER(iwp) :: j !<
1161	INTEGER(iwp) :: k !<
1162
1163	REAL(wp) :: sbt !<
1164
1165
1166	IF ( debug_output ) THEN
1167	WRITE( debug_string, * ) 'prognostic_equations_vector'
1168	CALL debug_message( debug_string, 'start' )
1169	ENDIF
1170	!
1171	!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1172	!-- step. The exchange of ghost points is required after this update of the
1173	!-- concentrations of aerosol number and mass
1174	IF ( salsa ) THEN
1175	IF ( time_since_reference_point >= skip_time_do_salsa ) THEN
1176	IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa ) &
1177	THEN
1178	CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1179	!$OMP PARALLEL PRIVATE (i,j)
1180	!$OMP DO
1181	!
1182	!-- Call salsa processes
1183	DO i = nxl, nxr
1184	DO j = nys, nyn
1185	CALL salsa_diagnostics( i, j )
1186	CALL salsa_driver( i, j, 3 )
1187	CALL salsa_diagnostics( i, j )
1188	ENDDO
1189	ENDDO
1190	!$OMP END PARALLEL
1191
1192	CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1193	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1194	!
1195	!-- Exchange ghost points and decycle if needed.
1196	DO ib = 1, nbins_aerosol
1197	CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
1198	CALL salsa_boundary_conds( aerosol_number(ib)%conc, &
1199	aerosol_number(ib)%init )
1200	DO ic = 1, ncomponents_mass
1201	icc = ( ic - 1 ) * nbins_aerosol + ib
1202	CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
1203	CALL salsa_boundary_conds( aerosol_mass(icc)%conc, &
1204	aerosol_mass(icc)%init )
1205	ENDDO
1206	ENDDO
1207	IF ( .NOT. salsa_gases_from_chem ) THEN
1208	DO ig = 1, ngases_salsa
1209	CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
1210	CALL salsa_boundary_conds( salsa_gas(ig)%conc, &
1211	salsa_gas(ig)%init )
1212	ENDDO
1213	ENDIF
1214	CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1215	last_salsa_time = time_since_reference_point
1216	ENDIF
1217	ENDIF
1218	ENDIF
1219
1220	!
1221	!-- Calculate non transport physics for all other modules
1222	CALL module_interface_non_transport_physics
1223	!
1224	!-- Module Inferface for exchange horiz after non_transport_physics but before
1225	!-- advection. Therefore, non_transport_physics must not run for ghost points.
1226	CALL module_interface_exchange_horiz()
1227
1228	!
1229	!-- u-velocity component
1230	CALL cpu_log( log_point(5), 'u-equation', 'start' )
1231
1232	!$ACC KERNELS PRESENT(tend)
1233	tend = 0.0_wp
1234	!$ACC END KERNELS
1235	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1236	IF ( ws_scheme_mom ) THEN
1237	CALL advec_u_ws
1238	ELSE
1239	CALL advec_u_pw
1240	ENDIF
1241	ELSE
1242	CALL advec_u_up
1243	ENDIF
1244	CALL diffusion_u
1245	CALL coriolis( 1 )
1246	IF ( sloping_surface .AND. .NOT. neutral ) THEN
1247	CALL buoyancy( pt, 1 )
1248	ENDIF
1249
1250	!
1251	!-- Drag by plant canopy
1252	IF ( plant_canopy ) CALL pcm_tendency( 1 )
1253
1254	!
1255	!-- External pressure gradient
1256	IF ( dp_external ) THEN
1257	DO i = nxlu, nxr
1258	DO j = nys, nyn
1259	DO k = dp_level_ind_b+1, nzt
1260	tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1261	ENDDO
1262	ENDDO
1263	ENDDO
1264	ENDIF
1265
1266	!
1267	!-- Nudging
1268	IF ( nudging ) CALL nudge( simulated_time, 'u' )
1269
1270	!
1271	!-- Effect of Stokes drift (in ocean mode only)
1272	IF ( stokes_force ) CALL stokes_drift_terms( 1 )
1273
1274	CALL module_interface_actions( 'u-tendency' )
1275
1276	!
1277	!-- Prognostic equation for u-velocity component
1278	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1279	!$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1280	!$ACC PRESENT(tsc(2:5)) &
1281	!$ACC PRESENT(u_p)
1282	DO i = nxlu, nxr
1283	DO j = nys, nyn
1284	DO k = nzb+1, nzt
1285	u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1286	tsc(3) * tu_m(k,j,i) ) &
1287	- tsc(5) * rdf(k) * &
1288	( u(k,j,i) - u_init(k) ) &
1289	) * MERGE( 1.0_wp, 0.0_wp, &
1290	BTEST( wall_flags_0(k,j,i), 1 ) &
1291	)
1292	ENDDO
1293	ENDDO
1294	ENDDO
1295
1296	!
1297	!-- Add turbulence generated by wave breaking (in ocean mode only)
1298	IF ( wave_breaking .AND. &
1299	intermediate_timestep_count == intermediate_timestep_count_max ) &
1300	THEN
1301	CALL wave_breaking_term( 1 )
1302	ENDIF
1303
1304	!
1305	!-- Calculate tendencies for the next Runge-Kutta step
1306	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1307	IF ( intermediate_timestep_count == 1 ) THEN
1308	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1309	!$ACC PRESENT(tend, tu_m)
1310	DO i = nxlu, nxr
1311	DO j = nys, nyn
1312	DO k = nzb+1, nzt
1313	tu_m(k,j,i) = tend(k,j,i)
1314	ENDDO
1315	ENDDO
1316	ENDDO
1317	ELSEIF ( intermediate_timestep_count < &
1318	intermediate_timestep_count_max ) THEN
1319	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1320	!$ACC PRESENT(tend, tu_m)
1321	DO i = nxlu, nxr
1322	DO j = nys, nyn
1323	DO k = nzb+1, nzt
1324	tu_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1325	+ 5.3125_wp * tu_m(k,j,i)
1326	ENDDO
1327	ENDDO
1328	ENDDO
1329	ENDIF
1330	ENDIF
1331
1332	CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1333
1334	!
1335	!-- v-velocity component
1336	CALL cpu_log( log_point(6), 'v-equation', 'start' )
1337
1338	!$ACC KERNELS PRESENT(tend)
1339	tend = 0.0_wp
1340	!$ACC END KERNELS
1341	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1342	IF ( ws_scheme_mom ) THEN
1343	CALL advec_v_ws
1344	ELSE
1345	CALL advec_v_pw
1346	END IF
1347	ELSE
1348	CALL advec_v_up
1349	ENDIF
1350	CALL diffusion_v
1351	CALL coriolis( 2 )
1352
1353	!
1354	!-- Drag by plant canopy
1355	IF ( plant_canopy ) CALL pcm_tendency( 2 )
1356
1357	!
1358	!-- External pressure gradient
1359	IF ( dp_external ) THEN
1360	DO i = nxl, nxr
1361	DO j = nysv, nyn
1362	DO k = dp_level_ind_b+1, nzt
1363	tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1364	ENDDO
1365	ENDDO
1366	ENDDO
1367	ENDIF
1368
1369	!
1370	!-- Nudging
1371	IF ( nudging ) CALL nudge( simulated_time, 'v' )
1372
1373	!
1374	!-- Effect of Stokes drift (in ocean mode only)
1375	IF ( stokes_force ) CALL stokes_drift_terms( 2 )
1376
1377	CALL module_interface_actions( 'v-tendency' )
1378
1379	!
1380	!-- Prognostic equation for v-velocity component
1381	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1382	!$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1383	!$ACC PRESENT(tsc(2:5)) &
1384	!$ACC PRESENT(v_p)
1385	DO i = nxl, nxr
1386	DO j = nysv, nyn
1387	DO k = nzb+1, nzt
1388	v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1389	tsc(3) * tv_m(k,j,i) ) &
1390	- tsc(5) * rdf(k) * &
1391	( v(k,j,i) - v_init(k) ) &
1392	) * MERGE( 1.0_wp, 0.0_wp, &
1393	BTEST( wall_flags_0(k,j,i), 2 )&
1394	)
1395	ENDDO
1396	ENDDO
1397	ENDDO
1398
1399	!
1400	!-- Add turbulence generated by wave breaking (in ocean mode only)
1401	IF ( wave_breaking .AND. &
1402	intermediate_timestep_count == intermediate_timestep_count_max ) &
1403	THEN
1404	CALL wave_breaking_term( 2 )
1405	ENDIF
1406
1407	!
1408	!-- Calculate tendencies for the next Runge-Kutta step
1409	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1410	IF ( intermediate_timestep_count == 1 ) THEN
1411	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1412	!$ACC PRESENT(tend, tv_m)
1413	DO i = nxl, nxr
1414	DO j = nysv, nyn
1415	DO k = nzb+1, nzt
1416	tv_m(k,j,i) = tend(k,j,i)
1417	ENDDO
1418	ENDDO
1419	ENDDO
1420	ELSEIF ( intermediate_timestep_count < &
1421	intermediate_timestep_count_max ) THEN
1422	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1423	!$ACC PRESENT(tend, tv_m)
1424	DO i = nxl, nxr
1425	DO j = nysv, nyn
1426	DO k = nzb+1, nzt
1427	tv_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1428	+ 5.3125_wp * tv_m(k,j,i)
1429	ENDDO
1430	ENDDO
1431	ENDDO
1432	ENDIF
1433	ENDIF
1434
1435	CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1436
1437	!
1438	!-- w-velocity component
1439	CALL cpu_log( log_point(7), 'w-equation', 'start' )
1440
1441	!$ACC KERNELS PRESENT(tend)
1442	tend = 0.0_wp
1443	!$ACC END KERNELS
1444	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1445	IF ( ws_scheme_mom ) THEN
1446	CALL advec_w_ws
1447	ELSE
1448	CALL advec_w_pw
1449	ENDIF
1450	ELSE
1451	CALL advec_w_up
1452	ENDIF
1453	CALL diffusion_w
1454	CALL coriolis( 3 )
1455
1456	IF ( .NOT. neutral ) THEN
1457	IF ( ocean_mode ) THEN
1458	CALL buoyancy( rho_ocean, 3 )
1459	ELSE
1460	IF ( .NOT. humidity ) THEN
1461	CALL buoyancy( pt, 3 )
1462	ELSE
1463	CALL buoyancy( vpt, 3 )
1464	ENDIF
1465	ENDIF
1466	ENDIF
1467
1468	!
1469	!-- Drag by plant canopy
1470	IF ( plant_canopy ) CALL pcm_tendency( 3 )
1471
1472	!
1473	!-- Effect of Stokes drift (in ocean mode only)
1474	IF ( stokes_force ) CALL stokes_drift_terms( 3 )
1475
1476	CALL module_interface_actions( 'w-tendency' )
1477
1478	!
1479	!-- Prognostic equation for w-velocity component
1480	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1481	!$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1482	!$ACC PRESENT(tsc(2:5)) &
1483	!$ACC PRESENT(w_p)
1484	DO i = nxl, nxr
1485	DO j = nys, nyn
1486	DO k = nzb+1, nzt-1
1487	w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
1488	tsc(3) * tw_m(k,j,i) ) &
1489	- tsc(5) * rdf(k) * w(k,j,i) &
1490	) * MERGE( 1.0_wp, 0.0_wp, &
1491	BTEST( wall_flags_0(k,j,i), 3 )&
1492	)
1493	ENDDO
1494	ENDDO
1495	ENDDO
1496
1497	!
1498	!-- Calculate tendencies for the next Runge-Kutta step
1499	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1500	IF ( intermediate_timestep_count == 1 ) THEN
1501	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1502	!$ACC PRESENT(tend, tw_m)
1503	DO i = nxl, nxr
1504	DO j = nys, nyn
1505	DO k = nzb+1, nzt-1
1506	tw_m(k,j,i) = tend(k,j,i)
1507	ENDDO
1508	ENDDO
1509	ENDDO
1510	ELSEIF ( intermediate_timestep_count < &
1511	intermediate_timestep_count_max ) THEN
1512	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1513	!$ACC PRESENT(tend, tw_m)
1514	DO i = nxl, nxr
1515	DO j = nys, nyn
1516	DO k = nzb+1, nzt-1
1517	tw_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1518	+ 5.3125_wp * tw_m(k,j,i)
1519	ENDDO
1520	ENDDO
1521	ENDDO
1522	ENDIF
1523	ENDIF
1524
1525	CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1526
1527
1528	!
1529	!-- If required, compute prognostic equation for potential temperature
1530	IF ( .NOT. neutral ) THEN
1531
1532	CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1533
1534	!
1535	!-- pt-tendency terms with communication
1536	sbt = tsc(2)
1537	IF ( scalar_advec == 'bc-scheme' ) THEN
1538
1539	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1540	!
1541	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1542	sbt = 1.0_wp
1543	ENDIF
1544	tend = 0.0_wp
1545	CALL advec_s_bc( pt, 'pt' )
1546
1547	ENDIF
1548
1549	!
1550	!-- pt-tendency terms with no communication
1551	IF ( scalar_advec /= 'bc-scheme' ) THEN
1552	!$ACC KERNELS PRESENT(tend)
1553	tend = 0.0_wp
1554	!$ACC END KERNELS
1555	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1556	IF ( ws_scheme_sca ) THEN
1557	CALL advec_s_ws( pt, 'pt' )
1558	ELSE
1559	CALL advec_s_pw( pt )
1560	ENDIF
1561	ELSE
1562	CALL advec_s_up( pt )
1563	ENDIF
1564	ENDIF
1565
1566	CALL diffusion_s( pt, &
1567	surf_def_h(0)%shf, surf_def_h(1)%shf, &
1568	surf_def_h(2)%shf, &
1569	surf_lsm_h%shf, surf_usm_h%shf, &
1570	surf_def_v(0)%shf, surf_def_v(1)%shf, &
1571	surf_def_v(2)%shf, surf_def_v(3)%shf, &
1572	surf_lsm_v(0)%shf, surf_lsm_v(1)%shf, &
1573	surf_lsm_v(2)%shf, surf_lsm_v(3)%shf, &
1574	surf_usm_v(0)%shf, surf_usm_v(1)%shf, &
1575	surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1576
1577	!
1578	!-- Consideration of heat sources within the plant canopy
1579	IF ( plant_canopy .AND. &
1580	(cthf /= 0.0_wp .OR. urban_surface .OR. land_surface) ) THEN
1581	CALL pcm_tendency( 4 )
1582	ENDIF
1583
1584	!
1585	!-- Large scale advection
1586	IF ( large_scale_forcing ) THEN
1587	CALL ls_advec( simulated_time, 'pt' )
1588	ENDIF
1589
1590	!
1591	!-- Nudging
1592	IF ( nudging ) CALL nudge( simulated_time, 'pt' )
1593
1594	!
1595	!-- If required compute influence of large-scale subsidence/ascent
1596	IF ( large_scale_subsidence .AND. &
1597	.NOT. use_subsidence_tendencies ) THEN
1598	CALL subsidence( tend, pt, pt_init, 2 )
1599	ENDIF
1600
1601	#if defined( __rrtmg )
1602	!
1603	!-- If required, add tendency due to radiative heating/cooling
1604	IF ( radiation .AND. &
1605	simulated_time > skip_time_do_radiation ) THEN
1606	CALL radiation_tendency ( tend )
1607	ENDIF
1608	#endif
1609
1610	CALL module_interface_actions( 'pt-tendency' )
1611
1612	!
1613	!-- Prognostic equation for potential temperature
1614	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1615	!$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1616	!$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1617	!$ACC PRESENT(tsc(3:5)) &
1618	!$ACC PRESENT(pt_p)
1619	DO i = nxl, nxr
1620	DO j = nys, nyn
1621	DO k = nzb+1, nzt
1622	pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1623	tsc(3) * tpt_m(k,j,i) ) &
1624	- tsc(5) * &
1625	( pt(k,j,i) - pt_init(k) ) *&
1626	( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1627	) &
1628	* MERGE( 1.0_wp, 0.0_wp, &
1629	BTEST( wall_flags_0(k,j,i), 0 ) &
1630	)
1631	ENDDO
1632	ENDDO
1633	ENDDO
1634	!
1635	!-- Calculate tendencies for the next Runge-Kutta step
1636	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1637	IF ( intermediate_timestep_count == 1 ) THEN
1638	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1639	!$ACC PRESENT(tend, tpt_m)
1640	DO i = nxl, nxr
1641	DO j = nys, nyn
1642	DO k = nzb+1, nzt
1643	tpt_m(k,j,i) = tend(k,j,i)
1644	ENDDO
1645	ENDDO
1646	ENDDO
1647	ELSEIF ( intermediate_timestep_count < &
1648	intermediate_timestep_count_max ) THEN
1649	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1650	!$ACC PRESENT(tend, tpt_m)
1651	DO i = nxl, nxr
1652	DO j = nys, nyn
1653	DO k = nzb+1, nzt
1654	tpt_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
1655	5.3125_wp * tpt_m(k,j,i)
1656	ENDDO
1657	ENDDO
1658	ENDDO
1659	ENDIF
1660	ENDIF
1661
1662	CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1663
1664	ENDIF
1665
1666	!
1667	!-- If required, compute prognostic equation for total water content
1668	IF ( humidity ) THEN
1669
1670	CALL cpu_log( log_point(29), 'q-equation', 'start' )
1671
1672	!
1673	!-- Scalar/q-tendency terms with communication
1674	sbt = tsc(2)
1675	IF ( scalar_advec == 'bc-scheme' ) THEN
1676
1677	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1678	!
1679	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1680	sbt = 1.0_wp
1681	ENDIF
1682	tend = 0.0_wp
1683	CALL advec_s_bc( q, 'q' )
1684
1685	ENDIF
1686
1687	!
1688	!-- Scalar/q-tendency terms with no communication
1689	IF ( scalar_advec /= 'bc-scheme' ) THEN
1690	tend = 0.0_wp
1691	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1692	IF ( ws_scheme_sca ) THEN
1693	CALL advec_s_ws( q, 'q' )
1694	ELSE
1695	CALL advec_s_pw( q )
1696	ENDIF
1697	ELSE
1698	CALL advec_s_up( q )
1699	ENDIF
1700	ENDIF
1701
1702	CALL diffusion_s( q, &
1703	surf_def_h(0)%qsws, surf_def_h(1)%qsws, &
1704	surf_def_h(2)%qsws, &
1705	surf_lsm_h%qsws, surf_usm_h%qsws, &
1706	surf_def_v(0)%qsws, surf_def_v(1)%qsws, &
1707	surf_def_v(2)%qsws, surf_def_v(3)%qsws, &
1708	surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws, &
1709	surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws, &
1710	surf_usm_v(0)%qsws, surf_usm_v(1)%qsws, &
1711	surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1712
1713	!
1714	!-- Sink or source of humidity due to canopy elements
1715	IF ( plant_canopy ) CALL pcm_tendency( 5 )
1716
1717	!
1718	!-- Large scale advection
1719	IF ( large_scale_forcing ) THEN
1720	CALL ls_advec( simulated_time, 'q' )
1721	ENDIF
1722
1723	!
1724	!-- Nudging
1725	IF ( nudging ) CALL nudge( simulated_time, 'q' )
1726
1727	!
1728	!-- If required compute influence of large-scale subsidence/ascent
1729	IF ( large_scale_subsidence .AND. &
1730	.NOT. use_subsidence_tendencies ) THEN
1731	CALL subsidence( tend, q, q_init, 3 )
1732	ENDIF
1733
1734	CALL module_interface_actions( 'q-tendency' )
1735
1736	!
1737	!-- Prognostic equation for total water content
1738	DO i = nxl, nxr
1739	DO j = nys, nyn
1740	DO k = nzb+1, nzt
1741	q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1742	tsc(3) * tq_m(k,j,i) ) &
1743	- tsc(5) * rdf_sc(k) * &
1744	( q(k,j,i) - q_init(k) ) &
1745	) * MERGE( 1.0_wp, 0.0_wp, &
1746	BTEST( wall_flags_0(k,j,i), 0 ) &
1747	)
1748	IF ( q_p(k,j,i) < 0.0_wp ) q_p(k,j,i) = 0.1_wp * q(k,j,i)
1749	ENDDO
1750	ENDDO
1751	ENDDO
1752
1753	!
1754	!-- Calculate tendencies for the next Runge-Kutta step
1755	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1756	IF ( intermediate_timestep_count == 1 ) THEN
1757	DO i = nxl, nxr
1758	DO j = nys, nyn
1759	DO k = nzb+1, nzt
1760	tq_m(k,j,i) = tend(k,j,i)
1761	ENDDO
1762	ENDDO
1763	ENDDO
1764	ELSEIF ( intermediate_timestep_count < &
1765	intermediate_timestep_count_max ) THEN
1766	DO i = nxl, nxr
1767	DO j = nys, nyn
1768	DO k = nzb+1, nzt
1769	tq_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1770	+ 5.3125_wp * tq_m(k,j,i)
1771	ENDDO
1772	ENDDO
1773	ENDDO
1774	ENDIF
1775	ENDIF
1776
1777	CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1778
1779	ENDIF
1780	!
1781	!-- If required, compute prognostic equation for scalar
1782	IF ( passive_scalar ) THEN
1783
1784	CALL cpu_log( log_point(66), 's-equation', 'start' )
1785
1786	!
1787	!-- Scalar/q-tendency terms with communication
1788	sbt = tsc(2)
1789	IF ( scalar_advec == 'bc-scheme' ) THEN
1790
1791	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1792	!
1793	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1794	sbt = 1.0_wp
1795	ENDIF
1796	tend = 0.0_wp
1797	CALL advec_s_bc( s, 's' )
1798
1799	ENDIF
1800
1801	!
1802	!-- Scalar/q-tendency terms with no communication
1803	IF ( scalar_advec /= 'bc-scheme' ) THEN
1804	tend = 0.0_wp
1805	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1806	IF ( ws_scheme_sca ) THEN
1807	CALL advec_s_ws( s, 's' )
1808	ELSE
1809	CALL advec_s_pw( s )
1810	ENDIF
1811	ELSE
1812	CALL advec_s_up( s )
1813	ENDIF
1814	ENDIF
1815
1816	CALL diffusion_s( s, &
1817	surf_def_h(0)%ssws, surf_def_h(1)%ssws, &
1818	surf_def_h(2)%ssws, &
1819	surf_lsm_h%ssws, surf_usm_h%ssws, &
1820	surf_def_v(0)%ssws, surf_def_v(1)%ssws, &
1821	surf_def_v(2)%ssws, surf_def_v(3)%ssws, &
1822	surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws, &
1823	surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws, &
1824	surf_usm_v(0)%ssws, surf_usm_v(1)%ssws, &
1825	surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1826
1827	!
1828	!-- Sink or source of humidity due to canopy elements
1829	IF ( plant_canopy ) CALL pcm_tendency( 7 )
1830
1831	!
1832	!-- Large scale advection. Not implemented for scalars so far.
1833	! IF ( large_scale_forcing ) THEN
1834	! CALL ls_advec( simulated_time, 'q' )
1835	! ENDIF
1836
1837	!
1838	!-- Nudging. Not implemented for scalars so far.
1839	! IF ( nudging ) CALL nudge( simulated_time, 'q' )
1840
1841	!
1842	!-- If required compute influence of large-scale subsidence/ascent.
1843	!-- Not implemented for scalars so far.
1844	IF ( large_scale_subsidence .AND. &
1845	.NOT. use_subsidence_tendencies .AND. &
1846	.NOT. large_scale_forcing ) THEN
1847	CALL subsidence( tend, s, s_init, 3 )
1848	ENDIF
1849
1850	CALL module_interface_actions( 's-tendency' )
1851
1852	!
1853	!-- Prognostic equation for total water content
1854	DO i = nxl, nxr
1855	DO j = nys, nyn
1856	DO k = nzb+1, nzt
1857	s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
1858	tsc(3) * ts_m(k,j,i) ) &
1859	- tsc(5) * rdf_sc(k) * &
1860	( s(k,j,i) - s_init(k) ) &
1861	) &
1862	* MERGE( 1.0_wp, 0.0_wp, &
1863	BTEST( wall_flags_0(k,j,i), 0 ) &
1864	)
1865	IF ( s_p(k,j,i) < 0.0_wp ) s_p(k,j,i) = 0.1_wp * s(k,j,i)
1866	ENDDO
1867	ENDDO
1868	ENDDO
1869
1870	!
1871	!-- Calculate tendencies for the next Runge-Kutta step
1872	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1873	IF ( intermediate_timestep_count == 1 ) THEN
1874	DO i = nxl, nxr
1875	DO j = nys, nyn
1876	DO k = nzb+1, nzt
1877	ts_m(k,j,i) = tend(k,j,i)
1878	ENDDO
1879	ENDDO
1880	ENDDO
1881	ELSEIF ( intermediate_timestep_count < &
1882	intermediate_timestep_count_max ) THEN
1883	DO i = nxl, nxr
1884	DO j = nys, nyn
1885	DO k = nzb+1, nzt
1886	ts_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
1887	+ 5.3125_wp * ts_m(k,j,i)
1888	ENDDO
1889	ENDDO
1890	ENDDO
1891	ENDIF
1892	ENDIF
1893
1894	CALL cpu_log( log_point(66), 's-equation', 'stop' )
1895
1896	ENDIF
1897
1898	!
1899	!-- Calculate prognostic equations for turbulence closure
1900	CALL tcm_prognostic_equations()
1901	!
1902	!-- Calculate prognostic equations for all other modules
1903	CALL module_interface_prognostic_equations()
1904
1905	IF ( debug_output ) THEN
1906	WRITE( debug_string, * ) 'prognostic_equations_vector'
1907	CALL debug_message( debug_string, 'end' )
1908	ENDIF
1909
1910	END SUBROUTINE prognostic_equations_vector
1911
1912
1913	END MODULE prognostic_equations_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |