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

Last change on this file since 1366 was 1366, checked in by boeske, 7 years ago

last commit documented

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