source: palm/trunk/SOURCE/production_e.f90 @ 2127

Last change on this file since 2127 was 2127, checked in by raasch, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 57.9 KB
Line 
1!> @file production_e.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: production_e.f90 2127 2017-01-20 16:06:40Z raasch $
27!
28! 2126 2017-01-20 15:54:21Z raasch
29! density in ocean case replaced by potential density
30!
31! 2118 2017-01-17 16:38:49Z raasch
32! OpenACC version of subroutine removed
33!
34! 2031 2016-10-21 15:11:58Z knoop
35! renamed variable rho to rho_ocean
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1873 2016-04-18 14:50:06Z maronga
41! Module renamed (removed _mod)
42!
43!
44! 1850 2016-04-08 13:29:27Z maronga
45! Module renamed
46!
47!
48! 1691 2015-10-26 16:17:44Z maronga
49! Renamed prandtl_layer to constant_flux_layer.
50!
51! 1682 2015-10-07 23:56:08Z knoop
52! Code annotations made doxygen readable
53!
54! 1374 2014-04-25 12:55:07Z raasch
55! nzb_s_outer removed from acc-present-list
56!
57! 1353 2014-04-08 15:21:23Z heinze
58! REAL constants provided with KIND-attribute
59!
60! 1342 2014-03-26 17:04:47Z kanani
61! REAL constants defined as wp-kind
62!
63! 1320 2014-03-20 08:40:49Z raasch
64! ONLY-attribute added to USE-statements,
65! kind-parameters added to all INTEGER and REAL declaration statements,
66! kinds are defined in new module kinds,
67! old module precision_kind is removed,
68! revision history before 2012 removed,
69! comment fields (!:) to be used for variable explanations added to
70! all variable declaration statements
71!
72! 1257 2013-11-08 15:18:40Z raasch
73! openacc loop and loop vector clauses removed, declare create moved after
74! the FORTRAN declaration statement
75!
76! 1179 2013-06-14 05:57:58Z raasch
77! use_reference renamed use_single_reference_value
78!
79! 1128 2013-04-12 06:19:32Z raasch
80! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
81! j_north
82!
83! 1036 2012-10-22 13:43:42Z raasch
84! code put under GPL (PALM 3.9)
85!
86! 1015 2012-09-27 09:23:24Z raasch
87! accelerator version (*_acc) added
88!
89! 1007 2012-09-19 14:30:36Z franke
90! Bugfix: calculation of buoyancy production has to consider the liquid water
91! mixing ratio in case of cloud droplets
92!
93! 940 2012-07-09 14:31:00Z raasch
94! TKE production by buoyancy can be switched off in case of runs with pure
95! neutral stratification
96!
97! Revision 1.1  1997/09/19 07:45:35  raasch
98! Initial revision
99!
100!
101! Description:
102! ------------
103!> Production terms (shear + buoyancy) of the TKE.
104!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
105!>          not considered well!
106!------------------------------------------------------------------------------!
107 MODULE production_e_mod
108 
109
110    USE wall_fluxes_mod,                                                       &
111        ONLY:  wall_fluxes_e
112
113    USE kinds
114
115    PRIVATE
116    PUBLIC production_e, production_e_init
117
118    LOGICAL, SAVE ::  first_call = .TRUE.  !<
119
120    REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0  !<
121    REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  v_0  !<
122
123    INTERFACE production_e
124       MODULE PROCEDURE production_e
125       MODULE PROCEDURE production_e_ij
126    END INTERFACE production_e
127   
128    INTERFACE production_e_init
129       MODULE PROCEDURE production_e_init
130    END INTERFACE production_e_init
131 
132 CONTAINS
133
134
135!------------------------------------------------------------------------------!
136! Description:
137! ------------
138!> Call for all grid points
139!------------------------------------------------------------------------------!
140    SUBROUTINE production_e
141
142       USE arrays_3d,                                                          &
143           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf,      &
144                  tend, tswst, u, v, vpt, w
145
146       USE cloud_parameters,                                                   &
147           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
148
149       USE control_parameters,                                                 &
150           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
151                  humidity, kappa, neutral, ocean, pt_reference,               &
152                  rho_reference, use_single_reference_value,                   &
153                  use_surface_fluxes, use_top_fluxes
154
155       USE grid_variables,                                                     &
156           ONLY:  ddx, dx, ddy, dy, wall_e_x, wall_e_y
157
158       USE indices,                                                            &
159           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner,                   &
160                   nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
161
162       IMPLICIT NONE
163
164       INTEGER(iwp) ::  i           !<
165       INTEGER(iwp) ::  j           !<
166       INTEGER(iwp) ::  k           !<
167
168       REAL(wp)     ::  def         !<
169       REAL(wp)     ::  dudx        !<
170       REAL(wp)     ::  dudy        !<
171       REAL(wp)     ::  dudz        !<
172       REAL(wp)     ::  dvdx        !<
173       REAL(wp)     ::  dvdy        !<
174       REAL(wp)     ::  dvdz        !<
175       REAL(wp)     ::  dwdx        !<
176       REAL(wp)     ::  dwdy        !<
177       REAL(wp)     ::  dwdz        !<
178       REAL(wp)     ::  k1          !<
179       REAL(wp)     ::  k2          !<
180       REAL(wp)     ::  km_neutral  !<
181       REAL(wp)     ::  theta       !<
182       REAL(wp)     ::  temp        !<
183
184!       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
185       REAL(wp), DIMENSION(nzb:nzt+1) ::  usvs  !<
186       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
187       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus  !<
188       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs  !<
189
190!
191!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
192!--    vertical walls, if neccessary
193!--    So far, results are slightly different from the ij-Version.
194!--    Therefore, ij-Version is called further below within the ij-loops.
195!       IF ( topography /= 'flat' )  THEN
196!          CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
197!          CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
198!          CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
199!          CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
200!       ENDIF
201
202
203       DO  i = nxl, nxr
204
205!
206!--       Calculate TKE production by shear
207          DO  j = nys, nyn
208             DO  k = nzb_diff_s_outer(j,i), nzt
209
210                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
211                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
212                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
213                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
214                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
215
216                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
217                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
218                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
219                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
220                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
221
222                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
223                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
224                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
225                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
226                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
227
228                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
229                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
230                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
231
232                IF ( def < 0.0_wp )  def = 0.0_wp
233
234                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
235
236             ENDDO
237          ENDDO
238
239          IF ( constant_flux_layer )  THEN
240
241!
242!--          Position beneath wall
243!--          (2) - Will allways be executed.
244!--          'bottom and wall: use u_0,v_0 and wall functions'
245             DO  j = nys, nyn
246
247                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
248                THEN
249
250                   k = nzb_diff_s_inner(j,i) - 1
251                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
252                   dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
253                                     u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
254                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
255                   dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
256                                     v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
257                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
258
259                   IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
260!
261!--                   Inconsistency removed: as the thermal stratification is
262!--                   not taken into account for the evaluation of the wall
263!--                   fluxes at vertical walls, the eddy viscosity km must not
264!--                   be used for the evaluation of the velocity gradients dudy
265!--                   and dwdy
266!--                   Note: The validity of the new method has not yet been
267!--                         shown, as so far no suitable data for a validation
268!--                         has been available
269                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
270                                          usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
271                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
272                                          wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
273                      km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
274                                   0.5_wp * dy 
275                      IF ( km_neutral > 0.0_wp )  THEN
276                         dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
277                         dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
278                      ELSE
279                         dudy = 0.0_wp
280                         dwdy = 0.0_wp
281                      ENDIF
282                   ELSE
283                      dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
284                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
285                      dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
286                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
287                   ENDIF
288
289                   IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
290!
291!--                   Inconsistency removed: as the thermal stratification is
292!--                   not taken into account for the evaluation of the wall
293!--                   fluxes at vertical walls, the eddy viscosity km must not
294!--                   be used for the evaluation of the velocity gradients dvdx
295!--                   and dwdx
296!--                   Note: The validity of the new method has not yet been
297!--                         shown, as so far no suitable data for a validation
298!--                         has been available
299                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
300                                          vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
301                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
302                                          wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
303                      km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * &
304                                   0.5_wp * dx
305                      IF ( km_neutral > 0.0_wp )  THEN
306                         dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
307                         dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
308                      ELSE
309                         dvdx = 0.0_wp
310                         dwdx = 0.0_wp
311                      ENDIF
312                   ELSE
313                      dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
314                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
315                      dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
316                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
317                   ENDIF
318
319                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
320                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
321                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
322
323                   IF ( def < 0.0_wp )  def = 0.0_wp
324
325                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
326
327
328!
329!--                (3) - will be executed only, if there is at least one level
330!--                between (2) and (4), i.e. the topography must have a
331!--                minimum height of 2 dz. Wall fluxes for this case have
332!--                already been calculated for (2).
333!--                'wall only: use wall functions'
334
335                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
336
337                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
338                      dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
339                                        u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
340                      dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
341                      dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
342                                        v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
343                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
344
345                      IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
346!
347!--                      Inconsistency removed: as the thermal stratification
348!--                      is not taken into account for the evaluation of the
349!--                      wall fluxes at vertical walls, the eddy viscosity km
350!--                      must not be used for the evaluation of the velocity
351!--                      gradients dudy and dwdy
352!--                      Note: The validity of the new method has not yet
353!--                            been shown, as so far no suitable data for a
354!--                            validation has been available
355                         km_neutral = kappa * ( usvs(k)**2 + & 
356                                                wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
357                         IF ( km_neutral > 0.0_wp )  THEN
358                            dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
359                            dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
360                         ELSE
361                            dudy = 0.0_wp
362                            dwdy = 0.0_wp
363                         ENDIF
364                      ELSE
365                         dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
366                                            u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
367                         dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
368                                            w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
369                      ENDIF
370
371                      IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
372!
373!--                      Inconsistency removed: as the thermal stratification
374!--                      is not taken into account for the evaluation of the
375!--                      wall fluxes at vertical walls, the eddy viscosity km
376!--                      must not be used for the evaluation of the velocity
377!--                      gradients dvdx and dwdx
378!--                      Note: The validity of the new method has not yet
379!--                            been shown, as so far no suitable data for a
380!--                            validation has been available
381                         km_neutral = kappa * ( vsus(k)**2 + & 
382                                                wsus(k)**2 )**0.25_wp * 0.5_wp * dx
383                         IF ( km_neutral > 0.0_wp )  THEN
384                            dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
385                            dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
386                         ELSE
387                            dvdx = 0.0_wp
388                            dwdx = 0.0_wp
389                         ENDIF
390                      ELSE
391                         dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
392                                            v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
393                         dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
394                                            w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
395                      ENDIF
396
397                      def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
398                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
399                           dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
400
401                      IF ( def < 0.0_wp )  def = 0.0_wp
402
403                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
404
405                   ENDDO
406
407                ENDIF
408
409             ENDDO
410
411!
412!--          (4) - will allways be executed.
413!--          'special case: free atmosphere' (as for case (0))
414             DO  j = nys, nyn
415
416                IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
417                THEN
418
419                   k = nzb_diff_s_outer(j,i)-1
420
421                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
422                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
423                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
424                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
425                                       u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
426
427                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
428                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
429                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
430                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
431                                       v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
432
433                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
434                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
435                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
436                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
437                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
438
439                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
440                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
441                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
442
443                   IF ( def < 0.0_wp )  def = 0.0_wp
444
445                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
446
447                ENDIF
448
449             ENDDO
450
451!
452!--          Position without adjacent wall
453!--          (1) - will allways be executed.
454!--          'bottom only: use u_0,v_0'
455             DO  j = nys, nyn
456
457                IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
458                THEN
459
460                   k = nzb_diff_s_inner(j,i)-1
461
462                   dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
463                   dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
464                                       u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
465                   dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
466                                       u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
467
468                   dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
469                                       v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
470                   dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
471                   dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
472                                       v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
473
474                   dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
475                                       w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
476                   dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
477                                       w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
478                   dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
479
480                   def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
481                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
482                         dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
483
484                   IF ( def < 0.0_wp )  def = 0.0_wp
485
486                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
487
488                ENDIF
489
490             ENDDO
491
492          ELSEIF ( use_surface_fluxes )  THEN
493
494             DO  j = nys, nyn
495
496                k = nzb_diff_s_outer(j,i)-1
497
498                dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
499                dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
500                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
501                dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
502                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
503
504                dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
505                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
506                dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
507                dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
508                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
509
510                dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
511                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
512                dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
513                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
514                dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
515
516                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
517                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
518                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
519
520                IF ( def < 0.0_wp )  def = 0.0_wp
521
522                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
523
524             ENDDO
525
526          ENDIF
527
528!
529!--       If required, calculate TKE production by buoyancy
530          IF ( .NOT. neutral )  THEN
531
532             IF ( .NOT. humidity )  THEN
533
534                IF ( use_single_reference_value )  THEN
535
536                   IF ( ocean )  THEN
537!
538!--                   So far in the ocean no special treatment of density flux
539!--                   in the bottom and top surface layer
540                      DO  j = nys, nyn
541                         DO  k = nzb_s_inner(j,i)+1, nzt
542                            tend(k,j,i) = tend(k,j,i) +                        &
543                                          kh(k,j,i) * g / rho_reference *      &
544                                          ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
545                                          dd2zu(k)
546                         ENDDO
547                      ENDDO
548
549                   ELSE
550
551                      DO  j = nys, nyn
552                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
553                            tend(k,j,i) = tend(k,j,i) -                   &
554                                          kh(k,j,i) * g / pt_reference *  &
555                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
556                                          dd2zu(k)
557                         ENDDO
558
559                         IF ( use_surface_fluxes )  THEN
560                            k = nzb_diff_s_inner(j,i)-1
561                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
562                                                        shf(j,i)
563                         ENDIF
564
565                         IF ( use_top_fluxes )  THEN
566                            k = nzt
567                            tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
568                                                        tswst(j,i)
569                         ENDIF
570                      ENDDO
571
572                   ENDIF
573
574                ELSE
575
576                   IF ( ocean )  THEN
577!
578!--                   So far in the ocean no special treatment of density flux
579!--                   in the bottom and top surface layer
580                      DO  j = nys, nyn
581                         DO  k = nzb_s_inner(j,i)+1, nzt
582                            tend(k,j,i) = tend(k,j,i) +                        &
583                                          kh(k,j,i) * g / prho(k,j,i) *        &
584                                          ( prho(k+1,j,i) - prho(k-1,j,i) ) *  &
585                                          dd2zu(k)
586                         ENDDO
587                      ENDDO
588
589                   ELSE
590
591                      DO  j = nys, nyn
592                         DO  k = nzb_diff_s_inner(j,i), nzt_diff
593                            tend(k,j,i) = tend(k,j,i) -                   &
594                                          kh(k,j,i) * g / pt(k,j,i) *     &
595                                          ( pt(k+1,j,i) - pt(k-1,j,i) ) * &
596                                          dd2zu(k)
597                         ENDDO
598
599                         IF ( use_surface_fluxes )  THEN
600                            k = nzb_diff_s_inner(j,i)-1
601                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
602                                                        shf(j,i)
603                         ENDIF
604
605                         IF ( use_top_fluxes )  THEN
606                            k = nzt
607                            tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
608                                                        tswst(j,i)
609                         ENDIF
610                      ENDDO
611
612                   ENDIF
613
614                ENDIF
615
616             ELSE
617
618                DO  j = nys, nyn
619
620                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
621
622                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
623                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
624                         k2 = 0.61_wp * pt(k,j,i)
625                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
626                                         g / vpt(k,j,i) *                      &
627                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
628                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
629                                         ) * dd2zu(k)
630                      ELSE IF ( cloud_physics )  THEN
631                         IF ( ql(k,j,i) == 0.0_wp )  THEN
632                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
633                            k2 = 0.61_wp * pt(k,j,i)
634                         ELSE
635                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
636                            temp  = theta * t_d_pt(k)
637                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
638                                          ( q(k,j,i) - ql(k,j,i) ) *          &
639                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
640                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
641                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
642                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
643                         ENDIF
644                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
645                                         g / vpt(k,j,i) *                      &
646                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
647                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
648                                         ) * dd2zu(k)
649                      ELSE IF ( cloud_droplets )  THEN
650                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
651                         k2 = 0.61_wp * pt(k,j,i)
652                         tend(k,j,i) = tend(k,j,i) -                          & 
653                                       kh(k,j,i) * g / vpt(k,j,i) *           &
654                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +  &
655                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -  &
656                                         pt(k,j,i) * ( ql(k+1,j,i) -          &
657                                         ql(k-1,j,i) ) ) * dd2zu(k)
658                      ENDIF
659
660                   ENDDO
661
662                ENDDO
663
664                IF ( use_surface_fluxes )  THEN
665
666                   DO  j = nys, nyn
667
668                      k = nzb_diff_s_inner(j,i)-1
669
670                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
671                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
672                         k2 = 0.61_wp * pt(k,j,i)
673                      ELSE IF ( cloud_physics )  THEN
674                         IF ( ql(k,j,i) == 0.0_wp )  THEN
675                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
676                            k2 = 0.61_wp * pt(k,j,i)
677                         ELSE
678                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
679                            temp  = theta * t_d_pt(k)
680                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
681                                          ( q(k,j,i) - ql(k,j,i) ) *           &
682                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
683                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
684                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
685                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
686                         ENDIF
687                      ELSE IF ( cloud_droplets )  THEN
688                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
689                         k2 = 0.61_wp * pt(k,j,i)
690                      ENDIF
691
692                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
693                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
694                   ENDDO
695
696                ENDIF
697
698                IF ( use_top_fluxes )  THEN
699
700                   DO  j = nys, nyn
701
702                      k = nzt
703
704                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
705                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
706                         k2 = 0.61_wp * pt(k,j,i)
707                      ELSE IF ( cloud_physics )  THEN
708                         IF ( ql(k,j,i) == 0.0_wp )  THEN
709                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
710                            k2 = 0.61_wp * pt(k,j,i)
711                         ELSE
712                            theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
713                            temp  = theta * t_d_pt(k)
714                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
715                                          ( q(k,j,i) - ql(k,j,i) ) *           &
716                                 ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /      &
717                                 ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *        &
718                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
719                            k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
720                         ENDIF
721                      ELSE IF ( cloud_droplets )  THEN
722                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
723                         k2 = 0.61_wp * pt(k,j,i)
724                      ENDIF
725
726                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
727                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
728                   ENDDO
729
730                ENDIF
731
732             ENDIF
733
734          ENDIF
735
736       ENDDO
737
738    END SUBROUTINE production_e
739
740
741!------------------------------------------------------------------------------!
742! Description:
743! ------------
744!> Call for grid point i,j
745!------------------------------------------------------------------------------!
746    SUBROUTINE production_e_ij( i, j )
747
748       USE arrays_3d,                                                          &
749           ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf,      &
750                  tend, tswst, u, v, vpt, w
751
752       USE cloud_parameters,                                                   &
753           ONLY:  l_d_cp, l_d_r, pt_d_t, t_d_pt
754
755       USE control_parameters,                                                 &
756           ONLY:  cloud_droplets, cloud_physics, constant_flux_layer, g,       &
757                  humidity, kappa, neutral, ocean, pt_reference,               &
758                  rho_reference, use_single_reference_value,                   &
759                  use_surface_fluxes, use_top_fluxes
760
761       USE grid_variables,                                                     &
762           ONLY:  ddx, dx, ddy, dy, wall_e_x, wall_e_y
763
764       USE indices,                                                            &
765           ONLY:  nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner,                   &
766                  nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
767
768       IMPLICIT NONE
769
770       INTEGER(iwp) ::  i           !<
771       INTEGER(iwp) ::  j           !<
772       INTEGER(iwp) ::  k           !<
773
774       REAL(wp)     ::  def         !<
775       REAL(wp)     ::  dudx        !<
776       REAL(wp)     ::  dudy        !<
777       REAL(wp)     ::  dudz        !<
778       REAL(wp)     ::  dvdx        !<
779       REAL(wp)     ::  dvdy        !<
780       REAL(wp)     ::  dvdz        !<
781       REAL(wp)     ::  dwdx        !<
782       REAL(wp)     ::  dwdy        !<
783       REAL(wp)     ::  dwdz        !<
784       REAL(wp)     ::  k1          !<
785       REAL(wp)     ::  k2          !<
786       REAL(wp)     ::  km_neutral  !<
787       REAL(wp)     ::  theta       !<
788       REAL(wp)     ::  temp        !<
789
790       REAL(wp), DIMENSION(nzb:nzt+1) ::  usvs  !<
791       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
792       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus  !<
793       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs  !<
794
795!
796!--    Calculate TKE production by shear
797       DO  k = nzb_diff_s_outer(j,i), nzt
798
799          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
800          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
801                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
802          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
803                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
804
805          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
806                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
807          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
808          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
809                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
810
811          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
812                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
813          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
814                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
815          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
816
817          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
818                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2    &
819                + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
820
821          IF ( def < 0.0_wp )  def = 0.0_wp
822
823          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
824
825       ENDDO
826
827       IF ( constant_flux_layer )  THEN
828
829          IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) )  THEN
830
831!
832!--          Position beneath wall
833!--          (2) - Will allways be executed.
834!--          'bottom and wall: use u_0,v_0 and wall functions'
835             k = nzb_diff_s_inner(j,i)-1
836
837             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
838             dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
839                               u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
840             dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
841             dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
842                               v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
843             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
844
845             IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
846!
847!--             Inconsistency removed: as the thermal stratification
848!--             is not taken into account for the evaluation of the
849!--             wall fluxes at vertical walls, the eddy viscosity km
850!--             must not be used for the evaluation of the velocity
851!--             gradients dudy and dwdy
852!--             Note: The validity of the new method has not yet
853!--                   been shown, as so far no suitable data for a
854!--                   validation has been available
855                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
856                                    usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
857                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
858                                    wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
859                km_neutral = kappa * ( usvs(k)**2 + wsvs(k)**2 )**0.25_wp * &
860                             0.5_wp * dy
861                IF ( km_neutral > 0.0_wp )  THEN
862                   dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
863                   dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
864                ELSE
865                   dudy = 0.0_wp
866                   dwdy = 0.0_wp
867                ENDIF
868             ELSE
869                dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
870                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
871                dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
872                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
873             ENDIF
874
875             IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
876!
877!--             Inconsistency removed: as the thermal stratification
878!--             is not taken into account for the evaluation of the
879!--             wall fluxes at vertical walls, the eddy viscosity km
880!--             must not be used for the evaluation of the velocity
881!--             gradients dvdx and dwdx
882!--             Note: The validity of the new method has not yet
883!--                   been shown, as so far no suitable data for a
884!--                   validation has been available
885                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
886                                    vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
887                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
888                                    wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
889                km_neutral = kappa * ( vsus(k)**2 + wsus(k)**2 )**0.25_wp * & 
890                             0.5_wp * dx
891                IF ( km_neutral > 0.0_wp )  THEN
892                   dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
893                   dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
894                ELSE
895                   dvdx = 0.0_wp
896                   dwdx = 0.0_wp
897                ENDIF
898             ELSE
899                dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
900                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
901                dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
902                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
903             ENDIF
904
905             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
906                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
907                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
908
909             IF ( def < 0.0_wp )  def = 0.0_wp
910
911             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
912
913!
914!--          (3) - will be executed only, if there is at least one level
915!--          between (2) and (4), i.e. the topography must have a
916!--          minimum height of 2 dz. Wall fluxes for this case have
917!--          already been calculated for (2).
918!--          'wall only: use wall functions'
919             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
920
921                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
922                dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
923                                  u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
924                dvdy =          ( v(k,j+1,i) - v(k,j,i)     ) * ddy
925                dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
926                                  v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
927                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
928
929                IF ( wall_e_y(j,i) /= 0.0_wp )  THEN
930!
931!--                Inconsistency removed: as the thermal stratification
932!--                is not taken into account for the evaluation of the
933!--                wall fluxes at vertical walls, the eddy viscosity km
934!--                must not be used for the evaluation of the velocity
935!--                gradients dudy and dwdy
936!--                Note: The validity of the new method has not yet
937!--                      been shown, as so far no suitable data for a
938!--                      validation has been available
939                   km_neutral = kappa * ( usvs(k)**2 + & 
940                                          wsvs(k)**2 )**0.25_wp * 0.5_wp * dy
941                   IF ( km_neutral > 0.0_wp )  THEN
942                      dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
943                      dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
944                   ELSE
945                      dudy = 0.0_wp
946                      dwdy = 0.0_wp
947                   ENDIF
948                ELSE
949                   dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
950                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
951                   dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
952                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
953                ENDIF
954
955                IF ( wall_e_x(j,i) /= 0.0_wp )  THEN
956!
957!--                Inconsistency removed: as the thermal stratification
958!--                is not taken into account for the evaluation of the
959!--                wall fluxes at vertical walls, the eddy viscosity km
960!--                must not be used for the evaluation of the velocity
961!--                gradients dvdx and dwdx
962!--                Note: The validity of the new method has not yet
963!--                      been shown, as so far no suitable data for a
964!--                      validation has been available
965                   km_neutral = kappa * ( vsus(k)**2 + & 
966                                          wsus(k)**2 )**0.25_wp * 0.5_wp * dx
967                   IF ( km_neutral > 0.0_wp )  THEN
968                      dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
969                      dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
970                   ELSE
971                      dvdx = 0.0_wp
972                      dwdx = 0.0_wp
973                   ENDIF
974                ELSE
975                   dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
976                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
977                   dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
978                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
979                ENDIF
980
981                def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
982                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
983                      dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
984
985                IF ( def < 0.0_wp )  def = 0.0_wp
986
987                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
988
989             ENDDO
990
991!
992!--          (4) - will allways be executed.
993!--          'special case: free atmosphere' (as for case (0))
994             k = nzb_diff_s_outer(j,i)-1
995
996             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
997             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
998                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
999             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1000                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1001
1002             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1003                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1004             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
1005             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1006                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1007
1008             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1009                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1010             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1011                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1012             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
1013
1014             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +        &
1015                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
1016                   dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
1017
1018             IF ( def < 0.0_wp )  def = 0.0_wp
1019
1020             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1021
1022          ELSE
1023
1024!
1025!--          Position without adjacent wall
1026!--          (1) - will allways be executed.
1027!--          'bottom only: use u_0,v_0'
1028             k = nzb_diff_s_inner(j,i)-1
1029
1030             dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
1031             dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1032                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1033             dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1034                                 u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
1035
1036             dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1037                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1038             dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
1039             dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1040                                 v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
1041
1042             dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1043                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1044             dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1045                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1046             dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
1047
1048             def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
1049                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
1050                   + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
1051
1052             IF ( def < 0.0_wp )  def = 0.0_wp
1053
1054             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1055
1056          ENDIF
1057
1058       ELSEIF ( use_surface_fluxes )  THEN
1059
1060          k = nzb_diff_s_outer(j,i)-1
1061
1062          dudx  =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
1063          dudy  = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1064                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1065          dudz  = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1066                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1067
1068          dvdx  = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1069                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1070          dvdy  =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
1071          dvdz  = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1072                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1073
1074          dwdx  = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1075                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1076          dwdy  = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1077                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1078          dwdz  =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
1079
1080          def = 2.0_wp * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
1081                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
1082                dvdz**2 + 2.0_wp * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
1083
1084          IF ( def < 0.0_wp )  def = 0.0_wp
1085
1086          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1087
1088       ENDIF
1089
1090!
1091!--    If required, calculate TKE production by buoyancy
1092       IF ( .NOT. neutral )  THEN
1093
1094          IF ( .NOT. humidity )  THEN
1095
1096             IF ( use_single_reference_value )  THEN
1097
1098                IF ( ocean )  THEN
1099!
1100!--                So far in the ocean no special treatment of density flux in
1101!--                the bottom and top surface layer
1102                   DO  k = nzb_s_inner(j,i)+1, nzt
1103                      tend(k,j,i) = tend(k,j,i) +                   &
1104                                    kh(k,j,i) * g / rho_reference * &
1105                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
1106                   ENDDO
1107
1108                ELSE
1109
1110                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
1111                      tend(k,j,i) = tend(k,j,i) -                  &
1112                                    kh(k,j,i) * g / pt_reference * &
1113                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1114                   ENDDO
1115
1116                   IF ( use_surface_fluxes )  THEN
1117                      k = nzb_diff_s_inner(j,i)-1
1118                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
1119                   ENDIF
1120
1121                   IF ( use_top_fluxes )  THEN
1122                      k = nzt
1123                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
1124                   ENDIF
1125
1126                ENDIF
1127
1128             ELSE
1129
1130                IF ( ocean )  THEN
1131!
1132!--                So far in the ocean no special treatment of density flux in
1133!--                the bottom and top surface layer
1134                   DO  k = nzb_s_inner(j,i)+1, nzt
1135                      tend(k,j,i) = tend(k,j,i) +                              &
1136                                    kh(k,j,i) * g / prho(k,j,i) *              &
1137                                    ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
1138                   ENDDO
1139
1140                ELSE
1141
1142                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
1143                      tend(k,j,i) = tend(k,j,i) -               &
1144                                    kh(k,j,i) * g / pt(k,j,i) * &
1145                                    ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1146                   ENDDO
1147
1148                   IF ( use_surface_fluxes )  THEN
1149                      k = nzb_diff_s_inner(j,i)-1
1150                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
1151                   ENDIF
1152
1153                   IF ( use_top_fluxes )  THEN
1154                      k = nzt
1155                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
1156                   ENDIF
1157
1158                ENDIF
1159
1160             ENDIF
1161
1162          ELSE
1163
1164             DO  k = nzb_diff_s_inner(j,i), nzt_diff
1165
1166                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
1167                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1168                   k2 = 0.61_wp * pt(k,j,i)
1169                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
1170                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1171                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
1172                                         ) * dd2zu(k)
1173                ELSE IF ( cloud_physics )  THEN
1174                   IF ( ql(k,j,i) == 0.0_wp )  THEN
1175                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1176                      k2 = 0.61_wp * pt(k,j,i)
1177                   ELSE
1178                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1179                      temp  = theta * t_d_pt(k)
1180                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
1181                                    ( q(k,j,i) - ql(k,j,i) ) *          &
1182                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
1183                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
1184                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1185                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1186                   ENDIF
1187                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
1188                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1189                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
1190                                         ) * dd2zu(k)
1191                ELSE IF ( cloud_droplets )  THEN
1192                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1193                   k2 = 0.61_wp * pt(k,j,i)
1194                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
1195                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
1196                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -    &
1197                                       pt(k,j,i) * ( ql(k+1,j,i) -           &
1198                                                     ql(k-1,j,i) ) ) * dd2zu(k)
1199                ENDIF
1200             ENDDO
1201
1202             IF ( use_surface_fluxes )  THEN
1203                k = nzb_diff_s_inner(j,i)-1
1204
1205                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
1206                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1207                   k2 = 0.61_wp * pt(k,j,i)
1208                ELSE IF ( cloud_physics )  THEN
1209                   IF ( ql(k,j,i) == 0.0_wp )  THEN
1210                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1211                      k2 = 0.61_wp * pt(k,j,i)
1212                   ELSE
1213                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1214                      temp  = theta * t_d_pt(k)
1215                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
1216                                    ( q(k,j,i) - ql(k,j,i) ) *          &
1217                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
1218                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
1219                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1220                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1221                   ENDIF
1222                ELSE IF ( cloud_droplets )  THEN
1223                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1224                   k2 = 0.61_wp * pt(k,j,i)
1225                ENDIF
1226
1227                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1228                                            ( k1* shf(j,i) + k2 * qsws(j,i) )
1229             ENDIF
1230
1231             IF ( use_top_fluxes )  THEN
1232                k = nzt
1233
1234                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
1235                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1236                   k2 = 0.61_wp * pt(k,j,i)
1237                ELSE IF ( cloud_physics )  THEN
1238                   IF ( ql(k,j,i) == 0.0_wp )  THEN
1239                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1240                      k2 = 0.61_wp * pt(k,j,i)
1241                   ELSE
1242                      theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1243                      temp  = theta * t_d_pt(k)
1244                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                 &
1245                                    ( q(k,j,i) - ql(k,j,i) ) *          &
1246                           ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /        &
1247                           ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp *          &
1248                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1249                      k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1250                   ENDIF
1251                ELSE IF ( cloud_droplets )  THEN
1252                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1253                   k2 = 0.61_wp * pt(k,j,i)
1254                ENDIF
1255
1256                tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1257                                            ( k1* tswst(j,i) + k2 * qswst(j,i) )
1258             ENDIF
1259
1260          ENDIF
1261
1262       ENDIF
1263
1264    END SUBROUTINE production_e_ij
1265
1266
1267!------------------------------------------------------------------------------!
1268! Description:
1269! ------------
1270!> @todo Missing subroutine description.
1271!------------------------------------------------------------------------------!
1272    SUBROUTINE production_e_init
1273
1274       USE arrays_3d,                                                          &
1275           ONLY:  kh, km, u, us, usws, v, vsws, zu
1276
1277       USE control_parameters,                                                 &
1278           ONLY:  constant_flux_layer, kappa
1279
1280       USE indices,                                                            &
1281           ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner,     &
1282                  nzb_v_inner
1283
1284       IMPLICIT NONE
1285
1286       INTEGER(iwp) ::  i   !<
1287       INTEGER(iwp) ::  j   !<
1288       INTEGER(iwp) ::  ku  !<
1289       INTEGER(iwp) ::  kv  !<
1290
1291       IF ( constant_flux_layer )  THEN
1292
1293          IF ( first_call )  THEN
1294             ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
1295             u_0 = 0.0_wp   ! just to avoid access of uninitialized memory
1296             v_0 = 0.0_wp   ! within exchange_horiz_2d
1297             first_call = .FALSE.
1298          ENDIF
1299
1300!
1301!--       Calculate a virtual velocity at the surface in a way that the
1302!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1303!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1304!--       production term at k=1 (see production_e_ij).
1305!--       The velocity gradient has to be limited in case of too small km
1306!--       (otherwise the timestep may be significantly reduced by large
1307!--       surface winds).
1308!--       Upper bounds are nxr+1 and nyn+1 because otherwise these values are
1309!--       not available in case of non-cyclic boundary conditions.
1310!--       WARNING: the exact analytical solution would require the determination
1311!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
1312          !$OMP PARALLEL DO PRIVATE( ku, kv )
1313          DO  i = nxl, nxr+1
1314             DO  j = nys, nyn+1
1315
1316                ku = nzb_u_inner(j,i)+1
1317                kv = nzb_v_inner(j,i)+1
1318
1319                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
1320                                 ( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) +    &
1321                                   1.0E-20_wp )
1322!                                  ( us(j,i) * kappa * zu(1) )
1323                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
1324                                 ( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) +    &
1325                                   1.0E-20_wp )
1326!                                  ( us(j,i) * kappa * zu(1) )
1327
1328                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
1329                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
1330                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
1331                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
1332
1333             ENDDO
1334          ENDDO
1335
1336          CALL exchange_horiz_2d( u_0 )
1337          CALL exchange_horiz_2d( v_0 )
1338
1339       ENDIF
1340
1341    END SUBROUTINE production_e_init
1342
1343 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.