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

Last change on this file since 3840 was 3840, checked in by knoop, 3 years ago

Moved ocean_prognostic_equations calls into module_interface

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