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

Last change on this file since 1375 was 1375, checked in by raasch, 7 years ago

last commit documented

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