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

Last change on this file since 3879 was 3879, checked in by knoop, 2 years ago

Moved loop over chem_species into chem_boundary_conds_decycle

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