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

Last change on this file since 3349 was 3337, checked in by kanani, 5 years ago

reintegrate branch resler to trunk

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