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

Last change on this file since 2001 was 2001, checked in by knoop, 8 years ago

last commit documented

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