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

Last change on this file since 1496 was 1496, checked in by maronga, 9 years ago

added beta version of a land surface model and a simple radiation model for clear sky conditions

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