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

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

Correction of "Former revisions" section

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