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

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

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

  • Property svn:keywords set to Id
File size: 67.2 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! update directives for prognostic quantities removed
23!
24! Former revisions:
25! -----------------
26! $Id: prognostic_equations.f90 1111 2013-03-08 23:54:10Z 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
1577!
1578!-- v-velocity component
1579    CALL cpu_log( log_point(6), 'v-equation', 'start' )
1580
1581    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1582       IF ( ws_scheme_mom )  THEN
1583          CALL advec_v_ws_acc
1584       ELSE
1585          tend = 0.0    ! to be removed later??
1586          CALL advec_v_pw
1587       END IF
1588    ELSE
1589       CALL advec_v_up
1590    ENDIF
1591    CALL diffusion_v_acc
1592    CALL coriolis_acc( 2 )
1593
1594!
1595!-- Drag by plant canopy
1596    IF ( plant_canopy )  CALL plant_canopy_model( 2 )
1597
1598!
1599!-- External pressure gradient
1600    IF ( dp_external )  THEN
1601       DO  i = nxl, nxr
1602          DO  j = nysv, nyn
1603             DO  k = dp_level_ind_b+1, nzt
1604                tend(k,j,i) = tend(k,j,i) - dpdxy(2) * dp_smooth_factor(k)
1605             ENDDO
1606          ENDDO
1607       ENDDO
1608    ENDIF
1609
1610    CALL user_actions( 'v-tendency' )
1611
1612!
1613!-- Prognostic equation for v-velocity component
1614    !$acc kernels present( nzb_v_inner, rdf, tend, tv_m, v, vg, v_p )
1615    !$acc loop
1616    DO  i = nxl, nxr
1617       DO  j = nysv, nyn
1618          !$acc loop vector( 32 )
1619          DO  k = 1, nzt
1620             IF ( k > nzb_v_inner(j,i) )  THEN
1621                v_p(k,j,i) = v(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
1622                                                  tsc(3) * tv_m(k,j,i) )       &
1623                                      - tsc(5) * rdf(k) * ( v(k,j,i) - vg(k) )
1624!
1625!--             Tendencies for the next Runge-Kutta step
1626                IF ( runge_step == 1 )  THEN
1627                   tv_m(k,j,i) = tend(k,j,i)
1628                ELSEIF ( runge_step == 2 )  THEN
1629                   tv_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tv_m(k,j,i)
1630                ENDIF
1631             ENDIF
1632          ENDDO
1633       ENDDO
1634    ENDDO
1635    !$acc end kernels
1636
1637    CALL cpu_log( log_point(6), 'v-equation', 'stop' )
1638
1639!
1640!-- w-velocity component
1641    CALL cpu_log( log_point(7), 'w-equation', 'start' )
1642
1643    IF ( timestep_scheme(1:5) == 'runge' )  THEN
1644       IF ( ws_scheme_mom )  THEN
1645          CALL advec_w_ws_acc
1646       ELSE
1647          tend = 0.0    ! to be removed later??
1648          CALL advec_w_pw
1649       ENDIF
1650    ELSE
1651       CALL advec_w_up
1652    ENDIF
1653    CALL diffusion_w_acc
1654    CALL coriolis_acc( 3 )
1655
1656    IF ( .NOT. neutral )  THEN
1657       IF ( ocean )  THEN
1658          CALL buoyancy( rho, rho_reference, 3, 64 )
1659       ELSE
1660          IF ( .NOT. humidity )  THEN
1661             CALL buoyancy_acc( pt, pt_reference, 3, 4 )
1662          ELSE
1663             CALL buoyancy( vpt, pt_reference, 3, 44 )
1664          ENDIF
1665       ENDIF
1666    ENDIF
1667
1668!
1669!-- Drag by plant canopy
1670    IF ( plant_canopy )  CALL plant_canopy_model( 3 )
1671
1672    CALL user_actions( 'w-tendency' )
1673
1674!
1675!-- Prognostic equation for w-velocity component
1676    !$acc kernels present( nzb_w_inner, rdf, tend, tw_m, w, w_p )
1677    !$acc loop
1678    DO  i = nxl, nxr
1679       DO  j = nys, nyn
1680          !$acc loop vector( 32 )
1681          DO  k = 1, nzt-1
1682             IF ( k > nzb_w_inner(j,i) )  THEN
1683                w_p(k,j,i) = w(k,j,i) + dt_3d * ( tsc(2) * tend(k,j,i) +       &
1684                                                  tsc(3) * tw_m(k,j,i) )       &
1685                                      - tsc(5) * rdf(k) * w(k,j,i)
1686   !
1687   !--          Tendencies for the next Runge-Kutta step
1688                IF ( runge_step == 1 )  THEN
1689                   tw_m(k,j,i) = tend(k,j,i)
1690                ELSEIF ( runge_step == 2 )  THEN
1691                   tw_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tw_m(k,j,i)
1692                ENDIF
1693             ENDIF
1694          ENDDO
1695       ENDDO
1696    ENDDO
1697    !$acc end kernels
1698
1699    CALL cpu_log( log_point(7), 'w-equation', 'stop' )
1700
1701
1702!
1703!-- If required, compute prognostic equation for potential temperature
1704    IF ( .NOT. neutral )  THEN
1705
1706       CALL cpu_log( log_point(13), 'pt-equation', 'start' )
1707
1708!
1709!--    pt-tendency terms with communication
1710       sbt = tsc(2)
1711       IF ( scalar_advec == 'bc-scheme' )  THEN
1712
1713          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1714!
1715!--          Bott-Chlond scheme always uses Euler time step. Thus:
1716             sbt = 1.0
1717          ENDIF
1718          tend = 0.0
1719          CALL advec_s_bc( pt, 'pt' )
1720
1721       ENDIF
1722
1723!
1724!--    pt-tendency terms with no communication
1725       IF ( scalar_advec /= 'bc-scheme' )  THEN
1726          tend = 0.0
1727          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1728             IF ( ws_scheme_sca )  THEN
1729                CALL advec_s_ws_acc( pt, 'pt' )
1730             ELSE
1731                tend = 0.0    ! to be removed later??
1732                CALL advec_s_pw( pt )
1733             ENDIF
1734          ELSE
1735             CALL advec_s_up( pt )
1736          ENDIF
1737       ENDIF
1738
1739       CALL diffusion_s_acc( pt, shf, tswst, wall_heatflux )
1740
1741!
1742!--    If required compute heating/cooling due to long wave radiation processes
1743       IF ( radiation )  THEN
1744          CALL calc_radiation
1745       ENDIF
1746
1747!
1748!--    If required compute impact of latent heat due to precipitation
1749       IF ( precipitation )  THEN
1750          CALL impact_of_latent_heat
1751       ENDIF
1752
1753!
1754!--    Consideration of heat sources within the plant canopy
1755       IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
1756          CALL plant_canopy_model( 4 )
1757       ENDIF
1758
1759!
1760!--    If required compute influence of large-scale subsidence/ascent
1761       IF ( large_scale_subsidence )  THEN
1762          CALL subsidence( tend, pt, pt_init )
1763       ENDIF
1764
1765       CALL user_actions( 'pt-tendency' )
1766
1767!
1768!--    Prognostic equation for potential temperature
1769       !$acc kernels present( nzb_s_inner, rdf_sc, ptdf_x, ptdf_y, pt_init ) &
1770       !$acc         present( tend, tpt_m, pt, pt_p )
1771       !$acc loop
1772       DO  i = nxl, nxr
1773          DO  j = nys, nyn
1774             !$acc loop vector( 32 )
1775             DO  k = 1, nzt
1776                IF ( k > nzb_s_inner(j,i) )  THEN
1777                   pt_p(k,j,i) = pt(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
1778                                                       tsc(3) * tpt_m(k,j,i) )    &
1779                                           - tsc(5) * ( pt(k,j,i) - pt_init(k) ) *&
1780                                             ( rdf_sc(k) + ptdf_x(i) + ptdf_y(j) )
1781!
1782!--                Tendencies for the next Runge-Kutta step
1783                   IF ( runge_step == 1 )  THEN
1784                      tpt_m(k,j,i) = tend(k,j,i)
1785                   ELSEIF ( runge_step == 2 )  THEN
1786                      tpt_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tpt_m(k,j,i)
1787                   ENDIF
1788                ENDIF
1789             ENDDO
1790          ENDDO
1791       ENDDO
1792       !$acc end kernels
1793
1794       CALL cpu_log( log_point(13), 'pt-equation', 'stop' )
1795
1796    ENDIF
1797
1798!
1799!-- If required, compute prognostic equation for salinity
1800    IF ( ocean )  THEN
1801
1802       CALL cpu_log( log_point(37), 'sa-equation', 'start' )
1803
1804!
1805!--    sa-tendency terms with communication
1806       sbt = tsc(2)
1807       IF ( scalar_advec == 'bc-scheme' )  THEN
1808
1809          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1810!
1811!--          Bott-Chlond scheme always uses Euler time step. Thus:
1812             sbt = 1.0
1813          ENDIF
1814          tend = 0.0
1815          CALL advec_s_bc( sa, 'sa' )
1816
1817       ENDIF
1818
1819!
1820!--    sa-tendency terms with no communication
1821       IF ( scalar_advec /= 'bc-scheme' )  THEN
1822          tend = 0.0
1823          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1824             IF ( ws_scheme_sca )  THEN
1825                 CALL advec_s_ws( sa, 'sa' )
1826             ELSE
1827                 CALL advec_s_pw( sa )
1828             ENDIF
1829          ELSE
1830             CALL advec_s_up( sa )
1831          ENDIF
1832       ENDIF
1833
1834       CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
1835
1836       CALL user_actions( 'sa-tendency' )
1837
1838!
1839!--    Prognostic equation for salinity
1840       DO  i = nxl, nxr
1841          DO  j = nys, nyn
1842             DO  k = nzb_s_inner(j,i)+1, nzt
1843                sa_p(k,j,i) = sa(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +        &
1844                                                    tsc(3) * tsa_m(k,j,i) )    &
1845                                        - tsc(5) * rdf_sc(k) *                 &
1846                                          ( sa(k,j,i) - sa_init(k) )
1847                IF ( sa_p(k,j,i) < 0.0 )  sa_p(k,j,i) = 0.1 * sa(k,j,i)
1848!
1849!--             Tendencies for the next Runge-Kutta step
1850                IF ( runge_step == 1 )  THEN
1851                   tsa_m(k,j,i) = tend(k,j,i)
1852                ELSEIF ( runge_step == 2 )  THEN
1853                   tsa_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tsa_m(k,j,i)
1854                ENDIF
1855             ENDDO
1856          ENDDO
1857       ENDDO
1858
1859       CALL cpu_log( log_point(37), 'sa-equation', 'stop' )
1860
1861!
1862!--    Calculate density by the equation of state for seawater
1863       CALL cpu_log( log_point(38), 'eqns-seawater', 'start' )
1864       CALL eqn_state_seawater
1865       CALL cpu_log( log_point(38), 'eqns-seawater', 'stop' )
1866
1867    ENDIF
1868
1869!
1870!-- If required, compute prognostic equation for total water content / scalar
1871    IF ( humidity  .OR.  passive_scalar )  THEN
1872
1873       CALL cpu_log( log_point(29), 'q/s-equation', 'start' )
1874
1875!
1876!--    Scalar/q-tendency terms with communication
1877       sbt = tsc(2)
1878       IF ( scalar_advec == 'bc-scheme' )  THEN
1879
1880          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1881!
1882!--          Bott-Chlond scheme always uses Euler time step. Thus:
1883             sbt = 1.0
1884          ENDIF
1885          tend = 0.0
1886          CALL advec_s_bc( q, 'q' )
1887
1888       ENDIF
1889
1890!
1891!--    Scalar/q-tendency terms with no communication
1892       IF ( scalar_advec /= 'bc-scheme' )  THEN
1893          tend = 0.0
1894          IF ( timestep_scheme(1:5) == 'runge' )  THEN
1895             IF ( ws_scheme_sca )  THEN
1896                CALL advec_s_ws( q, 'q' )
1897             ELSE
1898                CALL advec_s_pw( q )
1899             ENDIF
1900          ELSE
1901             CALL advec_s_up( q )
1902          ENDIF
1903       ENDIF
1904
1905       CALL diffusion_s( q, qsws, qswst, wall_qflux )
1906
1907!
1908!--    If required compute decrease of total water content due to
1909!--    precipitation
1910       IF ( precipitation )  THEN
1911          CALL calc_precipitation
1912       ENDIF
1913
1914!
1915!--    Sink or source of scalar concentration due to canopy elements
1916       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
1917
1918!
1919!--    If required compute influence of large-scale subsidence/ascent
1920       IF ( large_scale_subsidence )  THEN
1921         CALL subsidence( tend, q, q_init )
1922       ENDIF
1923
1924       CALL user_actions( 'q-tendency' )
1925
1926!
1927!--    Prognostic equation for total water content / scalar
1928       DO  i = nxl, nxr
1929          DO  j = nys, nyn
1930             DO  k = nzb_s_inner(j,i)+1, nzt
1931                q_p(k,j,i) = q(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
1932                                                  tsc(3) * tq_m(k,j,i) )       &
1933                                      - tsc(5) * rdf_sc(k) *                   &
1934                                        ( q(k,j,i) - q_init(k) )
1935                IF ( q_p(k,j,i) < 0.0 )  q_p(k,j,i) = 0.1 * q(k,j,i)
1936!
1937!--             Tendencies for the next Runge-Kutta step
1938                IF ( runge_step == 1 )  THEN
1939                   tq_m(k,j,i) = tend(k,j,i)
1940                ELSEIF ( runge_step == 2 )  THEN
1941                   tq_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * tq_m(k,j,i)
1942                ENDIF
1943             ENDDO
1944          ENDDO
1945       ENDDO
1946
1947       CALL cpu_log( log_point(29), 'q/s-equation', 'stop' )
1948
1949    ENDIF
1950
1951!
1952!-- If required, compute prognostic equation for turbulent kinetic
1953!-- energy (TKE)
1954    IF ( .NOT. constant_diffusion )  THEN
1955
1956       CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1957
1958!
1959!--    TKE-tendency terms with communication
1960       CALL production_e_init
1961
1962       sbt = tsc(2)
1963       IF ( .NOT. use_upstream_for_tke )  THEN
1964          IF ( scalar_advec == 'bc-scheme' )  THEN
1965
1966             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
1967!
1968!--             Bott-Chlond scheme always uses Euler time step. Thus:
1969                sbt = 1.0
1970             ENDIF
1971             tend = 0.0
1972             CALL advec_s_bc( e, 'e' )
1973
1974          ENDIF
1975       ENDIF
1976
1977!
1978!--    TKE-tendency terms with no communication
1979       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
1980          IF ( use_upstream_for_tke )  THEN
1981             tend = 0.0
1982             CALL advec_s_up( e )
1983          ELSE
1984             IF ( timestep_scheme(1:5) == 'runge' )  THEN
1985                IF ( ws_scheme_sca )  THEN
1986                   CALL advec_s_ws_acc( e, 'e' )
1987                ELSE
1988                   tend = 0.0    ! to be removed later??
1989                   CALL advec_s_pw( e )
1990                ENDIF
1991             ELSE
1992                tend = 0.0    ! to be removed later??
1993                CALL advec_s_up( e )
1994             ENDIF
1995          ENDIF
1996       ENDIF
1997
1998       IF ( .NOT. humidity )  THEN
1999          IF ( ocean )  THEN
2000             CALL diffusion_e( prho, prho_reference )
2001          ELSE
2002             CALL diffusion_e_acc( pt, pt_reference )
2003          ENDIF
2004       ELSE
2005          CALL diffusion_e( vpt, pt_reference )
2006       ENDIF
2007
2008       CALL production_e_acc
2009
2010!
2011!--    Additional sink term for flows through plant canopies
2012       IF ( plant_canopy )  CALL plant_canopy_model( 6 )
2013       CALL user_actions( 'e-tendency' )
2014
2015!
2016!--    Prognostic equation for TKE.
2017!--    Eliminate negative TKE values, which can occur due to numerical
2018!--    reasons in the course of the integration. In such cases the old TKE
2019!--    value is reduced by 90%.
2020       !$acc kernels present( e, e_p, nzb_s_inner, tend, te_m )
2021       !$acc loop
2022       DO  i = nxl, nxr
2023          DO  j = nys, nyn
2024             !$acc loop vector( 32 )
2025             DO  k = 1, nzt
2026                IF ( k > nzb_s_inner(j,i) )  THEN
2027                   e_p(k,j,i) = e(k,j,i) + dt_3d * ( sbt * tend(k,j,i) +          &
2028                                                     tsc(3) * te_m(k,j,i) )
2029                   IF ( e_p(k,j,i) < 0.0 )  e_p(k,j,i) = 0.1 * e(k,j,i)
2030!
2031!--                Tendencies for the next Runge-Kutta step
2032                   IF ( runge_step == 1 )  THEN
2033                      te_m(k,j,i) = tend(k,j,i)
2034                   ELSEIF ( runge_step == 2 )  THEN
2035                      te_m(k,j,i) = -9.5625 * tend(k,j,i) + 5.3125 * te_m(k,j,i)
2036                   ENDIF
2037                ENDIF
2038             ENDDO
2039          ENDDO
2040       ENDDO
2041       !$acc end kernels
2042
2043       CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
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.