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

Last change on this file since 3014 was 3014, checked in by maronga, 5 years ago

series of bugfixes

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