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

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

Implemented wtm_actions and moved respective module code into it

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