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

Last change on this file since 3768 was 3761, checked in by raasch, 2 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

  • Property svn:keywords set to Id
File size: 104.3 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 3761 2019-02-25 15:31:42Z raasch $
27! unused variable removed
28!
29! 3719 2019-02-06 13:10:18Z kanani
30! Cleaned up chemistry cpu measurements
31!
32! 3684 2019-01-20 20:20:58Z knoop
33! OpenACC port for SPEC
34!
35! 3589 2018-11-30 15:09:51Z suehring
36! Move the control parameter "salsa" from salsa_mod to control_parameters
37! (M. Kurppa)
38!
39! 3582 2018-11-29 19:16:36Z suehring
40! Implementation of a new aerosol module salsa.
41! Remove cpu-logs from i,j loop in cache version.
42!
43! 3458 2018-10-30 14:51:23Z kanani
44! remove duplicate USE chem_modules
45! from chemistry branch r3443, banzhafs, basit:
46! chem_depo call introduced
47! code added for decycling chemistry
48!
49! 3386 2018-10-19 16:28:22Z gronemeier
50! Renamed tcm_prognostic to tcm_prognostic_equations
51!
52! 3355 2018-10-16 14:03:34Z knoop
53! (from branch resler)
54! Fix for chemistry call
55!
56! 3302 2018-10-03 02:39:40Z raasch
57! Stokes drift + wave breaking term added
58!
59! 3298 2018-10-02 12:21:11Z kanani
60! Code added for decycling chemistry (basit)
61!
62! 3294 2018-10-01 02:37:10Z raasch
63! changes concerning modularization of ocean option
64!
65! 3274 2018-09-24 15:42:55Z knoop
66! Modularization of all bulk cloud physics code components
67!
68! 3241 2018-09-12 15:02:00Z raasch
69! omp_get_thread_num now declared in omp directive
70!
71! 3183 2018-07-27 14:25:55Z suehring
72! Remove unused variables from USE statements
73!
74! 3182 2018-07-27 13:36:03Z suehring
75! Revise recent bugfix for nesting
76!
77! 3021 2018-05-16 08:14:20Z maronga
78! Bugfix in IF clause for nesting
79!
80! 3014 2018-05-09 08:42:38Z maronga
81! Fixed a bug in the IF condition to call pcm_tendency in case of
82! potential temperature
83!
84! 2815 2018-02-19 11:29:57Z kanani
85! Rename chem_tendency to chem_prognostic_equations,
86! implement vector version for air chemistry
87!
88! 2766 2018-01-22 17:17:47Z kanani
89! Removed preprocessor directive __chem
90!
91! 2746 2018-01-15 12:06:04Z suehring
92! Move flag plant canopy to modules
93!
94! 2719 2018-01-02 09:02:06Z maronga
95! Bugfix for last change.
96!
97! 2718 2018-01-02 08:49:38Z maronga
98! Corrected "Former revisions" section
99!
100! 2696 2017-12-14 17:12:51Z kanani
101! - Change in file header (GPL part)
102! - Moved TKE equation to tcm_prognostic (TG)
103! - Added switch for chemical reactions (RF, FK)
104! - Implementation of chemistry module (RF, BK, FK)
105!
106! 2563 2017-10-19 15:36:10Z Giersch
107! Variable wind_turbine moved to module control_parameters
108!
109! 2320 2017-07-21 12:47:43Z suehring
110! Modularize large-scale forcing and nudging
111!
112! 2292 2017-06-20 09:51:42Z schwenkel
113! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
114! includes two more prognostic equations for cloud drop concentration (nc) 
115! and cloud water content (qc).
116!
117! 2261 2017-06-08 14:25:57Z raasch
118! bugfix for r2232: openmp directives removed
119!
120! 2233 2017-05-30 18:08:54Z suehring
121!
122! 2232 2017-05-30 17:47:52Z suehring
123! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
124! which is realized directly in diffusion_s now.
125!
126! 2192 2017-03-22 04:14:10Z raasch
127! Bugfix for misplaced and missing openMP directives from r2155
128!
129! 2155 2017-02-21 09:57:40Z hoffmann
130! Bugfix in the calculation of microphysical quantities on ghost points.
131!
132! 2118 2017-01-17 16:38:49Z raasch
133! OpenACC version of subroutine removed
134!
135! 2031 2016-10-21 15:11:58Z knoop
136! renamed variable rho to rho_ocean
137!
138! 2011 2016-09-19 17:29:57Z kanani
139! Flag urban_surface is now defined in module control_parameters.
140!
141! 2007 2016-08-24 15:47:17Z kanani
142! Added pt tendency calculation based on energy balance at urban surfaces
143! (new urban surface model)
144!
145! 2000 2016-08-20 18:09:15Z knoop
146! Forced header and separation lines into 80 columns
147!
148! 1976 2016-07-27 13:28:04Z maronga
149! Simplied calls to radiation model
150!
151! 1960 2016-07-12 16:34:24Z suehring
152! Separate humidity and passive scalar
153!
154! 1914 2016-05-26 14:44:07Z witha
155! Added calls for wind turbine model
156!
157! 1873 2016-04-18 14:50:06Z maronga
158! Module renamed (removed _mod)
159!
160! 1850 2016-04-08 13:29:27Z maronga
161! Module renamed
162!
163! 1826 2016-04-07 12:01:39Z maronga
164! Renamed canopy model calls.
165!
166! 1822 2016-04-07 07:49:42Z hoffmann
167! Kessler microphysics scheme moved to microphysics.
168!
169! 1757 2016-02-22 15:49:32Z maronga
170!
171! 1691 2015-10-26 16:17:44Z maronga
172! Added optional model spin-up without radiation / land surface model calls.
173! Formatting corrections.
174!
175! 1682 2015-10-07 23:56:08Z knoop
176! Code annotations made doxygen readable
177!
178! 1585 2015-04-30 07:05:52Z maronga
179! Added call for temperature tendency calculation due to radiative flux divergence
180!
181! 1517 2015-01-07 19:12:25Z hoffmann
182! advec_s_bc_mod addded, since advec_s_bc is now a module
183!
184! 1484 2014-10-21 10:53:05Z kanani
185! Changes due to new module structure of the plant canopy model:
186! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
187! Removed double-listing of use_upstream_for_tke in ONLY-list of module
188! control_parameters
189!
190! 1409 2014-05-23 12:11:32Z suehring
191! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
192! This ensures that left-hand side fluxes are also calculated for nxl in that
193! case, even though the solution at nxl is overwritten in boundary_conds()
194!
195! 1398 2014-05-07 11:15:00Z heinze
196! Rayleigh-damping for horizontal velocity components changed: instead of damping
197! against ug and vg, damping against u_init and v_init is used to allow for a
198! homogenized treatment in case of nudging
199!
200! 1380 2014-04-28 12:40:45Z heinze
201! Change order of calls for scalar prognostic quantities:
202! ls_advec -> nudging -> subsidence since initial profiles
203!
204! 1374 2014-04-25 12:55:07Z raasch
205! missing variables added to ONLY lists
206!
207! 1365 2014-04-22 15:03:56Z boeske
208! Calls of ls_advec for large scale advection added,
209! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
210! new argument ls_index added to the calls of subsidence
211! +ls_index
212!
213! 1361 2014-04-16 15:17:48Z hoffmann
214! Two-moment microphysics moved to the start of prognostic equations. This makes
215! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
216! Additionally, it is allowed to call the microphysics just once during the time
217! step (not at each sub-time step).
218!
219! Two-moment cloud physics added for vector and accelerator optimization.
220!
221! 1353 2014-04-08 15:21:23Z heinze
222! REAL constants provided with KIND-attribute
223!
224! 1337 2014-03-25 15:11:48Z heinze
225! Bugfix: REAL constants provided with KIND-attribute
226!
227! 1332 2014-03-25 11:59:43Z suehring
228! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
229!
230! 1330 2014-03-24 17:29:32Z suehring
231! In case of SGS-particle velocity advection of TKE is also allowed with
232! dissipative 5th-order scheme.
233!
234! 1320 2014-03-20 08:40:49Z raasch
235! ONLY-attribute added to USE-statements,
236! kind-parameters added to all INTEGER and REAL declaration statements,
237! kinds are defined in new module kinds,
238! old module precision_kind is removed,
239! revision history before 2012 removed,
240! comment fields (!:) to be used for variable explanations added to
241! all variable declaration statements
242!
243! 1318 2014-03-17 13:35:16Z raasch
244! module interfaces removed
245!
246! 1257 2013-11-08 15:18:40Z raasch
247! openacc loop vector clauses removed, independent clauses added
248!
249! 1246 2013-11-01 08:59:45Z heinze
250! enable nudging also for accelerator version
251!
252! 1241 2013-10-30 11:36:58Z heinze
253! usage of nudging enabled (so far not implemented for accelerator version)
254!
255! 1179 2013-06-14 05:57:58Z raasch
256! two arguments removed from routine buoyancy, ref_state updated on device
257!
258! 1128 2013-04-12 06:19:32Z raasch
259! those parts requiring global communication moved to time_integration,
260! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
261! j_north
262!
263! 1115 2013-03-26 18:16:16Z hoffmann
264! optimized cloud physics: calculation of microphysical tendencies transfered
265! to microphysics.f90; qr and nr are only calculated if precipitation is required
266!
267! 1111 2013-03-08 23:54:10Z raasch
268! update directives for prognostic quantities removed
269!
270! 1106 2013-03-04 05:31:38Z raasch
271! small changes in code formatting
272!
273! 1092 2013-02-02 11:24:22Z raasch
274! unused variables removed
275!
276! 1053 2012-11-13 17:11:03Z hoffmann
277! implementation of two new prognostic equations for rain drop concentration (nr)
278! and rain water content (qr)
279!
280! currently, only available for cache loop optimization
281!
282! 1036 2012-10-22 13:43:42Z raasch
283! code put under GPL (PALM 3.9)
284!
285! 1019 2012-09-28 06:46:45Z raasch
286! non-optimized version of prognostic_equations removed
287!
288! 1015 2012-09-27 09:23:24Z raasch
289! new branch prognostic_equations_acc
290! OpenACC statements added + code changes required for GPU optimization
291!
292! 1001 2012-09-13 14:08:46Z raasch
293! all actions concerning leapfrog- and upstream-spline-scheme removed
294!
295! 978 2012-08-09 08:28:32Z fricke
296! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
297! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
298! boundary in case of non-cyclic lateral boundaries
299! Bugfix: first thread index changes for WS-scheme at the inflow
300!
301! 940 2012-07-09 14:31:00Z raasch
302! temperature equation can be switched off
303!
304! Revision 1.1  2000/04/13 14:56:27  schroeter
305! Initial revision
306!
307!
308! Description:
309! ------------
310!> Solving the prognostic equations.
311!------------------------------------------------------------------------------!
312 MODULE prognostic_equations_mod
313
314
315    USE advec_s_bc_mod,                                                        &
316        ONLY:  advec_s_bc
317
318    USE advec_s_pw_mod,                                                        &
319        ONLY:  advec_s_pw
320
321    USE advec_s_up_mod,                                                        &
322        ONLY:  advec_s_up
323
324    USE advec_u_pw_mod,                                                        &
325        ONLY:  advec_u_pw
326
327    USE advec_u_up_mod,                                                        &
328        ONLY:  advec_u_up
329
330    USE advec_v_pw_mod,                                                        &
331        ONLY:  advec_v_pw
332
333    USE advec_v_up_mod,                                                        &
334        ONLY:  advec_v_up
335
336    USE advec_w_pw_mod,                                                        &
337        ONLY:  advec_w_pw
338
339    USE advec_w_up_mod,                                                        &
340        ONLY:  advec_w_up
341
342    USE advec_ws,                                                              &
343        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
344
345    USE arrays_3d,                                                             &
346        ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qc, &
347               diss_l_qr, diss_l_s, diss_l_sa, diss_s_e, diss_s_nc, diss_s_nr, &
348               diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, diss_s_sa, &
349               e, e_p, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,    &
350               flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, flux_l_nc, &
351               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
352               flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
353               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
354               ref_state, rho_ocean, s, s_init, s_p, tend, te_m, tnc_m,        &
355               tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tu_m, tv_m, tw_m, u,    &
356               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
357
358    USE bulk_cloud_model_mod,                                                  &
359        ONLY:  call_microphysics_at_all_substeps, bulk_cloud_model,            &
360               bcm_actions, microphysics_sat_adjust,                           &
361               microphysics_morrison, microphysics_seifert
362
363    USE buoyancy_mod,                                                          &
364        ONLY:  buoyancy
365           
366    USE chem_modules,                                                          &
367        ONLY:  call_chem_at_all_substeps, chem_gasphase_on, cs_name, do_depo
368
369    USE chem_photolysis_mod,                                                   &
370        ONLY:  photolysis_control
371
372    USE chemistry_model_mod,                                                   &
373        ONLY:  chem_boundary_conds, chem_depo, chem_integrate,                 &
374               chem_prognostic_equations, chem_species,                        &
375               nspec, nvar, spc_names
376
377    USE control_parameters,                                                    &
378        ONLY:  air_chemistry, constant_diffusion,                              &
379               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
380               humidity, intermediate_timestep_count,                          &
381               intermediate_timestep_count_max, large_scale_forcing,           &
382               large_scale_subsidence, neutral, nudging,                       &
383               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
384               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
385               timestep_scheme, tsc, use_subsidence_tendencies,                &
386               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
387               ws_scheme_sca, urban_surface, land_surface,                     &
388               time_since_reference_point, salsa
389
390    USE coriolis_mod,                                                          &
391        ONLY:  coriolis
392
393    USE cpulog,                                                                &
394        ONLY:  cpu_log, log_point, log_point_s
395
396    USE diffusion_s_mod,                                                       &
397        ONLY:  diffusion_s
398
399    USE diffusion_u_mod,                                                       &
400        ONLY:  diffusion_u
401
402    USE diffusion_v_mod,                                                       &
403        ONLY:  diffusion_v
404
405    USE diffusion_w_mod,                                                       &
406        ONLY:  diffusion_w
407
408    USE indices,                                                               &
409        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
410               nzb, nzt, wall_flags_0
411
412    USE kinds
413
414    USE lsf_nudging_mod,                                                       &
415        ONLY:  ls_advec, nudge
416
417    USE module_interface,                                                      &
418        ONLY:  module_interface_actions
419
420    USE ocean_mod,                                                             &
421        ONLY:  ocean_prognostic_equations, stokes_drift_terms, stokes_force,   &
422               wave_breaking, wave_breaking_term
423
424    USE plant_canopy_model_mod,                                                &
425        ONLY:  cthf, pcm_tendency
426
427    USE radiation_model_mod,                                                   &
428        ONLY:  radiation, radiation_tendency,                                  &
429               skip_time_do_radiation
430                         
431    USE salsa_mod,                                                             &
432        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, nbins, &
433               ncc_tot, ngast, salsa_boundary_conds, salsa_diagnostics,        &
434               salsa_driver, salsa_gas, salsa_gases_from_chem, salsa_tendency, &
435               skip_time_do_salsa
436               
437    USE salsa_util_mod,                                                        &
438        ONLY:  sums_salsa_ws_l                             
439   
440    USE statistics,                                                            &
441        ONLY:  hom
442
443    USE subsidence_mod,                                                        &
444        ONLY:  subsidence
445
446    USE surface_mod,                                                           &
447        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
448                surf_usm_v
449
450    USE turbulence_closure_mod,                                                &
451        ONLY:  tcm_prognostic_equations
452
453    USE wind_turbine_model_mod,                                                &
454        ONLY:  wtm_tendencies
455
456
457    PRIVATE
458    PUBLIC prognostic_equations_cache, prognostic_equations_vector
459
460    INTERFACE prognostic_equations_cache
461       MODULE PROCEDURE prognostic_equations_cache
462    END INTERFACE prognostic_equations_cache
463
464    INTERFACE prognostic_equations_vector
465       MODULE PROCEDURE prognostic_equations_vector
466    END INTERFACE prognostic_equations_vector
467
468
469 CONTAINS
470
471
472!------------------------------------------------------------------------------!
473! Description:
474! ------------
475!> Version with one optimized loop over all equations. It is only allowed to
476!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
477!>
478!> Here the calls of most subroutines are embedded in two DO loops over i and j,
479!> so communication between CPUs is not allowed (does not make sense) within
480!> these loops.
481!>
482!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
483!------------------------------------------------------------------------------!
484
485 SUBROUTINE prognostic_equations_cache
486
487
488    IMPLICIT NONE
489   
490    INTEGER(iwp) ::  b                   !< index for aerosol size bins (salsa)
491    INTEGER(iwp) ::  c                   !< index for chemical compounds (salsa)
492    INTEGER(iwp) ::  g                   !< index for gaseous compounds (salsa)
493    INTEGER(iwp) ::  i                   !<
494    INTEGER(iwp) ::  i_omp_start         !<
495    INTEGER(iwp) ::  j                   !<
496    INTEGER(iwp) ::  k                   !<
497!$  INTEGER(iwp) ::  omp_get_thread_num  !<
498    INTEGER(iwp) ::  tn = 0              !<
499
500    LOGICAL      ::  loop_start          !<
501    INTEGER(iwp) ::  lsp
502    INTEGER(iwp) ::  lsp_usr             !< lsp running index for chem spcs
503
504
505!
506!-- Time measurement can only be performed for the whole set of equations
507    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
508
509!
510!-- Calculation of chemical reactions. This is done outside of main loop,
511!-- since exchange of ghost points is required after this update of the
512!-- concentrations of chemical species                                   
513    IF ( air_chemistry )  THEN
514       lsp_usr = 1
515!
516!--    Chemical reactions and deposition
517       IF ( chem_gasphase_on ) THEN
518!
519!--       If required, calculate photolysis frequencies -
520!--       UNFINISHED: Why not before the intermediate timestep loop?
521          IF ( intermediate_timestep_count ==  1 )  THEN
522             CALL photolysis_control
523          ENDIF
524          DO  i = nxl, nxr
525             DO  j = nys, nyn
526
527                IF ( intermediate_timestep_count == 1 .OR.                        &
528                     call_chem_at_all_substeps )  THEN
529
530                   CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' ) 
531                   CALL chem_integrate (i,j)
532                   CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
533
534                   IF ( do_depo )  THEN
535                      CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
536                      CALL chem_depo(i,j)
537                      CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
538                   ENDIF
539                ENDIF
540             ENDDO
541          ENDDO
542       ENDIF
543!
544!--    Loop over chemical species       
545       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
546       DO  lsp = 1, nspec
547          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
548          lsp_usr = 1 
549          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
550             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
551
552                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
553                                          chem_species(lsp)%conc_pr_init ) 
554             
555             ENDIF
556             lsp_usr = lsp_usr +1
557          ENDDO
558
559
560       ENDDO
561       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
562
563    ENDIF       
564!
565!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
566!-- step. The exchange of ghost points is required after this update of the
567!-- concentrations of aerosol number and mass
568    IF ( salsa )  THEN
569       
570       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN                 
571          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
572          THEN
573             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
574             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
575             !$OMP DO
576!
577!--          Call salsa processes
578             DO  i = nxl, nxr
579                DO  j = nys, nyn
580                   CALL salsa_diagnostics( i, j )
581                   CALL salsa_driver( i, j, 3 )
582                   CALL salsa_diagnostics( i, j )
583                ENDDO
584             ENDDO
585             
586             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
587             
588             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
589!
590!--          Exchange ghost points and decycle if needed.
591             DO  b = 1, nbins
592                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
593                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
594                                           aerosol_number(b)%init )
595                DO  c = 1, ncc_tot
596                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
597                   CALL salsa_boundary_conds(                                  &
598                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
599                                           aerosol_mass((c-1)*nbins+b)%init )
600                ENDDO
601             ENDDO
602             
603             IF ( .NOT. salsa_gases_from_chem )  THEN
604                DO  g = 1, ngast
605                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
606                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
607                                              salsa_gas(g)%init )
608                ENDDO
609             ENDIF
610             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
611             
612             !$OMP END PARALLEL
613             last_salsa_time = time_since_reference_point
614             
615          ENDIF
616         
617       ENDIF
618       
619    ENDIF
620
621!
622!-- If required, calculate cloud microphysics
623    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
624         ( intermediate_timestep_count == 1  .OR.                              &
625           call_microphysics_at_all_substeps ) )                               &
626    THEN
627       !$OMP PARALLEL PRIVATE (i,j)
628       !$OMP DO
629       DO  i = nxlg, nxrg
630          DO  j = nysg, nyng
631             CALL bcm_actions( i, j )
632           ENDDO
633       ENDDO
634       !$OMP END PARALLEL
635    ENDIF
636
637!
638!-- Loop over all prognostic equations
639!-- b, c ang g added for SALSA
640    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
641
642    !$  tn = omp_get_thread_num()
643    loop_start = .TRUE.
644
645    !$OMP DO
646    DO  i = nxl, nxr
647
648!
649!--    Store the first loop index. It differs for each thread and is required
650!--    later in advec_ws
651       IF ( loop_start )  THEN
652          loop_start  = .FALSE.
653          i_omp_start = i
654       ENDIF
655
656       DO  j = nys, nyn
657!
658!--       Tendency terms for u-velocity component. Please note, in case of
659!--       non-cyclic boundary conditions the grid point i=0 is excluded from
660!--       the prognostic equations for the u-component.   
661          IF ( i >= nxlu )  THEN
662
663             tend(:,j,i) = 0.0_wp
664             IF ( timestep_scheme(1:5) == 'runge' )  THEN
665                IF ( ws_scheme_mom )  THEN
666                   CALL advec_u_ws( i, j, i_omp_start, tn )
667                ELSE
668                   CALL advec_u_pw( i, j )
669                ENDIF
670             ELSE
671                CALL advec_u_up( i, j )
672             ENDIF
673             CALL diffusion_u( i, j )
674             CALL coriolis( i, j, 1 )
675             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
676                CALL buoyancy( i, j, pt, 1 )
677             ENDIF
678
679!
680!--          Drag by plant canopy
681             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
682
683!
684!--          External pressure gradient
685             IF ( dp_external )  THEN
686                DO  k = dp_level_ind_b+1, nzt
687                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
688                ENDDO
689             ENDIF
690
691!
692!--          Nudging
693             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
694
695!
696!--          Effect of Stokes drift (in ocean mode only)
697             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
698
699!
700!--          Forces by wind turbines
701             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
702
703             CALL module_interface_actions( i, j, 'u-tendency' )
704!
705!--          Prognostic equation for u-velocity component
706             DO  k = nzb+1, nzt
707
708                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
709                                            ( tsc(2) * tend(k,j,i) +            &
710                                              tsc(3) * tu_m(k,j,i) )            &
711                                            - tsc(5) * rdf(k)                   &
712                                                     * ( u(k,j,i) - u_init(k) ) &
713                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
714                                                   BTEST( wall_flags_0(k,j,i), 1 )&
715                                                 )
716             ENDDO
717
718!
719!--          Add turbulence generated by wave breaking (in ocean mode only)
720             IF ( wave_breaking  .AND.                                         &
721               intermediate_timestep_count == intermediate_timestep_count_max )&
722             THEN
723                CALL wave_breaking_term( i, j, 1 )
724             ENDIF
725
726!
727!--          Calculate tendencies for the next Runge-Kutta step
728             IF ( timestep_scheme(1:5) == 'runge' )  THEN
729                IF ( intermediate_timestep_count == 1 )  THEN
730                   DO  k = nzb+1, nzt
731                      tu_m(k,j,i) = tend(k,j,i)
732                   ENDDO
733                ELSEIF ( intermediate_timestep_count < &
734                         intermediate_timestep_count_max )  THEN
735                   DO  k = nzb+1, nzt
736                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
737                                     + 5.3125_wp * tu_m(k,j,i)
738                   ENDDO
739                ENDIF
740             ENDIF
741
742          ENDIF
743!
744!--       Tendency terms for v-velocity component. Please note, in case of
745!--       non-cyclic boundary conditions the grid point j=0 is excluded from
746!--       the prognostic equations for the v-component. !--       
747          IF ( j >= nysv )  THEN
748
749             tend(:,j,i) = 0.0_wp
750             IF ( timestep_scheme(1:5) == 'runge' )  THEN
751                IF ( ws_scheme_mom )  THEN
752                    CALL advec_v_ws( i, j, i_omp_start, tn )
753                ELSE
754                    CALL advec_v_pw( i, j )
755                ENDIF
756             ELSE
757                CALL advec_v_up( i, j )
758             ENDIF
759             CALL diffusion_v( i, j )
760             CALL coriolis( i, j, 2 )
761
762!
763!--          Drag by plant canopy
764             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
765
766!
767!--          External pressure gradient
768             IF ( dp_external )  THEN
769                DO  k = dp_level_ind_b+1, nzt
770                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
771                ENDDO
772             ENDIF
773
774!
775!--          Nudging
776             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
777
778!
779!--          Effect of Stokes drift (in ocean mode only)
780             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
781
782!
783!--          Forces by wind turbines
784             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
785
786             CALL module_interface_actions( i, j, 'v-tendency' )
787!
788!--          Prognostic equation for v-velocity component
789             DO  k = nzb+1, nzt
790                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
791                                            ( tsc(2) * tend(k,j,i) +           &
792                                              tsc(3) * tv_m(k,j,i) )           &
793                                            - tsc(5) * rdf(k)                  &
794                                                     * ( v(k,j,i) - v_init(k) )&
795                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
796                                                   BTEST( wall_flags_0(k,j,i), 2 )&
797                                                 )
798             ENDDO
799
800!
801!--          Add turbulence generated by wave breaking (in ocean mode only)
802             IF ( wave_breaking  .AND.                                         &
803               intermediate_timestep_count == intermediate_timestep_count_max )&
804             THEN
805                CALL wave_breaking_term( i, j, 2 )
806             ENDIF
807
808!
809!--          Calculate tendencies for the next Runge-Kutta step
810             IF ( timestep_scheme(1:5) == 'runge' )  THEN
811                IF ( intermediate_timestep_count == 1 )  THEN
812                   DO  k = nzb+1, nzt
813                      tv_m(k,j,i) = tend(k,j,i)
814                   ENDDO
815                ELSEIF ( intermediate_timestep_count < &
816                         intermediate_timestep_count_max )  THEN
817                   DO  k = nzb+1, nzt
818                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
819                                     + 5.3125_wp * tv_m(k,j,i)
820                   ENDDO
821                ENDIF
822             ENDIF
823
824          ENDIF
825
826!
827!--       Tendency terms for w-velocity component
828          tend(:,j,i) = 0.0_wp
829          IF ( timestep_scheme(1:5) == 'runge' )  THEN
830             IF ( ws_scheme_mom )  THEN
831                CALL advec_w_ws( i, j, i_omp_start, tn )
832             ELSE
833                CALL advec_w_pw( i, j )
834             END IF
835          ELSE
836             CALL advec_w_up( i, j )
837          ENDIF
838          CALL diffusion_w( i, j )
839          CALL coriolis( i, j, 3 )
840
841          IF ( .NOT. neutral )  THEN
842             IF ( ocean_mode )  THEN
843                CALL buoyancy( i, j, rho_ocean, 3 )
844             ELSE
845                IF ( .NOT. humidity )  THEN
846                   CALL buoyancy( i, j, pt, 3 )
847                ELSE
848                   CALL buoyancy( i, j, vpt, 3 )
849                ENDIF
850             ENDIF
851          ENDIF
852
853!
854!--       Drag by plant canopy
855          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
856
857!
858!--       Effect of Stokes drift (in ocean mode only)
859          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
860
861!
862!--       Forces by wind turbines
863          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
864
865          CALL module_interface_actions( i, j, 'w-tendency' )
866!
867!--       Prognostic equation for w-velocity component
868          DO  k = nzb+1, nzt-1
869             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
870                                             ( tsc(2) * tend(k,j,i) +          &
871                                               tsc(3) * tw_m(k,j,i) )          &
872                                             - tsc(5) * rdf(k) * w(k,j,i)      &
873                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
874                                                BTEST( wall_flags_0(k,j,i), 3 )&
875                                              )
876          ENDDO
877
878!
879!--       Calculate tendencies for the next Runge-Kutta step
880          IF ( timestep_scheme(1:5) == 'runge' )  THEN
881             IF ( intermediate_timestep_count == 1 )  THEN
882                DO  k = nzb+1, nzt-1
883                   tw_m(k,j,i) = tend(k,j,i)
884                ENDDO
885             ELSEIF ( intermediate_timestep_count < &
886                      intermediate_timestep_count_max )  THEN
887                DO  k = nzb+1, nzt-1
888                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
889                                  + 5.3125_wp * tw_m(k,j,i)
890                ENDDO
891             ENDIF
892          ENDIF
893
894!
895!--       If required, compute prognostic equation for potential temperature
896          IF ( .NOT. neutral )  THEN
897!
898!--          Tendency terms for potential temperature
899             tend(:,j,i) = 0.0_wp
900             IF ( timestep_scheme(1:5) == 'runge' )  THEN
901                   IF ( ws_scheme_sca )  THEN
902                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
903                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
904                   ELSE
905                      CALL advec_s_pw( i, j, pt )
906                   ENDIF
907             ELSE
908                CALL advec_s_up( i, j, pt )
909             ENDIF
910             CALL diffusion_s( i, j, pt,                                       &
911                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
912                               surf_def_h(2)%shf,                              &
913                               surf_lsm_h%shf,    surf_usm_h%shf,              &
914                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
915                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
916                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
917                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
918                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
919                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
920
921!
922!--          Consideration of heat sources within the plant canopy
923             IF ( plant_canopy  .AND.                                          &
924                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
925                CALL pcm_tendency( i, j, 4 )
926             ENDIF
927
928!
929!--          Large scale advection
930             IF ( large_scale_forcing )  THEN
931                CALL ls_advec( i, j, simulated_time, 'pt' )
932             ENDIF
933
934!
935!--          Nudging
936             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
937
938!
939!--          If required, compute effect of large-scale subsidence/ascent
940             IF ( large_scale_subsidence  .AND.                                &
941                  .NOT. use_subsidence_tendencies )  THEN
942                CALL subsidence( i, j, tend, pt, pt_init, 2 )
943             ENDIF
944
945!
946!--          If required, add tendency due to radiative heating/cooling
947             IF ( radiation  .AND.                                             &
948                  simulated_time > skip_time_do_radiation )  THEN
949                CALL radiation_tendency ( i, j, tend )
950             ENDIF
951
952
953             CALL module_interface_actions( i, j, 'pt-tendency' )
954!
955!--          Prognostic equation for potential temperature
956             DO  k = nzb+1, nzt
957                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
958                                                  ( tsc(2) * tend(k,j,i) +     &
959                                                    tsc(3) * tpt_m(k,j,i) )    &
960                                                  - tsc(5)                     &
961                                                  * ( pt(k,j,i) - pt_init(k) ) &
962                                                  * ( rdf_sc(k) + ptdf_x(i)    &
963                                                                + ptdf_y(j) )  &
964                                          )                                    &
965                                       * MERGE( 1.0_wp, 0.0_wp,                &
966                                                BTEST( wall_flags_0(k,j,i), 0 )&
967                                              )
968             ENDDO
969
970!
971!--          Calculate tendencies for the next Runge-Kutta step
972             IF ( timestep_scheme(1:5) == 'runge' )  THEN
973                IF ( intermediate_timestep_count == 1 )  THEN
974                   DO  k = nzb+1, nzt
975                      tpt_m(k,j,i) = tend(k,j,i)
976                   ENDDO
977                ELSEIF ( intermediate_timestep_count < &
978                         intermediate_timestep_count_max )  THEN
979                   DO  k = nzb+1, nzt
980                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
981                                        5.3125_wp * tpt_m(k,j,i)
982                   ENDDO
983                ENDIF
984             ENDIF
985
986          ENDIF
987
988!
989!--       If required, compute prognostic equation for total water content
990          IF ( humidity )  THEN
991
992!
993!--          Tendency-terms for total water content / scalar
994             tend(:,j,i) = 0.0_wp
995             IF ( timestep_scheme(1:5) == 'runge' ) &
996             THEN
997                IF ( ws_scheme_sca )  THEN
998                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
999                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1000                ELSE
1001                   CALL advec_s_pw( i, j, q )
1002                ENDIF
1003             ELSE
1004                CALL advec_s_up( i, j, q )
1005             ENDIF
1006             CALL diffusion_s( i, j, q,                                        &
1007                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
1008                               surf_def_h(2)%qsws,                             &
1009                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
1010                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
1011                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
1012                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
1013                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
1014                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
1015                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1016
1017!
1018!--          Sink or source of humidity due to canopy elements
1019             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
1020
1021!
1022!--          Large scale advection
1023             IF ( large_scale_forcing )  THEN
1024                CALL ls_advec( i, j, simulated_time, 'q' )
1025             ENDIF
1026
1027!
1028!--          Nudging
1029             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
1030
1031!
1032!--          If required compute influence of large-scale subsidence/ascent
1033             IF ( large_scale_subsidence  .AND.                                &
1034                  .NOT. use_subsidence_tendencies )  THEN
1035                CALL subsidence( i, j, tend, q, q_init, 3 )
1036             ENDIF
1037
1038             CALL module_interface_actions( i, j, 'q-tendency' )
1039
1040!
1041!--          Prognostic equation for total water content / scalar
1042             DO  k = nzb+1, nzt
1043                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
1044                                                ( tsc(2) * tend(k,j,i) +       &
1045                                                  tsc(3) * tq_m(k,j,i) )       &
1046                                                - tsc(5) * rdf_sc(k) *         &
1047                                                      ( q(k,j,i) - q_init(k) ) &
1048                                        )                                      &
1049                                       * MERGE( 1.0_wp, 0.0_wp,                &
1050                                                BTEST( wall_flags_0(k,j,i), 0 )&
1051                                              )               
1052                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1053             ENDDO
1054
1055!
1056!--          Calculate tendencies for the next Runge-Kutta step
1057             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1058                IF ( intermediate_timestep_count == 1 )  THEN
1059                   DO  k = nzb+1, nzt
1060                      tq_m(k,j,i) = tend(k,j,i)
1061                   ENDDO
1062                ELSEIF ( intermediate_timestep_count < &
1063                         intermediate_timestep_count_max )  THEN
1064                   DO  k = nzb+1, nzt
1065                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1066                                       5.3125_wp * tq_m(k,j,i)
1067                   ENDDO
1068                ENDIF
1069             ENDIF
1070
1071!
1072!--          If required, calculate prognostic equations for cloud water content
1073!--          and cloud drop concentration
1074             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1075!
1076!--             Calculate prognostic equation for cloud water content
1077                tend(:,j,i) = 0.0_wp
1078                IF ( timestep_scheme(1:5) == 'runge' ) &
1079                THEN
1080                   IF ( ws_scheme_sca )  THEN
1081                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
1082                                       diss_s_qc, flux_l_qc, diss_l_qc, &
1083                                       i_omp_start, tn )
1084                   ELSE
1085                      CALL advec_s_pw( i, j, qc )
1086                   ENDIF
1087                ELSE
1088                   CALL advec_s_up( i, j, qc )
1089                ENDIF
1090                CALL diffusion_s( i, j, qc,                                   &
1091                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
1092                                  surf_def_h(2)%qcsws,                        &
1093                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      & 
1094                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
1095                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
1096                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
1097                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
1098                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
1099                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
1100
1101!
1102!--             Prognostic equation for cloud water content
1103                DO  k = nzb+1, nzt
1104                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
1105                                                      ( tsc(2) * tend(k,j,i) + &
1106                                                        tsc(3) * tqc_m(k,j,i) )&
1107                                                      - tsc(5) * rdf_sc(k)     &
1108                                                               * qc(k,j,i)     &
1109                                             )                                 &
1110                                       * MERGE( 1.0_wp, 0.0_wp,                &
1111                                                BTEST( wall_flags_0(k,j,i), 0 )&
1112                                              ) 
1113                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
1114                ENDDO
1115!
1116!--             Calculate tendencies for the next Runge-Kutta step
1117                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1118                   IF ( intermediate_timestep_count == 1 )  THEN
1119                      DO  k = nzb+1, nzt
1120                         tqc_m(k,j,i) = tend(k,j,i)
1121                      ENDDO
1122                   ELSEIF ( intermediate_timestep_count < &
1123                            intermediate_timestep_count_max )  THEN
1124                      DO  k = nzb+1, nzt
1125                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1126                                           5.3125_wp * tqc_m(k,j,i)
1127                      ENDDO
1128                   ENDIF
1129                ENDIF
1130
1131!
1132!--             Calculate prognostic equation for cloud drop concentration.
1133                tend(:,j,i) = 0.0_wp
1134                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1135                   IF ( ws_scheme_sca )  THEN
1136                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
1137                                    diss_s_nc, flux_l_nc, diss_l_nc, &
1138                                    i_omp_start, tn )
1139                   ELSE
1140                      CALL advec_s_pw( i, j, nc )
1141                   ENDIF
1142                ELSE
1143                   CALL advec_s_up( i, j, nc )
1144                ENDIF
1145                CALL diffusion_s( i, j, nc,                                    &
1146                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
1147                                  surf_def_h(2)%ncsws,                         &
1148                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
1149                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
1150                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
1151                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
1152                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
1153                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
1154                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
1155
1156!
1157!--             Prognostic equation for cloud drop concentration
1158                DO  k = nzb+1, nzt
1159                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
1160                                                      ( tsc(2) * tend(k,j,i) + &
1161                                                        tsc(3) * tnc_m(k,j,i) )&
1162                                                      - tsc(5) * rdf_sc(k)     &
1163                                                               * nc(k,j,i)     &
1164                                             )                                 &
1165                                       * MERGE( 1.0_wp, 0.0_wp,                &
1166                                                BTEST( wall_flags_0(k,j,i), 0 )&
1167                                              )
1168                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
1169                ENDDO
1170!
1171!--             Calculate tendencies for the next Runge-Kutta step
1172                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1173                   IF ( intermediate_timestep_count == 1 )  THEN
1174                      DO  k = nzb+1, nzt
1175                         tnc_m(k,j,i) = tend(k,j,i)
1176                      ENDDO
1177                   ELSEIF ( intermediate_timestep_count < &
1178                            intermediate_timestep_count_max )  THEN
1179                      DO  k = nzb+1, nzt
1180                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1181                                           5.3125_wp * tnc_m(k,j,i)
1182                      ENDDO
1183                   ENDIF
1184                ENDIF
1185
1186             ENDIF
1187!
1188!--          If required, calculate prognostic equations for rain water content
1189!--          and rain drop concentration
1190             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1191!
1192!--             Calculate prognostic equation for rain water content
1193                tend(:,j,i) = 0.0_wp
1194                IF ( timestep_scheme(1:5) == 'runge' ) &
1195                THEN
1196                   IF ( ws_scheme_sca )  THEN
1197                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
1198                                       diss_s_qr, flux_l_qr, diss_l_qr, &
1199                                       i_omp_start, tn )
1200                   ELSE
1201                      CALL advec_s_pw( i, j, qr )
1202                   ENDIF
1203                ELSE
1204                   CALL advec_s_up( i, j, qr )
1205                ENDIF
1206                CALL diffusion_s( i, j, qr,                                   &
1207                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
1208                                  surf_def_h(2)%qrsws,                        &
1209                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      & 
1210                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
1211                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
1212                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
1213                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
1214                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
1215                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1216
1217!
1218!--             Prognostic equation for rain water content
1219                DO  k = nzb+1, nzt
1220                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
1221                                                      ( tsc(2) * tend(k,j,i) + &
1222                                                        tsc(3) * tqr_m(k,j,i) )&
1223                                                      - tsc(5) * rdf_sc(k)     &
1224                                                               * qr(k,j,i)     &
1225                                             )                                 &
1226                                       * MERGE( 1.0_wp, 0.0_wp,                &
1227                                                BTEST( wall_flags_0(k,j,i), 0 )&
1228                                              ) 
1229                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
1230                ENDDO
1231!
1232!--             Calculate tendencies for the next Runge-Kutta step
1233                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1234                   IF ( intermediate_timestep_count == 1 )  THEN
1235                      DO  k = nzb+1, nzt
1236                         tqr_m(k,j,i) = tend(k,j,i)
1237                      ENDDO
1238                   ELSEIF ( intermediate_timestep_count < &
1239                            intermediate_timestep_count_max )  THEN
1240                      DO  k = nzb+1, nzt
1241                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1242                                           5.3125_wp * tqr_m(k,j,i)
1243                      ENDDO
1244                   ENDIF
1245                ENDIF
1246
1247!
1248!--             Calculate prognostic equation for rain drop concentration.
1249                tend(:,j,i) = 0.0_wp
1250                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1251                   IF ( ws_scheme_sca )  THEN
1252                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
1253                                    diss_s_nr, flux_l_nr, diss_l_nr, &
1254                                    i_omp_start, tn )
1255                   ELSE
1256                      CALL advec_s_pw( i, j, nr )
1257                   ENDIF
1258                ELSE
1259                   CALL advec_s_up( i, j, nr )
1260                ENDIF
1261                CALL diffusion_s( i, j, nr,                                    &
1262                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
1263                                  surf_def_h(2)%nrsws,                         &
1264                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
1265                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
1266                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
1267                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
1268                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
1269                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
1270                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1271
1272!
1273!--             Prognostic equation for rain drop concentration
1274                DO  k = nzb+1, nzt
1275                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
1276                                                      ( tsc(2) * tend(k,j,i) + &
1277                                                        tsc(3) * tnr_m(k,j,i) )&
1278                                                      - tsc(5) * rdf_sc(k)     &
1279                                                               * nr(k,j,i)     &
1280                                             )                                 &
1281                                       * MERGE( 1.0_wp, 0.0_wp,                &
1282                                                BTEST( wall_flags_0(k,j,i), 0 )&
1283                                              )
1284                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1285                ENDDO
1286!
1287!--             Calculate tendencies for the next Runge-Kutta step
1288                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1289                   IF ( intermediate_timestep_count == 1 )  THEN
1290                      DO  k = nzb+1, nzt
1291                         tnr_m(k,j,i) = tend(k,j,i)
1292                      ENDDO
1293                   ELSEIF ( intermediate_timestep_count < &
1294                            intermediate_timestep_count_max )  THEN
1295                      DO  k = nzb+1, nzt
1296                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1297                                           5.3125_wp * tnr_m(k,j,i)
1298                      ENDDO
1299                   ENDIF
1300                ENDIF
1301
1302             ENDIF
1303
1304          ENDIF
1305
1306!
1307!--       If required, compute prognostic equation for scalar
1308          IF ( passive_scalar )  THEN
1309!
1310!--          Tendency-terms for total water content / scalar
1311             tend(:,j,i) = 0.0_wp
1312             IF ( timestep_scheme(1:5) == 'runge' ) &
1313             THEN
1314                IF ( ws_scheme_sca )  THEN
1315                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1316                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1317                ELSE
1318                   CALL advec_s_pw( i, j, s )
1319                ENDIF
1320             ELSE
1321                CALL advec_s_up( i, j, s )
1322             ENDIF
1323             CALL diffusion_s( i, j, s,                                        &
1324                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1325                               surf_def_h(2)%ssws,                             &
1326                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1327                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1328                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1329                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1330                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1331                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1332                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1333
1334!
1335!--          Sink or source of scalar concentration due to canopy elements
1336             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1337
1338!
1339!--          Large scale advection, still need to be extended for scalars
1340!              IF ( large_scale_forcing )  THEN
1341!                 CALL ls_advec( i, j, simulated_time, 's' )
1342!              ENDIF
1343
1344!
1345!--          Nudging, still need to be extended for scalars
1346!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1347
1348!
1349!--          If required compute influence of large-scale subsidence/ascent.
1350!--          Note, the last argument is of no meaning in this case, as it is
1351!--          only used in conjunction with large_scale_forcing, which is to
1352!--          date not implemented for scalars.
1353             IF ( large_scale_subsidence  .AND.                                &
1354                  .NOT. use_subsidence_tendencies  .AND.                       &
1355                  .NOT. large_scale_forcing )  THEN
1356                CALL subsidence( i, j, tend, s, s_init, 3 )
1357             ENDIF
1358
1359             CALL module_interface_actions( i, j, 's-tendency' )
1360
1361!
1362!--          Prognostic equation for scalar
1363             DO  k = nzb+1, nzt
1364                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1365                                            ( tsc(2) * tend(k,j,i) +           &
1366                                              tsc(3) * ts_m(k,j,i) )           &
1367                                            - tsc(5) * rdf_sc(k)               &
1368                                                     * ( s(k,j,i) - s_init(k) )&
1369                                        )                                      &
1370                                       * MERGE( 1.0_wp, 0.0_wp,                &
1371                                                BTEST( wall_flags_0(k,j,i), 0 )&
1372                                              )
1373                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1374             ENDDO
1375
1376!
1377!--          Calculate tendencies for the next Runge-Kutta step
1378             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1379                IF ( intermediate_timestep_count == 1 )  THEN
1380                   DO  k = nzb+1, nzt
1381                      ts_m(k,j,i) = tend(k,j,i)
1382                   ENDDO
1383                ELSEIF ( intermediate_timestep_count < &
1384                         intermediate_timestep_count_max )  THEN
1385                   DO  k = nzb+1, nzt
1386                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1387                                       5.3125_wp * ts_m(k,j,i)
1388                   ENDDO
1389                ENDIF
1390             ENDIF
1391
1392          ENDIF
1393!
1394!--       Calculate prognostic equations for turbulence closure
1395          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1396
1397!
1398!--       Calculate prognostic equation for chemical quantites
1399          IF ( air_chemistry )  THEN
1400             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
1401!
1402!--          Loop over chemical species
1403             DO  lsp = 1, nvar                         
1404                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
1405                                     chem_species(lsp)%conc,                   &
1406                                     chem_species(lsp)%tconc_m,                &
1407                                     chem_species(lsp)%conc_pr_init,           &
1408                                     i, j, i_omp_start, tn, lsp,               &
1409                                     chem_species(lsp)%flux_s_cs,              &
1410                                     chem_species(lsp)%diss_s_cs,              &
1411                                     chem_species(lsp)%flux_l_cs,              &
1412                                     chem_species(lsp)%diss_l_cs )       
1413             ENDDO
1414         
1415          ENDIF   ! Chemical equations
1416!
1417!--       Calculate prognostic equations for the ocean
1418          IF ( ocean_mode )  THEN
1419             CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
1420          ENDIF
1421       
1422          IF ( salsa )  THEN 
1423!
1424!--          Loop over aerosol size bins: number and mass bins
1425             IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1426             
1427                DO  b = 1, nbins               
1428                   sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
1429                   CALL salsa_tendency( 'aerosol_number',                      &
1430                                         aerosol_number(b)%conc_p,             &
1431                                         aerosol_number(b)%conc,               &
1432                                         aerosol_number(b)%tconc_m,            &
1433                                         i, j, i_omp_start, tn, b, b,          &
1434                                         aerosol_number(b)%flux_s,             &
1435                                         aerosol_number(b)%diss_s,             &
1436                                         aerosol_number(b)%flux_l,             &
1437                                         aerosol_number(b)%diss_l,             &
1438                                         aerosol_number(b)%init )
1439                   aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
1440                   DO  c = 1, ncc_tot
1441                      sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
1442                      CALL salsa_tendency( 'aerosol_mass',                     &
1443                                            aerosol_mass((c-1)*nbins+b)%conc_p,&
1444                                            aerosol_mass((c-1)*nbins+b)%conc,  &
1445                                            aerosol_mass((c-1)*nbins+b)%tconc_m,&
1446                                            i, j, i_omp_start, tn, b, c,       &
1447                                            aerosol_mass((c-1)*nbins+b)%flux_s,&
1448                                            aerosol_mass((c-1)*nbins+b)%diss_s,&
1449                                            aerosol_mass((c-1)*nbins+b)%flux_l,&
1450                                            aerosol_mass((c-1)*nbins+b)%diss_l,&
1451                                            aerosol_mass((c-1)*nbins+b)%init )
1452                      aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
1453                   ENDDO
1454                ENDDO
1455                IF ( .NOT. salsa_gases_from_chem )  THEN
1456                   DO  g = 1, ngast
1457                      sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
1458                      CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,   &
1459                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
1460                                      i, j, i_omp_start, tn, g, g,             &
1461                                      salsa_gas(g)%flux_s, salsa_gas(g)%diss_s,&
1462                                      salsa_gas(g)%flux_l, salsa_gas(g)%diss_l,&
1463                                      salsa_gas(g)%init )
1464                      salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
1465                   ENDDO
1466                ENDIF
1467             
1468             ENDIF
1469             
1470          ENDIF       
1471
1472       ENDDO  ! loop over j
1473    ENDDO  ! loop over i
1474    !$OMP END PARALLEL
1475
1476
1477    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1478
1479
1480 END SUBROUTINE prognostic_equations_cache
1481
1482
1483!------------------------------------------------------------------------------!
1484! Description:
1485! ------------
1486!> Version for vector machines
1487!------------------------------------------------------------------------------!
1488
1489 SUBROUTINE prognostic_equations_vector
1490
1491
1492    IMPLICIT NONE
1493
1494    INTEGER(iwp) ::  b     !< index for aerosol size bins (salsa)
1495    INTEGER(iwp) ::  c     !< index for chemical compounds (salsa)
1496    INTEGER(iwp) ::  g     !< index for gaseous compounds (salsa)
1497    INTEGER(iwp) ::  i     !<
1498    INTEGER(iwp) ::  j     !<
1499    INTEGER(iwp) ::  k     !<
1500    INTEGER(iwp) ::  lsp   !< running index for chemical species
1501
1502    REAL(wp)     ::  sbt  !<
1503
1504!
1505!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1506!-- step. The exchange of ghost points is required after this update of the
1507!-- concentrations of aerosol number and mass
1508    IF ( salsa )  THEN
1509       
1510       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1511       
1512          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
1513          THEN
1514         
1515             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1516             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
1517             !$OMP DO
1518!
1519!--          Call salsa processes
1520             DO  i = nxl, nxr
1521                DO  j = nys, nyn
1522                   CALL salsa_diagnostics( i, j )
1523                   CALL salsa_driver( i, j, 3 )
1524                   CALL salsa_diagnostics( i, j )
1525                ENDDO
1526             ENDDO
1527             
1528             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1529             
1530             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1531!
1532!--          Exchange ghost points and decycle if needed.
1533             DO  b = 1, nbins
1534                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
1535                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
1536                                           aerosol_number(b)%init )
1537                DO  c = 1, ncc_tot
1538                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
1539                   CALL salsa_boundary_conds(                                  &
1540                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
1541                                           aerosol_mass((c-1)*nbins+b)%init )
1542                ENDDO
1543             ENDDO
1544             
1545             IF ( .NOT. salsa_gases_from_chem )  THEN
1546                DO  g = 1, ngast
1547                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
1548                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
1549                                              salsa_gas(g)%init )
1550                ENDDO
1551             ENDIF
1552             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1553             
1554             !$OMP END PARALLEL
1555             last_salsa_time = time_since_reference_point
1556             
1557          ENDIF
1558         
1559       ENDIF
1560       
1561    ENDIF
1562
1563!
1564!-- If required, calculate cloud microphysical impacts
1565    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
1566         ( intermediate_timestep_count == 1  .OR.                              &
1567           call_microphysics_at_all_substeps )                                 &
1568       )  THEN
1569       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1570       CALL bcm_actions
1571       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1572    ENDIF
1573
1574!
1575!-- u-velocity component
1576    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1577
1578    !$ACC KERNELS PRESENT(tend)
1579    tend = 0.0_wp
1580    !$ACC END KERNELS
1581    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1582       IF ( ws_scheme_mom )  THEN
1583          CALL advec_u_ws
1584       ELSE
1585          CALL advec_u_pw
1586       ENDIF
1587    ELSE
1588       CALL advec_u_up
1589    ENDIF
1590    CALL diffusion_u
1591    CALL coriolis( 1 )
1592    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1593       CALL buoyancy( pt, 1 )
1594    ENDIF
1595
1596!
1597!-- Drag by plant canopy
1598    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1599
1600!
1601!-- External pressure gradient
1602    IF ( dp_external )  THEN
1603       DO  i = nxlu, nxr
1604          DO  j = nys, nyn
1605             DO  k = dp_level_ind_b+1, nzt
1606                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1607             ENDDO
1608          ENDDO
1609       ENDDO
1610    ENDIF
1611
1612!
1613!-- Nudging
1614    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1615
1616!
1617!-- Effect of Stokes drift (in ocean mode only)
1618    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1619
1620!
1621!-- Forces by wind turbines
1622    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1623
1624    CALL module_interface_actions( 'u-tendency' )
1625
1626!
1627!-- Prognostic equation for u-velocity component
1628    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1629    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1630    !$ACC PRESENT(tsc(2:5)) &
1631    !$ACC PRESENT(u_p)
1632    DO  i = nxlu, nxr
1633       DO  j = nys, nyn
1634          DO  k = nzb+1, nzt
1635             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1636                                                 tsc(3) * tu_m(k,j,i) )          &
1637                                               - tsc(5) * rdf(k) *               &
1638                                                        ( u(k,j,i) - u_init(k) ) &
1639                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1640                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1641                                              )
1642          ENDDO
1643       ENDDO
1644    ENDDO
1645
1646!
1647!-- Add turbulence generated by wave breaking (in ocean mode only)
1648    IF ( wave_breaking  .AND.                                                  &
1649         intermediate_timestep_count == intermediate_timestep_count_max )      &
1650    THEN
1651       CALL wave_breaking_term( 1 )
1652    ENDIF
1653
1654!
1655!-- Calculate tendencies for the next Runge-Kutta step
1656    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1657       IF ( intermediate_timestep_count == 1 )  THEN
1658          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1659          !$ACC PRESENT(tend, tu_m)
1660          DO  i = nxlu, nxr
1661             DO  j = nys, nyn
1662                DO  k = nzb+1, nzt
1663                   tu_m(k,j,i) = tend(k,j,i)
1664                ENDDO
1665             ENDDO
1666          ENDDO
1667       ELSEIF ( intermediate_timestep_count < &
1668                intermediate_timestep_count_max )  THEN
1669          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1670          !$ACC PRESENT(tend, tu_m)
1671          DO  i = nxlu, nxr
1672             DO  j = nys, nyn
1673                DO  k = nzb+1, nzt
1674                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1675                                   + 5.3125_wp * tu_m(k,j,i)
1676                ENDDO
1677             ENDDO
1678          ENDDO
1679       ENDIF
1680    ENDIF
1681
1682    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1683
1684!
1685!-- v-velocity component
1686    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1687
1688    !$ACC KERNELS PRESENT(tend)
1689    tend = 0.0_wp
1690    !$ACC END KERNELS
1691    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1692       IF ( ws_scheme_mom )  THEN
1693          CALL advec_v_ws
1694       ELSE
1695          CALL advec_v_pw
1696       END IF
1697    ELSE
1698       CALL advec_v_up
1699    ENDIF
1700    CALL diffusion_v
1701    CALL coriolis( 2 )
1702
1703!
1704!-- Drag by plant canopy
1705    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1706
1707!
1708!-- External pressure gradient
1709    IF ( dp_external )  THEN
1710       DO  i = nxl, nxr
1711          DO  j = nysv, nyn
1712             DO  k = dp_level_ind_b+1, nzt
1713                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1714             ENDDO
1715          ENDDO
1716       ENDDO
1717    ENDIF
1718
1719!
1720!-- Nudging
1721    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1722
1723!
1724!-- Effect of Stokes drift (in ocean mode only)
1725    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1726
1727!
1728!-- Forces by wind turbines
1729    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1730
1731    CALL module_interface_actions( 'v-tendency' )
1732
1733!
1734!-- Prognostic equation for v-velocity component
1735    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1736    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1737    !$ACC PRESENT(tsc(2:5)) &
1738    !$ACC PRESENT(v_p)
1739    DO  i = nxl, nxr
1740       DO  j = nysv, nyn
1741          DO  k = nzb+1, nzt
1742             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1743                                                 tsc(3) * tv_m(k,j,i) )        &
1744                                               - tsc(5) * rdf(k) *             &
1745                                                      ( v(k,j,i) - v_init(k) ) &
1746                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1747                                                BTEST( wall_flags_0(k,j,i), 2 )&
1748                                              )
1749          ENDDO
1750       ENDDO
1751    ENDDO
1752
1753!
1754!-- Add turbulence generated by wave breaking (in ocean mode only)
1755    IF ( wave_breaking  .AND.                                                  &
1756         intermediate_timestep_count == intermediate_timestep_count_max )      &
1757    THEN
1758       CALL wave_breaking_term( 2 )
1759    ENDIF
1760
1761!
1762!-- Calculate tendencies for the next Runge-Kutta step
1763    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1764       IF ( intermediate_timestep_count == 1 )  THEN
1765          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1766          !$ACC PRESENT(tend, tv_m)
1767          DO  i = nxl, nxr
1768             DO  j = nysv, nyn
1769                DO  k = nzb+1, nzt
1770                   tv_m(k,j,i) = tend(k,j,i)
1771                ENDDO
1772             ENDDO
1773          ENDDO
1774       ELSEIF ( intermediate_timestep_count < &
1775                intermediate_timestep_count_max )  THEN
1776          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1777          !$ACC PRESENT(tend, tv_m)
1778          DO  i = nxl, nxr
1779             DO  j = nysv, nyn
1780                DO  k = nzb+1, nzt
1781                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1782                                  + 5.3125_wp * tv_m(k,j,i)
1783                ENDDO
1784             ENDDO
1785          ENDDO
1786       ENDIF
1787    ENDIF
1788
1789    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1790
1791!
1792!-- w-velocity component
1793    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1794
1795    !$ACC KERNELS PRESENT(tend)
1796    tend = 0.0_wp
1797    !$ACC END KERNELS
1798    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1799       IF ( ws_scheme_mom )  THEN
1800          CALL advec_w_ws
1801       ELSE
1802          CALL advec_w_pw
1803       ENDIF
1804    ELSE
1805       CALL advec_w_up
1806    ENDIF
1807    CALL diffusion_w
1808    CALL coriolis( 3 )
1809
1810    IF ( .NOT. neutral )  THEN
1811       IF ( ocean_mode )  THEN
1812          CALL buoyancy( rho_ocean, 3 )
1813       ELSE
1814          IF ( .NOT. humidity )  THEN
1815             CALL buoyancy( pt, 3 )
1816          ELSE
1817             CALL buoyancy( vpt, 3 )
1818          ENDIF
1819       ENDIF
1820    ENDIF
1821
1822!
1823!-- Drag by plant canopy
1824    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1825
1826!
1827!-- Effect of Stokes drift (in ocean mode only)
1828    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1829
1830!
1831!-- Forces by wind turbines
1832    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1833
1834    CALL module_interface_actions( 'w-tendency' )
1835
1836!
1837!-- Prognostic equation for w-velocity component
1838    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1839    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1840    !$ACC PRESENT(tsc(2:5)) &
1841    !$ACC PRESENT(w_p)
1842    DO  i = nxl, nxr
1843       DO  j = nys, nyn
1844          DO  k = nzb+1, nzt-1
1845             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1846                                                 tsc(3) * tw_m(k,j,i) )        &
1847                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1848                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1849                                                BTEST( wall_flags_0(k,j,i), 3 )&
1850                                              )
1851          ENDDO
1852       ENDDO
1853    ENDDO
1854
1855!
1856!-- Calculate tendencies for the next Runge-Kutta step
1857    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1858       IF ( intermediate_timestep_count == 1 )  THEN
1859          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1860          !$ACC PRESENT(tend, tw_m)
1861          DO  i = nxl, nxr
1862             DO  j = nys, nyn
1863                DO  k = nzb+1, nzt-1
1864                   tw_m(k,j,i) = tend(k,j,i)
1865                ENDDO
1866             ENDDO
1867          ENDDO
1868       ELSEIF ( intermediate_timestep_count < &
1869                intermediate_timestep_count_max )  THEN
1870          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1871          !$ACC PRESENT(tend, tw_m)
1872          DO  i = nxl, nxr
1873             DO  j = nys, nyn
1874                DO  k = nzb+1, nzt-1
1875                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1876                                  + 5.3125_wp * tw_m(k,j,i)
1877                ENDDO
1878             ENDDO
1879          ENDDO
1880       ENDIF
1881    ENDIF
1882
1883    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1884
1885
1886!
1887!-- If required, compute prognostic equation for potential temperature
1888    IF ( .NOT. neutral )  THEN
1889
1890       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1891
1892!
1893!--    pt-tendency terms with communication
1894       sbt = tsc(2)
1895       IF ( scalar_advec == 'bc-scheme' )  THEN
1896
1897          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1898!
1899!--          Bott-Chlond scheme always uses Euler time step. Thus:
1900             sbt = 1.0_wp
1901          ENDIF
1902          tend = 0.0_wp
1903          CALL advec_s_bc( pt, 'pt' )
1904
1905       ENDIF
1906
1907!
1908!--    pt-tendency terms with no communication
1909       IF ( scalar_advec /= 'bc-scheme' )  THEN
1910          !$ACC KERNELS PRESENT(tend)
1911          tend = 0.0_wp
1912          !$ACC END KERNELS
1913          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1914             IF ( ws_scheme_sca )  THEN
1915                CALL advec_s_ws( pt, 'pt' )
1916             ELSE
1917                CALL advec_s_pw( pt )
1918             ENDIF
1919          ELSE
1920             CALL advec_s_up( pt )
1921          ENDIF
1922       ENDIF
1923
1924       CALL diffusion_s( pt,                                                   &
1925                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1926                         surf_def_h(2)%shf,                                    &
1927                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1928                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1929                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1930                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1931                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1932                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1933                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1934
1935!
1936!--    Consideration of heat sources within the plant canopy
1937       IF ( plant_canopy  .AND.                                          &
1938            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1939          CALL pcm_tendency( 4 )
1940       ENDIF
1941
1942!
1943!--    Large scale advection
1944       IF ( large_scale_forcing )  THEN
1945          CALL ls_advec( simulated_time, 'pt' )
1946       ENDIF
1947
1948!
1949!--    Nudging
1950       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1951
1952!
1953!--    If required compute influence of large-scale subsidence/ascent
1954       IF ( large_scale_subsidence  .AND.                                      &
1955            .NOT. use_subsidence_tendencies )  THEN
1956          CALL subsidence( tend, pt, pt_init, 2 )
1957       ENDIF
1958
1959!
1960!--    If required, add tendency due to radiative heating/cooling
1961       IF ( radiation  .AND.                                                   &
1962            simulated_time > skip_time_do_radiation )  THEN
1963            CALL radiation_tendency ( tend )
1964       ENDIF
1965
1966       CALL module_interface_actions( 'pt-tendency' )
1967
1968!
1969!--    Prognostic equation for potential temperature
1970       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1971       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1972       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1973       !$ACC PRESENT(tsc(3:5)) &
1974       !$ACC PRESENT(pt_p)
1975       DO  i = nxl, nxr
1976          DO  j = nys, nyn
1977             DO  k = nzb+1, nzt
1978                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1979                                                   tsc(3) * tpt_m(k,j,i) )     &
1980                                                 - tsc(5) *                    &
1981                                                   ( pt(k,j,i) - pt_init(k) ) *&
1982                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1983                                          )                                    &
1984                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1985                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1986                                          )
1987             ENDDO
1988          ENDDO
1989       ENDDO
1990!
1991!--    Calculate tendencies for the next Runge-Kutta step
1992       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1993          IF ( intermediate_timestep_count == 1 )  THEN
1994             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1995             !$ACC PRESENT(tend, tpt_m)
1996             DO  i = nxl, nxr
1997                DO  j = nys, nyn
1998                   DO  k = nzb+1, nzt
1999                      tpt_m(k,j,i) = tend(k,j,i)
2000                   ENDDO
2001                ENDDO
2002             ENDDO
2003          ELSEIF ( intermediate_timestep_count < &
2004                   intermediate_timestep_count_max )  THEN
2005             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2006             !$ACC PRESENT(tend, tpt_m)
2007             DO  i = nxl, nxr
2008                DO  j = nys, nyn
2009                   DO  k = nzb+1, nzt
2010                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
2011                                        5.3125_wp * tpt_m(k,j,i)
2012                   ENDDO
2013                ENDDO
2014             ENDDO
2015          ENDIF
2016       ENDIF
2017
2018       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
2019
2020    ENDIF
2021
2022!
2023!-- If required, compute prognostic equation for total water content
2024    IF ( humidity )  THEN
2025
2026       CALL cpu_log( log_point(29), 'q-equation', 'start' )
2027
2028!
2029!--    Scalar/q-tendency terms with communication
2030       sbt = tsc(2)
2031       IF ( scalar_advec == 'bc-scheme' )  THEN
2032
2033          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2034!
2035!--          Bott-Chlond scheme always uses Euler time step. Thus:
2036             sbt = 1.0_wp
2037          ENDIF
2038          tend = 0.0_wp
2039          CALL advec_s_bc( q, 'q' )
2040
2041       ENDIF
2042
2043!
2044!--    Scalar/q-tendency terms with no communication
2045       IF ( scalar_advec /= 'bc-scheme' )  THEN
2046          tend = 0.0_wp
2047          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2048             IF ( ws_scheme_sca )  THEN
2049                CALL advec_s_ws( q, 'q' )
2050             ELSE
2051                CALL advec_s_pw( q )
2052             ENDIF
2053          ELSE
2054             CALL advec_s_up( q )
2055          ENDIF
2056       ENDIF
2057
2058       CALL diffusion_s( q,                                                    &
2059                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
2060                         surf_def_h(2)%qsws,                                   &
2061                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
2062                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
2063                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
2064                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
2065                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
2066                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
2067                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
2068
2069!
2070!--    Sink or source of humidity due to canopy elements
2071       IF ( plant_canopy ) CALL pcm_tendency( 5 )
2072
2073!
2074!--    Large scale advection
2075       IF ( large_scale_forcing )  THEN
2076          CALL ls_advec( simulated_time, 'q' )
2077       ENDIF
2078
2079!
2080!--    Nudging
2081       IF ( nudging )  CALL nudge( simulated_time, 'q' )
2082
2083!
2084!--    If required compute influence of large-scale subsidence/ascent
2085       IF ( large_scale_subsidence  .AND.                                      &
2086            .NOT. use_subsidence_tendencies )  THEN
2087         CALL subsidence( tend, q, q_init, 3 )
2088       ENDIF
2089
2090       CALL module_interface_actions( 'q-tendency' )
2091
2092!
2093!--    Prognostic equation for total water content
2094       DO  i = nxl, nxr
2095          DO  j = nys, nyn
2096             DO  k = nzb+1, nzt
2097                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2098                                                 tsc(3) * tq_m(k,j,i) )        &
2099                                               - tsc(5) * rdf_sc(k) *          &
2100                                                      ( q(k,j,i) - q_init(k) ) &
2101                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
2102                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
2103                                                 )
2104                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
2105             ENDDO
2106          ENDDO
2107       ENDDO
2108
2109!
2110!--    Calculate tendencies for the next Runge-Kutta step
2111       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2112          IF ( intermediate_timestep_count == 1 )  THEN
2113             DO  i = nxl, nxr
2114                DO  j = nys, nyn
2115                   DO  k = nzb+1, nzt
2116                      tq_m(k,j,i) = tend(k,j,i)
2117                   ENDDO
2118                ENDDO
2119             ENDDO
2120          ELSEIF ( intermediate_timestep_count < &
2121                   intermediate_timestep_count_max )  THEN
2122             DO  i = nxl, nxr
2123                DO  j = nys, nyn
2124                   DO  k = nzb+1, nzt
2125                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2126                                     + 5.3125_wp * tq_m(k,j,i)
2127                   ENDDO
2128                ENDDO
2129             ENDDO
2130          ENDIF
2131       ENDIF
2132
2133       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
2134
2135!
2136!--    If required, calculate prognostic equations for cloud water content
2137!--    and cloud drop concentration
2138       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
2139
2140          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
2141
2142!
2143!--       Calculate prognostic equation for cloud water content
2144          sbt = tsc(2)
2145          IF ( scalar_advec == 'bc-scheme' )  THEN
2146
2147             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2148!
2149!--             Bott-Chlond scheme always uses Euler time step. Thus:
2150                sbt = 1.0_wp
2151             ENDIF
2152             tend = 0.0_wp
2153             CALL advec_s_bc( qc, 'qc' )
2154
2155          ENDIF
2156
2157!
2158!--       qc-tendency terms with no communication
2159          IF ( scalar_advec /= 'bc-scheme' )  THEN
2160             tend = 0.0_wp
2161             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2162                IF ( ws_scheme_sca )  THEN
2163                   CALL advec_s_ws( qc, 'qc' )
2164                ELSE
2165                   CALL advec_s_pw( qc )
2166                ENDIF
2167             ELSE
2168                CALL advec_s_up( qc )
2169             ENDIF
2170          ENDIF
2171
2172          CALL diffusion_s( qc,                                                &
2173                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
2174                            surf_def_h(2)%qcsws,                               &
2175                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
2176                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
2177                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
2178                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
2179                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
2180                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
2181                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
2182
2183!
2184!--       Prognostic equation for cloud water content
2185          DO  i = nxl, nxr
2186             DO  j = nys, nyn
2187                DO  k = nzb+1, nzt
2188                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2189                                                      tsc(3) * tqc_m(k,j,i) )  &
2190                                                    - tsc(5) * rdf_sc(k) *     &
2191                                                               qc(k,j,i)       &
2192                                             )                                 &
2193                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2194                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2195                                          )
2196                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
2197                ENDDO
2198             ENDDO
2199          ENDDO
2200
2201!
2202!--       Calculate tendencies for the next Runge-Kutta step
2203          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2204             IF ( intermediate_timestep_count == 1 )  THEN
2205                DO  i = nxl, nxr
2206                   DO  j = nys, nyn
2207                      DO  k = nzb+1, nzt
2208                         tqc_m(k,j,i) = tend(k,j,i)
2209                      ENDDO
2210                   ENDDO
2211                ENDDO
2212             ELSEIF ( intermediate_timestep_count < &
2213                      intermediate_timestep_count_max )  THEN
2214                DO  i = nxl, nxr
2215                   DO  j = nys, nyn
2216                      DO  k = nzb+1, nzt
2217                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2218                                         + 5.3125_wp * tqc_m(k,j,i)
2219                      ENDDO
2220                   ENDDO
2221                ENDDO
2222             ENDIF
2223          ENDIF
2224
2225          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
2226
2227          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
2228!
2229!--       Calculate prognostic equation for cloud drop concentration
2230          sbt = tsc(2)
2231          IF ( scalar_advec == 'bc-scheme' )  THEN
2232
2233             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2234!
2235!--             Bott-Chlond scheme always uses Euler time step. Thus:
2236                sbt = 1.0_wp
2237             ENDIF
2238             tend = 0.0_wp
2239             CALL advec_s_bc( nc, 'nc' )
2240
2241          ENDIF
2242
2243!
2244!--       nc-tendency terms with no communication
2245          IF ( scalar_advec /= 'bc-scheme' )  THEN
2246             tend = 0.0_wp
2247             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2248                IF ( ws_scheme_sca )  THEN
2249                   CALL advec_s_ws( nc, 'nc' )
2250                ELSE
2251                   CALL advec_s_pw( nc )
2252                ENDIF
2253             ELSE
2254                CALL advec_s_up( nc )
2255             ENDIF
2256          ENDIF
2257
2258          CALL diffusion_s( nc,                                                &
2259                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
2260                            surf_def_h(2)%ncsws,                               &
2261                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
2262                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
2263                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
2264                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
2265                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
2266                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
2267                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
2268
2269!
2270!--       Prognostic equation for cloud drop concentration
2271          DO  i = nxl, nxr
2272             DO  j = nys, nyn
2273                DO  k = nzb+1, nzt
2274                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2275                                                      tsc(3) * tnc_m(k,j,i) )  &
2276                                                    - tsc(5) * rdf_sc(k) *     &
2277                                                               nc(k,j,i)       &
2278                                             )                                 &
2279                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2280                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2281                                          )
2282                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
2283                ENDDO
2284             ENDDO
2285          ENDDO
2286
2287!
2288!--       Calculate tendencies for the next Runge-Kutta step
2289          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2290             IF ( intermediate_timestep_count == 1 )  THEN
2291                DO  i = nxl, nxr
2292                   DO  j = nys, nyn
2293                      DO  k = nzb+1, nzt
2294                         tnc_m(k,j,i) = tend(k,j,i)
2295                      ENDDO
2296                   ENDDO
2297                ENDDO
2298             ELSEIF ( intermediate_timestep_count < &
2299                      intermediate_timestep_count_max )  THEN
2300                DO  i = nxl, nxr
2301                   DO  j = nys, nyn
2302                      DO  k = nzb+1, nzt
2303                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2304                                         + 5.3125_wp * tnc_m(k,j,i)
2305                      ENDDO
2306                   ENDDO
2307                ENDDO
2308             ENDIF
2309          ENDIF
2310
2311          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
2312
2313       ENDIF
2314!
2315!--    If required, calculate prognostic equations for rain water content
2316!--    and rain drop concentration
2317       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
2318
2319          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
2320
2321!
2322!--       Calculate prognostic equation for rain water content
2323          sbt = tsc(2)
2324          IF ( scalar_advec == 'bc-scheme' )  THEN
2325
2326             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2327!
2328!--             Bott-Chlond scheme always uses Euler time step. Thus:
2329                sbt = 1.0_wp
2330             ENDIF
2331             tend = 0.0_wp
2332             CALL advec_s_bc( qr, 'qr' )
2333
2334          ENDIF
2335
2336!
2337!--       qr-tendency terms with no communication
2338          IF ( scalar_advec /= 'bc-scheme' )  THEN
2339             tend = 0.0_wp
2340             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2341                IF ( ws_scheme_sca )  THEN
2342                   CALL advec_s_ws( qr, 'qr' )
2343                ELSE
2344                   CALL advec_s_pw( qr )
2345                ENDIF
2346             ELSE
2347                CALL advec_s_up( qr )
2348             ENDIF
2349          ENDIF
2350
2351          CALL diffusion_s( qr,                                                &
2352                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
2353                            surf_def_h(2)%qrsws,                               &
2354                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
2355                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
2356                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
2357                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
2358                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
2359                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
2360                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
2361
2362!
2363!--       Prognostic equation for rain water content
2364          DO  i = nxl, nxr
2365             DO  j = nys, nyn
2366                DO  k = nzb+1, nzt
2367                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2368                                                      tsc(3) * tqr_m(k,j,i) )  &
2369                                                    - tsc(5) * rdf_sc(k) *     &
2370                                                               qr(k,j,i)       &
2371                                             )                                 &
2372                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2373                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2374                                          )
2375                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
2376                ENDDO
2377             ENDDO
2378          ENDDO
2379
2380!
2381!--       Calculate tendencies for the next Runge-Kutta step
2382          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2383             IF ( intermediate_timestep_count == 1 )  THEN
2384                DO  i = nxl, nxr
2385                   DO  j = nys, nyn
2386                      DO  k = nzb+1, nzt
2387                         tqr_m(k,j,i) = tend(k,j,i)
2388                      ENDDO
2389                   ENDDO
2390                ENDDO
2391             ELSEIF ( intermediate_timestep_count < &
2392                      intermediate_timestep_count_max )  THEN
2393                DO  i = nxl, nxr
2394                   DO  j = nys, nyn
2395                      DO  k = nzb+1, nzt
2396                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2397                                         + 5.3125_wp * tqr_m(k,j,i)
2398                      ENDDO
2399                   ENDDO
2400                ENDDO
2401             ENDIF
2402          ENDIF
2403
2404          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
2405          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
2406
2407!
2408!--       Calculate prognostic equation for rain drop concentration
2409          sbt = tsc(2)
2410          IF ( scalar_advec == 'bc-scheme' )  THEN
2411
2412             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2413!
2414!--             Bott-Chlond scheme always uses Euler time step. Thus:
2415                sbt = 1.0_wp
2416             ENDIF
2417             tend = 0.0_wp
2418             CALL advec_s_bc( nr, 'nr' )
2419
2420          ENDIF
2421
2422!
2423!--       nr-tendency terms with no communication
2424          IF ( scalar_advec /= 'bc-scheme' )  THEN
2425             tend = 0.0_wp
2426             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2427                IF ( ws_scheme_sca )  THEN
2428                   CALL advec_s_ws( nr, 'nr' )
2429                ELSE
2430                   CALL advec_s_pw( nr )
2431                ENDIF
2432             ELSE
2433                CALL advec_s_up( nr )
2434             ENDIF
2435          ENDIF
2436
2437          CALL diffusion_s( nr,                                                &
2438                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
2439                            surf_def_h(2)%nrsws,                               &
2440                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
2441                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
2442                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
2443                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
2444                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
2445                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
2446                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
2447
2448!
2449!--       Prognostic equation for rain drop concentration
2450          DO  i = nxl, nxr
2451             DO  j = nys, nyn
2452                DO  k = nzb+1, nzt
2453                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2454                                                      tsc(3) * tnr_m(k,j,i) )  &
2455                                                    - tsc(5) * rdf_sc(k) *     &
2456                                                               nr(k,j,i)       &
2457                                             )                                 &
2458                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2459                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2460                                          )
2461                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
2462                ENDDO
2463             ENDDO
2464          ENDDO
2465
2466!
2467!--       Calculate tendencies for the next Runge-Kutta step
2468          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2469             IF ( intermediate_timestep_count == 1 )  THEN
2470                DO  i = nxl, nxr
2471                   DO  j = nys, nyn
2472                      DO  k = nzb+1, nzt
2473                         tnr_m(k,j,i) = tend(k,j,i)
2474                      ENDDO
2475                   ENDDO
2476                ENDDO
2477             ELSEIF ( intermediate_timestep_count < &
2478                      intermediate_timestep_count_max )  THEN
2479                DO  i = nxl, nxr
2480                   DO  j = nys, nyn
2481                      DO  k = nzb+1, nzt
2482                         tnr_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2483                                         + 5.3125_wp * tnr_m(k,j,i)
2484                      ENDDO
2485                   ENDDO
2486                ENDDO
2487             ENDIF
2488          ENDIF
2489
2490          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
2491
2492       ENDIF
2493
2494    ENDIF
2495!
2496!-- If required, compute prognostic equation for scalar
2497    IF ( passive_scalar )  THEN
2498
2499       CALL cpu_log( log_point(66), 's-equation', 'start' )
2500
2501!
2502!--    Scalar/q-tendency terms with communication
2503       sbt = tsc(2)
2504       IF ( scalar_advec == 'bc-scheme' )  THEN
2505
2506          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2507!
2508!--          Bott-Chlond scheme always uses Euler time step. Thus:
2509             sbt = 1.0_wp
2510          ENDIF
2511          tend = 0.0_wp
2512          CALL advec_s_bc( s, 's' )
2513
2514       ENDIF
2515
2516!
2517!--    Scalar/q-tendency terms with no communication
2518       IF ( scalar_advec /= 'bc-scheme' )  THEN
2519          tend = 0.0_wp
2520          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2521             IF ( ws_scheme_sca )  THEN
2522                CALL advec_s_ws( s, 's' )
2523             ELSE
2524                CALL advec_s_pw( s )
2525             ENDIF
2526          ELSE
2527             CALL advec_s_up( s )
2528          ENDIF
2529       ENDIF
2530
2531       CALL diffusion_s( s,                                                    &
2532                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
2533                         surf_def_h(2)%ssws,                                   &
2534                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
2535                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
2536                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
2537                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
2538                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
2539                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
2540                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2541
2542!
2543!--    Sink or source of humidity due to canopy elements
2544       IF ( plant_canopy ) CALL pcm_tendency( 7 )
2545
2546!
2547!--    Large scale advection. Not implemented for scalars so far.
2548!        IF ( large_scale_forcing )  THEN
2549!           CALL ls_advec( simulated_time, 'q' )
2550!        ENDIF
2551
2552!
2553!--    Nudging. Not implemented for scalars so far.
2554!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
2555
2556!
2557!--    If required compute influence of large-scale subsidence/ascent.
2558!--    Not implemented for scalars so far.
2559       IF ( large_scale_subsidence  .AND.                                      &
2560            .NOT. use_subsidence_tendencies  .AND.                             &
2561            .NOT. large_scale_forcing )  THEN
2562         CALL subsidence( tend, s, s_init, 3 )
2563       ENDIF
2564
2565       CALL module_interface_actions( 's-tendency' )
2566
2567!
2568!--    Prognostic equation for total water content
2569       DO  i = nxl, nxr
2570          DO  j = nys, nyn
2571             DO  k = nzb+1, nzt
2572                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2573                                                 tsc(3) * ts_m(k,j,i) )        &
2574                                               - tsc(5) * rdf_sc(k) *          &
2575                                                    ( s(k,j,i) - s_init(k) )   &
2576                                        )                                      &
2577                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2578                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2579                                          )
2580                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
2581             ENDDO
2582          ENDDO
2583       ENDDO
2584
2585!
2586!--    Calculate tendencies for the next Runge-Kutta step
2587       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2588          IF ( intermediate_timestep_count == 1 )  THEN
2589             DO  i = nxl, nxr
2590                DO  j = nys, nyn
2591                   DO  k = nzb+1, nzt
2592                      ts_m(k,j,i) = tend(k,j,i)
2593                   ENDDO
2594                ENDDO
2595             ENDDO
2596          ELSEIF ( intermediate_timestep_count < &
2597                   intermediate_timestep_count_max )  THEN
2598             DO  i = nxl, nxr
2599                DO  j = nys, nyn
2600                   DO  k = nzb+1, nzt
2601                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2602                                     + 5.3125_wp * ts_m(k,j,i)
2603                   ENDDO
2604                ENDDO
2605             ENDDO
2606          ENDIF
2607       ENDIF
2608
2609       CALL cpu_log( log_point(66), 's-equation', 'stop' )
2610
2611    ENDIF
2612
2613!
2614!-- Calculate prognostic equations for turbulence closure
2615    CALL tcm_prognostic_equations()
2616
2617!
2618!-- Calculate prognostic equation for chemical quantites
2619    IF ( air_chemistry )  THEN
2620       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2621!
2622!--    Loop over chemical species
2623       DO  lsp = 1, nvar                         
2624          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
2625                                          chem_species(lsp)%conc,              &
2626                                          chem_species(lsp)%tconc_m,           &
2627                                          chem_species(lsp)%conc_pr_init,      &
2628                                          lsp )
2629       ENDDO
2630
2631       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )             
2632    ENDIF   ! Chemicals equations
2633       
2634    IF ( salsa )  THEN
2635       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'start' )
2636!
2637!--          Loop over aerosol size bins: number and mass bins
2638       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
2639       
2640          DO  b = 1, nbins               
2641             sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
2642             CALL salsa_tendency( 'aerosol_number', aerosol_number(b)%conc_p,  &
2643                                   aerosol_number(b)%conc,                     &
2644                                   aerosol_number(b)%tconc_m,                  &
2645                                   b, b, aerosol_number(b)%init )
2646             aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
2647             DO  c = 1, ncc_tot
2648                sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
2649                CALL salsa_tendency( 'aerosol_mass',                           &
2650                                      aerosol_mass((c-1)*nbins+b)%conc_p,      &
2651                                      aerosol_mass((c-1)*nbins+b)%conc,        &
2652                                      aerosol_mass((c-1)*nbins+b)%tconc_m,     &
2653                                      b, c, aerosol_mass((c-1)*nbins+b)%init )
2654                aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
2655             ENDDO
2656          ENDDO
2657          IF ( .NOT. salsa_gases_from_chem )  THEN
2658             DO  g = 1, ngast
2659                sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
2660                CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,         &
2661                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
2662                                      g, g, salsa_gas(g)%init )
2663                salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
2664             ENDDO
2665          ENDIF
2666       
2667       ENDIF
2668       
2669       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'stop' )
2670    ENDIF 
2671
2672!
2673!-- Calculate prognostic equations for the ocean
2674    IF ( ocean_mode )  CALL ocean_prognostic_equations()
2675
2676 END SUBROUTINE prognostic_equations_vector
2677
2678
2679 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.