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

Last change on this file since 3929 was 3929, checked in by banzhafs, 3 years ago

Correct/complete module_interface introduction for chemistry model and bug fix in chem_depo subroutine

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