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

Last change on this file since 1318 was 1318, checked in by raasch, 10 years ago

former files/routines cpu_log and cpu_statistics combined to one module,
which also includes the former data module cpulog from the modules-file,
module interfaces removed

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