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

Last change on this file since 1242 was 1242, checked in by heinze, 8 years ago

Last commmit documented

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