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

Last change on this file since 1961 was 1961, checked in by suehring, 5 years ago

last commit documented

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