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

Last change on this file since 3833 was 3833, checked in by forkel, 3 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

  • 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 3833 2019-03-28 15:04:04Z forkel $
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
436    USE ocean_mod,                                                             &
437        ONLY:  ocean_prognostic_equations, stokes_drift_terms, stokes_force,   &
438               wave_breaking, wave_breaking_term
439
440    USE plant_canopy_model_mod,                                                &
441        ONLY:  cthf, pcm_tendency
442
443#if defined( __rrtmg )
444    USE radiation_model_mod,                                                   &
445        ONLY:  radiation, radiation_tendency,                                  &
446               skip_time_do_radiation
447#endif
448                         
449    USE salsa_mod,                                                             &
450        ONLY:  aerosol_mass, aerosol_number, dt_salsa, last_salsa_time, nbins, &
451               ncc_tot, ngast, salsa_boundary_conds, salsa_diagnostics,        &
452               salsa_driver, salsa_gas, salsa_gases_from_chem, salsa_tendency, &
453               skip_time_do_salsa
454               
455    USE salsa_util_mod,                                                        &
456        ONLY:  sums_salsa_ws_l                             
457   
458    USE statistics,                                                            &
459        ONLY:  hom
460
461    USE subsidence_mod,                                                        &
462        ONLY:  subsidence
463
464    USE surface_mod,                                                           &
465        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
466                surf_usm_v
467
468    USE turbulence_closure_mod,                                                &
469        ONLY:  tcm_prognostic_equations
470
471    USE wind_turbine_model_mod,                                                &
472        ONLY:  wtm_tendencies
473
474
475    PRIVATE
476    PUBLIC prognostic_equations_cache, prognostic_equations_vector
477
478    INTERFACE prognostic_equations_cache
479       MODULE PROCEDURE prognostic_equations_cache
480    END INTERFACE prognostic_equations_cache
481
482    INTERFACE prognostic_equations_vector
483       MODULE PROCEDURE prognostic_equations_vector
484    END INTERFACE prognostic_equations_vector
485
486
487 CONTAINS
488
489
490!------------------------------------------------------------------------------!
491! Description:
492! ------------
493!> Version with one optimized loop over all equations. It is only allowed to
494!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
495!>
496!> Here the calls of most subroutines are embedded in two DO loops over i and j,
497!> so communication between CPUs is not allowed (does not make sense) within
498!> these loops.
499!>
500!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
501!------------------------------------------------------------------------------!
502
503 SUBROUTINE prognostic_equations_cache
504
505
506    IMPLICIT NONE
507   
508    INTEGER(iwp) ::  b                   !< index for aerosol size bins (salsa)
509    INTEGER(iwp) ::  c                   !< index for chemical compounds (salsa)
510    INTEGER(iwp) ::  g                   !< index for gaseous compounds (salsa)
511    INTEGER(iwp) ::  i                   !<
512    INTEGER(iwp) ::  i_omp_start         !<
513    INTEGER(iwp) ::  j                   !<
514    INTEGER(iwp) ::  k                   !<
515!$  INTEGER(iwp) ::  omp_get_thread_num  !<
516    INTEGER(iwp) ::  tn = 0              !<
517
518    LOGICAL      ::  loop_start          !<
519    INTEGER(iwp) ::  lsp
520    INTEGER(iwp) ::  lsp_usr             !< lsp running index for chem spcs
521
522
523!
524!-- Time measurement can only be performed for the whole set of equations
525    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
526
527!
528!-- Calculation of chemical reactions. This is done outside of main loop,
529!-- since exchange of ghost points is required after this update of the
530!-- concentrations of chemical species                                   
531    IF ( air_chemistry )  THEN
532       lsp_usr = 1
533!
534!--    Chemical reactions and deposition
535       IF ( chem_gasphase_on ) THEN
536!
537!--       If required, calculate photolysis frequencies -
538!--       UNFINISHED: Why not before the intermediate timestep loop?
539          IF ( intermediate_timestep_count ==  1 )  THEN
540             CALL photolysis_control
541          ENDIF
542
543          IF ( intermediate_timestep_count == 1 .OR.                        &
544             call_chem_at_all_substeps )  THEN
545
546             CALL cpu_log( log_point_s(19), 'chem.reactions', 'start' )
547             !$OMP PARALLEL PRIVATE (i,j)
548             !$OMP DO schedule(static,1)
549             DO  i = nxl, nxr
550                DO  j = nys, nyn
551                   CALL chem_integrate (i,j)
552                ENDDO
553             ENDDO
554             !$OMP END PARALLEL
555             CALL cpu_log( log_point_s(19), 'chem.reactions', 'stop' )
556
557             IF ( deposition_dry )  THEN
558                CALL cpu_log( log_point_s(24), 'chem.deposition', 'start' )
559                DO  i = nxl, nxr
560                   DO  j = nys, nyn
561                      CALL chem_depo(i,j)
562                   ENDDO
563                ENDDO
564                CALL cpu_log( log_point_s(24), 'chem.deposition', 'stop' )
565             ENDIF
566          ENDIF
567       ENDIF
568!
569!--    Loop over chemical species       
570       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'start' )
571       DO  lsp = 1, nspec
572          CALL exchange_horiz( chem_species(lsp)%conc, nbgp )   
573          lsp_usr = 1 
574          DO WHILE ( TRIM( cs_name( lsp_usr ) ) /= 'novalue' )
575             IF ( TRIM(chem_species(lsp)%name) == TRIM(cs_name(lsp_usr)) )  THEN
576
577                CALL chem_boundary_conds( chem_species(lsp)%conc_p,                                 &
578                                          chem_species(lsp)%conc_pr_init ) 
579             
580             ENDIF
581             lsp_usr = lsp_usr +1
582          ENDDO
583
584
585       ENDDO
586       CALL cpu_log( log_point_s(84), 'chem.exch-horiz', 'stop' )
587
588    ENDIF       
589!
590!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
591!-- step. The exchange of ghost points is required after this update of the
592!-- concentrations of aerosol number and mass
593    IF ( salsa )  THEN
594       
595       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN                 
596          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
597          THEN
598             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
599             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
600             !$OMP DO
601!
602!--          Call salsa processes
603             DO  i = nxl, nxr
604                DO  j = nys, nyn
605                   CALL salsa_diagnostics( i, j )
606                   CALL salsa_driver( i, j, 3 )
607                   CALL salsa_diagnostics( i, j )
608                ENDDO
609             ENDDO
610             
611             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
612             
613             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
614!
615!--          Exchange ghost points and decycle if needed.
616             DO  b = 1, nbins
617                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
618                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
619                                           aerosol_number(b)%init )
620                DO  c = 1, ncc_tot
621                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
622                   CALL salsa_boundary_conds(                                  &
623                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
624                                           aerosol_mass((c-1)*nbins+b)%init )
625                ENDDO
626             ENDDO
627             
628             IF ( .NOT. salsa_gases_from_chem )  THEN
629                DO  g = 1, ngast
630                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
631                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
632                                              salsa_gas(g)%init )
633                ENDDO
634             ENDIF
635             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
636             
637             !$OMP END PARALLEL
638             last_salsa_time = time_since_reference_point
639             
640          ENDIF
641         
642       ENDIF
643       
644    ENDIF
645
646!
647!-- If required, calculate cloud microphysics
648    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
649         ( intermediate_timestep_count == 1  .OR.                              &
650           call_microphysics_at_all_substeps ) )                               &
651    THEN
652       !$OMP PARALLEL PRIVATE (i,j)
653       !$OMP DO
654       DO  i = nxlg, nxrg
655          DO  j = nysg, nyng
656             CALL bcm_actions( i, j )
657           ENDDO
658       ENDDO
659       !$OMP END PARALLEL
660    ENDIF
661
662!
663!-- Loop over all prognostic equations
664!-- b, c ang g added for SALSA
665    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn,b,c,g)
666
667    !$  tn = omp_get_thread_num()
668    loop_start = .TRUE.
669
670    !$OMP DO
671    DO  i = nxl, nxr
672
673!
674!--    Store the first loop index. It differs for each thread and is required
675!--    later in advec_ws
676       IF ( loop_start )  THEN
677          loop_start  = .FALSE.
678          i_omp_start = i
679       ENDIF
680
681       DO  j = nys, nyn
682!
683!--       Tendency terms for u-velocity component. Please note, in case of
684!--       non-cyclic boundary conditions the grid point i=0 is excluded from
685!--       the prognostic equations for the u-component.   
686          IF ( i >= nxlu )  THEN
687
688             tend(:,j,i) = 0.0_wp
689             IF ( timestep_scheme(1:5) == 'runge' )  THEN
690                IF ( ws_scheme_mom )  THEN
691                   CALL advec_u_ws( i, j, i_omp_start, tn )
692                ELSE
693                   CALL advec_u_pw( i, j )
694                ENDIF
695             ELSE
696                CALL advec_u_up( i, j )
697             ENDIF
698             CALL diffusion_u( i, j )
699             CALL coriolis( i, j, 1 )
700             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
701                CALL buoyancy( i, j, pt, 1 )
702             ENDIF
703
704!
705!--          Drag by plant canopy
706             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
707
708!
709!--          External pressure gradient
710             IF ( dp_external )  THEN
711                DO  k = dp_level_ind_b+1, nzt
712                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
713                ENDDO
714             ENDIF
715
716!
717!--          Nudging
718             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
719
720!
721!--          Effect of Stokes drift (in ocean mode only)
722             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
723
724!
725!--          Forces by wind turbines
726             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 1 )
727
728             CALL module_interface_actions( i, j, 'u-tendency' )
729!
730!--          Prognostic equation for u-velocity component
731             DO  k = nzb+1, nzt
732
733                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
734                                            ( tsc(2) * tend(k,j,i) +            &
735                                              tsc(3) * tu_m(k,j,i) )            &
736                                            - tsc(5) * rdf(k)                   &
737                                                     * ( u(k,j,i) - u_init(k) ) &
738                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
739                                                   BTEST( wall_flags_0(k,j,i), 1 )&
740                                                 )
741             ENDDO
742
743!
744!--          Add turbulence generated by wave breaking (in ocean mode only)
745             IF ( wave_breaking  .AND.                                         &
746               intermediate_timestep_count == intermediate_timestep_count_max )&
747             THEN
748                CALL wave_breaking_term( i, j, 1 )
749             ENDIF
750
751!
752!--          Calculate tendencies for the next Runge-Kutta step
753             IF ( timestep_scheme(1:5) == 'runge' )  THEN
754                IF ( intermediate_timestep_count == 1 )  THEN
755                   DO  k = nzb+1, nzt
756                      tu_m(k,j,i) = tend(k,j,i)
757                   ENDDO
758                ELSEIF ( intermediate_timestep_count < &
759                         intermediate_timestep_count_max )  THEN
760                   DO  k = nzb+1, nzt
761                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
762                                     + 5.3125_wp * tu_m(k,j,i)
763                   ENDDO
764                ENDIF
765             ENDIF
766
767          ENDIF
768!
769!--       Tendency terms for v-velocity component. Please note, in case of
770!--       non-cyclic boundary conditions the grid point j=0 is excluded from
771!--       the prognostic equations for the v-component. !--       
772          IF ( j >= nysv )  THEN
773
774             tend(:,j,i) = 0.0_wp
775             IF ( timestep_scheme(1:5) == 'runge' )  THEN
776                IF ( ws_scheme_mom )  THEN
777                    CALL advec_v_ws( i, j, i_omp_start, tn )
778                ELSE
779                    CALL advec_v_pw( i, j )
780                ENDIF
781             ELSE
782                CALL advec_v_up( i, j )
783             ENDIF
784             CALL diffusion_v( i, j )
785             CALL coriolis( i, j, 2 )
786
787!
788!--          Drag by plant canopy
789             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
790
791!
792!--          External pressure gradient
793             IF ( dp_external )  THEN
794                DO  k = dp_level_ind_b+1, nzt
795                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
796                ENDDO
797             ENDIF
798
799!
800!--          Nudging
801             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
802
803!
804!--          Effect of Stokes drift (in ocean mode only)
805             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
806
807!
808!--          Forces by wind turbines
809             IF ( wind_turbine )  CALL wtm_tendencies( i, j, 2 )
810
811             CALL module_interface_actions( i, j, 'v-tendency' )
812!
813!--          Prognostic equation for v-velocity component
814             DO  k = nzb+1, nzt
815                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
816                                            ( tsc(2) * tend(k,j,i) +           &
817                                              tsc(3) * tv_m(k,j,i) )           &
818                                            - tsc(5) * rdf(k)                  &
819                                                     * ( v(k,j,i) - v_init(k) )&
820                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
821                                                   BTEST( wall_flags_0(k,j,i), 2 )&
822                                                 )
823             ENDDO
824
825!
826!--          Add turbulence generated by wave breaking (in ocean mode only)
827             IF ( wave_breaking  .AND.                                         &
828               intermediate_timestep_count == intermediate_timestep_count_max )&
829             THEN
830                CALL wave_breaking_term( i, j, 2 )
831             ENDIF
832
833!
834!--          Calculate tendencies for the next Runge-Kutta step
835             IF ( timestep_scheme(1:5) == 'runge' )  THEN
836                IF ( intermediate_timestep_count == 1 )  THEN
837                   DO  k = nzb+1, nzt
838                      tv_m(k,j,i) = tend(k,j,i)
839                   ENDDO
840                ELSEIF ( intermediate_timestep_count < &
841                         intermediate_timestep_count_max )  THEN
842                   DO  k = nzb+1, nzt
843                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
844                                     + 5.3125_wp * tv_m(k,j,i)
845                   ENDDO
846                ENDIF
847             ENDIF
848
849          ENDIF
850
851!
852!--       Tendency terms for w-velocity component
853          tend(:,j,i) = 0.0_wp
854          IF ( timestep_scheme(1:5) == 'runge' )  THEN
855             IF ( ws_scheme_mom )  THEN
856                CALL advec_w_ws( i, j, i_omp_start, tn )
857             ELSE
858                CALL advec_w_pw( i, j )
859             END IF
860          ELSE
861             CALL advec_w_up( i, j )
862          ENDIF
863          CALL diffusion_w( i, j )
864          CALL coriolis( i, j, 3 )
865
866          IF ( .NOT. neutral )  THEN
867             IF ( ocean_mode )  THEN
868                CALL buoyancy( i, j, rho_ocean, 3 )
869             ELSE
870                IF ( .NOT. humidity )  THEN
871                   CALL buoyancy( i, j, pt, 3 )
872                ELSE
873                   CALL buoyancy( i, j, vpt, 3 )
874                ENDIF
875             ENDIF
876          ENDIF
877
878!
879!--       Drag by plant canopy
880          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
881
882!
883!--       Effect of Stokes drift (in ocean mode only)
884          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
885
886!
887!--       Forces by wind turbines
888          IF ( wind_turbine )  CALL wtm_tendencies( i, j, 3 )
889
890          CALL module_interface_actions( i, j, 'w-tendency' )
891!
892!--       Prognostic equation for w-velocity component
893          DO  k = nzb+1, nzt-1
894             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
895                                             ( tsc(2) * tend(k,j,i) +          &
896                                               tsc(3) * tw_m(k,j,i) )          &
897                                             - tsc(5) * rdf(k) * w(k,j,i)      &
898                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
899                                                BTEST( wall_flags_0(k,j,i), 3 )&
900                                              )
901          ENDDO
902
903!
904!--       Calculate tendencies for the next Runge-Kutta step
905          IF ( timestep_scheme(1:5) == 'runge' )  THEN
906             IF ( intermediate_timestep_count == 1 )  THEN
907                DO  k = nzb+1, nzt-1
908                   tw_m(k,j,i) = tend(k,j,i)
909                ENDDO
910             ELSEIF ( intermediate_timestep_count < &
911                      intermediate_timestep_count_max )  THEN
912                DO  k = nzb+1, nzt-1
913                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
914                                  + 5.3125_wp * tw_m(k,j,i)
915                ENDDO
916             ENDIF
917          ENDIF
918
919!
920!--       If required, compute prognostic equation for potential temperature
921          IF ( .NOT. neutral )  THEN
922!
923!--          Tendency terms for potential temperature
924             tend(:,j,i) = 0.0_wp
925             IF ( timestep_scheme(1:5) == 'runge' )  THEN
926                   IF ( ws_scheme_sca )  THEN
927                      CALL advec_s_ws( i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
928                                       flux_l_pt, diss_l_pt, i_omp_start, tn )
929                   ELSE
930                      CALL advec_s_pw( i, j, pt )
931                   ENDIF
932             ELSE
933                CALL advec_s_up( i, j, pt )
934             ENDIF
935             CALL diffusion_s( i, j, pt,                                       &
936                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
937                               surf_def_h(2)%shf,                              &
938                               surf_lsm_h%shf,    surf_usm_h%shf,              &
939                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
940                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
941                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
942                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
943                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
944                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
945
946!
947!--          Consideration of heat sources within the plant canopy
948             IF ( plant_canopy  .AND.                                          &
949                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
950                CALL pcm_tendency( i, j, 4 )
951             ENDIF
952
953!
954!--          Large scale advection
955             IF ( large_scale_forcing )  THEN
956                CALL ls_advec( i, j, simulated_time, 'pt' )
957             ENDIF
958
959!
960!--          Nudging
961             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
962
963!
964!--          If required, compute effect of large-scale subsidence/ascent
965             IF ( large_scale_subsidence  .AND.                                &
966                  .NOT. use_subsidence_tendencies )  THEN
967                CALL subsidence( i, j, tend, pt, pt_init, 2 )
968             ENDIF
969
970#if defined( __rrtmg )
971!
972!--          If required, add tendency due to radiative heating/cooling
973             IF ( radiation  .AND.                                             &
974                  simulated_time > skip_time_do_radiation )  THEN
975                CALL radiation_tendency ( i, j, tend )
976             ENDIF
977#endif
978
979             CALL module_interface_actions( i, j, 'pt-tendency' )
980!
981!--          Prognostic equation for potential temperature
982             DO  k = nzb+1, nzt
983                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
984                                                  ( tsc(2) * tend(k,j,i) +     &
985                                                    tsc(3) * tpt_m(k,j,i) )    &
986                                                  - tsc(5)                     &
987                                                  * ( pt(k,j,i) - pt_init(k) ) &
988                                                  * ( rdf_sc(k) + ptdf_x(i)    &
989                                                                + ptdf_y(j) )  &
990                                          )                                    &
991                                       * MERGE( 1.0_wp, 0.0_wp,                &
992                                                BTEST( wall_flags_0(k,j,i), 0 )&
993                                              )
994             ENDDO
995
996!
997!--          Calculate tendencies for the next Runge-Kutta step
998             IF ( timestep_scheme(1:5) == 'runge' )  THEN
999                IF ( intermediate_timestep_count == 1 )  THEN
1000                   DO  k = nzb+1, nzt
1001                      tpt_m(k,j,i) = tend(k,j,i)
1002                   ENDDO
1003                ELSEIF ( intermediate_timestep_count < &
1004                         intermediate_timestep_count_max )  THEN
1005                   DO  k = nzb+1, nzt
1006                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1007                                        5.3125_wp * tpt_m(k,j,i)
1008                   ENDDO
1009                ENDIF
1010             ENDIF
1011
1012          ENDIF
1013
1014!
1015!--       If required, compute prognostic equation for total water content
1016          IF ( humidity )  THEN
1017
1018!
1019!--          Tendency-terms for total water content / scalar
1020             tend(:,j,i) = 0.0_wp
1021             IF ( timestep_scheme(1:5) == 'runge' ) &
1022             THEN
1023                IF ( ws_scheme_sca )  THEN
1024                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
1025                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
1026                ELSE
1027                   CALL advec_s_pw( i, j, q )
1028                ENDIF
1029             ELSE
1030                CALL advec_s_up( i, j, q )
1031             ENDIF
1032             CALL diffusion_s( i, j, q,                                        &
1033                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
1034                               surf_def_h(2)%qsws,                             &
1035                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
1036                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
1037                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
1038                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
1039                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
1040                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
1041                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1042
1043!
1044!--          Sink or source of humidity due to canopy elements
1045             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
1046
1047!
1048!--          Large scale advection
1049             IF ( large_scale_forcing )  THEN
1050                CALL ls_advec( i, j, simulated_time, 'q' )
1051             ENDIF
1052
1053!
1054!--          Nudging
1055             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
1056
1057!
1058!--          If required compute influence of large-scale subsidence/ascent
1059             IF ( large_scale_subsidence  .AND.                                &
1060                  .NOT. use_subsidence_tendencies )  THEN
1061                CALL subsidence( i, j, tend, q, q_init, 3 )
1062             ENDIF
1063
1064             CALL module_interface_actions( i, j, 'q-tendency' )
1065
1066!
1067!--          Prognostic equation for total water content / scalar
1068             DO  k = nzb+1, nzt
1069                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
1070                                                ( tsc(2) * tend(k,j,i) +       &
1071                                                  tsc(3) * tq_m(k,j,i) )       &
1072                                                - tsc(5) * rdf_sc(k) *         &
1073                                                      ( q(k,j,i) - q_init(k) ) &
1074                                        )                                      &
1075                                       * MERGE( 1.0_wp, 0.0_wp,                &
1076                                                BTEST( wall_flags_0(k,j,i), 0 )&
1077                                              )               
1078                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1079             ENDDO
1080
1081!
1082!--          Calculate tendencies for the next Runge-Kutta step
1083             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1084                IF ( intermediate_timestep_count == 1 )  THEN
1085                   DO  k = nzb+1, nzt
1086                      tq_m(k,j,i) = tend(k,j,i)
1087                   ENDDO
1088                ELSEIF ( intermediate_timestep_count < &
1089                         intermediate_timestep_count_max )  THEN
1090                   DO  k = nzb+1, nzt
1091                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1092                                       5.3125_wp * tq_m(k,j,i)
1093                   ENDDO
1094                ENDIF
1095             ENDIF
1096
1097!
1098!--          If required, calculate prognostic equations for cloud water content
1099!--          and cloud drop concentration
1100             IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
1101!
1102!--             Calculate prognostic equation for cloud water content
1103                tend(:,j,i) = 0.0_wp
1104                IF ( timestep_scheme(1:5) == 'runge' ) &
1105                THEN
1106                   IF ( ws_scheme_sca )  THEN
1107                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
1108                                       diss_s_qc, flux_l_qc, diss_l_qc, &
1109                                       i_omp_start, tn )
1110                   ELSE
1111                      CALL advec_s_pw( i, j, qc )
1112                   ENDIF
1113                ELSE
1114                   CALL advec_s_up( i, j, qc )
1115                ENDIF
1116                CALL diffusion_s( i, j, qc,                                   &
1117                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
1118                                  surf_def_h(2)%qcsws,                        &
1119                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      & 
1120                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
1121                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
1122                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
1123                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
1124                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
1125                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
1126
1127!
1128!--             Prognostic equation for cloud water content
1129                DO  k = nzb+1, nzt
1130                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
1131                                                      ( tsc(2) * tend(k,j,i) + &
1132                                                        tsc(3) * tqc_m(k,j,i) )&
1133                                                      - tsc(5) * rdf_sc(k)     &
1134                                                               * qc(k,j,i)     &
1135                                             )                                 &
1136                                       * MERGE( 1.0_wp, 0.0_wp,                &
1137                                                BTEST( wall_flags_0(k,j,i), 0 )&
1138                                              ) 
1139                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
1140                ENDDO
1141!
1142!--             Calculate tendencies for the next Runge-Kutta step
1143                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1144                   IF ( intermediate_timestep_count == 1 )  THEN
1145                      DO  k = nzb+1, nzt
1146                         tqc_m(k,j,i) = tend(k,j,i)
1147                      ENDDO
1148                   ELSEIF ( intermediate_timestep_count < &
1149                            intermediate_timestep_count_max )  THEN
1150                      DO  k = nzb+1, nzt
1151                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1152                                           5.3125_wp * tqc_m(k,j,i)
1153                      ENDDO
1154                   ENDIF
1155                ENDIF
1156
1157!
1158!--             Calculate prognostic equation for cloud drop concentration.
1159                tend(:,j,i) = 0.0_wp
1160                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1161                   IF ( ws_scheme_sca )  THEN
1162                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
1163                                    diss_s_nc, flux_l_nc, diss_l_nc, &
1164                                    i_omp_start, tn )
1165                   ELSE
1166                      CALL advec_s_pw( i, j, nc )
1167                   ENDIF
1168                ELSE
1169                   CALL advec_s_up( i, j, nc )
1170                ENDIF
1171                CALL diffusion_s( i, j, nc,                                    &
1172                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
1173                                  surf_def_h(2)%ncsws,                         &
1174                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
1175                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
1176                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
1177                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
1178                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
1179                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
1180                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
1181
1182!
1183!--             Prognostic equation for cloud drop concentration
1184                DO  k = nzb+1, nzt
1185                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
1186                                                      ( tsc(2) * tend(k,j,i) + &
1187                                                        tsc(3) * tnc_m(k,j,i) )&
1188                                                      - tsc(5) * rdf_sc(k)     &
1189                                                               * nc(k,j,i)     &
1190                                             )                                 &
1191                                       * MERGE( 1.0_wp, 0.0_wp,                &
1192                                                BTEST( wall_flags_0(k,j,i), 0 )&
1193                                              )
1194                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
1195                ENDDO
1196!
1197!--             Calculate tendencies for the next Runge-Kutta step
1198                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1199                   IF ( intermediate_timestep_count == 1 )  THEN
1200                      DO  k = nzb+1, nzt
1201                         tnc_m(k,j,i) = tend(k,j,i)
1202                      ENDDO
1203                   ELSEIF ( intermediate_timestep_count < &
1204                            intermediate_timestep_count_max )  THEN
1205                      DO  k = nzb+1, nzt
1206                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1207                                           5.3125_wp * tnc_m(k,j,i)
1208                      ENDDO
1209                   ENDIF
1210                ENDIF
1211
1212             ENDIF
1213!
1214!--          If required, calculate prognostic equations for rain water content
1215!--          and rain drop concentration
1216             IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
1217!
1218!--             Calculate prognostic equation for rain water content
1219                tend(:,j,i) = 0.0_wp
1220                IF ( timestep_scheme(1:5) == 'runge' ) &
1221                THEN
1222                   IF ( ws_scheme_sca )  THEN
1223                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
1224                                       diss_s_qr, flux_l_qr, diss_l_qr, &
1225                                       i_omp_start, tn )
1226                   ELSE
1227                      CALL advec_s_pw( i, j, qr )
1228                   ENDIF
1229                ELSE
1230                   CALL advec_s_up( i, j, qr )
1231                ENDIF
1232                CALL diffusion_s( i, j, qr,                                   &
1233                                  surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,   &
1234                                  surf_def_h(2)%qrsws,                        &
1235                                  surf_lsm_h%qrsws,    surf_usm_h%qrsws,      & 
1236                                  surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,   &
1237                                  surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,   &
1238                                  surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,   &
1239                                  surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,   &
1240                                  surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,   &
1241                                  surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
1242
1243!
1244!--             Prognostic equation for rain water content
1245                DO  k = nzb+1, nzt
1246                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d *                         &
1247                                                      ( tsc(2) * tend(k,j,i) + &
1248                                                        tsc(3) * tqr_m(k,j,i) )&
1249                                                      - tsc(5) * rdf_sc(k)     &
1250                                                               * qr(k,j,i)     &
1251                                             )                                 &
1252                                       * MERGE( 1.0_wp, 0.0_wp,                &
1253                                                BTEST( wall_flags_0(k,j,i), 0 )&
1254                                              ) 
1255                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
1256                ENDDO
1257!
1258!--             Calculate tendencies for the next Runge-Kutta step
1259                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1260                   IF ( intermediate_timestep_count == 1 )  THEN
1261                      DO  k = nzb+1, nzt
1262                         tqr_m(k,j,i) = tend(k,j,i)
1263                      ENDDO
1264                   ELSEIF ( intermediate_timestep_count < &
1265                            intermediate_timestep_count_max )  THEN
1266                      DO  k = nzb+1, nzt
1267                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1268                                           5.3125_wp * tqr_m(k,j,i)
1269                      ENDDO
1270                   ENDIF
1271                ENDIF
1272
1273!
1274!--             Calculate prognostic equation for rain drop concentration.
1275                tend(:,j,i) = 0.0_wp
1276                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1277                   IF ( ws_scheme_sca )  THEN
1278                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
1279                                    diss_s_nr, flux_l_nr, diss_l_nr, &
1280                                    i_omp_start, tn )
1281                   ELSE
1282                      CALL advec_s_pw( i, j, nr )
1283                   ENDIF
1284                ELSE
1285                   CALL advec_s_up( i, j, nr )
1286                ENDIF
1287                CALL diffusion_s( i, j, nr,                                    &
1288                                  surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,    &
1289                                  surf_def_h(2)%nrsws,                         &
1290                                  surf_lsm_h%nrsws,    surf_usm_h%nrsws,       &
1291                                  surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,    &
1292                                  surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,    &
1293                                  surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,    &
1294                                  surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,    &
1295                                  surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,    &
1296                                  surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
1297
1298!
1299!--             Prognostic equation for rain drop concentration
1300                DO  k = nzb+1, nzt
1301                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d *                         &
1302                                                      ( tsc(2) * tend(k,j,i) + &
1303                                                        tsc(3) * tnr_m(k,j,i) )&
1304                                                      - tsc(5) * rdf_sc(k)     &
1305                                                               * nr(k,j,i)     &
1306                                             )                                 &
1307                                       * MERGE( 1.0_wp, 0.0_wp,                &
1308                                                BTEST( wall_flags_0(k,j,i), 0 )&
1309                                              )
1310                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
1311                ENDDO
1312!
1313!--             Calculate tendencies for the next Runge-Kutta step
1314                IF ( timestep_scheme(1:5) == 'runge' )  THEN
1315                   IF ( intermediate_timestep_count == 1 )  THEN
1316                      DO  k = nzb+1, nzt
1317                         tnr_m(k,j,i) = tend(k,j,i)
1318                      ENDDO
1319                   ELSEIF ( intermediate_timestep_count < &
1320                            intermediate_timestep_count_max )  THEN
1321                      DO  k = nzb+1, nzt
1322                         tnr_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
1323                                           5.3125_wp * tnr_m(k,j,i)
1324                      ENDDO
1325                   ENDIF
1326                ENDIF
1327
1328             ENDIF
1329
1330          ENDIF
1331
1332!
1333!--       If required, compute prognostic equation for scalar
1334          IF ( passive_scalar )  THEN
1335!
1336!--          Tendency-terms for total water content / scalar
1337             tend(:,j,i) = 0.0_wp
1338             IF ( timestep_scheme(1:5) == 'runge' ) &
1339             THEN
1340                IF ( ws_scheme_sca )  THEN
1341                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
1342                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
1343                ELSE
1344                   CALL advec_s_pw( i, j, s )
1345                ENDIF
1346             ELSE
1347                CALL advec_s_up( i, j, s )
1348             ENDIF
1349             CALL diffusion_s( i, j, s,                                        &
1350                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1351                               surf_def_h(2)%ssws,                             &
1352                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1353                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1354                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1355                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1356                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1357                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1358                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1359
1360!
1361!--          Sink or source of scalar concentration due to canopy elements
1362             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1363
1364!
1365!--          Large scale advection, still need to be extended for scalars
1366!              IF ( large_scale_forcing )  THEN
1367!                 CALL ls_advec( i, j, simulated_time, 's' )
1368!              ENDIF
1369
1370!
1371!--          Nudging, still need to be extended for scalars
1372!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1373
1374!
1375!--          If required compute influence of large-scale subsidence/ascent.
1376!--          Note, the last argument is of no meaning in this case, as it is
1377!--          only used in conjunction with large_scale_forcing, which is to
1378!--          date not implemented for scalars.
1379             IF ( large_scale_subsidence  .AND.                                &
1380                  .NOT. use_subsidence_tendencies  .AND.                       &
1381                  .NOT. large_scale_forcing )  THEN
1382                CALL subsidence( i, j, tend, s, s_init, 3 )
1383             ENDIF
1384
1385             CALL module_interface_actions( i, j, 's-tendency' )
1386
1387!
1388!--          Prognostic equation for scalar
1389             DO  k = nzb+1, nzt
1390                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1391                                            ( tsc(2) * tend(k,j,i) +           &
1392                                              tsc(3) * ts_m(k,j,i) )           &
1393                                            - tsc(5) * rdf_sc(k)               &
1394                                                     * ( s(k,j,i) - s_init(k) )&
1395                                        )                                      &
1396                                       * MERGE( 1.0_wp, 0.0_wp,                &
1397                                                BTEST( wall_flags_0(k,j,i), 0 )&
1398                                              )
1399                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1400             ENDDO
1401
1402!
1403!--          Calculate tendencies for the next Runge-Kutta step
1404             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1405                IF ( intermediate_timestep_count == 1 )  THEN
1406                   DO  k = nzb+1, nzt
1407                      ts_m(k,j,i) = tend(k,j,i)
1408                   ENDDO
1409                ELSEIF ( intermediate_timestep_count < &
1410                         intermediate_timestep_count_max )  THEN
1411                   DO  k = nzb+1, nzt
1412                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1413                                       5.3125_wp * ts_m(k,j,i)
1414                   ENDDO
1415                ENDIF
1416             ENDIF
1417
1418          ENDIF
1419!
1420!--       Calculate prognostic equations for turbulence closure
1421          CALL tcm_prognostic_equations( i, j, i_omp_start, tn )
1422
1423!
1424!--       Calculate prognostic equation for chemical quantites
1425          IF ( air_chemistry )  THEN
1426             !> TODO: remove time measurement since it slows down performance because it will be called extremely often
1427!
1428!--          Loop over chemical species
1429             DO  lsp = 1, nvar                         
1430                CALL chem_prognostic_equations( chem_species(lsp)%conc_p,      &
1431                                     chem_species(lsp)%conc,                   &
1432                                     chem_species(lsp)%tconc_m,                &
1433                                     chem_species(lsp)%conc_pr_init,           &
1434                                     i, j, i_omp_start, tn, lsp,               &
1435                                     chem_species(lsp)%flux_s_cs,              &
1436                                     chem_species(lsp)%diss_s_cs,              &
1437                                     chem_species(lsp)%flux_l_cs,              &
1438                                     chem_species(lsp)%diss_l_cs )       
1439             ENDDO
1440         
1441          ENDIF   ! Chemical equations
1442!
1443!--       Calculate prognostic equations for the ocean
1444          IF ( ocean_mode )  THEN
1445             CALL ocean_prognostic_equations( i, j, i_omp_start, tn )
1446          ENDIF
1447       
1448          IF ( salsa )  THEN 
1449!
1450!--          Loop over aerosol size bins: number and mass bins
1451             IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1452             
1453                DO  b = 1, nbins               
1454                   sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
1455                   CALL salsa_tendency( 'aerosol_number',                      &
1456                                         aerosol_number(b)%conc_p,             &
1457                                         aerosol_number(b)%conc,               &
1458                                         aerosol_number(b)%tconc_m,            &
1459                                         i, j, i_omp_start, tn, b, b,          &
1460                                         aerosol_number(b)%flux_s,             &
1461                                         aerosol_number(b)%diss_s,             &
1462                                         aerosol_number(b)%flux_l,             &
1463                                         aerosol_number(b)%diss_l,             &
1464                                         aerosol_number(b)%init )
1465                   aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
1466                   DO  c = 1, ncc_tot
1467                      sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
1468                      CALL salsa_tendency( 'aerosol_mass',                     &
1469                                            aerosol_mass((c-1)*nbins+b)%conc_p,&
1470                                            aerosol_mass((c-1)*nbins+b)%conc,  &
1471                                            aerosol_mass((c-1)*nbins+b)%tconc_m,&
1472                                            i, j, i_omp_start, tn, b, c,       &
1473                                            aerosol_mass((c-1)*nbins+b)%flux_s,&
1474                                            aerosol_mass((c-1)*nbins+b)%diss_s,&
1475                                            aerosol_mass((c-1)*nbins+b)%flux_l,&
1476                                            aerosol_mass((c-1)*nbins+b)%diss_l,&
1477                                            aerosol_mass((c-1)*nbins+b)%init )
1478                      aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
1479                   ENDDO
1480                ENDDO
1481                IF ( .NOT. salsa_gases_from_chem )  THEN
1482                   DO  g = 1, ngast
1483                      sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
1484                      CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,   &
1485                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
1486                                      i, j, i_omp_start, tn, g, g,             &
1487                                      salsa_gas(g)%flux_s, salsa_gas(g)%diss_s,&
1488                                      salsa_gas(g)%flux_l, salsa_gas(g)%diss_l,&
1489                                      salsa_gas(g)%init )
1490                      salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
1491                   ENDDO
1492                ENDIF
1493             
1494             ENDIF
1495             
1496          ENDIF       
1497
1498       ENDDO  ! loop over j
1499    ENDDO  ! loop over i
1500    !$OMP END PARALLEL
1501
1502
1503    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1504
1505
1506 END SUBROUTINE prognostic_equations_cache
1507
1508
1509!------------------------------------------------------------------------------!
1510! Description:
1511! ------------
1512!> Version for vector machines
1513!------------------------------------------------------------------------------!
1514
1515 SUBROUTINE prognostic_equations_vector
1516
1517
1518    IMPLICIT NONE
1519
1520    INTEGER(iwp) ::  b     !< index for aerosol size bins (salsa)
1521    INTEGER(iwp) ::  c     !< index for chemical compounds (salsa)
1522    INTEGER(iwp) ::  g     !< index for gaseous compounds (salsa)
1523    INTEGER(iwp) ::  i     !<
1524    INTEGER(iwp) ::  j     !<
1525    INTEGER(iwp) ::  k     !<
1526    INTEGER(iwp) ::  lsp   !< running index for chemical species
1527
1528    REAL(wp)     ::  sbt  !<
1529
1530!
1531!-- Run SALSA and aerosol dynamic processes. SALSA is run with a longer time
1532!-- step. The exchange of ghost points is required after this update of the
1533!-- concentrations of aerosol number and mass
1534    IF ( salsa )  THEN
1535       
1536       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
1537       
1538          IF ( ( time_since_reference_point - last_salsa_time ) >= dt_salsa )  &
1539          THEN
1540         
1541             CALL cpu_log( log_point_s(90), 'salsa processes ', 'start' )
1542             !$OMP PARALLEL PRIVATE (i,j,b,c,g)
1543             !$OMP DO
1544!
1545!--          Call salsa processes
1546             DO  i = nxl, nxr
1547                DO  j = nys, nyn
1548                   CALL salsa_diagnostics( i, j )
1549                   CALL salsa_driver( i, j, 3 )
1550                   CALL salsa_diagnostics( i, j )
1551                ENDDO
1552             ENDDO
1553             
1554             CALL cpu_log( log_point_s(90), 'salsa processes ', 'stop' )
1555             
1556             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'start' )
1557!
1558!--          Exchange ghost points and decycle if needed.
1559             DO  b = 1, nbins
1560                CALL exchange_horiz( aerosol_number(b)%conc, nbgp )
1561                CALL salsa_boundary_conds( aerosol_number(b)%conc_p,           &
1562                                           aerosol_number(b)%init )
1563                DO  c = 1, ncc_tot
1564                   CALL exchange_horiz( aerosol_mass((c-1)*nbins+b)%conc, nbgp )
1565                   CALL salsa_boundary_conds(                                  &
1566                                           aerosol_mass((c-1)*nbins+b)%conc_p, &
1567                                           aerosol_mass((c-1)*nbins+b)%init )
1568                ENDDO
1569             ENDDO
1570             
1571             IF ( .NOT. salsa_gases_from_chem )  THEN
1572                DO  g = 1, ngast
1573                   CALL exchange_horiz( salsa_gas(g)%conc, nbgp )
1574                   CALL salsa_boundary_conds( salsa_gas(g)%conc_p,             &
1575                                              salsa_gas(g)%init )
1576                ENDDO
1577             ENDIF
1578             CALL cpu_log( log_point_s(91), 'salsa exch-horiz ', 'stop' )
1579             
1580             !$OMP END PARALLEL
1581             last_salsa_time = time_since_reference_point
1582             
1583          ENDIF
1584         
1585       ENDIF
1586       
1587    ENDIF
1588
1589!
1590!-- If required, calculate cloud microphysical impacts
1591    IF ( bulk_cloud_model  .AND.  .NOT. microphysics_sat_adjust  .AND.         &
1592         ( intermediate_timestep_count == 1  .OR.                              &
1593           call_microphysics_at_all_substeps )                                 &
1594       )  THEN
1595       CALL cpu_log( log_point(51), 'microphysics', 'start' )
1596       CALL bcm_actions
1597       CALL cpu_log( log_point(51), 'microphysics', 'stop' )
1598    ENDIF
1599
1600!
1601!-- u-velocity component
1602    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1603
1604    !$ACC KERNELS PRESENT(tend)
1605    tend = 0.0_wp
1606    !$ACC END KERNELS
1607    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1608       IF ( ws_scheme_mom )  THEN
1609          CALL advec_u_ws
1610       ELSE
1611          CALL advec_u_pw
1612       ENDIF
1613    ELSE
1614       CALL advec_u_up
1615    ENDIF
1616    CALL diffusion_u
1617    CALL coriolis( 1 )
1618    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1619       CALL buoyancy( pt, 1 )
1620    ENDIF
1621
1622!
1623!-- Drag by plant canopy
1624    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1625
1626!
1627!-- External pressure gradient
1628    IF ( dp_external )  THEN
1629       DO  i = nxlu, nxr
1630          DO  j = nys, nyn
1631             DO  k = dp_level_ind_b+1, nzt
1632                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1633             ENDDO
1634          ENDDO
1635       ENDDO
1636    ENDIF
1637
1638!
1639!-- Nudging
1640    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1641
1642!
1643!-- Effect of Stokes drift (in ocean mode only)
1644    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1645
1646!
1647!-- Forces by wind turbines
1648    IF ( wind_turbine )  CALL wtm_tendencies( 1 )
1649
1650    CALL module_interface_actions( 'u-tendency' )
1651
1652!
1653!-- Prognostic equation for u-velocity component
1654    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1655    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1656    !$ACC PRESENT(tsc(2:5)) &
1657    !$ACC PRESENT(u_p)
1658    DO  i = nxlu, nxr
1659       DO  j = nys, nyn
1660          DO  k = nzb+1, nzt
1661             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1662                                                 tsc(3) * tu_m(k,j,i) )          &
1663                                               - tsc(5) * rdf(k) *               &
1664                                                        ( u(k,j,i) - u_init(k) ) &
1665                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1666                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1667                                              )
1668          ENDDO
1669       ENDDO
1670    ENDDO
1671
1672!
1673!-- Add turbulence generated by wave breaking (in ocean mode only)
1674    IF ( wave_breaking  .AND.                                                  &
1675         intermediate_timestep_count == intermediate_timestep_count_max )      &
1676    THEN
1677       CALL wave_breaking_term( 1 )
1678    ENDIF
1679
1680!
1681!-- Calculate tendencies for the next Runge-Kutta step
1682    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1683       IF ( intermediate_timestep_count == 1 )  THEN
1684          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1685          !$ACC PRESENT(tend, tu_m)
1686          DO  i = nxlu, nxr
1687             DO  j = nys, nyn
1688                DO  k = nzb+1, nzt
1689                   tu_m(k,j,i) = tend(k,j,i)
1690                ENDDO
1691             ENDDO
1692          ENDDO
1693       ELSEIF ( intermediate_timestep_count < &
1694                intermediate_timestep_count_max )  THEN
1695          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1696          !$ACC PRESENT(tend, tu_m)
1697          DO  i = nxlu, nxr
1698             DO  j = nys, nyn
1699                DO  k = nzb+1, nzt
1700                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1701                                   + 5.3125_wp * tu_m(k,j,i)
1702                ENDDO
1703             ENDDO
1704          ENDDO
1705       ENDIF
1706    ENDIF
1707
1708    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1709
1710!
1711!-- v-velocity component
1712    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1713
1714    !$ACC KERNELS PRESENT(tend)
1715    tend = 0.0_wp
1716    !$ACC END KERNELS
1717    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1718       IF ( ws_scheme_mom )  THEN
1719          CALL advec_v_ws
1720       ELSE
1721          CALL advec_v_pw
1722       END IF
1723    ELSE
1724       CALL advec_v_up
1725    ENDIF
1726    CALL diffusion_v
1727    CALL coriolis( 2 )
1728
1729!
1730!-- Drag by plant canopy
1731    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1732
1733!
1734!-- External pressure gradient
1735    IF ( dp_external )  THEN
1736       DO  i = nxl, nxr
1737          DO  j = nysv, nyn
1738             DO  k = dp_level_ind_b+1, nzt
1739                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1740             ENDDO
1741          ENDDO
1742       ENDDO
1743    ENDIF
1744
1745!
1746!-- Nudging
1747    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1748
1749!
1750!-- Effect of Stokes drift (in ocean mode only)
1751    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1752
1753!
1754!-- Forces by wind turbines
1755    IF ( wind_turbine )  CALL wtm_tendencies( 2 )
1756
1757    CALL module_interface_actions( 'v-tendency' )
1758
1759!
1760!-- Prognostic equation for v-velocity component
1761    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1762    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1763    !$ACC PRESENT(tsc(2:5)) &
1764    !$ACC PRESENT(v_p)
1765    DO  i = nxl, nxr
1766       DO  j = nysv, nyn
1767          DO  k = nzb+1, nzt
1768             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1769                                                 tsc(3) * tv_m(k,j,i) )        &
1770                                               - tsc(5) * rdf(k) *             &
1771                                                      ( v(k,j,i) - v_init(k) ) &
1772                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1773                                                BTEST( wall_flags_0(k,j,i), 2 )&
1774                                              )
1775          ENDDO
1776       ENDDO
1777    ENDDO
1778
1779!
1780!-- Add turbulence generated by wave breaking (in ocean mode only)
1781    IF ( wave_breaking  .AND.                                                  &
1782         intermediate_timestep_count == intermediate_timestep_count_max )      &
1783    THEN
1784       CALL wave_breaking_term( 2 )
1785    ENDIF
1786
1787!
1788!-- Calculate tendencies for the next Runge-Kutta step
1789    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1790       IF ( intermediate_timestep_count == 1 )  THEN
1791          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1792          !$ACC PRESENT(tend, tv_m)
1793          DO  i = nxl, nxr
1794             DO  j = nysv, nyn
1795                DO  k = nzb+1, nzt
1796                   tv_m(k,j,i) = tend(k,j,i)
1797                ENDDO
1798             ENDDO
1799          ENDDO
1800       ELSEIF ( intermediate_timestep_count < &
1801                intermediate_timestep_count_max )  THEN
1802          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1803          !$ACC PRESENT(tend, tv_m)
1804          DO  i = nxl, nxr
1805             DO  j = nysv, nyn
1806                DO  k = nzb+1, nzt
1807                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1808                                  + 5.3125_wp * tv_m(k,j,i)
1809                ENDDO
1810             ENDDO
1811          ENDDO
1812       ENDIF
1813    ENDIF
1814
1815    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1816
1817!
1818!-- w-velocity component
1819    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1820
1821    !$ACC KERNELS PRESENT(tend)
1822    tend = 0.0_wp
1823    !$ACC END KERNELS
1824    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1825       IF ( ws_scheme_mom )  THEN
1826          CALL advec_w_ws
1827       ELSE
1828          CALL advec_w_pw
1829       ENDIF
1830    ELSE
1831       CALL advec_w_up
1832    ENDIF
1833    CALL diffusion_w
1834    CALL coriolis( 3 )
1835
1836    IF ( .NOT. neutral )  THEN
1837       IF ( ocean_mode )  THEN
1838          CALL buoyancy( rho_ocean, 3 )
1839       ELSE
1840          IF ( .NOT. humidity )  THEN
1841             CALL buoyancy( pt, 3 )
1842          ELSE
1843             CALL buoyancy( vpt, 3 )
1844          ENDIF
1845       ENDIF
1846    ENDIF
1847
1848!
1849!-- Drag by plant canopy
1850    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1851
1852!
1853!-- Effect of Stokes drift (in ocean mode only)
1854    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1855
1856!
1857!-- Forces by wind turbines
1858    IF ( wind_turbine )  CALL wtm_tendencies( 3 )
1859
1860    CALL module_interface_actions( 'w-tendency' )
1861
1862!
1863!-- Prognostic equation for w-velocity component
1864    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1865    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1866    !$ACC PRESENT(tsc(2:5)) &
1867    !$ACC PRESENT(w_p)
1868    DO  i = nxl, nxr
1869       DO  j = nys, nyn
1870          DO  k = nzb+1, nzt-1
1871             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1872                                                 tsc(3) * tw_m(k,j,i) )        &
1873                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1874                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1875                                                BTEST( wall_flags_0(k,j,i), 3 )&
1876                                              )
1877          ENDDO
1878       ENDDO
1879    ENDDO
1880
1881!
1882!-- Calculate tendencies for the next Runge-Kutta step
1883    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1884       IF ( intermediate_timestep_count == 1 )  THEN
1885          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1886          !$ACC PRESENT(tend, tw_m)
1887          DO  i = nxl, nxr
1888             DO  j = nys, nyn
1889                DO  k = nzb+1, nzt-1
1890                   tw_m(k,j,i) = tend(k,j,i)
1891                ENDDO
1892             ENDDO
1893          ENDDO
1894       ELSEIF ( intermediate_timestep_count < &
1895                intermediate_timestep_count_max )  THEN
1896          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1897          !$ACC PRESENT(tend, tw_m)
1898          DO  i = nxl, nxr
1899             DO  j = nys, nyn
1900                DO  k = nzb+1, nzt-1
1901                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1902                                  + 5.3125_wp * tw_m(k,j,i)
1903                ENDDO
1904             ENDDO
1905          ENDDO
1906       ENDIF
1907    ENDIF
1908
1909    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1910
1911
1912!
1913!-- If required, compute prognostic equation for potential temperature
1914    IF ( .NOT. neutral )  THEN
1915
1916       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1917
1918!
1919!--    pt-tendency terms with communication
1920       sbt = tsc(2)
1921       IF ( scalar_advec == 'bc-scheme' )  THEN
1922
1923          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1924!
1925!--          Bott-Chlond scheme always uses Euler time step. Thus:
1926             sbt = 1.0_wp
1927          ENDIF
1928          tend = 0.0_wp
1929          CALL advec_s_bc( pt, 'pt' )
1930
1931       ENDIF
1932
1933!
1934!--    pt-tendency terms with no communication
1935       IF ( scalar_advec /= 'bc-scheme' )  THEN
1936          !$ACC KERNELS PRESENT(tend)
1937          tend = 0.0_wp
1938          !$ACC END KERNELS
1939          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1940             IF ( ws_scheme_sca )  THEN
1941                CALL advec_s_ws( pt, 'pt' )
1942             ELSE
1943                CALL advec_s_pw( pt )
1944             ENDIF
1945          ELSE
1946             CALL advec_s_up( pt )
1947          ENDIF
1948       ENDIF
1949
1950       CALL diffusion_s( pt,                                                   &
1951                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1952                         surf_def_h(2)%shf,                                    &
1953                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1954                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1955                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1956                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1957                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1958                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1959                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1960
1961!
1962!--    Consideration of heat sources within the plant canopy
1963       IF ( plant_canopy  .AND.                                          &
1964            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1965          CALL pcm_tendency( 4 )
1966       ENDIF
1967
1968!
1969!--    Large scale advection
1970       IF ( large_scale_forcing )  THEN
1971          CALL ls_advec( simulated_time, 'pt' )
1972       ENDIF
1973
1974!
1975!--    Nudging
1976       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1977
1978!
1979!--    If required compute influence of large-scale subsidence/ascent
1980       IF ( large_scale_subsidence  .AND.                                      &
1981            .NOT. use_subsidence_tendencies )  THEN
1982          CALL subsidence( tend, pt, pt_init, 2 )
1983       ENDIF
1984
1985#if defined( __rrtmg )
1986!
1987!--    If required, add tendency due to radiative heating/cooling
1988       IF ( radiation  .AND.                                                   &
1989            simulated_time > skip_time_do_radiation )  THEN
1990            CALL radiation_tendency ( tend )
1991       ENDIF
1992#endif
1993
1994       CALL module_interface_actions( 'pt-tendency' )
1995
1996!
1997!--    Prognostic equation for potential temperature
1998       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1999       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
2000       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
2001       !$ACC PRESENT(tsc(3:5)) &
2002       !$ACC PRESENT(pt_p)
2003       DO  i = nxl, nxr
2004          DO  j = nys, nyn
2005             DO  k = nzb+1, nzt
2006                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
2007                                                   tsc(3) * tpt_m(k,j,i) )     &
2008                                                 - tsc(5) *                    &
2009                                                   ( pt(k,j,i) - pt_init(k) ) *&
2010                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
2011                                          )                                    &
2012                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2013                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2014                                          )
2015             ENDDO
2016          ENDDO
2017       ENDDO
2018!
2019!--    Calculate tendencies for the next Runge-Kutta step
2020       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2021          IF ( intermediate_timestep_count == 1 )  THEN
2022             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2023             !$ACC PRESENT(tend, tpt_m)
2024             DO  i = nxl, nxr
2025                DO  j = nys, nyn
2026                   DO  k = nzb+1, nzt
2027                      tpt_m(k,j,i) = tend(k,j,i)
2028                   ENDDO
2029                ENDDO
2030             ENDDO
2031          ELSEIF ( intermediate_timestep_count < &
2032                   intermediate_timestep_count_max )  THEN
2033             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2034             !$ACC PRESENT(tend, tpt_m)
2035             DO  i = nxl, nxr
2036                DO  j = nys, nyn
2037                   DO  k = nzb+1, nzt
2038                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
2039                                        5.3125_wp * tpt_m(k,j,i)
2040                   ENDDO
2041                ENDDO
2042             ENDDO
2043          ENDIF
2044       ENDIF
2045
2046       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
2047
2048    ENDIF
2049
2050!
2051!-- If required, compute prognostic equation for total water content
2052    IF ( humidity )  THEN
2053
2054       CALL cpu_log( log_point(29), 'q-equation', 'start' )
2055
2056!
2057!--    Scalar/q-tendency terms with communication
2058       sbt = tsc(2)
2059       IF ( scalar_advec == 'bc-scheme' )  THEN
2060
2061          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2062!
2063!--          Bott-Chlond scheme always uses Euler time step. Thus:
2064             sbt = 1.0_wp
2065          ENDIF
2066          tend = 0.0_wp
2067          CALL advec_s_bc( q, 'q' )
2068
2069       ENDIF
2070
2071!
2072!--    Scalar/q-tendency terms with no communication
2073       IF ( scalar_advec /= 'bc-scheme' )  THEN
2074          tend = 0.0_wp
2075          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2076             IF ( ws_scheme_sca )  THEN
2077                CALL advec_s_ws( q, 'q' )
2078             ELSE
2079                CALL advec_s_pw( q )
2080             ENDIF
2081          ELSE
2082             CALL advec_s_up( q )
2083          ENDIF
2084       ENDIF
2085
2086       CALL diffusion_s( q,                                                    &
2087                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
2088                         surf_def_h(2)%qsws,                                   &
2089                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
2090                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
2091                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
2092                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
2093                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
2094                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
2095                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
2096
2097!
2098!--    Sink or source of humidity due to canopy elements
2099       IF ( plant_canopy ) CALL pcm_tendency( 5 )
2100
2101!
2102!--    Large scale advection
2103       IF ( large_scale_forcing )  THEN
2104          CALL ls_advec( simulated_time, 'q' )
2105       ENDIF
2106
2107!
2108!--    Nudging
2109       IF ( nudging )  CALL nudge( simulated_time, 'q' )
2110
2111!
2112!--    If required compute influence of large-scale subsidence/ascent
2113       IF ( large_scale_subsidence  .AND.                                      &
2114            .NOT. use_subsidence_tendencies )  THEN
2115         CALL subsidence( tend, q, q_init, 3 )
2116       ENDIF
2117
2118       CALL module_interface_actions( 'q-tendency' )
2119
2120!
2121!--    Prognostic equation for total water content
2122       DO  i = nxl, nxr
2123          DO  j = nys, nyn
2124             DO  k = nzb+1, nzt
2125                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2126                                                 tsc(3) * tq_m(k,j,i) )        &
2127                                               - tsc(5) * rdf_sc(k) *          &
2128                                                      ( q(k,j,i) - q_init(k) ) &
2129                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
2130                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
2131                                                 )
2132                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
2133             ENDDO
2134          ENDDO
2135       ENDDO
2136
2137!
2138!--    Calculate tendencies for the next Runge-Kutta step
2139       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2140          IF ( intermediate_timestep_count == 1 )  THEN
2141             DO  i = nxl, nxr
2142                DO  j = nys, nyn
2143                   DO  k = nzb+1, nzt
2144                      tq_m(k,j,i) = tend(k,j,i)
2145                   ENDDO
2146                ENDDO
2147             ENDDO
2148          ELSEIF ( intermediate_timestep_count < &
2149                   intermediate_timestep_count_max )  THEN
2150             DO  i = nxl, nxr
2151                DO  j = nys, nyn
2152                   DO  k = nzb+1, nzt
2153                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2154                                     + 5.3125_wp * tq_m(k,j,i)
2155                   ENDDO
2156                ENDDO
2157             ENDDO
2158          ENDIF
2159       ENDIF
2160
2161       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
2162
2163!
2164!--    If required, calculate prognostic equations for cloud water content
2165!--    and cloud drop concentration
2166       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
2167
2168          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
2169
2170!
2171!--       Calculate prognostic equation for cloud water content
2172          sbt = tsc(2)
2173          IF ( scalar_advec == 'bc-scheme' )  THEN
2174
2175             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2176!
2177!--             Bott-Chlond scheme always uses Euler time step. Thus:
2178                sbt = 1.0_wp
2179             ENDIF
2180             tend = 0.0_wp
2181             CALL advec_s_bc( qc, 'qc' )
2182
2183          ENDIF
2184
2185!
2186!--       qc-tendency terms with no communication
2187          IF ( scalar_advec /= 'bc-scheme' )  THEN
2188             tend = 0.0_wp
2189             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2190                IF ( ws_scheme_sca )  THEN
2191                   CALL advec_s_ws( qc, 'qc' )
2192                ELSE
2193                   CALL advec_s_pw( qc )
2194                ENDIF
2195             ELSE
2196                CALL advec_s_up( qc )
2197             ENDIF
2198          ENDIF
2199
2200          CALL diffusion_s( qc,                                                &
2201                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
2202                            surf_def_h(2)%qcsws,                               &
2203                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
2204                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
2205                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
2206                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
2207                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
2208                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
2209                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
2210
2211!
2212!--       Prognostic equation for cloud water content
2213          DO  i = nxl, nxr
2214             DO  j = nys, nyn
2215                DO  k = nzb+1, nzt
2216                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2217                                                      tsc(3) * tqc_m(k,j,i) )  &
2218                                                    - tsc(5) * rdf_sc(k) *     &
2219                                                               qc(k,j,i)       &
2220                                             )                                 &
2221                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2222                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2223                                          )
2224                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
2225                ENDDO
2226             ENDDO
2227          ENDDO
2228
2229!
2230!--       Calculate tendencies for the next Runge-Kutta step
2231          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2232             IF ( intermediate_timestep_count == 1 )  THEN
2233                DO  i = nxl, nxr
2234                   DO  j = nys, nyn
2235                      DO  k = nzb+1, nzt
2236                         tqc_m(k,j,i) = tend(k,j,i)
2237                      ENDDO
2238                   ENDDO
2239                ENDDO
2240             ELSEIF ( intermediate_timestep_count < &
2241                      intermediate_timestep_count_max )  THEN
2242                DO  i = nxl, nxr
2243                   DO  j = nys, nyn
2244                      DO  k = nzb+1, nzt
2245                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2246                                         + 5.3125_wp * tqc_m(k,j,i)
2247                      ENDDO
2248                   ENDDO
2249                ENDDO
2250             ENDIF
2251          ENDIF
2252
2253          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
2254
2255          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
2256!
2257!--       Calculate prognostic equation for cloud drop concentration
2258          sbt = tsc(2)
2259          IF ( scalar_advec == 'bc-scheme' )  THEN
2260
2261             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2262!
2263!--             Bott-Chlond scheme always uses Euler time step. Thus:
2264                sbt = 1.0_wp
2265             ENDIF
2266             tend = 0.0_wp
2267             CALL advec_s_bc( nc, 'nc' )
2268
2269          ENDIF
2270
2271!
2272!--       nc-tendency terms with no communication
2273          IF ( scalar_advec /= 'bc-scheme' )  THEN
2274             tend = 0.0_wp
2275             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2276                IF ( ws_scheme_sca )  THEN
2277                   CALL advec_s_ws( nc, 'nc' )
2278                ELSE
2279                   CALL advec_s_pw( nc )
2280                ENDIF
2281             ELSE
2282                CALL advec_s_up( nc )
2283             ENDIF
2284          ENDIF
2285
2286          CALL diffusion_s( nc,                                                &
2287                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
2288                            surf_def_h(2)%ncsws,                               &
2289                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             & 
2290                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
2291                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
2292                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
2293                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
2294                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
2295                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
2296
2297!
2298!--       Prognostic equation for cloud drop concentration
2299          DO  i = nxl, nxr
2300             DO  j = nys, nyn
2301                DO  k = nzb+1, nzt
2302                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2303                                                      tsc(3) * tnc_m(k,j,i) )  &
2304                                                    - tsc(5) * rdf_sc(k) *     &
2305                                                               nc(k,j,i)       &
2306                                             )                                 &
2307                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2308                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2309                                          )
2310                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
2311                ENDDO
2312             ENDDO
2313          ENDDO
2314
2315!
2316!--       Calculate tendencies for the next Runge-Kutta step
2317          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2318             IF ( intermediate_timestep_count == 1 )  THEN
2319                DO  i = nxl, nxr
2320                   DO  j = nys, nyn
2321                      DO  k = nzb+1, nzt
2322                         tnc_m(k,j,i) = tend(k,j,i)
2323                      ENDDO
2324                   ENDDO
2325                ENDDO
2326             ELSEIF ( intermediate_timestep_count < &
2327                      intermediate_timestep_count_max )  THEN
2328                DO  i = nxl, nxr
2329                   DO  j = nys, nyn
2330                      DO  k = nzb+1, nzt
2331                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2332                                         + 5.3125_wp * tnc_m(k,j,i)
2333                      ENDDO
2334                   ENDDO
2335                ENDDO
2336             ENDIF
2337          ENDIF
2338
2339          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
2340
2341       ENDIF
2342!
2343!--    If required, calculate prognostic equations for rain water content
2344!--    and rain drop concentration
2345       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
2346
2347          CALL cpu_log( log_point(52), 'qr-equation', 'start' )
2348
2349!
2350!--       Calculate prognostic equation for rain water content
2351          sbt = tsc(2)
2352          IF ( scalar_advec == 'bc-scheme' )  THEN
2353
2354             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2355!
2356!--             Bott-Chlond scheme always uses Euler time step. Thus:
2357                sbt = 1.0_wp
2358             ENDIF
2359             tend = 0.0_wp
2360             CALL advec_s_bc( qr, 'qr' )
2361
2362          ENDIF
2363
2364!
2365!--       qr-tendency terms with no communication
2366          IF ( scalar_advec /= 'bc-scheme' )  THEN
2367             tend = 0.0_wp
2368             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2369                IF ( ws_scheme_sca )  THEN
2370                   CALL advec_s_ws( qr, 'qr' )
2371                ELSE
2372                   CALL advec_s_pw( qr )
2373                ENDIF
2374             ELSE
2375                CALL advec_s_up( qr )
2376             ENDIF
2377          ENDIF
2378
2379          CALL diffusion_s( qr,                                                &
2380                            surf_def_h(0)%qrsws, surf_def_h(1)%qrsws,          &
2381                            surf_def_h(2)%qrsws,                               &
2382                            surf_lsm_h%qrsws,    surf_usm_h%qrsws,             &
2383                            surf_def_v(0)%qrsws, surf_def_v(1)%qrsws,          &
2384                            surf_def_v(2)%qrsws, surf_def_v(3)%qrsws,          &
2385                            surf_lsm_v(0)%qrsws, surf_lsm_v(1)%qrsws,          &
2386                            surf_lsm_v(2)%qrsws, surf_lsm_v(3)%qrsws,          &
2387                            surf_usm_v(0)%qrsws, surf_usm_v(1)%qrsws,          &
2388                            surf_usm_v(2)%qrsws, surf_usm_v(3)%qrsws )
2389
2390!
2391!--       Prognostic equation for rain water content
2392          DO  i = nxl, nxr
2393             DO  j = nys, nyn
2394                DO  k = nzb+1, nzt
2395                   qr_p(k,j,i) = qr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2396                                                      tsc(3) * tqr_m(k,j,i) )  &
2397                                                    - tsc(5) * rdf_sc(k) *     &
2398                                                               qr(k,j,i)       &
2399                                             )                                 &
2400                                    * MERGE( 1.0_wp, 0.0_wp,                   &
2401                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2402                                          )
2403                   IF ( qr_p(k,j,i) < 0.0_wp )  qr_p(k,j,i) = 0.0_wp
2404                ENDDO
2405             ENDDO
2406          ENDDO
2407
2408!
2409!--       Calculate tendencies for the next Runge-Kutta step
2410          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2411             IF ( intermediate_timestep_count == 1 )  THEN
2412                DO  i = nxl, nxr
2413                   DO  j = nys, nyn
2414                      DO  k = nzb+1, nzt
2415                         tqr_m(k,j,i) = tend(k,j,i)
2416                      ENDDO
2417                   ENDDO
2418                ENDDO
2419             ELSEIF ( intermediate_timestep_count < &
2420                      intermediate_timestep_count_max )  THEN
2421                DO  i = nxl, nxr
2422                   DO  j = nys, nyn
2423                      DO  k = nzb+1, nzt
2424                         tqr_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
2425                                         + 5.3125_wp * tqr_m(k,j,i)
2426                      ENDDO
2427                   ENDDO
2428                ENDDO
2429             ENDIF
2430          ENDIF
2431
2432          CALL cpu_log( log_point(52), 'qr-equation', 'stop' )
2433          CALL cpu_log( log_point(53), 'nr-equation', 'start' )
2434
2435!
2436!--       Calculate prognostic equation for rain drop concentration
2437          sbt = tsc(2)
2438          IF ( scalar_advec == 'bc-scheme' )  THEN
2439
2440             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2441!
2442!--             Bott-Chlond scheme always uses Euler time step. Thus:
2443                sbt = 1.0_wp
2444             ENDIF
2445             tend = 0.0_wp
2446             CALL advec_s_bc( nr, 'nr' )
2447
2448          ENDIF
2449
2450!
2451!--       nr-tendency terms with no communication
2452          IF ( scalar_advec /= 'bc-scheme' )  THEN
2453             tend = 0.0_wp
2454             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2455                IF ( ws_scheme_sca )  THEN
2456                   CALL advec_s_ws( nr, 'nr' )
2457                ELSE
2458                   CALL advec_s_pw( nr )
2459                ENDIF
2460             ELSE
2461                CALL advec_s_up( nr )
2462             ENDIF
2463          ENDIF
2464
2465          CALL diffusion_s( nr,                                                &
2466                            surf_def_h(0)%nrsws, surf_def_h(1)%nrsws,          &
2467                            surf_def_h(2)%nrsws,                               &
2468                            surf_lsm_h%nrsws,    surf_usm_h%nrsws,             & 
2469                            surf_def_v(0)%nrsws, surf_def_v(1)%nrsws,          &
2470                            surf_def_v(2)%nrsws, surf_def_v(3)%nrsws,          &
2471                            surf_lsm_v(0)%nrsws, surf_lsm_v(1)%nrsws,          &
2472                            surf_lsm_v(2)%nrsws, surf_lsm_v(3)%nrsws,          &
2473                            surf_usm_v(0)%nrsws, surf_usm_v(1)%nrsws,          &
2474                            surf_usm_v(2)%nrsws, surf_usm_v(3)%nrsws )
2475
2476!
2477!--       Prognostic equation for rain drop concentration
2478          DO  i = nxl, nxr
2479             DO  j = nys, nyn
2480                DO  k = nzb+1, nzt
2481                   nr_p(k,j,i) = nr(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
2482                                                      tsc(3) * tnr_m(k,j,i) )  &
2483                                                    - tsc(5) * rdf_sc(k) *     &
2484                                                               nr(k,j,i)       &
2485                                             )                                 &
2486                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2487                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2488                                          )
2489                   IF ( nr_p(k,j,i) < 0.0_wp )  nr_p(k,j,i) = 0.0_wp
2490                ENDDO
2491             ENDDO
2492          ENDDO
2493
2494!
2495!--       Calculate tendencies for the next Runge-Kutta step
2496          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2497             IF ( intermediate_timestep_count == 1 )  THEN
2498                DO  i = nxl, nxr
2499                   DO  j = nys, nyn
2500                      DO  k = nzb+1, nzt
2501                         tnr_m(k,j,i) = tend(k,j,i)
2502                      ENDDO
2503                   ENDDO
2504                ENDDO
2505             ELSEIF ( intermediate_timestep_count < &
2506                      intermediate_timestep_count_max )  THEN
2507                DO  i = nxl, nxr
2508                   DO  j = nys, nyn
2509                      DO  k = nzb+1, nzt
2510                         tnr_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
2511                                         + 5.3125_wp * tnr_m(k,j,i)
2512                      ENDDO
2513                   ENDDO
2514                ENDDO
2515             ENDIF
2516          ENDIF
2517
2518          CALL cpu_log( log_point(53), 'nr-equation', 'stop' )
2519
2520       ENDIF
2521
2522    ENDIF
2523!
2524!-- If required, compute prognostic equation for scalar
2525    IF ( passive_scalar )  THEN
2526
2527       CALL cpu_log( log_point(66), 's-equation', 'start' )
2528
2529!
2530!--    Scalar/q-tendency terms with communication
2531       sbt = tsc(2)
2532       IF ( scalar_advec == 'bc-scheme' )  THEN
2533
2534          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2535!
2536!--          Bott-Chlond scheme always uses Euler time step. Thus:
2537             sbt = 1.0_wp
2538          ENDIF
2539          tend = 0.0_wp
2540          CALL advec_s_bc( s, 's' )
2541
2542       ENDIF
2543
2544!
2545!--    Scalar/q-tendency terms with no communication
2546       IF ( scalar_advec /= 'bc-scheme' )  THEN
2547          tend = 0.0_wp
2548          IF ( timestep_scheme(1:5) == 'runge' )  THEN
2549             IF ( ws_scheme_sca )  THEN
2550                CALL advec_s_ws( s, 's' )
2551             ELSE
2552                CALL advec_s_pw( s )
2553             ENDIF
2554          ELSE
2555             CALL advec_s_up( s )
2556          ENDIF
2557       ENDIF
2558
2559       CALL diffusion_s( s,                                                    &
2560                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
2561                         surf_def_h(2)%ssws,                                   &
2562                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
2563                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
2564                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
2565                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
2566                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
2567                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
2568                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
2569
2570!
2571!--    Sink or source of humidity due to canopy elements
2572       IF ( plant_canopy ) CALL pcm_tendency( 7 )
2573
2574!
2575!--    Large scale advection. Not implemented for scalars so far.
2576!        IF ( large_scale_forcing )  THEN
2577!           CALL ls_advec( simulated_time, 'q' )
2578!        ENDIF
2579
2580!
2581!--    Nudging. Not implemented for scalars so far.
2582!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
2583
2584!
2585!--    If required compute influence of large-scale subsidence/ascent.
2586!--    Not implemented for scalars so far.
2587       IF ( large_scale_subsidence  .AND.                                      &
2588            .NOT. use_subsidence_tendencies  .AND.                             &
2589            .NOT. large_scale_forcing )  THEN
2590         CALL subsidence( tend, s, s_init, 3 )
2591       ENDIF
2592
2593       CALL module_interface_actions( 's-tendency' )
2594
2595!
2596!--    Prognostic equation for total water content
2597       DO  i = nxl, nxr
2598          DO  j = nys, nyn
2599             DO  k = nzb+1, nzt
2600                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2601                                                 tsc(3) * ts_m(k,j,i) )        &
2602                                               - tsc(5) * rdf_sc(k) *          &
2603                                                    ( s(k,j,i) - s_init(k) )   &
2604                                        )                                      &
2605                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2606                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2607                                          )
2608                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
2609             ENDDO
2610          ENDDO
2611       ENDDO
2612
2613!
2614!--    Calculate tendencies for the next Runge-Kutta step
2615       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2616          IF ( intermediate_timestep_count == 1 )  THEN
2617             DO  i = nxl, nxr
2618                DO  j = nys, nyn
2619                   DO  k = nzb+1, nzt
2620                      ts_m(k,j,i) = tend(k,j,i)
2621                   ENDDO
2622                ENDDO
2623             ENDDO
2624          ELSEIF ( intermediate_timestep_count < &
2625                   intermediate_timestep_count_max )  THEN
2626             DO  i = nxl, nxr
2627                DO  j = nys, nyn
2628                   DO  k = nzb+1, nzt
2629                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2630                                     + 5.3125_wp * ts_m(k,j,i)
2631                   ENDDO
2632                ENDDO
2633             ENDDO
2634          ENDIF
2635       ENDIF
2636
2637       CALL cpu_log( log_point(66), 's-equation', 'stop' )
2638
2639    ENDIF
2640
2641!
2642!-- Calculate prognostic equations for turbulence closure
2643    CALL tcm_prognostic_equations()
2644
2645!
2646!-- Calculate prognostic equation for chemical quantites
2647    IF ( air_chemistry )  THEN
2648       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'start' )
2649!
2650!--    Loop over chemical species
2651       DO  lsp = 1, nvar                         
2652          CALL chem_prognostic_equations( chem_species(lsp)%conc_p,            &
2653                                          chem_species(lsp)%conc,              &
2654                                          chem_species(lsp)%tconc_m,           &
2655                                          chem_species(lsp)%conc_pr_init,      &
2656                                          lsp )
2657       ENDDO
2658
2659       CALL cpu_log( log_point_s(25), 'chem.advec+diff+prog', 'stop' )             
2660    ENDIF   ! Chemicals equations
2661       
2662    IF ( salsa )  THEN
2663       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'start' )
2664!
2665!--          Loop over aerosol size bins: number and mass bins
2666       IF ( time_since_reference_point >= skip_time_do_salsa )  THEN
2667       
2668          DO  b = 1, nbins               
2669             sums_salsa_ws_l = aerosol_number(b)%sums_ws_l
2670             CALL salsa_tendency( 'aerosol_number', aerosol_number(b)%conc_p,  &
2671                                   aerosol_number(b)%conc,                     &
2672                                   aerosol_number(b)%tconc_m,                  &
2673                                   b, b, aerosol_number(b)%init )
2674             aerosol_number(b)%sums_ws_l = sums_salsa_ws_l
2675             DO  c = 1, ncc_tot
2676                sums_salsa_ws_l = aerosol_mass((c-1)*nbins+b)%sums_ws_l
2677                CALL salsa_tendency( 'aerosol_mass',                           &
2678                                      aerosol_mass((c-1)*nbins+b)%conc_p,      &
2679                                      aerosol_mass((c-1)*nbins+b)%conc,        &
2680                                      aerosol_mass((c-1)*nbins+b)%tconc_m,     &
2681                                      b, c, aerosol_mass((c-1)*nbins+b)%init )
2682                aerosol_mass((c-1)*nbins+b)%sums_ws_l = sums_salsa_ws_l
2683             ENDDO
2684          ENDDO
2685          IF ( .NOT. salsa_gases_from_chem )  THEN
2686             DO  g = 1, ngast
2687                sums_salsa_ws_l = salsa_gas(g)%sums_ws_l
2688                CALL salsa_tendency( 'salsa_gas', salsa_gas(g)%conc_p,         &
2689                                      salsa_gas(g)%conc, salsa_gas(g)%tconc_m, &
2690                                      g, g, salsa_gas(g)%init )
2691                salsa_gas(g)%sums_ws_l = sums_salsa_ws_l
2692             ENDDO
2693          ENDIF
2694       
2695       ENDIF
2696       
2697       CALL cpu_log( log_point_s(92), 'salsa advec+diff+prog ', 'stop' )
2698    ENDIF 
2699
2700!
2701!-- Calculate prognostic equations for the ocean
2702    IF ( ocean_mode )  CALL ocean_prognostic_equations()
2703
2704 END SUBROUTINE prognostic_equations_vector
2705
2706
2707 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.