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

Last change on this file since 1207 was 1182, checked in by raasch, 11 years ago

last commit documented, rc-file for example run updated

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