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

Last change on this file since 3738 was 3719, checked in by kanani, 6 years ago

Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)

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