Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |