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

Last change on this file since 3421 was 3386, checked in by gronemeier, 5 years ago

renamed tcm_prognostic to tcm_prognostic_equations

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