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

Last change on this file since 3870 was 3870, checked in by knoop, 5 years ago

Moving prognostic equations of bcm into bulk_cloud_model_mod

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