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

Last change on this file since 3779 was 3771, checked in by raasch, 6 years ago

rrtmg preprocessor for directives moved/added, save attribute added to temporary pointers to avoid compiler warnings about outlived pointer targets, statement added to avoid compiler warning about unused variable

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