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

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

all OpenACC directives and related parts removed from the code

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