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

Last change on this file since 3837 was 3837, checked in by knoop, 3 years ago

Added module_interface for prognostic_equations

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