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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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