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

Last change on this file since 966 was 941, checked in by raasch, 12 years ago

last commit documented

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