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

Last change on this file since 1140 was 1132, checked in by raasch, 11 years ago

r1028 documented

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