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

Last change on this file since 3712 was 3684, checked in by knoop, 5 years ago

Enabled module_interface_actions in time_integration and prognostic_equations

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