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

Last change on this file since 3458 was 3458, checked in by kanani, 3 years ago

Reintegrated fixes/changes from branch chemistry

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