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

Last change on this file since 3887 was 3887, checked in by schwenkel, 2 years ago

bugfix for chemistry_model_mod via introducing module_interface_exchange_horiz

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