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

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