source: palm/tags/release-5.0/SOURCE/prognostic_equations.f90 @ 4012

Last change on this file since 4012 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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