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

Last change on this file since 1362 was 1362, checked in by hoffmann, 8 years ago

last commit documented

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