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

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

last commit documented

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