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

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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