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

Last change on this file since 3022 was 3022, checked in by suehring, 3 years ago

Revise recent bugfix in nested runs at left and south boundary; bugfix in advection of u in case of OpenMP parallelization; bugfix in plant transpiration

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