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

Last change on this file since 3611 was 3589, checked in by suehring, 5 years ago

Remove erroneous UTF encoding; last commit documented

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