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

Last change on this file since 1258 was 1258, checked in by raasch, 8 years ago

last commit documented

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