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

Last change on this file since 3864 was 3864, checked in by monakurppa, 2 years ago

major changes in salsa: data input, format and performance

  • Time-dependent emissions enabled: lod=1 for yearly PM emissions that are normalised depending on the time, and lod=2 for preprocessed emissions (similar to the chemistry module).
  • Additionally, 'uniform' emissions allowed. This emission is set constant on all horisontal upward facing surfaces and it is created based on parameters surface_aerosol_flux, aerosol_flux_dpg/sigmag/mass_fracs_a/mass_fracs_b.
  • All emissions are now implemented as surface fluxes! No 3D sources anymore.
  • Update the emission information by calling salsa_emission_update if skip_time_do_salsa >= time_since_reference_point and next_aero_emission_update <= time_since_reference_point
  • Aerosol background concentrations read from PIDS_DYNAMIC. The vertical grid must match the one applied in the model.
  • Gas emissions and background concentrations can be also read in in salsa_mod if the chemistry module is not applied.
  • In deposition, information on the land use type can be now imported from the land use model
  • Use SI units in PARIN, i.e. n_lognorm given in #/m3 and dpg in metres.
  • Apply 100 character line limit
  • Change all variable names from capital to lowercase letter
  • Change real exponents to integer if possible. If not, precalculate the value of exponent
  • Rename in1a to start_subrange_1a, fn2a to end_subrange_1a etc.
  • Rename nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and ngast --> ngases_salsa
  • Rename ibc to index_bc, idu to index_du etc.
  • Renamed loop indices b, c and sg to ib, ic and ig
  • run_salsa subroutine removed
  • Corrected a bud in salsa_driver: falsely applied ino instead of inh
  • Call salsa_tendency within salsa_prognostic_equations which is called in module_interface_mod instead of prognostic_equations_mod
  • Removed tailing white spaces and unused variables
  • Change error message to start by PA instead of SA
  • Property svn:keywords set to Id
