source: palm/tags/release-3.5/SOURCE/production_e.f90 @ 3999

Last change on this file since 3999 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

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