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

Last change on this file since 3373 was 3355, checked in by knoop, 5 years ago

removed calc_radiation.f90 and related cloud_top_radiation namelist parameter (functionality replaced by radiation model)

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