File size: 100.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 3864 2019-04-05 09:01:56Z monakurppa $
27! Modifications made for salsa:
28! - salsa_prognostic_equations moved to salsa_mod (and the call to
29!   module_interface_mod)
30! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
31!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
32!
33! 3840 2019-03-29 10:35:52Z knoop
34! added USE chem_gasphase_mod for nvar, nspec and spc_names
35!
36! 3820 2019-03-27 11:53:41Z forkel
37! renamed do_depo to deposition_dry (ecc)
38!
39! 3797 2019-03-15 11:15:38Z forkel
40! Call chem_integegrate in OpenMP loop   (ketelsen)
41!
42!
43! 3771 2019-02-28 12:19:33Z raasch
44! preprocessor directivs fro rrtmg added
45!
46! 3761 2019-02-25 15:31:42Z raasch
47! unused variable removed
48!
49! 3719 2019-02-06 13:10:18Z kanani
50! Cleaned up chemistry cpu measurements
51!
52! 3684 2019-01-20 20:20:58Z knoop
53! OpenACC port for SPEC
54!
55! 3589 2018-11-30 15:09:51Z suehring
56! Move the control parameter "salsa" from salsa_mod to control_parameters
57! (M. Kurppa)
58!
59! 3582 2018-11-29 19:16:36Z suehring
60! Implementation of a new aerosol module salsa.
61! Remove cpu-logs from i,j loop in cache version.
62!
63! 3458 2018-10-30 14:51:23Z kanani
64! remove duplicate USE chem_modules
65! from chemistry branch r3443, banzhafs, basit:
66! chem_depo call introduced
67! code added for decycling chemistry
68!
69! 3386 2018-10-19 16:28:22Z gronemeier
70! Renamed tcm_prognostic to tcm_prognostic_equations
71!
72! 3355 2018-10-16 14:03:34Z knoop
73! (from branch resler)
74! Fix for chemistry call
75!
76! 3302 2018-10-03 02:39:40Z raasch
77! Stokes drift + wave breaking term added
78!
79! 3298 2018-10-02 12:21:11Z kanani
80! Code added for decycling chemistry (basit)
81!
82! 3294 2018-10-01 02:37:10Z raasch
83! changes concerning modularization of ocean option
84!
85! 3274 2018-09-24 15:42:55Z knoop
86! Modularization of all bulk cloud physics code components
87!
88! 3241 2018-09-12 15:02:00Z raasch
89! omp_get_thread_num now declared in omp directive
90!
91! 3183 2018-07-27 14:25:55Z suehring
92! Remove unused variables from USE statements
93!
94! 3182 2018-07-27 13:36:03Z suehring
95! Revise recent bugfix for nesting
96!
97! 3021 2018-05-16 08:14:20Z maronga
98! Bugfix in IF clause for nesting
99!
100! 3014 2018-05-09 08:42:38Z maronga
101! Fixed a bug in the IF condition to call pcm_tendency in case of
102! potential temperature
103!
104! 2815 2018-02-19 11:29:57Z kanani
105! Rename chem_tendency to chem_prognostic_equations,
106! implement vector version for air chemistry
107!
108! 2766 2018-01-22 17:17:47Z kanani
109! Removed preprocessor directive __chem
110!
111! 2746 2018-01-15 12:06:04Z suehring
112! Move flag plant canopy to modules
113!
114! 2719 2018-01-02 09:02:06Z maronga
115! Bugfix for last change.
116!
117! 2718 2018-01-02 08:49:38Z maronga
118! Corrected "Former revisions" section
119!
120! 2696 2017-12-14 17:12:51Z kanani
121! - Change in file header (GPL part)
122! - Moved TKE equation to tcm_prognostic (TG)
123! - Added switch for chemical reactions (RF, FK)
124! - Implementation of chemistry module (RF, BK, FK)
125!
126! 2563 2017-10-19 15:36:10Z Giersch
127! Variable wind_turbine moved to module control_parameters
128!
129! 2320 2017-07-21 12:47:43Z suehring
130! Modularize large-scale forcing and nudging
131!
132! 2292 2017-06-20 09:51:42Z schwenkel
133! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
134! includes two more prognostic equations for cloud drop concentration (nc) 
135! and cloud water content (qc).
136!
137! 2261 2017-06-08 14:25:57Z raasch
138! bugfix for r2232: openmp directives removed
139!
140! 2233 2017-05-30 18:08:54Z suehring
141!
142! 2232 2017-05-30 17:47:52Z suehring
143! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
144! which is realized directly in diffusion_s now.
145!
146! 2192 2017-03-22 04:14:10Z raasch
147! Bugfix for misplaced and missing openMP directives from r2155
148!
149! 2155 2017-02-21 09:57:40Z hoffmann
150! Bugfix in the calculation of microphysical quantities on ghost points.
151!
152! 2118 2017-01-17 16:38:49Z raasch
153! OpenACC version of subroutine removed
154!
155! 2031 2016-10-21 15:11:58Z knoop
156! renamed variable rho to rho_ocean
157!
158! 2011 2016-09-19 17:29:57Z kanani
159! Flag urban_surface is now defined in module control_parameters.
160!
161! 2007 2016-08-24 15:47:17Z kanani
162! Added pt tendency calculation based on energy balance at urban surfaces
163! (new urban surface model)
164!
165! 2000 2016-08-20 18:09:15Z knoop
166! Forced header and separation lines into 80 columns
167!
168! 1976 2016-07-27 13:28:04Z maronga
169! Simplied calls to radiation model
170!
171! 1960 2016-07-12 16:34:24Z suehring
172! Separate humidity and passive scalar
173!
174! 1914 2016-05-26 14:44:07Z witha
175! Added calls for wind turbine model
176!
177! 1873 2016-04-18 14:50:06Z maronga
178! Module renamed (removed _mod)
179!
180! 1850 2016-04-08 13:29:27Z maronga
181! Module renamed
182!
183! 1826 2016-04-07 12:01:39Z maronga
184! Renamed canopy model calls.
185!
186! 1822 2016-04-07 07:49:42Z hoffmann
187! Kessler microphysics scheme moved to microphysics.
188!
189! 1757 2016-02-22 15:49:32Z maronga
190!
191! 1691 2015-10-26 16:17:44Z maronga
192! Added optional model spin-up without radiation / land surface model calls.
193! Formatting corrections.
194!
195! 1682 2015-10-07 23:56:08Z knoop
196! Code annotations made doxygen readable
197!
198! 1585 2015-04-30 07:05:52Z maronga
199! Added call for temperature tendency calculation due to radiative flux divergence
200!
201! 1517 2015-01-07 19:12:25Z hoffmann
202! advec_s_bc_mod addded, since advec_s_bc is now a module
203!
204! 1484 2014-10-21 10:53:05Z kanani
205! Changes due to new module structure of the plant canopy model:
206! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
207! Removed double-listing of use_upstream_for_tke in ONLY-list of module
208! control_parameters
209!
210! 1409 2014-05-23 12:11:32Z suehring
211! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
212! This ensures that left-hand side fluxes are also calculated for nxl in that
213! case, even though the solution at nxl is overwritten in boundary_conds()
214!
215! 1398 2014-05-07 11:15:00Z heinze
216! Rayleigh-damping for horizontal velocity components changed: instead of damping
217! against ug and vg, damping against u_init and v_init is used to allow for a
218! homogenized treatment in case of nudging
219!
220! 1380 2014-04-28 12:40:45Z heinze
221! Change order of calls for scalar prognostic quantities:
222! ls_advec -> nudging -> subsidence since initial profiles
223!
224! 1374 2014-04-25 12:55:07Z raasch
225! missing variables added to ONLY lists
226!
227! 1365 2014-04-22 15:03:56Z boeske
228! Calls of ls_advec for large scale advection added,
229! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
230! new argument ls_index added to the calls of subsidence
231! +ls_index
232!
233! 1361 2014-04-16 15:17:48Z hoffmann
234! Two-moment microphysics moved to the start of prognostic equations. This makes
235! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
236! Additionally, it is allowed to call the microphysics just once during the time
237! step (not at each sub-time step).
238!
239! Two-moment cloud physics added for vector and accelerator optimization.
240!
241! 1353 2014-04-08 15:21:23Z heinze
242! REAL constants provided with KIND-attribute
243!
244! 1337 2014-03-25 15:11:48Z heinze
245! Bugfix: REAL constants provided with KIND-attribute
246!
247! 1332 2014-03-25 11:59:43Z suehring
248! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
249!
250! 1330 2014-03-24 17:29:32Z suehring
251! In case of SGS-particle velocity advection of TKE is also allowed with
252! dissipative 5th-order scheme.
253!
254! 1320 2014-03-20 08:40:49Z raasch
255! ONLY-attribute added to USE-statements,
256! kind-parameters added to all INTEGER and REAL declaration statements,
257! kinds are defined in new module kinds,
258! old module precision_kind is removed,
259! revision history before 2012 removed,
260! comment fields (!:) to be used for variable explanations added to
261! all variable declaration statements
262!
263! 1318 2014-03-17 13:35:16Z raasch
264! module interfaces removed
265!
266! 1257 2013-11-08 15:18:40Z raasch
267! openacc loop vector clauses removed, independent clauses added
268!
269! 1246 2013-11-01 08:59:45Z heinze
270! enable nudging also for accelerator version
271!
272! 1241 2013-10-30 11:36:58Z heinze
273! usage of nudging enabled (so far not implemented for accelerator version)
274!
275! 1179 2013-06-14 05:57:58Z raasch
276! two arguments removed from routine buoyancy, ref_state updated on device
277!
278! 1128 2013-04-12 06:19:32Z raasch
279! those parts requiring global communication moved to time_integration,
280! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
281! j_north
282!
283! 1115 2013-03-26 18:16:16Z hoffmann
284! optimized cloud physics: calculation of microphysical tendencies transfered
285! to microphysics.f90; qr and nr are only calculated if precipitation is required
286!
287! 1111 2013-03-08 23:54:10Z raasch
288! update directives for prognostic quantities removed
289!
290! 1106 2013-03-04 05:31:38Z raasch
291! small changes in code formatting
292!
293! 1092 2013-02-02 11:24:22Z raasch
294! unused variables removed
295!
296! 1053 2012-11-13 17:11:03Z hoffmann
297! implementation of two new prognostic equations for rain drop concentration (nr)
298! and rain water content (qr)
299!
300! currently, only available for cache loop optimization
301!
302! 1036 2012-10-22 13:43:42Z raasch
303! code put under GPL (PALM 3.9)
304!
305! 1019 2012-09-28 06:46:45Z raasch
306! non-optimized version of prognostic_equations removed
307!
308! 1015 2012-09-27 09:23:24Z raasch
309! new branch prognostic_equations_acc
310! OpenACC statements added + code changes required for GPU optimization
311!
312! 1001 2012-09-13 14:08:46Z raasch
313! all actions concerning leapfrog- and upstream-spline-scheme removed
314!
315! 978 2012-08-09 08:28:32Z fricke
316! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
317! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
318! boundary in case of non-cyclic lateral boundaries
319! Bugfix: first thread index changes for WS-scheme at the inflow
320!
321! 940 2012-07-09 14:31:00Z raasch
322! temperature equation can be switched off
323!
324! Revision 1.1  2000/04/13 14:56:27  schroeter
325! Initial revision
326!
327!
328! Description:
329! ------------
330!> Solving the prognostic equations.
331!------------------------------------------------------------------------------!
332 MODULE prognostic_equations_mod
333
334
335    USE advec_s_bc_mod,                                                        &
336        ONLY:  advec_s_bc
337
338    USE advec_s_pw_mod,                                                        &
339        ONLY:  advec_s_pw
340
341    USE advec_s_up_mod,                                                        &
342        ONLY:  advec_s_up
343
344    USE advec_u_pw_mod,                                                        &
345        ONLY:  advec_u_pw
346
347    USE advec_u_up_mod,                                                        &
348        ONLY:  advec_u_up
349
350    USE advec_v_pw_mod,                                                        &
351        ONLY:  advec_v_pw
352
353    USE advec_v_up_mod,                                                        &
354        ONLY:  advec_v_up
355
356    USE advec_w_pw_mod,                                                        &
357        ONLY:  advec_w_pw
358
359    USE advec_w_up_mod,                                                        &
360        ONLY:  advec_w_up
361
362    USE advec_ws,                                                              &
363        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
364
365    USE arrays_3d,                                                             &
366        ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qc, &
367               diss_l_qr, diss_l_s, diss_l_sa, diss_s_e, diss_s_nc, diss_s_nr, &
368               diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, diss_s_sa, &
369               e, e_p, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,    &
370               flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, flux_l_nc, &
371               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
372               flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
373               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
374               ref_state, rho_ocean, s, s_init, s_p, tend, te_m, tnc_m,        &
375               tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tu_m, tv_m, tw_m, u,    &
376               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
377
378    USE bulk_cloud_model_mod,                                                  &
379        ONLY:  call_microphysics_at_all_substeps, bulk_cloud_model,            &
380               bcm_actions, microphysics_sat_adjust,                           &
381               microphysics_morrison, microphysics_seifert
382
383    USE buoyancy_mod,                                                          &
384        ONLY:  buoyancy
385
386    USE chem_modules,                                                          &
387        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name,           &
388               deposition_dry
389
390    USE chem_gasphase_mod,                                                     &
391        ONLY:  nspec, nvar, spc_names
392
393    USE chem_photolysis_mod,                                                   &
394        ONLY:  photolysis_control
395
396    USE chemistry_model_mod,                                                   &
397        ONLY:  chem_boundary_conds, chem_depo, chem_integrate,                 &
398               chem_prognostic_equations, chem_species
399
400    USE control_parameters,                                                    &
401        ONLY:  air_chemistry, constant_diffusion,                              &
402               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
403               humidity, intermediate_timestep_count,                          &
404               intermediate_timestep_count_max, large_scale_forcing,           &
405               large_scale_subsidence, neutral, nudging,                       &
406               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
407               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
408               timestep_scheme, tsc, use_subsidence_tendencies,                &
409               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
410               ws_scheme_sca, urban_surface, land_surface,                     &
411               time_since_reference_point, salsa
412
413    USE coriolis_mod,                                                          &
414        ONLY:  coriolis
415
416    USE cpulog,                                                                &
417        ONLY:  cpu_log, log_point, log_point_s
418
419    USE diffusion_s_mod,                                                       &
420        ONLY:  diffusion_s
421
422    USE diffusion_u_mod,                                                       &
423        ONLY:  diffusion_u
424
425    USE diffusion_v_mod,                                                       &
426        ONLY:  diffusion_v
427
428    USE diffusion_w_mod,                                                       &
429        ONLY:  diffusion_w
430
431    USE indices,                                                               &
432        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
433               nzb, nzt, wall_flags_0
434
435    USE kinds
436
437    USE lsf_nudging_mod,                                                       &
438        ONLY:  ls_advec, nudge
439
440    USE module_interface,                                                      &
441        ONLY:  module_interface_actions, &
442               module_interface_prognostic_equations
443
444    USE ocean_mod,                                                             &
445        ONLY:  stokes_drift_terms, stokes_force,   &
446               wave_breaking, wave_breaking_term
447
448    USE plant_canopy_model_mod,                                                &
449        ONLY:  cthf, pcm_tendency
450
451#if defined( __rrtmg )
452    USE radiation_model_mod,                                                   &
453        ONLY:  radiation, radiation_tendency,                                  &
454               skip_time_do_radiation
455#endif
456
457    USE salsa_mod,                                                             &
458        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time,        &
459               nbins_aerosol, ncomponents_mass, ngases_salsa,                  &
460               salsa_boundary_conds, salsa_diagnostics, salsa_driver,          &
461               salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
462
463    USE salsa_util_mod,                                                        &
464        ONLY:  sums_salsa_ws_l
465
466    USE statistics,                                                            &
467        ONLY:  hom
468
469    USE subsidence_mod,                                                        &
470        ONLY:  subsidence
471
472    USE surface_mod,                                                           &
473        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
474                surf_usm_v
475
476    USE turbulence_closure_mod,                                                &
477        ONLY:  tcm_prognostic_equations
478
479    USE wind_turbine_model_mod,                                                &
480        ONLY:  wtm_tendencies
481
482
483    PRIVATE
484    PUBLIC prognostic_equations_cache, prognostic_equations_vector
485
486    INTERFACE prognostic_equations_cache
487       MODULE PROCEDURE prognostic_equations_cache
488    END INTERFACE prognostic_equations_cache
489
490    INTERFACE prognostic_equations_vector
491       MODULE PROCEDURE prognostic_equations_vector
492    END INTERFACE prognostic_equations_vector
493
494
495 CONTAINS
496
497
498!------------------------------------------------------------------------------!
499! Description:
500! ------------
501!> Version with one optimized loop over all equations. It is only allowed to
502!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
503!>
504!> Here the calls of most subroutines are embedded in two DO loops over i and j,
505!> so communication between CPUs is not allowed (does not make sense) within
506!> these loops.
507!>
508!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
509!------------------------------------------------------------------------------!
510
511 SUBROUTINE prognostic_equations_cache
512
513
514    IMPLICIT NONE
515
516    INTEGER(iwp) ::  i                   !<
517    INTEGER(iwp) ::  i_omp_start         !<
518    INTEGER(iwp) ::  ib                  !< index for aerosol size bins (salsa)
519    INTEGER(iwp) ::  ic                  !< index for chemical compounds (salsa)
520    INTEGER(iwp) ::  icc                 !< additional index for chemical compounds (salsa)
521    INTEGER(iwp) ::  ig                  !< index for gaseous compounds (salsa)
522    INTEGER(iwp) ::  j                   !<
523    INTEGER(iwp) ::  k                   !<
524!$  INTEGER(iwp) ::  omp_get_thread_num  !<
525    INTEGER(iwp) ::  tn = 0              !<
526
527    LOGICAL      ::  loop_start          !<
528    INTEGER(iwp) ::  lsp
529    INTEGER(iwp) ::  lsp_usr             !< lsp running index for chem spcs
530
531
532!
533!-- Time measurement can only be performed for the whole set of equations
534    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
535
536!
537!-- Calculation of chemical reactions. This is done outside of main loop,
538!-- since exchange of ghost points is required after this update of the
539!-- concentrations of chemical species                                   
540    IF ( air_chemistry )  THEN
541       lsp_usr = 1
542!
543!--    Chemical reactions and deposition
544       IF ( chem_gasphase_on ) THEN
545!
546!--       If required, calculate photolysis frequencies -
547!--       UNFINISHED: Why not before the intermediate timestep loop?
548          IF ( intermediate_timestep_count ==  1 )  THEN
549             CALL photolysis_control
550          ENDIF
551
552          IF ( intermediate_timestep_count == 1 .OR.                        &
553             call_chem_at_all_substeps )  THEN
554
555             CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
556             !$OMP PARALLEL PRIVATE (i,j)
557             !$OMP DO schedule(static,1)
558             DO  i = nxl, nxr
559                DO  j = nys, nyn
560                   CALL chem_integrate (i,j)
561                ENDDO
562             ENDDO
563             !$OMP END PARALLEL
564             CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
565
566             IF ( deposition_dry )  THEN
567                CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
568                DO  i = nxl, nxr
569                   DO  j = nys, nyn
570                      CALL chem_depo(i,j)
571                   ENDDO
572                ENDDO
573                CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
574             ENDIF
575          ENDIF
576       ENDIF
577!
578!--    Loop over chemical species       
579       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
580       DO  lsp = 1, nspec
581          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
582          lsp_usr = 1 
583          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
584             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
585
586                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
587                                          chem_species(lsp)%conc_pr_init ) 
588             
589             ENDIF
590             lsp_usr = lsp_usr +1
591          ENDDO
592
593
594       ENDDO
595       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
596
597    ENDIF
598!
599!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
600!-- step. The exchange of ghost points is required after this update of the
601!-- concentrations of aerosol number and mass
602    IF ( salsa )  THEN
603
604       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
605          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
606          THEN
607             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
608             !$OMP PARALLEL PRIVATE (i,j,ib,ic,icc,ig)
609             !$OMP DO
610!
611!--          Call salsa processes
612             DO  i = nxl, nxr
613                DO  j = nys, nyn
614                   CALL salsa_diagnostics( i, j )
615                   CALL salsa_driver( i, j, 3 )
616                   CALL salsa_diagnostics( i, j )
617                ENDDO
618             ENDDO
619
620             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
621             
622             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
623!
624!--          Exchange ghost points and decycle if needed.
625             DO  ib = 1, nbins_aerosol
626                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
627                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
628                                           aerosol_number(ib)%init )
629                DO  ic = 1, ncomponents_mass
630                   icc = ( ic - 1 ) * nbins_aerosol + ib
631                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
632                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
633                                              aerosol_mass(icc)%init )
634                ENDDO
635             ENDDO
636
637             IF ( .NOT. salsa_gases_from_chem )  THEN
638                DO  ig = 1, ngases_salsa
639                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
640                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
641                                              salsa_gas(ig)%init )
642                ENDDO
643             ENDIF
644             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
645             
646             !$OMP END PARALLEL
647             last_salsa_time = time_since_reference_point
648
649          ENDIF
650
651       ENDIF
652
653    ENDIF
654
655!
656!-- If required, calculate cloud microphysics
657    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
658         ( intermediate_timestep_count == 1  .OR.                              &
659           call_microphysics_at_all_substeps ) )                               &
660    THEN
661       !$OMP PARALLEL PRIVATE (i,j)
662       !$OMP DO
663       DO  i = nxlg, nxrg
664          DO  j = nysg, nyng
665             CALL bcm_actions( i, j )
666           ENDDO
667       ENDDO
668       !$OMP END PARALLEL
669    ENDIF
670
671!
672!-- Loop over all prognostic equations
673!-- b, c ang g added for SALSA
674    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
675
676    !$  tn = omp_get_thread_num()
677    loop_start = .TRUE.
678
679    !$OMP DO
680    DO  i = nxl, nxr
681
682!
683!--    Store the first loop index. It differs for each thread and is required
684!--    later in advec_ws
685       IF ( loop_start )  THEN
686          loop_start  = .FALSE.
687          i_omp_start = i
688       ENDIF
689
690       DO  j = nys, nyn
691!
692!--       Tendency terms for u-velocity component. Please note, in case of
693!--       non-cyclic boundary conditions the grid point i=0 is excluded from
694!--       the prognostic equations for the u-component.   
695          IF ( i >= nxlu )  THEN
696
697             tend(:,j,i) = 0.0_wp
698             IF ( timestep_scheme(1:5) == 'runge' )  THEN
699                IF ( ws_scheme_mom )  THEN
700                   CALL advec_u_ws( i, j, i_omp_start, tn )
701                ELSE
702                   CALL advec_u_pw( i, j )
703                ENDIF
704             ELSE
705                CALL advec_u_up( i, j )
706             ENDIF
707             CALL diffusion_u( i, j )
708             CALL coriolis( i, j, 1 )
709             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
710                CALL buoyancy( i, j, pt, 1 )
711             ENDIF
712
713!
714!--          Drag by plant canopy
715             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
716
717!
718!--          External pressure gradient
719             IF ( dp_external )  THEN
720                DO  k = dp_level_ind_b+1, nzt
721                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
722                ENDDO
723             ENDIF
724
725!
726!--          Nudging
727             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
728
729!
730!--          Effect of Stokes drift (in ocean mode only)
731             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
732
733!
734!--          Forces by wind turbines
735             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
736
737             CALL module_interface_actions( i, j, 'u-tendency' )
738!
739!--          Prognostic equation for u-velocity component
740             DO  k = nzb+1, nzt
741
742                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
743                                            ( tsc(2) * tend(k,j,i) +            &
744                                              tsc(3) * tu_m(k,j,i) )            &
745                                            - tsc(5) * rdf(k)                   &
746                                                     * ( u(k,j,i) - u_init(k) ) &
747                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
748                                                   BTEST( wall_flags_0(k,j,i), 1 )&
749                                                 )
750             ENDDO
751
752!
753!--          Add turbulence generated by wave breaking (in ocean mode only)
754             IF ( wave_breaking  .AND.                                         &
755               intermediate_timestep_count == intermediate_timestep_count_max )&
756             THEN
757                CALL wave_breaking_term( i, j, 1 )
758             ENDIF
759
760!
761!--          Calculate tendencies for the next Runge-Kutta step
762             IF ( timestep_scheme(1:5) == 'runge' )  THEN
763                IF ( intermediate_timestep_count == 1 )  THEN
764                   DO  k = nzb+1, nzt
765                      tu_m(k,j,i) = tend(k,j,i)
766                   ENDDO
767                ELSEIF ( intermediate_timestep_count < &
768                         intermediate_timestep_count_max )  THEN
769                   DO  k = nzb+1, nzt
770                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
771                                     + 5.3125_wp * tu_m(k,j,i)
772                   ENDDO
773                ENDIF
774             ENDIF
775
776          ENDIF
777!
778!--       Tendency terms for v-velocity component. Please note, in case of
779!--       non-cyclic boundary conditions the grid point j=0 is excluded from
780!--       the prognostic equations for the v-component. !--       
781          IF ( j >= nysv )  THEN
782
783             tend(:,j,i) = 0.0_wp
784             IF ( timestep_scheme(1:5) == 'runge' )  THEN
785                IF ( ws_scheme_mom )  THEN
786                    CALL advec_v_ws( i, j, i_omp_start, tn )
787                ELSE
788                    CALL advec_v_pw( i, j )
789                ENDIF
790             ELSE
791                CALL advec_v_up( i, j )
792             ENDIF
793             CALL diffusion_v( i, j )
794             CALL coriolis( i, j, 2 )
795
796!
797!--          Drag by plant canopy
798             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
799
800!
801!--          External pressure gradient
802             IF ( dp_external )  THEN
803                DO  k = dp_level_ind_b+1, nzt
804                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
805                ENDDO
806             ENDIF
807
808!
809!--          Nudging
810             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
811
812!
813!--          Effect of Stokes drift (in ocean mode only)
814             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
815
816!
817!--          Forces by wind turbines
818             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
819
820             CALL module_interface_actions( i, j, 'v-tendency' )
821!
822!--          Prognostic equation for v-velocity component
823             DO  k = nzb+1, nzt
824                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
825                                            ( tsc(2) * tend(k,j,i) +           &
826                                              tsc(3) * tv_m(k,j,i) )           &
827                                            - tsc(5) * rdf(k)                  &
828                                                     * ( v(k,j,i) - v_init(k) )&
829                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
830                                                   BTEST( wall_flags_0(k,j,i), 2 )&
831                                                 )
832             ENDDO
833
834!
835!--          Add turbulence generated by wave breaking (in ocean mode only)
836             IF ( wave_breaking  .AND.                                         &
837               intermediate_timestep_count == intermediate_timestep_count_max )&
838             THEN
839                CALL wave_breaking_term( i, j, 2 )
840             ENDIF
841
842!
843!--          Calculate tendencies for the next Runge-Kutta step
844             IF ( timestep_scheme(1:5) == 'runge' )  THEN
845                IF ( intermediate_timestep_count == 1 )  THEN
846                   DO  k = nzb+1, nzt
847                      tv_m(k,j,i) = tend(k,j,i)
848                   ENDDO
849                ELSEIF ( intermediate_timestep_count < &
850                         intermediate_timestep_count_max )  THEN
851                   DO  k = nzb+1, nzt
852                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
853                                     + 5.3125_wp * tv_m(k,j,i)
854                   ENDDO
855                ENDIF
856             ENDIF
857
858          ENDIF
859
860!
861!--       Tendency terms for w-velocity component
862          tend(:,j,i) = 0.0_wp
863          IF ( timestep_scheme(1:5) == 'runge' )  THEN
864             IF ( ws_scheme_mom )  THEN
865                CALL advec_w_ws( i, j, i_omp_start, tn )
866             ELSE
867                CALL advec_w_pw( i, j )
868             END IF
869          ELSE
870             CALL advec_w_up( i, j )
871          ENDIF
872          CALL diffusion_w( i, j )
873          CALL coriolis( i, j, 3 )
874
875          IF ( .NOT. neutral )  THEN
876             IF ( ocean_mode )  THEN
877                CALL buoyancy( i, j, rho_ocean, 3 )
878             ELSE
879                IF ( .NOT. humidity )  THEN
880                   CALL buoyancy( i, j, pt, 3 )
881                ELSE
882                   CALL buoyancy( i, j, vpt, 3 )
883                ENDIF
884             ENDIF
885          ENDIF
886
887!
888!--       Drag by plant canopy
889          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
890
891!
892!--       Effect of Stokes drift (in ocean mode only)
893          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
894
895!
896!--       Forces by wind turbines
897          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
898
899          CALL module_interface_actions( i, j, 'w-tendency' )
900!
901!--       Prognostic equation for w-velocity component
902          DO  k = nzb+1, nzt-1
903             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
904                                             ( tsc(2) * tend(k,j,i) +          &
905                                               tsc(3) * tw_m(k,j,i) )          &
906                                             - tsc(5) * rdf(k) * w(k,j,i)      &
907                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
908                                                BTEST( wall_flags_0(k,j,i), 3 )&
909                                              )
910          ENDDO
911
912!
913!--       Calculate tendencies for the next Runge-Kutta step
914          IF ( timestep_scheme(1:5) == 'runge' )  THEN
915             IF ( intermediate_timestep_count == 1 )  THEN
916                DO  k = nzb+1, nzt-1
917                   tw_m(k,j,i) = tend(k,j,i)
918                ENDDO
919             ELSEIF ( intermediate_timestep_count < &
920                      intermediate_timestep_count_max )  THEN
921                DO  k = nzb+1, nzt-1
922                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
923                                  + 5.3125_wp * tw_m(k,j,i)
924                ENDDO
925             ENDIF
926          ENDIF
927
928!
929!--       If required, compute prognostic equation for potential temperature
930          IF ( .NOT. neutral )  THEN
931!
932!--          Tendency terms for potential temperature
933             tend(:,j,i) = 0.0_wp
934             IF ( timestep_scheme(1:5) == 'runge' )  THEN
935                   IF ( ws_scheme_sca )  THEN
936                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
937                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
938                   ELSE
939                      CALL advec_s_pw( i, j, pt )
940                   ENDIF
941             ELSE
942                CALL advec_s_up( i, j, pt )
943             ENDIF
944             CALL diffusion_s( i, j, pt,                                       &
945                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
946                               surf_def_h(2)%shf,                              &
947                               surf_lsm_h%shf,    surf_usm_h%shf,              &
948                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
949                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
950                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
951                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
952                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
953                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
954
955!
956!--          Consideration of heat sources within the plant canopy
957             IF ( plant_canopy  .AND.                                          &
958                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
959                CALL pcm_tendency( i, j, 4 )
960             ENDIF
961
962!
963!--          Large scale advection
964             IF ( large_scale_forcing )  THEN
965                CALL ls_advec( i, j, simulated_time, 'pt' )
966             ENDIF
967
968!
969!--          Nudging
970             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
971
972!
973!--          If required, compute effect of large-scale subsidence/ascent
974             IF ( large_scale_subsidence  .AND.                                &
975                  .NOT. use_subsidence_tendencies )  THEN
976                CALL subsidence( i, j, tend, pt, pt_init, 2 )
977             ENDIF
978
979#if defined( __rrtmg )
980!
981!--          If required, add tendency due to radiative heating/cooling
982             IF ( radiation  .AND.                                             &
983                  simulated_time > skip_time_do_radiation )  THEN
984                CALL radiation_tendency ( i, j, tend )
985             ENDIF
986#endif
987
988             CALL module_interface_actions( i, j, 'pt-tendency' )
989!
990!--          Prognostic equation for potential temperature
991             DO  k = nzb+1, nzt
992                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
993                                                  ( tsc(2) * tend(k,j,i) +     &
994                                                    tsc(3) * tpt_m(k,j,i) )    &
995                                                  - tsc(5)                     &
996                                                  * ( pt(k,j,i) - pt_init(k) ) &
997                                                  * ( rdf_sc(k) + ptdf_x(i)    &
998                                                                + ptdf_y(j) )  &
999                                          )                                    &
1000                                       * MERGE( 1.0_wp, 0.0_wp,                &
1001                                                BTEST( wall_flags_0(k,j,i), 0 )&
1002                                              )
1003             ENDDO
1004
1005!
1006!--          Calculate tendencies for the next Runge-Kutta step
1007             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1008                IF ( intermediate_timestep_count == 1 )  THEN
1009                   DO  k = nzb+1, nzt
1010                      tpt_m(k,j,i) = tend(k,j,i)
1011                   ENDDO
1012                ELSEIF ( intermediate_timestep_count < &
1013                         intermediate_timestep_count_max )  THEN
1014                   DO  k = nzb+1, nzt
1015                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1016                                        5.3125_wp * tpt_m(k,j,i)
1017                   ENDDO
1018                ENDIF
1019             ENDIF
1020
1021          ENDIF
1022
1023!
1024!--       If required, compute prognostic equation for total water content
1025          IF ( humidity )  THEN
1026
1027!
1028!--          Tendency-terms for total water content / scalar
1029             tend(:,j,i) = 0.0_wp
1030             IF ( timestep_scheme(1:5) == 'runge' ) &
1031             THEN
1032                IF ( ws_scheme_sca )  THEN
1033                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
1034                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1035                ELSE
1036                   CALL advec_s_pw( i, j, q )
1037                ENDIF
1038             ELSE
1039                CALL advec_s_up( i, j, q )
1040             ENDIF
1041             CALL diffusion_s( i, j, q,                                        &
1042                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
1043                               surf_def_h(2)%qsws,                             &
1044                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
1045                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
1046                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
1047                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
1048                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
1049                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
1050                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1051
1052!
1053!--          Sink or source of humidity due to canopy elements
1054             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
1055
1056!
1057!--          Large scale advection
1058             IF ( large_scale_forcing )  THEN
1059                CALL ls_advec( i, j, simulated_time, 'q' )
1060             ENDIF
1061
1062!
1063!--          Nudging
1064             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
1065
1066!
1067!--          If required compute influence of large-scale subsidence/ascent
1068             IF ( large_scale_subsidence  .AND.                                &
1069                  .NOT. use_subsidence_tendencies )  THEN
1070                CALL subsidence( i, j, tend, q, q_init, 3 )
1071             ENDIF
1072
1073             CALL module_interface_actions( i, j, 'q-tendency' )
1074
1075!
1076!--          Prognostic equation for total water content / scalar
1077             DO  k = nzb+1, nzt
1078                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
1079                                                ( tsc(2) * tend(k,j,i) +       &
1080                                                  tsc(3) * tq_m(k,j,i) )       &
1081                                                - tsc(5) * rdf_sc(k) *         &
1082                                                      ( q(k,j,i) - q_init(k) ) &
1083                                        )                                      &
1084                                       * MERGE( 1.0_wp, 0.0_wp,                &
1085                                                BTEST( wall_flags_0(k,j,i), 0 )&
1086                                              )               
1087                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1088             ENDDO
1089
1090!
1091!--          Calculate tendencies for the next Runge-Kutta step
1092             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1093                IF ( intermediate_timestep_count == 1 )  THEN
1094                   DO  k = nzb+1, nzt
1095                      tq_m(k,j,i) = tend(k,j,i)
1096                   ENDDO
1097                ELSEIF ( intermediate_timestep_count < &
1098                         intermediate_timestep_count_max )  THEN
1099                   DO  k = nzb+1, nzt
1100                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1101                                       5.3125_wp * tq_m(k,j,i)
1102                   ENDDO
1103                ENDIF
1104             ENDIF
1105
1106!
1107!--          If required, calculate prognostic equations for cloud water content
1108!--          and cloud drop concentration
1109             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1110!
1111!--             Calculate prognostic equation for cloud water content
1112                tend(:,j,i) = 0.0_wp
1113                IF ( timestep_scheme(1:5) == 'runge' ) &
1114                THEN
1115                   IF ( ws_scheme_sca )  THEN
1116                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
1117                                       diss_s_qc, flux_l_qc, diss_l_qc, &
1118                                       i_omp_start, tn )
1119                   ELSE
1120                      CALL advec_s_pw( i, j, qc )
1121                   ENDIF
1122                ELSE
1123                   CALL advec_s_up( i, j, qc )
1124                ENDIF
1125                CALL diffusion_s( i, j, qc,                                   &
1126                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
1127                                  surf_def_h(2)%qcsws,                        &
1128                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      & 
1129                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
1130                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
1131                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
1132                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
1133                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
1134                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
1135
1136!
1137!--             Prognostic equation for cloud water content
1138                DO  k = nzb+1, nzt
1139                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
1140                                                      ( tsc(2) * tend(k,j,i) + &
1141                                                        tsc(3) * tqc_m(k,j,i) )&
1142                                                      - tsc(5) * rdf_sc(k)     &
1143                                                               * qc(k,j,i)     &
1144                                             )                                 &
1145                                       * MERGE( 1.0_wp, 0.0_wp,                &
1146                                                BTEST( wall_flags_0(k,j,i), 0 )&
1147                                              ) 
1148                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
1149                ENDDO
1150!
1151!--             Calculate tendencies for the next Runge-Kutta step
1152                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1153                   IF ( intermediate_timestep_count == 1 )  THEN
1154                      DO  k = nzb+1, nzt
1155                         tqc_m(k,j,i) = tend(k,j,i)
1156                      ENDDO
1157                   ELSEIF ( intermediate_timestep_count < &
1158                            intermediate_timestep_count_max )  THEN
1159                      DO  k = nzb+1, nzt
1160                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1161                                           5.3125_wp * tqc_m(k,j,i)
1162                      ENDDO
1163                   ENDIF
1164                ENDIF
1165
1166!
1167!--             Calculate prognostic equation for cloud drop concentration.
1168                tend(:,j,i) = 0.0_wp
1169                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1170                   IF ( ws_scheme_sca )  THEN
1171                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
1172                                    diss_s_nc, flux_l_nc, diss_l_nc, &
1173                                    i_omp_start, tn )
1174                   ELSE
1175                      CALL advec_s_pw( i, j, nc )
1176                   ENDIF
1177                ELSE
1178                   CALL advec_s_up( i, j, nc )
1179                ENDIF
1180                CALL diffusion_s( i, j, nc,                                    &
1181                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
1182                                  surf_def_h(2)%ncsws,                         &
1183                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
1184                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
1185                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
1186                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
1187                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
1188                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
1189                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
1190
1191!
1192!--             Prognostic equation for cloud drop concentration
1193                DO  k = nzb+1, nzt
1194                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
1195                                                      ( tsc(2) * tend(k,j,i) + &
1196                                                        tsc(3) * tnc_m(k,j,i) )&
1197                                                      - tsc(5) * rdf_sc(k)     &
1198                                                               * nc(k,j,i)     &
1199                                             )                                 &
1200                                       * MERGE( 1.0_wp, 0.0_wp,                &
1201                                                BTEST( wall_flags_0(k,j,i), 0 )&
1202                                              )
1203                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
1204                ENDDO
1205!
1206!--             Calculate tendencies for the next Runge-Kutta step
1207                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1208                   IF ( intermediate_timestep_count == 1 )  THEN
1209                      DO  k = nzb+1, nzt
1210                         tnc_m(k,j,i) = tend(k,j,i)
1211                      ENDDO
1212                   ELSEIF ( intermediate_timestep_count < &
1213                            intermediate_timestep_count_max )  THEN
1214                      DO  k = nzb+1, nzt
1215                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1216                                           5.3125_wp * tnc_m(k,j,i)
1217                      ENDDO
1218                   ENDIF
1219                ENDIF
1220
1221             ENDIF
1222!
1223!--          If required, calculate prognostic equations for rain water content
1224!--          and rain drop concentration
1225             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1226!
1227!--             Calculate prognostic equation for rain water content
1228                tend(:,j,i) = 0.0_wp
1229                IF ( timestep_scheme(1:5) == 'runge' ) &
1230                THEN
1231                   IF ( ws_scheme_sca )  THEN
1232                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
1233                                       diss_s_qr, flux_l_qr, diss_l_qr, &
1234                                       i_omp_start, tn )
1235                   ELSE
1236                      CALL advec_s_pw( i, j, qr )
1237                   ENDIF
1238                ELSE
1239                   CALL advec_s_up( i, j, qr )
1240                ENDIF
1241                CALL diffusion_s( i, j, qr,                                   &
1242                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
1243                                  surf_def_h(2)%qrsws,                        &
1244                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      & 
1245                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
1246                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
1247                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
1248                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
1249                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
1250                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1251
1252!
1253!--             Prognostic equation for rain water content
1254                DO  k = nzb+1, nzt
1255                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
1256                                                      ( tsc(2) * tend(k,j,i) + &
1257                                                        tsc(3) * tqr_m(k,j,i) )&
1258                                                      - tsc(5) * rdf_sc(k)     &
1259                                                               * qr(k,j,i)     &
1260                                             )                                 &
1261                                       * MERGE( 1.0_wp, 0.0_wp,                &
1262                                                BTEST( wall_flags_0(k,j,i), 0 )&
1263                                              ) 
1264                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
1265                ENDDO
1266!
1267!--             Calculate tendencies for the next Runge-Kutta step
1268                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1269                   IF ( intermediate_timestep_count == 1 )  THEN
1270                      DO  k = nzb+1, nzt
1271                         tqr_m(k,j,i) = tend(k,j,i)
1272                      ENDDO
1273                   ELSEIF ( intermediate_timestep_count < &
1274                            intermediate_timestep_count_max )  THEN
1275                      DO  k = nzb+1, nzt
1276                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1277                                           5.3125_wp * tqr_m(k,j,i)
1278                      ENDDO
1279                   ENDIF
1280                ENDIF
1281
1282!
1283!--             Calculate prognostic equation for rain drop concentration.
1284                tend(:,j,i) = 0.0_wp
1285                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1286                   IF ( ws_scheme_sca )  THEN
1287                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
1288                                    diss_s_nr, flux_l_nr, diss_l_nr, &
1289                                    i_omp_start, tn )
1290                   ELSE
1291                      CALL advec_s_pw( i, j, nr )
1292                   ENDIF
1293                ELSE
1294                   CALL advec_s_up( i, j, nr )
1295                ENDIF
1296                CALL diffusion_s( i, j, nr,                                    &
1297                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
1298                                  surf_def_h(2)%nrsws,                         &
1299                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
1300                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
1301                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
1302                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
1303                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
1304                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
1305                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1306
1307!
1308!--             Prognostic equation for rain drop concentration
1309                DO  k = nzb+1, nzt
1310                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
1311                                                      ( tsc(2) * tend(k,j,i) + &
1312                                                        tsc(3) * tnr_m(k,j,i) )&
1313                                                      - tsc(5) * rdf_sc(k)     &
1314                                                               * nr(k,j,i)     &
1315                                             )                                 &
1316                                       * MERGE( 1.0_wp, 0.0_wp,                &
1317                                                BTEST( wall_flags_0(k,j,i), 0 )&
1318                                              )
1319                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1320                ENDDO
1321!
1322!--             Calculate tendencies for the next Runge-Kutta step
1323                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1324                   IF ( intermediate_timestep_count == 1 )  THEN
1325                      DO  k = nzb+1, nzt
1326                         tnr_m(k,j,i) = tend(k,j,i)
1327                      ENDDO
1328                   ELSEIF ( intermediate_timestep_count < &
1329                            intermediate_timestep_count_max )  THEN
1330                      DO  k = nzb+1, nzt
1331                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1332                                           5.3125_wp * tnr_m(k,j,i)
1333                      ENDDO
1334                   ENDIF
1335                ENDIF
1336
1337             ENDIF
1338
1339          ENDIF
1340
1341!
1342!--       If required, compute prognostic equation for scalar
1343          IF ( passive_scalar )  THEN
1344!
1345!--          Tendency-terms for total water content / scalar
1346             tend(:,j,i) = 0.0_wp
1347             IF ( timestep_scheme(1:5) == 'runge' ) &
1348             THEN
1349                IF ( ws_scheme_sca )  THEN
1350                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1351                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1352                ELSE
1353                   CALL advec_s_pw( i, j, s )
1354                ENDIF
1355             ELSE
1356                CALL advec_s_up( i, j, s )
1357             ENDIF
1358             CALL diffusion_s( i, j, s,                                        &
1359                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1360                               surf_def_h(2)%ssws,                             &
1361                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1362                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1363                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1364                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1365                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1366                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1367                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1368
1369!
1370!--          Sink or source of scalar concentration due to canopy elements
1371             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1372
1373!
1374!--          Large scale advection, still need to be extended for scalars
1375!              IF ( large_scale_forcing )  THEN
1376!                 CALL ls_advec( i, j, simulated_time, 's' )
1377!              ENDIF
1378
1379!
1380!--          Nudging, still need to be extended for scalars
1381!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1382
1383!
1384!--          If required compute influence of large-scale subsidence/ascent.
1385!--          Note, the last argument is of no meaning in this case, as it is
1386!--          only used in conjunction with large_scale_forcing, which is to
1387!--          date not implemented for scalars.
1388             IF ( large_scale_subsidence  .AND.                                &
1389                  .NOT. use_subsidence_tendencies  .AND.                       &
1390                  .NOT. large_scale_forcing )  THEN
1391                CALL subsidence( i, j, tend, s, s_init, 3 )
1392             ENDIF
1393
1394             CALL module_interface_actions( i, j, 's-tendency' )
1395
1396!
1397!--          Prognostic equation for scalar
1398             DO  k = nzb+1, nzt
1399                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1400                                            ( tsc(2) * tend(k,j,i) +           &
1401                                              tsc(3) * ts_m(k,j,i) )           &
1402                                            - tsc(5) * rdf_sc(k)               &
1403                                                     * ( s(k,j,i) - s_init(k) )&
1404                                        )                                      &
1405                                       * MERGE( 1.0_wp, 0.0_wp,                &
1406                                                BTEST( wall_flags_0(k,j,i), 0 )&
1407                                              )
1408                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1409             ENDDO
1410
1411!
1412!--          Calculate tendencies for the next Runge-Kutta step
1413             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1414                IF ( intermediate_timestep_count == 1 )  THEN
1415                   DO  k = nzb+1, nzt
1416                      ts_m(k,j,i) = tend(k,j,i)
1417                   ENDDO
1418                ELSEIF ( intermediate_timestep_count < &
1419                         intermediate_timestep_count_max )  THEN
1420                   DO  k = nzb+1, nzt
1421                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1422                                       5.3125_wp * ts_m(k,j,i)
1423                   ENDDO
1424                ENDIF
1425             ENDIF
1426
1427          ENDIF
1428!
1429!--       Calculate prognostic equations for turbulence closure
1430          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1431!
1432!--       Calculate prognostic equations for all other modules
1433          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1434
1435!
1436!--       Calculate prognostic equation for chemical quantites
1437          IF ( air_chemistry )  THEN
1438             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
1439!
1440!--          Loop over chemical species
1441             DO  lsp = 1, nvar                         
1442                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
1443                                     chem_species(lsp)%conc,                   &
1444                                     chem_species(lsp)%tconc_m,                &
1445                                     chem_species(lsp)%conc_pr_init,           &
1446                                     i, j, i_omp_start, tn, lsp,               &
1447                                     chem_species(lsp)%flux_s_cs,              &
1448                                     chem_species(lsp)%diss_s_cs,              &
1449                                     chem_species(lsp)%flux_l_cs,              &
1450                                     chem_species(lsp)%diss_l_cs )       
1451             ENDDO
1452
1453          ENDIF   ! Chemical equations
1454
1455       ENDDO  ! loop over j
1456    ENDDO  ! loop over i
1457    !$OMP END PARALLEL
1458
1459
1460    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1461
1462
1463 END SUBROUTINE prognostic_equations_cache
1464
1465
1466!------------------------------------------------------------------------------!
1467! Description:
1468! ------------
1469!> Version for vector machines
1470!------------------------------------------------------------------------------!
1471
1472 SUBROUTINE prognostic_equations_vector
1473
1474
1475    IMPLICIT NONE
1476
1477    INTEGER(iwp) ::  i     !<
1478    INTEGER(iwp) ::  ib    !< index for aerosol size bins (salsa)
1479    INTEGER(iwp) ::  ic    !< index for chemical compounds (salsa)
1480    INTEGER(iwp) ::  icc   !< additional index for chemical compounds (salsa)
1481    INTEGER(iwp) ::  ig    !< index for gaseous compounds (salsa)
1482    INTEGER(iwp) ::  j     !<
1483    INTEGER(iwp) ::  k     !<
1484    INTEGER(iwp) ::  lsp   !< running index for chemical species
1485
1486    REAL(wp)     ::  sbt  !<
1487
1488!
1489!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1490!-- step. The exchange of ghost points is required after this update of the
1491!-- concentrations of aerosol number and mass
1492    IF ( salsa )  THEN
1493       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1494          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
1495          THEN
1496             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1497             !$OMP PARALLEL PRIVATE (i,j,ib,ic,icc,ig)
1498             !$OMP DO
1499!
1500!--          Call salsa processes
1501             DO  i = nxl, nxr
1502                DO  j = nys, nyn
1503                   CALL salsa_diagnostics( i, j )
1504                   CALL salsa_driver( i, j, 3 )
1505                   CALL salsa_diagnostics( i, j )
1506                ENDDO
1507             ENDDO
1508
1509             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1510             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1511!
1512!--          Exchange ghost points and decycle if needed.
1513             DO  ib = 1, nbins_aerosol
1514                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
1515                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
1516                                           aerosol_number(ib)%init )
1517                DO  ic = 1, ncomponents_mass
1518                   icc = ( ic - 1 ) * nbins_aerosol + ib
1519                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
1520                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
1521                                              aerosol_mass(icc)%init )
1522                ENDDO
1523             ENDDO
1524             IF ( .NOT. salsa_gases_from_chem )  THEN
1525                DO  ig = 1, ngases_salsa
1526                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
1527                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
1528                                              salsa_gas(ig)%init )
1529                ENDDO
1530             ENDIF
1531             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1532             !$OMP END PARALLEL
1533             last_salsa_time = time_since_reference_point
1534          ENDIF
1535       ENDIF
1536    ENDIF
1537
1538!
1539!-- If required, calculate cloud microphysical impacts
1540    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
1541         ( intermediate_timestep_count == 1  .OR.                              &
1542           call_microphysics_at_all_substeps )                                 &
1543       )  THEN
1544       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1545       CALL bcm_actions
1546       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1547    ENDIF
1548
1549!
1550!-- u-velocity component
1551    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1552
1553    !$ACC KERNELS PRESENT(tend)
1554    tend = 0.0_wp
1555    !$ACC END KERNELS
1556    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1557       IF ( ws_scheme_mom )  THEN
1558          CALL advec_u_ws
1559       ELSE
1560          CALL advec_u_pw
1561       ENDIF
1562    ELSE
1563       CALL advec_u_up
1564    ENDIF
1565    CALL diffusion_u
1566    CALL coriolis( 1 )
1567    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1568       CALL buoyancy( pt, 1 )
1569    ENDIF
1570
1571!
1572!-- Drag by plant canopy
1573    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1574
1575!
1576!-- External pressure gradient
1577    IF ( dp_external )  THEN
1578       DO  i = nxlu, nxr
1579          DO  j = nys, nyn
1580             DO  k = dp_level_ind_b+1, nzt
1581                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1582             ENDDO
1583          ENDDO
1584       ENDDO
1585    ENDIF
1586
1587!
1588!-- Nudging
1589    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1590
1591!
1592!-- Effect of Stokes drift (in ocean mode only)
1593    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1594
1595!
1596!-- Forces by wind turbines
1597    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1598
1599    CALL module_interface_actions( 'u-tendency' )
1600
1601!
1602!-- Prognostic equation for u-velocity component
1603    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1604    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1605    !$ACC PRESENT(tsc(2:5)) &
1606    !$ACC PRESENT(u_p)
1607    DO  i = nxlu, nxr
1608       DO  j = nys, nyn
1609          DO  k = nzb+1, nzt
1610             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1611                                                 tsc(3) * tu_m(k,j,i) )          &
1612                                               - tsc(5) * rdf(k) *               &
1613                                                        ( u(k,j,i) - u_init(k) ) &
1614                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1615                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1616                                              )
1617          ENDDO
1618       ENDDO
1619    ENDDO
1620
1621!
1622!-- Add turbulence generated by wave breaking (in ocean mode only)
1623    IF ( wave_breaking  .AND.                                                  &
1624         intermediate_timestep_count == intermediate_timestep_count_max )      &
1625    THEN
1626       CALL wave_breaking_term( 1 )
1627    ENDIF
1628
1629!
1630!-- Calculate tendencies for the next Runge-Kutta step
1631    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1632       IF ( intermediate_timestep_count == 1 )  THEN
1633          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1634          !$ACC PRESENT(tend, tu_m)
1635          DO  i = nxlu, nxr
1636             DO  j = nys, nyn
1637                DO  k = nzb+1, nzt
1638                   tu_m(k,j,i) = tend(k,j,i)
1639                ENDDO
1640             ENDDO
1641          ENDDO
1642       ELSEIF ( intermediate_timestep_count < &
1643                intermediate_timestep_count_max )  THEN
1644          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1645          !$ACC PRESENT(tend, tu_m)
1646          DO  i = nxlu, nxr
1647             DO  j = nys, nyn
1648                DO  k = nzb+1, nzt
1649                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1650                                   + 5.3125_wp * tu_m(k,j,i)
1651                ENDDO
1652             ENDDO
1653          ENDDO
1654       ENDIF
1655    ENDIF
1656
1657    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1658
1659!
1660!-- v-velocity component
1661    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1662
1663    !$ACC KERNELS PRESENT(tend)
1664    tend = 0.0_wp
1665    !$ACC END KERNELS
1666    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1667       IF ( ws_scheme_mom )  THEN
1668          CALL advec_v_ws
1669       ELSE
1670          CALL advec_v_pw
1671       END IF
1672    ELSE
1673       CALL advec_v_up
1674    ENDIF
1675    CALL diffusion_v
1676    CALL coriolis( 2 )
1677
1678!
1679!-- Drag by plant canopy
1680    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1681
1682!
1683!-- External pressure gradient
1684    IF ( dp_external )  THEN
1685       DO  i = nxl, nxr
1686          DO  j = nysv, nyn
1687             DO  k = dp_level_ind_b+1, nzt
1688                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1689             ENDDO
1690          ENDDO
1691       ENDDO
1692    ENDIF
1693
1694!
1695!-- Nudging
1696    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1697
1698!
1699!-- Effect of Stokes drift (in ocean mode only)
1700    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1701
1702!
1703!-- Forces by wind turbines
1704    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1705
1706    CALL module_interface_actions( 'v-tendency' )
1707
1708!
1709!-- Prognostic equation for v-velocity component
1710    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1711    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1712    !$ACC PRESENT(tsc(2:5)) &
1713    !$ACC PRESENT(v_p)
1714    DO  i = nxl, nxr
1715       DO  j = nysv, nyn
1716          DO  k = nzb+1, nzt
1717             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1718                                                 tsc(3) * tv_m(k,j,i) )        &
1719                                               - tsc(5) * rdf(k) *             &
1720                                                      ( v(k,j,i) - v_init(k) ) &
1721                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1722                                                BTEST( wall_flags_0(k,j,i), 2 )&
1723                                              )
1724          ENDDO
1725       ENDDO
1726    ENDDO
1727
1728!
1729!-- Add turbulence generated by wave breaking (in ocean mode only)
1730    IF ( wave_breaking  .AND.                                                  &
1731         intermediate_timestep_count == intermediate_timestep_count_max )      &
1732    THEN
1733       CALL wave_breaking_term( 2 )
1734    ENDIF
1735
1736!
1737!-- Calculate tendencies for the next Runge-Kutta step
1738    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1739       IF ( intermediate_timestep_count == 1 )  THEN
1740          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1741          !$ACC PRESENT(tend, tv_m)
1742          DO  i = nxl, nxr
1743             DO  j = nysv, nyn
1744                DO  k = nzb+1, nzt
1745                   tv_m(k,j,i) = tend(k,j,i)
1746                ENDDO
1747             ENDDO
1748          ENDDO
1749       ELSEIF ( intermediate_timestep_count < &
1750                intermediate_timestep_count_max )  THEN
1751          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1752          !$ACC PRESENT(tend, tv_m)
1753          DO  i = nxl, nxr
1754             DO  j = nysv, nyn
1755                DO  k = nzb+1, nzt
1756                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1757                                  + 5.3125_wp * tv_m(k,j,i)
1758                ENDDO
1759             ENDDO
1760          ENDDO
1761       ENDIF
1762    ENDIF
1763
1764    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1765
1766!
1767!-- w-velocity component
1768    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1769
1770    !$ACC KERNELS PRESENT(tend)
1771    tend = 0.0_wp
1772    !$ACC END KERNELS
1773    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1774       IF ( ws_scheme_mom )  THEN
1775          CALL advec_w_ws
1776       ELSE
1777          CALL advec_w_pw
1778       ENDIF
1779    ELSE
1780       CALL advec_w_up
1781    ENDIF
1782    CALL diffusion_w
1783    CALL coriolis( 3 )
1784
1785    IF ( .NOT. neutral )  THEN
1786       IF ( ocean_mode )  THEN
1787          CALL buoyancy( rho_ocean, 3 )
1788       ELSE
1789          IF ( .NOT. humidity )  THEN
1790             CALL buoyancy( pt, 3 )
1791          ELSE
1792             CALL buoyancy( vpt, 3 )
1793          ENDIF
1794       ENDIF
1795    ENDIF
1796
1797!
1798!-- Drag by plant canopy
1799    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1800
1801!
1802!-- Effect of Stokes drift (in ocean mode only)
1803    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1804
1805!
1806!-- Forces by wind turbines
1807    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1808
1809    CALL module_interface_actions( 'w-tendency' )
1810
1811!
1812!-- Prognostic equation for w-velocity component
1813    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1814    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1815    !$ACC PRESENT(tsc(2:5)) &
1816    !$ACC PRESENT(w_p)
1817    DO  i = nxl, nxr
1818       DO  j = nys, nyn
1819          DO  k = nzb+1, nzt-1
1820             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1821                                                 tsc(3) * tw_m(k,j,i) )        &
1822                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1823                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1824                                                BTEST( wall_flags_0(k,j,i), 3 )&
1825                                              )
1826          ENDDO
1827       ENDDO
1828    ENDDO
1829
1830!
1831!-- Calculate tendencies for the next Runge-Kutta step
1832    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1833       IF ( intermediate_timestep_count == 1 )  THEN
1834          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1835          !$ACC PRESENT(tend, tw_m)
1836          DO  i = nxl, nxr
1837             DO  j = nys, nyn
1838                DO  k = nzb+1, nzt-1
1839                   tw_m(k,j,i) = tend(k,j,i)
1840                ENDDO
1841             ENDDO
1842          ENDDO
1843       ELSEIF ( intermediate_timestep_count < &
1844                intermediate_timestep_count_max )  THEN
1845          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1846          !$ACC PRESENT(tend, tw_m)
1847          DO  i = nxl, nxr
1848             DO  j = nys, nyn
1849                DO  k = nzb+1, nzt-1
1850                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1851                                  + 5.3125_wp * tw_m(k,j,i)
1852                ENDDO
1853             ENDDO
1854          ENDDO
1855       ENDIF
1856    ENDIF
1857
1858    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1859
1860
1861!
1862!-- If required, compute prognostic equation for potential temperature
1863    IF ( .NOT. neutral )  THEN
1864
1865       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1866
1867!
1868!--    pt-tendency terms with communication
1869       sbt = tsc(2)
1870       IF ( scalar_advec == 'bc-scheme' )  THEN
1871
1872          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1873!
1874!--          Bott-Chlond scheme always uses Euler time step. Thus:
1875             sbt = 1.0_wp
1876          ENDIF
1877          tend = 0.0_wp
1878          CALL advec_s_bc( pt, 'pt' )
1879
1880       ENDIF
1881
1882!
1883!--    pt-tendency terms with no communication
1884       IF ( scalar_advec /= 'bc-scheme' )  THEN
1885          !$ACC KERNELS PRESENT(tend)
1886          tend = 0.0_wp
1887          !$ACC END KERNELS
1888          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1889             IF ( ws_scheme_sca )  THEN
1890                CALL advec_s_ws( pt, 'pt' )
1891             ELSE
1892                CALL advec_s_pw( pt )
1893             ENDIF
1894          ELSE
1895             CALL advec_s_up( pt )
1896          ENDIF
1897       ENDIF
1898
1899       CALL diffusion_s( pt,                                                   &
1900                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1901                         surf_def_h(2)%shf,                                    &
1902                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1903                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1904                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1905                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1906                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1907                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1908                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1909
1910!
1911!--    Consideration of heat sources within the plant canopy
1912       IF ( plant_canopy  .AND.                                          &
1913            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1914          CALL pcm_tendency( 4 )
1915       ENDIF
1916
1917!
1918!--    Large scale advection
1919       IF ( large_scale_forcing )  THEN
1920          CALL ls_advec( simulated_time, 'pt' )
1921       ENDIF
1922
1923!
1924!--    Nudging
1925       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1926
1927!
1928!--    If required compute influence of large-scale subsidence/ascent
1929       IF ( large_scale_subsidence  .AND.                                      &
1930            .NOT. use_subsidence_tendencies )  THEN
1931          CALL subsidence( tend, pt, pt_init, 2 )
1932       ENDIF
1933
1934#if defined( __rrtmg )
1935!
1936!--    If required, add tendency due to radiative heating/cooling
1937       IF ( radiation  .AND.                                                   &
1938            simulated_time > skip_time_do_radiation )  THEN
1939            CALL radiation_tendency ( tend )
1940       ENDIF
1941#endif
1942
1943       CALL module_interface_actions( 'pt-tendency' )
1944
1945!
1946!--    Prognostic equation for potential temperature
1947       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1948       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1949       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1950       !$ACC PRESENT(tsc(3:5)) &
1951       !$ACC PRESENT(pt_p)
1952       DO  i = nxl, nxr
1953          DO  j = nys, nyn
1954             DO  k = nzb+1, nzt
1955                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1956                                                   tsc(3) * tpt_m(k,j,i) )     &
1957                                                 - tsc(5) *                    &
1958                                                   ( pt(k,j,i) - pt_init(k) ) *&
1959                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1960                                          )                                    &
1961                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1962                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1963                                          )
1964             ENDDO
1965          ENDDO
1966       ENDDO
1967!
1968!--    Calculate tendencies for the next Runge-Kutta step
1969       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1970          IF ( intermediate_timestep_count == 1 )  THEN
1971             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1972             !$ACC PRESENT(tend, tpt_m)
1973             DO  i = nxl, nxr
1974                DO  j = nys, nyn
1975                   DO  k = nzb+1, nzt
1976                      tpt_m(k,j,i) = tend(k,j,i)
1977                   ENDDO
1978                ENDDO
1979             ENDDO
1980          ELSEIF ( intermediate_timestep_count < &
1981                   intermediate_timestep_count_max )  THEN
1982             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1983             !$ACC PRESENT(tend, tpt_m)
1984             DO  i = nxl, nxr
1985                DO  j = nys, nyn
1986                   DO  k = nzb+1, nzt
1987                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1988                                        5.3125_wp * tpt_m(k,j,i)
1989                   ENDDO
1990                ENDDO
1991             ENDDO
1992          ENDIF
1993       ENDIF
1994
1995       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1996
1997    ENDIF
1998
1999!
2000!-- If required, compute prognostic equation for total water content
2001    IF ( humidity )  THEN
2002
2003       CALL cpu_log( log_point(29), 'q-equation', 'start' )
2004
2005!
2006!--    Scalar/q-tendency terms with communication
2007       sbt = tsc(2)
2008       IF ( scalar_advec == 'bc-scheme' )  THEN
2009
2010          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2011!
2012!--          Bott-Chlond scheme always uses Euler time step. Thus:
2013             sbt = 1.0_wp
2014          ENDIF
2015          tend = 0.0_wp
2016          CALL advec_s_bc( q, 'q' )
2017
2018       ENDIF
2019
2020!
2021!--    Scalar/q-tendency terms with no communication
2022       IF ( scalar_advec /= 'bc-scheme' )  THEN
2023          tend = 0.0_wp
2024          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2025             IF ( ws_scheme_sca )  THEN
2026                CALL advec_s_ws( q, 'q' )
2027             ELSE
2028                CALL advec_s_pw( q )
2029             ENDIF
2030          ELSE
2031             CALL advec_s_up( q )
2032          ENDIF
2033       ENDIF
2034
2035       CALL diffusion_s( q,                                                    &
2036                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
2037                         surf_def_h(2)%qsws,                                   &
2038                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
2039                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
2040                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
2041                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
2042                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
2043                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
2044                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
2045
2046!
2047!--    Sink or source of humidity due to canopy elements
2048       IF ( plant_canopy ) CALL pcm_tendency( 5 )
2049
2050!
2051!--    Large scale advection
2052       IF ( large_scale_forcing )  THEN
2053          CALL ls_advec( simulated_time, 'q' )
2054       ENDIF
2055
2056!
2057!--    Nudging
2058       IF ( nudging )  CALL nudge( simulated_time, 'q' )
2059
2060!
2061!--    If required compute influence of large-scale subsidence/ascent
2062       IF ( large_scale_subsidence  .AND.                                      &
2063            .NOT. use_subsidence_tendencies )  THEN
2064         CALL subsidence( tend, q, q_init, 3 )
2065       ENDIF
2066
2067       CALL module_interface_actions( 'q-tendency' )
2068
2069!
2070!--    Prognostic equation for total water content
2071       DO  i = nxl, nxr
2072          DO  j = nys, nyn
2073             DO  k = nzb+1, nzt
2074                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2075                                                 tsc(3) * tq_m(k,j,i) )        &
2076                                               - tsc(5) * rdf_sc(k) *          &
2077                                                      ( q(k,j,i) - q_init(k) ) &
2078                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
2079                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
2080                                                 )
2081                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
2082             ENDDO
2083          ENDDO
2084       ENDDO
2085
2086!
2087!--    Calculate tendencies for the next Runge-Kutta step
2088       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2089          IF ( intermediate_timestep_count == 1 )  THEN
2090             DO  i = nxl, nxr
2091                DO  j = nys, nyn
2092                   DO  k = nzb+1, nzt
2093                      tq_m(k,j,i) = tend(k,j,i)
2094                   ENDDO
2095                ENDDO
2096             ENDDO
2097          ELSEIF ( intermediate_timestep_count < &
2098                   intermediate_timestep_count_max )  THEN
2099             DO  i = nxl, nxr
2100                DO  j = nys, nyn
2101                   DO  k = nzb+1, nzt
2102                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2103                                     + 5.3125_wp * tq_m(k,j,i)
2104                   ENDDO
2105                ENDDO
2106             ENDDO
2107          ENDIF
2108       ENDIF
2109
2110       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
2111
2112!
2113!--    If required, calculate prognostic equations for cloud water content
2114!--    and cloud drop concentration
2115       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
2116
2117          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
2118
2119!
2120!--       Calculate prognostic equation for cloud water content
2121          sbt = tsc(2)
2122          IF ( scalar_advec == 'bc-scheme' )  THEN
2123
2124             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2125!
2126!--             Bott-Chlond scheme always uses Euler time step. Thus:
2127                sbt = 1.0_wp
2128             ENDIF
2129             tend = 0.0_wp
2130             CALL advec_s_bc( qc, 'qc' )
2131
2132          ENDIF
2133
2134!
2135!--       qc-tendency terms with no communication
2136          IF ( scalar_advec /= 'bc-scheme' )  THEN
2137             tend = 0.0_wp
2138             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2139                IF ( ws_scheme_sca )  THEN
2140                   CALL advec_s_ws( qc, 'qc' )
2141                ELSE
2142                   CALL advec_s_pw( qc )
2143                ENDIF
2144             ELSE
2145                CALL advec_s_up( qc )
2146             ENDIF
2147          ENDIF
2148
2149          CALL diffusion_s( qc,                                                &
2150                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
2151                            surf_def_h(2)%qcsws,                               &
2152                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
2153                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
2154                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
2155                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
2156                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
2157                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
2158                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
2159
2160!
2161!--       Prognostic equation for cloud water content
2162          DO  i = nxl, nxr
2163             DO  j = nys, nyn
2164                DO  k = nzb+1, nzt
2165                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2166                                                      tsc(3) * tqc_m(k,j,i) )  &
2167                                                    - tsc(5) * rdf_sc(k) *     &
2168                                                               qc(k,j,i)       &
2169                                             )                                 &
2170                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2171                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2172                                          )
2173                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
2174                ENDDO
2175             ENDDO
2176          ENDDO
2177
2178!
2179!--       Calculate tendencies for the next Runge-Kutta step
2180          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2181             IF ( intermediate_timestep_count == 1 )  THEN
2182                DO  i = nxl, nxr
2183                   DO  j = nys, nyn
2184                      DO  k = nzb+1, nzt
2185                         tqc_m(k,j,i) = tend(k,j,i)
2186                      ENDDO
2187                   ENDDO
2188                ENDDO
2189             ELSEIF ( intermediate_timestep_count < &
2190                      intermediate_timestep_count_max )  THEN
2191                DO  i = nxl, nxr
2192                   DO  j = nys, nyn
2193                      DO  k = nzb+1, nzt
2194                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2195                                         + 5.3125_wp * tqc_m(k,j,i)
2196                      ENDDO
2197                   ENDDO
2198                ENDDO
2199             ENDIF
2200          ENDIF
2201
2202          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
2203
2204          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
2205!
2206!--       Calculate prognostic equation for cloud drop concentration
2207          sbt = tsc(2)
2208          IF ( scalar_advec == 'bc-scheme' )  THEN
2209
2210             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2211!
2212!--             Bott-Chlond scheme always uses Euler time step. Thus:
2213                sbt = 1.0_wp
2214             ENDIF
2215             tend = 0.0_wp
2216             CALL advec_s_bc( nc, 'nc' )
2217
2218          ENDIF
2219
2220!
2221!--       nc-tendency terms with no communication
2222          IF ( scalar_advec /= 'bc-scheme' )  THEN
2223             tend = 0.0_wp
2224             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2225                IF ( ws_scheme_sca )  THEN
2226                   CALL advec_s_ws( nc, 'nc' )
2227                ELSE
2228                   CALL advec_s_pw( nc )
2229                ENDIF
2230             ELSE
2231                CALL advec_s_up( nc )
2232             ENDIF
2233          ENDIF
2234
2235          CALL diffusion_s( nc,                                                &
2236                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
2237                            surf_def_h(2)%ncsws,                               &
2238                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
2239                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
2240                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
2241                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
2242                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
2243                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
2244                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
2245
2246!
2247!--       Prognostic equation for cloud drop concentration
2248          DO  i = nxl, nxr
2249             DO  j = nys, nyn
2250                DO  k = nzb+1, nzt
2251                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2252                                                      tsc(3) * tnc_m(k,j,i) )  &
2253                                                    - tsc(5) * rdf_sc(k) *     &
2254                                                               nc(k,j,i)       &
2255                                             )                                 &
2256                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2257                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2258                                          )
2259                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
2260                ENDDO
2261             ENDDO
2262          ENDDO
2263
2264!
2265!--       Calculate tendencies for the next Runge-Kutta step
2266          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2267             IF ( intermediate_timestep_count == 1 )  THEN
2268                DO  i = nxl, nxr
2269                   DO  j = nys, nyn
2270                      DO  k = nzb+1, nzt
2271                         tnc_m(k,j,i) = tend(k,j,i)
2272                      ENDDO
2273                   ENDDO
2274                ENDDO
2275             ELSEIF ( intermediate_timestep_count < &
2276                      intermediate_timestep_count_max )  THEN
2277                DO  i = nxl, nxr
2278                   DO  j = nys, nyn
2279                      DO  k = nzb+1, nzt
2280                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2281                                         + 5.3125_wp * tnc_m(k,j,i)
2282                      ENDDO
2283                   ENDDO
2284                ENDDO
2285             ENDIF
2286          ENDIF
2287
2288          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
2289
2290       ENDIF
2291!
2292!--    If required, calculate prognostic equations for rain water content
2293!--    and rain drop concentration
2294       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
2295
2296          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
2297
2298!
2299!--       Calculate prognostic equation for rain water content
2300          sbt = tsc(2)
2301          IF ( scalar_advec == 'bc-scheme' )  THEN
2302
2303             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2304!
2305!--             Bott-Chlond scheme always uses Euler time step. Thus:
2306                sbt = 1.0_wp
2307             ENDIF
2308             tend = 0.0_wp
2309             CALL advec_s_bc( qr, 'qr' )
2310
2311          ENDIF
2312
2313!
2314!--       qr-tendency terms with no communication
2315          IF ( scalar_advec /= 'bc-scheme' )  THEN
2316             tend = 0.0_wp
2317             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2318                IF ( ws_scheme_sca )  THEN
2319                   CALL advec_s_ws( qr, 'qr' )
2320                ELSE
2321                   CALL advec_s_pw( qr )
2322                ENDIF
2323             ELSE
2324                CALL advec_s_up( qr )
2325             ENDIF
2326          ENDIF
2327
2328          CALL diffusion_s( qr,                                                &
2329                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
2330                            surf_def_h(2)%qrsws,                               &
2331                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
2332                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
2333                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
2334                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
2335                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
2336                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
2337                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
2338
2339!
2340!--       Prognostic equation for rain water content
2341          DO  i = nxl, nxr
2342             DO  j = nys, nyn
2343                DO  k = nzb+1, nzt
2344                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2345                                                      tsc(3) * tqr_m(k,j,i) )  &
2346                                                    - tsc(5) * rdf_sc(k) *     &
2347                                                               qr(k,j,i)       &
2348                                             )                                 &
2349                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2350                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2351                                          )
2352                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
2353                ENDDO
2354             ENDDO
2355          ENDDO
2356
2357!
2358!--       Calculate tendencies for the next Runge-Kutta step
2359          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2360             IF ( intermediate_timestep_count == 1 )  THEN
2361                DO  i = nxl, nxr
2362                   DO  j = nys, nyn
2363                      DO  k = nzb+1, nzt
2364                         tqr_m(k,j,i) = tend(k,j,i)
2365                      ENDDO
2366                   ENDDO
2367                ENDDO
2368             ELSEIF ( intermediate_timestep_count < &
2369                      intermediate_timestep_count_max )  THEN
2370                DO  i = nxl, nxr
2371                   DO  j = nys, nyn
2372                      DO  k = nzb+1, nzt
2373                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2374                                         + 5.3125_wp * tqr_m(k,j,i)
2375                      ENDDO
2376                   ENDDO
2377                ENDDO
2378             ENDIF
2379          ENDIF
2380
2381          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
2382          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
2383
2384!
2385!--       Calculate prognostic equation for rain drop concentration
2386          sbt = tsc(2)
2387          IF ( scalar_advec == 'bc-scheme' )  THEN
2388
2389             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2390!
2391!--             Bott-Chlond scheme always uses Euler time step. Thus:
2392                sbt = 1.0_wp
2393             ENDIF
2394             tend = 0.0_wp
2395             CALL advec_s_bc( nr, 'nr' )
2396
2397          ENDIF
2398
2399!
2400!--       nr-tendency terms with no communication
2401          IF ( scalar_advec /= 'bc-scheme' )  THEN
2402             tend = 0.0_wp
2403             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2404                IF ( ws_scheme_sca )  THEN
2405                   CALL advec_s_ws( nr, 'nr' )
2406                ELSE
2407                   CALL advec_s_pw( nr )
2408                ENDIF
2409             ELSE
2410                CALL advec_s_up( nr )
2411             ENDIF
2412          ENDIF
2413
2414          CALL diffusion_s( nr,                                                &
2415                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
2416                            surf_def_h(2)%nrsws,                               &
2417                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
2418                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
2419                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
2420                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
2421                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
2422                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
2423                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
2424
2425!
2426!--       Prognostic equation for rain drop concentration
2427          DO  i = nxl, nxr
2428             DO  j = nys, nyn
2429                DO  k = nzb+1, nzt
2430                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2431                                                      tsc(3) * tnr_m(k,j,i) )  &
2432                                                    - tsc(5) * rdf_sc(k) *     &
2433                                                               nr(k,j,i)       &
2434                                             )                                 &
2435                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2436                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2437                                          )
2438                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
2439                ENDDO
2440             ENDDO
2441          ENDDO
2442
2443!
2444!--       Calculate tendencies for the next Runge-Kutta step
2445          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2446             IF ( intermediate_timestep_count == 1 )  THEN
2447                DO  i = nxl, nxr
2448                   DO  j = nys, nyn
2449                      DO  k = nzb+1, nzt
2450                         tnr_m(k,j,i) = tend(k,j,i)
2451                      ENDDO
2452                   ENDDO
2453                ENDDO
2454             ELSEIF ( intermediate_timestep_count < &
2455                      intermediate_timestep_count_max )  THEN
2456                DO  i = nxl, nxr
2457                   DO  j = nys, nyn
2458                      DO  k = nzb+1, nzt
2459                         tnr_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2460                                         + 5.3125_wp * tnr_m(k,j,i)
2461                      ENDDO
2462                   ENDDO
2463                ENDDO
2464             ENDIF
2465          ENDIF
2466
2467          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
2468
2469       ENDIF
2470
2471    ENDIF
2472!
2473!-- If required, compute prognostic equation for scalar
2474    IF ( passive_scalar )  THEN
2475
2476       CALL cpu_log( log_point(66), 's-equation', 'start' )
2477
2478!
2479!--    Scalar/q-tendency terms with communication
2480       sbt = tsc(2)
2481       IF ( scalar_advec == 'bc-scheme' )  THEN
2482
2483          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2484!
2485!--          Bott-Chlond scheme always uses Euler time step. Thus:
2486             sbt = 1.0_wp
2487          ENDIF
2488          tend = 0.0_wp
2489          CALL advec_s_bc( s, 's' )
2490
2491       ENDIF
2492
2493!
2494!--    Scalar/q-tendency terms with no communication
2495       IF ( scalar_advec /= 'bc-scheme' )  THEN
2496          tend = 0.0_wp
2497          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2498             IF ( ws_scheme_sca )  THEN
2499                CALL advec_s_ws( s, 's' )
2500             ELSE
2501                CALL advec_s_pw( s )
2502             ENDIF
2503          ELSE
2504             CALL advec_s_up( s )
2505          ENDIF
2506       ENDIF
2507
2508       CALL diffusion_s( s,                                                    &
2509                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
2510                         surf_def_h(2)%ssws,                                   &
2511                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
2512                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
2513                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
2514                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
2515                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
2516                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
2517                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2518
2519!
2520!--    Sink or source of humidity due to canopy elements
2521       IF ( plant_canopy ) CALL pcm_tendency( 7 )
2522
2523!
2524!--    Large scale advection. Not implemented for scalars so far.
2525!        IF ( large_scale_forcing )  THEN
2526!           CALL ls_advec( simulated_time, 'q' )
2527!        ENDIF
2528
2529!
2530!--    Nudging. Not implemented for scalars so far.
2531!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
2532
2533!
2534!--    If required compute influence of large-scale subsidence/ascent.
2535!--    Not implemented for scalars so far.
2536       IF ( large_scale_subsidence  .AND.                                      &
2537            .NOT. use_subsidence_tendencies  .AND.                             &
2538            .NOT. large_scale_forcing )  THEN
2539         CALL subsidence( tend, s, s_init, 3 )
2540       ENDIF
2541
2542       CALL module_interface_actions( 's-tendency' )
2543
2544!
2545!--    Prognostic equation for total water content
2546       DO  i = nxl, nxr
2547          DO  j = nys, nyn
2548             DO  k = nzb+1, nzt
2549                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2550                                                 tsc(3) * ts_m(k,j,i) )        &
2551                                               - tsc(5) * rdf_sc(k) *          &
2552                                                    ( s(k,j,i) - s_init(k) )   &
2553                                        )                                      &
2554                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2555                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2556                                          )
2557                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
2558             ENDDO
2559          ENDDO
2560       ENDDO
2561
2562!
2563!--    Calculate tendencies for the next Runge-Kutta step
2564       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2565          IF ( intermediate_timestep_count == 1 )  THEN
2566             DO  i = nxl, nxr
2567                DO  j = nys, nyn
2568                   DO  k = nzb+1, nzt
2569                      ts_m(k,j,i) = tend(k,j,i)
2570                   ENDDO
2571                ENDDO
2572             ENDDO
2573          ELSEIF ( intermediate_timestep_count < &
2574                   intermediate_timestep_count_max )  THEN
2575             DO  i = nxl, nxr
2576                DO  j = nys, nyn
2577                   DO  k = nzb+1, nzt
2578                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2579                                     + 5.3125_wp * ts_m(k,j,i)
2580                   ENDDO
2581                ENDDO
2582             ENDDO
2583          ENDIF
2584       ENDIF
2585
2586       CALL cpu_log( log_point(66), 's-equation', 'stop' )
2587
2588    ENDIF
2589
2590!
2591!-- Calculate prognostic equations for turbulence closure
2592    CALL tcm_prognostic_equations()
2593!
2594!-- Calculate prognostic equations for all other modules
2595    CALL module_interface_prognostic_equations()
2596
2597!
2598!-- Calculate prognostic equation for chemical quantites
2599    IF ( air_chemistry )  THEN
2600       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2601!
2602!--    Loop over chemical species
2603       DO  lsp = 1, nvar                         
2604          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
2605                                          chem_species(lsp)%conc,              &
2606                                          chem_species(lsp)%tconc_m,           &
2607                                          chem_species(lsp)%conc_pr_init,      &
2608                                          lsp )
2609       ENDDO
2610
2611       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )             
2612    ENDIF   ! Chemicals equations
2613
2614 END SUBROUTINE prognostic_equations_vector
2615
2616
2617 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.