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

Last change on this file since 3911 was 3899, checked in by monakurppa, 5 years ago

corrected the OpenMP implementation for salsa and some minor bugs in salsa_mod

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