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

Last change on this file since 1092 was 1092, checked in by raasch, 9 years ago

unused variables remove from several routines

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