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

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

Including last commit, salsa dependency for advec_ws removed

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