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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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