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

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

Added chem_non_transport_physics module interface to chemistry_model_mod and moved respective calls into it

  • Property svn:keywords set to Id
File size: 70.9 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 3878 2019-04-08 19:35:54Z 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 nvar, 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:  call_chem_at_all_substeps, chem_gasphase_on, cs_name,           &
391               deposition_dry, chem_species
392
393    USE chem_gasphase_mod,                                                     &
394        ONLY:  nspec, nvar, spc_names
395
396    USE chemistry_model_mod,                                                   &
397        ONLY:  chem_boundary_conds, chem_prognostic_equations
398
399    USE control_parameters,                                                    &
400        ONLY:  air_chemistry, constant_diffusion,                              &
401               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
402               humidity, intermediate_timestep_count,                          &
403               intermediate_timestep_count_max, large_scale_forcing,           &
404               large_scale_subsidence, neutral, nudging,                       &
405               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
406               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
407               timestep_scheme, tsc, use_subsidence_tendencies,                &
408               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
409               ws_scheme_sca, urban_surface, land_surface,                     &
410               time_since_reference_point, salsa
411
412    USE coriolis_mod,                                                          &
413        ONLY:  coriolis
414
415    USE cpulog,                                                                &
416        ONLY:  cpu_log, log_point, log_point_s
417
418    USE diffusion_s_mod,                                                       &
419        ONLY:  diffusion_s
420
421    USE diffusion_u_mod,                                                       &
422        ONLY:  diffusion_u
423
424    USE diffusion_v_mod,                                                       &
425        ONLY:  diffusion_v
426
427    USE diffusion_w_mod,                                                       &
428        ONLY:  diffusion_w
429
430    USE indices,                                                               &
431        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
432               nzb, nzt, wall_flags_0
433
434    USE kinds
435
436    USE lsf_nudging_mod,                                                       &
437        ONLY:  ls_advec, nudge
438
439    USE module_interface,                                                      &
440        ONLY:  module_interface_actions, &
441               module_interface_non_transport_physics, &
442               module_interface_prognostic_equations
443
444    USE ocean_mod,                                                             &
445        ONLY:  stokes_drift_terms, stokes_force,   &
446               wave_breaking, wave_breaking_term
447
448    USE plant_canopy_model_mod,                                                &
449        ONLY:  cthf, pcm_tendency
450
451#if defined( __rrtmg )
452    USE radiation_model_mod,                                                   &
453        ONLY:  radiation, radiation_tendency,                                  &
454               skip_time_do_radiation
455#endif
456
457    USE salsa_mod,                                                             &
458        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time,        &
459               nbins_aerosol, ncomponents_mass, ngases_salsa,                  &
460               salsa_boundary_conds, salsa_diagnostics, salsa_driver,          &
461               salsa_gas, salsa_gases_from_chem, skip_time_do_salsa
462
463    USE statistics,                                                            &
464        ONLY:  hom
465
466    USE subsidence_mod,                                                        &
467        ONLY:  subsidence
468
469    USE surface_mod,                                                           &
470        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
471                surf_usm_v
472
473    USE turbulence_closure_mod,                                                &
474        ONLY:  tcm_prognostic_equations
475
476    IMPLICIT NONE
477
478    PRIVATE
479    PUBLIC prognostic_equations_cache, prognostic_equations_vector
480
481    INTERFACE prognostic_equations_cache
482       MODULE PROCEDURE prognostic_equations_cache
483    END INTERFACE prognostic_equations_cache
484
485    INTERFACE prognostic_equations_vector
486       MODULE PROCEDURE prognostic_equations_vector
487    END INTERFACE prognostic_equations_vector
488
489
490 CONTAINS
491
492
493!------------------------------------------------------------------------------!
494! Description:
495! ------------
496!> Version with one optimized loop over all equations. It is only allowed to
497!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
498!>
499!> Here the calls of most subroutines are embedded in two DO loops over i and j,
500!> so communication between CPUs is not allowed (does not make sense) within
501!> these loops.
502!>
503!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
504!------------------------------------------------------------------------------!
505
506 SUBROUTINE prognostic_equations_cache
507
508
509    INTEGER(iwp) ::  i                   !<
510    INTEGER(iwp) ::  i_omp_start         !<
511    INTEGER(iwp) ::  ib                  !< index for aerosol size bins (salsa)
512    INTEGER(iwp) ::  ic                  !< index for chemical compounds (salsa)
513    INTEGER(iwp) ::  icc                 !< additional index for chemical compounds (salsa)
514    INTEGER(iwp) ::  ig                  !< index for gaseous compounds (salsa)
515    INTEGER(iwp) ::  j                   !<
516    INTEGER(iwp) ::  k                   !<
517!$  INTEGER(iwp) ::  omp_get_thread_num  !<
518    INTEGER(iwp) ::  tn = 0              !<
519
520    LOGICAL      ::  loop_start          !<
521    INTEGER(iwp) ::  lsp
522    INTEGER(iwp) ::  lsp_usr             !< lsp running index for chem spcs
523
524
525!
526!-- Time measurement can only be performed for the whole set of equations
527    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
528
529    !$OMP PARALLEL PRIVATE (i,j)
530    !$OMP DO
531    DO  i = nxlg, nxrg
532       DO  j = nysg, nyng
533!
534!--       Calculate non transport physics for all other modules
535          CALL module_interface_non_transport_physics( i, j )
536       ENDDO
537    ENDDO
538    !$OMP END PARALLEL
539
540    IF ( air_chemistry )  THEN
541!
542!--    Loop over chemical species       
543       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
544       DO  lsp = 1, nspec
545          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
546          lsp_usr = 1 
547          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
548             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
549
550                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
551                                          chem_species(lsp)%conc_pr_init ) 
552             
553             ENDIF
554             lsp_usr = lsp_usr +1
555          ENDDO
556
557
558       ENDDO
559       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
560
561    ENDIF
562!
563!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
564!-- step. The exchange of ghost points is required after this update of the
565!-- concentrations of aerosol number and mass
566    IF ( salsa )  THEN
567
568       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
569          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
570          THEN
571             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
572             !$OMP PARALLEL PRIVATE (i,j,ib,ic,icc,ig)
573             !$OMP DO
574!
575!--          Call salsa processes
576             DO  i = nxl, nxr
577                DO  j = nys, nyn
578                   CALL salsa_diagnostics( i, j )
579                   CALL salsa_driver( i, j, 3 )
580                   CALL salsa_diagnostics( i, j )
581                ENDDO
582             ENDDO
583
584             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
585             
586             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
587!
588!--          Exchange ghost points and decycle if needed.
589             DO  ib = 1, nbins_aerosol
590                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
591                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
592                                           aerosol_number(ib)%init )
593                DO  ic = 1, ncomponents_mass
594                   icc = ( ic - 1 ) * nbins_aerosol + ib
595                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
596                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
597                                              aerosol_mass(icc)%init )
598                ENDDO
599             ENDDO
600
601             IF ( .NOT. salsa_gases_from_chem )  THEN
602                DO  ig = 1, ngases_salsa
603                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
604                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
605                                              salsa_gas(ig)%init )
606                ENDDO
607             ENDIF
608             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
609             
610             !$OMP END PARALLEL
611             last_salsa_time = time_since_reference_point
612
613          ENDIF
614
615       ENDIF
616
617    ENDIF
618
619!
620!-- Loop over all prognostic equations
621!-- b, c ang g added for SALSA
622    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
623
624    !$  tn = omp_get_thread_num()
625    loop_start = .TRUE.
626
627    !$OMP DO
628    DO  i = nxl, nxr
629
630!
631!--    Store the first loop index. It differs for each thread and is required
632!--    later in advec_ws
633       IF ( loop_start )  THEN
634          loop_start  = .FALSE.
635          i_omp_start = i
636       ENDIF
637
638       DO  j = nys, nyn
639!
640!--       Tendency terms for u-velocity component. Please note, in case of
641!--       non-cyclic boundary conditions the grid point i=0 is excluded from
642!--       the prognostic equations for the u-component.   
643          IF ( i >= nxlu )  THEN
644
645             tend(:,j,i) = 0.0_wp
646             IF ( timestep_scheme(1:5) == 'runge' )  THEN
647                IF ( ws_scheme_mom )  THEN
648                   CALL advec_u_ws( i, j, i_omp_start, tn )
649                ELSE
650                   CALL advec_u_pw( i, j )
651                ENDIF
652             ELSE
653                CALL advec_u_up( i, j )
654             ENDIF
655             CALL diffusion_u( i, j )
656             CALL coriolis( i, j, 1 )
657             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
658                CALL buoyancy( i, j, pt, 1 )
659             ENDIF
660
661!
662!--          Drag by plant canopy
663             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
664
665!
666!--          External pressure gradient
667             IF ( dp_external )  THEN
668                DO  k = dp_level_ind_b+1, nzt
669                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
670                ENDDO
671             ENDIF
672
673!
674!--          Nudging
675             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
676
677!
678!--          Effect of Stokes drift (in ocean mode only)
679             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
680
681             CALL module_interface_actions( i, j, 'u-tendency' )
682!
683!--          Prognostic equation for u-velocity component
684             DO  k = nzb+1, nzt
685
686                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
687                                            ( tsc(2) * tend(k,j,i) +            &
688                                              tsc(3) * tu_m(k,j,i) )            &
689                                            - tsc(5) * rdf(k)                   &
690                                                     * ( u(k,j,i) - u_init(k) ) &
691                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
692                                                   BTEST( wall_flags_0(k,j,i), 1 )&
693                                                 )
694             ENDDO
695
696!
697!--          Add turbulence generated by wave breaking (in ocean mode only)
698             IF ( wave_breaking  .AND.                                         &
699               intermediate_timestep_count == intermediate_timestep_count_max )&
700             THEN
701                CALL wave_breaking_term( i, j, 1 )
702             ENDIF
703
704!
705!--          Calculate tendencies for the next Runge-Kutta step
706             IF ( timestep_scheme(1:5) == 'runge' )  THEN
707                IF ( intermediate_timestep_count == 1 )  THEN
708                   DO  k = nzb+1, nzt
709                      tu_m(k,j,i) = tend(k,j,i)
710                   ENDDO
711                ELSEIF ( intermediate_timestep_count < &
712                         intermediate_timestep_count_max )  THEN
713                   DO  k = nzb+1, nzt
714                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
715                                     + 5.3125_wp * tu_m(k,j,i)
716                   ENDDO
717                ENDIF
718             ENDIF
719
720          ENDIF
721!
722!--       Tendency terms for v-velocity component. Please note, in case of
723!--       non-cyclic boundary conditions the grid point j=0 is excluded from
724!--       the prognostic equations for the v-component. !--       
725          IF ( j >= nysv )  THEN
726
727             tend(:,j,i) = 0.0_wp
728             IF ( timestep_scheme(1:5) == 'runge' )  THEN
729                IF ( ws_scheme_mom )  THEN
730                    CALL advec_v_ws( i, j, i_omp_start, tn )
731                ELSE
732                    CALL advec_v_pw( i, j )
733                ENDIF
734             ELSE
735                CALL advec_v_up( i, j )
736             ENDIF
737             CALL diffusion_v( i, j )
738             CALL coriolis( i, j, 2 )
739
740!
741!--          Drag by plant canopy
742             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
743
744!
745!--          External pressure gradient
746             IF ( dp_external )  THEN
747                DO  k = dp_level_ind_b+1, nzt
748                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
749                ENDDO
750             ENDIF
751
752!
753!--          Nudging
754             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
755
756!
757!--          Effect of Stokes drift (in ocean mode only)
758             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
759
760             CALL module_interface_actions( i, j, 'v-tendency' )
761!
762!--          Prognostic equation for v-velocity component
763             DO  k = nzb+1, nzt
764                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
765                                            ( tsc(2) * tend(k,j,i) +           &
766                                              tsc(3) * tv_m(k,j,i) )           &
767                                            - tsc(5) * rdf(k)                  &
768                                                     * ( v(k,j,i) - v_init(k) )&
769                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
770                                                   BTEST( wall_flags_0(k,j,i), 2 )&
771                                                 )
772             ENDDO
773
774!
775!--          Add turbulence generated by wave breaking (in ocean mode only)
776             IF ( wave_breaking  .AND.                                         &
777               intermediate_timestep_count == intermediate_timestep_count_max )&
778             THEN
779                CALL wave_breaking_term( i, j, 2 )
780             ENDIF
781
782!
783!--          Calculate tendencies for the next Runge-Kutta step
784             IF ( timestep_scheme(1:5) == 'runge' )  THEN
785                IF ( intermediate_timestep_count == 1 )  THEN
786                   DO  k = nzb+1, nzt
787                      tv_m(k,j,i) = tend(k,j,i)
788                   ENDDO
789                ELSEIF ( intermediate_timestep_count < &
790                         intermediate_timestep_count_max )  THEN
791                   DO  k = nzb+1, nzt
792                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
793                                     + 5.3125_wp * tv_m(k,j,i)
794                   ENDDO
795                ENDIF
796             ENDIF
797
798          ENDIF
799
800!
801!--       Tendency terms for w-velocity component
802          tend(:,j,i) = 0.0_wp
803          IF ( timestep_scheme(1:5) == 'runge' )  THEN
804             IF ( ws_scheme_mom )  THEN
805                CALL advec_w_ws( i, j, i_omp_start, tn )
806             ELSE
807                CALL advec_w_pw( i, j )
808             END IF
809          ELSE
810             CALL advec_w_up( i, j )
811          ENDIF
812          CALL diffusion_w( i, j )
813          CALL coriolis( i, j, 3 )
814
815          IF ( .NOT. neutral )  THEN
816             IF ( ocean_mode )  THEN
817                CALL buoyancy( i, j, rho_ocean, 3 )
818             ELSE
819                IF ( .NOT. humidity )  THEN
820                   CALL buoyancy( i, j, pt, 3 )
821                ELSE
822                   CALL buoyancy( i, j, vpt, 3 )
823                ENDIF
824             ENDIF
825          ENDIF
826
827!
828!--       Drag by plant canopy
829          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
830
831!
832!--       Effect of Stokes drift (in ocean mode only)
833          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
834
835          CALL module_interface_actions( i, j, 'w-tendency' )
836!
837!--       Prognostic equation for w-velocity component
838          DO  k = nzb+1, nzt-1
839             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
840                                             ( tsc(2) * tend(k,j,i) +          &
841                                               tsc(3) * tw_m(k,j,i) )          &
842                                             - tsc(5) * rdf(k) * w(k,j,i)      &
843                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
844                                                BTEST( wall_flags_0(k,j,i), 3 )&
845                                              )
846          ENDDO
847
848!
849!--       Calculate tendencies for the next Runge-Kutta step
850          IF ( timestep_scheme(1:5) == 'runge' )  THEN
851             IF ( intermediate_timestep_count == 1 )  THEN
852                DO  k = nzb+1, nzt-1
853                   tw_m(k,j,i) = tend(k,j,i)
854                ENDDO
855             ELSEIF ( intermediate_timestep_count < &
856                      intermediate_timestep_count_max )  THEN
857                DO  k = nzb+1, nzt-1
858                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
859                                  + 5.3125_wp * tw_m(k,j,i)
860                ENDDO
861             ENDIF
862          ENDIF
863
864!
865!--       If required, compute prognostic equation for potential temperature
866          IF ( .NOT. neutral )  THEN
867!
868!--          Tendency terms for potential temperature
869             tend(:,j,i) = 0.0_wp
870             IF ( timestep_scheme(1:5) == 'runge' )  THEN
871                   IF ( ws_scheme_sca )  THEN
872                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
873                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
874                   ELSE
875                      CALL advec_s_pw( i, j, pt )
876                   ENDIF
877             ELSE
878                CALL advec_s_up( i, j, pt )
879             ENDIF
880             CALL diffusion_s( i, j, pt,                                       &
881                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
882                               surf_def_h(2)%shf,                              &
883                               surf_lsm_h%shf,    surf_usm_h%shf,              &
884                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
885                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
886                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
887                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
888                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
889                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
890
891!
892!--          Consideration of heat sources within the plant canopy
893             IF ( plant_canopy  .AND.                                          &
894                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
895                CALL pcm_tendency( i, j, 4 )
896             ENDIF
897
898!
899!--          Large scale advection
900             IF ( large_scale_forcing )  THEN
901                CALL ls_advec( i, j, simulated_time, 'pt' )
902             ENDIF
903
904!
905!--          Nudging
906             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
907
908!
909!--          If required, compute effect of large-scale subsidence/ascent
910             IF ( large_scale_subsidence  .AND.                                &
911                  .NOT. use_subsidence_tendencies )  THEN
912                CALL subsidence( i, j, tend, pt, pt_init, 2 )
913             ENDIF
914
915#if defined( __rrtmg )
916!
917!--          If required, add tendency due to radiative heating/cooling
918             IF ( radiation  .AND.                                             &
919                  simulated_time > skip_time_do_radiation )  THEN
920                CALL radiation_tendency ( i, j, tend )
921             ENDIF
922#endif
923
924             CALL module_interface_actions( i, j, 'pt-tendency' )
925!
926!--          Prognostic equation for potential temperature
927             DO  k = nzb+1, nzt
928                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
929                                                  ( tsc(2) * tend(k,j,i) +     &
930                                                    tsc(3) * tpt_m(k,j,i) )    &
931                                                  - tsc(5)                     &
932                                                  * ( pt(k,j,i) - pt_init(k) ) &
933                                                  * ( rdf_sc(k) + ptdf_x(i)    &
934                                                                + ptdf_y(j) )  &
935                                          )                                    &
936                                       * MERGE( 1.0_wp, 0.0_wp,                &
937                                                BTEST( wall_flags_0(k,j,i), 0 )&
938                                              )
939             ENDDO
940
941!
942!--          Calculate tendencies for the next Runge-Kutta step
943             IF ( timestep_scheme(1:5) == 'runge' )  THEN
944                IF ( intermediate_timestep_count == 1 )  THEN
945                   DO  k = nzb+1, nzt
946                      tpt_m(k,j,i) = tend(k,j,i)
947                   ENDDO
948                ELSEIF ( intermediate_timestep_count < &
949                         intermediate_timestep_count_max )  THEN
950                   DO  k = nzb+1, nzt
951                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
952                                        5.3125_wp * tpt_m(k,j,i)
953                   ENDDO
954                ENDIF
955             ENDIF
956
957          ENDIF
958
959!
960!--       If required, compute prognostic equation for total water content
961          IF ( humidity )  THEN
962
963!
964!--          Tendency-terms for total water content / scalar
965             tend(:,j,i) = 0.0_wp
966             IF ( timestep_scheme(1:5) == 'runge' ) &
967             THEN
968                IF ( ws_scheme_sca )  THEN
969                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
970                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
971                ELSE
972                   CALL advec_s_pw( i, j, q )
973                ENDIF
974             ELSE
975                CALL advec_s_up( i, j, q )
976             ENDIF
977             CALL diffusion_s( i, j, q,                                        &
978                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
979                               surf_def_h(2)%qsws,                             &
980                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
981                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
982                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
983                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
984                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
985                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
986                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
987
988!
989!--          Sink or source of humidity due to canopy elements
990             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
991
992!
993!--          Large scale advection
994             IF ( large_scale_forcing )  THEN
995                CALL ls_advec( i, j, simulated_time, 'q' )
996             ENDIF
997
998!
999!--          Nudging
1000             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
1001
1002!
1003!--          If required compute influence of large-scale subsidence/ascent
1004             IF ( large_scale_subsidence  .AND.                                &
1005                  .NOT. use_subsidence_tendencies )  THEN
1006                CALL subsidence( i, j, tend, q, q_init, 3 )
1007             ENDIF
1008
1009             CALL module_interface_actions( i, j, 'q-tendency' )
1010
1011!
1012!--          Prognostic equation for total water content / scalar
1013             DO  k = nzb+1, nzt
1014                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
1015                                                ( tsc(2) * tend(k,j,i) +       &
1016                                                  tsc(3) * tq_m(k,j,i) )       &
1017                                                - tsc(5) * rdf_sc(k) *         &
1018                                                      ( q(k,j,i) - q_init(k) ) &
1019                                        )                                      &
1020                                       * MERGE( 1.0_wp, 0.0_wp,                &
1021                                                BTEST( wall_flags_0(k,j,i), 0 )&
1022                                              )               
1023                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1024             ENDDO
1025
1026!
1027!--          Calculate tendencies for the next Runge-Kutta step
1028             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1029                IF ( intermediate_timestep_count == 1 )  THEN
1030                   DO  k = nzb+1, nzt
1031                      tq_m(k,j,i) = tend(k,j,i)
1032                   ENDDO
1033                ELSEIF ( intermediate_timestep_count < &
1034                         intermediate_timestep_count_max )  THEN
1035                   DO  k = nzb+1, nzt
1036                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1037                                       5.3125_wp * tq_m(k,j,i)
1038                   ENDDO
1039                ENDIF
1040             ENDIF
1041
1042          ENDIF
1043
1044!
1045!--       If required, compute prognostic equation for scalar
1046          IF ( passive_scalar )  THEN
1047!
1048!--          Tendency-terms for total water content / scalar
1049             tend(:,j,i) = 0.0_wp
1050             IF ( timestep_scheme(1:5) == 'runge' ) &
1051             THEN
1052                IF ( ws_scheme_sca )  THEN
1053                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1054                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1055                ELSE
1056                   CALL advec_s_pw( i, j, s )
1057                ENDIF
1058             ELSE
1059                CALL advec_s_up( i, j, s )
1060             ENDIF
1061             CALL diffusion_s( i, j, s,                                        &
1062                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1063                               surf_def_h(2)%ssws,                             &
1064                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1065                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1066                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1067                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1068                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1069                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1070                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1071
1072!
1073!--          Sink or source of scalar concentration due to canopy elements
1074             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1075
1076!
1077!--          Large scale advection, still need to be extended for scalars
1078!              IF ( large_scale_forcing )  THEN
1079!                 CALL ls_advec( i, j, simulated_time, 's' )
1080!              ENDIF
1081
1082!
1083!--          Nudging, still need to be extended for scalars
1084!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1085
1086!
1087!--          If required compute influence of large-scale subsidence/ascent.
1088!--          Note, the last argument is of no meaning in this case, as it is
1089!--          only used in conjunction with large_scale_forcing, which is to
1090!--          date not implemented for scalars.
1091             IF ( large_scale_subsidence  .AND.                                &
1092                  .NOT. use_subsidence_tendencies  .AND.                       &
1093                  .NOT. large_scale_forcing )  THEN
1094                CALL subsidence( i, j, tend, s, s_init, 3 )
1095             ENDIF
1096
1097             CALL module_interface_actions( i, j, 's-tendency' )
1098
1099!
1100!--          Prognostic equation for scalar
1101             DO  k = nzb+1, nzt
1102                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1103                                            ( tsc(2) * tend(k,j,i) +           &
1104                                              tsc(3) * ts_m(k,j,i) )           &
1105                                            - tsc(5) * rdf_sc(k)               &
1106                                                     * ( s(k,j,i) - s_init(k) )&
1107                                        )                                      &
1108                                       * MERGE( 1.0_wp, 0.0_wp,                &
1109                                                BTEST( wall_flags_0(k,j,i), 0 )&
1110                                              )
1111                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1112             ENDDO
1113
1114!
1115!--          Calculate tendencies for the next Runge-Kutta step
1116             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1117                IF ( intermediate_timestep_count == 1 )  THEN
1118                   DO  k = nzb+1, nzt
1119                      ts_m(k,j,i) = tend(k,j,i)
1120                   ENDDO
1121                ELSEIF ( intermediate_timestep_count < &
1122                         intermediate_timestep_count_max )  THEN
1123                   DO  k = nzb+1, nzt
1124                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1125                                       5.3125_wp * ts_m(k,j,i)
1126                   ENDDO
1127                ENDIF
1128             ENDIF
1129
1130          ENDIF
1131!
1132!--       Calculate prognostic equations for turbulence closure
1133          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1134!
1135!--       Calculate prognostic equations for all other modules
1136          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1137
1138!
1139!--       Calculate prognostic equation for chemical quantites
1140          IF ( air_chemistry )  THEN
1141             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
1142!
1143!--          Loop over chemical species
1144             DO  lsp = 1, nvar                         
1145                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
1146                                     chem_species(lsp)%conc,                   &
1147                                     chem_species(lsp)%tconc_m,                &
1148                                     chem_species(lsp)%conc_pr_init,           &
1149                                     i, j, i_omp_start, tn, lsp,               &
1150                                     chem_species(lsp)%flux_s_cs,              &
1151                                     chem_species(lsp)%diss_s_cs,              &
1152                                     chem_species(lsp)%flux_l_cs,              &
1153                                     chem_species(lsp)%diss_l_cs )       
1154             ENDDO
1155
1156          ENDIF   ! Chemical equations
1157
1158       ENDDO  ! loop over j
1159    ENDDO  ! loop over i
1160    !$OMP END PARALLEL
1161
1162
1163    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1164
1165
1166 END SUBROUTINE prognostic_equations_cache
1167
1168
1169!------------------------------------------------------------------------------!
1170! Description:
1171! ------------
1172!> Version for vector machines
1173!------------------------------------------------------------------------------!
1174
1175 SUBROUTINE prognostic_equations_vector
1176
1177
1178    INTEGER(iwp) ::  i     !<
1179    INTEGER(iwp) ::  ib    !< index for aerosol size bins (salsa)
1180    INTEGER(iwp) ::  ic    !< index for chemical compounds (salsa)
1181    INTEGER(iwp) ::  icc   !< additional index for chemical compounds (salsa)
1182    INTEGER(iwp) ::  ig    !< index for gaseous compounds (salsa)
1183    INTEGER(iwp) ::  j     !<
1184    INTEGER(iwp) ::  k     !<
1185    INTEGER(iwp) ::  lsp   !< running index for chemical species
1186
1187    REAL(wp)     ::  sbt  !<
1188
1189!
1190!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1191!-- step. The exchange of ghost points is required after this update of the
1192!-- concentrations of aerosol number and mass
1193    IF ( salsa )  THEN
1194       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1195          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
1196          THEN
1197             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1198             !$OMP PARALLEL PRIVATE (i,j,ib,ic,icc,ig)
1199             !$OMP DO
1200!
1201!--          Call salsa processes
1202             DO  i = nxl, nxr
1203                DO  j = nys, nyn
1204                   CALL salsa_diagnostics( i, j )
1205                   CALL salsa_driver( i, j, 3 )
1206                   CALL salsa_diagnostics( i, j )
1207                ENDDO
1208             ENDDO
1209
1210             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1211             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1212!
1213!--          Exchange ghost points and decycle if needed.
1214             DO  ib = 1, nbins_aerosol
1215                CALL exchange_horiz( aerosol_number(ib)%conc, nbgp )
1216                CALL salsa_boundary_conds( aerosol_number(ib)%conc,            &
1217                                           aerosol_number(ib)%init )
1218                DO  ic = 1, ncomponents_mass
1219                   icc = ( ic - 1 ) * nbins_aerosol + ib
1220                   CALL exchange_horiz( aerosol_mass(icc)%conc, nbgp )
1221                   CALL salsa_boundary_conds( aerosol_mass(icc)%conc,          &
1222                                              aerosol_mass(icc)%init )
1223                ENDDO
1224             ENDDO
1225             IF ( .NOT. salsa_gases_from_chem )  THEN
1226                DO  ig = 1, ngases_salsa
1227                   CALL exchange_horiz( salsa_gas(ig)%conc, nbgp )
1228                   CALL salsa_boundary_conds( salsa_gas(ig)%conc,              &
1229                                              salsa_gas(ig)%init )
1230                ENDDO
1231             ENDIF
1232             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1233             !$OMP END PARALLEL
1234             last_salsa_time = time_since_reference_point
1235          ENDIF
1236       ENDIF
1237    ENDIF
1238
1239!
1240!-- Calculate non transport physics for all other modules
1241    CALL module_interface_non_transport_physics
1242
1243!
1244!-- u-velocity component
1245    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1246
1247    !$ACC KERNELS PRESENT(tend)
1248    tend = 0.0_wp
1249    !$ACC END KERNELS
1250    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1251       IF ( ws_scheme_mom )  THEN
1252          CALL advec_u_ws
1253       ELSE
1254          CALL advec_u_pw
1255       ENDIF
1256    ELSE
1257       CALL advec_u_up
1258    ENDIF
1259    CALL diffusion_u
1260    CALL coriolis( 1 )
1261    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1262       CALL buoyancy( pt, 1 )
1263    ENDIF
1264
1265!
1266!-- Drag by plant canopy
1267    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1268
1269!
1270!-- External pressure gradient
1271    IF ( dp_external )  THEN
1272       DO  i = nxlu, nxr
1273          DO  j = nys, nyn
1274             DO  k = dp_level_ind_b+1, nzt
1275                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1276             ENDDO
1277          ENDDO
1278       ENDDO
1279    ENDIF
1280
1281!
1282!-- Nudging
1283    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1284
1285!
1286!-- Effect of Stokes drift (in ocean mode only)
1287    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1288
1289    CALL module_interface_actions( 'u-tendency' )
1290
1291!
1292!-- Prognostic equation for u-velocity component
1293    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1294    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1295    !$ACC PRESENT(tsc(2:5)) &
1296    !$ACC PRESENT(u_p)
1297    DO  i = nxlu, nxr
1298       DO  j = nys, nyn
1299          DO  k = nzb+1, nzt
1300             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1301                                                 tsc(3) * tu_m(k,j,i) )          &
1302                                               - tsc(5) * rdf(k) *               &
1303                                                        ( u(k,j,i) - u_init(k) ) &
1304                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1305                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1306                                              )
1307          ENDDO
1308       ENDDO
1309    ENDDO
1310
1311!
1312!-- Add turbulence generated by wave breaking (in ocean mode only)
1313    IF ( wave_breaking  .AND.                                                  &
1314         intermediate_timestep_count == intermediate_timestep_count_max )      &
1315    THEN
1316       CALL wave_breaking_term( 1 )
1317    ENDIF
1318
1319!
1320!-- Calculate tendencies for the next Runge-Kutta step
1321    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1322       IF ( intermediate_timestep_count == 1 )  THEN
1323          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1324          !$ACC PRESENT(tend, tu_m)
1325          DO  i = nxlu, nxr
1326             DO  j = nys, nyn
1327                DO  k = nzb+1, nzt
1328                   tu_m(k,j,i) = tend(k,j,i)
1329                ENDDO
1330             ENDDO
1331          ENDDO
1332       ELSEIF ( intermediate_timestep_count < &
1333                intermediate_timestep_count_max )  THEN
1334          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1335          !$ACC PRESENT(tend, tu_m)
1336          DO  i = nxlu, nxr
1337             DO  j = nys, nyn
1338                DO  k = nzb+1, nzt
1339                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1340                                   + 5.3125_wp * tu_m(k,j,i)
1341                ENDDO
1342             ENDDO
1343          ENDDO
1344       ENDIF
1345    ENDIF
1346
1347    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1348
1349!
1350!-- v-velocity component
1351    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1352
1353    !$ACC KERNELS PRESENT(tend)
1354    tend = 0.0_wp
1355    !$ACC END KERNELS
1356    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1357       IF ( ws_scheme_mom )  THEN
1358          CALL advec_v_ws
1359       ELSE
1360          CALL advec_v_pw
1361       END IF
1362    ELSE
1363       CALL advec_v_up
1364    ENDIF
1365    CALL diffusion_v
1366    CALL coriolis( 2 )
1367
1368!
1369!-- Drag by plant canopy
1370    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1371
1372!
1373!-- External pressure gradient
1374    IF ( dp_external )  THEN
1375       DO  i = nxl, nxr
1376          DO  j = nysv, nyn
1377             DO  k = dp_level_ind_b+1, nzt
1378                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1379             ENDDO
1380          ENDDO
1381       ENDDO
1382    ENDIF
1383
1384!
1385!-- Nudging
1386    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1387
1388!
1389!-- Effect of Stokes drift (in ocean mode only)
1390    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1391
1392    CALL module_interface_actions( 'v-tendency' )
1393
1394!
1395!-- Prognostic equation for v-velocity component
1396    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1397    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1398    !$ACC PRESENT(tsc(2:5)) &
1399    !$ACC PRESENT(v_p)
1400    DO  i = nxl, nxr
1401       DO  j = nysv, nyn
1402          DO  k = nzb+1, nzt
1403             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1404                                                 tsc(3) * tv_m(k,j,i) )        &
1405                                               - tsc(5) * rdf(k) *             &
1406                                                      ( v(k,j,i) - v_init(k) ) &
1407                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1408                                                BTEST( wall_flags_0(k,j,i), 2 )&
1409                                              )
1410          ENDDO
1411       ENDDO
1412    ENDDO
1413
1414!
1415!-- Add turbulence generated by wave breaking (in ocean mode only)
1416    IF ( wave_breaking  .AND.                                                  &
1417         intermediate_timestep_count == intermediate_timestep_count_max )      &
1418    THEN
1419       CALL wave_breaking_term( 2 )
1420    ENDIF
1421
1422!
1423!-- Calculate tendencies for the next Runge-Kutta step
1424    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1425       IF ( intermediate_timestep_count == 1 )  THEN
1426          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1427          !$ACC PRESENT(tend, tv_m)
1428          DO  i = nxl, nxr
1429             DO  j = nysv, nyn
1430                DO  k = nzb+1, nzt
1431                   tv_m(k,j,i) = tend(k,j,i)
1432                ENDDO
1433             ENDDO
1434          ENDDO
1435       ELSEIF ( intermediate_timestep_count < &
1436                intermediate_timestep_count_max )  THEN
1437          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1438          !$ACC PRESENT(tend, tv_m)
1439          DO  i = nxl, nxr
1440             DO  j = nysv, nyn
1441                DO  k = nzb+1, nzt
1442                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1443                                  + 5.3125_wp * tv_m(k,j,i)
1444                ENDDO
1445             ENDDO
1446          ENDDO
1447       ENDIF
1448    ENDIF
1449
1450    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1451
1452!
1453!-- w-velocity component
1454    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1455
1456    !$ACC KERNELS PRESENT(tend)
1457    tend = 0.0_wp
1458    !$ACC END KERNELS
1459    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1460       IF ( ws_scheme_mom )  THEN
1461          CALL advec_w_ws
1462       ELSE
1463          CALL advec_w_pw
1464       ENDIF
1465    ELSE
1466       CALL advec_w_up
1467    ENDIF
1468    CALL diffusion_w
1469    CALL coriolis( 3 )
1470
1471    IF ( .NOT. neutral )  THEN
1472       IF ( ocean_mode )  THEN
1473          CALL buoyancy( rho_ocean, 3 )
1474       ELSE
1475          IF ( .NOT. humidity )  THEN
1476             CALL buoyancy( pt, 3 )
1477          ELSE
1478             CALL buoyancy( vpt, 3 )
1479          ENDIF
1480       ENDIF
1481    ENDIF
1482
1483!
1484!-- Drag by plant canopy
1485    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1486
1487!
1488!-- Effect of Stokes drift (in ocean mode only)
1489    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1490
1491    CALL module_interface_actions( 'w-tendency' )
1492
1493!
1494!-- Prognostic equation for w-velocity component
1495    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1496    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1497    !$ACC PRESENT(tsc(2:5)) &
1498    !$ACC PRESENT(w_p)
1499    DO  i = nxl, nxr
1500       DO  j = nys, nyn
1501          DO  k = nzb+1, nzt-1
1502             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1503                                                 tsc(3) * tw_m(k,j,i) )        &
1504                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1505                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1506                                                BTEST( wall_flags_0(k,j,i), 3 )&
1507                                              )
1508          ENDDO
1509       ENDDO
1510    ENDDO
1511
1512!
1513!-- Calculate tendencies for the next Runge-Kutta step
1514    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1515       IF ( intermediate_timestep_count == 1 )  THEN
1516          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1517          !$ACC PRESENT(tend, tw_m)
1518          DO  i = nxl, nxr
1519             DO  j = nys, nyn
1520                DO  k = nzb+1, nzt-1
1521                   tw_m(k,j,i) = tend(k,j,i)
1522                ENDDO
1523             ENDDO
1524          ENDDO
1525       ELSEIF ( intermediate_timestep_count < &
1526                intermediate_timestep_count_max )  THEN
1527          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1528          !$ACC PRESENT(tend, tw_m)
1529          DO  i = nxl, nxr
1530             DO  j = nys, nyn
1531                DO  k = nzb+1, nzt-1
1532                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1533                                  + 5.3125_wp * tw_m(k,j,i)
1534                ENDDO
1535             ENDDO
1536          ENDDO
1537       ENDIF
1538    ENDIF
1539
1540    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1541
1542
1543!
1544!-- If required, compute prognostic equation for potential temperature
1545    IF ( .NOT. neutral )  THEN
1546
1547       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1548
1549!
1550!--    pt-tendency terms with communication
1551       sbt = tsc(2)
1552       IF ( scalar_advec == 'bc-scheme' )  THEN
1553
1554          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1555!
1556!--          Bott-Chlond scheme always uses Euler time step. Thus:
1557             sbt = 1.0_wp
1558          ENDIF
1559          tend = 0.0_wp
1560          CALL advec_s_bc( pt, 'pt' )
1561
1562       ENDIF
1563
1564!
1565!--    pt-tendency terms with no communication
1566       IF ( scalar_advec /= 'bc-scheme' )  THEN
1567          !$ACC KERNELS PRESENT(tend)
1568          tend = 0.0_wp
1569          !$ACC END KERNELS
1570          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1571             IF ( ws_scheme_sca )  THEN
1572                CALL advec_s_ws( pt, 'pt' )
1573             ELSE
1574                CALL advec_s_pw( pt )
1575             ENDIF
1576          ELSE
1577             CALL advec_s_up( pt )
1578          ENDIF
1579       ENDIF
1580
1581       CALL diffusion_s( pt,                                                   &
1582                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1583                         surf_def_h(2)%shf,                                    &
1584                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1585                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1586                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1587                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1588                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1589                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1590                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1591
1592!
1593!--    Consideration of heat sources within the plant canopy
1594       IF ( plant_canopy  .AND.                                          &
1595            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1596          CALL pcm_tendency( 4 )
1597       ENDIF
1598
1599!
1600!--    Large scale advection
1601       IF ( large_scale_forcing )  THEN
1602          CALL ls_advec( simulated_time, 'pt' )
1603       ENDIF
1604
1605!
1606!--    Nudging
1607       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1608
1609!
1610!--    If required compute influence of large-scale subsidence/ascent
1611       IF ( large_scale_subsidence  .AND.                                      &
1612            .NOT. use_subsidence_tendencies )  THEN
1613          CALL subsidence( tend, pt, pt_init, 2 )
1614       ENDIF
1615
1616#if defined( __rrtmg )
1617!
1618!--    If required, add tendency due to radiative heating/cooling
1619       IF ( radiation  .AND.                                                   &
1620            simulated_time > skip_time_do_radiation )  THEN
1621            CALL radiation_tendency ( tend )
1622       ENDIF
1623#endif
1624
1625       CALL module_interface_actions( 'pt-tendency' )
1626
1627!
1628!--    Prognostic equation for potential temperature
1629       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1630       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1631       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1632       !$ACC PRESENT(tsc(3:5)) &
1633       !$ACC PRESENT(pt_p)
1634       DO  i = nxl, nxr
1635          DO  j = nys, nyn
1636             DO  k = nzb+1, nzt
1637                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1638                                                   tsc(3) * tpt_m(k,j,i) )     &
1639                                                 - tsc(5) *                    &
1640                                                   ( pt(k,j,i) - pt_init(k) ) *&
1641                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1642                                          )                                    &
1643                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1644                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1645                                          )
1646             ENDDO
1647          ENDDO
1648       ENDDO
1649!
1650!--    Calculate tendencies for the next Runge-Kutta step
1651       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1652          IF ( intermediate_timestep_count == 1 )  THEN
1653             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1654             !$ACC PRESENT(tend, tpt_m)
1655             DO  i = nxl, nxr
1656                DO  j = nys, nyn
1657                   DO  k = nzb+1, nzt
1658                      tpt_m(k,j,i) = tend(k,j,i)
1659                   ENDDO
1660                ENDDO
1661             ENDDO
1662          ELSEIF ( intermediate_timestep_count < &
1663                   intermediate_timestep_count_max )  THEN
1664             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1665             !$ACC PRESENT(tend, tpt_m)
1666             DO  i = nxl, nxr
1667                DO  j = nys, nyn
1668                   DO  k = nzb+1, nzt
1669                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1670                                        5.3125_wp * tpt_m(k,j,i)
1671                   ENDDO
1672                ENDDO
1673             ENDDO
1674          ENDIF
1675       ENDIF
1676
1677       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1678
1679    ENDIF
1680
1681!
1682!-- If required, compute prognostic equation for total water content
1683    IF ( humidity )  THEN
1684
1685       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1686
1687!
1688!--    Scalar/q-tendency terms with communication
1689       sbt = tsc(2)
1690       IF ( scalar_advec == 'bc-scheme' )  THEN
1691
1692          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1693!
1694!--          Bott-Chlond scheme always uses Euler time step. Thus:
1695             sbt = 1.0_wp
1696          ENDIF
1697          tend = 0.0_wp
1698          CALL advec_s_bc( q, 'q' )
1699
1700       ENDIF
1701
1702!
1703!--    Scalar/q-tendency terms with no communication
1704       IF ( scalar_advec /= 'bc-scheme' )  THEN
1705          tend = 0.0_wp
1706          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1707             IF ( ws_scheme_sca )  THEN
1708                CALL advec_s_ws( q, 'q' )
1709             ELSE
1710                CALL advec_s_pw( q )
1711             ENDIF
1712          ELSE
1713             CALL advec_s_up( q )
1714          ENDIF
1715       ENDIF
1716
1717       CALL diffusion_s( q,                                                    &
1718                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1719                         surf_def_h(2)%qsws,                                   &
1720                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1721                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1722                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1723                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1724                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1725                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1726                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1727
1728!
1729!--    Sink or source of humidity due to canopy elements
1730       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1731
1732!
1733!--    Large scale advection
1734       IF ( large_scale_forcing )  THEN
1735          CALL ls_advec( simulated_time, 'q' )
1736       ENDIF
1737
1738!
1739!--    Nudging
1740       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1741
1742!
1743!--    If required compute influence of large-scale subsidence/ascent
1744       IF ( large_scale_subsidence  .AND.                                      &
1745            .NOT. use_subsidence_tendencies )  THEN
1746         CALL subsidence( tend, q, q_init, 3 )
1747       ENDIF
1748
1749       CALL module_interface_actions( 'q-tendency' )
1750
1751!
1752!--    Prognostic equation for total water content
1753       DO  i = nxl, nxr
1754          DO  j = nys, nyn
1755             DO  k = nzb+1, nzt
1756                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1757                                                 tsc(3) * tq_m(k,j,i) )        &
1758                                               - tsc(5) * rdf_sc(k) *          &
1759                                                      ( q(k,j,i) - q_init(k) ) &
1760                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1761                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1762                                                 )
1763                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1764             ENDDO
1765          ENDDO
1766       ENDDO
1767
1768!
1769!--    Calculate tendencies for the next Runge-Kutta step
1770       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1771          IF ( intermediate_timestep_count == 1 )  THEN
1772             DO  i = nxl, nxr
1773                DO  j = nys, nyn
1774                   DO  k = nzb+1, nzt
1775                      tq_m(k,j,i) = tend(k,j,i)
1776                   ENDDO
1777                ENDDO
1778             ENDDO
1779          ELSEIF ( intermediate_timestep_count < &
1780                   intermediate_timestep_count_max )  THEN
1781             DO  i = nxl, nxr
1782                DO  j = nys, nyn
1783                   DO  k = nzb+1, nzt
1784                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1785                                     + 5.3125_wp * tq_m(k,j,i)
1786                   ENDDO
1787                ENDDO
1788             ENDDO
1789          ENDIF
1790       ENDIF
1791
1792       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1793
1794    ENDIF
1795!
1796!-- If required, compute prognostic equation for scalar
1797    IF ( passive_scalar )  THEN
1798
1799       CALL cpu_log( log_point(66), 's-equation', 'start' )
1800
1801!
1802!--    Scalar/q-tendency terms with communication
1803       sbt = tsc(2)
1804       IF ( scalar_advec == 'bc-scheme' )  THEN
1805
1806          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1807!
1808!--          Bott-Chlond scheme always uses Euler time step. Thus:
1809             sbt = 1.0_wp
1810          ENDIF
1811          tend = 0.0_wp
1812          CALL advec_s_bc( s, 's' )
1813
1814       ENDIF
1815
1816!
1817!--    Scalar/q-tendency terms with no communication
1818       IF ( scalar_advec /= 'bc-scheme' )  THEN
1819          tend = 0.0_wp
1820          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1821             IF ( ws_scheme_sca )  THEN
1822                CALL advec_s_ws( s, 's' )
1823             ELSE
1824                CALL advec_s_pw( s )
1825             ENDIF
1826          ELSE
1827             CALL advec_s_up( s )
1828          ENDIF
1829       ENDIF
1830
1831       CALL diffusion_s( s,                                                    &
1832                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
1833                         surf_def_h(2)%ssws,                                   &
1834                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
1835                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
1836                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
1837                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
1838                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
1839                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
1840                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1841
1842!
1843!--    Sink or source of humidity due to canopy elements
1844       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1845
1846!
1847!--    Large scale advection. Not implemented for scalars so far.
1848!        IF ( large_scale_forcing )  THEN
1849!           CALL ls_advec( simulated_time, 'q' )
1850!        ENDIF
1851
1852!
1853!--    Nudging. Not implemented for scalars so far.
1854!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1855
1856!
1857!--    If required compute influence of large-scale subsidence/ascent.
1858!--    Not implemented for scalars so far.
1859       IF ( large_scale_subsidence  .AND.                                      &
1860            .NOT. use_subsidence_tendencies  .AND.                             &
1861            .NOT. large_scale_forcing )  THEN
1862         CALL subsidence( tend, s, s_init, 3 )
1863       ENDIF
1864
1865       CALL module_interface_actions( 's-tendency' )
1866
1867!
1868!--    Prognostic equation for total water content
1869       DO  i = nxl, nxr
1870          DO  j = nys, nyn
1871             DO  k = nzb+1, nzt
1872                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1873                                                 tsc(3) * ts_m(k,j,i) )        &
1874                                               - tsc(5) * rdf_sc(k) *          &
1875                                                    ( s(k,j,i) - s_init(k) )   &
1876                                        )                                      &
1877                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1878                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1879                                          )
1880                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1881             ENDDO
1882          ENDDO
1883       ENDDO
1884
1885!
1886!--    Calculate tendencies for the next Runge-Kutta step
1887       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1888          IF ( intermediate_timestep_count == 1 )  THEN
1889             DO  i = nxl, nxr
1890                DO  j = nys, nyn
1891                   DO  k = nzb+1, nzt
1892                      ts_m(k,j,i) = tend(k,j,i)
1893                   ENDDO
1894                ENDDO
1895             ENDDO
1896          ELSEIF ( intermediate_timestep_count < &
1897                   intermediate_timestep_count_max )  THEN
1898             DO  i = nxl, nxr
1899                DO  j = nys, nyn
1900                   DO  k = nzb+1, nzt
1901                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1902                                     + 5.3125_wp * ts_m(k,j,i)
1903                   ENDDO
1904                ENDDO
1905             ENDDO
1906          ENDIF
1907       ENDIF
1908
1909       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1910
1911    ENDIF
1912
1913!
1914!-- Calculate prognostic equations for turbulence closure
1915    CALL tcm_prognostic_equations()
1916!
1917!-- Calculate prognostic equations for all other modules
1918    CALL module_interface_prognostic_equations()
1919
1920!
1921!-- Calculate prognostic equation for chemical quantites
1922    IF ( air_chemistry )  THEN
1923       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
1924!
1925!--    Loop over chemical species
1926       DO  lsp = 1, nvar                         
1927          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
1928                                          chem_species(lsp)%conc,              &
1929                                          chem_species(lsp)%tconc_m,           &
1930                                          chem_species(lsp)%conc_pr_init,      &
1931                                          lsp )
1932       ENDDO
1933
1934       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )             
1935    ENDIF   ! Chemicals equations
1936
1937 END SUBROUTINE prognostic_equations_vector
1938
1939
1940 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.