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

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

Added chem_actions module interface to chemistry_model_mod and moved photolysis_control call into it

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