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

Last change on this file since 1054 was 1054, checked in by hoffmann, 9 years ago

last commit documented

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