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

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

New:
---

Porting of FFT-solver for serial runs to GPU using CUDA FFT,
preprocessor lines in transpose routines rearranged, so that routines can also
be used in serial (non-parallel) mode,
transpositions also carried out in serial mode, routines fftx, fftxp replaced
by calls of fft_x, fft_x replaced by fft_x_1d in the 1D-decomposition routines
(Makefile, Makefile_check, fft_xy, poisfft, poisfft_hybrid, transpose, new: cuda_fft_interfaces)

--stdin argument for mpiexec on lckyuh, -y and -Y settings output to header (mrun)

Changed:


Module array_kind renamed precision_kind
(check_open, data_output_3d, fft_xy, modules, user_data_output_3d)

some format changes for coupled atmosphere-ocean runs (header)
small changes in code formatting (microphysics, prognostic_equations)

Errors:


bugfix: default value (0) assigned to coupling_start_time (modules)
bugfix: initial time for preruns of coupled runs is output as -coupling_start_time (data_output_profiles)

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