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

Last change on this file since 2119 was 2119, checked in by raasch, 5 years ago

last commit documented

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