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

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

Implemented non_transport_physics module interfaces

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