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

Last change on this file since 4131 was 4110, checked in by suehring, 5 years ago

last changes documented

  • Property svn:keywords set to Id
File size: 67.8 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: prognostic_equations.f90 4110 2019-07-22 17:05:21Z monakurppa $
27! pass integer flag array to WS scalar advection routine which is now necessary
28! as the flags may differ for scalars, e.g. pt can be cyclic while chemical
29! species may be non-cyclic. Further, pass boundary flags.
30!
31! 4109 2019-07-22 17:00:34Z suehring
32! Application of monotonic flux limiter for the vertical scalar advection
33! up to the topography top (only for the cache-optimized version at the
34! moment). Please note, at the moment the limiter is only applied for passive
35! scalars.
36!
37! 4048 2019-06-21 21:00:21Z knoop
38! Moved tcm_prognostic_equations to module_interface
39!
40! 3987 2019-05-22 09:52:13Z kanani
41! Introduce alternative switch for debug output during timestepping
42!
43! 3956 2019-05-07 12:32:52Z monakurppa
44! Removed salsa calls.
45!
46! 3931 2019-04-24 16:34:28Z schwenkel
47! Correct/complete module_interface introduction for chemistry model
48!
49! 3899 2019-04-16 14:05:27Z monakurppa
50! Corrections in the OpenMP version of salsa
51!
52! 3887 2019 -04-12 08:47:41Z schwenkel
53! Implicit Bugfix for chemistry model, loop for non_transport_physics over
54! ghost points is avoided. Instead introducing module_interface_exchange_horiz.
55!
56! 3885 2019-04-11 11:29:34Z kanani
57! Changes related to global restructuring of location messages and introduction
58! of additional debug messages
59!
60! 3881 2019-04-10 09:31:22Z suehring
61! Bugfix in OpenMP directive
62!
63! 3880 2019-04-08 21:43:02Z knoop
64! Moved wtm_tendencies to module_interface_actions
65!
66! 3874 2019-04-08 16:53:48Z knoop
67! Added non_transport_physics module interfaces and moved bcm code into it
68!
69! 3872 2019-04-08 15:03:06Z knoop
70! Moving prognostic equations of bcm into bulk_cloud_model_mod
71!
72! 3864 2019-04-05 09:01:56Z monakurppa
73! Modifications made for salsa:
74! - salsa_prognostic_equations moved to salsa_mod (and the call to
75!   module_interface_mod)
76! - Renamed nbins --> nbins_aerosol, ncc_tot --> ncomponents_mass and
77!   ngast --> ngases_salsa and loop indices b, c and sg to ib, ic and ig
78!
79! 3840 2019-03-29 10:35:52Z knoop
80! added USE chem_gasphase_mod for nspec, nspec and spc_names
81!
82! 3820 2019-03-27 11:53:41Z forkel
83! renamed do_depo to deposition_dry (ecc)
84!
85! 3797 2019-03-15 11:15:38Z forkel
86! Call chem_integegrate in OpenMP loop   (ketelsen)
87!
88!
89! 3771 2019-02-28 12:19:33Z raasch
90! preprocessor directivs fro rrtmg added
91!
92! 3761 2019-02-25 15:31:42Z raasch
93! unused variable removed
94!
95! 3719 2019-02-06 13:10:18Z kanani
96! Cleaned up chemistry cpu measurements
97!
98! 3684 2019-01-20 20:20:58Z knoop
99! OpenACC port for SPEC
100!
101! 3589 2018-11-30 15:09:51Z suehring
102! Move the control parameter "salsa" from salsa_mod to control_parameters
103! (M. Kurppa)
104!
105! 3582 2018-11-29 19:16:36Z suehring
106! Implementation of a new aerosol module salsa.
107! Remove cpu-logs from i,j loop in cache version.
108!
109! 3458 2018-10-30 14:51:23Z kanani
110! remove duplicate USE chem_modules
111! from chemistry branch r3443, banzhafs, basit:
112! chem_depo call introduced
113! code added for decycling chemistry
114!
115! 3386 2018-10-19 16:28:22Z gronemeier
116! Renamed tcm_prognostic to tcm_prognostic_equations
117!
118! 3355 2018-10-16 14:03:34Z knoop
119! (from branch resler)
120! Fix for chemistry call
121!
122! 3302 2018-10-03 02:39:40Z raasch
123! Stokes drift + wave breaking term added
124!
125! 3298 2018-10-02 12:21:11Z kanani
126! Code added for decycling chemistry (basit)
127!
128! 3294 2018-10-01 02:37:10Z raasch
129! changes concerning modularization of ocean option
130!
131! 3274 2018-09-24 15:42:55Z knoop
132! Modularization of all bulk cloud physics code components
133!
134! 3241 2018-09-12 15:02:00Z raasch
135! omp_get_thread_num now declared in omp directive
136!
137! 3183 2018-07-27 14:25:55Z suehring
138! Remove unused variables from USE statements
139!
140! 3182 2018-07-27 13:36:03Z suehring
141! Revise recent bugfix for nesting
142!
143! 3021 2018-05-16 08:14:20Z maronga
144! Bugfix in IF clause for nesting
145!
146! 3014 2018-05-09 08:42:38Z maronga
147! Fixed a bug in the IF condition to call pcm_tendency in case of
148! potential temperature
149!
150! 2815 2018-02-19 11:29:57Z kanani
151! Rename chem_tendency to chem_prognostic_equations,
152! implement vector version for air chemistry
153!
154! 2766 2018-01-22 17:17:47Z kanani
155! Removed preprocessor directive __chem
156!
157! 2746 2018-01-15 12:06:04Z suehring
158! Move flag plant canopy to modules
159!
160! 2719 2018-01-02 09:02:06Z maronga
161! Bugfix for last change.
162!
163! 2718 2018-01-02 08:49:38Z maronga
164! Corrected "Former revisions" section
165!
166! 2696 2017-12-14 17:12:51Z kanani
167! - Change in file header (GPL part)
168! - Moved TKE equation to tcm_prognostic (TG)
169! - Added switch for chemical reactions (RF, FK)
170! - Implementation of chemistry module (RF, BK, FK)
171!
172! 2563 2017-10-19 15:36:10Z Giersch
173! Variable wind_turbine moved to module control_parameters
174!
175! 2320 2017-07-21 12:47:43Z suehring
176! Modularize large-scale forcing and nudging
177!
178! 2292 2017-06-20 09:51:42Z schwenkel
179! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
180! includes two more prognostic equations for cloud drop concentration (nc) 
181! and cloud water content (qc).
182!
183! 2261 2017-06-08 14:25:57Z raasch
184! bugfix for r2232: openmp directives removed
185!
186! 2233 2017-05-30 18:08:54Z suehring
187!
188! 2232 2017-05-30 17:47:52Z suehring
189! Adjutst to new surface-type structure. Remove call for usm_wall_heat_flux,
190! which is realized directly in diffusion_s now.
191!
192! 2192 2017-03-22 04:14:10Z raasch
193! Bugfix for misplaced and missing openMP directives from r2155
194!
195! 2155 2017-02-21 09:57:40Z hoffmann
196! Bugfix in the calculation of microphysical quantities on ghost points.
197!
198! 2118 2017-01-17 16:38:49Z raasch
199! OpenACC version of subroutine removed
200!
201! 2031 2016-10-21 15:11:58Z knoop
202! renamed variable rho to rho_ocean
203!
204! 2011 2016-09-19 17:29:57Z kanani
205! Flag urban_surface is now defined in module control_parameters.
206!
207! 2007 2016-08-24 15:47:17Z kanani
208! Added pt tendency calculation based on energy balance at urban surfaces
209! (new urban surface model)
210!
211! 2000 2016-08-20 18:09:15Z knoop
212! Forced header and separation lines into 80 columns
213!
214! 1976 2016-07-27 13:28:04Z maronga
215! Simplied calls to radiation model
216!
217! 1960 2016-07-12 16:34:24Z suehring
218! Separate humidity and passive scalar
219!
220! 1914 2016-05-26 14:44:07Z witha
221! Added calls for wind turbine model
222!
223! 1873 2016-04-18 14:50:06Z maronga
224! Module renamed (removed _mod)
225!
226! 1850 2016-04-08 13:29:27Z maronga
227! Module renamed
228!
229! 1826 2016-04-07 12:01:39Z maronga
230! Renamed canopy model calls.
231!
232! 1822 2016-04-07 07:49:42Z hoffmann
233! Kessler microphysics scheme moved to microphysics.
234!
235! 1757 2016-02-22 15:49:32Z maronga
236!
237! 1691 2015-10-26 16:17:44Z maronga
238! Added optional model spin-up without radiation / land surface model calls.
239! Formatting corrections.
240!
241! 1682 2015-10-07 23:56:08Z knoop
242! Code annotations made doxygen readable
243!
244! 1585 2015-04-30 07:05:52Z maronga
245! Added call for temperature tendency calculation due to radiative flux divergence
246!
247! 1517 2015-01-07 19:12:25Z hoffmann
248! advec_s_bc_mod addded, since advec_s_bc is now a module
249!
250! 1484 2014-10-21 10:53:05Z kanani
251! Changes due to new module structure of the plant canopy model:
252! parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
253! Removed double-listing of use_upstream_for_tke in ONLY-list of module
254! control_parameters
255!
256! 1409 2014-05-23 12:11:32Z suehring
257! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
258! This ensures that left-hand side fluxes are also calculated for nxl in that
259! case, even though the solution at nxl is overwritten in boundary_conds()
260!
261! 1398 2014-05-07 11:15:00Z heinze
262! Rayleigh-damping for horizontal velocity components changed: instead of damping
263! against ug and vg, damping against u_init and v_init is used to allow for a
264! homogenized treatment in case of nudging
265!
266! 1380 2014-04-28 12:40:45Z heinze
267! Change order of calls for scalar prognostic quantities:
268! ls_advec -> nudging -> subsidence since initial profiles
269!
270! 1374 2014-04-25 12:55:07Z raasch
271! missing variables added to ONLY lists
272!
273! 1365 2014-04-22 15:03:56Z boeske
274! Calls of ls_advec for large scale advection added,
275! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
276! new argument ls_index added to the calls of subsidence
277! +ls_index
278!
279! 1361 2014-04-16 15:17:48Z hoffmann
280! Two-moment microphysics moved to the start of prognostic equations. This makes
281! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant.
282! Additionally, it is allowed to call the microphysics just once during the time
283! step (not at each sub-time step).
284!
285! Two-moment cloud physics added for vector and accelerator optimization.
286!
287! 1353 2014-04-08 15:21:23Z heinze
288! REAL constants provided with KIND-attribute
289!
290! 1337 2014-03-25 15:11:48Z heinze
291! Bugfix: REAL constants provided with KIND-attribute
292!
293! 1332 2014-03-25 11:59:43Z suehring
294! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
295!
296! 1330 2014-03-24 17:29:32Z suehring
297! In case of SGS-particle velocity advection of TKE is also allowed with
298! dissipative 5th-order scheme.
299!
300! 1320 2014-03-20 08:40:49Z raasch
301! ONLY-attribute added to USE-statements,
302! kind-parameters added to all INTEGER and REAL declaration statements,
303! kinds are defined in new module kinds,
304! old module precision_kind is removed,
305! revision history before 2012 removed,
306! comment fields (!:) to be used for variable explanations added to
307! all variable declaration statements
308!
309! 1318 2014-03-17 13:35:16Z raasch
310! module interfaces removed
311!
312! 1257 2013-11-08 15:18:40Z raasch
313! openacc loop vector clauses removed, independent clauses added
314!
315! 1246 2013-11-01 08:59:45Z heinze
316! enable nudging also for accelerator version
317!
318! 1241 2013-10-30 11:36:58Z heinze
319! usage of nudging enabled (so far not implemented for accelerator version)
320!
321! 1179 2013-06-14 05:57:58Z raasch
322! two arguments removed from routine buoyancy, ref_state updated on device
323!
324! 1128 2013-04-12 06:19:32Z raasch
325! those parts requiring global communication moved to time_integration,
326! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
327! j_north
328!
329! 1115 2013-03-26 18:16:16Z hoffmann
330! optimized cloud physics: calculation of microphysical tendencies transfered
331! to microphysics.f90; qr and nr are only calculated if precipitation is required
332!
333! 1111 2013-03-08 23:54:10Z raasch
334! update directives for prognostic quantities removed
335!
336! 1106 2013-03-04 05:31:38Z raasch
337! small changes in code formatting
338!
339! 1092 2013-02-02 11:24:22Z raasch
340! unused variables removed
341!
342! 1053 2012-11-13 17:11:03Z hoffmann
343! implementation of two new prognostic equations for rain drop concentration (nr)
344! and rain water content (qr)
345!
346! currently, only available for cache loop optimization
347!
348! 1036 2012-10-22 13:43:42Z raasch
349! code put under GPL (PALM 3.9)
350!
351! 1019 2012-09-28 06:46:45Z raasch
352! non-optimized version of prognostic_equations removed
353!
354! 1015 2012-09-27 09:23:24Z raasch
355! new branch prognostic_equations_acc
356! OpenACC statements added + code changes required for GPU optimization
357!
358! 1001 2012-09-13 14:08:46Z raasch
359! all actions concerning leapfrog- and upstream-spline-scheme removed
360!
361! 978 2012-08-09 08:28:32Z fricke
362! km_damp_x and km_damp_y removed in calls of diffusion_u and diffusion_v
363! add ptdf_x, ptdf_y for damping the potential temperature at the inflow
364! boundary in case of non-cyclic lateral boundaries
365! Bugfix: first thread index changes for WS-scheme at the inflow
366!
367! 940 2012-07-09 14:31:00Z raasch
368! temperature equation can be switched off
369!
370! Revision 1.1  2000/04/13 14:56:27  schroeter
371! Initial revision
372!
373!
374! Description:
375! ------------
376!> Solving the prognostic equations.
377!------------------------------------------------------------------------------!
378 MODULE prognostic_equations_mod
379
380    USE advec_s_bc_mod,                                                        &
381        ONLY:  advec_s_bc
382
383    USE advec_s_pw_mod,                                                        &
384        ONLY:  advec_s_pw
385
386    USE advec_s_up_mod,                                                        &
387        ONLY:  advec_s_up
388
389    USE advec_u_pw_mod,                                                        &
390        ONLY:  advec_u_pw
391
392    USE advec_u_up_mod,                                                        &
393        ONLY:  advec_u_up
394
395    USE advec_v_pw_mod,                                                        &
396        ONLY:  advec_v_pw
397
398    USE advec_v_up_mod,                                                        &
399        ONLY:  advec_v_up
400
401    USE advec_w_pw_mod,                                                        &
402        ONLY:  advec_w_pw
403
404    USE advec_w_up_mod,                                                        &
405        ONLY:  advec_w_up
406
407    USE advec_ws,                                                              &
408        ONLY:  advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws
409
410    USE arrays_3d,                                                             &
411        ONLY:  diss_l_e, diss_l_pt, diss_l_q,                                  &
412               diss_l_s, diss_l_sa, diss_s_e,                                  &
413               diss_s_pt, diss_s_q, diss_s_s, diss_s_sa,                       &
414               e, e_p, flux_s_e, flux_s_pt, flux_s_q,                          &
415               flux_s_s, flux_s_sa, flux_l_e,                                  &
416               flux_l_pt, flux_l_q, flux_l_s,                                  &
417               flux_l_sa, pt, ptdf_x, ptdf_y, pt_init,                         &
418               pt_p, prho, q, q_init, q_p, rdf, rdf_sc,                        &
419               ref_state, rho_ocean, s, s_init, s_p, tend, te_m,               &
420               tpt_m, tq_m, ts_m, tu_m, tv_m, tw_m, u,                         &
421               ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
422
423    USE buoyancy_mod,                                                          &
424        ONLY:  buoyancy
425
426    USE control_parameters,                                                    &
427        ONLY:  bc_dirichlet_l,                                                 &
428               bc_dirichlet_n,                                                 &
429               bc_dirichlet_r,                                                 &
430               bc_dirichlet_s,                                                 &
431               bc_radiation_l,                                                 &
432               bc_radiation_n,                                                 &
433               bc_radiation_r,                                                 &
434               bc_radiation_s,                                                 &
435               constant_diffusion,                                             &
436               debug_output_timestep,                                          &
437               dp_external, dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d,    &
438               humidity, intermediate_timestep_count,                          &
439               intermediate_timestep_count_max, large_scale_forcing,           &
440               large_scale_subsidence,                                         &
441               monotonic_limiter_z,                                            &
442               neutral, nudging,                                               &
443               ocean_mode, passive_scalar, plant_canopy, pt_reference,         &
444               scalar_advec, scalar_advec, simulated_time, sloping_surface,    &
445               timestep_scheme, tsc, use_subsidence_tendencies,                &
446               use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
447               ws_scheme_sca, urban_surface, land_surface,                     &
448               time_since_reference_point, salsa
449
450    USE coriolis_mod,                                                          &
451        ONLY:  coriolis
452
453    USE cpulog,                                                                &
454        ONLY:  cpu_log, log_point, log_point_s
455
456    USE diffusion_s_mod,                                                       &
457        ONLY:  diffusion_s
458
459    USE diffusion_u_mod,                                                       &
460        ONLY:  diffusion_u
461
462    USE diffusion_v_mod,                                                       &
463        ONLY:  diffusion_v
464
465    USE diffusion_w_mod,                                                       &
466        ONLY:  diffusion_w
467
468    USE indices,                                                               &
469        ONLY:  advc_flags_s,                                                   &
470               nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
471               nzb, nzt, wall_flags_0
472
473    USE kinds
474
475    USE lsf_nudging_mod,                                                       &
476        ONLY:  ls_advec, nudge
477
478    USE module_interface,                                                      &
479        ONLY:  module_interface_actions, &
480               module_interface_non_advective_processes, &
481               module_interface_exchange_horiz, &
482               module_interface_prognostic_equations
483
484    USE ocean_mod,                                                             &
485        ONLY:  stokes_drift_terms, stokes_force,   &
486               wave_breaking, wave_breaking_term
487
488    USE plant_canopy_model_mod,                                                &
489        ONLY:  cthf, pcm_tendency
490
491#if defined( __rrtmg )
492    USE radiation_model_mod,                                                   &
493        ONLY:  radiation, radiation_tendency,                                  &
494               skip_time_do_radiation
495#endif
496
497    USE statistics,                                                            &
498        ONLY:  hom
499
500    USE subsidence_mod,                                                        &
501        ONLY:  subsidence
502
503    USE surface_mod,                                                           &
504        ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h,    &
505                surf_usm_v
506
507    IMPLICIT NONE
508
509    PRIVATE
510    PUBLIC prognostic_equations_cache, prognostic_equations_vector
511
512    INTERFACE prognostic_equations_cache
513       MODULE PROCEDURE prognostic_equations_cache
514    END INTERFACE prognostic_equations_cache
515
516    INTERFACE prognostic_equations_vector
517       MODULE PROCEDURE prognostic_equations_vector
518    END INTERFACE prognostic_equations_vector
519
520
521 CONTAINS
522
523
524!------------------------------------------------------------------------------!
525! Description:
526! ------------
527!> Version with one optimized loop over all equations. It is only allowed to
528!> be called for the Wicker and Skamarock or Piascek-Williams advection scheme.
529!>
530!> Here the calls of most subroutines are embedded in two DO loops over i and j,
531!> so communication between CPUs is not allowed (does not make sense) within
532!> these loops.
533!>
534!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
535!------------------------------------------------------------------------------!
536
537 SUBROUTINE prognostic_equations_cache
538
539
540    INTEGER(iwp) ::  i                   !<
541    INTEGER(iwp) ::  i_omp_start         !<
542    INTEGER(iwp) ::  j                   !<
543    INTEGER(iwp) ::  k                   !<
544!$  INTEGER(iwp) ::  omp_get_thread_num  !<
545    INTEGER(iwp) ::  tn = 0              !<
546
547    LOGICAL      ::  loop_start          !<
548
549
550    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'start' )
551!
552!-- Time measurement can only be performed for the whole set of equations
553    CALL cpu_log( log_point(32), 'all progn.equations', 'start' )
554
555    !$OMP PARALLEL PRIVATE (i,j)
556    !$OMP DO
557    DO  i = nxl, nxr
558       DO  j = nys, nyn
559!
560!--       Calculate non advective processes for all other modules
561          CALL module_interface_non_advective_processes( i, j )
562       ENDDO
563    ENDDO
564!
565!-- Module Inferface for exchange horiz after non_advective_processes but before
566!-- advection. Therefore, non_advective_processes must not run for ghost points.
567    !$OMP END PARALLEL
568    CALL module_interface_exchange_horiz()
569!
570!-- Loop over all prognostic equations
571    !$OMP PARALLEL PRIVATE (i,i_omp_start,j,k,loop_start,tn)
572
573    !$  tn = omp_get_thread_num()
574    loop_start = .TRUE.
575
576    !$OMP DO
577    DO  i = nxl, nxr
578
579!
580!--    Store the first loop index. It differs for each thread and is required
581!--    later in advec_ws
582       IF ( loop_start )  THEN
583          loop_start  = .FALSE.
584          i_omp_start = i
585       ENDIF
586
587       DO  j = nys, nyn
588!
589!--       Tendency terms for u-velocity component. Please note, in case of
590!--       non-cyclic boundary conditions the grid point i=0 is excluded from
591!--       the prognostic equations for the u-component.
592          IF ( i >= nxlu )  THEN
593
594             tend(:,j,i) = 0.0_wp
595             IF ( timestep_scheme(1:5) == 'runge' )  THEN
596                IF ( ws_scheme_mom )  THEN
597                   CALL advec_u_ws( i, j, i_omp_start, tn )
598                ELSE
599                   CALL advec_u_pw( i, j )
600                ENDIF
601             ELSE
602                CALL advec_u_up( i, j )
603             ENDIF
604             CALL diffusion_u( i, j )
605             CALL coriolis( i, j, 1 )
606             IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
607                CALL buoyancy( i, j, pt, 1 )
608             ENDIF
609
610!
611!--          Drag by plant canopy
612             IF ( plant_canopy )  CALL pcm_tendency( i, j, 1 )
613
614!
615!--          External pressure gradient
616             IF ( dp_external )  THEN
617                DO  k = dp_level_ind_b+1, nzt
618                   tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
619                ENDDO
620             ENDIF
621
622!
623!--          Nudging
624             IF ( nudging )  CALL nudge( i, j, simulated_time, 'u' )
625
626!
627!--          Effect of Stokes drift (in ocean mode only)
628             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 1 )
629
630             CALL module_interface_actions( i, j, 'u-tendency' )
631!
632!--          Prognostic equation for u-velocity component
633             DO  k = nzb+1, nzt
634
635                u_p(k,j,i) = u(k,j,i) + ( dt_3d *                               &
636                                            ( tsc(2) * tend(k,j,i) +            &
637                                              tsc(3) * tu_m(k,j,i) )            &
638                                            - tsc(5) * rdf(k)                   &
639                                                     * ( u(k,j,i) - u_init(k) ) &
640                                        ) * MERGE( 1.0_wp, 0.0_wp,              &
641                                                   BTEST( wall_flags_0(k,j,i), 1 )&
642                                                 )
643             ENDDO
644
645!
646!--          Add turbulence generated by wave breaking (in ocean mode only)
647             IF ( wave_breaking  .AND.                                         &
648               intermediate_timestep_count == intermediate_timestep_count_max )&
649             THEN
650                CALL wave_breaking_term( i, j, 1 )
651             ENDIF
652
653!
654!--          Calculate tendencies for the next Runge-Kutta step
655             IF ( timestep_scheme(1:5) == 'runge' )  THEN
656                IF ( intermediate_timestep_count == 1 )  THEN
657                   DO  k = nzb+1, nzt
658                      tu_m(k,j,i) = tend(k,j,i)
659                   ENDDO
660                ELSEIF ( intermediate_timestep_count < &
661                         intermediate_timestep_count_max )  THEN
662                   DO  k = nzb+1, nzt
663                      tu_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                &
664                                     + 5.3125_wp * tu_m(k,j,i)
665                   ENDDO
666                ENDIF
667             ENDIF
668
669          ENDIF
670!
671!--       Tendency terms for v-velocity component. Please note, in case of
672!--       non-cyclic boundary conditions the grid point j=0 is excluded from
673!--       the prognostic equations for the v-component. !--       
674          IF ( j >= nysv )  THEN
675
676             tend(:,j,i) = 0.0_wp
677             IF ( timestep_scheme(1:5) == 'runge' )  THEN
678                IF ( ws_scheme_mom )  THEN
679                    CALL advec_v_ws( i, j, i_omp_start, tn )
680                ELSE
681                    CALL advec_v_pw( i, j )
682                ENDIF
683             ELSE
684                CALL advec_v_up( i, j )
685             ENDIF
686             CALL diffusion_v( i, j )
687             CALL coriolis( i, j, 2 )
688
689!
690!--          Drag by plant canopy
691             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
692
693!
694!--          External pressure gradient
695             IF ( dp_external )  THEN
696                DO  k = dp_level_ind_b+1, nzt
697                   tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
698                ENDDO
699             ENDIF
700
701!
702!--          Nudging
703             IF ( nudging )  CALL nudge( i, j, simulated_time, 'v' )
704
705!
706!--          Effect of Stokes drift (in ocean mode only)
707             IF ( stokes_force )  CALL stokes_drift_terms( i, j, 2 )
708
709             CALL module_interface_actions( i, j, 'v-tendency' )
710!
711!--          Prognostic equation for v-velocity component
712             DO  k = nzb+1, nzt
713                v_p(k,j,i) = v(k,j,i) + ( dt_3d *                              &
714                                            ( tsc(2) * tend(k,j,i) +           &
715                                              tsc(3) * tv_m(k,j,i) )           &
716                                            - tsc(5) * rdf(k)                  &
717                                                     * ( v(k,j,i) - v_init(k) )&
718                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
719                                                   BTEST( wall_flags_0(k,j,i), 2 )&
720                                                 )
721             ENDDO
722
723!
724!--          Add turbulence generated by wave breaking (in ocean mode only)
725             IF ( wave_breaking  .AND.                                         &
726               intermediate_timestep_count == intermediate_timestep_count_max )&
727             THEN
728                CALL wave_breaking_term( i, j, 2 )
729             ENDIF
730
731!
732!--          Calculate tendencies for the next Runge-Kutta step
733             IF ( timestep_scheme(1:5) == 'runge' )  THEN
734                IF ( intermediate_timestep_count == 1 )  THEN
735                   DO  k = nzb+1, nzt
736                      tv_m(k,j,i) = tend(k,j,i)
737                   ENDDO
738                ELSEIF ( intermediate_timestep_count < &
739                         intermediate_timestep_count_max )  THEN
740                   DO  k = nzb+1, nzt
741                      tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
742                                     + 5.3125_wp * tv_m(k,j,i)
743                   ENDDO
744                ENDIF
745             ENDIF
746
747          ENDIF
748
749!
750!--       Tendency terms for w-velocity component
751          tend(:,j,i) = 0.0_wp
752          IF ( timestep_scheme(1:5) == 'runge' )  THEN
753             IF ( ws_scheme_mom )  THEN
754                CALL advec_w_ws( i, j, i_omp_start, tn )
755             ELSE
756                CALL advec_w_pw( i, j )
757             END IF
758          ELSE
759             CALL advec_w_up( i, j )
760          ENDIF
761          CALL diffusion_w( i, j )
762          CALL coriolis( i, j, 3 )
763
764          IF ( .NOT. neutral )  THEN
765             IF ( ocean_mode )  THEN
766                CALL buoyancy( i, j, rho_ocean, 3 )
767             ELSE
768                IF ( .NOT. humidity )  THEN
769                   CALL buoyancy( i, j, pt, 3 )
770                ELSE
771                   CALL buoyancy( i, j, vpt, 3 )
772                ENDIF
773             ENDIF
774          ENDIF
775
776!
777!--       Drag by plant canopy
778          IF ( plant_canopy )  CALL pcm_tendency( i, j, 3 )
779
780!
781!--       Effect of Stokes drift (in ocean mode only)
782          IF ( stokes_force )  CALL stokes_drift_terms( i, j, 3 )
783
784          CALL module_interface_actions( i, j, 'w-tendency' )
785!
786!--       Prognostic equation for w-velocity component
787          DO  k = nzb+1, nzt-1
788             w_p(k,j,i) = w(k,j,i) + ( dt_3d *                                 &
789                                             ( tsc(2) * tend(k,j,i) +          &
790                                               tsc(3) * tw_m(k,j,i) )          &
791                                             - tsc(5) * rdf(k) * w(k,j,i)      &
792                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
793                                                BTEST( wall_flags_0(k,j,i), 3 )&
794                                              )
795          ENDDO
796
797!
798!--       Calculate tendencies for the next Runge-Kutta step
799          IF ( timestep_scheme(1:5) == 'runge' )  THEN
800             IF ( intermediate_timestep_count == 1 )  THEN
801                DO  k = nzb+1, nzt-1
802                   tw_m(k,j,i) = tend(k,j,i)
803                ENDDO
804             ELSEIF ( intermediate_timestep_count < &
805                      intermediate_timestep_count_max )  THEN
806                DO  k = nzb+1, nzt-1
807                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
808                                  + 5.3125_wp * tw_m(k,j,i)
809                ENDDO
810             ENDIF
811          ENDIF
812
813!
814!--       If required, compute prognostic equation for potential temperature
815          IF ( .NOT. neutral )  THEN
816!
817!--          Tendency terms for potential temperature
818             tend(:,j,i) = 0.0_wp
819             IF ( timestep_scheme(1:5) == 'runge' )  THEN
820                   IF ( ws_scheme_sca )  THEN
821                      CALL advec_s_ws( advc_flags_s,                           &
822                                       i, j, pt, 'pt', flux_s_pt, diss_s_pt,   &
823                                       flux_l_pt, diss_l_pt, i_omp_start, tn,  &
824                                       bc_dirichlet_l  .OR.  bc_radiation_l,   &
825                                       bc_dirichlet_n  .OR.  bc_radiation_n,   &
826                                       bc_dirichlet_r  .OR.  bc_radiation_r,   &
827                                       bc_dirichlet_s  .OR.  bc_radiation_s )
828                   ELSE
829                      CALL advec_s_pw( i, j, pt )
830                   ENDIF
831             ELSE
832                CALL advec_s_up( i, j, pt )
833             ENDIF
834             CALL diffusion_s( i, j, pt,                                       &
835                               surf_def_h(0)%shf, surf_def_h(1)%shf,           &
836                               surf_def_h(2)%shf,                              &
837                               surf_lsm_h%shf,    surf_usm_h%shf,              &
838                               surf_def_v(0)%shf, surf_def_v(1)%shf,           &
839                               surf_def_v(2)%shf, surf_def_v(3)%shf,           &
840                               surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,           &
841                               surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,           &
842                               surf_usm_v(0)%shf, surf_usm_v(1)%shf,           &
843                               surf_usm_v(2)%shf, surf_usm_v(3)%shf )
844
845!
846!--          Consideration of heat sources within the plant canopy
847             IF ( plant_canopy  .AND.                                          &
848                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
849                CALL pcm_tendency( i, j, 4 )
850             ENDIF
851
852!
853!--          Large scale advection
854             IF ( large_scale_forcing )  THEN
855                CALL ls_advec( i, j, simulated_time, 'pt' )
856             ENDIF
857
858!
859!--          Nudging
860             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
861
862!
863!--          If required, compute effect of large-scale subsidence/ascent
864             IF ( large_scale_subsidence  .AND.                                &
865                  .NOT. use_subsidence_tendencies )  THEN
866                CALL subsidence( i, j, tend, pt, pt_init, 2 )
867             ENDIF
868
869#if defined( __rrtmg )
870!
871!--          If required, add tendency due to radiative heating/cooling
872             IF ( radiation  .AND.                                             &
873                  simulated_time > skip_time_do_radiation )  THEN
874                CALL radiation_tendency ( i, j, tend )
875             ENDIF
876#endif
877
878             CALL module_interface_actions( i, j, 'pt-tendency' )
879!
880!--          Prognostic equation for potential temperature
881             DO  k = nzb+1, nzt
882                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d *                            &
883                                                  ( tsc(2) * tend(k,j,i) +     &
884                                                    tsc(3) * tpt_m(k,j,i) )    &
885                                                  - tsc(5)                     &
886                                                  * ( pt(k,j,i) - pt_init(k) ) &
887                                                  * ( rdf_sc(k) + ptdf_x(i)    &
888                                                                + ptdf_y(j) )  &
889                                          )                                    &
890                                       * MERGE( 1.0_wp, 0.0_wp,                &
891                                                BTEST( wall_flags_0(k,j,i), 0 )&
892                                              )
893             ENDDO
894
895!
896!--          Calculate tendencies for the next Runge-Kutta step
897             IF ( timestep_scheme(1:5) == 'runge' )  THEN
898                IF ( intermediate_timestep_count == 1 )  THEN
899                   DO  k = nzb+1, nzt
900                      tpt_m(k,j,i) = tend(k,j,i)
901                   ENDDO
902                ELSEIF ( intermediate_timestep_count < &
903                         intermediate_timestep_count_max )  THEN
904                   DO  k = nzb+1, nzt
905                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
906                                        5.3125_wp * tpt_m(k,j,i)
907                   ENDDO
908                ENDIF
909             ENDIF
910
911          ENDIF
912
913!
914!--       If required, compute prognostic equation for total water content
915          IF ( humidity )  THEN
916
917!
918!--          Tendency-terms for total water content / scalar
919             tend(:,j,i) = 0.0_wp
920             IF ( timestep_scheme(1:5) == 'runge' )                            &
921             THEN
922                IF ( ws_scheme_sca )  THEN
923                   CALL advec_s_ws( advc_flags_s,                              &
924                                    i, j, q, 'q', flux_s_q,                    &
925                                    diss_s_q, flux_l_q, diss_l_q,              &
926                                    i_omp_start, tn,                           &
927                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
928                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
929                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
930                                    bc_dirichlet_s  .OR.  bc_radiation_s )
931                ELSE
932                   CALL advec_s_pw( i, j, q )
933                ENDIF
934             ELSE
935                CALL advec_s_up( i, j, q )
936             ENDIF
937             CALL diffusion_s( i, j, q,                                        &
938                               surf_def_h(0)%qsws, surf_def_h(1)%qsws,         &
939                               surf_def_h(2)%qsws,                             &
940                               surf_lsm_h%qsws,    surf_usm_h%qsws,            &
941                               surf_def_v(0)%qsws, surf_def_v(1)%qsws,         &
942                               surf_def_v(2)%qsws, surf_def_v(3)%qsws,         &
943                               surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,         &
944                               surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,         &
945                               surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,         &
946                               surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
947
948!
949!--          Sink or source of humidity due to canopy elements
950             IF ( plant_canopy )  CALL pcm_tendency( i, j, 5 )
951
952!
953!--          Large scale advection
954             IF ( large_scale_forcing )  THEN
955                CALL ls_advec( i, j, simulated_time, 'q' )
956             ENDIF
957
958!
959!--          Nudging
960             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
961
962!
963!--          If required compute influence of large-scale subsidence/ascent
964             IF ( large_scale_subsidence  .AND.                                &
965                  .NOT. use_subsidence_tendencies )  THEN
966                CALL subsidence( i, j, tend, q, q_init, 3 )
967             ENDIF
968
969             CALL module_interface_actions( i, j, 'q-tendency' )
970
971!
972!--          Prognostic equation for total water content / scalar
973             DO  k = nzb+1, nzt
974                q_p(k,j,i) = q(k,j,i) + ( dt_3d *                              &
975                                                ( tsc(2) * tend(k,j,i) +       &
976                                                  tsc(3) * tq_m(k,j,i) )       &
977                                                - tsc(5) * rdf_sc(k) *         &
978                                                      ( q(k,j,i) - q_init(k) ) &
979                                        )                                      &
980                                       * MERGE( 1.0_wp, 0.0_wp,                &
981                                                BTEST( wall_flags_0(k,j,i), 0 )&
982                                              )               
983                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
984             ENDDO
985
986!
987!--          Calculate tendencies for the next Runge-Kutta step
988             IF ( timestep_scheme(1:5) == 'runge' )  THEN
989                IF ( intermediate_timestep_count == 1 )  THEN
990                   DO  k = nzb+1, nzt
991                      tq_m(k,j,i) = tend(k,j,i)
992                   ENDDO
993                ELSEIF ( intermediate_timestep_count < &
994                         intermediate_timestep_count_max )  THEN
995                   DO  k = nzb+1, nzt
996                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
997                                       5.3125_wp * tq_m(k,j,i)
998                   ENDDO
999                ENDIF
1000             ENDIF
1001
1002          ENDIF
1003
1004!
1005!--       If required, compute prognostic equation for scalar
1006          IF ( passive_scalar )  THEN
1007!
1008!--          Tendency-terms for total water content / scalar
1009             tend(:,j,i) = 0.0_wp
1010             IF ( timestep_scheme(1:5) == 'runge' )                            &
1011             THEN
1012                IF ( ws_scheme_sca )  THEN
1013!
1014!--                For scalar advection apply monotonic flux limiter near
1015!--                topography.
1016                   CALL advec_s_ws( advc_flags_s,                              &
1017                                    i, j, s, 's', flux_s_s,                    &
1018                                    diss_s_s, flux_l_s, diss_l_s, i_omp_start, &
1019                                    tn,                                        &
1020                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
1021                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
1022                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
1023                                    bc_dirichlet_s  .OR.  bc_radiation_s,      &
1024                                    monotonic_limiter_z )
1025                ELSE
1026                   CALL advec_s_pw( i, j, s )
1027                ENDIF
1028             ELSE
1029                CALL advec_s_up( i, j, s )
1030             ENDIF
1031             CALL diffusion_s( i, j, s,                                        &
1032                               surf_def_h(0)%ssws, surf_def_h(1)%ssws,         &
1033                               surf_def_h(2)%ssws,                             &
1034                               surf_lsm_h%ssws,    surf_usm_h%ssws,            &
1035                               surf_def_v(0)%ssws, surf_def_v(1)%ssws,         &
1036                               surf_def_v(2)%ssws, surf_def_v(3)%ssws,         &
1037                               surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,         &
1038                               surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,         &
1039                               surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,         &
1040                               surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1041
1042!
1043!--          Sink or source of scalar concentration due to canopy elements
1044             IF ( plant_canopy )  CALL pcm_tendency( i, j, 7 )
1045
1046!
1047!--          Large scale advection, still need to be extended for scalars
1048!              IF ( large_scale_forcing )  THEN
1049!                 CALL ls_advec( i, j, simulated_time, 's' )
1050!              ENDIF
1051
1052!
1053!--          Nudging, still need to be extended for scalars
1054!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
1055
1056!
1057!--          If required compute influence of large-scale subsidence/ascent.
1058!--          Note, the last argument is of no meaning in this case, as it is
1059!--          only used in conjunction with large_scale_forcing, which is to
1060!--          date not implemented for scalars.
1061             IF ( large_scale_subsidence  .AND.                                &
1062                  .NOT. use_subsidence_tendencies  .AND.                       &
1063                  .NOT. large_scale_forcing )  THEN
1064                CALL subsidence( i, j, tend, s, s_init, 3 )
1065             ENDIF
1066
1067             CALL module_interface_actions( i, j, 's-tendency' )
1068
1069!
1070!--          Prognostic equation for scalar
1071             DO  k = nzb+1, nzt
1072                s_p(k,j,i) = s(k,j,i) + (  dt_3d *                             &
1073                                            ( tsc(2) * tend(k,j,i) +           &
1074                                              tsc(3) * ts_m(k,j,i) )           &
1075                                            - tsc(5) * rdf_sc(k)               &
1076                                                     * ( s(k,j,i) - s_init(k) )&
1077                                        )                                      &
1078                                       * MERGE( 1.0_wp, 0.0_wp,                &
1079                                                BTEST( wall_flags_0(k,j,i), 0 )&
1080                                              )
1081                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1082             ENDDO
1083
1084!
1085!--          Calculate tendencies for the next Runge-Kutta step
1086             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1087                IF ( intermediate_timestep_count == 1 )  THEN
1088                   DO  k = nzb+1, nzt
1089                      ts_m(k,j,i) = tend(k,j,i)
1090                   ENDDO
1091                ELSEIF ( intermediate_timestep_count < &
1092                         intermediate_timestep_count_max )  THEN
1093                   DO  k = nzb+1, nzt
1094                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +               &
1095                                       5.3125_wp * ts_m(k,j,i)
1096                   ENDDO
1097                ENDIF
1098             ENDIF
1099
1100          ENDIF
1101!
1102!--       Calculate prognostic equations for all other modules
1103          CALL module_interface_prognostic_equations( i, j, i_omp_start, tn )
1104
1105       ENDDO  ! loop over j
1106    ENDDO  ! loop over i
1107    !$OMP END PARALLEL
1108
1109
1110    CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
1111
1112    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_cache', 'end' )
1113
1114 END SUBROUTINE prognostic_equations_cache
1115
1116
1117!------------------------------------------------------------------------------!
1118! Description:
1119! ------------
1120!> Version for vector machines
1121!------------------------------------------------------------------------------!
1122
1123 SUBROUTINE prognostic_equations_vector
1124
1125
1126    INTEGER(iwp) ::  i     !<
1127    INTEGER(iwp) ::  j     !<
1128    INTEGER(iwp) ::  k     !<
1129
1130    REAL(wp)     ::  sbt  !<
1131
1132
1133    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'start' )
1134!
1135!-- Calculate non advective processes for all other modules
1136    CALL module_interface_non_advective_processes
1137!
1138!-- Module Inferface for exchange horiz after non_advective_processes but before
1139!-- advection. Therefore, non_advective_processes must not run for ghost points.
1140    CALL module_interface_exchange_horiz()
1141!
1142!-- u-velocity component
1143    CALL cpu_log( log_point(5), 'u-equation', 'start' )
1144
1145    !$ACC KERNELS PRESENT(tend)
1146    tend = 0.0_wp
1147    !$ACC END KERNELS
1148    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1149       IF ( ws_scheme_mom )  THEN
1150          CALL advec_u_ws
1151       ELSE
1152          CALL advec_u_pw
1153       ENDIF
1154    ELSE
1155       CALL advec_u_up
1156    ENDIF
1157    CALL diffusion_u
1158    CALL coriolis( 1 )
1159    IF ( sloping_surface  .AND.  .NOT. neutral )  THEN
1160       CALL buoyancy( pt, 1 )
1161    ENDIF
1162
1163!
1164!-- Drag by plant canopy
1165    IF ( plant_canopy )  CALL pcm_tendency( 1 )
1166
1167!
1168!-- External pressure gradient
1169    IF ( dp_external )  THEN
1170       DO  i = nxlu, nxr
1171          DO  j = nys, nyn
1172             DO  k = dp_level_ind_b+1, nzt
1173                tend(k,j,i) = tend(k,j,i) - dpdxy(1) * dp_smooth_factor(k)
1174             ENDDO
1175          ENDDO
1176       ENDDO
1177    ENDIF
1178
1179!
1180!-- Nudging
1181    IF ( nudging )  CALL nudge( simulated_time, 'u' )
1182
1183!
1184!-- Effect of Stokes drift (in ocean mode only)
1185    IF ( stokes_force )  CALL stokes_drift_terms( 1 )
1186
1187    CALL module_interface_actions( 'u-tendency' )
1188
1189!
1190!-- Prognostic equation for u-velocity component
1191    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1192    !$ACC PRESENT(u, tend, tu_m, u_init, rdf, wall_flags_0) &
1193    !$ACC PRESENT(tsc(2:5)) &
1194    !$ACC PRESENT(u_p)
1195    DO  i = nxlu, nxr
1196       DO  j = nys, nyn
1197          DO  k = nzb+1, nzt
1198             u_p(k,j,i) = u(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +          &
1199                                                 tsc(3) * tu_m(k,j,i) )          &
1200                                               - tsc(5) * rdf(k) *               &
1201                                                        ( u(k,j,i) - u_init(k) ) &
1202                                     ) * MERGE( 1.0_wp, 0.0_wp,                  &
1203                                                 BTEST( wall_flags_0(k,j,i), 1 ) &
1204                                              )
1205          ENDDO
1206       ENDDO
1207    ENDDO
1208
1209!
1210!-- Add turbulence generated by wave breaking (in ocean mode only)
1211    IF ( wave_breaking  .AND.                                                  &
1212         intermediate_timestep_count == intermediate_timestep_count_max )      &
1213    THEN
1214       CALL wave_breaking_term( 1 )
1215    ENDIF
1216
1217!
1218!-- Calculate tendencies for the next Runge-Kutta step
1219    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1220       IF ( intermediate_timestep_count == 1 )  THEN
1221          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1222          !$ACC PRESENT(tend, tu_m)
1223          DO  i = nxlu, nxr
1224             DO  j = nys, nyn
1225                DO  k = nzb+1, nzt
1226                   tu_m(k,j,i) = tend(k,j,i)
1227                ENDDO
1228             ENDDO
1229          ENDDO
1230       ELSEIF ( intermediate_timestep_count < &
1231                intermediate_timestep_count_max )  THEN
1232          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1233          !$ACC PRESENT(tend, tu_m)
1234          DO  i = nxlu, nxr
1235             DO  j = nys, nyn
1236                DO  k = nzb+1, nzt
1237                   tu_m(k,j,i) =    -9.5625_wp * tend(k,j,i)                   &
1238                                   + 5.3125_wp * tu_m(k,j,i)
1239                ENDDO
1240             ENDDO
1241          ENDDO
1242       ENDIF
1243    ENDIF
1244
1245    CALL cpu_log( log_point(5), 'u-equation', 'stop' )
1246
1247!
1248!-- v-velocity component
1249    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1250
1251    !$ACC KERNELS PRESENT(tend)
1252    tend = 0.0_wp
1253    !$ACC END KERNELS
1254    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1255       IF ( ws_scheme_mom )  THEN
1256          CALL advec_v_ws
1257       ELSE
1258          CALL advec_v_pw
1259       END IF
1260    ELSE
1261       CALL advec_v_up
1262    ENDIF
1263    CALL diffusion_v
1264    CALL coriolis( 2 )
1265
1266!
1267!-- Drag by plant canopy
1268    IF ( plant_canopy )  CALL pcm_tendency( 2 )
1269
1270!
1271!-- External pressure gradient
1272    IF ( dp_external )  THEN
1273       DO  i = nxl, nxr
1274          DO  j = nysv, nyn
1275             DO  k = dp_level_ind_b+1, nzt
1276                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1277             ENDDO
1278          ENDDO
1279       ENDDO
1280    ENDIF
1281
1282!
1283!-- Nudging
1284    IF ( nudging )  CALL nudge( simulated_time, 'v' )
1285
1286!
1287!-- Effect of Stokes drift (in ocean mode only)
1288    IF ( stokes_force )  CALL stokes_drift_terms( 2 )
1289
1290    CALL module_interface_actions( 'v-tendency' )
1291
1292!
1293!-- Prognostic equation for v-velocity component
1294    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1295    !$ACC PRESENT(v, tend, tv_m, v_init, rdf, wall_flags_0) &
1296    !$ACC PRESENT(tsc(2:5)) &
1297    !$ACC PRESENT(v_p)
1298    DO  i = nxl, nxr
1299       DO  j = nysv, nyn
1300          DO  k = nzb+1, nzt
1301             v_p(k,j,i) = v(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1302                                                 tsc(3) * tv_m(k,j,i) )        &
1303                                               - tsc(5) * rdf(k) *             &
1304                                                      ( v(k,j,i) - v_init(k) ) &
1305                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1306                                                BTEST( wall_flags_0(k,j,i), 2 )&
1307                                              )
1308          ENDDO
1309       ENDDO
1310    ENDDO
1311
1312!
1313!-- Add turbulence generated by wave breaking (in ocean mode only)
1314    IF ( wave_breaking  .AND.                                                  &
1315         intermediate_timestep_count == intermediate_timestep_count_max )      &
1316    THEN
1317       CALL wave_breaking_term( 2 )
1318    ENDIF
1319
1320!
1321!-- Calculate tendencies for the next Runge-Kutta step
1322    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1323       IF ( intermediate_timestep_count == 1 )  THEN
1324          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1325          !$ACC PRESENT(tend, tv_m)
1326          DO  i = nxl, nxr
1327             DO  j = nysv, nyn
1328                DO  k = nzb+1, nzt
1329                   tv_m(k,j,i) = tend(k,j,i)
1330                ENDDO
1331             ENDDO
1332          ENDDO
1333       ELSEIF ( intermediate_timestep_count < &
1334                intermediate_timestep_count_max )  THEN
1335          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1336          !$ACC PRESENT(tend, tv_m)
1337          DO  i = nxl, nxr
1338             DO  j = nysv, nyn
1339                DO  k = nzb+1, nzt
1340                   tv_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1341                                  + 5.3125_wp * tv_m(k,j,i)
1342                ENDDO
1343             ENDDO
1344          ENDDO
1345       ENDIF
1346    ENDIF
1347
1348    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1349
1350!
1351!-- w-velocity component
1352    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1353
1354    !$ACC KERNELS PRESENT(tend)
1355    tend = 0.0_wp
1356    !$ACC END KERNELS
1357    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1358       IF ( ws_scheme_mom )  THEN
1359          CALL advec_w_ws
1360       ELSE
1361          CALL advec_w_pw
1362       ENDIF
1363    ELSE
1364       CALL advec_w_up
1365    ENDIF
1366    CALL diffusion_w
1367    CALL coriolis( 3 )
1368
1369    IF ( .NOT. neutral )  THEN
1370       IF ( ocean_mode )  THEN
1371          CALL buoyancy( rho_ocean, 3 )
1372       ELSE
1373          IF ( .NOT. humidity )  THEN
1374             CALL buoyancy( pt, 3 )
1375          ELSE
1376             CALL buoyancy( vpt, 3 )
1377          ENDIF
1378       ENDIF
1379    ENDIF
1380
1381!
1382!-- Drag by plant canopy
1383    IF ( plant_canopy )  CALL pcm_tendency( 3 )
1384
1385!
1386!-- Effect of Stokes drift (in ocean mode only)
1387    IF ( stokes_force )  CALL stokes_drift_terms( 3 )
1388
1389    CALL module_interface_actions( 'w-tendency' )
1390
1391!
1392!-- Prognostic equation for w-velocity component
1393    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1394    !$ACC PRESENT(w, tend, tw_m, v_init, rdf, wall_flags_0) &
1395    !$ACC PRESENT(tsc(2:5)) &
1396    !$ACC PRESENT(w_p)
1397    DO  i = nxl, nxr
1398       DO  j = nys, nyn
1399          DO  k = nzb+1, nzt-1
1400             w_p(k,j,i) = w(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +        &
1401                                                 tsc(3) * tw_m(k,j,i) )        &
1402                                               - tsc(5) * rdf(k) * w(k,j,i)    &
1403                                     ) * MERGE( 1.0_wp, 0.0_wp,                &
1404                                                BTEST( wall_flags_0(k,j,i), 3 )&
1405                                              )
1406          ENDDO
1407       ENDDO
1408    ENDDO
1409
1410!
1411!-- Calculate tendencies for the next Runge-Kutta step
1412    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1413       IF ( intermediate_timestep_count == 1 )  THEN
1414          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1415          !$ACC PRESENT(tend, tw_m)
1416          DO  i = nxl, nxr
1417             DO  j = nys, nyn
1418                DO  k = nzb+1, nzt-1
1419                   tw_m(k,j,i) = tend(k,j,i)
1420                ENDDO
1421             ENDDO
1422          ENDDO
1423       ELSEIF ( intermediate_timestep_count < &
1424                intermediate_timestep_count_max )  THEN
1425          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1426          !$ACC PRESENT(tend, tw_m)
1427          DO  i = nxl, nxr
1428             DO  j = nys, nyn
1429                DO  k = nzb+1, nzt-1
1430                   tw_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                    &
1431                                  + 5.3125_wp * tw_m(k,j,i)
1432                ENDDO
1433             ENDDO
1434          ENDDO
1435       ENDIF
1436    ENDIF
1437
1438    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1439
1440
1441!
1442!-- If required, compute prognostic equation for potential temperature
1443    IF ( .NOT. neutral )  THEN
1444
1445       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1446
1447!
1448!--    pt-tendency terms with communication
1449       sbt = tsc(2)
1450       IF ( scalar_advec == 'bc-scheme' )  THEN
1451
1452          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1453!
1454!--          Bott-Chlond scheme always uses Euler time step. Thus:
1455             sbt = 1.0_wp
1456          ENDIF
1457          tend = 0.0_wp
1458          CALL advec_s_bc( pt, 'pt' )
1459
1460       ENDIF
1461
1462!
1463!--    pt-tendency terms with no communication
1464       IF ( scalar_advec /= 'bc-scheme' )  THEN
1465          !$ACC KERNELS PRESENT(tend)
1466          tend = 0.0_wp
1467          !$ACC END KERNELS
1468          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1469             IF ( ws_scheme_sca )  THEN
1470                CALL advec_s_ws( advc_flags_s, pt, 'pt',                       &
1471                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1472                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1473                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1474                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1475             ELSE
1476                CALL advec_s_pw( pt )
1477             ENDIF
1478          ELSE
1479             CALL advec_s_up( pt )
1480          ENDIF
1481       ENDIF
1482
1483       CALL diffusion_s( pt,                                                   &
1484                         surf_def_h(0)%shf, surf_def_h(1)%shf,                 &
1485                         surf_def_h(2)%shf,                                    &
1486                         surf_lsm_h%shf,    surf_usm_h%shf,                    &
1487                         surf_def_v(0)%shf, surf_def_v(1)%shf,                 &
1488                         surf_def_v(2)%shf, surf_def_v(3)%shf,                 &
1489                         surf_lsm_v(0)%shf, surf_lsm_v(1)%shf,                 &
1490                         surf_lsm_v(2)%shf, surf_lsm_v(3)%shf,                 &
1491                         surf_usm_v(0)%shf, surf_usm_v(1)%shf,                 &
1492                         surf_usm_v(2)%shf, surf_usm_v(3)%shf )
1493
1494!
1495!--    Consideration of heat sources within the plant canopy
1496       IF ( plant_canopy  .AND.                                          &
1497            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
1498          CALL pcm_tendency( 4 )
1499       ENDIF
1500
1501!
1502!--    Large scale advection
1503       IF ( large_scale_forcing )  THEN
1504          CALL ls_advec( simulated_time, 'pt' )
1505       ENDIF
1506
1507!
1508!--    Nudging
1509       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
1510
1511!
1512!--    If required compute influence of large-scale subsidence/ascent
1513       IF ( large_scale_subsidence  .AND.                                      &
1514            .NOT. use_subsidence_tendencies )  THEN
1515          CALL subsidence( tend, pt, pt_init, 2 )
1516       ENDIF
1517
1518#if defined( __rrtmg )
1519!
1520!--    If required, add tendency due to radiative heating/cooling
1521       IF ( radiation  .AND.                                                   &
1522            simulated_time > skip_time_do_radiation )  THEN
1523            CALL radiation_tendency ( tend )
1524       ENDIF
1525#endif
1526
1527       CALL module_interface_actions( 'pt-tendency' )
1528
1529!
1530!--    Prognostic equation for potential temperature
1531       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1532       !$ACC PRESENT(pt, tend, tpt_m, wall_flags_0) &
1533       !$ACC PRESENT(pt_init, rdf_sc, ptdf_x, ptdf_y) &
1534       !$ACC PRESENT(tsc(3:5)) &
1535       !$ACC PRESENT(pt_p)
1536       DO  i = nxl, nxr
1537          DO  j = nys, nyn
1538             DO  k = nzb+1, nzt
1539                pt_p(k,j,i) = pt(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +      &
1540                                                   tsc(3) * tpt_m(k,j,i) )     &
1541                                                 - tsc(5) *                    &
1542                                                   ( pt(k,j,i) - pt_init(k) ) *&
1543                                          ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )&
1544                                          )                                    &
1545                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1546                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1547                                          )
1548             ENDDO
1549          ENDDO
1550       ENDDO
1551!
1552!--    Calculate tendencies for the next Runge-Kutta step
1553       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1554          IF ( intermediate_timestep_count == 1 )  THEN
1555             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1556             !$ACC PRESENT(tend, tpt_m)
1557             DO  i = nxl, nxr
1558                DO  j = nys, nyn
1559                   DO  k = nzb+1, nzt
1560                      tpt_m(k,j,i) = tend(k,j,i)
1561                   ENDDO
1562                ENDDO
1563             ENDDO
1564          ELSEIF ( intermediate_timestep_count < &
1565                   intermediate_timestep_count_max )  THEN
1566             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
1567             !$ACC PRESENT(tend, tpt_m)
1568             DO  i = nxl, nxr
1569                DO  j = nys, nyn
1570                   DO  k = nzb+1, nzt
1571                      tpt_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +              &
1572                                        5.3125_wp * tpt_m(k,j,i)
1573                   ENDDO
1574                ENDDO
1575             ENDDO
1576          ENDIF
1577       ENDIF
1578
1579       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1580
1581    ENDIF
1582
1583!
1584!-- If required, compute prognostic equation for total water content
1585    IF ( humidity )  THEN
1586
1587       CALL cpu_log( log_point(29), 'q-equation', 'start' )
1588
1589!
1590!--    Scalar/q-tendency terms with communication
1591       sbt = tsc(2)
1592       IF ( scalar_advec == 'bc-scheme' )  THEN
1593
1594          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1595!
1596!--          Bott-Chlond scheme always uses Euler time step. Thus:
1597             sbt = 1.0_wp
1598          ENDIF
1599          tend = 0.0_wp
1600          CALL advec_s_bc( q, 'q' )
1601
1602       ENDIF
1603
1604!
1605!--    Scalar/q-tendency terms with no communication
1606       IF ( scalar_advec /= 'bc-scheme' )  THEN
1607          tend = 0.0_wp
1608          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1609             IF ( ws_scheme_sca )  THEN
1610                CALL advec_s_ws( advc_flags_s, q, 'q',                         &
1611                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1612                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1613                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1614                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1615             ELSE
1616                CALL advec_s_pw( q )
1617             ENDIF
1618          ELSE
1619             CALL advec_s_up( q )
1620          ENDIF
1621       ENDIF
1622
1623       CALL diffusion_s( q,                                                    &
1624                         surf_def_h(0)%qsws, surf_def_h(1)%qsws,               &
1625                         surf_def_h(2)%qsws,                                   &
1626                         surf_lsm_h%qsws,    surf_usm_h%qsws,                  &
1627                         surf_def_v(0)%qsws, surf_def_v(1)%qsws,               &
1628                         surf_def_v(2)%qsws, surf_def_v(3)%qsws,               &
1629                         surf_lsm_v(0)%qsws, surf_lsm_v(1)%qsws,               &
1630                         surf_lsm_v(2)%qsws, surf_lsm_v(3)%qsws,               &
1631                         surf_usm_v(0)%qsws, surf_usm_v(1)%qsws,               &
1632                         surf_usm_v(2)%qsws, surf_usm_v(3)%qsws )
1633
1634!
1635!--    Sink or source of humidity due to canopy elements
1636       IF ( plant_canopy ) CALL pcm_tendency( 5 )
1637
1638!
1639!--    Large scale advection
1640       IF ( large_scale_forcing )  THEN
1641          CALL ls_advec( simulated_time, 'q' )
1642       ENDIF
1643
1644!
1645!--    Nudging
1646       IF ( nudging )  CALL nudge( simulated_time, 'q' )
1647
1648!
1649!--    If required compute influence of large-scale subsidence/ascent
1650       IF ( large_scale_subsidence  .AND.                                      &
1651            .NOT. use_subsidence_tendencies )  THEN
1652         CALL subsidence( tend, q, q_init, 3 )
1653       ENDIF
1654
1655       CALL module_interface_actions( 'q-tendency' )
1656
1657!
1658!--    Prognostic equation for total water content
1659       DO  i = nxl, nxr
1660          DO  j = nys, nyn
1661             DO  k = nzb+1, nzt
1662                q_p(k,j,i) = q(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1663                                                 tsc(3) * tq_m(k,j,i) )        &
1664                                               - tsc(5) * rdf_sc(k) *          &
1665                                                      ( q(k,j,i) - q_init(k) ) &
1666                                        ) * MERGE( 1.0_wp, 0.0_wp,             &
1667                                                   BTEST( wall_flags_0(k,j,i), 0 ) &
1668                                                 )
1669                IF ( q_p(k,j,i) < 0.0_wp )  q_p(k,j,i) = 0.1_wp * q(k,j,i)
1670             ENDDO
1671          ENDDO
1672       ENDDO
1673
1674!
1675!--    Calculate tendencies for the next Runge-Kutta step
1676       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1677          IF ( intermediate_timestep_count == 1 )  THEN
1678             DO  i = nxl, nxr
1679                DO  j = nys, nyn
1680                   DO  k = nzb+1, nzt
1681                      tq_m(k,j,i) = tend(k,j,i)
1682                   ENDDO
1683                ENDDO
1684             ENDDO
1685          ELSEIF ( intermediate_timestep_count < &
1686                   intermediate_timestep_count_max )  THEN
1687             DO  i = nxl, nxr
1688                DO  j = nys, nyn
1689                   DO  k = nzb+1, nzt
1690                      tq_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1691                                     + 5.3125_wp * tq_m(k,j,i)
1692                   ENDDO
1693                ENDDO
1694             ENDDO
1695          ENDIF
1696       ENDIF
1697
1698       CALL cpu_log( log_point(29), 'q-equation', 'stop' )
1699
1700    ENDIF
1701!
1702!-- If required, compute prognostic equation for scalar
1703    IF ( passive_scalar )  THEN
1704
1705       CALL cpu_log( log_point(66), 's-equation', 'start' )
1706
1707!
1708!--    Scalar/q-tendency terms with communication
1709       sbt = tsc(2)
1710       IF ( scalar_advec == 'bc-scheme' )  THEN
1711
1712          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1713!
1714!--          Bott-Chlond scheme always uses Euler time step. Thus:
1715             sbt = 1.0_wp
1716          ENDIF
1717          tend = 0.0_wp
1718          CALL advec_s_bc( s, 's' )
1719
1720       ENDIF
1721
1722!
1723!--    Scalar/q-tendency terms with no communication
1724       IF ( scalar_advec /= 'bc-scheme' )  THEN
1725          tend = 0.0_wp
1726          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1727             IF ( ws_scheme_sca )  THEN
1728                CALL advec_s_ws( advc_flags_s, s, 's',                         &
1729                                 bc_dirichlet_l  .OR.  bc_radiation_l,         &
1730                                 bc_dirichlet_n  .OR.  bc_radiation_n,         &
1731                                 bc_dirichlet_r  .OR.  bc_radiation_r,         &
1732                                 bc_dirichlet_s  .OR.  bc_radiation_s )
1733             ELSE
1734                CALL advec_s_pw( s )
1735             ENDIF
1736          ELSE
1737             CALL advec_s_up( s )
1738          ENDIF
1739       ENDIF
1740
1741       CALL diffusion_s( s,                                                    &
1742                         surf_def_h(0)%ssws, surf_def_h(1)%ssws,               &
1743                         surf_def_h(2)%ssws,                                   &
1744                         surf_lsm_h%ssws,    surf_usm_h%ssws,                  &
1745                         surf_def_v(0)%ssws, surf_def_v(1)%ssws,               &
1746                         surf_def_v(2)%ssws, surf_def_v(3)%ssws,               &
1747                         surf_lsm_v(0)%ssws, surf_lsm_v(1)%ssws,               &
1748                         surf_lsm_v(2)%ssws, surf_lsm_v(3)%ssws,               &
1749                         surf_usm_v(0)%ssws, surf_usm_v(1)%ssws,               &
1750                         surf_usm_v(2)%ssws, surf_usm_v(3)%ssws )
1751
1752!
1753!--    Sink or source of humidity due to canopy elements
1754       IF ( plant_canopy ) CALL pcm_tendency( 7 )
1755
1756!
1757!--    Large scale advection. Not implemented for scalars so far.
1758!        IF ( large_scale_forcing )  THEN
1759!           CALL ls_advec( simulated_time, 'q' )
1760!        ENDIF
1761
1762!
1763!--    Nudging. Not implemented for scalars so far.
1764!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
1765
1766!
1767!--    If required compute influence of large-scale subsidence/ascent.
1768!--    Not implemented for scalars so far.
1769       IF ( large_scale_subsidence  .AND.                                      &
1770            .NOT. use_subsidence_tendencies  .AND.                             &
1771            .NOT. large_scale_forcing )  THEN
1772         CALL subsidence( tend, s, s_init, 3 )
1773       ENDIF
1774
1775       CALL module_interface_actions( 's-tendency' )
1776
1777!
1778!--    Prognostic equation for total water content
1779       DO  i = nxl, nxr
1780          DO  j = nys, nyn
1781             DO  k = nzb+1, nzt
1782                s_p(k,j,i) = s(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
1783                                                 tsc(3) * ts_m(k,j,i) )        &
1784                                               - tsc(5) * rdf_sc(k) *          &
1785                                                    ( s(k,j,i) - s_init(k) )   &
1786                                        )                                      &
1787                                   * MERGE( 1.0_wp, 0.0_wp,                    &
1788                                             BTEST( wall_flags_0(k,j,i), 0 )   &
1789                                          )
1790                IF ( s_p(k,j,i) < 0.0_wp )  s_p(k,j,i) = 0.1_wp * s(k,j,i)
1791             ENDDO
1792          ENDDO
1793       ENDDO
1794
1795!
1796!--    Calculate tendencies for the next Runge-Kutta step
1797       IF ( timestep_scheme(1:5) == 'runge' )  THEN
1798          IF ( intermediate_timestep_count == 1 )  THEN
1799             DO  i = nxl, nxr
1800                DO  j = nys, nyn
1801                   DO  k = nzb+1, nzt
1802                      ts_m(k,j,i) = tend(k,j,i)
1803                   ENDDO
1804                ENDDO
1805             ENDDO
1806          ELSEIF ( intermediate_timestep_count < &
1807                   intermediate_timestep_count_max )  THEN
1808             DO  i = nxl, nxr
1809                DO  j = nys, nyn
1810                   DO  k = nzb+1, nzt
1811                      ts_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
1812                                     + 5.3125_wp * ts_m(k,j,i)
1813                   ENDDO
1814                ENDDO
1815             ENDDO
1816          ENDIF
1817       ENDIF
1818
1819       CALL cpu_log( log_point(66), 's-equation', 'stop' )
1820
1821    ENDIF
1822!
1823!-- Calculate prognostic equations for all other modules
1824    CALL module_interface_prognostic_equations()
1825
1826    IF ( debug_output_timestep )  CALL debug_message( 'prognostic_equations_vector', 'end' )
1827
1828 END SUBROUTINE prognostic_equations_vector
1829
1830
1831 END MODULE prognostic_equations_mod
Note: See TracBrowser for help on using the repository browser.