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

Last change on this file since 1053 was 1053, checked in by hoffmann, 11 years ago

two-moment cloud physics implemented

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