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

Last change on this file since 1365 was 1365, checked in by boeske, 7 years ago

large scale forcing enabled

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