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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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