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

Last change on this file since 2007 was 2007, checked in by kanani, 5 years ago

changes in the course of urban surface model implementation

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