Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

Last change on this file since 1007 was 1007, checked in by franke, 12 years ago

Bugfixes

Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes

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