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

Last change on this file since 3467 was 3467, checked in by suehring, 3 years ago

Branch salsa @3446 re-integrated into trunk

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