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

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

last commit documented

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