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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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