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

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago

last commit documented

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