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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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