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

Last change on this file since 1380 was 1380, checked in by heinze, 10 years ago

Upper boundary conditions for pt and q in case of nudging adjusted

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