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

Last change on this file since 1985 was 1977, checked in by maronga, 8 years ago

last commit documented

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