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

Last change on this file since 2728 was 2719, checked in by maronga, 6 years ago

bugfix for last commit

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