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

Last change on this file since 2629 was 2563, checked in by Giersch, 7 years ago

Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data

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