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

Last change on this file since 124 was 124, checked in by raasch, 14 years ago

further bugfixes for the ocean part

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