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

Last change on this file since 3312 was 3302, checked in by raasch, 5 years ago

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

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